00001
00002
00003 #ifndef DUNE_FMATRIXEIGENVALUES_HH
00004 #define DUNE_FMATRIXEIGENVALUES_HH
00005
00010 #include <iostream>
00011 #include <cmath>
00012 #include <cassert>
00013
00014 #include <dune/common/exceptions.hh>
00015 #include <dune/common/fvector.hh>
00016 #include <dune/common/fmatrix.hh>
00017
00018 namespace Dune {
00019
00025 namespace FMatrixHelp {
00026
00027
00028 extern void eigenValuesLapackCall(
00029 const char* jobz, const char* uplo, const long
00030 int* n, double* a, const long int* lda, double* w,
00031 double* work, const long int* lwork, long int* info);
00032
00033 extern void eigenValuesNonsymLapackCall(
00034 const char* jobvl, const char* jobvr, const long
00035 int* n, double* a, const long int* lda, double* wr, double* wi, double* vl,
00036 const long int* ldvl, double* vr, const long int* ldvr, double* work,
00037 const long int* lwork, const long int* info);
00038
00044 template <typename K>
00045 static void eigenValues(const FieldMatrix<K, 1, 1>& matrix,
00046 FieldVector<K, 1>& eigenvalues)
00047 {
00048 eigenvalues[0] = matrix[0][0];
00049 }
00050
00056 template <typename K>
00057 static void eigenValues(const FieldMatrix<K, 2, 2>& matrix,
00058 FieldVector<K, 2>& eigenvalues)
00059 {
00060 const K detM = matrix[0][0] * matrix[1][1] - matrix[1][0] * matrix[0][1];
00061 const K p = 0.5 * (matrix[0][0] + matrix [1][1]);
00062 K q = p * p - detM;
00063 if( q < 0 && q > -1e-14 ) q = 0;
00064 if (q < 0)
00065 {
00066 std::cout << matrix << std::endl;
00067
00068 DUNE_THROW(MathError, "Complex eigenvalue detected (which this implementation cannot handle).");
00069 }
00070
00071
00072 q = std :: sqrt(q);
00073
00074
00075 eigenvalues[0] = p - q;
00076 eigenvalues[1] = p + q;
00077 }
00078
00091 template <typename K>
00092 static void eigenValues(const FieldMatrix<K, 3, 3>& matrix,
00093 FieldVector<K, 3>& eigenvalues)
00094 {
00095 K p1 = matrix[0][1]*matrix[0][1] + matrix[0][2]*matrix[0][2] + matrix[1][2]*matrix[1][2];
00096
00097 if (p1 <= 1e-8)
00098 {
00099
00100 eigenvalues[0] = matrix[0][0];
00101 eigenvalues[1] = matrix[1][1];
00102 eigenvalues[2] = matrix[2][2];
00103 }
00104 else
00105 {
00106
00107 K q = 0;
00108 for (int i=0; i<3; i++)
00109 q += matrix[i][i]/3.0;
00110
00111 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;
00112 K p = std::sqrt(p2 / 6);
00113
00114 FieldMatrix<K,3,3> B;
00115 for (int i=0; i<3; i++)
00116 for (int j=0; j<3; j++)
00117 B[i][j] = (1/p) * (matrix[i][j] - q*(i==j));
00118
00119 K r = B.determinant() / 2.0;
00120
00121
00122
00123 K phi;
00124 if (r <= -1)
00125 phi = M_PI / 3.0;
00126 else if (r >= 1)
00127 phi = 0;
00128 else
00129 phi = std::acos(r) / 3;
00130
00131
00132 eigenvalues[2] = q + 2 * p * cos(phi);
00133 eigenvalues[0] = q + 2 * p * cos(phi + (2*M_PI/3));
00134 eigenvalues[1] = 3 * q - eigenvalues[0] - eigenvalues[2];
00135 }
00136 }
00137
00145 template <int dim, typename K>
00146 static void eigenValues(const FieldMatrix<K, dim, dim>& matrix,
00147 FieldVector<K, dim>& eigenvalues)
00148 {
00149 {
00150 const long int N = dim ;
00151 const char jobz = 'n';
00152 const char uplo = 'u';
00153
00154
00155 const long int w = N * N ;
00156
00157
00158 double matrixVector[dim * dim];
00159
00160
00161 int row = 0;
00162 for(int i=0; i<dim; ++i)
00163 {
00164 for(int j=0; j<dim; ++j, ++row)
00165 {
00166 matrixVector[ row ] = matrix[ i ][ j ];
00167 }
00168 }
00169
00170
00171 double workSpace[dim * dim];
00172
00173
00174 long int info = 0;
00175
00176
00177 eigenValuesLapackCall(&jobz, &uplo, &N, &matrixVector[0], &N,
00178 &eigenvalues[0], &workSpace[0], &w, &info);
00179
00180 if( info != 0 )
00181 {
00182 std::cerr << "For matrix " << matrix << " eigenvalue calculation failed! " << std::endl;
00183 DUNE_THROW(InvalidStateException,"eigenValues: Eigenvalue calculation failed!");
00184 }
00185 }
00186 }
00194 template <int dim, typename K, class C>
00195 static void eigenValuesNonSym(const FieldMatrix<K, dim, dim>& matrix,
00196 FieldVector<C, dim>& eigenValues)
00197 {
00198 {
00199 const long int N = dim ;
00200 const char jobvl = 'n';
00201 const char jobvr = 'n';
00202
00203
00204 double matrixVector[dim * dim];
00205
00206
00207 int row = 0;
00208 for(int i=0; i<dim; ++i)
00209 {
00210 for(int j=0; j<dim; ++j, ++row)
00211 {
00212 matrixVector[ row ] = matrix[ i ][ j ];
00213 }
00214 }
00215
00216
00217 double eigenR[dim];
00218 double eigenI[dim];
00219 double work[3*dim];
00220
00221
00222 long int info = 0;
00223 long int lwork = 3*dim;
00224
00225
00226 eigenValuesNonsymLapackCall(&jobvl, &jobvr, &N, &matrixVector[0], &N,
00227 &eigenR[0], &eigenI[0], 0, &N, 0, &N, &work[0],
00228 &lwork, &info);
00229
00230 if( info != 0 )
00231 {
00232 std::cerr << "For matrix " << matrix << " eigenvalue calculation failed! " << std::endl;
00233 DUNE_THROW(InvalidStateException,"eigenValues: Eigenvalue calculation failed!");
00234 }
00235 for (int i=0; i<N; ++i) {
00236 eigenValues[i].real = eigenR[i];
00237 eigenValues[i].imag = eigenI[i];
00238 }
00239 }
00240
00241 }
00242
00243 }
00244
00247 }
00248 #endif