From: kronbichler Date: Tue, 11 Nov 2008 18:49:50 +0000 (+0000) Subject: Added the functionality to calculate the inverse in LAPACKFullMatrix class. X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=d8263a85763511fbf512254845ec21f33fc28af2;p=dealii-svn.git Added the functionality to calculate the inverse in LAPACKFullMatrix class. git-svn-id: https://svn.dealii.org/trunk@17547 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/configure b/deal.II/configure index 33b7a60eaa..8adaaed733 100755 --- a/deal.II/configure +++ b/deal.II/configure @@ -15420,7 +15420,7 @@ done -for ac_func in dgetrf_ sgetrf_ dgetrs_ sgetrs_ dstev_ sstev_ +for ac_func in dgetrf_ sgetrf_ dgetri_ sgetri_ 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 fd9b758729..03371f18f1 100644 --- a/deal.II/configure.in +++ b/deal.II/configure.in @@ -542,7 +542,7 @@ if test "x$with_blas" != "x" -a "x$with_blas" != "xno" ; then fi AC_CHECK_FUNCS([dgemv_ sgemv_ dgeev_ sgeev_ dgeevx_ sgeevx_ dgesvd_ sgesvd_]) -AC_CHECK_FUNCS([dgetrf_ sgetrf_ dgetrs_ sgetrs_ dstev_ sstev_]) +AC_CHECK_FUNCS([dgetrf_ sgetrf_ dgetri_ sgetri_ dgetrs_ sgetrs_ dstev_ sstev_]) dnl ------------------------------------------------------------- dnl set include paths of several libraries diff --git a/deal.II/doc/news/changes.h b/deal.II/doc/news/changes.h index 5586bff1fe..844155aa66 100644 --- a/deal.II/doc/news/changes.h +++ b/deal.II/doc/news/changes.h @@ -322,6 +322,16 @@ inconvenience this causes.

lac

    +
  1. +

    + New: The class LAPACKFullMatrix can now invert full matrices using + the (optimized) LAPACK functions getrf and getri. The speedup over + the FullMatrix::gauss_jordan() function is a factor of two for matrices + with 100 rows and columns, and grows with matrix size. +
    + (Martin Kronbichler 2008/11/11) +

    +
  2. Fixed: The BlockMatrixBase::clear() function that is used by all other diff --git a/deal.II/lac/include/lac/full_matrix.h b/deal.II/lac/include/lac/full_matrix.h index 976398c612..dd89409b5c 100644 --- a/deal.II/lac/include/lac/full_matrix.h +++ b/deal.II/lac/include/lac/full_matrix.h @@ -19,6 +19,7 @@ #include #include #include +#include #include #include @@ -330,6 +331,16 @@ class FullMatrix : public Table<2,number> FullMatrix & operator = (const IdentityMatrix &id); + /** + * Assignment operator for a + * LapackFullMatrix. The calling matrix + * must be of the same size as the + * LAPACK matrix. + */ + template + FullMatrix & + operator = (const LAPACKFullMatrix&); + /** * Assignment from different * matrix classes. This diff --git a/deal.II/lac/include/lac/full_matrix.templates.h b/deal.II/lac/include/lac/full_matrix.templates.h index aafc65b128..73af2bc426 100644 --- a/deal.II/lac/include/lac/full_matrix.templates.h +++ b/deal.II/lac/include/lac/full_matrix.templates.h @@ -108,6 +108,22 @@ FullMatrix::operator = (const IdentityMatrix &id) +template +template +FullMatrix & +FullMatrix::operator = (const LAPACKFullMatrix& M) +{ + Assert (this->m() == M.n_rows(), ExcDimensionMismatch(this->m(), M.n_rows())); + Assert (this->n() == M.n_cols(), ExcDimensionMismatch(this->n(), M.n_rows())); + for (unsigned int i=0;im();++i) + for (unsigned int j=0;jn();++j) + (*this)(i,j) = M(i,j); + + return *this; +} + + + template bool FullMatrix::all_zero () const diff --git a/deal.II/lac/include/lac/lapack_full_matrix.h b/deal.II/lac/include/lac/lapack_full_matrix.h index 0f238cef55..4576f23625 100644 --- a/deal.II/lac/include/lac/lapack_full_matrix.h +++ b/deal.II/lac/include/lac/lapack_full_matrix.h @@ -2,7 +2,7 @@ // $Id$ // Version: $Name$ // -// Copyright (C) 2005, 2006 by the deal.II authors +// Copyright (C) 2005, 2006, 2008 by the deal.II authors // // This file is subject to QPL and may not be distributed // without copyright and license information. Please refer @@ -231,6 +231,14 @@ class LAPACKFullMatrix : public TransposeTable */ void compute_lu_factorization (); + /** + * Invert the matrix by first computing + * an LU factorization with the LAPACK + * function Xgetrf and then building + * the actual inverse using Xgetri. + */ + void invert (); + /** * Solve the linear system with * right hand side given by @@ -370,6 +378,13 @@ class LAPACKFullMatrix : public TransposeTable */ std::vector ipiv; + /** + * Workspace for calculating the + * inverse matrix from an LU + * factorization. + */ + std::vector iwork; + /** * Real parts of * eigenvalues. Filled by diff --git a/deal.II/lac/include/lac/lapack_support.h b/deal.II/lac/include/lac/lapack_support.h index dcfab53d44..034eb9f4e6 100644 --- a/deal.II/lac/include/lac/lapack_support.h +++ b/deal.II/lac/include/lac/lapack_support.h @@ -2,7 +2,7 @@ // $Id$ // Version: $Name$ // -// Copyright (C) 2005, 2006 by the deal.II authors +// Copyright (C) 2005, 2006, 2008 by the deal.II authors // // This file is subject to QPL and may not be distributed // without copyright and license information. Please refer @@ -34,6 +34,8 @@ namespace LAPACKSupport { /// Contents is actually a matrix. matrix, + /// Contents is the inverse of a matrix. + inverse_matrix, /// Contents is an LU decomposition. lu, /// Eigenvalue vector is filled @@ -51,6 +53,8 @@ namespace LAPACKSupport { case matrix: return "matrix"; + case inverse_matrix: + return "inverse matrix"; case lu: return "lu decomposition"; case eigenvalues: diff --git a/deal.II/lac/include/lac/lapack_templates.h b/deal.II/lac/include/lac/lapack_templates.h index cf14636785..a44da9470e 100644 --- a/deal.II/lac/include/lac/lapack_templates.h +++ b/deal.II/lac/include/lac/lapack_templates.h @@ -5,7 +5,7 @@ // This file was automatically generated from blas.h.in // See blastemplates in the deal.II contrib directory // -// Copyright (C) 2005, 2006, 2007 by the deal authors +// Copyright (C) 2005, 2006, 2007, 2008 by the deal authors // // This file is subject to QPL and may not be distributed // without copyright and license information. Please refer @@ -49,6 +49,11 @@ 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); +// Invert matrix from LU factorization +void dgetri_ (const int* n, double* A, const int* lda, + int* ipiv, double* iwork, const int* lwork, int* info); +void sgetri_ (const int* n, float* A, const int* lda, + int* ipiv, float* iwork, const int* lwork, 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, @@ -213,6 +218,36 @@ getrf (const int*, const int*, float*, const int*, int*, int*) #endif +#ifdef HAVE_DGETRI_ +inline void +getri (const int* n, double* A, const int* lda, int* ipiv, double* iwork, const int* lwork, int* info) +{ + dgetri_ (n,A,lda,ipiv,iwork,lwork,info); +} +#else +inline void +getri (const int*, double*, const int*, int*, double*, const int*, int*) +{ + LAPACKSupport::ExcMissing("dgetrf"); +} +#endif + + +#ifdef HAVE_SGETRI_ +inline void +getri (const int* n, float* A, const int* lda, int* ipiv, float* iwork, const int* lwork, int* info) +{ + sgetri_ (n,A,lda,ipiv,iwork,lwork,info); +} +#else +inline void +getri (const int*, float*, const int*, int*, float*, const 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) diff --git a/deal.II/lac/source/full_matrix.inst.in b/deal.II/lac/source/full_matrix.inst.in index e1e0e89a83..4b527ba37d 100644 --- a/deal.II/lac/source/full_matrix.inst.in +++ b/deal.II/lac/source/full_matrix.inst.in @@ -24,6 +24,9 @@ for (S : REAL_SCALARS) for (S1, S2 : REAL_SCALARS) { + template + FullMatrix& FullMatrix::operator = (const LAPACKFullMatrix&); + template void FullMatrix::fill ( const FullMatrix&, unsigned, unsigned, unsigned, unsigned); diff --git a/deal.II/lac/source/lapack_full_matrix.cc b/deal.II/lac/source/lapack_full_matrix.cc index 4a8a01ff1a..20e8f44230 100644 --- a/deal.II/lac/source/lapack_full_matrix.cc +++ b/deal.II/lac/source/lapack_full_matrix.cc @@ -2,7 +2,7 @@ // $Id$ // Version: $Name$ // -// Copyright (C) 2005, 2006 by the deal.II authors +// Copyright (C) 2005, 2006, 2008 by the deal.II authors // // This file is subject to QPL and may not be distributed // without copyright and license information. Please refer @@ -175,6 +175,38 @@ LAPACKFullMatrix::compute_lu_factorization() } +template +void +LAPACKFullMatrix::invert() +{ + Assert(state == matrix || state == lu, + ExcState(state)); + const int mm = this->n_rows(); + const int nn = this->n_cols(); + Assert (nn == mm, ExcNotQuadratic()); + + number* values = const_cast (this->data()); + ipiv.resize(mm); + int info; + + if (state == matrix) + { + getrf(&mm, &nn, values, &mm, &ipiv[0], &info); + + Assert(info >= 0, ExcInternalError()); + Assert(info == 0, LACExceptions::ExcSingular()); + } + + iwork.resize (mm); + getri(&mm, values, &mm, &ipiv[0], &iwork[0], &mm, &info); + + Assert(info >= 0, ExcInternalError()); + Assert(info == 0, LACExceptions::ExcSingular()); + + state = inverse_matrix; +} + + template void LAPACKFullMatrix::apply_lu_factorization(Vector& v,