dune-common
2.2.0
|
00001 #ifndef DUNE_FMATRIXEIGENVALUES_HH 00002 #define DUNE_FMATRIXEIGENVALUES_HH 00003 00008 #include <iostream> 00009 #include <cmath> 00010 #include <cassert> 00011 00012 #include <dune/common/exceptions.hh> 00013 #include <dune/common/fvector.hh> 00014 #include <dune/common/fmatrix.hh> 00015 00016 namespace Dune { 00017 00023 namespace FMatrixHelp { 00024 00025 // defined in fmatrixev.cc 00026 extern void eigenValuesLapackCall( 00027 const char* jobz, const char* uplo, const long 00028 int* n, double* a, const long int* lda, double* w, 00029 double* work, const long int* lwork, long int* info); 00030 00036 template <typename K> 00037 static void eigenValues(const FieldMatrix<K, 1, 1>& matrix, 00038 FieldVector<K, 1>& eigenvalues) 00039 { 00040 eigenvalues[0] = matrix[0][0]; 00041 } 00042 00048 template <typename K> 00049 static void eigenValues(const FieldMatrix<K, 2, 2>& matrix, 00050 FieldVector<K, 2>& eigenvalues) 00051 { 00052 const K detM = matrix[0][0] * matrix[1][1] - matrix[1][0] * matrix[0][1]; 00053 const K p = 0.5 * (matrix[0][0] + matrix [1][1]); 00054 K q = p * p - detM; 00055 if( q < 0 && q > -1e-14 ) q = 0; 00056 if (p < 0 || q < 0) 00057 { 00058 std::cout << p << " p | q " << q << "\n"; 00059 std::cout << matrix << std::endl; 00060 std::cout << "something went wrong in Eigenvalues for matrix!" << std::endl; 00061 assert(false); 00062 abort(); 00063 } 00064 00065 // get square root 00066 q = std :: sqrt(q); 00067 00068 // store eigenvalues in ascending order 00069 eigenvalues[0] = p - q; 00070 eigenvalues[1] = p + q; 00071 } 00072 00080 template <int dim, typename K> 00081 static void eigenValues(const FieldMatrix<K, dim, dim>& matrix, 00082 FieldVector<K, dim>& eigenvalues) 00083 { 00084 { 00085 const long int N = dim ; 00086 const char jobz = 'n'; // only calculate eigenvalues 00087 const char uplo = 'u'; // use upper triangular matrix 00088 00089 // length of matrix vector 00090 const long int w = N * N ; 00091 00092 // matrix to put into dsyev 00093 double matrixVector[dim * dim]; 00094 00095 // copy matrix 00096 int row = 0; 00097 for(int i=0; i<dim; ++i) 00098 { 00099 for(int j=0; j<dim; ++j, ++row) 00100 { 00101 matrixVector[ row ] = matrix[ i ][ j ]; 00102 } 00103 } 00104 00105 // working memory 00106 double workSpace[dim * dim]; 00107 00108 // return value information 00109 long int info = 0; 00110 00111 // call LAPACK routine (see fmatrixev.cc) 00112 eigenValuesLapackCall(&jobz, &uplo, &N, &matrixVector[0], &N, 00113 &eigenvalues[0], &workSpace[0], &w, &info); 00114 00115 if( info != 0 ) 00116 { 00117 std::cerr << "For matrix " << matrix << " eigenvalue calculation failed! " << std::endl; 00118 DUNE_THROW(InvalidStateException,"eigenValues: Eigenvalue calculation failed!"); 00119 } 00120 } 00121 } 00122 00123 } // end namespace FMatrixHelp 00124 00127 } // end namespace Dune 00128 #endif