dune-common  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 
16 #include <dune/common/fvector.hh>
17 #include <dune/common/fmatrix.hh>
18 #include <dune/common/math.hh>
19 
20 namespace 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  {
374  std::swap(eigenVectors[0], eigenVectors[1]);
375  }
376  if (eigenValues[1] > eigenValues[2])
377  {
379  std::swap(eigenVectors[1], eigenVectors[2]);
380  }
381  if (eigenValues[0] > eigenValues[1])
382  {
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,
520  {
521  Impl::EVDummy<K,dim> dummy;
522  Impl::eigenValuesVectorsImpl<Impl::Jobs::OnlyEigenvalues>(matrix, eigenValues, dummy);
523  }
524 
533  template <int dim, typename K>
534  static void eigenValuesVectors(const FieldMatrix<K, dim, dim>& matrix,
536  FieldMatrix<K, dim, dim>& eigenVectors)
537  {
538  Impl::eigenValuesVectorsImpl<Impl::Jobs::EigenvaluesEigenvectors>(matrix, eigenValues, eigenVectors);
539  }
540 
548  template <int dim, typename K>
549  static void eigenValuesLapack(const FieldMatrix<K, dim, dim>& matrix,
551  {
552  Impl::EVDummy<K,dim> dummy;
553  Impl::eigenValuesVectorsLapackImpl<Impl::Jobs::EigenvaluesEigenvectors>(matrix, eigenValues, dummy);
554  }
555 
564  template <int dim, typename K>
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>
580  static void eigenValuesNonSym(const FieldMatrix<K, dim, dim>& matrix,
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 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:11
void swap(T &v1, T &v2, bool mask)
Definition: simd.hh:470
auto max(const AlignedNumber< T, align > &a, const AlignedNumber< T, align > &b)
Definition: debugalign.hh:412
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
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 for dummy implementations.
Definition: exceptions.hh:261
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