From: kanschat Date: Sun, 1 Oct 2006 23:01:21 +0000 (+0000) Subject: LU decomposition added to LAPACKFullMatrix X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=cbb9ebeaffa26d4a4f49023adf4caa4a63eef593;p=dealii-svn.git LU decomposition added to LAPACKFullMatrix some items in LAPACKFullMatrix fixed git-svn-id: https://svn.dealii.org/trunk@13973 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/base/include/base/config.h.in b/deal.II/base/include/base/config.h.in index e1016ea246..0f7d6ba18d 100644 --- a/deal.II/base/include/base/config.h.in +++ b/deal.II/base/include/base/config.h.in @@ -221,6 +221,12 @@ /* Define to 1 if you have the `dgesvd_' function. */ #undef HAVE_DGESVD_ +/* Define to 1 if you have the `dgetrf_' function. */ +#undef HAVE_DGETRF_ + +/* Define to 1 if you have the `dgetrs_' function. */ +#undef HAVE_DGETRS_ + /* Define to 1 if you have the `dstev_' function. */ #undef HAVE_DSTEV_ @@ -280,6 +286,12 @@ /* Define to 1 if you have the `sgesvd_' function. */ #undef HAVE_SGESVD_ +/* Define to 1 if you have the `sgetrf_' function. */ +#undef HAVE_SGETRF_ + +/* Define to 1 if you have the `sgetrs_' function. */ +#undef HAVE_SGETRS_ + /* Define to 1 if you have the `sstev_' function. */ #undef HAVE_SSTEV_ diff --git a/deal.II/configure b/deal.II/configure index b8492dedbe..82b3612a3a 100755 --- a/deal.II/configure +++ b/deal.II/configure @@ -12067,9 +12067,114 @@ fi; +for ac_func in dgemv_ sgemv_ dgeev_ sgeev_ dgeevx_ sgeevx_ dgesvd_ sgesvd_ +do +as_ac_var=`echo "ac_cv_func_$ac_func" | $as_tr_sh` +echo "$as_me:$LINENO: checking for $ac_func" >&5 +echo $ECHO_N "checking for $ac_func... $ECHO_C" >&6 +if eval "test \"\${$as_ac_var+set}\" = set"; then + echo $ECHO_N "(cached) $ECHO_C" >&6 +else + cat >conftest.$ac_ext <<_ACEOF +/* confdefs.h. */ +_ACEOF +cat confdefs.h >>conftest.$ac_ext +cat >>conftest.$ac_ext <<_ACEOF +/* end confdefs.h. */ +/* Define $ac_func to an innocuous variant, in case declares $ac_func. + For example, HP-UX 11i declares gettimeofday. */ +#define $ac_func innocuous_$ac_func + +/* System header to define __stub macros and hopefully few prototypes, + which can conflict with char $ac_func (); below. + Prefer to if __STDC__ is defined, since + exists even on freestanding compilers. */ + +#ifdef __STDC__ +# include +#else +# include +#endif + +#undef $ac_func + +/* Override any gcc2 internal prototype to avoid an error. */ +#ifdef __cplusplus +extern "C" +{ +#endif +/* We use char because int might match the return type of a gcc2 + builtin and then its argument prototype would still apply. */ +char $ac_func (); +/* The GNU C library defines this for functions which it implements + to always fail with ENOSYS. Some functions are actually named + something starting with __ and the normal name is an alias. */ +#if defined (__stub_$ac_func) || defined (__stub___$ac_func) +choke me +#else +char (*f) () = $ac_func; +#endif +#ifdef __cplusplus +} +#endif + +int +main () +{ +return f != $ac_func; + ; + return 0; +} +_ACEOF +rm -f conftest.$ac_objext conftest$ac_exeext +if { (eval echo "$as_me:$LINENO: \"$ac_link\"") >&5 + (eval $ac_link) 2>conftest.er1 + ac_status=$? + grep -v '^ *+' conftest.er1 >conftest.err + rm -f conftest.er1 + cat conftest.err >&5 + echo "$as_me:$LINENO: \$? = $ac_status" >&5 + (exit $ac_status); } && + { ac_try='test -z "$ac_cxx_werror_flag" + || test ! -s conftest.err' + { (eval echo "$as_me:$LINENO: \"$ac_try\"") >&5 + (eval $ac_try) 2>&5 + ac_status=$? + echo "$as_me:$LINENO: \$? = $ac_status" >&5 + (exit $ac_status); }; } && + { ac_try='test -s conftest$ac_exeext' + { (eval echo "$as_me:$LINENO: \"$ac_try\"") >&5 + (eval $ac_try) 2>&5 + ac_status=$? + echo "$as_me:$LINENO: \$? = $ac_status" >&5 + (exit $ac_status); }; }; then + eval "$as_ac_var=yes" +else + echo "$as_me: failed program was:" >&5 +sed 's/^/| /' conftest.$ac_ext >&5 + +eval "$as_ac_var=no" +fi +rm -f conftest.err conftest.$ac_objext \ + conftest$ac_exeext conftest.$ac_ext +fi +echo "$as_me:$LINENO: result: `eval echo '${'$as_ac_var'}'`" >&5 +echo "${ECHO_T}`eval echo '${'$as_ac_var'}'`" >&6 +if test `eval echo '${'$as_ac_var'}'` = yes; then + cat >>confdefs.h <<_ACEOF +#define `echo "HAVE_$ac_func" | $as_tr_cpp` 1 +_ACEOF + +fi +done + + + + + -for ac_func in dgemv_ sgemv_ dgeev_ sgeev_ dgeevx_ sgeevx_ dgesvd_ sgesvd_ dstev_ sstev_ +for ac_func in dgetrf_ sgetrf_ dgetrs_ sgetrs_ dstev_ sstev_ do as_ac_var=`echo "ac_cv_func_$ac_func" | $as_tr_sh` echo "$as_me:$LINENO: checking for $ac_func" >&5 diff --git a/deal.II/configure.in b/deal.II/configure.in index ce8c038823..1ee4c2c5ce 100644 --- a/deal.II/configure.in +++ b/deal.II/configure.in @@ -510,7 +510,8 @@ AC_ARG_WITH(lapack, argument is given, use -llapack. Default is not to use -llapack.], DEAL_II_WITH_LAPACK($withval)) -AC_CHECK_FUNCS([dgemv_ sgemv_ dgeev_ sgeev_ dgeevx_ sgeevx_ dgesvd_ sgesvd_ dstev_ sstev_]) +AC_CHECK_FUNCS([dgemv_ sgemv_ dgeev_ sgeev_ dgeevx_ sgeevx_ dgesvd_ sgesvd_]) +AC_CHECK_FUNCS([dgetrf_ sgetrf_ dgetrs_ sgetrs_ dstev_ sstev_]) dnl ------------------------------------------------------------- dnl set include paths of several libraries diff --git a/deal.II/contrib/blastemplates/lapack_templates.h.in b/deal.II/contrib/blastemplates/lapack_templates.h.in index 1b23d0d64f..89c8578f2b 100644 --- a/deal.II/contrib/blastemplates/lapack_templates.h.in +++ b/deal.II/contrib/blastemplates/lapack_templates.h.in @@ -4,6 +4,15 @@ void dgemv_ (const char* trans, const int* m, const int* n, const double* alpha, const double* A, const int* lda, const double* x, const int* incx, const double* b, double* y, const int* incy); + +// Compute LU factorization +void dgetrf_ (const int* m, const int* n, double* A, + const int* lda, int* ipiv, int* info); +// Apply forward/backward substitution to LU factorization +void dgetrs_ (const char* trans, const int* n, const int* nrhs, + const double* A, const int* lda, const int* ipiv, + double* b, const int* ldb, int* info); + // Compute eigenvalues and vectors void dgeev_ (const char* jobvl, const char* jobvr, const int* n, double* A, const int* lda, diff --git a/deal.II/lac/include/lac/exceptions.h b/deal.II/lac/include/lac/exceptions.h index 0e63732d8c..d286ad19fe 100644 --- a/deal.II/lac/include/lac/exceptions.h +++ b/deal.II/lac/include/lac/exceptions.h @@ -2,7 +2,7 @@ // $Id$ // Version: $Name$ // -// Copyright (C) 2004, 2005 by the deal.II authors +// Copyright (C) 2004, 2005, 2006 by the deal.II authors // // This file is subject to QPL and may not be distributed // without copyright and license information. Please refer @@ -28,11 +28,17 @@ namespace LACExceptions */ DeclException0 (ExcNotQuadratic); + /** + * The operation cannot be finished + * since the matrix is singular. + */ + DeclException0 (ExcSingular); + /** * Block indices of two block * objects are different. */ - DeclException0(ExcDifferentBlockIndices); + DeclException0 (ExcDifferentBlockIndices); /** * An error of a PETSc function was diff --git a/deal.II/lac/include/lac/lapack_full_matrix.h b/deal.II/lac/include/lac/lapack_full_matrix.h index 945fe88b71..6f00e2b46e 100644 --- a/deal.II/lac/include/lac/lapack_full_matrix.h +++ b/deal.II/lac/include/lac/lapack_full_matrix.h @@ -15,24 +15,33 @@ #include +#include #include #include +#include #include #include // forward declarations template class Vector; +template class BlockVector; template class FullMatrix; - -/*! @addtogroup Matrix1 - *@{ - */ +template class VectorMemory; /** - * A variant of FullMatrix using LAPACK functions whereever possible. + * A variant of FullMatrix using LAPACK functions whereever + * possible. In order to do this, the matrix is stored in transposed + * order. The element access functions hide this fact by reverting the + * transposition. * + * @note In order to perform LAPACK functions, the class contains a lot of + * auxiliary data in the private section. The names of these data + * vectors are usually the names chosen for the arguments in the + * LAPACK documentation. + * + * @ingroup Matrix1 * @author Guido Kanschat, 2005 */ template @@ -136,13 +145,19 @@ class LAPACKFullMatrix : public TransposeTable * maximum size possible, * determined either by the size * of this or src. + * + * The final two arguments allow + * to enter a multiple of the + * source or its transpose. */ template void fill (const MATRIX &src, const unsigned int dst_offset_i = 0, const unsigned int dst_offset_j = 0, const unsigned int src_offset_i = 0, - const unsigned int src_offset_j = 0); + const unsigned int src_offset_j = 0, + const number factor = 1., + const bool transpose = false); /** * Matrix-vector-multiplication. @@ -208,6 +223,25 @@ class LAPACKFullMatrix : public TransposeTable void Tvmult_add (Vector &w, const Vector &v) const; + /** + * Compute the LU factorization + * of the matrix using LAPACK + * function Xgetrf. + */ + void compute_lu_factorization (); + + /** + * Solve the linear system with + * right hand side given by + * applying forward/backward + * substitution to the previously + * computed LU + * factorization. Uses LAPACK + * function Xgetrs. + */ + void apply_lu_factorization (Vector& v, + const bool transposed) const; + /** * Compute eigenvalues of the * matrix. After this routine has @@ -326,6 +360,14 @@ class LAPACKFullMatrix : public TransposeTable * some LAPACK functions. */ mutable std::vector work; + + /** + * The vector storing the + * permutations applied for + * pivoting in the + * LU-factorization. + */ + std::vector ipiv; /** * Real parts of @@ -354,7 +396,35 @@ class LAPACKFullMatrix : public TransposeTable std::vector vr; }; -/**@}*/ + + +/** + * A preconditioner based on the LU-factorization of LAPACKFullMatrix. + * + * @ingroup Preconditioners + * @author Guido Kanschat, 2006 + */ +template +class LUPrecondition + : + public Subscriptor +{ + public: + void initialize(const LAPACKFullMatrix&); + void initialize(const LAPACKFullMatrix&, + VectorMemory >&); + void vmult(Vector&, const Vector&) const; + void Tvmult(Vector&, const Vector&) const; + void vmult(BlockVector&, + const BlockVector&) const; + void Tvmult(BlockVector&, + const BlockVector&) const; + private: +// SmartPointer > matrix; +// SmartPointer > mem; +}; + + template template @@ -380,19 +450,21 @@ LAPACKFullMatrix::fill ( const unsigned int dst_offset_i, const unsigned int dst_offset_j, const unsigned int src_offset_i, - const unsigned int src_offset_j) + const unsigned int src_offset_j, + const number factor, + const bool transpose) { const typename MATRIX::const_iterator end = M.end(); for (typename MATRIX::const_iterator entry = M.begin(src_offset_i); entry != end; ++entry) { - const unsigned int i = entry->row(); - const unsigned int j = entry->column(); - + const unsigned int i = transpose ? entry->column() : entry->row(); + const unsigned int j = transpose ? entry->row() : entry->column(); + const unsigned int dst_i=dst_offset_i+i-src_offset_i; const unsigned int dst_j=dst_offset_j+j-src_offset_j; if (dst_in_rows() && dst_jn_cols()) - (*this)(dst_i, dst_j) = entry->value(); + (*this)(dst_i, dst_j) = factor * entry->value(); } state = LAPACKSupport::matrix; diff --git a/deal.II/lac/include/lac/lapack_templates.h b/deal.II/lac/include/lac/lapack_templates.h index d6803b9eb0..3d772cbbc1 100644 --- a/deal.II/lac/include/lac/lapack_templates.h +++ b/deal.II/lac/include/lac/lapack_templates.h @@ -31,6 +31,18 @@ void sgemv_ (const char* trans, const int* m, const int* n, const float* alpha, const float* A, const int* lda, const float* x, const int* incx, const float* b, float* y, const int* incy); +// Compute LU factorization +void dgetrf_ (const int* m, const int* n, double* A, + const int* lda, int* ipiv, int* info); +void sgetrf_ (const int* m, const int* n, float* A, + const int* lda, int* ipiv, int* info); +// Apply forward/backward substitution to LU factorization +void dgetrs_ (const char* trans, const int* n, const int* nrhs, + const double* A, const int* lda, const int* ipiv, + double* b, const int* ldb, int* info); +void sgetrs_ (const char* trans, const int* n, const int* nrhs, + const float* A, const int* lda, const int* ipiv, + float* b, const int* ldb, int* info); // Compute eigenvalues and vectors void dgeev_ (const char* jobvl, const char* jobvr, const int* n, double* A, const int* lda, @@ -121,13 +133,73 @@ gemv (const char* trans, const int* m, const int* n, const float* alpha, const f } #else inline void -gemv (const char*, const int*, const int*, const float*, const float*, const int*, const float*, const int*, const float*, float*, const int*) +gemv (const char*, const int*, const int*, const double*, const double*, const int*, const double*, const int*, const double*, double*, const int*) { LAPACKSupport::ExcMissing("sgemv"); } #endif +#ifdef HAVE_DGETRF_ +inline void +getrf (const int* m, const int* n, double* A, const int* lda, int* ipiv, int* info) +{ + dgetrf_ (m,n,A,lda,ipiv,info); +} +#else +inline void +getrf (const int*, const int*, double*, const int*, int*, int*) +{ + LAPACKSupport::ExcMissing("dgetrf"); +} +#endif + + +#ifdef HAVE_SGETRF_ +inline void +getrf (const int* m, const int* n, float* A, const int* lda, int* ipiv, int* info) +{ + sgetrf_ (m,n,A,lda,ipiv,info); +} +#else +inline void +getrf (const int*, const int*, double*, const int*, int*, int*) +{ + LAPACKSupport::ExcMissing("sgetrf"); +} +#endif + + +#ifdef HAVE_DGETRS_ +inline void +getrs (const char* trans, const int* n, const int* nrhs, const double* A, const int* lda, const int* ipiv, double* b, const int* ldb, int* info) +{ + dgetrs_ (trans,n,nrhs,A,lda,ipiv,b,ldb,info); +} +#else +inline void +getrs (const char*, const int*, const int*, const double*, const int*, const int*, double*, const int*, int*) +{ + LAPACKSupport::ExcMissing("dgetrs"); +} +#endif + + +#ifdef HAVE_SGETRS_ +inline void +getrs (const char* trans, const int* n, const int* nrhs, const float* A, const int* lda, const int* ipiv, float* b, const int* ldb, int* info) +{ + sgetrs_ (trans,n,nrhs,A,lda,ipiv,b,ldb,info); +} +#else +inline void +getrs (const char*, const int*, const int*, const double*, const int*, const int*, double*, const int*, int*) +{ + LAPACKSupport::ExcMissing("sgetrs"); +} +#endif + + #ifdef HAVE_DGEEV_ inline void geev (const char* jobvl, const char* jobvr, const int* n, double* A, const int* lda, double* lambda_re, double* lambda_im, double* vl, const int* ldvl, double* vr, const int* ldva, double* work, const int* lwork, int* info) @@ -151,7 +223,7 @@ geev (const char* jobvl, const char* jobvr, const int* n, float* A, const int* l } #else inline void -geev (const char*, const char*, const int*, float*, const int*, float*, float*, float*, const int*, float*, const int*, float*, const int*, int*) +geev (const char*, const char*, const int*, double*, const int*, double*, double*, double*, const int*, double*, const int*, double*, const int*, int*) { LAPACKSupport::ExcMissing("sgeev"); } @@ -181,7 +253,7 @@ geevx (const char* balanc, const char* jobvl, const char* jobvr, const char* sen } #else inline void -geevx (const char*, const char*, const char*, const char*, const int*, float*, const int*, float*, float*, float*, const int*, float*, const int*, int*, int*, float*, float*, float*, float*, float*, const int*, int*, int*) +geevx (const char*, const char*, const char*, const char*, const int*, double*, const int*, double*, double*, double*, const int*, double*, const int*, int*, int*, double*, double*, double*, double*, double*, const int*, int*, int*) { LAPACKSupport::ExcMissing("sgeevx"); } @@ -211,7 +283,7 @@ gesvd (int* jobu, int* jobvt, const int* n, const int* m, float* A, const int* l } #else inline void -gesvd (int*, int*, const int*, const int*, float*, const int*, float*, float*, const int*, float*, const int*, float*, const int*, int*) +gesvd (int*, int*, const int*, const int*, double*, const int*, double*, double*, const int*, double*, const int*, double*, const int*, int*) { LAPACKSupport::ExcMissing("sgesvd"); } @@ -241,7 +313,7 @@ stev (const char* jobz, const int* n, float* d, float* e, float* z, const int* l } #else inline void -stev (const char*, const int*, float*, float*, float*, const int*, float*, int*) +stev (const char*, const int*, double*, double*, double*, const int*, double*, int*) { LAPACKSupport::ExcMissing("sstev"); } diff --git a/deal.II/lac/source/lapack_full_matrix.cc b/deal.II/lac/source/lapack_full_matrix.cc index 599bd43d55..e65f1440cb 100644 --- a/deal.II/lac/source/lapack_full_matrix.cc +++ b/deal.II/lac/source/lapack_full_matrix.cc @@ -25,7 +25,8 @@ using namespace LAPACKSupport; template LAPACKFullMatrix::LAPACKFullMatrix(const unsigned int n) : - TransposeTable (n,n) + TransposeTable (n,n), + state(matrix) {} @@ -34,17 +35,17 @@ LAPACKFullMatrix::LAPACKFullMatrix( const unsigned int m, const unsigned int n) : - TransposeTable (m,n) + TransposeTable (m,n), + state(matrix) {} template LAPACKFullMatrix::LAPACKFullMatrix(const LAPACKFullMatrix &M) : - TransposeTable (M) -{ - state = LAPACKSupport::matrix; -} + TransposeTable (M), + state(matrix) +{} template @@ -95,7 +96,7 @@ LAPACKFullMatrix::vmult ( const Vector &v, const bool adding) const { - Assert (state == matrix, ExcInvalidState()); + Assert (state == matrix, ExcState(state)); const int mm = this->n_rows(); const int nn = this->n_cols(); @@ -113,7 +114,7 @@ LAPACKFullMatrix::Tvmult ( const Vector &v, const bool adding) const { - Assert (state == matrix, ExcInvalidState()); + Assert (state == matrix, ExcState(state)); const int mm = this->n_rows(); const int nn = this->n_cols(); @@ -152,12 +153,53 @@ LAPACKFullMatrix::Tvmult ( #ifdef HAVE_LIBLAPACK +template +void +LAPACKFullMatrix::compute_lu_factorization() +{ + Assert(state == matrix, ExcState(state)); + const int mm = this->n_rows(); + const int nn = this->n_cols(); + number* values = const_cast (this->data()); + ipiv.resize(mm); + int info; + getrf(&mm, &nn, values, &mm, &ipiv[0], &info); + + Assert(info >= 0, ExcInternalError()); + Assert(info == 0, LACExceptions::ExcSingular()); + + state = lu; +} + + +template +void +LAPACKFullMatrix::apply_lu_factorization(Vector& v, + const bool transposed) const +{ + Assert(state == lu, ExcState(state)); + Assert(this->n_rows() == this->n_cols(), + LACExceptions::ExcNotQuadratic()); + + const char* trans = transposed ? &T : &N; + const int nn = this->n_cols(); + const number* values = this->data(); + int info; + + getrs(trans, &nn, &one, values, &nn, &ipiv[0], + v.begin(), &nn, &info); + + Assert(info == 0, ExcInternalError()); +} + + template void LAPACKFullMatrix::compute_eigenvalues( const bool right, const bool left) { + Assert(state == matrix, ExcState(state)); const int nn = this->n_cols(); wr.resize(nn); wi.resize(nn); @@ -232,6 +274,21 @@ LAPACKFullMatrix::compute_eigenvalues( #else +template +LAPACKFullMatrix::compute_lu_factorization() +{ +Assert(false, ExcNeedsLAPACK()); + +} + + +template +LAPACKFullMatrix::apply_lu_factorization(Vector&, bool) +{ + Assert(false, ExcNeedsLAPACK()); +} + + template void LAPACKFullMatrix::compute_eigenvalues(const bool /*right*/, @@ -243,6 +300,7 @@ LAPACKFullMatrix::compute_eigenvalues(const bool /*right*/, #endif + // template // LAPACKFullMatrix::() // {}