dune-common  2.3.1
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 "dynmatrix.hh"
7 
16 namespace Dune {
17 
18  namespace DynamicMatrixHelp {
19 
20  // defined in fmatrixev_ext.cpp
21  extern void eigenValuesNonsymLapackCall(
22  const char* jobvl, const char* jobvr, const long
23  int* n, double* a, const long int* lda, double* wr, double* wi, double* vl,
24  const long int* ldvl, double* vr, const long int* ldvr, double* work,
25  const long int* lwork, const long int* info);
26 
34  template <typename K, class C>
35  static void eigenValuesNonSym(const DynamicMatrix<K>& matrix,
37  {
38  {
39  const long int N = matrix.rows();
40  const char jobvl = 'n';
41  const char jobvr = 'n';
42 
43 
44  // matrix to put into dgeev
45  double matrixVector[N * N];
46 
47  // copy matrix
48  int row = 0;
49  for(int i=0; i<N; ++i)
50  {
51  for(int j=0; j<N; ++j, ++row)
52  {
53  matrixVector[ row ] = matrix[ i ][ j ];
54  }
55  }
56 
57  // working memory
58  double eigenR[N];
59  double eigenI[N];
60  double work[3*N];
61 
62  // return value information
63  long int info = 0;
64  long int lwork = 3*N;
65 
66  // call LAPACK routine (see fmatrixev_ext.cc)
67  eigenValuesNonsymLapackCall(&jobvl, &jobvr, &N, &matrixVector[0], &N,
68  &eigenR[0], &eigenI[0], 0, &N, 0, &N, &work[0],
69  &lwork, &info);
70 
71  if( info != 0 )
72  {
73  std::cerr << "For matrix " << matrix << " eigenvalue calculation failed! " << std::endl;
74  DUNE_THROW(InvalidStateException,"eigenValues: Eigenvalue calculation failed!");
75  }
76  for (int i=0; i<N; ++i)
77  eigenValues[i] = std::complex<double>(eigenR[i], eigenI[i]);
78  }
79  }
80 
81  }
82 
83 }
85 #endif