3#ifndef DUNE_FMATRIXEIGENVALUES_HH
4#define DUNE_FMATRIXEIGENVALUES_HH
25 namespace FMatrixHelp {
28 extern void eigenValuesLapackCall(
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);
33 extern void eigenValuesNonsymLapackCall(
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;
97 K p1 = matrix[0][1]*matrix[0][1] + matrix[0][2]*matrix[0][2] + matrix[1][2]*matrix[1][2];
102 eigenvalues[0] = matrix[0][0];
103 eigenvalues[1] = matrix[1][1];
104 eigenvalues[2] = matrix[2][2];
110 for (
int i=0; i<3; i++)
111 q += matrix[i][i]/3.0;
113 K p2 = (matrix[0][0] - q)*(matrix[0][0] - q) + (matrix[1][1] - q)*(matrix[1][1] - q) + (matrix[2][2] - q)*(matrix[2][2] - q) + 2 * p1;
114 K p = std::sqrt(p2 / 6);
117 for (
int i=0; i<3; i++)
118 for (
int j=0; j<3; j++)
119 B[i][j] = (1/p) * (matrix[i][j] - q*(i==j));
131 phi = std::acos(r) / 3;
134 eigenvalues[2] = q + 2 * p * cos(phi);
135 eigenvalues[0] = q + 2 * p * cos(phi + (2*M_PI/3));
136 eigenvalues[1] = 3 * q - eigenvalues[0] - eigenvalues[2];
147 template <
int dim,
typename K>
152 const long int N = dim ;
153 const char jobz =
'n';
154 const char uplo =
'u';
157 const long int w = N * N ;
160 double matrixVector[dim * dim];
164 for(
int i=0; i<dim; ++i)
166 for(
int j=0; j<dim; ++j, ++row)
168 matrixVector[ row ] = matrix[ i ][ j ];
173 double workSpace[dim * dim];
179 eigenValuesLapackCall(&jobz, &uplo, &N, &matrixVector[0], &N,
180 &eigenvalues[0], &workSpace[0], &w, &info);
184 std::cerr <<
"For matrix " << matrix <<
" eigenvalue calculation failed! " << std::endl;
196 template <
int dim,
typename K,
class C>
201 const long int N = dim ;
202 const char jobvl =
'n';
203 const char jobvr =
'n';
206 double matrixVector[dim * dim];
210 for(
int i=0; i<dim; ++i)
212 for(
int j=0; j<dim; ++j, ++row)
214 matrixVector[ row ] = matrix[ i ][ j ];
225 long int lwork = 3*dim;
228 eigenValuesNonsymLapackCall(&jobvl, &jobvr, &N, &matrixVector[0], &N,
229 &eigenR[0], &eigenI[0], 0, &N, 0, &N, &work[0],
234 std::cerr <<
"For matrix " << matrix <<
" eigenvalue calculation failed! " << std::endl;
237 for (
int i=0; i<N; ++i) {
field_type determinant() const
calculates the determinant of this matrix
A dense n x m matrix.
Definition: fmatrix.hh:68
vector space out of a tensor product of fields.
Definition: fvector.hh:93
Default exception if a function was called while the object is not in a valid state for that function...
Definition: exceptions.hh:279
A few common exception classes.
Implements a matrix constructed from a given type representing a field and compile-time given number ...
static void eigenValuesNonSym(const FieldMatrix< K, dim, dim > &matrix, FieldVector< C, dim > &eigenValues)
calculates the eigenvalues of a symmetric field matrix
Definition: fmatrixev.hh:197
static void eigenValues(const FieldMatrix< K, 1, 1 > &matrix, FieldVector< K, 1 > &eigenvalues)
calculates the eigenvalues of a symmetric field matrix
Definition: fmatrixev.hh:45
Implements a vector constructed from a given type representing a field and a compile-time given size.
#define DUNE_THROW(E, m)
Definition: exceptions.hh:216
Dune namespace.
Definition: alignment.hh:11