DUNE-FEM (unstable)

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// SPDX-FileCopyrightInfo: Copyright © DUNE Project contributors, see file LICENSE.md in module root
4// SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception
5#ifndef DUNE_FMATRIXEIGENVALUES_HH
6#define DUNE_FMATRIXEIGENVALUES_HH
7
12#include <algorithm>
13#include <iostream>
14#include <cmath>
15#include <cassert>
16
17#include <dune-common-config.hh> // HAVE_LAPACK
21#include <dune/common/math.hh>
22
23namespace Dune {
24
30 namespace FMatrixHelp {
31
32#if HAVE_LAPACK
33 // defined in fmatrixev.cc
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);
38
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);
44
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);
49
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);
55
56#endif
57
58 namespace Impl {
59 //internal tag to activate/disable code for eigenvector calculation at compile time
60 enum Jobs { OnlyEigenvalues=0, EigenvaluesEigenvectors=1 };
61
62 //internal dummy used if only eigenvalues are to be calculated
63 template<typename K, int dim>
64 using EVDummy = FieldMatrix<K, dim, dim>;
65
66 //compute the cross-product of two vectors
67 template<typename K>
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]};
70 }
71
72 template <typename K>
73 static void eigenValues2dImpl(const FieldMatrix<K, 2, 2>& matrix,
74 FieldVector<K, 2>& eigenvalues)
75 {
76 using std::sqrt;
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;
81 if (q < 0)
82 {
83 std::cout << matrix << std::endl;
84 // Complex eigenvalues are either caused by non-symmetric matrices or by round-off errors
85 DUNE_THROW(MathError, "Complex eigenvalue detected (which this implementation cannot handle).");
86 }
87
88 // get square root
89 q = sqrt(q);
90
91 // store eigenvalues in ascending order
92 eigenvalues[0] = p - q;
93 eigenvalues[1] = p + q;
94 }
95
96 /*
97 This implementation was adapted from the pseudo-code (Python?) implementation found on
98 http://en.wikipedia.org/wiki/Eigenvalue_algorithm (retrieved late August 2014).
99 Wikipedia claims to have taken it from
100 Smith, Oliver K. (April 1961), Eigenvalues of a symmetric 3 × 3 matrix.,
101 Communications of the ACM 4 (4): 168, doi:10.1145/355578.366316
102 */
103 template <typename K>
104 static K eigenValues3dImpl(const FieldMatrix<K, 3, 3>& matrix,
105 FieldVector<K, 3>& eigenvalues)
106 {
107 using std::sqrt;
108 using std::acos;
109 using real_type = typename FieldTraits<K>::real_type;
110 const K pi = MathematicalConstants<K>::pi();
111 K p1 = matrix[0][1]*matrix[0][1] + matrix[0][2]*matrix[0][2] + matrix[1][2]*matrix[1][2];
112
113 if (p1 <= std::numeric_limits<K>::epsilon()) {
114 // A is diagonal.
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());
119
120 return 0.0;
121 }
122 else
123 {
124 // q = trace(A)/3
125 K q = 0;
126 for (int i=0; i<3; i++)
127 q += matrix[i][i] / 3.0;
128
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;
130 K p = sqrt(p2 / 6);
131 // B = (1 / p) * (A - q * I); // I is the identity matrix
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));
136
137 K r = B.determinant() / 2.0;
138
139 /*In exact arithmetic for a symmetric matrix -1 <= r <= 1
140 but computation error can leave it slightly outside this range.
141 acos(z) function requires |z| <= 1, but will fail silently
142 and return NaN if the input is larger than 1 in magnitude.
143 Thus r is clamped to [-1,1].*/
144 using std::clamp;
145 r = clamp<K>(r, -1.0, 1.0);
146 K phi = acos(r) / 3.0;
147
148 // the eigenvalues satisfy eig[2] <= eig[1] <= eig[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]; // since trace(matrix) = eig1 + eig2 + eig3
152
153 return r;
154 }
155 }
156
157 //see https://www.geometrictools.com/Documentation/RobustEigenSymmetric3x3.pdf
158 //Robustly compute a right-handed orthonormal set {u, v, evec0}.
159 template<typename K>
160 void orthoComp(const FieldVector<K,3>& evec0, FieldVector<K,3>& u, FieldVector<K,3>& v) {
161 using std::abs;
162 if(abs(evec0[0]) > abs(evec0[1])) {
163 //The component of maximum absolute value is either evec0[0] or evec0[2].
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]});
167 }
168 else {
169 //The component of maximum absolute value is either evec0[1] or evec0[2].
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]});
173 }
174 v = crossProduct(evec0, u);
175 }
176
177 //see https://www.geometrictools.com/Documentation/RobustEigenSymmetric3x3.pdf
178 template<typename K>
179 void eig0(const FieldMatrix<K,3,3>& matrix, K eval0, FieldVector<K,3>& evec0) {
180 /* Compute a unit-length eigenvector for eigenvalue[i0]. The
181 matrix is rank 2, so two of the rows are linearly independent.
182 For a robust computation of the eigenvector, select the two
183 rows whose cross product has largest length of all pairs of
184 rows. */
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};
189
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();
196
197 auto dmax = d0 ;
198 int imax = 0;
199 if(d1>dmax) {
200 dmax = d1;
201 imax = 1;
202 }
203 if(d2>dmax)
204 imax = 2;
205
206 if(imax == 0)
207 evec0 = r0xr1 / d0;
208 else if(imax == 1)
209 evec0 = r0xr2 / d1;
210 else
211 evec0 = r1xr2 / d2;
212 }
213
214 //see https://www.geometrictools.com/Documentation/RobustEigenSymmetric3x3.pdf
215 template<typename K>
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>;
218
219 //Robustly compute a right-handed orthonormal set {u, v, evec0}.
220 Vector u,v;
221 orthoComp(evec0, u, v);
222
223 /* Let e be eval1 and let E be a corresponding eigenvector which
224 is a solution to the linear system (A - e*I)*E = 0. The matrix
225 (A - e*I) is 3x3, not invertible (so infinitely many
226 solutions), and has rank 2 when eval1 and eval are different.
227 It has rank 1 when eval1 and eval2 are equal. Numerically, it
228 is difficult to compute robustly the rank of a matrix. Instead,
229 the 3x3 linear system is reduced to a 2x2 system as follows.
230 Define the 3x2 matrix J = [u,v] whose columns are the u and v
231 computed previously. Define the 2x1 vector X = J*E. The 2x2
232 system is 0 = M * X = (J^T * (A - e*I) * J) * X where J^T is
233 the transpose of J and M = J^T * (A - e*I) * J is a 2x2 matrix.
234 The system may be written as
235 +- -++- -+ +- -+
236 | U^T*A*U - e U^T*A*V || x0 | = e * | x0 |
237 | V^T*A*U V^T*A*V - e || x1 | | x1 |
238 +- -++ -+ +- -+
239 where X has row entries x0 and x1. */
240
241 Vector Au, Av;
242 matrix.mv(u, Au);
243 matrix.mv(v, Av);
244
245 auto m00 = u.dot(Au) - eval1;
246 auto m01 = u.dot(Av);
247 auto m11 = v.dot(Av) - eval1;
248
249 /* For robustness, choose the largest-length row of M to compute
250 the eigenvector. The 2-tuple of coefficients of U and V in the
251 assignments to eigenvector[1] lies on a circle, and U and V are
252 unit length and perpendicular, so eigenvector[1] is unit length
253 (within numerical tolerance). */
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) {
262 m01 /= m00;
263 m00 = 1.0 / sqrt(1.0 + m01*m01);
264 m01 *= m00;
265 }
266 else {
267 m00 /= m01;
268 m01 = 1.0 / sqrt(1.0 + m00*m00);
269 m00 *= m01;
270 }
271 evec1 = m01*u - m00*v;
272 }
273 else
274 evec1 = u;
275 }
276 else {
277 auto maxAbsComp = max(absM11, absM01);
278 if(maxAbsComp > 0.0) {
279 if(absM11 >= absM01) {
280 m01 /= m11;
281 m11 = 1.0 / sqrt(1.0 + m01*m01);
282 m01 *= m11;
283 }
284 else {
285 m11 /= m01;
286 m01 = 1.0 / sqrt(1.0 + m11*m11);
287 m11 *= m01;
288 }
289 evec1 = m11*u - m01*v;
290 }
291 else
292 evec1 = u;
293 }
294 }
295
296 // 1d specialization
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)
301 {
302 eigenValues[0] = matrix[0][0];
303 if constexpr(Tag==EigenvaluesEigenvectors)
304 eigenVectors[0] = {1.0};
305 }
306
307
308 // 2d specialization
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)
313 {
314 // Compute eigen values
315 Impl::eigenValues2dImpl(matrix, eigenValues);
316
317 // Compute eigenvectors by exploiting the Cayley–Hamilton theorem.
318 // If λ_1, λ_2 are the eigenvalues, then (A - λ_1I )(A - λ_2I ) = (A - λ_2I )(A - λ_1I ) = 0,
319 // so the columns of (A - λ_2I ) are annihilated by (A - λ_1I ) and vice versa.
320 // Assuming neither matrix is zero, the columns of each must include eigenvectors
321 // for the other eigenvalue. (If either matrix is zero, then A is a multiple of the
322 // identity and any non-zero vector is an eigenvector.)
323 // From: https://en.wikipedia.org/wiki/Eigenvalue_algorithm#2x2_matrices
324 if constexpr(Tag==EigenvaluesEigenvectors) {
325
326 // Special casing for multiples of the identity
327 FieldMatrix<K,2,2> temp = matrix;
328 temp[0][0] -= eigenValues[0];
329 temp[1][1] -= eigenValues[0];
330 if(temp.infinity_norm() <= 1e-14) {
331 eigenVectors[0] = {1.0, 0.0};
332 eigenVectors[1] = {0.0, 1.0};
333 }
334 else {
335 // The columns of A - λ_2I are eigenvectors for λ_1, or zero.
336 // Take the column with the larger norm to avoid zero columns.
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();
340
341 // The columns of A - λ_1I are eigenvectors for λ_2, or zero.
342 // Take the column with the larger norm to avoid zero columns.
343 ev0 = {matrix[0][0]-eigenValues[0], matrix[1][0]};
344 ev1 = {matrix[0][1], matrix[1][1]-eigenValues[0]};
345 eigenVectors[1] = (ev0.two_norm2() >= ev1.two_norm2()) ? ev0/ev0.two_norm() : ev1/ev1.two_norm();
346 }
347 }
348 }
349
350 // 3d specialization
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)
355 {
356 using Vector = FieldVector<K,3>;
357 using Matrix = FieldMatrix<K,3,3>;
358
359 //compute eigenvalues
360 /* Precondition the matrix by factoring out the maximum absolute
361 value of the components. This guards against floating-point
362 overflow when computing the eigenvalues.*/
363 using std::isnormal;
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);
367
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())
371 {
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}};
374
375 // Use bubble sort to jointly sort eigenvalues and eigenvectors
376 // such that eigenvalues are ascending
377 if (eigenValues[0] > eigenValues[1])
378 {
379 std::swap(eigenValues[0], eigenValues[1]);
380 std::swap(eigenVectors[0], eigenVectors[1]);
381 }
382 if (eigenValues[1] > eigenValues[2])
383 {
384 std::swap(eigenValues[1], eigenValues[2]);
385 std::swap(eigenVectors[1], eigenVectors[2]);
386 }
387 if (eigenValues[0] > eigenValues[1])
388 {
389 std::swap(eigenValues[0], eigenValues[1]);
390 std::swap(eigenVectors[0], eigenVectors[1]);
391 }
392 }
393 else {
394 /*Compute the eigenvectors so that the set
395 [evec[0], evec[1], evec[2]] is right handed and
396 orthonormal. */
397
398 Matrix evec(0.0);
399 Vector eval(eigenValues);
400 if(r >= 0) {
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]);
404 }
405 else {
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]);
409 }
410 //sort eval/evec-pairs in ascending order
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){
418 eigenValues[i] = pairs[i].first;
419 eigenVectors[i] = pairs[i].second;
420 }
421 }
422 }
423 //The preconditioning scaled the matrix, which scales the eigenvalues. Revert the scaling.
424 eigenValues *= maxAbsElement;
425 }
426
427 // forwarding to LAPACK with corresponding tag
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)
432 {
433 {
434#if HAVE_LAPACK
435 /*Lapack uses a proprietary tag to determine whether both eigenvalues and
436 -vectors ('v') or only eigenvalues ('n') should be calculated */
437 const char jobz = "nv"[Tag];
438
439 const long int N = dim ;
440 const char uplo = 'u'; // use upper triangular matrix
441
442 // length of matrix vector, LWORK >= max(1,3*N-1)
443 const long int lwork = 3*N -1 ;
444
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>;
447
448 // matrix to put into dsyev
449 LapackNumType matrixVector[dim * dim];
450
451 // copy matrix
452 int row = 0;
453 for(int i=0; i<dim; ++i)
454 {
455 for(int j=0; j<dim; ++j, ++row)
456 {
457 matrixVector[ row ] = matrix[ i ][ j ];
458 }
459 }
460
461 // working memory
462 LapackNumType workSpace[lwork];
463
464 // return value information
465 long int info = 0;
466 LapackNumType* ev;
467 if constexpr (isKLapackType){
468 ev = &eigenValues[0];
469 }else{
470 ev = new LapackNumType[dim];
471 }
472
473 // call LAPACK routine (see fmatrixev.cc)
474 eigenValuesLapackCall(&jobz, &uplo, &N, &matrixVector[0], &N,
475 ev, &workSpace[0], &lwork, &info);
476
477 if constexpr (!isKLapackType){
478 for(size_t i=0;i<dim;++i)
479 eigenValues[i] = ev[i];
480 delete[] ev;
481 }
482
483 // restore eigenvectors matrix
484 if (Tag==Jobs::EigenvaluesEigenvectors){
485 row = 0;
486 for(int i=0; i<dim; ++i)
487 {
488 for(int j=0; j<dim; ++j, ++row)
489 {
490 eigenVectors[ i ][ j ] = matrixVector[ row ];
491 }
492 }
493 }
494
495 if( info != 0 )
496 {
497 std::cerr << "For matrix " << matrix << " eigenvalue calculation failed! " << std::endl;
498 DUNE_THROW(InvalidStateException,"eigenValues: Eigenvalue calculation failed!");
499 }
500#else
501 DUNE_THROW(NotImplemented,"LAPACK not found!");
502#endif
503 }
504 }
505
506 // generic specialization
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)
511 {
512 eigenValuesVectorsLapackImpl<Tag>(matrix,eigenValues,eigenVectors);
513 }
514 } //namespace Impl
515
523 template <int dim, typename K>
524 static void eigenValues(const FieldMatrix<K, dim, dim>& matrix,
525 FieldVector<K ,dim>& eigenValues)
526 {
528 Impl::eigenValuesVectorsImpl<Impl::Jobs::OnlyEigenvalues>(matrix, eigenValues, dummy);
529 }
530
539 template <int dim, typename K>
541 FieldVector<K ,dim>& eigenValues,
542 FieldMatrix<K, dim, dim>& eigenVectors)
543 {
544 Impl::eigenValuesVectorsImpl<Impl::Jobs::EigenvaluesEigenvectors>(matrix, eigenValues, eigenVectors);
545 }
546
554 template <int dim, typename K>
556 FieldVector<K, dim>& eigenValues)
557 {
559 Impl::eigenValuesVectorsLapackImpl<Impl::Jobs::EigenvaluesEigenvectors>(matrix, eigenValues, dummy);
560 }
561
570 template <int dim, typename K>
572 FieldVector<K, dim>& eigenValues,
573 FieldMatrix<K, dim, dim>& eigenVectors)
574 {
575 Impl::eigenValuesVectorsLapackImpl<Impl::Jobs::EigenvaluesEigenvectors>(matrix, eigenValues, eigenVectors);
576 }
577
585 template <int dim, typename K, class C>
587 FieldVector<C, dim>& eigenValues)
588 {
589#if HAVE_LAPACK
590 {
591 const long int N = dim ;
592 const char jobvl = 'n';
593 const char jobvr = 'n';
594
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>;
597
598 // matrix to put into dgeev
599 LapackNumType matrixVector[dim * dim];
600
601 // copy matrix
602 int row = 0;
603 for(int i=0; i<dim; ++i)
604 {
605 for(int j=0; j<dim; ++j, ++row)
606 {
607 matrixVector[ row ] = matrix[ i ][ j ];
608 }
609 }
610
611 // working memory
612 LapackNumType eigenR[dim];
613 LapackNumType eigenI[dim];
614 LapackNumType work[3*dim];
615
616 // return value information
617 long int info = 0;
618 const long int lwork = 3*dim;
619
620 // call LAPACK routine (see fmatrixev_ext.cc)
621 eigenValuesNonsymLapackCall(&jobvl, &jobvr, &N, &matrixVector[0], &N,
622 &eigenR[0], &eigenI[0], nullptr, &N, nullptr, &N, &work[0],
623 &lwork, &info);
624
625 if( info != 0 )
626 {
627 std::cerr << "For matrix " << matrix << " eigenvalue calculation failed! " << std::endl;
628 DUNE_THROW(InvalidStateException,"eigenValues: Eigenvalue calculation failed!");
629 }
630 for (int i=0; i<N; ++i) {
631 eigenValues[i].real = eigenR[i];
632 eigenValues[i].imag = eigenI[i];
633 }
634 }
635#else
636 DUNE_THROW(NotImplemented,"LAPACK not found!");
637#endif
638 }
639 } // end namespace FMatrixHelp
640
643} // end namespace Dune
644#endif
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
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 21, 23:30, 2024)