From 277cecf26510a667ba08d3823e77bc612e3d72eb Mon Sep 17 00:00:00 2001 From: Guido Kanschat Date: Sun, 7 Nov 2010 23:11:12 +0000 Subject: [PATCH] implement SVD git-svn-id: https://svn.dealii.org/trunk@22628 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/doc/news/changes.h | 6 + deal.II/include/deal.II/lac/full_matrix.h | 5 +- .../deal.II/lac/full_matrix.templates.h | 1 + .../include/deal.II/lac/lapack_full_matrix.h | 97 +++++++- deal.II/include/deal.II/lac/lapack_support.h | 25 +- .../include/deal.II/lac/lapack_templates.h.in | 10 + .../deal.II/lac/sparse_matrix.templates.h | 2 +- deal.II/source/lac/full_matrix.cc | 4 +- deal.II/source/lac/lapack_full_matrix.cc | 234 ++++++++++++------ 9 files changed, 294 insertions(+), 90 deletions(-) diff --git a/deal.II/doc/news/changes.h b/deal.II/doc/news/changes.h index c24b26546b..db4b4f91b5 100644 --- a/deal.II/doc/news/changes.h +++ b/deal.II/doc/news/changes.h @@ -296,6 +296,12 @@ through DoFHandler::get_tria() and DoFHandler::get_fe(), respectively.
    +

  1. New: The class LAPACKFullMatrix now has functions to compute the + singular value decomposition of a matrix and its inverse. +
    + (GK 2010/11/7) +

    2. +

  2. New: The classes RelaxationBlockSOR and RelaxationBlockSSOR implement overlapping Schwarz relaxation methods. Additionally, their base class RelaxationBlock and the helper class BlockList have diff --git a/deal.II/include/deal.II/lac/full_matrix.h b/deal.II/include/deal.II/lac/full_matrix.h index a62e83af85..5195577946 100644 --- a/deal.II/include/deal.II/lac/full_matrix.h +++ b/deal.II/include/deal.II/lac/full_matrix.h @@ -19,18 +19,17 @@ #include #include #include -#include - #include #include +#include DEAL_II_NAMESPACE_OPEN // forward declarations template class Vector; - +template class LAPACKFullMatrix; template class Tensor; diff --git a/deal.II/include/deal.II/lac/full_matrix.templates.h b/deal.II/include/deal.II/lac/full_matrix.templates.h index c45f2763f6..7db3267f9f 100644 --- a/deal.II/include/deal.II/lac/full_matrix.templates.h +++ b/deal.II/include/deal.II/lac/full_matrix.templates.h @@ -20,6 +20,7 @@ #include #include #include +#include #include #include diff --git a/deal.II/include/deal.II/lac/lapack_full_matrix.h b/deal.II/include/deal.II/lac/lapack_full_matrix.h index 52506732b0..d80c651a2d 100644 --- a/deal.II/include/deal.II/lac/lapack_full_matrix.h +++ b/deal.II/include/deal.II/lac/lapack_full_matrix.h @@ -341,6 +341,60 @@ class LAPACKFullMatrix : public TransposeTable std::vector > & eigenvectors, const int itype = 1); + /** + * Compute the singular value + * decomposition of the + * matrix using LAPACK function + * Xgesdd. + * + * Requires that the #state is + * LAPACKSupport::matrix, fills + * the data members #wr, #svd_u, + * and #svd_vt, and leaves the + * object in the #state + * LAPACKSupport::svd. + */ + void compute_svd(); + /** + * Compute the inverse of the + * matrix by singular value + * decomposition. + * + * Requires that #state is either + * LAPACKSupport::matrix or + * LAPACKSupport::svd. In the + * first case, this function + * calls compute_svd(). After + * this function, the object will + * have the #state + * LAPACKSupport::inverse_svd. + * + * For a singular value + * decomposition, the inverse is + * simply computed by replacing + * all singular values by their + * reciprocal values. If the + * matrix does not have maximal + * rank, singular values 0 are + * not touched, thus computing + * the minimal norm right inverse + * of the matrix. + * + * The parameter + * threshold determines, + * when a singular value should + * be considered zero. It is the + * ratio of the smallest to the + * largest nonzero singular + * value + * smax. Thus, + * the inverses of all singular + * values less than + * smax/threshold + * will be set to zero. + */ + void compute_inverse_svd (const double threshold = 0.); + /** * Retrieve eigenvalue after * compute_eigenvalues() was @@ -349,6 +403,15 @@ class LAPACKFullMatrix : public TransposeTable std::complex eigenvalue (const unsigned int i) const; + /** + * Retrieve singular values after + * compute_svd() or + * compute_inverse_svd() was + * called. + */ + number + singular_value (const unsigned int i) const; + /** * Print the matrix and allow * formatting of entries. @@ -426,6 +489,10 @@ class LAPACKFullMatrix : public TransposeTable * permutations applied for * pivoting in the * LU-factorization. + * + * Also used as the scratch array + * IWORK for LAPACK functions + * needing it. */ std::vector ipiv; @@ -462,15 +529,18 @@ class LAPACKFullMatrix : public TransposeTable */ std::vector vr; - /** - * The matrix U in the singular value decomposition - * USVT. - */ + /** + * The matrix U in the + * singular value decomposition + * USVT. + */ boost::shared_ptr > svd_u; - /** - * The matrix VT in the singular value decomposition - * USVT. - */ + /** + * The matrix + * VT in the + * singular value decomposition + * USVT. + */ boost::shared_ptr > svd_vt; }; @@ -562,6 +632,17 @@ LAPACKFullMatrix::eigenvalue (const unsigned int i) const } +template +number +LAPACKFullMatrix::singular_value (const unsigned int i) const +{ + Assert (state == LAPACKSupport::svd || state == LAPACKSupport::inverse_svd, LAPACKSupport::ExcState(state)); + AssertIndexRange(i,wr.size()); + + return wr[i]; +} + + DEAL_II_NAMESPACE_CLOSE diff --git a/deal.II/include/deal.II/lac/lapack_support.h b/deal.II/include/deal.II/lac/lapack_support.h index 6aa4b9ee05..34de869a10 100644 --- a/deal.II/include/deal.II/lac/lapack_support.h +++ b/deal.II/include/deal.II/lac/lapack_support.h @@ -2,7 +2,7 @@ // $Id$ // Version: $Name$ // -// Copyright (C) 2005, 2006, 2008 by the deal.II authors +// Copyright (C) 2005, 2006, 2008, 2010 by the deal.II authors // // This file is subject to QPL and may not be distributed // without copyright and license information. Please refer @@ -42,6 +42,8 @@ namespace LAPACKSupport eigenvalues, /// Matrix contains singular value decomposition, svd, + /// Matrix is the inverse of a singular value decomposition + inverse_svd, /// Contents is something useless. unusable = 0x8000 }; @@ -61,6 +63,10 @@ namespace LAPACKSupport return "lu decomposition"; case eigenvalues: return "eigenvalues"; + case svd: + return "svd"; + case inverse_svd: + return "inverse_svd"; case unusable: return "unusable"; default: @@ -92,19 +98,23 @@ namespace LAPACKSupport /** * Character constant. */ - static const char V = 'V'; + static const char A = 'A'; /** * Character constant. */ - static const char T = 'T'; + static const char N = 'N'; /** * Character constant. */ - static const char N = 'N'; + static const char T = 'T'; /** * Character constant. */ static const char U = 'U'; + /** + * Character constant. + */ + static const char V = 'V'; /** * Integer constant. */ @@ -114,6 +124,13 @@ namespace LAPACKSupport */ static const int one = 1; + /** + * A LAPACK function returned an + * error code. + */ + DeclException2(ExcErrorCode, char*, int, + << "The function " << arg1 << " returned with an error code " << arg2); + /** * Exception thrown when a matrix * is not in a suitable state for diff --git a/deal.II/include/deal.II/lac/lapack_templates.h.in b/deal.II/include/deal.II/lac/lapack_templates.h.in index 631796654d..76e78e7216 100644 --- a/deal.II/include/deal.II/lac/lapack_templates.h.in +++ b/deal.II/include/deal.II/lac/lapack_templates.h.in @@ -78,6 +78,16 @@ void dsygv_ (const int* itype, const char* jobz, const char* uplo, const int* ldb, double* w, double* work, const int* lwork, int* info); +// Compute singular value decomposition using divide and conquer +void dgesdd_ (const char* jobz, + const int* m, const int* n, double* A, const int* lda, + double* s, + double* u, const int* ldu, + double* vt, const int* ldvt, + double* work, const int* lwork, + int* iwork, + int* info); + // Compute singular value decomposition void dgesvd_ (int* jobu, int* jobvt, const int* n, const int* m, double* A, const int* lda, diff --git a/deal.II/include/deal.II/lac/sparse_matrix.templates.h b/deal.II/include/deal.II/lac/sparse_matrix.templates.h index f107bdc6c9..24b6b8dd1d 100644 --- a/deal.II/include/deal.II/lac/sparse_matrix.templates.h +++ b/deal.II/include/deal.II/lac/sparse_matrix.templates.h @@ -24,7 +24,7 @@ #include #include #include - +#include // we only need output streams, but older compilers did not provide // them in a separate include file diff --git a/deal.II/source/lac/full_matrix.cc b/deal.II/source/lac/full_matrix.cc index 9af29628ec..f6c9a2e2c7 100644 --- a/deal.II/source/lac/full_matrix.cc +++ b/deal.II/source/lac/full_matrix.cc @@ -2,7 +2,7 @@ // $Id$ // Version: $Name$ // -// Copyright (C) 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008 by the deal.II authors +// Copyright (C) 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2010 by the deal.II authors // // This file is subject to QPL and may not be distributed // without copyright and license information. Please refer @@ -12,10 +12,10 @@ //--------------------------------------------------------------------------- -#include #include #include #include +#include DEAL_II_NAMESPACE_OPEN diff --git a/deal.II/source/lac/lapack_full_matrix.cc b/deal.II/source/lac/lapack_full_matrix.cc index f051a05374..407a63e014 100644 --- a/deal.II/source/lac/lapack_full_matrix.cc +++ b/deal.II/source/lac/lapack_full_matrix.cc @@ -90,7 +90,6 @@ LAPACKFullMatrix::operator = (const double d) return *this; } -#ifdef HAVE_LIBBLAS template void @@ -98,16 +97,55 @@ LAPACKFullMatrix::vmult ( Vector &w, const Vector &v, const bool adding) const -{ - Assert ((state == matrix) || (state == inverse_matrix), - ExcState(state)); - +{ const int mm = this->n_rows(); const int nn = this->n_cols(); const number alpha = 1.; const number beta = (adding ? 1. : 0.); - - gemv("N", &mm, &nn, &alpha, this->data(), &mm, v.val, &one, &beta, w.val, &one); + const number null = 0.; + + switch (state) + { + case matrix: + case inverse_matrix: + { + AssertDimension(v.size(), this->n_cols()); + AssertDimension(w.size(), this->n_rows()); + + gemv("N", &mm, &nn, &alpha, this->data(), &mm, v.val, &one, &beta, w.val, &one); + break; + } + case svd: + { + AssertDimension(v.size(), this->n_cols()); + AssertDimension(w.size(), this->n_rows()); + // Compute V^T v + work.resize(std::max(mm,nn)); + gemv("N", &nn, &nn, &alpha, svd_vt->data(), &nn, v.val, &one, &null, &work[0], &one); + // Multiply by singular values + for (unsigned int i=0;idata(), &mm, &work[0], &one, &beta, w.val, &one); + break; + } + case inverse_svd: + { + AssertDimension(w.size(), this->n_cols()); + AssertDimension(v.size(), this->n_rows()); + // Compute U^T v + work.resize(std::max(mm,nn)); + gemv("T", &mm, &mm, &alpha, svd_u->data(), &mm, v.val, &one, &null, &work[0], &one); + // Multiply by singular values + for (unsigned int i=0;idata(), &nn, &work[0], &one, &beta, w.val, &one); + break; + } + default: + Assert (false, ExcState(state)); + } } @@ -118,61 +156,148 @@ LAPACKFullMatrix::Tvmult ( const Vector &v, const bool adding) const { - Assert (state == matrix, ExcState(state)); - const int mm = this->n_rows(); const int nn = this->n_cols(); const number alpha = 1.; const number beta = (adding ? 1. : 0.); - - gemv("T", &mm, &nn, &alpha, this->data(), &mm, v.val, &one, &beta, w.val, &one); + const number null = 0.; + + switch (state) + { + case matrix: + case inverse_matrix: + { + AssertDimension(w.size(), this->n_cols()); + AssertDimension(v.size(), this->n_rows()); + + gemv("T", &mm, &nn, &alpha, this->data(), &mm, v.val, &one, &beta, w.val, &one); + break; + } + case svd: + { + AssertDimension(w.size(), this->n_cols()); + AssertDimension(v.size(), this->n_rows()); + + // Compute U^T v + work.resize(std::max(mm,nn)); + gemv("T", &mm, &mm, &alpha, svd_u->data(), &mm, v.val, &one, &null, &work[0], &one); + // Multiply by singular values + for (unsigned int i=0;idata(), &nn, &work[0], &one, &beta, w.val, &one); + break; + case inverse_svd: + { + AssertDimension(v.size(), this->n_cols()); + AssertDimension(w.size(), this->n_rows()); + + // Compute V^T v + work.resize(std::max(mm,nn)); + gemv("N", &nn, &nn, &alpha, svd_vt->data(), &nn, v.val, &one, &null, &work[0], &one); + // Multiply by singular values + for (unsigned int i=0;idata(), &mm, &work[0], &one, &beta, w.val, &one); + break; + } + } + default: + Assert (false, ExcState(state)); + } } -#else - template void -LAPACKFullMatrix::vmult ( - Vector &, - const Vector &, - const bool ) const +LAPACKFullMatrix::compute_lu_factorization() { - Assert(false, ExcNeedsBLAS()); -} + 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()); -template -void -LAPACKFullMatrix::Tvmult ( - Vector &, - const Vector &, - const bool ) const -{ - Assert(false, ExcNeedsBLAS()); + state = lu; } -#endif - -#ifdef HAVE_LIBLAPACK - template void -LAPACKFullMatrix::compute_lu_factorization() +LAPACKFullMatrix::compute_svd() { Assert(state == matrix, ExcState(state)); + state = LAPACKSupport::unusable; + const int mm = this->n_rows(); const int nn = this->n_cols(); number* values = const_cast (this->data()); - ipiv.resize(mm); + wr.resize(std::max(mm,nn)); + std::fill(wr.begin(), wr.end(), 0.); + ipiv.resize(8*mm); + + svd_u = boost::shared_ptr >(new LAPACKFullMatrix(mm,mm)); + svd_vt = boost::shared_ptr >(new LAPACKFullMatrix(nn,nn)); + number* mu = const_cast (svd_u->data()); + number* mvt = const_cast (svd_vt->data()); int info; - getrf(&mm, &nn, values, &mm, &ipiv[0], &info); - Assert(info >= 0, ExcInternalError()); - Assert(info == 0, LACExceptions::ExcSingular()); + // see comment on this #if + // below. Another reason to love Petsc +#ifndef DEAL_II_LIBLAPACK_NOQUERYMODE - state = lu; + // First determine optimal + // workspace size + work.resize(1); + int lwork = -1; + gesdd(&LAPACKSupport::A, &mm, &nn, values, &mm, + &wr[0], mu, &mm, mvt, &nn, + &work[0], &lwork, &ipiv[0], &info); + Assert (info==0, LAPACKSupport::ExcErrorCode("gesdd", info)); + // Now resize work array and + lwork = work[0] + .5; +#else + int lwork = 3*std::min(mm,nn) + + std::max(std::max(mm,nn),4*std::min(mm,nn)*std::min(mm,nn)+4*std::min(mm,nn)); +#endif + work.resize(lwork); + // Do the actual SVD. + gesdd(&LAPACKSupport::A, &mm, &nn, values, &mm, + &wr[0], mu, &mm, mvt, &nn, + &work[0], &lwork, &ipiv[0], &info); + Assert (info==0, LAPACKSupport::ExcErrorCode("gesdd", info)); + + work.resize(0); + ipiv.resize(0); + + state = LAPACKSupport::svd; +} + + +template +void +LAPACKFullMatrix::compute_inverse_svd(const double threshold) +{ + if (state == LAPACKSupport::matrix) + compute_svd(); + + Assert (state==LAPACKSupport::svd, ExcState(state)); + + const double lim = wr[0]*threshold; + for (unsigned int i=0;i lim) + wr[i] = 1./wr[i]; + else + wr[i] = 0.; + } + state = LAPACKSupport::inverse_svd; } @@ -403,41 +528,6 @@ LAPACKFullMatrix::compute_generalized_eigenvalues_symmetric ( state = LAPACKSupport::State(eigenvalues | unusable); } -#else - -template -void -LAPACKFullMatrix::compute_lu_factorization() -{ -Assert(false, ExcNeedsLAPACK()); - -} - - -template -void -LAPACKFullMatrix::apply_lu_factorization(Vector&, bool) const -{ - Assert(false, ExcNeedsLAPACK()); -} - - -template -void -LAPACKFullMatrix::compute_eigenvalues(const bool /*right*/, - const bool /*left*/) -{ - Assert(false, ExcNeedsLAPACK()); -} - - -#endif - - -// template -// LAPACKFullMatrix::() -// {} - // template // LAPACKFullMatrix::() -- 2.39.5