1 #ifndef DUNE_FMATRIXEIGENVALUES_HH
2 #define DUNE_FMATRIXEIGENVALUES_HH
23 namespace FMatrixHelp {
27 const char* jobz,
const char* uplo,
const long
28 int* n,
double* a,
const long int* lda,
double* w,
29 double* work,
const long int* lwork,
long int* info);
40 eigenvalues[0] = matrix[0][0];
52 const K detM = matrix[0][0] * matrix[1][1] - matrix[1][0] * matrix[0][1];
53 const K p = 0.5 * (matrix[0][0] + matrix [1][1]);
55 if( q < 0 && q > -1e-14 ) q = 0;
58 std::cout << p <<
" p | q " << q <<
"\n";
59 std::cout << matrix << std::endl;
60 std::cout <<
"something went wrong in Eigenvalues for matrix!" << std::endl;
69 eigenvalues[0] = p - q;
70 eigenvalues[1] = p + q;
80 template <
int dim,
typename K>
85 const long int N = dim ;
86 const char jobz =
'n';
87 const char uplo =
'u';
90 const long int w = N * N ;
93 double matrixVector[dim * dim];
97 for(
int i=0; i<dim; ++i)
99 for(
int j=0; j<dim; ++j, ++row)
101 matrixVector[ row ] = matrix[ i ][ j ];
106 double workSpace[dim * dim];
113 &eigenvalues[0], &workSpace[0], &w, &info);
117 std::cerr <<
"For matrix " << matrix <<
" eigenvalue calculation failed! " << std::endl;