]> https://gitweb.dealii.org/ - dealii.git/commitdiff
ScaLAPACKMatrix::invert() for non-symmetric matrices 6080/head
authorBenjamin Brands <benjamin.brands@fau.de>
Tue, 20 Mar 2018 22:34:54 +0000 (23:34 +0100)
committerBenjamin Brands <benjamin.brands@fau.de>
Wed, 21 Mar 2018 12:42:45 +0000 (13:42 +0100)
include/deal.II/lac/scalapack.h
include/deal.II/lac/scalapack.templates.h
source/lac/scalapack.cc
tests/scalapack/scalapack_03_b.cc [new file with mode: 0644]
tests/scalapack/scalapack_03_b.mpirun=1.output [new file with mode: 0644]
tests/scalapack/scalapack_03_b.mpirun=10.output [new file with mode: 0644]
tests/scalapack/scalapack_03_b.mpirun=4.output [new file with mode: 0644]

index 19dca92b6e52bc74cd64e5f4e81c9610c4344d30..0356234084b5851f63e9283d30fe83d235aad591 100644 (file)
@@ -387,9 +387,21 @@ public:
   void compute_cholesky_factorization ();
 
   /**
-   * Invert the matrix by first computing a Cholesky factorization and then
-   * building the actual inverse using <code>pXpotri</code>. The inverse is stored
-   * in this object.
+   * Compute the LU factorization of the matrix using ScaLAPACK
+   * function <code>pXgetrf</code> and partial pivoting with row interchanges.
+   * The result of the factorization is stored in this object.
+   */
+  void compute_lu_factorization ();
+
+  /**
+   * 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>.
+   *
+   * If a Cholesky or LU factorization has been applied previously,
+   * <code>pXpotri</code> or <code>pXgetri</code> are called directly.
+   *
+   * The inverse is stored in this object.
    */
   void invert();
 
@@ -695,6 +707,12 @@ private:
    */
   mutable std::vector<int> iwork;
 
+  /**
+   * Integer array holding pivoting information required
+   * by ScaLAPACK's matrix factorization routines.
+   */
+  std::vector<int> ipiv;
+
   /**
    * A character to define where elements are stored in case
    * ScaLAPACK operations support this.
index 824c7f198526e772aed053f12058fe5c2a815885..00db91dee499e37678d478b5ee054a77f7a66339 100644 (file)
@@ -192,6 +192,29 @@ extern "C"
                 const int *DESCA,
                 int *INFO);
 
+  /**
+   * Computes an LU factorization of a general distributed matrix sub( A )
+   * using partial pivoting with row interchanges.
+   *
+   * http://www.netlib.org/scalapack/explore-html/df/dfe/pdgetrf_8f_source.html
+   * https://www.ibm.com/support/knowledgecenter/en/SSNR5K_4.2.0/com.ibm.cluster.pessl.v4r2.pssl100.doc/am6gr_lgetrf.htm
+   */
+  void pdgetrf_(const int *m,
+                const int *n,
+                double *A,
+                const int *IA,
+                const int *JA,
+                const int *DESCA,
+                int *ipiv,
+                int *INFO);
+  void psgetrf_(const int *m,
+                const int *n,
+                float *A,
+                const int *IA,
+                const int *JA,
+                const int *DESCA,
+                int *ipiv,
+                int *INFO);
 
   /**
    * Compute the inverse of a real symmetric positive definite
@@ -218,6 +241,38 @@ extern "C"
                 const int *DESCA,
                 int *INFO);
 
+  /**
+   * PDGETRI computes the inverse of a distributed matrix using the LU
+   * factorization computed by PDGETRF. This method inverts U and then
+   * computes the inverse of sub( A ) = A(IA:IA+N-1,JA:JA+N-1) denoted
+   * InvA by solving the system InvA*L = inv(U) for InvA.
+   *
+   * http://www.netlib.org/scalapack/explore-html/d3/df3/pdgetri_8f_source.html
+   * https://www.ibm.com/support/knowledgecenter/SSNR5K_4.2.0/com.ibm.cluster.pessl.v4r2.pssl100.doc/am6gr_lgetri.htm
+   */
+  void pdgetri_(const int *N,
+                double *A,
+                const int *IA,
+                const int *JA,
+                const int *DESCA,
+                const int *ipiv,
+                double *work,
+                int *lwork,
+                int *iwork,
+                int *liwork,
+                int *info);
+  void psgetri_(const int *N,
+                float *A,
+                const int *IA,
+                const int *JA,
+                const int *DESCA,
+                const int *ipiv,
+                float *work,
+                int *lwork,
+                int *iwork,
+                int *liwork,
+                int *info);
+
   /**
    * Estimate the reciprocal of the condition number (in the
    * l1-norm) of a real symmetric positive definite distributed matrix
@@ -856,6 +911,44 @@ inline void ppotrf(const char *UPLO,
 }
 
 
+template <typename number>
+inline void pgetrf(const int * /*m*/,
+                   const int * /*n*/,
+                   number * /*A*/,
+                   const int * /*IA*/,
+                   const int * /*JA*/,
+                   const int * /*DESCA*/,
+                   int * /*ipiv*/,
+                   int * /*INFO*/)
+{
+  Assert (false, dealii::ExcNotImplemented());
+}
+
+inline void pgetrf(const int *m,
+                   const int *n,
+                   double *A,
+                   const int *IA,
+                   const int *JA,
+                   const int *DESCA,
+                   int *ipiv,
+                   int *INFO)
+{
+  pdgetrf_(m,n,A,IA,JA,DESCA,ipiv,INFO);
+}
+
+inline void pgetrf(const int *m,
+                   const int *n,
+                   float *A,
+                   const int *IA,
+                   const int *JA,
+                   const int *DESCA,
+                   int *ipiv,
+                   int *INFO)
+{
+  psgetrf_(m,n,A,IA,JA,DESCA,ipiv,INFO);
+}
+
+
 template <typename number>
 inline void ppotri(const char * /*UPLO*/,
                    const int * /*N*/,
@@ -891,6 +984,53 @@ inline void ppotri(const char *UPLO,
 }
 
 
+template <typename number>
+inline void pgetri(const int * /*N*/,
+                   number * /*A*/,
+                   const int * /*IA*/,
+                   const int * /*JA*/,
+                   const int * /*DESCA*/,
+                   const int * /*ipiv*/,
+                   number * /*work*/,
+                   int * /*lwork*/,
+                   int * /*iwork*/,
+                   int * /*liwork*/,
+                   int * /*info*/)
+{
+  Assert (false, dealii::ExcNotImplemented());
+}
+
+inline void pgetri(const int *N,
+                   double *A,
+                   const int *IA,
+                   const int *JA,
+                   const int *DESCA,
+                   const int *ipiv,
+                   double *work,
+                   int *lwork,
+                   int *iwork,
+                   int *liwork,
+                   int *info)
+{
+  pdgetri_(N,A,IA,JA,DESCA,ipiv,work,lwork,iwork,liwork,info);
+}
+
+inline void pgetri(const int *N,
+                   float *A,
+                   const int *IA,
+                   const int *JA,
+                   const int *DESCA,
+                   const int *ipiv,
+                   float *work,
+                   int *lwork,
+                   int *iwork,
+                   int *liwork,
+                   int *info)
+{
+  psgetri_(N,A,IA,JA,DESCA,ipiv,work,lwork,iwork,liwork,info);
+}
+
+
 template <typename number>
 inline void ppocon(const char * /*uplo*/,
                    const int * /*N*/,
index 6ea55a6f9646d746405f7c44355db62f36765d31..4c978a8d08d3a37a05fd5d6209f797a2110ecf36 100644 (file)
@@ -656,8 +656,10 @@ void ScaLAPACKMatrix<NumberType>::TmTmult(ScaLAPACKMatrix<NumberType> &C,
 template <typename NumberType>
 void ScaLAPACKMatrix<NumberType>::compute_cholesky_factorization()
 {
-  Assert (n_columns == n_rows,
-          ExcMessage("Cholesky factorization can be applied to SPD matrices only."));
+  Assert (n_columns==n_rows && property == LAPACKSupport::Property::symmetric,
+          ExcMessage("Cholesky factorization can be applied to symmetric matrices only."));
+  Assert (state == LAPACKSupport::matrix,
+          ExcMessage("Matrix has to be in Matrix state before calling this function."));
 
   if (grid->mpi_process_is_active)
     {
@@ -667,26 +669,85 @@ void ScaLAPACKMatrix<NumberType>::compute_cholesky_factorization()
       ppotrf(&uplo,&n_columns,A_loc,&submatrix_row,&submatrix_column,descriptor,&info);
       AssertThrow (info==0, LAPACKSupport::ExcErrorCode("ppotrf", info));
     }
-  property = (uplo=='L' ? LAPACKSupport::lower_triangular : LAPACKSupport::upper_triangular);
   state = LAPACKSupport::cholesky;
+  property = (uplo=='L' ? LAPACKSupport::lower_triangular : LAPACKSupport::upper_triangular);
 }
 
 
 
 template <typename NumberType>
-void ScaLAPACKMatrix<NumberType>::invert()
+void ScaLAPACKMatrix<NumberType>::compute_lu_factorization()
 {
-  if (state == LAPACKSupport::matrix)
-    compute_cholesky_factorization();
+  Assert (state == LAPACKSupport::matrix,
+          ExcMessage("Matrix has to be in Matrix state before calling this function."));
 
   if (grid->mpi_process_is_active)
     {
       int info = 0;
       NumberType *A_loc = &this->values[0];
-      ppotri (&uplo,&n_columns, A_loc, &submatrix_row, &submatrix_column, descriptor,&info);
-      AssertThrow (info==0, LAPACKSupport::ExcErrorCode("ppotri", info));
+
+      const int iarow = indxg2p_(&submatrix_row, &row_block_size, &(grid->this_process_row), &first_process_row, &(grid->n_process_rows));
+      const int mp = numroc_(&n_rows, &row_block_size, &(grid->this_process_row), &iarow, &(grid->n_process_rows));
+      ipiv.resize(mp+row_block_size);
+
+      pgetrf(&n_rows,&n_columns,A_loc,&submatrix_row,&submatrix_column,descriptor,ipiv.data(),&info);
+      AssertThrow (info==0, LAPACKSupport::ExcErrorCode("pgetrf", info));
+    }
+  state = LAPACKSupport::State::lu;
+  property = LAPACKSupport::Property::general;
+}
+
+
+
+template <typename NumberType>
+void ScaLAPACKMatrix<NumberType>::invert()
+{
+  // Check whether matrix is symmetric and save flag.
+  // If a Cholesky factorization has been applied previously,
+  // the original matrix was symmetric.
+  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];
+
+      if (is_symmetric)
+        {
+          ppotri(&uplo, &n_columns, A_loc, &submatrix_row, &submatrix_column, descriptor, &info);
+          AssertThrow (info==0, LAPACKSupport::ExcErrorCode("ppotri", info));
+        }
+      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::inverse_matrix;
+  state = LAPACKSupport::State::inverse_matrix;
 }
 
 
diff --git a/tests/scalapack/scalapack_03_b.cc b/tests/scalapack/scalapack_03_b.cc
new file mode 100644 (file)
index 0000000..44fed3f
--- /dev/null
@@ -0,0 +1,130 @@
+// ---------------------------------------------------------------------
+//
+// 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 "../tests.h"
+#include "../lapack/create_matrix.h"
+
+/*
+ * test inverse via LU decomposition in ScaLAPACK vs FullMatrix for general matrices
+ *
+ * relevant tests:
+ * full_matrix/full_matrix_invert_01
+ * full_matrix/full_matrix_56
+ * full_matrix/full_matrix_53
+ * full_matrix/full_matrix_42
+ * full_matrix/full_matrix_41
+ */
+
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/utilities.h>
+#include <deal.II/base/conditional_ostream.h>
+#include <deal.II/base/timer.h>
+#include <deal.II/base/multithread_info.h>
+
+#include <deal.II/lac/vector.h>
+
+#include <deal.II/lac/scalapack.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 SPD matrices of requested size:
+  FullMatrix<NumberType> full_in(size), inverse(size), full_out(size), diff(size), prod1(size), prod2(size), one(size);
+  create_spd (full_in);
+
+  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::general);
+
+  pcout << size << " " << block_size << " " << grid->get_process_grid_rows() << " " << grid->get_process_grid_columns() << std::endl;
+
+  // add a skew-symmetric matrix to the s.p.d matrix
+  // the positive definiteness is inherited from the symmetric part
+  for (unsigned int i=1; i<size; ++i)
+    for (unsigned int j=i+1; j<size; ++j)
+      {
+        full_in(i,j) *= 1.3;
+        full_in(j,i) *= -1.3;
+      }
+  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_b.mpirun=1.output b/tests/scalapack/scalapack_03_b.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_b.mpirun=10.output b/tests/scalapack/scalapack_03_b.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_b.mpirun=4.output b/tests/scalapack/scalapack_03_b.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.