]> https://gitweb.dealii.org/ - dealii.git/commitdiff
implement SVD
authorGuido Kanschat <dr.guido.kanschat@gmail.com>
Sun, 7 Nov 2010 23:11:12 +0000 (23:11 +0000)
committerGuido Kanschat <dr.guido.kanschat@gmail.com>
Sun, 7 Nov 2010 23:11:12 +0000 (23:11 +0000)
git-svn-id: https://svn.dealii.org/trunk@22628 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/doc/news/changes.h
deal.II/include/deal.II/lac/full_matrix.h
deal.II/include/deal.II/lac/full_matrix.templates.h
deal.II/include/deal.II/lac/lapack_full_matrix.h
deal.II/include/deal.II/lac/lapack_support.h
deal.II/include/deal.II/lac/lapack_templates.h.in
deal.II/include/deal.II/lac/sparse_matrix.templates.h
deal.II/source/lac/full_matrix.cc
deal.II/source/lac/lapack_full_matrix.cc

index c24b26546b77506a7e004403c7dc6db3a23cc0fc..db4b4f91b56a38d3039b9ee9afd1f477307e8a3b 100644 (file)
@@ -296,6 +296,12 @@ through DoFHandler::get_tria() and DoFHandler::get_fe(), respectively.
 
 <ol>
 
+  <li><p>New: The class LAPACKFullMatrix now has functions to compute the
+  singular value decomposition of a matrix and its inverse.
+  <br>
+  (GK 2010/11/7)
+  </p></li>
+  
   <li><p>New: The classes RelaxationBlockSOR and RelaxationBlockSSOR
   implement overlapping Schwarz relaxation methods. Additionally,
   their base class RelaxationBlock and the helper class BlockList have
index a62e83af8562159909728d11a29185b8086cc2cb..5195577946cfd381c88642e4b11995779b1f9e84 100644 (file)
 #include <base/table.h>
 #include <lac/exceptions.h>
 #include <lac/identity_matrix.h>
-#include <lac/lapack_full_matrix.h>
-
 
 #include <vector>
 #include <iomanip>
+#include <cstring>
 
 DEAL_II_NAMESPACE_OPEN
 
 
 // forward declarations
 template <typename number> class Vector;
-
+template <typename number> class LAPACKFullMatrix;
 template <int rank, int dim> class Tensor;
 
 
index c45f2763f62aeab016c665580f6ce0ad01354b98..7db3267f9f6cf81308d9276889f49b9cfb3ffb46 100644 (file)
@@ -20,6 +20,7 @@
 #include <base/template_constraints.h>
 #include <lac/vector.h>
 #include <lac/full_matrix.h>
+#include <lac/lapack_full_matrix.h>
 #include <lac/lapack_templates.h>
 
 #include <vector>
index 52506732b0ffc0c719909245d026a890ef4b9494..d80c651a2d230427dc8a2385f9ae5d645b4b0e5d 100644 (file)
@@ -341,6 +341,60 @@ class LAPACKFullMatrix : public TransposeTable<number>
                               std::vector<Vector<number> > & 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
+                                     * <tt>threshold</tt> determines,
+                                     * when a singular value should
+                                     * be considered zero. It is the
+                                     * ratio of the smallest to the
+                                     * largest nonzero singular
+                                     * value
+                                     * <i>s</i><sub>max</sub>. Thus,
+                                     * the inverses of all singular
+                                     * values less than
+                                     * <i>s</i><sub>max</sub>/<tt>threshold</tt>
+                                     * 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<number>
     std::complex<number>
     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<number>
                                      * permutations applied for
                                      * pivoting in the
                                      * LU-factorization.
+                                     *
+                                     * Also used as the scratch array
+                                     * IWORK for LAPACK functions
+                                     * needing it.
                                      */
     std::vector<int> ipiv;
 
@@ -462,15 +529,18 @@ class LAPACKFullMatrix : public TransposeTable<number>
                                      */
     std::vector<number> vr;
     
-    /**
-     * The matrix <i>U</i> in the singular value decomposition
-     * <i>USV<sup>T</sup></i>.
-     */
+                                    /**
+                                     * The matrix <i>U</i> in the
+                                     * singular value decomposition
+                                     * <i>USV<sup>T</sup></i>.
+                                     */
     boost::shared_ptr<LAPACKFullMatrix<number> > svd_u;
-    /**
-     * The matrix <i>V<sup>T</sup></i> in the singular value decomposition
-     * <i>USV<sup>T</sup></i>.
-     */
+                                    /**
+                                     * The matrix
+                                     * <i>V<sup>T</sup></i> in the
+                                     * singular value decomposition
+                                     * <i>USV<sup>T</sup></i>.
+                                     */
     boost::shared_ptr<LAPACKFullMatrix<number> > svd_vt;
 };
 
@@ -562,6 +632,17 @@ LAPACKFullMatrix<number>::eigenvalue (const unsigned int i) const
 }
 
 
+template <typename number>
+number
+LAPACKFullMatrix<number>::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
 
index 6aa4b9ee0591905026ffcebbabbed27471d17b2b..34de869a1093e9a71b89f46a2350cbd7cd212d5a 100644 (file)
@@ -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
index 631796654dd8309ce05cfa5d7a903f5c10807f7c..76e78e7216a5295526ddcebbf41e213874662faf 100644 (file)
@@ -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,
index f107bdc6c9aa50591f1d7e22a63c0c1913b7a05d..24b6b8dd1d5c440e8af30cbda4f0267956d798a1 100644 (file)
@@ -24,7 +24,7 @@
 #include <lac/vector.h>
 #include <lac/full_matrix.h>
 #include <lac/compressed_simple_sparsity_pattern.h>
-
+#include <lac/vector_memory.h>
 
 // we only need output streams, but older compilers did not provide
 // them in a separate include file
index 9af29628ec7e70d9045ea213e8cb2828bc1ce22b..f6c9a2e2c754cc892a21b81c951a688edca992cb 100644 (file)
@@ -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
 //---------------------------------------------------------------------------
 
 
-#include <lac/full_matrix.templates.h>
 #include <base/logstream.h>
 #include <base/tensor.h>
 #include <base/tensor_base.h>
+#include <lac/full_matrix.templates.h>
 
 DEAL_II_NAMESPACE_OPEN
 
index f051a053740e6a09cf39ad9b53db1e8699a009d1..407a63e01495b37c3b8c2a5ffc135a91581d7a57 100644 (file)
@@ -90,7 +90,6 @@ LAPACKFullMatrix<number>::operator = (const double d)
   return *this;
 }
 
-#ifdef HAVE_LIBBLAS
 
 template <typename number>
 void
@@ -98,16 +97,55 @@ LAPACKFullMatrix<number>::vmult (
   Vector<number>       &w,
   const Vector<number> &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;i<wr.size();++i)
+         work[i] *= wr[i];
+                                        // Multiply with U
+       gemv("N", &mm, &mm, &alpha, svd_u->data(), &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;i<wr.size();++i)
+         work[i] *= wr[i];
+                                        // Multiply with V
+       gemv("T", &nn, &nn, &alpha, svd_vt->data(), &nn, &work[0], &one, &beta, w.val, &one);
+       break;
+      }
+      default:
+           Assert (false, ExcState(state));
+    }
 }
 
 
@@ -118,61 +156,148 @@ LAPACKFullMatrix<number>::Tvmult (
   const Vector<number> &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;i<wr.size();++i)
+         work[i] *= wr[i];
+                                        // Multiply with V
+       gemv("T", &nn, &nn, &alpha, svd_vt->data(), &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;i<wr.size();++i)
+         work[i] *= wr[i];
+                                        // Multiply with U
+       gemv("N", &mm, &mm, &alpha, svd_u->data(), &mm, &work[0], &one, &beta, w.val, &one);
+       break;
+      }
+      }
+      default:
+           Assert (false, ExcState(state));
+    }
 }
 
-#else
-
 
 template <typename number>
 void
-LAPACKFullMatrix<number>::vmult (
-  Vector<number>       &,
-  const Vector<number> &,
-  const bool            ) const
+LAPACKFullMatrix<number>::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<number*> (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 <typename number>
-void
-LAPACKFullMatrix<number>::Tvmult (
-  Vector<number>       &,
-  const Vector<number> &,
-  const bool            ) const
-{
-  Assert(false, ExcNeedsBLAS());
+  state = lu;
 }
 
 
-#endif
-
-#ifdef HAVE_LIBLAPACK
-
 template <typename number>
 void
-LAPACKFullMatrix<number>::compute_lu_factorization()
+LAPACKFullMatrix<number>::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<number*> (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<LAPACKFullMatrix<number> >(new LAPACKFullMatrix<number>(mm,mm));
+  svd_vt = boost::shared_ptr<LAPACKFullMatrix<number> >(new LAPACKFullMatrix<number>(nn,nn));
+  number* mu  = const_cast<number*> (svd_u->data());
+  number* mvt = const_cast<number*> (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 <typename number>
+void
+LAPACKFullMatrix<number>::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<wr.size();++i)
+    {
+      if (wr[i] > lim)
+       wr[i] = 1./wr[i];
+      else
+       wr[i] = 0.;
+    }
+  state = LAPACKSupport::inverse_svd;
 }
 
 
@@ -403,41 +528,6 @@ LAPACKFullMatrix<number>::compute_generalized_eigenvalues_symmetric (
   state = LAPACKSupport::State(eigenvalues | unusable);
 }
 
-#else
-
-template <typename number>
-void
-LAPACKFullMatrix<number>::compute_lu_factorization()
-{
-Assert(false, ExcNeedsLAPACK());
-
-}
-
-
-template <typename number>
-void
-LAPACKFullMatrix<number>::apply_lu_factorization(Vector<number>&, bool) const
-{
-  Assert(false, ExcNeedsLAPACK());
-}
-
-
-template <typename number>
-void
-LAPACKFullMatrix<number>::compute_eigenvalues(const bool /*right*/,
-                                             const bool /*left*/)
-{
-  Assert(false, ExcNeedsLAPACK());
-}
-
-
-#endif
-
-
-// template <typename number>
-// LAPACKFullMatrix<number>::()
-// {}
-
 
 // template <typename number>
 // LAPACKFullMatrix<number>::()

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.