From: kronbichler Date: Wed, 19 Nov 2008 13:25:13 +0000 (+0000) Subject: Invert full matrices using Lapack functions if these are detected. Fixed a constructo... X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=9b7feb9e7dca474fd03fc60768c4dd80d4cc3d51;p=dealii-svn.git Invert full matrices using Lapack functions if these are detected. Fixed a constructor problem in Trilinos block vector classes. git-svn-id: https://svn.dealii.org/trunk@17641 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/lac/include/lac/full_matrix.h b/deal.II/lac/include/lac/full_matrix.h index dd89409b5c..ffbebca87c 100644 --- a/deal.II/lac/include/lac/full_matrix.h +++ b/deal.II/lac/include/lac/full_matrix.h @@ -2,7 +2,7 @@ // $Id$ // Version: $Name$ // -// Copyright (C) 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007 by the deal.II authors +// Copyright (C) 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008 by the deal.II authors // // This file is subject to QPL and may not be distributed // without copyright and license information. Please refer @@ -906,7 +906,15 @@ class FullMatrix : public Table<2,number> * well-behaved for positive * definite matrices, but be * aware of round-off errors in - * the indefinite case. + * the indefinite case. + * + * In case deal.II was configured with + * LAPACK, the functions Xgetrf and + * Xgetri build an LU factorization and + * invert the matrix upon that + * factorization, providing best + * performance up to matrices with a + * few hundreds rows and columns. * * The numerical effort to invert * an n x n matrix is of the @@ -1145,6 +1153,26 @@ class FullMatrix : public Table<2,number> //@} friend class Accessor; + + private: + +#ifdef HAVE_LIBLAPACK + /** + * The vector storing the + * permutations applied for + * pivoting in the + * LU-factorization. + */ + std::vector ipiv; + + /** + * Workspace for calculating the + * inverse matrix from an LU + * factorization. + */ + std::vector inv_work; +#endif + }; /**@}*/ diff --git a/deal.II/lac/include/lac/full_matrix.templates.h b/deal.II/lac/include/lac/full_matrix.templates.h index 73af2bc426..08f678ce9a 100644 --- a/deal.II/lac/include/lac/full_matrix.templates.h +++ b/deal.II/lac/include/lac/full_matrix.templates.h @@ -2,7 +2,7 @@ // $Id$ // Version: $Name$ // -// Copyright (C) 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007 by the deal.II authors +// Copyright (C) 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008 by the deal.II authors // // This file is subject to QPL and may not be distributed // without copyright and license information. Please refer diff --git a/deal.II/lac/include/lac/lapack_full_matrix.h b/deal.II/lac/include/lac/lapack_full_matrix.h index 4576f23625..e84ba90d6e 100644 --- a/deal.II/lac/include/lac/lapack_full_matrix.h +++ b/deal.II/lac/include/lac/lapack_full_matrix.h @@ -383,7 +383,7 @@ class LAPACKFullMatrix : public TransposeTable * inverse matrix from an LU * factorization. */ - std::vector iwork; + std::vector inv_work; /** * Real parts of diff --git a/deal.II/lac/include/lac/lapack_templates.h b/deal.II/lac/include/lac/lapack_templates.h index a44da9470e..89fcebf77f 100644 --- a/deal.II/lac/include/lac/lapack_templates.h +++ b/deal.II/lac/include/lac/lapack_templates.h @@ -51,9 +51,9 @@ 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); + int* ipiv, double* inv_work, const int* lwork, int* info); void sgetri_ (const int* n, float* A, const int* lda, - int* ipiv, float* iwork, const int* lwork, int* info); + int* ipiv, float* inv_work, 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, @@ -218,32 +218,32 @@ getrf (const int*, const int*, float*, const int*, int*, int*) #endif -#ifdef HAVE_DGETRI_ +#ifdef HAVE_DGETRF_ inline void -getri (const int* n, double* A, const int* lda, int* ipiv, double* iwork, const int* lwork, int* info) +getri (const int* n, double* A, const int* lda, int* ipiv, double* inv_work, const int* lwork, int* info) { - dgetri_ (n,A,lda,ipiv,iwork,lwork,info); + dgetri_ (n,A,lda,ipiv,inv_work,lwork,info); } #else inline void getri (const int*, double*, const int*, int*, double*, const int*, int*) { - LAPACKSupport::ExcMissing("dgetrf"); + LAPACKSupport::ExcMissing("dgetri"); } #endif -#ifdef HAVE_SGETRI_ +#ifdef HAVE_SGETRF_ inline void -getri (const int* n, float* A, const int* lda, int* ipiv, float* iwork, const int* lwork, int* info) +getri (const int* n, float* A, const int* lda, int* ipiv, float* inv_work, const int* lwork, int* info) { - sgetri_ (n,A,lda,ipiv,iwork,lwork,info); + sgetri_ (n,A,lda,ipiv,inv_work,lwork,info); } #else inline void getri (const int*, float*, const int*, int*, float*, const int*, int*) { - LAPACKSupport::ExcMissing("sgetrf"); + LAPACKSupport::ExcMissing("sgetri"); } #endif diff --git a/deal.II/lac/include/lac/trilinos_block_vector.h b/deal.II/lac/include/lac/trilinos_block_vector.h index 34534dbb59..ce144e0079 100644 --- a/deal.II/lac/include/lac/trilinos_block_vector.h +++ b/deal.II/lac/include/lac/trilinos_block_vector.h @@ -836,7 +836,11 @@ namespace TrilinosWrappers : BlockVectorBase () { - reinit(v); + this->components.resize (v.n_blocks()); + this->block_indices = v.block_indices; + + for (unsigned int i=0; in_blocks(); ++i) + this->components[i] = v.components[i]; } diff --git a/deal.II/lac/source/full_matrix.cc b/deal.II/lac/source/full_matrix.cc index 0fa3339cdb..e5c70d1090 100644 --- a/deal.II/lac/source/full_matrix.cc +++ b/deal.II/lac/source/full_matrix.cc @@ -13,10 +13,135 @@ #include +#include #include DEAL_II_NAMESPACE_OPEN + + // Need to explicitly state the Lapack + // inversion since it only works with + // floats and doubles in case LAPACK was + // detected by configure. +#ifdef HAVE_LIBLAPACK + +template <> +void +FullMatrix::gauss_jordan () +{ + Assert (!this->empty(), ExcEmptyMatrix()); + Assert (this->n_cols() == this->n_rows(), ExcNotQuadratic()); + + // In case we have the LAPACK functions + // getrf and getri detected at configure, + // we use these algorithms for inversion + // since they provide better performance + // than the deal.II native functions. + // + // Note that BLAS/LAPACK stores matrix + // elements column-wise (i.e., all values in + // one column, then all in the next, etc.), + // whereas the FullMatrix stores them + // row-wise. + // We ignore that difference, and give our + // row-wise data to LAPACK, + // let LAPACK build the inverse of the + // transpose matrix, and read the result as + // if it were row-wise again. In other words, + // we just got ((A^T)^{-1})^T, which is + // A^{-1}. + + const int nn = this->n(); + float* values = const_cast (this->data()); + ipiv.resize(nn); + int info; + + // Use the LAPACK function getrf for + // calculating the LU factorization. + getrf(&nn, &nn, values, &nn, &ipiv[0], &info); + + Assert(info >= 0, ExcInternalError()); + Assert(info == 0, LACExceptions::ExcSingular()); + + inv_work.resize (nn); + // Use the LAPACK function getri for + // calculating the actual inverse using + // the LU factorization. + getri(&nn, values, &nn, &ipiv[0], &inv_work[0], &nn, &info); + + Assert(info >= 0, ExcInternalError()); + Assert(info == 0, LACExceptions::ExcSingular()); +} + +template <> +void +FullMatrix::gauss_jordan () +{ + Assert (!this->empty(), ExcEmptyMatrix()); + Assert (this->n_cols() == this->n_rows(), ExcNotQuadratic()); + + // In case we have the LAPACK functions + // getrf and getri detected at configure, + // we use these algorithms for inversion + // since they provide better performance + // than the deal.II native functions. + // + // Note that BLAS/LAPACK stores matrix + // elements column-wise (i.e., all values in + // one column, then all in the next, etc.), + // whereas the FullMatrix stores them + // row-wise. + // We ignore that difference, and give our + // row-wise data to LAPACK, + // let LAPACK build the inverse of the + // transpose matrix, and read the result as + // if it were row-wise again. In other words, + // we just got ((A^T)^{-1})^T, which is + // A^{-1}. + + const int nn = this->n(); + double* values = const_cast (this->data()); + ipiv.resize(nn); + int info; + + // Use the LAPACK function getrf for + // calculating the LU factorization. + getrf(&nn, &nn, values, &nn, &ipiv[0], &info); + + Assert(info >= 0, ExcInternalError()); + Assert(info == 0, LACExceptions::ExcSingular()); + + inv_work.resize (nn); + // Use the LAPACK function getri for + // calculating the actual inverse using + // the LU factorization. + std::cout << " blas"; + getri(&nn, values, &nn, &ipiv[0], &inv_work[0], &nn, &info); + std::cout << " 2 mal" << std::endl; + + Assert(info >= 0, ExcInternalError()); + Assert(info == 0, LACExceptions::ExcSingular()); +} + + // ... and now the usual instantiations + // of gauss_jordan() and all the rest. +template void FullMatrix::gauss_jordan (); +template void FullMatrix >::gauss_jordan (); +template void FullMatrix >::gauss_jordan (); +template void FullMatrix >::gauss_jordan (); + +#else + +template void FullMatrix::gauss_jordan (); +template void FullMatrix::gauss_jordan (); +template void FullMatrix::gauss_jordan (); +template void FullMatrix >::gauss_jordan (); +template void FullMatrix >::gauss_jordan (); +template void FullMatrix >::gauss_jordan (); + +#endif + + #include "full_matrix.inst" @@ -50,4 +175,5 @@ TEMPL_OP_EQ(std::complex,std::complex); #undef TEMPL_OP_EQ + DEAL_II_NAMESPACE_CLOSE diff --git a/deal.II/lac/source/lapack_full_matrix.cc b/deal.II/lac/source/lapack_full_matrix.cc index 20e8f44230..386b85703f 100644 --- a/deal.II/lac/source/lapack_full_matrix.cc +++ b/deal.II/lac/source/lapack_full_matrix.cc @@ -197,8 +197,8 @@ LAPACKFullMatrix::invert() Assert(info == 0, LACExceptions::ExcSingular()); } - iwork.resize (mm); - getri(&mm, values, &mm, &ipiv[0], &iwork[0], &mm, &info); + inv_work.resize (mm); + getri(&mm, values, &mm, &ipiv[0], &inv_work[0], &mm, &info); Assert(info >= 0, ExcInternalError()); Assert(info == 0, LACExceptions::ExcSingular()); diff --git a/deal.II/lac/source/trilinos_block_vector.cc b/deal.II/lac/source/trilinos_block_vector.cc index 0d64f03947..369b99ecc9 100644 --- a/deal.II/lac/source/trilinos_block_vector.cc +++ b/deal.II/lac/source/trilinos_block_vector.cc @@ -231,9 +231,7 @@ namespace TrilinosWrappers components.resize(n_blocks()); for (unsigned int i=0;i