]> https://gitweb.dealii.org/ - dealii.git/commitdiff
lapack: add solve() for triangular matrices or factorized (LU/Cholesky) states
authorDenis Davydov <davydden@gmail.com>
Wed, 6 Dec 2017 15:00:36 +0000 (16:00 +0100)
committerDenis Davydov <davydden@gmail.com>
Wed, 6 Dec 2017 15:00:36 +0000 (16:00 +0100)
include/deal.II/lac/lapack_full_matrix.h
include/deal.II/lac/lapack_templates.h
source/lac/lapack_full_matrix.cc
tests/lapack/full_matrix_24.cc [new file with mode: 0644]
tests/lapack/full_matrix_24.output [new file with mode: 0644]

index 83660f69b990108de365684ef7fc9020099cd07b..daf294a1f1406ed8d264b2c0aa33f31ebd85e591 100644 (file)
@@ -484,6 +484,23 @@ public:
    */
   void invert ();
 
+  /**
+   * Solve the linear system with right hand side. The matrix should
+   * be either triangular or LU factorization should be previously computed.
+   *
+   * The flag transposed indicates whether the solution of the transposed
+   * system is to be performed.
+   */
+  void solve(Vector<number> &v,
+             const bool transposed) const;
+
+  /**
+   * Same as above but for multiple right hand sides (as many as there
+   * are columns in the matrix @p B).
+   */
+  void solve(LAPACKFullMatrix<number> &B,
+             const bool transposed) const;
+
   /**
    * Solve the linear system with right hand side given by applying
    * forward/backward substitution to the previously computed LU
@@ -491,7 +508,10 @@ public:
    *
    * The flag transposed indicates whether the solution of the transposed
    * system is to be performed.
+   *
+   * @deprecated use solve() instead.
    */
+  DEAL_II_DEPRECATED
   void apply_lu_factorization (Vector<number> &v,
                                const bool      transposed) const;
 
@@ -503,7 +523,10 @@ public:
    *
    * The flag transposed indicates whether the solution of the transposed
    * system is to be performed.
+   *
+   * @deprecated use solve() instead.
    */
+  DEAL_II_DEPRECATED
   void apply_lu_factorization (LAPACKFullMatrix<number> &B,
                                const bool                transposed) const;
 
index a48ee3739af10429e6507047e463ee8f23b1d3ce..8696a5d382996b8471051d3f56afe9ed2745718d 100644 (file)
@@ -76,6 +76,15 @@ extern "C"
                 const int *lda, int *info);
   void spotrf_ (const char *uplo, const int *n, float *A,
                 const int *lda, int *info);
+// Apply forward/backward substitution to Cholesky factorization
+  void dpotrs_ (const char *uplo, const int *n, const int *nrhs,
+                const double *A, const int *lda,
+                double *B, const int *ldb,
+                int *info);
+  void spotrs_ (const char *uplo, const int *n, const int *nrhs,
+                const float *A, const int *lda,
+                float *B, const int *ldb,
+                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,
@@ -819,6 +828,57 @@ getrs (const char *, const int *, const int *, const float *, const int *, const
 #endif
 
 
+
+///  Template wrapper for LAPACK functions dpotrs and spotrs
+template <typename number>
+inline void
+potrs (const char *, const int *, const int *,
+       const number *, const int *,
+       number *, const int *,
+       int *)
+{
+  Assert (false, ExcNotImplemented());
+}
+
+#ifdef DEAL_II_WITH_LAPACK
+inline void
+potrs (const char *uplo, const int *n, const int *nrhs,
+       const double *A, const int *lda,
+       double *B, const int *ldb,
+       int *info)
+{
+  dpotrs_(uplo,n,nrhs,A,lda,B,ldb,info);
+}
+inline void
+potrs (const char *uplo, const int *n, const int *nrhs,
+       const float *A, const int *lda,
+       float *B, const int *ldb,
+       int *info)
+{
+  spotrs_(uplo,n,nrhs,A,lda,B,ldb,info);
+}
+#else
+inline void
+potrs (const char *, const int *, const int *,
+       const double *, const int *,
+       double *, const int *,
+       int *)
+{
+  Assert (false, LAPACKSupport::ExcMissing("dpotrs"));
+}
+inline void
+potrs (const char *, const int *, const int *,
+       const float *, const int *,
+       flaot *, const int *,
+       int *)
+{
+  Assert (false, LAPACKSupport::ExcMissing("spotrs"));
+}
+#endif
+
+
+
+
 /// Template wrapper for LAPACK functions dgetri and sgetri
 template <typename number1, typename number2>
 inline void
index 3cfc34d7417164a5ca3e083844f847d26de837d6..a12e833979f44a50b3576f4f04f30fef4ad42489 100644 (file)
@@ -888,50 +888,123 @@ LAPACKFullMatrix<number>::invert()
 }
 
 
+
 template <typename number>
 void
-LAPACKFullMatrix<number>::apply_lu_factorization(Vector<number> &v,
-                                                 const bool transposed) const
+LAPACKFullMatrix<number>::solve(Vector<number> &v,
+                                const bool transposed) const
 {
-  Assert(state == lu, ExcState(state));
   Assert(this->n_rows() == this->n_cols(),
          LACExceptions::ExcNotQuadratic());
   AssertDimension(this->n_rows(), v.size());
-
   const char *trans = transposed ? &T : &N;
   const int nn = this->n_cols();
   const number *values = &this->values[0];
+  const int n_rhs = 1;
   int info = 0;
 
-  getrs(trans, &nn, &one, values, &nn, ipiv.data(),
-        v.begin(), &nn, &info);
+  if (state == lu)
+    {
+      getrs(trans, &nn, &n_rhs, values, &nn, ipiv.data(),
+            v.begin(), &nn, &info);
+    }
+  else if (state == cholesky)
+    {
+      potrs(&LAPACKSupport::L, &nn, &n_rhs,
+            values, &nn, v.begin(), &nn, &info);
+    }
+  else if (property == upper_triangular ||
+           property == lower_triangular)
+    {
+      const char uplo = (property == upper_triangular ? LAPACKSupport::U : LAPACKSupport::L);
+
+      const int lda = nn;
+      const int ldb = nn;
+      trtrs (&uplo, trans, "N",
+             &nn, &n_rhs,
+             values, &lda,
+             v.begin(), &ldb, &info);
+    }
+  else
+    {
+      Assert(false,
+             ExcMessage("The matrix has to be either factorized or triangular."));
+    }
 
   Assert(info == 0, ExcInternalError());
 }
 
 
+
 template <typename number>
 void
-LAPACKFullMatrix<number>::apply_lu_factorization(LAPACKFullMatrix<number> &B,
-                                                 const bool transposed) const
+LAPACKFullMatrix<number>::solve(LAPACKFullMatrix<number> &B,
+                                const bool transposed) const
 {
-  Assert(state == lu, ExcState(state));
-  Assert(B.state == matrix, ExcState(state));
-  Assert(this->n_rows() == this->n_cols(), LACExceptions::ExcNotQuadratic());
-  AssertDimension(this->n_rows(), B.n_rows());
+  Assert(B.state == matrix, ExcState(B.state));
 
+  Assert(this->n_rows() == this->n_cols(),
+         LACExceptions::ExcNotQuadratic());
+  AssertDimension(this->n_rows(), B.n_rows());
   const char *trans = transposed ? &T : &N;
   const int nn = this->n_cols();
-  const int kk = B.n_cols();
   const number *values = &this->values[0];
+  const int n_rhs = B.n_cols();
   int info = 0;
 
-  getrs(trans, &nn, &kk, values, &nn, ipiv.data(), &B.values[0], &nn, &info);
+  if (state == lu)
+    {
+      getrs(trans, &nn, &n_rhs, values, &nn, ipiv.data(),
+            &B.values[0], &nn, &info);
+    }
+  else if (state == cholesky)
+    {
+      potrs(&LAPACKSupport::L, &nn, &n_rhs,
+            values, &nn, &B.values[0], &nn, &info);
+    }
+  else if (property == upper_triangular ||
+           property == lower_triangular)
+    {
+      const char uplo = (property == upper_triangular ? LAPACKSupport::U : LAPACKSupport::L);
+
+      const int lda = nn;
+      const int ldb = nn;
+      trtrs (&uplo, trans, "N",
+             &nn, &n_rhs,
+             values, &lda,
+             &B.values[0], &ldb, &info);
+    }
+  else
+    {
+      Assert(false,
+             ExcMessage("The matrix has to be either factorized or triangular."));
+    }
 
   Assert(info == 0, ExcInternalError());
 }
 
 
+
+template <typename number>
+void
+LAPACKFullMatrix<number>::apply_lu_factorization(Vector<number> &v,
+                                                 const bool transposed) const
+{
+  solve(v,transposed);
+}
+
+
+
+template <typename number>
+void
+LAPACKFullMatrix<number>::apply_lu_factorization(LAPACKFullMatrix<number> &B,
+                                                 const bool transposed) const
+{
+  solve(B,transposed);
+}
+
+
+
 template <typename number>
 number
 LAPACKFullMatrix<number>::determinant() const
diff --git a/tests/lapack/full_matrix_24.cc b/tests/lapack/full_matrix_24.cc
new file mode 100644 (file)
index 0000000..72a63c0
--- /dev/null
@@ -0,0 +1,90 @@
+// ---------------------------------------------------------------------
+//
+// 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::solve() for triangular matrices
+
+/* MWE for size=3 in Octave:
+R = [1,2,3; 0, 5, 6; 0, 0, 9]
+x = [2; -7; 1]
+
+> R\x
+ans =
+
+   4.73333
+  -1.53333
+   0.11111
+
+> R'\x
+ans =
+
+   2.00000
+  -2.20000
+   0.91111
+*/
+
+#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 = 3;
+  LAPACKFullMatrix<NumberType> M(size);
+  M.set_property(LAPACKSupport::upper_triangular);
+
+  M = 0.;
+  unsigned int counter = 1;
+  for (unsigned int i = 0; i < size; ++i)
+    for (unsigned int j = 0; j < size; ++j)
+      {
+        if (j >= i)
+          M(i,j) = counter;
+
+        counter++;
+      }
+
+  Vector<NumberType> x(size), y(size);
+  x[0] = 2;
+  x[1] = -7;
+  x[2] = 1;
+
+  y = x;
+  M.solve(y, false);
+  y.print(deallog.get_file_stream(), 6, false);
+
+  y = x;
+  M.solve(y, true);
+  y.print(deallog.get_file_stream(), 6, false);
+}
+
+
+int main()
+{
+  const std::string logname = "output";
+  std::ofstream logfile(logname.c_str());
+  logfile.precision(3);
+  deallog.attach(logfile);
+
+  test<double>();
+
+}
diff --git a/tests/lapack/full_matrix_24.output b/tests/lapack/full_matrix_24.output
new file mode 100644 (file)
index 0000000..551292b
--- /dev/null
@@ -0,0 +1,3 @@
+
+4.733333 -1.533333 0.111111 
+2.000000 -2.200000 0.911111 

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.