Dune Core Modules (2.8.0)

fmatrixev.hh
Go to the documentation of this file.
1// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
2// vi: set et ts=4 sw=2 sts=2:
3#ifndef DUNE_FMATRIXEIGENVALUES_HH
4#define DUNE_FMATRIXEIGENVALUES_HH
5
10#include <algorithm>
11#include <iostream>
12#include <cmath>
13#include <cassert>
14
18#include <dune/common/math.hh>
19
20namespace Dune {
21
27 namespace FMatrixHelp {
28
29#if HAVE_LAPACK
30 // defined in fmatrixev.cc
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);
35
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);
41
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);
46
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);
52
53#endif
54
55 namespace Impl {
56 //internal tag to activate/disable code for eigenvector calculation at compile time
57 enum Jobs { OnlyEigenvalues=0, EigenvaluesEigenvectors=1 };
58
59 //internal dummy used if only eigenvalues are to be calculated
60 template<typename K, int dim>
61 using EVDummy = FieldMatrix<K, dim, dim>;
62
63 //compute the cross-product of two vectors
64 template<typename K>
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]};
67 }
68
69 template <typename K>
70 static void eigenValues2dImpl(const FieldMatrix<K, 2, 2>& matrix,
71 FieldVector<K, 2>& eigenvalues)
72 {
73 using std::sqrt;
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;
78 if (q < 0)
79 {
80 std::cout << matrix << std::endl;
81 // Complex eigenvalues are either caused by non-symmetric matrices or by round-off errors
82 DUNE_THROW(MathError, "Complex eigenvalue detected (which this implementation cannot handle).");
83 }
84
85 // get square root
86 q = sqrt(q);
87
88 // store eigenvalues in ascending order
89 eigenvalues[0] = p - q;
90 eigenvalues[1] = p + q;
91 }
92
93 /*
94 This implementation was adapted from the pseudo-code (Python?) implementation found on
95 http://en.wikipedia.org/wiki/Eigenvalue_algorithm (retrieved late August 2014).
96 Wikipedia claims to have taken it from
97 Smith, Oliver K. (April 1961), Eigenvalues of a symmetric 3 × 3 matrix.,
98 Communications of the ACM 4 (4): 168, doi:10.1145/355578.366316
99 */
100 template <typename K>
101 static K eigenValues3dImpl(const FieldMatrix<K, 3, 3>& matrix,
102 FieldVector<K, 3>& eigenvalues)
103 {
104 using std::sqrt;
105 using std::acos;
106 using real_type = typename FieldTraits<K>::real_type;
107 const K pi = MathematicalConstants<K>::pi();
108 K p1 = matrix[0][1]*matrix[0][1] + matrix[0][2]*matrix[0][2] + matrix[1][2]*matrix[1][2];
109
110 if (p1 <= std::numeric_limits<K>::epsilon()) {
111 // A is diagonal.
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());
116
117 return 0.0;
118 }
119 else
120 {
121 // q = trace(A)/3
122 K q = 0;
123 for (int i=0; i<3; i++)
124 q += matrix[i][i] / 3.0;
125
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;
127 K p = sqrt(p2 / 6);
128 // B = (1 / p) * (A - q * I); // I is the identity matrix
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));
133
134 K r = B.determinant() / 2.0;
135
136 /*In exact arithmetic for a symmetric matrix -1 <= r <= 1
137 but computation error can leave it slightly outside this range.
138 acos(z) function requires |z| <= 1, but will fail silently
139 and return NaN if the input is larger than 1 in magnitude.
140 Thus r is clamped to [-1,1].*/
141 r = std::min<K>(std::max<K>(r, -1.0), 1.0);
142 K phi = acos(r) / 3.0;
143
144 // the eigenvalues satisfy eig[2] <= eig[1] <= eig[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]; // since trace(matrix) = eig1 + eig2 + eig3
148
149 return r;
150 }
151 }
152
153 //see https://www.geometrictools.com/Documentation/RobustEigenSymmetric3x3.pdf
154 //Robustly compute a right-handed orthonormal set {u, v, evec0}.
155 template<typename K>
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])) {
158 //The component of maximum absolute value is either evec0[0] or evec0[2].
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]});
162 }
163 else {
164 //The component of maximum absolute value is either evec0[1] or evec0[2].
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]});
168 }
169 v = crossProduct(evec0, u);
170 }
171
172 //see https://www.geometrictools.com/Documentation/RobustEigenSymmetric3x3.pdf
173 template<typename K>
174 void eig0(const FieldMatrix<K,3,3>& matrix, K eval0, FieldVector<K,3>& evec0) {
175 /* Compute a unit-length eigenvector for eigenvalue[i0]. The
176 matrix is rank 2, so two of the rows are linearly independent.
177 For a robust computation of the eigenvector, select the two
178 rows whose cross product has largest length of all pairs of
179 rows. */
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};
184
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();
191
192 auto dmax = d0 ;
193 int imax = 0;
194 if(d1>dmax) {
195 dmax = d1;
196 imax = 1;
197 }
198 if(d2>dmax)
199 imax = 2;
200
201 if(imax == 0)
202 evec0 = r0xr1 / d0;
203 else if(imax == 1)
204 evec0 = r0xr2 / d1;
205 else
206 evec0 = r1xr2 / d2;
207 }
208
209 //see https://www.geometrictools.com/Documentation/RobustEigenSymmetric3x3.pdf
210 template<typename K>
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>;
213
214 //Robustly compute a right-handed orthonormal set {u, v, evec0}.
215 Vector u,v;
216 orthoComp(evec0, u, v);
217
218 /* Let e be eval1 and let E be a corresponding eigenvector which
219 is a solution to the linear system (A - e*I)*E = 0. The matrix
220 (A - e*I) is 3x3, not invertible (so infinitely many
221 solutions), and has rank 2 when eval1 and eval are different.
222 It has rank 1 when eval1 and eval2 are equal. Numerically, it
223 is difficult to compute robustly the rank of a matrix. Instead,
224 the 3x3 linear system is reduced to a 2x2 system as follows.
225 Define the 3x2 matrix J = [u,v] whose columns are the u and v
226 computed previously. Define the 2x1 vector X = J*E. The 2x2
227 system is 0 = M * X = (J^T * (A - e*I) * J) * X where J^T is
228 the transpose of J and M = J^T * (A - e*I) * J is a 2x2 matrix.
229 The system may be written as
230 +- -++- -+ +- -+
231 | U^T*A*U - e U^T*A*V || x0 | = e * | x0 |
232 | V^T*A*U V^T*A*V - e || x1 | | x1 |
233 +- -++ -+ +- -+
234 where X has row entries x0 and x1. */
235
236 Vector Au, Av;
237 matrix.mv(u, Au);
238 matrix.mv(v, Av);
239
240 auto m00 = u.dot(Au) - eval1;
241 auto m01 = u.dot(Av);
242 auto m11 = v.dot(Av) - eval1;
243
244 /* For robustness, choose the largest-length row of M to compute
245 the eigenvector. The 2-tuple of coefficients of U and V in the
246 assignments to eigenvector[1] lies on a circle, and U and V are
247 unit length and perpendicular, so eigenvector[1] is unit length
248 (within numerical tolerance). */
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) {
256 m01 /= m00;
257 m00 = 1.0 / std::sqrt(1.0 + m01*m01);
258 m01 *= m00;
259 }
260 else {
261 m00 /= m01;
262 m01 = 1.0 / std::sqrt(1.0 + m00*m00);
263 m00 *= m01;
264 }
265 evec1 = m01*u - m00*v;
266 }
267 else
268 evec1 = u;
269 }
270 else {
271 auto maxAbsComp = std::max(absM11, absM01);
272 if(maxAbsComp > 0.0) {
273 if(absM11 >= absM01) {
274 m01 /= m11;
275 m11 = 1.0 / std::sqrt(1.0 + m01*m01);
276 m01 *= m11;
277 }
278 else {
279 m11 /= m01;
280 m01 = 1.0 / std::sqrt(1.0 + m11*m11);
281 m11 *= m01;
282 }
283 evec1 = m11*u - m01*v;
284 }
285 else
286 evec1 = u;
287 }
288 }
289
290 // 1d specialization
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)
295 {
296 eigenValues[0] = matrix[0][0];
297 if constexpr(Tag==EigenvaluesEigenvectors)
298 eigenVectors[0] = {1.0};
299 }
300
301
302 // 2d specialization
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)
307 {
308 // Compute eigen values
309 Impl::eigenValues2dImpl(matrix, eigenValues);
310
311 // Compute eigenvectors by exploiting the Cayley–Hamilton theorem.
312 // If λ_1, λ_2 are the eigenvalues, then (A - λ_1I )(A - λ_2I ) = (A - λ_2I )(A - λ_1I ) = 0,
313 // so the columns of (A - λ_2I ) are annihilated by (A - λ_1I ) and vice versa.
314 // Assuming neither matrix is zero, the columns of each must include eigenvectors
315 // for the other eigenvalue. (If either matrix is zero, then A is a multiple of the
316 // identity and any non-zero vector is an eigenvector.)
317 // From: https://en.wikipedia.org/wiki/Eigenvalue_algorithm#2x2_matrices
318 if constexpr(Tag==EigenvaluesEigenvectors) {
319
320 // Special casing for multiples of the identity
321 FieldMatrix<K,2,2> temp = matrix;
322 temp[0][0] -= eigenValues[0];
323 temp[1][1] -= eigenValues[0];
324 if(temp.infinity_norm() <= 1e-14) {
325 eigenVectors[0] = {1.0, 0.0};
326 eigenVectors[1] = {0.0, 1.0};
327 }
328 else {
329 // The columns of A - λ_2I are eigenvectors for λ_1, or zero.
330 // Take the column with the larger norm to avoid zero columns.
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();
334
335 // The columns of A - λ_1I are eigenvectors for λ_2, or zero.
336 // Take the column with the larger norm to avoid zero columns.
337 ev0 = {matrix[0][0]-eigenValues[0], matrix[1][0]};
338 ev1 = {matrix[0][1], matrix[1][1]-eigenValues[0]};
339 eigenVectors[1] = (ev0.two_norm2() >= ev1.two_norm2()) ? ev0/ev0.two_norm() : ev1/ev1.two_norm();
340 }
341 }
342 }
343
344 // 3d specialization
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)
349 {
350 using Vector = FieldVector<K,3>;
351 using Matrix = FieldMatrix<K,3,3>;
352
353 //compute eigenvalues
354 /* Precondition the matrix by factoring out the maximum absolute
355 value of the components. This guards against floating-point
356 overflow when computing the eigenvalues.*/
357 using std::isnormal;
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);
361
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())
365 {
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}};
368
369 // Use bubble sort to jointly sort eigenvalues and eigenvectors
370 // such that eigenvalues are ascending
371 if (eigenValues[0] > eigenValues[1])
372 {
373 std::swap(eigenValues[0], eigenValues[1]);
374 std::swap(eigenVectors[0], eigenVectors[1]);
375 }
376 if (eigenValues[1] > eigenValues[2])
377 {
378 std::swap(eigenValues[1], eigenValues[2]);
379 std::swap(eigenVectors[1], eigenVectors[2]);
380 }
381 if (eigenValues[0] > eigenValues[1])
382 {
383 std::swap(eigenValues[0], eigenValues[1]);
384 std::swap(eigenVectors[0], eigenVectors[1]);
385 }
386 }
387 else {
388 /*Compute the eigenvectors so that the set
389 [evec[0], evec[1], evec[2]] is right handed and
390 orthonormal. */
391
392 Matrix evec(0.0);
393 Vector eval(eigenValues);
394 if(r >= 0) {
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]);
398 }
399 else {
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]);
403 }
404 //sort eval/evec-pairs in ascending order
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){
412 eigenValues[i] = pairs[i].first;
413 eigenVectors[i] = pairs[i].second;
414 }
415 }
416 }
417 //The preconditioning scaled the matrix, which scales the eigenvalues. Revert the scaling.
418 eigenValues *= maxAbsElement;
419 }
420
421 // forwarding to LAPACK with corresponding tag
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)
426 {
427 {
428#if HAVE_LAPACK
429 /*Lapack uses a proprietary tag to determine whether both eigenvalues and
430 -vectors ('v') or only eigenvalues ('n') should be calculated */
431 const char jobz = "nv"[Tag];
432
433 const long int N = dim ;
434 const char uplo = 'u'; // use upper triangular matrix
435
436 // length of matrix vector, LWORK >= max(1,3*N-1)
437 const long int lwork = 3*N -1 ;
438
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>;
441
442 // matrix to put into dsyev
443 LapackNumType matrixVector[dim * dim];
444
445 // copy matrix
446 int row = 0;
447 for(int i=0; i<dim; ++i)
448 {
449 for(int j=0; j<dim; ++j, ++row)
450 {
451 matrixVector[ row ] = matrix[ i ][ j ];
452 }
453 }
454
455 // working memory
456 LapackNumType workSpace[lwork];
457
458 // return value information
459 long int info = 0;
460 LapackNumType* ev;
461 if constexpr (isKLapackType){
462 ev = &eigenValues[0];
463 }else{
464 ev = new LapackNumType[dim];
465 }
466
467 // call LAPACK routine (see fmatrixev.cc)
468 eigenValuesLapackCall(&jobz, &uplo, &N, &matrixVector[0], &N,
469 ev, &workSpace[0], &lwork, &info);
470
471 if constexpr (!isKLapackType){
472 for(size_t i=0;i<dim;++i)
473 eigenValues[i] = ev[i];
474 delete[] ev;
475 }
476
477 // restore eigenvectors matrix
478 if (Tag==Jobs::EigenvaluesEigenvectors){
479 row = 0;
480 for(int i=0; i<dim; ++i)
481 {
482 for(int j=0; j<dim; ++j, ++row)
483 {
484 eigenVectors[ i ][ j ] = matrixVector[ row ];
485 }
486 }
487 }
488
489 if( info != 0 )
490 {
491 std::cerr << "For matrix " << matrix << " eigenvalue calculation failed! " << std::endl;
492 DUNE_THROW(InvalidStateException,"eigenValues: Eigenvalue calculation failed!");
493 }
494#else
495 DUNE_THROW(NotImplemented,"LAPACK not found!");
496#endif
497 }
498 }
499
500 // generic specialization
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)
505 {
506 eigenValuesVectorsLapackImpl<Tag>(matrix,eigenValues,eigenVectors);
507 }
508 } //namespace Impl
509
517 template <int dim, typename K>
518 static void eigenValues(const FieldMatrix<K, dim, dim>& matrix,
519 FieldVector<K ,dim>& eigenValues)
520 {
522 Impl::eigenValuesVectorsImpl<Impl::Jobs::OnlyEigenvalues>(matrix, eigenValues, dummy);
523 }
524
533 template <int dim, typename K>
535 FieldVector<K ,dim>& eigenValues,
536 FieldMatrix<K, dim, dim>& eigenVectors)
537 {
538 Impl::eigenValuesVectorsImpl<Impl::Jobs::EigenvaluesEigenvectors>(matrix, eigenValues, eigenVectors);
539 }
540
548 template <int dim, typename K>
550 FieldVector<K, dim>& eigenValues)
551 {
553 Impl::eigenValuesVectorsLapackImpl<Impl::Jobs::EigenvaluesEigenvectors>(matrix, eigenValues, dummy);
554 }
555
564 template <int dim, typename K>
566 FieldVector<K, dim>& eigenValues,
567 FieldMatrix<K, dim, dim>& eigenVectors)
568 {
569 Impl::eigenValuesVectorsLapackImpl<Impl::Jobs::EigenvaluesEigenvectors>(matrix, eigenValues, eigenVectors);
570 }
571
579 template <int dim, typename K, class C>
581 FieldVector<C, dim>& eigenValues)
582 {
583#if HAVE_LAPACK
584 {
585 const long int N = dim ;
586 const char jobvl = 'n';
587 const char jobvr = 'n';
588
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>;
591
592 // matrix to put into dgeev
593 LapackNumType matrixVector[dim * dim];
594
595 // copy matrix
596 int row = 0;
597 for(int i=0; i<dim; ++i)
598 {
599 for(int j=0; j<dim; ++j, ++row)
600 {
601 matrixVector[ row ] = matrix[ i ][ j ];
602 }
603 }
604
605 // working memory
606 LapackNumType eigenR[dim];
607 LapackNumType eigenI[dim];
608 LapackNumType work[3*dim];
609
610 // return value information
611 long int info = 0;
612 const long int lwork = 3*dim;
613
614 // call LAPACK routine (see fmatrixev_ext.cc)
615 eigenValuesNonsymLapackCall(&jobvl, &jobvr, &N, &matrixVector[0], &N,
616 &eigenR[0], &eigenI[0], nullptr, &N, nullptr, &N, &work[0],
617 &lwork, &info);
618
619 if( info != 0 )
620 {
621 std::cerr << "For matrix " << matrix << " eigenvalue calculation failed! " << std::endl;
622 DUNE_THROW(InvalidStateException,"eigenValues: Eigenvalue calculation failed!");
623 }
624 for (int i=0; i<N; ++i) {
625 eigenValues[i].real = eigenR[i];
626 eigenValues[i].imag = eigenI[i];
627 }
628 }
629#else
630 DUNE_THROW(NotImplemented,"LAPACK not found!");
631#endif
632 }
633 } // end namespace FMatrixHelp
634
637} // end namespace Dune
638#endif
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
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 12, 23:30, 2024)