3 #ifndef DUNE_FMATRIXEIGENVALUES_HH
4 #define DUNE_FMATRIXEIGENVALUES_HH
25 namespace FMatrixHelp {
29 const char* jobz,
const char* uplo,
const long
30 int* n,
double* a,
const long int* lda,
double* w,
31 double* work,
const long int* lwork,
long int* info);
34 const char* jobvl,
const char* jobvr,
const long
35 int* n,
double* a,
const long int* lda,
double* wr,
double* wi,
double* vl,
36 const long int* ldvl,
double* vr,
const long int* ldvr,
double* work,
37 const long int* lwork,
const long int* info);
48 eigenvalues[0] = matrix[0][0];
60 const K detM = matrix[0][0] * matrix[1][1] - matrix[1][0] * matrix[0][1];
61 const K p = 0.5 * (matrix[0][0] + matrix [1][1]);
63 if( q < 0 && q > -1e-14 ) q = 0;
66 std::cout << p <<
" p | q " << q <<
"\n";
67 std::cout << matrix << std::endl;
68 std::cout <<
"something went wrong in Eigenvalues for matrix!" << std::endl;
77 eigenvalues[0] = p - q;
78 eigenvalues[1] = p + q;
88 template <
int dim,
typename K>
93 const long int N = dim ;
94 const char jobz =
'n';
95 const char uplo =
'u';
98 const long int w = N * N ;
101 double matrixVector[dim * dim];
105 for(
int i=0; i<dim; ++i)
107 for(
int j=0; j<dim; ++j, ++row)
109 matrixVector[ row ] = matrix[ i ][ j ];
114 double workSpace[dim * dim];
121 &eigenvalues[0], &workSpace[0], &w, &info);
125 std::cerr <<
"For matrix " << matrix <<
" eigenvalue calculation failed! " << std::endl;
137 template <
int dim,
typename K,
class C>
142 const long int N = dim ;
143 const char jobvl =
'n';
144 const char jobvr =
'n';
147 double matrixVector[dim * dim];
151 for(
int i=0; i<dim; ++i)
153 for(
int j=0; j<dim; ++j, ++row)
155 matrixVector[ row ] = matrix[ i ][ j ];
166 long int lwork = 3*dim;
170 &eigenR[0], &eigenI[0], 0, &N, 0, &N, &work[0],
175 std::cerr <<
"For matrix " << matrix <<
" eigenvalue calculation failed! " << std::endl;
178 for (
int i=0; i<N; ++i) {
179 eigenValues[i].real = eigenR[i];
180 eigenValues[i].imag = eigenI[i];