]> https://gitweb.dealii.org/ - dealii.git/commitdiff
extend LAPACKFullMatrix to do Cholesky, inverse of SPD, norms and condition number...
authorDenis Davydov <davydden@gmail.com>
Tue, 26 Sep 2017 13:40:07 +0000 (15:40 +0200)
committerDenis Davydov <davydden@gmail.com>
Wed, 27 Sep 2017 07:14:32 +0000 (09:14 +0200)
15 files changed:
include/deal.II/lac/lapack_full_matrix.h
include/deal.II/lac/lapack_support.h
include/deal.II/lac/lapack_templates.h
source/lac/lapack_full_matrix.cc
tests/lapack/create_matrix.h [new file with mode: 0644]
tests/lapack/full_matrix_15.cc [new file with mode: 0644]
tests/lapack/full_matrix_15.output [new file with mode: 0644]
tests/lapack/full_matrix_16.cc [new file with mode: 0644]
tests/lapack/full_matrix_16.output [new file with mode: 0644]
tests/lapack/full_matrix_17.cc [new file with mode: 0644]
tests/lapack/full_matrix_17.output [new file with mode: 0644]
tests/lapack/full_matrix_18.cc [new file with mode: 0644]
tests/lapack/full_matrix_18.output [new file with mode: 0644]
tests/lapack/full_matrix_19.cc [new file with mode: 0644]
tests/lapack/full_matrix_19.output [new file with mode: 0644]

index 34d9247560727bde43bfbcc89e6b51c1497a6137..593794810bed34be22ba76c5a97b89d01e42c5cf 100644 (file)
@@ -46,7 +46,7 @@ template <typename number> class SparseMatrix;
  * usually the names chosen for the arguments in the LAPACK documentation.
  *
  * @ingroup Matrix1
- * @author Guido Kanschat, 2005
+ * @author Guido Kanschat, 2005, Denis Davydov, 2017
  */
 template <typename number>
 class LAPACKFullMatrix : public TransposeTable<number>
@@ -157,6 +157,11 @@ public:
   void reinit (const size_type rows,
                const size_type cols);
 
+  /**
+   * Assign @p property to this matrix.
+   */
+  void set_property(const LAPACKSupport::Property property);
+
   /**
    * Return the dimension of the codomain (or range) space.
    *
@@ -404,6 +409,23 @@ public:
    */
   void compute_lu_factorization ();
 
+  /**
+   * Compute the Cholesky factorization of the matrix using LAPACK function Xpotrf.
+   *
+   * @note The factorization is stored in the lower diagonal.
+   */
+  void compute_cholesky_factorization ();
+
+  /**
+   * Estimate reciprocal of the condition number (in 1-norm) of a SPD matrix
+   * using Cholesky factorization. This function can only be called if the matrix
+   * is already factorized.
+   *
+   * @param a_norm Is the 1-norm of the matrix before calling Cholesky
+   * factorization.
+   */
+  number reciprocal_condition_number(const number a_norm) const;
+
   /**
    * Compute the determinant of a matrix. As it requires the LU factorization of
    * the matrix, this function can only be called after
@@ -412,8 +434,23 @@ public:
   number determinant () const;
 
   /**
-   * Invert the matrix by first computing an LU factorization with the LAPACK
-   * function Xgetrf and then building the actual inverse using Xgetri.
+   * Compute $L_1$ norm.
+   */
+  number l1_norm() const;
+
+  /**
+   * Compute $L_\infty$ norm.
+   */
+  number linfty_norm() const;
+
+  /**
+   * Compute Frobenius norm
+   */
+  number frobenius_norm() const;
+
+  /**
+   * Invert the matrix by first computing an LU/Cholesky factorization with the LAPACK
+   * function Xgetrf/Xpotrf and then building the actual inverse using Xgetri/Xpotri.
    */
   void invert ();
 
@@ -613,6 +650,10 @@ public:
                         const double        threshold   = 0.) const;
 
 private:
+  /**
+   * Internal function to compute various norms.
+   */
+  number norm(const char type) const;
 
   /**
    * Since LAPACK operations notoriously change the meaning of the matrix
@@ -621,10 +662,10 @@ private:
   LAPACKSupport::State state;
 
   /**
-   * Additional properties of the matrix which may help to select more
+   * Additional property of the matrix which may help to select more
    * efficient LAPACK functions.
    */
-  LAPACKSupport::Properties properties;
+  LAPACKSupport::Property property;
 
   /**
    * The working array used for some LAPACK functions.
index 53bbb000a209cb119965675bea5bb660eed1f627..6b4032999a48c55e74539eb4af7faec004077d91 100644 (file)
@@ -85,7 +85,7 @@ namespace LAPACKSupport
    * A matrix can have certain features allowing for optimization, but hard to
    * test. These are listed here.
    */
-  enum Properties
+  enum Property
   {
     /// No special properties
     general = 0,
@@ -101,6 +101,28 @@ namespace LAPACKSupport
     hessenberg = 8
   };
 
+  /**
+   * Function printing the name of a Property.
+   */
+  inline const char *property_name(const Property s)
+  {
+    switch (s)
+      {
+      case general:
+        return "general";
+      case symmetric:
+        return "symmetric";
+      case upper_triangular:
+        return "upper triangular";
+      case lower_triangular:
+        return "lower triangular";
+      case diagonal:
+        return "diagonal";
+      case hessenberg:
+        return "Hessenberg";
+      }
+  }
+
   /**
    * Character constant.
    */
@@ -149,6 +171,15 @@ namespace LAPACKSupport
                  << "The function cannot be called while the matrix is in state "
                  << state_name(arg1));
 
+  /**
+   * Exception thrown when a matrix does not have suitable properties for an
+   * operation.
+   */
+  DeclException1(ExcProperty, Property,
+                 << "The function cannot be called while the "
+                 << property_name(arg1)
+                 << " matrix.");
+
   /**
    * This exception is thrown if a certain LAPACK function is not available
    * because no LAPACK installation was detected during configuration.
index 8f9398638165191b43a9048dd000d19d02da32a0..c110ff5a8896f1064a150ad959eff2e150afc820 100644 (file)
@@ -65,6 +65,33 @@ extern "C"
                 int *ipiv, double *inv_work, const int *lwork, int *info);
   void sgetri_ (const int *n, float *A, const int *lda,
                 int *ipiv, float *inv_work, const int *lwork, int *info);
+// Compute Cholesky factorization of SPD
+  void dpotrf_ (const char *uplo, const int *n, double *A,
+                const int *lda, int *info);
+  void spotrf_ (const char *uplo, const int *n, float *A,
+                const int *lda, int *info);
+// Estimate the reciprocal of the condition number in 1-norm from Cholesky
+  void dpocon_ (const char *uplo, const int *n, const double *A,
+                const int *lda, const double *anorm, double *rcond,
+                double *work, int *iwork, int *info);
+  void spocon_ (const char *uplo, const int *n, const float *A,
+                const int *lda, const float *anorm, float *rcond,
+                float *work, int *iwork, int *info);
+// Computes the inverse from Cholesky
+  void dpotri_ (const char *uplo, const int *n, double *A,
+                const int *lda, int *info);
+  void spotri_ (const char *uplo, const int *n, float *A,
+                const int *lda, int *info);
+// Norms
+  double dlansy_ (const char *norm, const char *uplo, const int *n, const double *A,
+                  const int *lda, double *work);
+  float  slansy_ (const char *norm, const char *uplo, const int *n, const float *A,
+                  const int *lda, float *work);
+  double dlange_ (const char *norm, const int *m, const int *n, const double *A,
+                  const int *lda, double *work);
+  float slange_ (const char *norm, const int *m, const int *n, const float *A,
+                 const int *lda, float *work);
+
 // Compute QR factorization (Householder)
   void dgeqrf_ (const int *m, const int *n, double *A,
                 const int *lda, double *tau, double *work,
@@ -365,6 +392,156 @@ gemm (const char *, const char *, const int *, const int *, const int *, const f
 #endif
 
 
+
+/// Template wrapper for potrf
+template <typename number1>
+inline void
+potrf (const char *uplo, const int *n, number1 *A, const int *lda, int *info)
+{
+  Assert (false, ExcNotImplemented());
+}
+
+inline void
+potrf (const char *uplo, const int *n, double *A, const int *lda, int *info)
+{
+#ifdef DEAL_II_WITH_LAPACK
+  dpotrf_ (uplo,n,A,lda,info);
+#else
+  Assert (false, LAPACKSupport::ExcMissing("dpotrf"));
+#endif
+}
+
+inline void
+potrf (const char *uplo, const int *n, float *A, const int *lda, int *info)
+{
+#ifdef DEAL_II_WITH_LAPACK
+  spotrf_ (uplo,n,A,lda,info);
+#else
+  Assert (false, LAPACKSupport::ExcMissing("spotrf"));
+#endif
+}
+
+
+
+/// Template wrapper for pocon
+template <typename number1>
+inline void
+pocon (const char *uplo, const int *n, const number1 *A, const int *lda, const number1 *anorm, number1 *rcond, number1 *work, int *iwork, int *info)
+{
+  Assert (false, ExcNotImplemented());
+}
+
+inline void
+pocon (const char *uplo, const int *n, const double *A, const int *lda, const double *anorm, double *rcond, double *work, int *iwork, int *info)
+{
+#ifdef DEAL_II_WITH_LAPACK
+  dpocon_ (uplo,n,A,lda,anorm,rcond,work,iwork,info);
+#else
+  Assert (false, LAPACKSupport::ExcMissing("dpocon"));
+#endif
+}
+
+inline void
+pocon (const char *uplo, const int *n, const float *A, const int *lda, const float *anorm, float *rcond, float *work, int *iwork, int *info)
+{
+#ifdef DEAL_II_WITH_LAPACK
+  spocon_ (uplo,n,A,lda,anorm,rcond,work,iwork,info);
+#else
+  Assert (false, LAPACKSupport::ExcMissing("dpocon"));
+#endif
+}
+
+/// Template wrapper for potri
+template <typename number1>
+inline void
+potri(const char *uplo, const int *n, number1 *A, const int *lda, int *info)
+{
+  Assert (false, ExcNotImplemented());
+}
+
+inline void
+potri(const char *uplo, const int *n, double *A, const int *lda, int *info)
+{
+#ifdef DEAL_II_WITH_LAPACK
+  dpotri_(uplo,n,A,lda,info);
+#else
+  Assert (false, LAPACKSupport::ExcMissing("dpotri"));
+#endif
+}
+
+inline void
+potri(const char *uplo, const int *n, float *A, const int *lda, int *info)
+{
+#ifdef DEAL_II_WITH_LAPACK
+  spotri_(uplo,n,A,lda,info);
+#else
+  Assert (false, LAPACKSupport::ExcMissing("spotri"));
+#endif
+}
+
+
+/// Template wrapper for lansy
+template <typename number>
+inline
+number lansy (const char *norm, const char *uplo, const int *n, const number *A, const int *lda, number *work)
+{
+  Assert (false, ExcNotImplemented());
+}
+
+inline
+double lansy (const char *norm, const char *uplo, const int *n, const double *A, const int *lda, double *work)
+{
+#ifdef DEAL_II_WITH_LAPACK
+  return dlansy_(norm,uplo,n,A,lda,work);
+#else
+  Assert (false, LAPACKSupport::ExcMissing("lansy"));
+  return 0.;
+#endif
+}
+
+inline
+float lansy (const char *norm, const char *uplo, const int *n, const float *A, const int *lda, float *work)
+{
+#ifdef DEAL_II_WITH_LAPACK
+  return slansy_(norm,uplo,n,A,lda,work);
+#else
+  Assert (false, LAPACKSupport::ExcMissing("lansy"));
+  return 0.;
+#endif
+}
+
+
+
+/// Template wrapper for lange
+template <typename number>
+inline
+number lange (const char *norm, const int *m, const int *n, const number *A, const int *lda, number *work)
+{
+  Assert (false, ExcNotImplemented());
+}
+
+inline
+double lange (const char *norm, const int *m, const int *n, const double *A, const int *lda, double *work)
+{
+#ifdef DEAL_II_WITH_LAPACK
+  return dlange_(norm,m,n,A,lda,work);
+#else
+  Assert (false, LAPACKSupport::ExcMissing("lange"));
+  return 0.;
+#endif
+}
+
+inline
+float lange (const char *norm, const int *m, const int *n, const float *A, const int *lda, float *work)
+{
+#ifdef DEAL_II_WITH_LAPACK
+  return slange_(norm,m,n,A,lda,work);
+#else
+  Assert (false, LAPACKSupport::ExcMissing("lange"));
+  return 0.;
+#endif
+}
+
 /// Template wrapper for LAPACK functions dgetrf and sgetrf
 template <typename number1>
 inline void
index ea9a3b84b6754e9892d812ad0da0ac8a6b21270b..763daac138aa40bab95e6293fdab2cd7eb29e51b 100644 (file)
@@ -33,7 +33,8 @@ template <typename number>
 LAPACKFullMatrix<number>::LAPACKFullMatrix (const size_type n)
   :
   TransposeTable<number> (n,n),
-  state (matrix)
+  state (matrix),
+  property(general)
 {}
 
 
@@ -42,7 +43,8 @@ LAPACKFullMatrix<number>::LAPACKFullMatrix (const size_type m,
                                             const size_type n)
   :
   TransposeTable<number> (m, n),
-  state (matrix)
+  state (matrix),
+  property(general)
 {}
 
 
@@ -50,7 +52,8 @@ template <typename number>
 LAPACKFullMatrix<number>::LAPACKFullMatrix (const LAPACKFullMatrix &M)
   :
   TransposeTable<number> (M),
-  state (matrix)
+  state (matrix),
+  property(general)
 {}
 
 
@@ -59,7 +62,8 @@ LAPACKFullMatrix<number> &
 LAPACKFullMatrix<number>::operator = (const LAPACKFullMatrix<number> &M)
 {
   TransposeTable<number>::operator=(M);
-  state = LAPACKSupport::matrix;
+  state = M.state;
+  property = M.property;
   return *this;
 }
 
@@ -95,6 +99,7 @@ LAPACKFullMatrix<number>::operator = (const FullMatrix<number2> &M)
       (*this)(i,j) = M(i,j);
 
   state = LAPACKSupport::matrix;
+  property = LAPACKSupport::general;
   return *this;
 }
 
@@ -111,6 +116,7 @@ LAPACKFullMatrix<number>::operator = (const SparseMatrix<number2> &M)
       (*this)(i,j) = M.el(i,j);
 
   state = LAPACKSupport::matrix;
+  property = LAPACKSupport::general;
   return *this;
 }
 
@@ -514,6 +520,127 @@ LAPACKFullMatrix<number>::compute_lu_factorization()
 }
 
 
+
+template <typename number>
+void
+LAPACKFullMatrix<number>::set_property(const Property p)
+{
+  property = p;
+}
+
+
+
+template <typename number>
+number LAPACKFullMatrix<number>::l1_norm() const
+{
+  const char type('O');
+  return norm(type);
+}
+
+
+
+template <typename number>
+number LAPACKFullMatrix<number>::linfty_norm() const
+{
+  const char type('I');
+  return norm(type);
+}
+
+
+
+template <typename number>
+number LAPACKFullMatrix<number>::frobenius_norm() const
+{
+  const char type('F');
+  return norm(type);
+}
+
+
+
+template <typename number>
+number LAPACKFullMatrix<number>::norm(const char type) const
+{
+  Assert (state == LAPACKSupport::matrix ||
+          state == LAPACKSupport::inverse_matrix,
+          ExcMessage("norms can be called in matrix state only."));
+
+  const int N = this->n_cols();
+  const int M = this->n_rows();
+  const number *values = &this->values[0];
+  std::vector<number> w;
+  if (property == symmetric)
+    {
+      const int lda = std::max(1,N);
+      const int lwork = (type == 'I' || type == 'O') ?
+                        std::max(1,N) :
+                        0;
+      w.resize(lwork);
+      return lansy (&type, &LAPACKSupport::L, &N, values, &lda, &w[0]);
+    }
+  else
+    {
+      const int lda = std::max(1,M);
+      const int lwork = (type == 'I') ?
+                        std::max(1,M) :
+                        0;
+      w.resize(lwork);
+      return lange (&type, &M, &N, values, &lda, &w[0]);
+    }
+}
+
+
+template <typename number>
+void
+LAPACKFullMatrix<number>::compute_cholesky_factorization()
+{
+  Assert(state == matrix, ExcState(state));
+  Assert(property == symmetric, ExcProperty(property));
+  state = LAPACKSupport::unusable;
+
+  const int mm = this->n_rows();
+  const int nn = this->n_cols();
+  Assert (mm == nn, ExcDimensionMismatch(mm,nn));
+
+  number *values = &this->values[0];
+  int info = 0;
+  const int lda = std::max(1,nn);
+  potrf (&LAPACKSupport::L, &nn, values, &lda, &info);
+
+  // info < 0 : the info-th argument had an illegal value
+  Assert(info >= 0, ExcInternalError());
+
+  state = cholesky;
+  AssertThrow(info == 0, LACExceptions::ExcSingular());
+}
+
+
+
+template <typename number>
+number
+LAPACKFullMatrix<number>::reciprocal_condition_number(const number a_norm) const
+{
+  Assert(state == cholesky, ExcState(state));
+  number rcond = 0.;
+
+  const int N = this->n_rows();
+  const number *values = &this->values[0];
+  int info = 0;
+  const int lda = std::max(1,N);
+  std::vector<number> w(3*N);
+  std::vector<int> iw(N);
+
+  // use the same uplo as in Cholesky
+  pocon (&LAPACKSupport::L, &N, values, &lda,
+         &a_norm, &rcond,
+         &w[0], &iw[0], &info);
+
+  Assert(info >= 0, ExcInternalError());
+
+  return rcond;
+}
+
+
+
 template <typename number>
 void
 LAPACKFullMatrix<number>::compute_svd()
@@ -584,26 +711,36 @@ template <typename number>
 void
 LAPACKFullMatrix<number>::invert()
 {
-  Assert(state == matrix || state == lu,
+  Assert(state == matrix || state == lu || state == cholesky,
          ExcState(state));
   const int mm = this->n_rows();
   const int nn = this->n_cols();
   Assert (nn == mm, ExcNotQuadratic());
 
   number *values = const_cast<number *> (&this->values[0]);
-  ipiv.resize(mm);
   int info = 0;
 
-  if (state == matrix)
+  if (property != symmetric)
     {
-      getrf(&mm, &nn, values, &mm, &ipiv[0], &info);
+      if (state == matrix)
+        compute_lu_factorization();
 
-      Assert(info >= 0, ExcInternalError());
-      AssertThrow(info == 0, LACExceptions::ExcSingular());
+      ipiv.resize(mm);
+      inv_work.resize (mm);
+      getri(&mm, values, &mm, &ipiv[0], &inv_work[0], &mm, &info);
+    }
+  else
+    {
+      if (state == matrix)
+        compute_cholesky_factorization();
+
+      const int lda = std::max(1,nn);
+      potri(&LAPACKSupport::L, &nn, values,&lda,&info);
+      // inverse is stored in lower diagonal, set the upper diagonal appropriately:
+      for (int i=0; i < nn; ++i)
+        for (int j=i+1; j < nn; ++j)
+          this->el(i,j) = this->el(j,i);
     }
-
-  inv_work.resize (mm);
-  getri(&mm, values, &mm, &ipiv[0], &inv_work[0], &mm, &info);
 
   Assert(info >= 0, ExcInternalError());
   AssertThrow(info == 0, LACExceptions::ExcSingular());
@@ -1016,7 +1153,8 @@ LAPACKFullMatrix<number>::print_formatted (
   Assert ((!this->empty()) || (this->n_cols()+this->n_rows()==0),
           ExcInternalError());
   Assert (state == LAPACKSupport::matrix ||
-          state == LAPACKSupport::inverse_matrix,
+          state == LAPACKSupport::inverse_matrix ||
+          state == LAPACKSupport::cholesky,
           ExcState(state));
 
   // set output format, but store old
diff --git a/tests/lapack/create_matrix.h b/tests/lapack/create_matrix.h
new file mode 100644 (file)
index 0000000..2e14784
--- /dev/null
@@ -0,0 +1,52 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2017 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+// auxiliary function to create a SPD FullMatrix
+
+#include "../tests.h"
+
+template <typename FullMatrix>
+void create_spd (FullMatrix &A)
+{
+  const unsigned int size = A.n();
+  Assert (size == A.m(), ExcDimensionMismatch(size,A.m()));
+
+  for (unsigned int i = 0; i < size; ++i)
+    for (unsigned int j = i; j < size; ++j)
+      {
+        const double val = (double)rand()/RAND_MAX;
+        Assert (val >= 0. && val <= 1.,
+                ExcInternalError());
+        if (i==j)
+          // since A(i,j) < 1 and
+          // a symmetric diagonally dominant matrix is SPD
+          A(i,j) = val + size;
+        else
+          {
+            A(i,j) = val;
+            A(j,i) = val;
+          }
+      }
+}
+
+
+
+template <typename FullMatrix>
+void create_random (FullMatrix &A)
+{
+  for (unsigned int i = 0; i < A.m(); ++i)
+    for (unsigned int j = 0; j < A.n(); ++j)
+      A(i,j) = (double)rand()/RAND_MAX;
+}
diff --git a/tests/lapack/full_matrix_15.cc b/tests/lapack/full_matrix_15.cc
new file mode 100644 (file)
index 0000000..faec1b3
--- /dev/null
@@ -0,0 +1,78 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2005 - 2015 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+// test LAPACKFullMatrix::compute_cholesky_factorization by comparing with FullMatrix
+
+#include "../tests.h"
+#include "create_matrix.h"
+#include <deal.II/lac/lapack_full_matrix.h>
+#include <deal.II/lac/full_matrix.h>
+#include <deal.II/lac/vector.h>
+
+#include <iostream>
+
+
+template <typename NumberType>
+void
+test(const unsigned int size)
+{
+  // Full matrix:
+  FullMatrix<NumberType> F(size), C(size);
+  create_spd(F);
+  C.cholesky(F);
+
+  // Lapack:
+  LAPACKFullMatrix<NumberType> M(size);
+  M = F;
+  M.set_property(LAPACKSupport::symmetric);
+  M.compute_cholesky_factorization();
+  // factorization is stored in the lower diagonal part
+  for (unsigned int i = 0; i < size; ++i)
+    for (unsigned int j=i+1; j < size; ++j)
+      M(i,j) = 0.;
+
+  FullMatrix<NumberType> diff(size);
+  diff = M;
+  diff.add(-1., C);
+
+  const NumberType error = diff.frobenius_norm();
+  deallog << size << " : " << diff.frobenius_norm() << std::endl;
+  if (false)
+    {
+      std::cout << "Lapack:" << std::endl;
+      M.print_formatted(std::cout);
+      std::cout << "FullMatrix:" << std::endl;
+      C.print_formatted(std::cout);
+      std::cout << std::flush;
+      AssertThrow(false, ExcInternalError());
+    }
+}
+
+
+int main()
+{
+  const std::string logname = "output";
+  std::ofstream logfile(logname.c_str());
+  logfile.precision(3);
+  deallog.attach(logfile);
+
+  const std::vector<unsigned int> sizes = {{1,3,11,17,32,64,200,391}};
+  for (const auto &s : sizes)
+    {
+      // test<float>(s);
+      test<double>(s);
+    }
+}
diff --git a/tests/lapack/full_matrix_15.output b/tests/lapack/full_matrix_15.output
new file mode 100644 (file)
index 0000000..5e55595
--- /dev/null
@@ -0,0 +1,9 @@
+
+DEAL::1 : 0.00000
+DEAL::3 : 0.00000
+DEAL::11 : 0.00000
+DEAL::17 : 0.00000
+DEAL::32 : 0.00000
+DEAL::64 : 0.00000
+DEAL::200 : 0.00000
+DEAL::391 : 0.00000
diff --git a/tests/lapack/full_matrix_16.cc b/tests/lapack/full_matrix_16.cc
new file mode 100644 (file)
index 0000000..a8fcaae
--- /dev/null
@@ -0,0 +1,61 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2005 - 2015 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+// test LAPACKFullMatrix:: norms
+
+#include "../tests.h"
+#include "create_matrix.h"
+#include <deal.II/lac/lapack_full_matrix.h>
+#include <deal.II/lac/full_matrix.h>
+#include <deal.II/lac/vector.h>
+
+#include <iostream>
+
+
+template <typename NumberType>
+void
+test(const unsigned int size)
+{
+  // Full matrix:
+  FullMatrix<NumberType> F(size);
+  create_spd(F);
+
+  // Lapack:
+  LAPACKFullMatrix<NumberType> M(size);
+  M = F;
+  M.set_property(LAPACKSupport::symmetric);
+
+  deallog << "l1:        " << (F.l1_norm()        - M.l1_norm())        << std::endl
+          << "linfty:    " << (F.linfty_norm()    - M.linfty_norm())    << std::endl
+          << "frobenius: " << (F.frobenius_norm() - M.frobenius_norm()) << std::endl;
+}
+
+
+int main()
+{
+  const std::string logname = "output";
+  std::ofstream logfile(logname.c_str());
+  logfile.precision(3);
+  deallog.attach(logfile);
+
+  const std::vector<unsigned int> sizes = {{1,3,11,17,32,64,200,391}};
+  for (const auto &s : sizes)
+    {
+      deallog << "size=" << s << std::endl;
+      // test<float>(s);
+      test<double>(s);
+    }
+}
diff --git a/tests/lapack/full_matrix_16.output b/tests/lapack/full_matrix_16.output
new file mode 100644 (file)
index 0000000..4fc3f3b
--- /dev/null
@@ -0,0 +1,33 @@
+
+DEAL::size=1
+DEAL::l1:        0.00000
+DEAL::linfty:    0.00000
+DEAL::frobenius: 0.00000
+DEAL::size=3
+DEAL::l1:        0.00000
+DEAL::linfty:    0.00000
+DEAL::frobenius: 0.00000
+DEAL::size=11
+DEAL::l1:        0.00000
+DEAL::linfty:    0.00000
+DEAL::frobenius: 0.00000
+DEAL::size=17
+DEAL::l1:        0.00000
+DEAL::linfty:    0.00000
+DEAL::frobenius: 0.00000
+DEAL::size=32
+DEAL::l1:        0.00000
+DEAL::linfty:    0.00000
+DEAL::frobenius: 0.00000
+DEAL::size=64
+DEAL::l1:        0.00000
+DEAL::linfty:    0.00000
+DEAL::frobenius: 0.00000
+DEAL::size=200
+DEAL::l1:        0.00000
+DEAL::linfty:    0.00000
+DEAL::frobenius: 0.00000
+DEAL::size=391
+DEAL::l1:        0.00000
+DEAL::linfty:    0.00000
+DEAL::frobenius: 0.00000
diff --git a/tests/lapack/full_matrix_17.cc b/tests/lapack/full_matrix_17.cc
new file mode 100644 (file)
index 0000000..f1c9065
--- /dev/null
@@ -0,0 +1,66 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2005 - 2015 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+// test LAPACKFullMatrix::reciprocal_condition_number() by comparing estimated
+// value to 1 / ( ||A||_1 * ||A^{-1}||_1 )
+
+#include "../tests.h"
+#include "create_matrix.h"
+#include <deal.II/lac/lapack_full_matrix.h>
+#include <deal.II/lac/full_matrix.h>
+#include <deal.II/lac/vector.h>
+
+#include <iostream>
+
+
+template <typename NumberType>
+void
+test(const unsigned int size)
+{
+  // Full matrix:
+  FullMatrix<NumberType> F(size), invF(size);
+  create_spd(F);
+  invF.invert(F);
+  const double l1 = F.l1_norm();
+  const double inv_l1 = invF.l1_norm();
+
+  // Lapack:
+  LAPACKFullMatrix<NumberType> M(size);
+  M = F;
+  M.set_property(LAPACKSupport::symmetric);
+  const double la_l1 = M.l1_norm();
+  M.compute_cholesky_factorization();
+  const double rcond = M.reciprocal_condition_number(la_l1);
+
+  deallog << 1./(l1*inv_l1) << " " << rcond << std::endl;
+}
+
+
+int main()
+{
+  const std::string logname = "output";
+  std::ofstream logfile(logname.c_str());
+  logfile.precision(3);
+  deallog.attach(logfile);
+
+  const std::vector<unsigned int> sizes = {{1,3,11,17,32,64,200,391}};
+  for (const auto &s : sizes)
+    {
+      deallog << "size=" << s << std::endl;
+      // test<float>(s);
+      test<double>(s);
+    }
+}
diff --git a/tests/lapack/full_matrix_17.output b/tests/lapack/full_matrix_17.output
new file mode 100644 (file)
index 0000000..a0247b2
--- /dev/null
@@ -0,0 +1,17 @@
+
+DEAL::size=1
+DEAL::1.00000 1.00000
+DEAL::size=3
+DEAL::0.482025 0.563918
+DEAL::size=11
+DEAL::0.425003 0.466425
+DEAL::size=17
+DEAL::0.426863 0.467023
+DEAL::size=32
+DEAL::0.441550 0.493985
+DEAL::size=64
+DEAL::0.452016 0.459357
+DEAL::size=200
+DEAL::0.467209 0.490028
+DEAL::size=391
+DEAL::0.468957 0.486322
diff --git a/tests/lapack/full_matrix_18.cc b/tests/lapack/full_matrix_18.cc
new file mode 100644 (file)
index 0000000..45a68bc
--- /dev/null
@@ -0,0 +1,76 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2005 - 2015 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+// test LAPACKFullMatrix::invert() wth SPD Cholesky by comparing to FullMatrix
+
+#include "../tests.h"
+#include "create_matrix.h"
+#include <deal.II/lac/lapack_full_matrix.h>
+#include <deal.II/lac/full_matrix.h>
+#include <deal.II/lac/vector.h>
+
+#include <iostream>
+
+
+template <typename NumberType>
+void
+test(const unsigned int size)
+{
+  // Full matrix:
+  FullMatrix<NumberType> F(size), invF(size);
+  create_spd(F);
+  invF.invert(F);
+
+  // Lapack:
+  LAPACKFullMatrix<NumberType> M(size);
+  M = F;
+  M.set_property(LAPACKSupport::symmetric);
+  M.compute_cholesky_factorization();
+  M.invert();
+
+  FullMatrix<NumberType> diff(size);
+  diff = M;
+  diff.add(-1., invF);
+  const double error = diff.frobenius_norm();
+  deallog << error << std::endl;
+
+  if (error > 1e-10)
+    {
+      std::cout << "Lapack:" << std::endl;
+      M.print_formatted(std::cout);
+      std::cout << "FullMatrix:" << std::endl;
+      invF.print_formatted(std::cout);
+      std::cout << std::flush;
+      AssertThrow(false, ExcInternalError());
+    }
+}
+
+
+int main()
+{
+  const std::string logname = "output";
+  std::ofstream logfile(logname.c_str());
+  logfile.precision(3);
+  deallog.attach(logfile);
+
+  const std::vector<unsigned int> sizes = {{1,3,11,17,32,64,200,391}};
+  for (const auto &s : sizes)
+    {
+      deallog << "size=" << s << std::endl;
+      // test<float>(s);
+      test<double>(s);
+    }
+}
diff --git a/tests/lapack/full_matrix_18.output b/tests/lapack/full_matrix_18.output
new file mode 100644 (file)
index 0000000..7204301
--- /dev/null
@@ -0,0 +1,17 @@
+
+DEAL::size=1
+DEAL::0.00000
+DEAL::size=3
+DEAL::0.00000
+DEAL::size=11
+DEAL::0.00000
+DEAL::size=17
+DEAL::0.00000
+DEAL::size=32
+DEAL::0.00000
+DEAL::size=64
+DEAL::0.00000
+DEAL::size=200
+DEAL::0.00000
+DEAL::size=391
+DEAL::0.00000
diff --git a/tests/lapack/full_matrix_19.cc b/tests/lapack/full_matrix_19.cc
new file mode 100644 (file)
index 0000000..c9df4fc
--- /dev/null
@@ -0,0 +1,60 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2005 - 2015 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+// test LAPACKFullMatrix:: norms for non-symmetric matrices
+
+#include "../tests.h"
+#include "create_matrix.h"
+#include <deal.II/lac/lapack_full_matrix.h>
+#include <deal.II/lac/full_matrix.h>
+#include <deal.II/lac/vector.h>
+
+#include <iostream>
+
+
+template <typename NumberType>
+void
+test(const unsigned int size)
+{
+  // Full matrix:
+  FullMatrix<NumberType> F(size);
+  create_random(F);
+
+  // Lapack:
+  LAPACKFullMatrix<NumberType> M(size);
+  M = F;
+
+  deallog << "l1:        " << (F.l1_norm()        - M.l1_norm())        << std::endl
+          << "linfty:    " << (F.linfty_norm()    - M.linfty_norm())    << std::endl
+          << "frobenius: " << (F.frobenius_norm() - M.frobenius_norm()) << std::endl;
+}
+
+
+int main()
+{
+  const std::string logname = "output";
+  std::ofstream logfile(logname.c_str());
+  logfile.precision(3);
+  deallog.attach(logfile);
+
+  const std::vector<unsigned int> sizes = {{1,3,11,17,32,64,200,391}};
+  for (const auto &s : sizes)
+    {
+      deallog << "size=" << s << std::endl;
+      // test<float>(s);
+      test<double>(s);
+    }
+}
diff --git a/tests/lapack/full_matrix_19.output b/tests/lapack/full_matrix_19.output
new file mode 100644 (file)
index 0000000..4fc3f3b
--- /dev/null
@@ -0,0 +1,33 @@
+
+DEAL::size=1
+DEAL::l1:        0.00000
+DEAL::linfty:    0.00000
+DEAL::frobenius: 0.00000
+DEAL::size=3
+DEAL::l1:        0.00000
+DEAL::linfty:    0.00000
+DEAL::frobenius: 0.00000
+DEAL::size=11
+DEAL::l1:        0.00000
+DEAL::linfty:    0.00000
+DEAL::frobenius: 0.00000
+DEAL::size=17
+DEAL::l1:        0.00000
+DEAL::linfty:    0.00000
+DEAL::frobenius: 0.00000
+DEAL::size=32
+DEAL::l1:        0.00000
+DEAL::linfty:    0.00000
+DEAL::frobenius: 0.00000
+DEAL::size=64
+DEAL::l1:        0.00000
+DEAL::linfty:    0.00000
+DEAL::frobenius: 0.00000
+DEAL::size=200
+DEAL::l1:        0.00000
+DEAL::linfty:    0.00000
+DEAL::frobenius: 0.00000
+DEAL::size=391
+DEAL::l1:        0.00000
+DEAL::linfty:    0.00000
+DEAL::frobenius: 0.00000

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.