5#ifndef DUNE_FMATRIXEIGENVALUES_HH
6#define DUNE_FMATRIXEIGENVALUES_HH
17#include <dune-common-config.hh>
30 namespace FMatrixHelp {
34 extern void eigenValuesLapackCall(
35 const char* jobz,
const char* uplo,
const long
36 int* n,
double* a,
const long int* lda,
double* w,
37 double* work,
const long int* lwork,
long int* info);
39 extern void eigenValuesNonsymLapackCall(
40 const char* jobvl,
const char* jobvr,
const long
41 int* n,
double* a,
const long int* lda,
double* wr,
double* wi,
double* vl,
42 const long int* ldvl,
double* vr,
const long int* ldvr,
double* work,
43 const long int* lwork,
long int* info);
45 extern void eigenValuesLapackCall(
46 const char* jobz,
const char* uplo,
const long
47 int* n,
float* a,
const long int* lda,
float* w,
48 float* work,
const long int* lwork,
long int* info);
50 extern void eigenValuesNonsymLapackCall(
51 const char* jobvl,
const char* jobvr,
const long
52 int* n,
float* a,
const long int* lda,
float* wr,
float* wi,
float* vl,
53 const long int* ldvl,
float* vr,
const long int* ldvr,
float* work,
54 const long int* lwork,
long int* info);
60 enum Jobs { OnlyEigenvalues=0, EigenvaluesEigenvectors=1 };
63 template<
typename K,
int dim>
64 using EVDummy = FieldMatrix<K, dim, dim>;
68 inline FieldVector<K,3> crossProduct(
const FieldVector<K,3>& vec0,
const FieldVector<K,3>& vec1) {
69 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]};
73 static void eigenValues2dImpl(
const FieldMatrix<K, 2, 2>& matrix,
74 FieldVector<K, 2>& eigenvalues)
77 const K p = 0.5 * (matrix[0][0] + matrix [1][1]);
78 const K p2 = p - matrix[1][1];
79 K q = p2 * p2 + matrix[1][0] * matrix[0][1];
80 if( q < 0 && q > -1e-14 ) q = 0;
83 std::cout << matrix << std::endl;
85 DUNE_THROW(MathError,
"Complex eigenvalue detected (which this implementation cannot handle).");
92 eigenvalues[0] = p - q;
93 eigenvalues[1] = p + q;
103 template <
typename K>
104 static K eigenValues3dImpl(
const FieldMatrix<K, 3, 3>& matrix,
105 FieldVector<K, 3>& eigenvalues)
109 using real_type =
typename FieldTraits<K>::real_type;
111 K p1 = matrix[0][1]*matrix[0][1] + matrix[0][2]*matrix[0][2] + matrix[1][2]*matrix[1][2];
113 if (p1 <= std::numeric_limits<K>::epsilon()) {
115 eigenvalues[0] = matrix[0][0];
116 eigenvalues[1] = matrix[1][1];
117 eigenvalues[2] = matrix[2][2];
118 std::sort(eigenvalues.begin(), eigenvalues.end());
126 for (
int i=0; i<3; i++)
127 q += matrix[i][i] / 3.0;
129 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;
132 FieldMatrix<K,3,3> B;
133 for (
int i=0; i<3; i++)
134 for (
int j=0; j<3; j++)
135 B[i][j] = (real_type(1.0)/p) * (matrix[i][j] - q*(i==j));
137 K r = B.determinant() / 2.0;
145 r = clamp<K>(r, -1.0, 1.0);
146 K phi = acos(r) / 3.0;
149 eigenvalues[2] = q + 2 * p * cos(phi);
150 eigenvalues[0] = q + 2 * p * cos(phi + (2*pi/3));
151 eigenvalues[1] = 3 * q - eigenvalues[0] - eigenvalues[2];
160 void orthoComp(
const FieldVector<K,3>& evec0, FieldVector<K,3>& u, FieldVector<K,3>& v) {
162 if(abs(evec0[0]) > abs(evec0[1])) {
164 FieldVector<K,2> temp = {evec0[0], evec0[2]};
165 auto L = 1.0 / temp.two_norm();
166 u = L * FieldVector<K,3>({-evec0[2], 0.0, evec0[0]});
170 FieldVector<K,2> temp = {evec0[1], evec0[2]};
171 auto L = 1.0 / temp.two_norm();
172 u = L * FieldVector<K,3>({0.0, evec0[2], -evec0[1]});
174 v = crossProduct(evec0, u);
179 void eig0(
const FieldMatrix<K,3,3>& matrix, K eval0, FieldVector<K,3>& evec0) {
185 using Vector = FieldVector<K,3>;
186 Vector row0 = {matrix[0][0]-eval0, matrix[0][1], matrix[0][2]};
187 Vector row1 = {matrix[1][0], matrix[1][1]-eval0, matrix[1][2]};
188 Vector row2 = {matrix[2][0], matrix[2][1], matrix[2][2]-eval0};
190 Vector r0xr1 = crossProduct(row0, row1);
191 Vector r0xr2 = crossProduct(row0, row2);
192 Vector r1xr2 = crossProduct(row1, row2);
193 auto d0 = r0xr1.two_norm();
194 auto d1 = r0xr2.two_norm();
195 auto d2 = r1xr2.two_norm();
216 void eig1(
const FieldMatrix<K,3,3>& matrix,
const FieldVector<K,3>& evec0, FieldVector<K,3>& evec1, K eval1) {
217 using Vector = FieldVector<K,3>;
221 orthoComp(evec0, u, v);
245 auto m00 = u.dot(Au) - eval1;
246 auto m01 = u.dot(Av);
247 auto m11 = v.dot(Av) - eval1;
254 using std::abs, std::sqrt, std::max;
255 auto absM00 = abs(m00);
256 auto absM01 = abs(m01);
257 auto absM11 = abs(m11);
258 if(absM00 >= absM11) {
259 auto maxAbsComp =
max(absM00, absM01);
260 if(maxAbsComp > 0.0) {
261 if(absM00 >= absM01) {
263 m00 = 1.0 / sqrt(1.0 + m01*m01);
268 m01 = 1.0 / sqrt(1.0 + m00*m00);
271 evec1 = m01*u - m00*v;
277 auto maxAbsComp =
max(absM11, absM01);
278 if(maxAbsComp > 0.0) {
279 if(absM11 >= absM01) {
281 m11 = 1.0 / sqrt(1.0 + m01*m01);
286 m01 = 1.0 / sqrt(1.0 + m11*m11);
289 evec1 = m11*u - m01*v;
297 template<Jobs Tag,
typename K>
298 static void eigenValuesVectorsImpl(
const FieldMatrix<K, 1, 1>& matrix,
299 FieldVector<K, 1>& eigenValues,
300 FieldMatrix<K, 1, 1>& eigenVectors)
303 if constexpr(Tag==EigenvaluesEigenvectors)
304 eigenVectors[0] = {1.0};
309 template <Jobs Tag,
typename K>
310 static void eigenValuesVectorsImpl(
const FieldMatrix<K, 2, 2>& matrix,
311 FieldVector<K, 2>& eigenValues,
312 FieldMatrix<K, 2, 2>& eigenVectors)
315 Impl::eigenValues2dImpl(matrix, eigenValues);
324 if constexpr(Tag==EigenvaluesEigenvectors) {
327 FieldMatrix<K,2,2> temp = matrix;
330 if(temp.infinity_norm() <= 1e-14) {
331 eigenVectors[0] = {1.0, 0.0};
332 eigenVectors[1] = {0.0, 1.0};
337 FieldVector<K,2> ev0 = {matrix[0][0]-
eigenValues[1], matrix[1][0]};
338 FieldVector<K,2> ev1 = {matrix[0][1], matrix[1][1]-
eigenValues[1]};
339 eigenVectors[0] = (ev0.two_norm2() >= ev1.two_norm2()) ? ev0/ev0.two_norm() : ev1/ev1.two_norm();
345 eigenVectors[1] = (ev0.two_norm2() >= ev1.two_norm2()) ? ev0/ev0.two_norm() : ev1/ev1.two_norm();
351 template <Jobs Tag,
typename K>
352 static void eigenValuesVectorsImpl(
const FieldMatrix<K, 3, 3>& matrix,
353 FieldVector<K, 3>& eigenValues,
354 FieldMatrix<K, 3, 3>& eigenVectors)
356 using Vector = FieldVector<K,3>;
357 using Matrix = FieldMatrix<K,3,3>;
364 K maxAbsElement = (isnormal(matrix.infinity_norm())) ? matrix.infinity_norm() : K(1.0);
365 Matrix scaledMatrix = matrix / maxAbsElement;
366 K r = Impl::eigenValues3dImpl(scaledMatrix, eigenValues);
368 if constexpr(Tag==EigenvaluesEigenvectors) {
369 K offDiagNorm = Vector{scaledMatrix[0][1],scaledMatrix[0][2],scaledMatrix[1][2]}.two_norm2();
370 if (offDiagNorm <= std::numeric_limits<K>::epsilon())
372 eigenValues = {scaledMatrix[0][0], scaledMatrix[1][1], scaledMatrix[2][2]};
373 eigenVectors = {{1.0, 0.0, 0.0}, {0.0, 1.0, 0.0}, {0.0, 0.0, 1.0}};
377 if (eigenValues[0] > eigenValues[1])
379 std::swap(eigenValues[0], eigenValues[1]);
380 std::swap(eigenVectors[0], eigenVectors[1]);
382 if (eigenValues[1] > eigenValues[2])
384 std::swap(eigenValues[1], eigenValues[2]);
385 std::swap(eigenVectors[1], eigenVectors[2]);
387 if (eigenValues[0] > eigenValues[1])
389 std::swap(eigenValues[0], eigenValues[1]);
390 std::swap(eigenVectors[0], eigenVectors[1]);
399 Vector eval(eigenValues);
401 Impl::eig0(scaledMatrix, eval[2], evec[2]);
402 Impl::eig1(scaledMatrix, evec[2], evec[1], eval[1]);
403 evec[0] = Impl::crossProduct(evec[1], evec[2]);
406 Impl::eig0(scaledMatrix, eval[0], evec[0]);
407 Impl::eig1(scaledMatrix, evec[0], evec[1], eval[1]);
408 evec[2] = Impl::crossProduct(evec[0], evec[1]);
411 using EVPair = std::pair<K, Vector>;
412 std::vector<EVPair> pairs;
413 for(std::size_t i=0; i<=2; ++i)
414 pairs.push_back(EVPair(eval[i], evec[i]));
415 auto comp = [](EVPair x, EVPair y){
return x.first < y.first; };
416 std::sort(pairs.begin(), pairs.end(), comp);
417 for(std::size_t i=0; i<=2; ++i){
419 eigenVectors[i] = pairs[i].second;
428 template <Jobs Tag,
int dim,
typename K>
429 static void eigenValuesVectorsLapackImpl(
const FieldMatrix<K, dim, dim>& matrix,
430 FieldVector<K, dim>& eigenValues,
431 FieldMatrix<K, dim, dim>& eigenVectors)
437 const char jobz =
"nv"[Tag];
439 const long int N = dim ;
440 const char uplo =
'u';
443 const long int lwork = 3*N -1 ;
445 constexpr bool isKLapackType = std::is_same_v<K,double> || std::is_same_v<K,float>;
446 using LapackNumType = std::conditional_t<isKLapackType, K, double>;
449 LapackNumType matrixVector[dim * dim];
453 for(
int i=0; i<dim; ++i)
455 for(
int j=0; j<dim; ++j, ++row)
457 matrixVector[ row ] = matrix[ i ][ j ];
462 LapackNumType workSpace[lwork];
467 if constexpr (isKLapackType){
470 ev =
new LapackNumType[dim];
474 eigenValuesLapackCall(&jobz, &uplo, &N, &matrixVector[0], &N,
475 ev, &workSpace[0], &lwork, &info);
477 if constexpr (!isKLapackType){
478 for(
size_t i=0;i<dim;++i)
479 eigenValues[i] = ev[i];
484 if (Tag==Jobs::EigenvaluesEigenvectors){
486 for(
int i=0; i<dim; ++i)
488 for(
int j=0; j<dim; ++j, ++row)
490 eigenVectors[ i ][ j ] = matrixVector[ row ];
497 std::cerr <<
"For matrix " << matrix <<
" eigenvalue calculation failed! " << std::endl;
498 DUNE_THROW(InvalidStateException,
"eigenValues: Eigenvalue calculation failed!");
501 DUNE_THROW(NotImplemented,
"LAPACK not found!");
507 template <Jobs Tag,
int dim,
typename K>
508 static void eigenValuesVectorsImpl(
const FieldMatrix<K, dim, dim>& matrix,
509 FieldVector<K, dim>& eigenValues,
510 FieldMatrix<K, dim, dim>& eigenVectors)
512 eigenValuesVectorsLapackImpl<Tag>(matrix,eigenValues,eigenVectors);
523 template <
int dim,
typename K>
528 Impl::eigenValuesVectorsImpl<Impl::Jobs::OnlyEigenvalues>(matrix,
eigenValues, dummy);
539 template <
int dim,
typename K>
544 Impl::eigenValuesVectorsImpl<Impl::Jobs::EigenvaluesEigenvectors>(matrix,
eigenValues, eigenVectors);
554 template <
int dim,
typename K>
559 Impl::eigenValuesVectorsLapackImpl<Impl::Jobs::EigenvaluesEigenvectors>(matrix,
eigenValues, dummy);
570 template <
int dim,
typename K>
575 Impl::eigenValuesVectorsLapackImpl<Impl::Jobs::EigenvaluesEigenvectors>(matrix,
eigenValues, eigenVectors);
585 template <
int dim,
typename K,
class C>
591 const long int N = dim ;
592 const char jobvl =
'n';
593 const char jobvr =
'n';
595 constexpr bool isKLapackType = std::is_same_v<K,double> || std::is_same_v<K,float>;
596 using LapackNumType = std::conditional_t<isKLapackType, K, double>;
599 LapackNumType matrixVector[dim * dim];
603 for(
int i=0; i<dim; ++i)
605 for(
int j=0; j<dim; ++j, ++row)
607 matrixVector[ row ] = matrix[ i ][ j ];
612 LapackNumType eigenR[dim];
613 LapackNumType eigenI[dim];
614 LapackNumType work[3*dim];
618 const long int lwork = 3*dim;
621 eigenValuesNonsymLapackCall(&jobvl, &jobvr, &N, &matrixVector[0], &N,
622 &eigenR[0], &eigenI[0],
nullptr, &N,
nullptr, &N, &work[0],
627 std::cerr <<
"For matrix " << matrix <<
" eigenvalue calculation failed! " << std::endl;
630 for (
int i=0; i<N; ++i) {
A dense n x m matrix.
Definition: fmatrix.hh:117
vector space out of a tensor product of fields.
Definition: fvector.hh:91
Default exception if a function was called while the object is not in a valid state for that function...
Definition: exceptions.hh:281
Default exception for dummy implementations.
Definition: exceptions.hh:263
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:586
static void eigenValues(const FieldMatrix< K, dim, dim > &matrix, FieldVector< K, dim > &eigenValues)
calculates the eigenvalues of a symmetric field matrix
Definition: fmatrixev.hh:524
static void eigenValuesLapack(const FieldMatrix< K, dim, dim > &matrix, FieldVector< K, dim > &eigenValues)
calculates the eigenvalues of a symmetric field matrix
Definition: fmatrixev.hh:555
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:540
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:571
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:218
constexpr auto max
Function object that returns the greater of the given values.
Definition: hybridutilities.hh:484
Some useful basic math stuff.
Dune namespace.
Definition: alignedallocator.hh:13
static const Field pi()
Archimedes' constant.
Definition: math.hh:48