From: Benjamin Brands Date: Tue, 20 Mar 2018 22:34:54 +0000 (+0100) Subject: ScaLAPACKMatrix::invert() for non-symmetric matrices X-Git-Tag: v9.0.0-rc1~275^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=refs%2Fpull%2F6080%2Fhead;p=dealii.git ScaLAPACKMatrix::invert() for non-symmetric matrices --- diff --git a/include/deal.II/lac/scalapack.h b/include/deal.II/lac/scalapack.h index 19dca92b6e..0356234084 100644 --- a/include/deal.II/lac/scalapack.h +++ b/include/deal.II/lac/scalapack.h @@ -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 pXpotri. The inverse is stored - * in this object. + * Compute the LU factorization of the matrix using ScaLAPACK + * function pXgetrf 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 pXpotri or pXgetri. + * + * If a Cholesky or LU factorization has been applied previously, + * pXpotri or pXgetri are called directly. + * + * The inverse is stored in this object. */ void invert(); @@ -695,6 +707,12 @@ private: */ mutable std::vector iwork; + /** + * Integer array holding pivoting information required + * by ScaLAPACK's matrix factorization routines. + */ + std::vector ipiv; + /** * A character to define where elements are stored in case * ScaLAPACK operations support this. diff --git a/include/deal.II/lac/scalapack.templates.h b/include/deal.II/lac/scalapack.templates.h index 824c7f1985..00db91dee4 100644 --- a/include/deal.II/lac/scalapack.templates.h +++ b/include/deal.II/lac/scalapack.templates.h @@ -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 +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 inline void ppotri(const char * /*UPLO*/, const int * /*N*/, @@ -891,6 +984,53 @@ inline void ppotri(const char *UPLO, } +template +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 inline void ppocon(const char * /*uplo*/, const int * /*N*/, diff --git a/source/lac/scalapack.cc b/source/lac/scalapack.cc index 6ea55a6f96..4c978a8d08 100644 --- a/source/lac/scalapack.cc +++ b/source/lac/scalapack.cc @@ -656,8 +656,10 @@ void ScaLAPACKMatrix::TmTmult(ScaLAPACKMatrix &C, template void ScaLAPACKMatrix::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::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 -void ScaLAPACKMatrix::invert() +void ScaLAPACKMatrix::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 +void ScaLAPACKMatrix::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 index 0000000000..44fed3fcbb --- /dev/null +++ b/tests/scalapack/scalapack_03_b.cc @@ -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 +#include +#include +#include +#include + +#include + +#include + +#include +#include + +template +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 full_in(size), inverse(size), full_out(size), diff(size), prod1(size), prod2(size), one(size); + create_spd (full_in); + + std::shared_ptr grid = std::make_shared(mpi_communicator,size,size,block_size,block_size); + ScaLAPACKMatrix 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::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 sizes = {{32,64,120,320,640}}; + const std::vector blocks = {{32,64}}; + + for (const auto &s : sizes) + for (const auto &b : blocks) + if (b <= s) + test(s,b); + + for (const auto &s : sizes) + for (const auto &b : blocks) + if (b <= s) + test(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 index 0000000000..21fc377823 --- /dev/null +++ b/tests/scalapack/scalapack_03_b.mpirun=1.output @@ -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 index 0000000000..384c6b7acd --- /dev/null +++ b/tests/scalapack/scalapack_03_b.mpirun=10.output @@ -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 index 0000000000..827afa216b --- /dev/null +++ b/tests/scalapack/scalapack_03_b.mpirun=4.output @@ -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