Skip to content
Snippets Groups Projects
Verified Commit 876989bc authored by Nils-Arne Dreier's avatar Nils-Arne Dreier
Browse files

add support for computing eigenvalues and vectors for FieldMatrices 

with field_type float
parent 27059e0f
Branches
Tags
1 merge request!787add eigenvalues and vectors computation for FieldMatrices for field_type float
Pipeline #25603 passed
Stage: test
Hide whitespace changes
Inline Side-by-side
dune/common/fmatrixev.cc
+ 46
− 0
View file @ 876989bc
......@@ -18,9 +18,11 @@
// symmetric matrices
#define DSYEV_FORTRAN FC_FUNC (dsyev, DSYEV)
#define SSYEV_FORTRAN FC_FUNC (ssyev, SSYEV)
// nonsymmetric matrices
#define DGEEV_FORTRAN FC_FUNC (dgeev, DGEEV)
#define SGEEV_FORTRAN FC_FUNC (sgeev, SGEEV)
// dsyev declaration (in liblapack)
extern "C" {
......@@ -90,6 +92,9 @@ extern "C" {
extern void DSYEV_FORTRAN(const char* jobz, const char* uplo, const long
int* n, double* a, const long int* lda, double* w,
double* work, const long int* lwork, long int* info);
extern void SSYEV_FORTRAN(const char* jobz, const char* uplo, const long
int* n, float* a, const long int* lda, float* w,
float* work, const long int* lwork, long int* info);
/*
*
......@@ -195,6 +200,10 @@ extern "C" {
int* n, double* a, const long int* lda, double* wr, double* wi, double* vl,
const long int* ldvl, double* vr, const long int* ldvr, double* work,
const long int* lwork, long int* info);
extern void SGEEV_FORTRAN(const char* jobvl, const char* jobvr, const long
int* n, float* a, const long int* lda, float* wr, float* wi, float* vl,
const long int* ldvl, float* vr, const long int* ldvr, float* work,
const long int* lwork, long int* info);
} // end extern C
#endif
......@@ -220,6 +229,23 @@ namespace Dune {
#endif
}
void eigenValuesLapackCall(
const char* jobz, const char* uplo, const long
int* n, float* a, const long int* lda, float* w,
float* work, const long int* lwork, long int* info)
{
#if HAVE_LAPACK
// call LAPACK dsyev
SSYEV_FORTRAN(jobz, uplo, n, a, lda, w, work, lwork, info);
#else
// silence unused variable warnings
DUNE_UNUSED_PARAMETER(jobz), DUNE_UNUSED_PARAMETER(uplo), DUNE_UNUSED_PARAMETER(n);
DUNE_UNUSED_PARAMETER(a), DUNE_UNUSED_PARAMETER(lda), DUNE_UNUSED_PARAMETER(w);
DUNE_UNUSED_PARAMETER(work), DUNE_UNUSED_PARAMETER(lwork), DUNE_UNUSED_PARAMETER(info);
DUNE_THROW(NotImplemented,"eigenValuesLapackCall: LAPACK not found!");
#endif
}
void eigenValuesNonsymLapackCall(
const char* jobvl, const char* jobvr, const long
int* n, double* a, const long int* lda, double* wr, double* wi, double* vl,
......@@ -240,6 +266,26 @@ namespace Dune {
#endif
}
void eigenValuesNonsymLapackCall(
const char* jobvl, const char* jobvr, const long
int* n, float* a, const long int* lda, float* wr, float* wi, float* vl,
const long int* ldvl, float* vr, const long int* ldvr, float* work,
const long int* lwork, long int* info)
{
#if HAVE_LAPACK
// call LAPACK dgeev
SGEEV_FORTRAN(jobvl, jobvr, n, a, lda, wr, wi, vl, ldvl, vr, ldvr,
work, lwork, info);
#else
// silence unused variable warnings
DUNE_UNUSED_PARAMETER(jobvl), DUNE_UNUSED_PARAMETER(jobvr), DUNE_UNUSED_PARAMETER(n);
DUNE_UNUSED_PARAMETER(a), DUNE_UNUSED_PARAMETER(lda), DUNE_UNUSED_PARAMETER(wr), DUNE_UNUSED_PARAMETER(wi);
DUNE_UNUSED_PARAMETER(vl), DUNE_UNUSED_PARAMETER(ldvl), DUNE_UNUSED_PARAMETER(vr);
DUNE_UNUSED_PARAMETER(ldvr), DUNE_UNUSED_PARAMETER(work), DUNE_UNUSED_PARAMETER(lwork), DUNE_UNUSED_PARAMETER(info);
DUNE_THROW(NotImplemented,"eigenValuesNonsymLapackCall: LAPACK not found!");
#endif
}
} // end namespace FMatrixHelp
} // end namespace Dune
......
dune/common/fmatrixev.hh
+ 18
− 7
View file @ 876989bc
......@@ -38,6 +38,17 @@ namespace Dune {
const long int* ldvl, double* vr, const long int* ldvr, double* work,
const long int* lwork, long int* info);
extern void eigenValuesLapackCall(
const char* jobz, const char* uplo, const long
int* n, float* a, const long int* lda, float* w,
float* work, const long int* lwork, long int* info);
extern void eigenValuesNonsymLapackCall(
const char* jobvl, const char* jobvr, const long
int* n, float* a, const long int* lda, float* wr, float* wi, float* vl,
const long int* ldvl, float* vr, const long int* ldvr, float* work,
const long int* lwork, long int* info);
namespace Impl {
//internal tag to activate/disable code for eigenvector calculation at compile time
enum Jobs { OnlyEigenvalues=0, EigenvaluesEigenvectors=1 };
......@@ -124,7 +135,7 @@ namespace Dune {
acos(z) function requires |z| <= 1, but will fail silently
and return NaN if the input is larger than 1 in magnitude.
Thus r is clamped to [-1,1].*/
r = std::min(std::max(r, -1.0), 1.0);
r = std::min<K>(std::max<K>(r, -1.0), 1.0);
K phi = acos(r) / 3.0;
// the eigenvalues satisfy eig[2] <= eig[1] <= eig[0]
......@@ -385,7 +396,7 @@ namespace Dune {
const long int lwork = 3*N -1 ;
// matrix to put into dsyev
double matrixVector[dim * dim];
K matrixVector[dim * dim];
// copy matrix
int row = 0;
......@@ -398,7 +409,7 @@ namespace Dune {
}
// working memory
double workSpace[lwork];
K workSpace[lwork];
// return value information
long int info = 0;
......@@ -516,7 +527,7 @@ namespace Dune {
const char jobvr = 'n';
// matrix to put into dgeev
double matrixVector[dim * dim];
K matrixVector[dim * dim];
// copy matrix
int row = 0;
......@@ -529,9 +540,9 @@ namespace Dune {
}
// working memory
double eigenR[dim];
double eigenI[dim];
double work[3*dim];
K eigenR[dim];
K eigenI[dim];
K work[3*dim];
// return value information
long int info = 0;
......
dune/common/test/eigenvaluestest.cc
+ 20
− 13
View file @ 876989bc
......@@ -125,13 +125,13 @@ void testSymmetricFieldMatrix()
{
FieldVector<field_type, dim> Av;
testMatrix.mv(eigenVectors[j], Av);
if((Av - eigenValues[j]*eigenVectors[j]).two_norm() > 1e-8)
if((Av - eigenValues[j]*eigenVectors[j]).two_norm() > dim*std::sqrt(std::numeric_limits<field_type>::epsilon()))
DUNE_THROW(MathError, "Vector computed by FMatrixHelp::eigenValuesVectors is not an eigenvector");
}
// Make sure the eigenvectors have unit length
for(auto& ev : eigenVectors) {
if(std::abs(ev.two_norm())-1 > 1e-14)
if(std::abs(ev.two_norm())-1 > dim*std::numeric_limits<field_type>::epsilon())
DUNE_THROW(MathError, "Vector computed by FMatrixHelp::eigenValuesVectors does not have unit length");
}
......@@ -221,32 +221,33 @@ void checkMatrixWithLAPACK(FieldMatrix<field_type, dim, dim> matrix)
}
}
template<class FT>
void checkMultiplicity()
{
//--2d--
//repeated eigenvalue (x2)
checkMatrixWithReference<double,2>({{1, 0},{0, 1}}, {{1,0}, {0,1}}, {1, 1});
checkMatrixWithReference<FT,2>({{1, 0},{0, 1}}, {{1,0}, {0,1}}, {1, 1});
//eigenvalues with same magnitude (x2)
checkMatrixWithReference<double,2>({{0, 1}, {1, 0}}, {{1,-1}, {1,1}}, {-1, 1});
checkMatrixWithReference<FT,2>({{0, 1}, {1, 0}}, {{1,-1}, {1,1}}, {-1, 1});
//--3d--
//repeated eigenvalue (x3)
checkMatrixWithReference<double,3>({{ 1, 0, 0},
checkMatrixWithReference<FT,3>({{ 1, 0, 0},
{ 0, 1, 0},
{ 0, 0, 1}},
{{1,0,0}, {0,1,0}, {0,0,1}},
{1, 1, 1});
//eigenvalues with same magnitude (x2)
checkMatrixWithReference<double,3>({{ 0, 1, 0},
checkMatrixWithReference<FT,3>({{ 0, 1, 0},
{ 1, 0, 0},
{ 0, 0, 5}},
{{-1,1,0}, {1,1,0}, {0,0,1}},
{-1, 1, 5});
//repeated eigenvalue (x2)
checkMatrixWithReference<double,3>({{ 3, -2, 0},
checkMatrixWithReference<FT,3>({{ 3, -2, 0},
{ -2, 3, 0},
{ 0, 0, 5}},
{{1,1,0}, {0,0,1}, {1,-1,0}},
......@@ -254,11 +255,11 @@ void checkMultiplicity()
//repeat tests with LAPACK (if found)
#if HAVE_LAPACK
checkMatrixWithLAPACK<double,2>({{1, 0}, {0, 1}});
checkMatrixWithLAPACK<double,2>({{0, 1}, {1, 0}});
checkMatrixWithLAPACK<double,3>({{1,0,0}, {0,1,0}, {0,0,1}});
checkMatrixWithLAPACK<double,3>({{0,1,0}, {1,0,0}, {0,0,5}});
checkMatrixWithLAPACK<double,3>({{3,-2,0}, {-2,3,0}, {0,0,5}});
checkMatrixWithLAPACK<FT,2>({{1, 0}, {0, 1}});
checkMatrixWithLAPACK<FT,2>({{0, 1}, {1, 0}});
checkMatrixWithLAPACK<FT,3>({{1,0,0}, {0,1,0}, {0,0,1}});
checkMatrixWithLAPACK<FT,3>({{0,1,0}, {1,0,0}, {0,0,5}});
checkMatrixWithLAPACK<FT,3>({{3,-2,0}, {-2,3,0}, {0,0,5}});
#endif
}
......@@ -267,6 +268,7 @@ int main()
{
#if HAVE_LAPACK
testRosserMatrix<double>();
testRosserMatrix<float>();
#else
std::cout << "WARNING: eigenvaluetest needs LAPACK, test disabled" << std::endl;
#endif // HAVE_LAPACK
......@@ -275,12 +277,17 @@ int main()
#if HAVE_LAPACK
testSymmetricFieldMatrix<double,4>();
testSymmetricFieldMatrix<double,200>();
testSymmetricFieldMatrix<float,4>();
testSymmetricFieldMatrix<float,200>();
#endif // HAVE_LAPACK
testSymmetricFieldMatrix<double,2>();
testSymmetricFieldMatrix<double,3>();
testSymmetricFieldMatrix<float,2>();
testSymmetricFieldMatrix<float,3>();
checkMultiplicity();
checkMultiplicity<double>();
checkMultiplicity<float>();
return 0;
}
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment