3 #ifndef DUNE_FMATRIXEIGENVALUES_HH
4 #define DUNE_FMATRIXEIGENVALUES_HH
26 namespace FMatrixHelp {
30 const char* jobz,
const char* uplo,
const long
31 int* n,
double* a,
const long int* lda,
double* w,
32 double* work,
const long int* lwork,
long int* info);
35 const char* jobvl,
const char* jobvr,
const long
36 int* n,
double* a,
const long int* lda,
double* wr,
double* wi,
double* vl,
37 const long int* ldvl,
double* vr,
const long int* ldvr,
double* work,
38 const long int* lwork,
long int* info);
49 eigenvalues[0] = matrix[0][0];
62 const K detM = matrix[0][0] * matrix[1][1] - matrix[1][0] * matrix[0][1];
63 const K p = 0.5 * (matrix[0][0] + matrix [1][1]);
65 if( q < 0 && q > -1e-14 ) q = 0;
68 std::cout << matrix << std::endl;
70 DUNE_THROW(
MathError,
"Complex eigenvalue detected (which this implementation cannot handle).");
77 eigenvalues[0] = p - q;
78 eigenvalues[1] = p + q;
101 K p1 = matrix[0][1]*matrix[0][1] + matrix[0][2]*matrix[0][2] + matrix[1][2]*matrix[1][2];
106 eigenvalues[0] = matrix[0][0];
107 eigenvalues[1] = matrix[1][1];
108 eigenvalues[2] = matrix[2][2];
114 for (
int i=0; i<3; i++)
115 q += matrix[i][i]/3.0;
117 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;
121 for (
int i=0; i<3; i++)
122 for (
int j=0; j<3; j++)
123 B[i][j] = (real_type(1.0)/p) * (matrix[i][j] - q*(i==j));
138 eigenvalues[2] = q + 2 * p * cos(phi);
139 eigenvalues[0] = q + 2 * p * cos(phi + (2*pi/3));
140 eigenvalues[1] = 3 * q - eigenvalues[0] - eigenvalues[2];
153 template <
int dim,
typename K>
160 const long int N = dim ;
161 const char uplo =
'u';
164 const long int lwork = 3*N -1 ;
167 double matrixVector[dim * dim];
171 for(
int i=0; i<dim; ++i)
173 for(
int j=0; j<dim; ++j, ++row)
175 matrixVector[ row ] = matrix[ i ][ j ];
180 double workSpace[lwork];
192 for(
int i=0; i<dim; ++i)
194 for(
int j=0; j<dim; ++j, ++row)
196 eigenVectors[ i ][ j ] = matrixVector[ row ];
203 std::cerr <<
"For matrix " << matrix <<
" eigenvalue calculation failed! " << std::endl;
216 template <
int dim,
typename K>
233 template <
int dim,
typename K>
249 template <
int dim,
typename K,
class C>
254 const long int N = dim ;
255 const char jobvl =
'n';
256 const char jobvr =
'n';
259 double matrixVector[dim * dim];
263 for(
int i=0; i<dim; ++i)
265 for(
int j=0; j<dim; ++j, ++row)
267 matrixVector[ row ] = matrix[ i ][ j ];
278 const long int lwork = 3*dim;
282 &eigenR[0], &eigenI[0],
nullptr, &N,
nullptr, &N, &work[0],
287 std::cerr <<
"For matrix " << matrix <<
" eigenvalue calculation failed! " << std::endl;
290 for (
int i=0; i<N; ++i) {
A few common exception classes.
Implements a matrix constructed from a given type representing a field and compile-time given number ...
Implements a vector constructed from a given type representing a field and a compile-time given size.
Some useful basic math stuff.
#define DUNE_THROW(E, m)
Definition: exceptions.hh:216
Dune namespace.
Definition: alignedallocator.hh:14
void eigenValuesNonsymLapackCall(const char *jobvl, const char *jobvr, const long int *n, double *a, const long int *lda, double *wr, double *wi, double *vl, const long int *ldvl, double *vr, const long int *ldvr, double *work, const long int *lwork, long int *info)
Definition: fmatrixev.cc:223
static void eigenValuesNonSym(const FieldMatrix< K, dim, dim > &matrix, FieldVector< C, dim > &eigenValues)
calculates the eigenvalues of a non-symmetric field matrix
Definition: fmatrixev.hh:250
static void eigenValues(const FieldMatrix< K, 1, 1 > &matrix, FieldVector< K, 1 > &eigenvalues)
calculates the eigenvalues of a symmetric field matrix
Definition: fmatrixev.hh:46
static void eigenValuesVectors(const FieldMatrix< K, dim, dim > &matrix, FieldVector< K, dim > &eigenValues, FieldMatrix< K, dim, dim > &eigenVectors)
calculates the eigenvalues and eigenvectors of a symmetric field matrix
Definition: fmatrixev.hh:234
void eigenValuesLapackCall(const char *jobz, const char *uplo, const long int *n, double *a, const long int *lda, double *w, double *work, const long int *lwork, long int *info)
Definition: fmatrixev.cc:206
static void eigenValuesVectorsLapack(const FieldMatrix< K, dim, dim > &matrix, FieldVector< K, dim > &eigenValues, FieldMatrix< K, dim, dim > &eigenVectors, const char &jobz)
calculates the eigenvalues and/or eigenvectors of a symmetric field matrix
Definition: fmatrixev.hh:154
field_type determinant(bool doPivoting=true) const
calculates the determinant of this matrix
A dense n x m matrix.
Definition: fmatrix.hh:69
vector space out of a tensor product of fields.
Definition: fvector.hh:96
Default exception class for mathematical errors.
Definition: exceptions.hh:239
Default exception if a function was called while the object is not in a valid state for that function...
Definition: exceptions.hh:279
T real_type
export the type representing the real type of the field
Definition: ftraits.hh:28
static const Field pi()
Archimedes' constant.
Definition: math.hh:46