dune-common  2.2.0
fmatrixev.hh
Go to the documentation of this file.
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