1 // -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
2 // vi: set et ts=4 sw=2 sts=2:
3 #ifndef DUNE_DYNMATRIXEIGENVALUES_HH
4 #define DUNE_DYNMATRIXEIGENVALUES_HH
5 
6 #include <algorithm>
7 #include <memory>
8 
9 #include "dynmatrix.hh"
10 #include "fmatrixev.hh"
11 
12 /*!
13    \file
14    \brief utility functions to compute eigenvalues for
15    dense matrices.
16    \addtogroup DenseMatVec
17     @{
18  */
19 
20 namespace Dune {
21 
22   namespace DynamicMatrixHelp {
23 
24 #if HAVE_LAPACK
25     using Dune::FMatrixHelp::eigenValuesNonsymLapackCall;
26 #endif
27 
28     /** \brief calculates the eigenvalues of a symmetric field matrix
29         \param[in]  matrix matrix eigenvalues are calculated for
30         \param[out] eigenValues FieldVector that contains eigenvalues in
31                     ascending order
32         \param[out] eigenVectors (optional) list of right eigenvectors
33 
34         \note LAPACK::dgeev is used to calculate the eigen values
35      */
36     template <typename K, class C>
eigenValuesNonSym(const DynamicMatrix<K> & matrix,DynamicVector<C> & eigenValues,std::vector<DynamicVector<K>> * eigenVectors=nullptr)37     static void eigenValuesNonSym(const DynamicMatrix<K>& matrix,
38                                   DynamicVector<C>& eigenValues,
39                                   std::vector<DynamicVector<K>>* eigenVectors = nullptr
40       )
41     {
42 
43 #if HAVE_LAPACK
44       {
45         const long int N = matrix.rows();
46         const char jobvl = 'n';
47         const char jobvr = eigenVectors ? 'v' : 'n';
48 
49 
50         // matrix to put into dgeev
51         auto matrixVector = std::make_unique<double[]>(N*N);
52 
53         // copy matrix
54         int row = 0;
55         for(int i=0; i<N; ++i)
56         {
57           for(int j=0; j<N; ++j, ++row)
58           {
59             matrixVector[ row ] = matrix[ i ][ j ];
60           }
61         }
62 
63         // working memory
64         auto eigenR = std::make_unique<double[]>(N);
65         auto eigenI = std::make_unique<double[]>(N);
66 
67         const long int lwork = eigenVectors ? 4*N : 3*N;
68         auto work = std::make_unique<double[]>(lwork);
69         auto vr = eigenVectors ? std::make_unique<double[]>(N*N) : std::unique_ptr<double[]>{};
70 
71         // return value information
72         long int info = 0;
73 
74         // call LAPACK routine (see fmatrixev_ext.cc)
75         eigenValuesNonsymLapackCall(&jobvl, &jobvr, &N, matrixVector.get(), &N,
76                                     eigenR.get(), eigenI.get(), nullptr, &N, vr.get(), &N, work.get(),
77                                     &lwork, &info);
78 
79         if( info != 0 )
80         {
81           std::cerr << "For matrix " << matrix << " eigenvalue calculation failed! " << std::endl;
82           DUNE_THROW(InvalidStateException,"eigenValues: Eigenvalue calculation failed!");
83         }
84 
85         eigenValues.resize(N);
86         for (int i=0; i<N; ++i)
87           eigenValues[i] = std::complex<double>(eigenR[i], eigenI[i]);
88 
89         if (eigenVectors) {
90           eigenVectors->resize(N);
91           for (int i = 0; i < N; ++i) {
92             auto& v = (*eigenVectors)[i];
93             v.resize(N);
94             std::copy(vr.get() + N*i, vr.get() + N*(i+1), &v[0]);
95           }
96         }
97       }
98 #else // #if HAVE_LAPACK
99       DUNE_THROW(NotImplemented,"LAPACK not found!");
100 #endif
101     }
102   }
103 
104 }
105 /** @} */
106 #endif
107