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 6 /** \file 7 * \brief Eigenvalue computations for the FieldMatrix class 8 */ 9 10 #include <algorithm> 11 #include <iostream> 12 #include <cmath> 13 #include <cassert> 14 15 #include <dune/common/exceptions.hh> 16 #include <dune/common/fvector.hh> 17 #include <dune/common/fmatrix.hh> 18 #include <dune/common/math.hh> 19 20 namespace Dune { 21 22 /** 23 @addtogroup DenseMatVec 24 @{ 25 */ 26 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> crossProduct(const FieldVector<K,3> & vec0,const FieldVector<K,3> & vec1)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> eigenValues2dImpl(const FieldMatrix<K,2,2> & matrix,FieldVector<K,2> & eigenvalues)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> eigenValues3dImpl(const FieldMatrix<K,3,3> & matrix,FieldVector<K,3> & eigenvalues)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> orthoComp(const FieldVector<K,3> & evec0,FieldVector<K,3> & u,FieldVector<K,3> & v)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> eig0(const FieldMatrix<K,3,3> & matrix,K eval0,FieldVector<K,3> & evec0)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> eig1(const FieldMatrix<K,3,3> & matrix,const FieldVector<K,3> & evec0,FieldVector<K,3> & evec1,K eval1)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> eigenValuesVectorsImpl(const FieldMatrix<K,1,1> & matrix,FieldVector<K,1> & eigenValues,FieldMatrix<K,1,1> & eigenVectors)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> eigenValuesVectorsImpl(const FieldMatrix<K,2,2> & matrix,FieldVector<K,2> & eigenValues,FieldMatrix<K,2,2> & eigenVectors)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> eigenValuesVectorsImpl(const FieldMatrix<K,3,3> & matrix,FieldVector<K,3> & eigenValues,FieldMatrix<K,3,3> & eigenVectors)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> eigenValuesVectorsLapackImpl(const FieldMatrix<K,dim,dim> & matrix,FieldVector<K,dim> & eigenValues,FieldMatrix<K,dim,dim> & eigenVectors)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> eigenValuesVectorsImpl(const FieldMatrix<K,dim,dim> & matrix,FieldVector<K,dim> & eigenValues,FieldMatrix<K,dim,dim> & eigenVectors)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 510 /** \brief calculates the eigenvalues of a symmetric field matrix 511 \param[in] matrix matrix eigenvalues are calculated for 512 \param[out] eigenValues FieldVector that contains eigenvalues in 513 ascending order 514 515 \note specializations for dim=1,2,3 exist, for dim>3 LAPACK::dsyev is used 516 */ 517 template <int dim, typename K> eigenValues(const FieldMatrix<K,dim,dim> & matrix,FieldVector<K,dim> & eigenValues)518 static void eigenValues(const FieldMatrix<K, dim, dim>& matrix, 519 FieldVector<K ,dim>& eigenValues) 520 { 521 Impl::EVDummy<K,dim> dummy; 522 Impl::eigenValuesVectorsImpl<Impl::Jobs::OnlyEigenvalues>(matrix, eigenValues, dummy); 523 } 524 525 /** \brief calculates the eigenvalues and eigenvectors of a symmetric field matrix 526 \param[in] matrix matrix eigenvalues are calculated for 527 \param[out] eigenValues FieldVector that contains eigenvalues in 528 ascending order 529 \param[out] eigenVectors FieldMatrix that contains the eigenvectors 530 531 \note specializations for dim=1,2,3 exist, for dim>3 LAPACK::dsyev is used 532 */ 533 template <int dim, typename K> eigenValuesVectors(const FieldMatrix<K,dim,dim> & matrix,FieldVector<K,dim> & eigenValues,FieldMatrix<K,dim,dim> & eigenVectors)534 static void eigenValuesVectors(const FieldMatrix<K, dim, dim>& matrix, 535 FieldVector<K ,dim>& eigenValues, 536 FieldMatrix<K, dim, dim>& eigenVectors) 537 { 538 Impl::eigenValuesVectorsImpl<Impl::Jobs::EigenvaluesEigenvectors>(matrix, eigenValues, eigenVectors); 539 } 540 541 /** \brief calculates the eigenvalues of a symmetric field matrix 542 \param[in] matrix matrix eigenvalues are calculated for 543 \param[out] eigenValues FieldVector that contains eigenvalues in 544 ascending order 545 546 \note LAPACK::dsyev is used to calculate the eigenvalues 547 */ 548 template <int dim, typename K> eigenValuesLapack(const FieldMatrix<K,dim,dim> & matrix,FieldVector<K,dim> & eigenValues)549 static void eigenValuesLapack(const FieldMatrix<K, dim, dim>& matrix, 550 FieldVector<K, dim>& eigenValues) 551 { 552 Impl::EVDummy<K,dim> dummy; 553 Impl::eigenValuesVectorsLapackImpl<Impl::Jobs::EigenvaluesEigenvectors>(matrix, eigenValues, dummy); 554 } 555 556 /** \brief calculates the eigenvalues and -vectors of a symmetric field matrix 557 \param[in] matrix matrix eigenvalues are calculated for 558 \param[out] eigenValues FieldVector that contains eigenvalues in 559 ascending order 560 \param[out] eigenVectors FieldMatrix that contains the eigenvectors 561 562 \note LAPACK::dsyev is used to calculate the eigenvalues and -vectors 563 */ 564 template <int dim, typename K> eigenValuesVectorsLapack(const FieldMatrix<K,dim,dim> & matrix,FieldVector<K,dim> & eigenValues,FieldMatrix<K,dim,dim> & eigenVectors)565 static void eigenValuesVectorsLapack(const FieldMatrix<K, dim, dim>& matrix, 566 FieldVector<K, dim>& eigenValues, 567 FieldMatrix<K, dim, dim>& eigenVectors) 568 { 569 Impl::eigenValuesVectorsLapackImpl<Impl::Jobs::EigenvaluesEigenvectors>(matrix, eigenValues, eigenVectors); 570 } 571 572 /** \brief calculates the eigenvalues of a non-symmetric field matrix 573 \param[in] matrix matrix eigenvalues are calculated for 574 \param[out] eigenValues FieldVector that contains eigenvalues in 575 ascending order 576 577 \note LAPACK::dgeev is used to calculate the eigenvalues 578 */ 579 template <int dim, typename K, class C> eigenValuesNonSym(const FieldMatrix<K,dim,dim> & matrix,FieldVector<C,dim> & eigenValues)580 static void eigenValuesNonSym(const FieldMatrix<K, dim, dim>& matrix, 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 635 /** @} end documentation */ 636 637 } // end namespace Dune 638 #endif 639