dune-common  2.2.0
fmatrixev.hh
Go to the documentation of this file.
1 #ifndef DUNE_FMATRIXEIGENVALUES_HH
2 #define DUNE_FMATRIXEIGENVALUES_HH
3 
8 #include <iostream>
9 #include <cmath>
10 #include <cassert>
11 
13 #include <dune/common/fvector.hh>
14 #include <dune/common/fmatrix.hh>
15 
16 namespace Dune {
17 
23 namespace FMatrixHelp {
24 
25 // defined in fmatrixev.cc
26 extern void eigenValuesLapackCall(
27  const char* jobz, const char* uplo, const long
28  int* n, double* a, const long int* lda, double* w,
29  double* work, const long int* lwork, long int* info);
30 
36 template <typename K>
37 static void eigenValues(const FieldMatrix<K, 1, 1>& matrix,
38  FieldVector<K, 1>& eigenvalues)
39 {
40  eigenvalues[0] = matrix[0][0];
41 }
42 
48 template <typename K>
49 static void eigenValues(const FieldMatrix<K, 2, 2>& matrix,
50  FieldVector<K, 2>& eigenvalues)
51 {
52  const K detM = matrix[0][0] * matrix[1][1] - matrix[1][0] * matrix[0][1];
53  const K p = 0.5 * (matrix[0][0] + matrix [1][1]);
54  K q = p * p - detM;
55  if( q < 0 && q > -1e-14 ) q = 0;
56  if (p < 0 || q < 0)
57  {
58  std::cout << p << " p | q " << q << "\n";
59  std::cout << matrix << std::endl;
60  std::cout << "something went wrong in Eigenvalues for matrix!" << std::endl;
61  assert(false);
62  abort();
63  }
64 
65  // get square root
66  q = std :: sqrt(q);
67 
68  // store eigenvalues in ascending order
69  eigenvalues[0] = p - q;
70  eigenvalues[1] = p + q;
71 }
72 
80 template <int dim, typename K>
81 static void eigenValues(const FieldMatrix<K, dim, dim>& matrix,
82  FieldVector<K, dim>& eigenvalues)
83 {
84  {
85  const long int N = dim ;
86  const char jobz = 'n'; // only calculate eigenvalues
87  const char uplo = 'u'; // use upper triangular matrix
88 
89  // length of matrix vector
90  const long int w = N * N ;
91 
92  // matrix to put into dsyev
93  double matrixVector[dim * dim];
94 
95  // copy matrix
96  int row = 0;
97  for(int i=0; i<dim; ++i)
98  {
99  for(int j=0; j<dim; ++j, ++row)
100  {
101  matrixVector[ row ] = matrix[ i ][ j ];
102  }
103  }
104 
105  // working memory
106  double workSpace[dim * dim];
107 
108  // return value information
109  long int info = 0;
110 
111  // call LAPACK routine (see fmatrixev.cc)
112  eigenValuesLapackCall(&jobz, &uplo, &N, &matrixVector[0], &N,
113  &eigenvalues[0], &workSpace[0], &w, &info);
114 
115  if( info != 0 )
116  {
117  std::cerr << "For matrix " << matrix << " eigenvalue calculation failed! " << std::endl;
118  DUNE_THROW(InvalidStateException,"eigenValues: Eigenvalue calculation failed!");
119  }
120  }
121 }
122 
123 } // end namespace FMatrixHelp
124 
127 } // end namespace Dune
128 #endif