dune-common  2.8.0
dynmatrixev.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_DYNMATRIXEIGENVALUES_HH
4 #define DUNE_DYNMATRIXEIGENVALUES_HH
5 
6 #include <algorithm>
7 #include <memory>
8 
9 #include "dynmatrix.hh"
10 #include "fmatrixev.hh"
11 
20 namespace Dune {
21 
22  namespace DynamicMatrixHelp {
23 
24 #if HAVE_LAPACK
25  using Dune::FMatrixHelp::eigenValuesNonsymLapackCall;
26 #endif
27 
36  template <typename K, class C>
37  static void eigenValuesNonSym(const DynamicMatrix<K>& matrix,
39  std::vector<DynamicVector<K>>* eigenVectors = nullptr
40  )
41  {
42 
43 #if HAVE_LAPACK
44  {
45  const long int N = matrix.rows();
46  const char jobvl = 'n';
47  const char jobvr = eigenVectors ? 'v' : 'n';
48 
49 
50  // matrix to put into dgeev
51  auto matrixVector = std::make_unique<double[]>(N*N);
52 
53  // copy matrix
54  int row = 0;
55  for(int i=0; i<N; ++i)
56  {
57  for(int j=0; j<N; ++j, ++row)
58  {
59  matrixVector[ row ] = matrix[ i ][ j ];
60  }
61  }
62 
63  // working memory
64  auto eigenR = std::make_unique<double[]>(N);
65  auto eigenI = std::make_unique<double[]>(N);
66 
67  const long int lwork = eigenVectors ? 4*N : 3*N;
68  auto work = std::make_unique<double[]>(lwork);
69  auto vr = eigenVectors ? std::make_unique<double[]>(N*N) : std::unique_ptr<double[]>{};
70 
71  // return value information
72  long int info = 0;
73 
74  // call LAPACK routine (see fmatrixev_ext.cc)
75  eigenValuesNonsymLapackCall(&jobvl, &jobvr, &N, matrixVector.get(), &N,
76  eigenR.get(), eigenI.get(), nullptr, &N, vr.get(), &N, work.get(),
77  &lwork, &info);
78 
79  if( info != 0 )
80  {
81  std::cerr << "For matrix " << matrix << " eigenvalue calculation failed! " << std::endl;
82  DUNE_THROW(InvalidStateException,"eigenValues: Eigenvalue calculation failed!");
83  }
84 
85  eigenValues.resize(N);
86  for (int i=0; i<N; ++i)
87  eigenValues[i] = std::complex<double>(eigenR[i], eigenI[i]);
88 
89  if (eigenVectors) {
90  eigenVectors->resize(N);
91  for (int i = 0; i < N; ++i) {
92  auto& v = (*eigenVectors)[i];
93  v.resize(N);
94  std::copy(vr.get() + N*i, vr.get() + N*(i+1), &v[0]);
95  }
96  }
97  }
98 #else // #if HAVE_LAPACK
99  DUNE_THROW(NotImplemented,"LAPACK not found!");
100 #endif
101  }
102  }
103 
104 }
106 #endif
This file implements a dense matrix with dynamic numbers of rows and columns.
Eigenvalue computations for the FieldMatrix class.
#define DUNE_THROW(E, m)
Definition: exceptions.hh:216
Dune namespace.
Definition: alignedallocator.hh:11
static void eigenValuesNonSym(const DynamicMatrix< K > &matrix, DynamicVector< C > &eigenValues, std::vector< DynamicVector< K >> *eigenVectors=nullptr)
calculates the eigenvalues of a symmetric field matrix
Definition: dynmatrixev.hh:37
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
constexpr size_type rows() const
number of rows
Definition: densematrix.hh:737
Construct a matrix with a dynamic size.
Definition: dynmatrix.hh:59
Construct a vector with a dynamic size.
Definition: dynvector.hh:57
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