00001
00002
00003 #ifndef DUNE_DYNMATRIXEIGENVALUES_HH
00004 #define DUNE_DYNMATRIXEIGENVALUES_HH
00005
00006 #include <memory>
00007
00008 #include <dune/common/std/memory.hh>
00009
00010 #include "dynmatrix.hh"
00011
00020 namespace Dune {
00021
00022 namespace DynamicMatrixHelp {
00023
00024
00025 extern void eigenValuesNonsymLapackCall(
00026 const char* jobvl, const char* jobvr, const long
00027 int* n, double* a, const long int* lda, double* wr, double* wi, double* vl,
00028 const long int* ldvl, double* vr, const long int* ldvr, double* work,
00029 const long int* lwork, const long int* info);
00030
00038 template <typename K, class C>
00039 static void eigenValuesNonSym(const DynamicMatrix<K>& matrix,
00040 DynamicVector<C>& eigenValues)
00041 {
00042 {
00043 const long int N = matrix.rows();
00044 const char jobvl = 'n';
00045 const char jobvr = 'n';
00046
00047
00048
00049 std::unique_ptr<double[]> matrixVector = Std::make_unique<double[]>(N*N);
00050
00051
00052 int row = 0;
00053 for(int i=0; i<N; ++i)
00054 {
00055 for(int j=0; j<N; ++j, ++row)
00056 {
00057 matrixVector[ row ] = matrix[ i ][ j ];
00058 }
00059 }
00060
00061
00062 std::unique_ptr<double[]> eigenR = Std::make_unique<double[]>(N);
00063 std::unique_ptr<double[]> eigenI = Std::make_unique<double[]>(N);
00064 std::unique_ptr<double[]> work = Std::make_unique<double[]>(3*N);
00065
00066
00067 long int info = 0;
00068 long int lwork = 3*N;
00069
00070
00071 eigenValuesNonsymLapackCall(&jobvl, &jobvr, &N, &matrixVector[0], &N,
00072 &eigenR[0], &eigenI[0], 0, &N, 0, &N, &work[0],
00073 &lwork, &info);
00074
00075 if( info != 0 )
00076 {
00077 std::cerr << "For matrix " << matrix << " eigenvalue calculation failed! " << std::endl;
00078 DUNE_THROW(InvalidStateException,"eigenValues: Eigenvalue calculation failed!");
00079 }
00080
00081 eigenValues.resize(N);
00082 for (int i=0; i<N; ++i)
00083 eigenValues[i] = std::complex<double>(eigenR[i], eigenI[i]);
00084 }
00085 }
00086
00087 }
00088
00089 }
00091 #endif