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 << matrix << std::endl;
68 DUNE_THROW(
MathError,
"Complex eigenvalue detected (which this implementation cannot handle).");
75 eigenvalues[0] = p - q;
76 eigenvalues[1] = p + q;
95 K p1 = matrix[0][1]*matrix[0][1] + matrix[0][2]*matrix[0][2] + matrix[1][2]*matrix[1][2];
100 eigenvalues[0] = matrix[0][0];
101 eigenvalues[1] = matrix[1][1];
102 eigenvalues[2] = matrix[2][2];
108 for (
int i=0; i<3; i++)
109 q += matrix[i][i]/3.0;
111 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;
112 K p = std::sqrt(p2 / 6);
115 for (
int i=0; i<3; i++)
116 for (
int j=0; j<3; j++)
117 B[i][j] = (1/p) * (matrix[i][j] - q*(i==j));
129 phi = std::acos(r) / 3;
132 eigenvalues[2] = q + 2 * p * cos(phi);
133 eigenvalues[0] = q + 2 * p * cos(phi + (2*M_PI/3));
134 eigenvalues[1] = 3 * q - eigenvalues[0] - eigenvalues[2];
145 template <
int dim,
typename K>
150 const long int N = dim ;
151 const char jobz =
'n';
152 const char uplo =
'u';
155 const long int w = N * N ;
158 double matrixVector[dim * dim];
162 for(
int i=0; i<dim; ++i)
164 for(
int j=0; j<dim; ++j, ++row)
166 matrixVector[ row ] = matrix[ i ][ j ];
171 double workSpace[dim * dim];
177 eigenValuesLapackCall(&jobz, &uplo, &N, &matrixVector[0], &N,
178 &eigenvalues[0], &workSpace[0], &w, &info);
182 std::cerr <<
"For matrix " << matrix <<
" eigenvalue calculation failed! " << std::endl;
194 template <
int dim,
typename K,
class C>
199 const long int N = dim ;
200 const char jobvl =
'n';
201 const char jobvr =
'n';
204 double matrixVector[dim * dim];
208 for(
int i=0; i<dim; ++i)
210 for(
int j=0; j<dim; ++j, ++row)
212 matrixVector[ row ] = matrix[ i ][ j ];
223 long int lwork = 3*dim;
226 eigenValuesNonsymLapackCall(&jobvl, &jobvr, &N, &matrixVector[0], &N,
227 &eigenR[0], &eigenI[0], 0, &N, 0, &N, &work[0],
232 std::cerr <<
"For matrix " << matrix <<
" eigenvalue calculation failed! " << std::endl;
235 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
Default exception class for mathematical errors.
Definition: exceptions.hh:239
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:195
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