3#ifndef DUNE_FMATRIXEIGENVALUES_HH
4#define DUNE_FMATRIXEIGENVALUES_HH
27 namespace FMatrixHelp {
31 extern void eigenValuesLapackCall(
32 const char* jobz,
const char* uplo,
const long
33 int* n,
double* a,
const long int* lda,
double* w,
34 double* work,
const long int* lwork,
long int* info);
36 extern void eigenValuesNonsymLapackCall(
37 const char* jobvl,
const char* jobvr,
const long
38 int* n,
double* a,
const long int* lda,
double* wr,
double* wi,
double* vl,
39 const long int* ldvl,
double* vr,
const long int* ldvr,
double* work,
40 const long int* lwork,
long int* info);
42 extern void eigenValuesLapackCall(
43 const char* jobz,
const char* uplo,
const long
44 int* n,
float* a,
const long int* lda,
float* w,
45 float* work,
const long int* lwork,
long int* info);
47 extern void eigenValuesNonsymLapackCall(
48 const char* jobvl,
const char* jobvr,
const long
49 int* n,
float* a,
const long int* lda,
float* wr,
float* wi,
float* vl,
50 const long int* ldvl,
float* vr,
const long int* ldvr,
float* work,
51 const long int* lwork,
long int* info);
57 enum Jobs { OnlyEigenvalues=0, EigenvaluesEigenvectors=1 };
60 template<
typename K,
int dim>
61 using EVDummy = FieldMatrix<K, dim, dim>;
65 inline FieldVector<K,3> crossProduct(
const FieldVector<K,3>& vec0,
const FieldVector<K,3>& vec1) {
66 return {vec0[1]*vec1[2] - vec0[2]*vec1[1], vec0[2]*vec1[0] - vec0[0]*vec1[2], vec0[0]*vec1[1] - vec0[1]*vec1[0]};
70 static void eigenValues2dImpl(
const FieldMatrix<K, 2, 2>& matrix,
71 FieldVector<K, 2>& eigenvalues)
74 const K p = 0.5 * (matrix[0][0] + matrix [1][1]);
75 const K p2 = p - matrix[1][1];
76 K q = p2 * p2 + matrix[1][0] * matrix[0][1];
77 if( q < 0 && q > -1e-14 ) q = 0;
80 std::cout << matrix << std::endl;
82 DUNE_THROW(MathError,
"Complex eigenvalue detected (which this implementation cannot handle).");
89 eigenvalues[0] = p - q;
90 eigenvalues[1] = p + q;
100 template <
typename K>
101 static K eigenValues3dImpl(
const FieldMatrix<K, 3, 3>& matrix,
102 FieldVector<K, 3>& eigenvalues)
106 using real_type =
typename FieldTraits<K>::real_type;
108 K p1 = matrix[0][1]*matrix[0][1] + matrix[0][2]*matrix[0][2] + matrix[1][2]*matrix[1][2];
110 if (p1 <= std::numeric_limits<K>::epsilon()) {
112 eigenvalues[0] = matrix[0][0];
113 eigenvalues[1] = matrix[1][1];
114 eigenvalues[2] = matrix[2][2];
115 std::sort(eigenvalues.begin(), eigenvalues.end());
123 for (
int i=0; i<3; i++)
124 q += matrix[i][i] / 3.0;
126 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.0 * p1;
129 FieldMatrix<K,3,3> B;
130 for (
int i=0; i<3; i++)
131 for (
int j=0; j<3; j++)
132 B[i][j] = (real_type(1.0)/p) * (matrix[i][j] - q*(i==j));
134 K r = B.determinant() / 2.0;
141 r = std::min<K>(std::max<K>(r, -1.0), 1.0);
142 K phi = acos(r) / 3.0;
145 eigenvalues[2] = q + 2 * p * cos(phi);
146 eigenvalues[0] = q + 2 * p * cos(phi + (2*pi/3));
147 eigenvalues[1] = 3 * q - eigenvalues[0] - eigenvalues[2];
156 void orthoComp(
const FieldVector<K,3>& evec0, FieldVector<K,3>& u, FieldVector<K,3>& v) {
157 if(std::abs(evec0[0]) > std::abs(evec0[1])) {
159 FieldVector<K,2> temp = {evec0[0], evec0[2]};
160 auto L = 1.0 / temp.two_norm();
161 u = L * FieldVector<K,3>({-evec0[2], 0.0, evec0[0]});
165 FieldVector<K,2> temp = {evec0[1], evec0[2]};
166 auto L = 1.0 / temp.two_norm();
167 u = L * FieldVector<K,3>({0.0, evec0[2], -evec0[1]});
169 v = crossProduct(evec0, u);
174 void eig0(
const FieldMatrix<K,3,3>& matrix, K eval0, FieldVector<K,3>& evec0) {
180 using Vector = FieldVector<K,3>;
181 Vector row0 = {matrix[0][0]-eval0, matrix[0][1], matrix[0][2]};
182 Vector row1 = {matrix[1][0], matrix[1][1]-eval0, matrix[1][2]};
183 Vector row2 = {matrix[2][0], matrix[2][1], matrix[2][2]-eval0};
185 Vector r0xr1 = crossProduct(row0, row1);
186 Vector r0xr2 = crossProduct(row0, row2);
187 Vector r1xr2 = crossProduct(row1, row2);
188 auto d0 = r0xr1.two_norm();
189 auto d1 = r0xr2.two_norm();
190 auto d2 = r1xr2.two_norm();
211 void eig1(
const FieldMatrix<K,3,3>& matrix,
const FieldVector<K,3>& evec0, FieldVector<K,3>& evec1, K eval1) {
212 using Vector = FieldVector<K,3>;
216 orthoComp(evec0, u, v);
240 auto m00 = u.dot(Au) - eval1;
241 auto m01 = u.dot(Av);
242 auto m11 = v.dot(Av) - eval1;
249 auto absM00 = std::abs(m00);
250 auto absM01 = std::abs(m01);
251 auto absM11 = std::abs(m11);
252 if(absM00 >= absM11) {
253 auto maxAbsComp =
std::max(absM00, absM01);
254 if(maxAbsComp > 0.0) {
255 if(absM00 >= absM01) {
257 m00 = 1.0 / std::sqrt(1.0 + m01*m01);
262 m01 = 1.0 / std::sqrt(1.0 + m00*m00);
265 evec1 = m01*u - m00*v;
271 auto maxAbsComp =
std::max(absM11, absM01);
272 if(maxAbsComp > 0.0) {
273 if(absM11 >= absM01) {
275 m11 = 1.0 / std::sqrt(1.0 + m01*m01);
280 m01 = 1.0 / std::sqrt(1.0 + m11*m11);
283 evec1 = m11*u - m01*v;
291 template<Jobs Tag,
typename K>
292 static void eigenValuesVectorsImpl(
const FieldMatrix<K, 1, 1>& matrix,
293 FieldVector<K, 1>& eigenValues,
294 FieldMatrix<K, 1, 1>& eigenVectors)
297 if constexpr(Tag==EigenvaluesEigenvectors)
298 eigenVectors[0] = {1.0};
303 template <Jobs Tag,
typename K>
304 static void eigenValuesVectorsImpl(
const FieldMatrix<K, 2, 2>& matrix,
305 FieldVector<K, 2>& eigenValues,
306 FieldMatrix<K, 2, 2>& eigenVectors)
309 Impl::eigenValues2dImpl(matrix, eigenValues);
318 if constexpr(Tag==EigenvaluesEigenvectors) {
321 FieldMatrix<K,2,2> temp = matrix;
324 if(temp.infinity_norm() <= 1e-14) {
325 eigenVectors[0] = {1.0, 0.0};
326 eigenVectors[1] = {0.0, 1.0};
331 FieldVector<K,2> ev0 = {matrix[0][0]-
eigenValues[1], matrix[1][0]};
332 FieldVector<K,2> ev1 = {matrix[0][1], matrix[1][1]-
eigenValues[1]};
333 eigenVectors[0] = (ev0.two_norm2() >= ev1.two_norm2()) ? ev0/ev0.two_norm() : ev1/ev1.two_norm();
339 eigenVectors[1] = (ev0.two_norm2() >= ev1.two_norm2()) ? ev0/ev0.two_norm() : ev1/ev1.two_norm();
345 template <Jobs Tag,
typename K>
346 static void eigenValuesVectorsImpl(
const FieldMatrix<K, 3, 3>& matrix,
347 FieldVector<K, 3>& eigenValues,
348 FieldMatrix<K, 3, 3>& eigenVectors)
350 using Vector = FieldVector<K,3>;
351 using Matrix = FieldMatrix<K,3,3>;
358 K maxAbsElement = (isnormal(matrix.infinity_norm())) ? matrix.infinity_norm() : K(1.0);
359 Matrix scaledMatrix = matrix / maxAbsElement;
360 K r = Impl::eigenValues3dImpl(scaledMatrix, eigenValues);
362 if constexpr(Tag==EigenvaluesEigenvectors) {
363 K offDiagNorm = Vector{scaledMatrix[0][1],scaledMatrix[0][2],scaledMatrix[1][2]}.two_norm2();
364 if (offDiagNorm <= std::numeric_limits<K>::epsilon())
366 eigenValues = {scaledMatrix[0][0], scaledMatrix[1][1], scaledMatrix[2][2]};
367 eigenVectors = {{1.0, 0.0, 0.0}, {0.0, 1.0, 0.0}, {0.0, 0.0, 1.0}};
371 if (eigenValues[0] > eigenValues[1])
373 std::swap(eigenValues[0], eigenValues[1]);
374 std::swap(eigenVectors[0], eigenVectors[1]);
376 if (eigenValues[1] > eigenValues[2])
378 std::swap(eigenValues[1], eigenValues[2]);
379 std::swap(eigenVectors[1], eigenVectors[2]);
381 if (eigenValues[0] > eigenValues[1])
383 std::swap(eigenValues[0], eigenValues[1]);
384 std::swap(eigenVectors[0], eigenVectors[1]);
393 Vector eval(eigenValues);
395 Impl::eig0(scaledMatrix, eval[2], evec[2]);
396 Impl::eig1(scaledMatrix, evec[2], evec[1], eval[1]);
397 evec[0] = Impl::crossProduct(evec[1], evec[2]);
400 Impl::eig0(scaledMatrix, eval[0], evec[0]);
401 Impl::eig1(scaledMatrix, evec[0], evec[1], eval[1]);
402 evec[2] = Impl::crossProduct(evec[0], evec[1]);
405 using EVPair = std::pair<K, Vector>;
406 std::vector<EVPair> pairs;
407 for(std::size_t i=0; i<=2; ++i)
408 pairs.push_back(EVPair(eval[i], evec[i]));
409 auto comp = [](EVPair x, EVPair y){
return x.first < y.first; };
410 std::sort(pairs.begin(), pairs.end(), comp);
411 for(std::size_t i=0; i<=2; ++i){
413 eigenVectors[i] = pairs[i].second;
422 template <Jobs Tag,
int dim,
typename K>
423 static void eigenValuesVectorsLapackImpl(
const FieldMatrix<K, dim, dim>& matrix,
424 FieldVector<K, dim>& eigenValues,
425 FieldMatrix<K, dim, dim>& eigenVectors)
431 const char jobz =
"nv"[Tag];
433 const long int N = dim ;
434 const char uplo =
'u';
437 const long int lwork = 3*N -1 ;
439 constexpr bool isKLapackType = std::is_same_v<K,double> || std::is_same_v<K,float>;
440 using LapackNumType = std::conditional_t<isKLapackType, K, double>;
443 LapackNumType matrixVector[dim * dim];
447 for(
int i=0; i<dim; ++i)
449 for(
int j=0; j<dim; ++j, ++row)
451 matrixVector[ row ] = matrix[ i ][ j ];
456 LapackNumType workSpace[lwork];
461 if constexpr (isKLapackType){
464 ev =
new LapackNumType[dim];
468 eigenValuesLapackCall(&jobz, &uplo, &N, &matrixVector[0], &N,
469 ev, &workSpace[0], &lwork, &info);
471 if constexpr (!isKLapackType){
472 for(
size_t i=0;i<dim;++i)
473 eigenValues[i] = ev[i];
478 if (Tag==Jobs::EigenvaluesEigenvectors){
480 for(
int i=0; i<dim; ++i)
482 for(
int j=0; j<dim; ++j, ++row)
484 eigenVectors[ i ][ j ] = matrixVector[ row ];
491 std::cerr <<
"For matrix " << matrix <<
" eigenvalue calculation failed! " << std::endl;
492 DUNE_THROW(InvalidStateException,
"eigenValues: Eigenvalue calculation failed!");
495 DUNE_THROW(NotImplemented,
"LAPACK not found!");
501 template <Jobs Tag,
int dim,
typename K>
502 static void eigenValuesVectorsImpl(
const FieldMatrix<K, dim, dim>& matrix,
503 FieldVector<K, dim>& eigenValues,
504 FieldMatrix<K, dim, dim>& eigenVectors)
506 eigenValuesVectorsLapackImpl<Tag>(matrix,eigenValues,eigenVectors);
517 template <
int dim,
typename K>
522 Impl::eigenValuesVectorsImpl<Impl::Jobs::OnlyEigenvalues>(matrix,
eigenValues, dummy);
533 template <
int dim,
typename K>
538 Impl::eigenValuesVectorsImpl<Impl::Jobs::EigenvaluesEigenvectors>(matrix,
eigenValues, eigenVectors);
548 template <
int dim,
typename K>
553 Impl::eigenValuesVectorsLapackImpl<Impl::Jobs::EigenvaluesEigenvectors>(matrix,
eigenValues, dummy);
564 template <
int dim,
typename K>
569 Impl::eigenValuesVectorsLapackImpl<Impl::Jobs::EigenvaluesEigenvectors>(matrix,
eigenValues, eigenVectors);
579 template <
int dim,
typename K,
class C>
585 const long int N = dim ;
586 const char jobvl =
'n';
587 const char jobvr =
'n';
589 constexpr bool isKLapackType = std::is_same_v<K,double> || std::is_same_v<K,float>;
590 using LapackNumType = std::conditional_t<isKLapackType, K, double>;
593 LapackNumType matrixVector[dim * dim];
597 for(
int i=0; i<dim; ++i)
599 for(
int j=0; j<dim; ++j, ++row)
601 matrixVector[ row ] = matrix[ i ][ j ];
606 LapackNumType eigenR[dim];
607 LapackNumType eigenI[dim];
608 LapackNumType work[3*dim];
612 const long int lwork = 3*dim;
615 eigenValuesNonsymLapackCall(&jobvl, &jobvr, &N, &matrixVector[0], &N,
616 &eigenR[0], &eigenI[0],
nullptr, &N,
nullptr, &N, &work[0],
621 std::cerr <<
"For matrix " << matrix <<
" eigenvalue calculation failed! " << std::endl;
624 for (
int i=0; i<N; ++i) {
A dense n x m matrix.
Definition: fmatrix.hh:69
vector space out of a tensor product of fields.
Definition: fvector.hh:95
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 for dummy implementations.
Definition: exceptions.hh:261
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 non-symmetric field matrix
Definition: fmatrixev.hh:580
static void eigenValues(const FieldMatrix< K, dim, dim > &matrix, FieldVector< K, dim > &eigenValues)
calculates the eigenvalues of a symmetric field matrix
Definition: fmatrixev.hh:518
static void eigenValuesLapack(const FieldMatrix< K, dim, dim > &matrix, FieldVector< K, dim > &eigenValues)
calculates the eigenvalues of a symmetric field matrix
Definition: fmatrixev.hh:549
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:534
static void eigenValuesVectorsLapack(const FieldMatrix< K, dim, dim > &matrix, FieldVector< K, dim > &eigenValues, FieldMatrix< K, dim, dim > &eigenVectors)
calculates the eigenvalues and -vectors of a symmetric field matrix
Definition: fmatrixev.hh:565
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
auto max(ADLTag< 0 >, const V &v1, const V &v2)
implements binary Simd::max()
Definition: defaults.hh:79
Some useful basic math stuff.
Dune namespace.
Definition: alignedallocator.hh:11
static const Field pi()
Archimedes' constant.
Definition: math.hh:46