Dune Core Modules (2.5.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 <iostream>
11#include <cmath>
12#include <cassert>
13
17
18namespace Dune {
19
25 namespace FMatrixHelp {
26
27 // defined in fmatrixev.cc
28 extern void eigenValuesLapackCall(
29 const char* jobz, const char* uplo, const long
30 int* n, double* a, const long int* lda, double* w,
31 double* work, const long int* lwork, long int* info);
32
33 extern void eigenValuesNonsymLapackCall(
34 const char* jobvl, const char* jobvr, const long
35 int* n, double* a, const long int* lda, double* wr, double* wi, double* vl,
36 const long int* ldvl, double* vr, const long int* ldvr, double* work,
37 const long int* lwork, const long int* info);
38
44 template <typename K>
45 static void eigenValues(const FieldMatrix<K, 1, 1>& matrix,
46 FieldVector<K, 1>& eigenvalues)
47 {
48 eigenvalues[0] = matrix[0][0];
49 }
50
56 template <typename K>
57 static void eigenValues(const FieldMatrix<K, 2, 2>& matrix,
58 FieldVector<K, 2>& eigenvalues)
59 {
60 const K detM = matrix[0][0] * matrix[1][1] - matrix[1][0] * matrix[0][1];
61 const K p = 0.5 * (matrix[0][0] + matrix [1][1]);
62 K q = p * p - detM;
63 if( q < 0 && q > -1e-14 ) q = 0;
64 if (p < 0 || q < 0)
65 {
66 std::cout << p << " p | q " << q << "\n";
67 std::cout << matrix << std::endl;
68 std::cout << "something went wrong in Eigenvalues for matrix!" << std::endl;
69 assert(false);
70 abort();
71 }
72
73 // get square root
74 q = std :: sqrt(q);
75
76 // store eigenvalues in ascending order
77 eigenvalues[0] = p - q;
78 eigenvalues[1] = p + q;
79 }
80
93 template <typename K>
94 static void eigenValues(const FieldMatrix<K, 3, 3>& matrix,
95 FieldVector<K, 3>& eigenvalues)
96 {
97 K p1 = matrix[0][1]*matrix[0][1] + matrix[0][2]*matrix[0][2] + matrix[1][2]*matrix[1][2];
98
99 if (p1 <= 1e-8)
100 {
101 // A is diagonal.
102 eigenvalues[0] = matrix[0][0];
103 eigenvalues[1] = matrix[1][1];
104 eigenvalues[2] = matrix[2][2];
105 }
106 else
107 {
108 // q = trace(A)/3
109 K q = 0;
110 for (int i=0; i<3; i++)
111 q += matrix[i][i]/3.0;
112
113 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 * p1;
114 K p = std::sqrt(p2 / 6);
115 // B = (1 / p) * (A - q * I); // I is the identity matrix
117 for (int i=0; i<3; i++)
118 for (int j=0; j<3; j++)
119 B[i][j] = (1/p) * (matrix[i][j] - q*(i==j));
120
121 K r = B.determinant() / 2.0;
122
123 // In exact arithmetic for a symmetric matrix -1 <= r <= 1
124 // but computation error can leave it slightly outside this range.
125 K phi;
126 if (r <= -1)
127 phi = M_PI / 3.0;
128 else if (r >= 1)
129 phi = 0;
130 else
131 phi = std::acos(r) / 3;
132
133 // the eigenvalues satisfy eig[2] <= eig[1] <= eig[0]
134 eigenvalues[2] = q + 2 * p * cos(phi);
135 eigenvalues[0] = q + 2 * p * cos(phi + (2*M_PI/3));
136 eigenvalues[1] = 3 * q - eigenvalues[0] - eigenvalues[2]; // since trace(matrix) = eig1 + eig2 + eig3
137 }
138 }
139
147 template <int dim, typename K>
148 static void eigenValues(const FieldMatrix<K, dim, dim>& matrix,
149 FieldVector<K, dim>& eigenvalues)
150 {
151 {
152 const long int N = dim ;
153 const char jobz = 'n'; // only calculate eigenvalues
154 const char uplo = 'u'; // use upper triangular matrix
155
156 // length of matrix vector
157 const long int w = N * N ;
158
159 // matrix to put into dsyev
160 double matrixVector[dim * dim];
161
162 // copy matrix
163 int row = 0;
164 for(int i=0; i<dim; ++i)
165 {
166 for(int j=0; j<dim; ++j, ++row)
167 {
168 matrixVector[ row ] = matrix[ i ][ j ];
169 }
170 }
171
172 // working memory
173 double workSpace[dim * dim];
174
175 // return value information
176 long int info = 0;
177
178 // call LAPACK routine (see fmatrixev.cc)
179 eigenValuesLapackCall(&jobz, &uplo, &N, &matrixVector[0], &N,
180 &eigenvalues[0], &workSpace[0], &w, &info);
181
182 if( info != 0 )
183 {
184 std::cerr << "For matrix " << matrix << " eigenvalue calculation failed! " << std::endl;
185 DUNE_THROW(InvalidStateException,"eigenValues: Eigenvalue calculation failed!");
186 }
187 }
188 }
196 template <int dim, typename K, class C>
198 FieldVector<C, dim>& eigenValues)
199 {
200 {
201 const long int N = dim ;
202 const char jobvl = 'n';
203 const char jobvr = 'n';
204
205 // matrix to put into dgeev
206 double matrixVector[dim * dim];
207
208 // copy matrix
209 int row = 0;
210 for(int i=0; i<dim; ++i)
211 {
212 for(int j=0; j<dim; ++j, ++row)
213 {
214 matrixVector[ row ] = matrix[ i ][ j ];
215 }
216 }
217
218 // working memory
219 double eigenR[dim];
220 double eigenI[dim];
221 double work[3*dim];
222
223 // return value information
224 long int info = 0;
225 long int lwork = 3*dim;
226
227 // call LAPACK routine (see fmatrixev_ext.cc)
228 eigenValuesNonsymLapackCall(&jobvl, &jobvr, &N, &matrixVector[0], &N,
229 &eigenR[0], &eigenI[0], 0, &N, 0, &N, &work[0],
230 &lwork, &info);
231
232 if( info != 0 )
233 {
234 std::cerr << "For matrix " << matrix << " eigenvalue calculation failed! " << std::endl;
235 DUNE_THROW(InvalidStateException,"eigenValues: Eigenvalue calculation failed!");
236 }
237 for (int i=0; i<N; ++i) {
238 eigenValues[i].real = eigenR[i];
239 eigenValues[i].imag = eigenI[i];
240 }
241 }
242
243 }
244
245 } // end namespace FMatrixHelp
246
249} // end namespace Dune
250#endif
field_type determinant() const
calculates the determinant of this matrix
A dense n x m matrix.
Definition: fmatrix.hh:68
vector space out of a tensor product of fields.
Definition: fvector.hh:93
Default exception if a function was called while the object is not in a valid state for that function...
Definition: exceptions.hh:279
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 symmetric field matrix
Definition: fmatrixev.hh:197
static void eigenValues(const FieldMatrix< K, 1, 1 > &matrix, FieldVector< K, 1 > &eigenvalues)
calculates the eigenvalues of a symmetric field matrix
Definition: fmatrixev.hh:45
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
Dune namespace.
Definition: alignment.hh:11
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 23, 23:29, 2024)