]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Extend ScaLAPACKMatrix::invert() to use pXtrtri for inversion of triangular matrices. 6715/head
authorSambit Das <dsambit@nyx1231.arc-ts.umich.edu>
Thu, 7 Jun 2018 06:59:20 +0000 (02:59 -0400)
committerSambit Das <dsambit@nyx1231.arc-ts.umich.edu>
Thu, 7 Jun 2018 06:59:20 +0000 (02:59 -0400)
doc/news/changes/minor/20180606SambitDas [new file with mode: 0644]
include/deal.II/lac/scalapack.h
include/deal.II/lac/scalapack.templates.h
source/lac/scalapack.cc
tests/lapack/create_matrix.h
tests/scalapack/scalapack_03_c.cc [new file with mode: 0644]
tests/scalapack/scalapack_03_c.mpirun=1.output [new file with mode: 0644]
tests/scalapack/scalapack_03_c.mpirun=10.output [new file with mode: 0644]
tests/scalapack/scalapack_03_c.mpirun=4.output [new file with mode: 0644]

diff --git a/doc/news/changes/minor/20180606SambitDas b/doc/news/changes/minor/20180606SambitDas
new file mode 100644 (file)
index 0000000..ac62bd9
--- /dev/null
@@ -0,0 +1,3 @@
+Improved: Extend ScaLAPACKMatrix::invert() to use pXtrtri for inversion of triangular matrices.
+<br>
+(Sambit Das 2018/06/06)
index b5ae7283cd1b5774df5dcb3a87c1a79dd3a71e04..fd8f0d40c24f630f743f7a6f80055c117bfcb71e 100644 (file)
@@ -473,7 +473,8 @@ public:
    * Invert the matrix by first computing a Cholesky for symmetric matrices
    * or a LU factorization for general matrices and then
    * building the actual inverse using <code>pXpotri</code> or
-   * <code>pXgetri</code>.
+   * <code>pXgetri</code>. If the matrix is triangular, the LU factorization
+   * step is skipped, and <code>pXtrtri</code> is used directly.
    *
    * If a Cholesky or LU factorization has been applied previously,
    * <code>pXpotri</code> or <code>pXgetri</code> are called directly.
index 46d6748bf54bf317602a35aa70e65e7f9a060eba..fa47c6ceee4f5d4f4a94c021510e31a4652e05aa 100644 (file)
@@ -271,6 +271,34 @@ extern "C"
            int *      liwork,
            int *      info);
 
+
+  /**
+   * PDTRTRI computes the inverse of a upper or lower triangular
+   * distributed matrix sub( A ) = A(IA:IA+N-1,JA:JA+N-1).
+   *
+   * http://www.netlib.org/scalapack/explore-html/d9/dc0/pdtrtri_8f_source.html
+   * https://www.ibm.com/support/knowledgecenter/SSNR5K_4.2.0/com.ibm.cluster.pessl.v4r2.pssl100.doc/am6gr_lpdtri.htm
+   * https://software.intel.com/en-us/mkl-developer-reference-c-p-trtri
+   */
+  void
+  pdtrtri_(const char *UPLO,
+           const char *DIAG,
+           const int * N,
+           double *    A,
+           const int * IA,
+           const int * JA,
+           const int * DESCA,
+           int *       INFO);
+  void
+  pstrtri_(const char *UPLO,
+           const char *DIAG,
+           const int * N,
+           float *     A,
+           const int * IA,
+           const int * JA,
+           const int * DESCA,
+           int *       INFO);
+
   /**
    * Estimate the reciprocal of the condition number (in the
    * l1-norm) of a real symmetric positive definite distributed matrix
@@ -1117,6 +1145,45 @@ pgetri(const int *N,
   psgetri_(N, A, IA, JA, DESCA, ipiv, work, lwork, iwork, liwork, info);
 }
 
+template <typename number>
+inline void
+ptrtri(const char * /*UPLO*/,
+       const char * /*DIAG*/,
+       const int * /*N*/,
+       number * /*A*/,
+       const int * /*IA*/,
+       const int * /*JA*/,
+       const int * /*DESCA*/,
+       int * /*INFO*/)
+{
+  Assert(false, dealii::ExcNotImplemented());
+}
+
+inline void
+ptrtri(const char *UPLO,
+       const char *DIAG,
+       const int * N,
+       double *    A,
+       const int * IA,
+       const int * JA,
+       const int * DESCA,
+       int *       INFO)
+{
+  pdtrtri_(UPLO, DIAG, N, A, IA, JA, DESCA, INFO);
+}
+
+inline void
+ptrtri(const char *UPLO,
+       const char *DIAG,
+       const int * N,
+       float *     A,
+       const int * IA,
+       const int * JA,
+       const int * DESCA,
+       int *       INFO)
+{
+  pstrtri_(UPLO, DIAG, N, A, IA, JA, DESCA, INFO);
+}
 
 template <typename number>
 inline void
index 24c5d765af43fc11d91132cf6fb42497cbbd9f14..4e7badf55b71cd91f5622598fe39f9e8785a6a84 100644 (file)
@@ -309,9 +309,9 @@ ScaLAPACKMatrix<NumberType>::copy_to(FullMatrix<NumberType> &matrix) const
   Assert(n_columns == int(matrix.n()),
          ExcDimensionMismatch(n_columns, matrix.n()));
 
+  matrix = 0.;
   if (grid->mpi_process_is_active)
     {
-      matrix = 0.;
       for (int i = 0; i < n_local_rows; ++i)
         {
           const int glob_i = global_row(i);
@@ -327,16 +327,29 @@ ScaLAPACKMatrix<NumberType>::copy_to(FullMatrix<NumberType> &matrix) const
   // we could move the following lines under the main loop above,
   // but they would be dependent on glob_i and glob_j, which
   // won't make it much prettier
-  if (property == LAPACKSupport::lower_triangular)
-    for (unsigned int i = 0; i < matrix.n(); ++i)
-      for (unsigned int j = i + 1; j < matrix.m(); ++j)
-        matrix(i, j) =
-          (state == LAPACKSupport::inverse_matrix ? matrix(j, i) : 0.);
-  else if (property == LAPACKSupport::upper_triangular)
-    for (unsigned int i = 0; i < matrix.n(); ++i)
-      for (unsigned int j = 0; j < i; ++j)
-        matrix(i, j) =
-          (state == LAPACKSupport::inverse_matrix ? matrix(j, i) : 0.);
+  if (state == LAPACKSupport::cholesky)
+    {
+      if (property == LAPACKSupport::lower_triangular)
+        for (unsigned int i = 0; i < matrix.n(); ++i)
+          for (unsigned int j = i + 1; j < matrix.m(); ++j)
+            matrix(i, j) = 0.;
+      else if (property == LAPACKSupport::upper_triangular)
+        for (unsigned int i = 0; i < matrix.n(); ++i)
+          for (unsigned int j = 0; j < i; ++j)
+            matrix(i, j) = 0.;
+    }
+  else if (property == LAPACKSupport::symmetric &&
+           state == LAPACKSupport::inverse_matrix)
+    {
+      if (uplo == 'L')
+        for (unsigned int i = 0; i < matrix.n(); ++i)
+          for (unsigned int j = i + 1; j < matrix.m(); ++j)
+            matrix(i, j) = matrix(j, i);
+      else if (uplo == 'U')
+        for (unsigned int i = 0; i < matrix.n(); ++i)
+          for (unsigned int j = 0; j < i; ++j)
+            matrix(i, j) = matrix(j, i);
+    }
 }
 
 
@@ -943,69 +956,105 @@ ScaLAPACKMatrix<NumberType>::invert()
   const bool is_symmetric = (property == LAPACKSupport::symmetric ||
                              state == LAPACKSupport::State::cholesky);
 
-  // Matrix is neither in Cholesky nor LU state.
-  // Compute the required factorizations based on the property of the matrix.
-  if (!(state == LAPACKSupport::State::lu ||
-        state == LAPACKSupport::State::cholesky))
-    {
-      if (is_symmetric)
-        compute_cholesky_factorization();
-      else
-        compute_lu_factorization();
-    }
-  if (grid->mpi_process_is_active)
-    {
-      int         info  = 0;
-      NumberType *A_loc = &this->values[0];
+  // Check whether matrix is triangular and is in an unfactorized state.
+  const bool is_triangular = (property == LAPACKSupport::upper_triangular ||
+                              property == LAPACKSupport::lower_triangular) &&
+                             (state == LAPACKSupport::State::matrix ||
+                              state == LAPACKSupport::State::inverse_matrix);
 
-      if (is_symmetric)
+  if (is_triangular)
+    {
+      if (grid->mpi_process_is_active)
         {
-          ppotri(&uplo,
+          const char uploTriangular =
+            property == LAPACKSupport::upper_triangular ? 'U' : 'L';
+          const char  diag  = 'N';
+          int         info  = 0;
+          NumberType *A_loc = &this->values[0];
+          ptrtri(&uploTriangular,
+                 &diag,
                  &n_columns,
                  A_loc,
                  &submatrix_row,
                  &submatrix_column,
                  descriptor,
                  &info);
-          AssertThrow(info == 0, LAPACKSupport::ExcErrorCode("ppotri", info));
+          AssertThrow(info == 0, LAPACKSupport::ExcErrorCode("ptrtri", info));
+          // The inversion is stored in the same part as the triangular matrix,
+          // so we don't need to re-set the property here.
         }
-      else
+    }
+  else
+    {
+      // Matrix is neither in Cholesky nor LU state.
+      // Compute the required factorizations based on the property of the
+      // matrix.
+      if (!(state == LAPACKSupport::State::lu ||
+            state == LAPACKSupport::State::cholesky))
         {
-          int lwork = -1, liwork = -1;
-          work.resize(1);
-          iwork.resize(1);
-
-          pgetri(&n_columns,
-                 A_loc,
-                 &submatrix_row,
-                 &submatrix_column,
-                 descriptor,
-                 ipiv.data(),
-                 work.data(),
-                 &lwork,
-                 iwork.data(),
-                 &liwork,
-                 &info);
-
-          AssertThrow(info == 0, LAPACKSupport::ExcErrorCode("pgetri", info));
-          lwork  = work[0];
-          liwork = iwork[0];
-          work.resize(lwork);
-          iwork.resize(liwork);
-
-          pgetri(&n_columns,
-                 A_loc,
-                 &submatrix_row,
-                 &submatrix_column,
-                 descriptor,
-                 ipiv.data(),
-                 work.data(),
-                 &lwork,
-                 iwork.data(),
-                 &liwork,
-                 &info);
+          if (is_symmetric)
+            compute_cholesky_factorization();
+          else
+            compute_lu_factorization();
+        }
+      if (grid->mpi_process_is_active)
+        {
+          int         info  = 0;
+          NumberType *A_loc = &this->values[0];
 
-          AssertThrow(info == 0, LAPACKSupport::ExcErrorCode("pgetri", info));
+          if (is_symmetric)
+            {
+              ppotri(&uplo,
+                     &n_columns,
+                     A_loc,
+                     &submatrix_row,
+                     &submatrix_column,
+                     descriptor,
+                     &info);
+              AssertThrow(info == 0,
+                          LAPACKSupport::ExcErrorCode("ppotri", info));
+              property = LAPACKSupport::Property::symmetric;
+            }
+          else
+            {
+              int lwork = -1, liwork = -1;
+              work.resize(1);
+              iwork.resize(1);
+
+              pgetri(&n_columns,
+                     A_loc,
+                     &submatrix_row,
+                     &submatrix_column,
+                     descriptor,
+                     ipiv.data(),
+                     work.data(),
+                     &lwork,
+                     iwork.data(),
+                     &liwork,
+                     &info);
+
+              AssertThrow(info == 0,
+                          LAPACKSupport::ExcErrorCode("pgetri", info));
+              lwork  = work[0];
+              liwork = iwork[0];
+              work.resize(lwork);
+              iwork.resize(liwork);
+
+              pgetri(&n_columns,
+                     A_loc,
+                     &submatrix_row,
+                     &submatrix_column,
+                     descriptor,
+                     ipiv.data(),
+                     work.data(),
+                     &lwork,
+                     iwork.data(),
+                     &liwork,
+                     &info);
+
+              AssertThrow(info == 0,
+                          LAPACKSupport::ExcErrorCode("pgetri", info));
+            }
         }
     }
   state = LAPACKSupport::State::inverse_matrix;
index d51dca4ffa45e88bab7eef32f4c54e9c568323eb..c3fc9aeb116e88e4eb8781154a3280d7ce59f982 100644 (file)
@@ -41,7 +41,23 @@ create_spd(FullMatrix &A)
       }
 }
 
-
+// create random invertible lower triangular matrix
+template <typename FullMatrix>
+void
+create_random_lt(FullMatrix &A)
+{
+  const unsigned int size = A.n();
+  Assert(size == A.m(), ExcDimensionMismatch(size, A.m()));
+  A = 0.;
+  for (unsigned int i = 0; i < size; ++i)
+    for (unsigned int j = 0; j <= i; ++j)
+      {
+        if (i == j)
+          A(i, j) = (typename FullMatrix::value_type)(1.);
+        else
+          A(i, j) = random_value<typename FullMatrix::value_type>();
+      }
+}
 
 template <typename FullMatrix>
 void
diff --git a/tests/scalapack/scalapack_03_c.cc b/tests/scalapack/scalapack_03_c.cc
new file mode 100644 (file)
index 0000000..47da17e
--- /dev/null
@@ -0,0 +1,125 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2017 - 2018 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.
+//
+// ---------------------------------------------------------------------
+
+#include "../lapack/create_matrix.h"
+#include "../tests.h"
+
+/*
+ * test inverse of triangular matrix using pXtrtri in ScaLAPACK vs FullMatrix
+ */
+
+#include <deal.II/base/conditional_ostream.h>
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/multithread_info.h>
+#include <deal.II/base/timer.h>
+#include <deal.II/base/utilities.h>
+
+#include <deal.II/lac/scalapack.h>
+#include <deal.II/lac/vector.h>
+
+#include <fstream>
+#include <iostream>
+
+template <typename NumberType>
+void
+test(const unsigned int size, const unsigned int block_size)
+{
+  MPI_Comm           mpi_communicator(MPI_COMM_WORLD);
+  const unsigned int n_mpi_processes(
+    Utilities::MPI::n_mpi_processes(mpi_communicator));
+  const unsigned int this_mpi_process(
+    Utilities::MPI::this_mpi_process(mpi_communicator));
+
+  ConditionalOStream pcout(std::cout, (this_mpi_process == 0));
+
+  // Create random lower triangular matrices of requested size:
+  FullMatrix<NumberType> full_in(size), inverse(size), full_out(size),
+    diff(size), prod1(size), prod2(size), one(size);
+
+  std::shared_ptr<Utilities::MPI::ProcessGrid> grid =
+    std::make_shared<Utilities::MPI::ProcessGrid>(
+      mpi_communicator, size, size, block_size, block_size);
+  ScaLAPACKMatrix<NumberType> scalapack_matrix(
+    size, grid, block_size, LAPACKSupport::Property::lower_triangular);
+
+  pcout << size << " " << block_size << " " << grid->get_process_grid_rows()
+        << " " << grid->get_process_grid_columns() << std::endl;
+
+  create_random_lt(full_in);
+
+  one = 0.;
+  for (unsigned int i = 0; i < size; ++i)
+    one(i, i) = 1.;
+  // invert via Lapack
+  inverse.invert(full_in);
+  inverse.mmult(prod1, full_in);
+  prod1.add(-1., one);
+  const NumberType lapack_error = prod1.linfty_norm();
+
+  // estimated condition number from 1-norm:
+  const NumberType k   = full_in.l1_norm() * inverse.l1_norm();
+  const NumberType tol = k * 1000 * std::numeric_limits<NumberType>::epsilon();
+
+  // invert via ScaLAPACK
+  scalapack_matrix = full_in;
+  scalapack_matrix.invert();
+  scalapack_matrix.copy_to(full_out);
+  full_out.mmult(prod2, full_in);
+  prod2.add(-1., one);
+  const NumberType error = prod2.linfty_norm();
+
+  if (error > tol && this_mpi_process == 0)
+    {
+      diff = 0;
+      diff.add(1., inverse);
+      diff.add(-1., full_out);
+
+      std::cout << "Norm of the error " << error
+                << " is more than the threshold " << tol
+                << " . Norm of the A^{-1}*A using Lapack is " << lapack_error
+                << std::endl
+                << "===== Expected to have:" << std::endl;
+      inverse.print_formatted(std::cout);
+      std::cout << "===== But got:" << std::endl;
+      full_out.print_formatted(std::cout);
+      std::cout << "===== Difference:" << std::endl;
+      diff.print_formatted(std::cout);
+      AssertThrow(false, dealii::ExcInternalError());
+    }
+  else
+    pcout << "Ok" << std::endl;
+}
+
+
+
+int
+main(int argc, char **argv)
+{
+  Utilities::MPI::MPI_InitFinalize mpi_initialization(
+    argc, argv, numbers::invalid_unsigned_int);
+
+  const std::vector<unsigned int> sizes  = {{32, 64, 120, 320, 640}};
+  const std::vector<unsigned int> blocks = {{32, 64}};
+
+  for (const auto &s : sizes)
+    for (const auto &b : blocks)
+      if (b <= s)
+        test<float>(s, b);
+
+  for (const auto &s : sizes)
+    for (const auto &b : blocks)
+      if (b <= s)
+        test<double>(s, b);
+}
diff --git a/tests/scalapack/scalapack_03_c.mpirun=1.output b/tests/scalapack/scalapack_03_c.mpirun=1.output
new file mode 100644 (file)
index 0000000..21fc377
--- /dev/null
@@ -0,0 +1,36 @@
+32 32 1 1
+Ok
+64 32 1 1
+Ok
+64 64 1 1
+Ok
+120 32 1 1
+Ok
+120 64 1 1
+Ok
+320 32 1 1
+Ok
+320 64 1 1
+Ok
+640 32 1 1
+Ok
+640 64 1 1
+Ok
+32 32 1 1
+Ok
+64 32 1 1
+Ok
+64 64 1 1
+Ok
+120 32 1 1
+Ok
+120 64 1 1
+Ok
+320 32 1 1
+Ok
+320 64 1 1
+Ok
+640 32 1 1
+Ok
+640 64 1 1
+Ok
diff --git a/tests/scalapack/scalapack_03_c.mpirun=10.output b/tests/scalapack/scalapack_03_c.mpirun=10.output
new file mode 100644 (file)
index 0000000..384c6b7
--- /dev/null
@@ -0,0 +1,36 @@
+32 32 1 1
+Ok
+64 32 2 2
+Ok
+64 64 1 1
+Ok
+120 32 3 3
+Ok
+120 64 2 2
+Ok
+320 32 3 3
+Ok
+320 64 3 3
+Ok
+640 32 3 3
+Ok
+640 64 3 3
+Ok
+32 32 1 1
+Ok
+64 32 2 2
+Ok
+64 64 1 1
+Ok
+120 32 3 3
+Ok
+120 64 2 2
+Ok
+320 32 3 3
+Ok
+320 64 3 3
+Ok
+640 32 3 3
+Ok
+640 64 3 3
+Ok
diff --git a/tests/scalapack/scalapack_03_c.mpirun=4.output b/tests/scalapack/scalapack_03_c.mpirun=4.output
new file mode 100644 (file)
index 0000000..827afa2
--- /dev/null
@@ -0,0 +1,36 @@
+32 32 1 1
+Ok
+64 32 2 2
+Ok
+64 64 1 1
+Ok
+120 32 2 2
+Ok
+120 64 2 2
+Ok
+320 32 2 2
+Ok
+320 64 2 2
+Ok
+640 32 2 2
+Ok
+640 64 2 2
+Ok
+32 32 1 1
+Ok
+64 32 2 2
+Ok
+64 64 1 1
+Ok
+120 32 2 2
+Ok
+120 64 2 2
+Ok
+320 32 2 2
+Ok
+320 64 2 2
+Ok
+640 32 2 2
+Ok
+640 64 2 2
+Ok

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.