From 3909180bf8e49e5ebaf9e0a72acd8237b0156b29 Mon Sep 17 00:00:00 2001 From: Benjamin Brands Date: Mon, 19 Mar 2018 15:07:55 +0100 Subject: [PATCH] add wrappers for LAPACK routines --- include/deal.II/lac/lapack_templates.h | 197 +++++++++++++++++++++++++ 1 file changed, 197 insertions(+) diff --git a/include/deal.II/lac/lapack_templates.h b/include/deal.II/lac/lapack_templates.h index bd894eaafc..cab2707058 100644 --- a/include/deal.II/lac/lapack_templates.h +++ b/include/deal.II/lac/lapack_templates.h @@ -221,6 +221,25 @@ extern "C" dealii::types::blas_int *m, float *w, float *z, const dealii::types::blas_int *ldz, float *work, const dealii::types::blas_int *lwork, dealii::types::blas_int *iwork, dealii::types::blas_int *ifail, dealii::types::blas_int *info); + // Same functionality as dsyev_ but using MRRR algorithm and with more options: + // E.g. compute only eigenvalues in a specific dealii::types::blas_interval, + // Compute only eigenvalues with a specific index. + void dsyevr_(const char *jobz, const char *range, + const char *uplo, const dealii::types::blas_int *n, double *A, const dealii::types::blas_int *lda, + const double *vl, const double *vu, + const dealii::types::blas_int *il, const dealii::types::blas_int *iu, const double *abstol, + dealii::types::blas_int *m, double *w, double *z, + const dealii::types::blas_int *ldz, dealii::types::blas_int *isuppz, + double *work, dealii::types::blas_int *lwork, dealii::types::blas_int *iwork, + dealii::types::blas_int *liwork, dealii::types::blas_int *info); + void ssyevr_(const char *jobz, const char *range, + const char *uplo, const dealii::types::blas_int *n, float *A, const dealii::types::blas_int *lda, + const float *vl, const float *vu, + const dealii::types::blas_int *il, const dealii::types::blas_int *iu, const float *abstol, + dealii::types::blas_int *m, float *w, float *z, + const dealii::types::blas_int *ldz, dealii::types::blas_int *isuppz, + float *work, dealii::types::blas_int *lwork, dealii::types::blas_int *iwork, + dealii::types::blas_int *liwork, dealii::types::blas_int *info); // Generalized eigenvalues and eigenvectors of // 1: A*x = lambda*B*x; 2: A*B*x = lambda*x; 3: B*A*x = lambda*x // A and B are symmetric and B is definite @@ -336,6 +355,9 @@ extern "C" float *A, const dealii::types::blas_int *lda, dealii::types::blas_int *info); +// dlamch and slamch help determining machine precision + double dlamch_(const char *chmach); + float slamch_(const char *chmach); } DEAL_II_NAMESPACE_OPEN @@ -1326,6 +1348,145 @@ syevx (const char *, const char *, const char *, const types::blas_int *, float #endif +// Template wrapper for LAPACK functions dsyevr and ssyevr +template +inline void +syevr(const char * /*jobz*/, + const char * /*range*/, + const char * /*uplo*/, + const types::blas_int * /*n*/, + number * /*A*/, + const types::blas_int * /*lda*/, + const number * /*vl*/, + const number * /*vu*/, + const types::blas_int * /*il*/, + const types::blas_int * /*iu*/, + const number * /*abstol*/, + types::blas_int * /*m*/, + number * /*w*/, + number * /*z*/, + const types::blas_int * /*ldz*/, + types::blas_int * /*isuppz*/, + number * /*work*/, + types::blas_int * /*lwork*/, + types::blas_int * /*iwork*/, + types::blas_int * /*liwork*/, + types::blas_int * /*info*/) +{ + Assert(false, ExcNotImplemented()); +} + + +#ifdef DEAL_II_WITH_LAPACK +inline void +syevr(const char *jobz, + const char *range, + const char *uplo, + const types::blas_int *n, + double *A, + const types::blas_int *lda, + const double *vl, + const double *vu, + const types::blas_int *il, + const types::blas_int *iu, + const double *abstol, + types::blas_int *m, + double *w, + double *z, + const types::blas_int *ldz, + types::blas_int *isuppz, + double *work, + types::blas_int *lwork, + types::blas_int *iwork, + types::blas_int *liwork, + types::blas_int *info) +{ + dsyevr_(jobz,range,uplo,n,A,lda,vl,vu,il,iu,abstol,m,w,z,ldz,isuppz,work,lwork,iwork,liwork,info); +} +#else +inline void +syevr(const char * /*jobz*/, + const char * /*range*/, + const char * /*uplo*/, + const types::blas_int * /*n*/, + double * /*A*/, + const types::blas_int * /*lda*/, + const double * /*vl*/, + const double * /*vu*/, + const types::blas_int * /*il*/, + const types::blas_int * /*iu*/, + const double * /*abstol*/, + types::blas_int * /*m*/, + double * /*w*/, + double * /*z*/, + const types::blas_int * /*ldz*/, + types::blas_int * /*isuppz*/, + double * /*work*/, + types::blas_int * /*lwork*/, + types::blas_int * /*iwork*/, + types::blas_int * /*liwork*/, + types::blas_int * /*info*/) +{ + Assert (false, LAPACKSupport::ExcMissing("dsyevr")); +} +#endif + + +#ifdef DEAL_II_WITH_LAPACK +inline void +syevr(const char *jobz, + const char *range, + const char *uplo, + const types::blas_int *n, + float *A, + const types::blas_int *lda, + const float *vl, + const float *vu, + const types::blas_int *il, + const types::blas_int *iu, + const float *abstol, + types::blas_int *m, + float *w, + float *z, + const types::blas_int *ldz, + types::blas_int *isuppz, + float *work, + types::blas_int *lwork, + types::blas_int *iwork, + types::blas_int *liwork, + types::blas_int *info) +{ + ssyevr_(jobz,range,uplo,n,A,lda,vl,vu,il,iu,abstol,m,w,z,ldz,isuppz,work,lwork,iwork,liwork,info); +} +#else +inline void +syevr(const char * /*jobz*/, + const char * /*range*/, + const char * /*uplo*/, + const types::blas_int * /*n*/, + float * /*A*/, + const types::blas_int * /*lda*/, + const float * /*vl*/, + const float * /*vu*/, + const types::blas_int * /*il*/, + const types::blas_int * /*iu*/, + const float * /*abstol*/, + types::blas_int * /*m*/, + float * /*w*/, + float * /*z*/, + const types::blas_int * /*ldz*/, + types::blas_int * /*isuppz*/, + float * /*work*/, + types::blas_int * /*lwork*/, + types::blas_int * /*iwork*/, + types::blas_int * /*liwork*/, + types::blas_int * /*info*/) +{ + Assert (false, LAPACKSupport::ExcMissing("ssyevr")); +} +#endif + + /// Template wrapper for LAPACK functions dsygv and ssygv template inline void @@ -1620,6 +1781,42 @@ lascl(const char *, const types::blas_int *, const types::blas_int *, const floa #endif +/// Template wrapper for LAPACK functions dlamch and slamch +template +inline void +lamch(const char * /*chmach*/, number & /*precision*/) +{ + Assert(false, ExcNotImplemented()); +} + +#ifdef DEAL_II_WITH_LAPACK +inline void +lamch(const char *chmach, double &precision) +{ + precision = dlamch_(chmach); +} +#else +inline void +lamch(const char * /*chmach*/, double & /*precision*/) +{ + Assert(false, LAPACKSupport::ExcMissing("dlamch")); +} +#endif + +#ifdef DEAL_II_WITH_LAPACK +inline void +lamch(const char *chmach, float &precision) +{ + precision = slamch_(chmach); +} +#else +inline void +lamch(const char * /*chmach*/, float & /*precision*/) +{ + Assert(false, LAPACKSupport::ExcMissing("slamch")); +} +#endif + DEAL_II_NAMESPACE_CLOSE #endif -- 2.39.5