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