From b023304b10a95dd845208b156388628b51aed295 Mon Sep 17 00:00:00 2001 From: Sambit Das Date: Thu, 7 Jun 2018 02:59:20 -0400 Subject: [PATCH] Extend ScaLAPACKMatrix::invert() to use pXtrtri for inversion of triangular matrices. --- doc/news/changes/minor/20180606SambitDas | 3 + include/deal.II/lac/scalapack.h | 3 +- include/deal.II/lac/scalapack.templates.h | 67 +++++++ source/lac/scalapack.cc | 175 +++++++++++------- tests/lapack/create_matrix.h | 18 +- tests/scalapack/scalapack_03_c.cc | 125 +++++++++++++ .../scalapack/scalapack_03_c.mpirun=1.output | 36 ++++ .../scalapack/scalapack_03_c.mpirun=10.output | 36 ++++ .../scalapack/scalapack_03_c.mpirun=4.output | 36 ++++ 9 files changed, 434 insertions(+), 65 deletions(-) create mode 100644 doc/news/changes/minor/20180606SambitDas create mode 100644 tests/scalapack/scalapack_03_c.cc create mode 100644 tests/scalapack/scalapack_03_c.mpirun=1.output create mode 100644 tests/scalapack/scalapack_03_c.mpirun=10.output create mode 100644 tests/scalapack/scalapack_03_c.mpirun=4.output diff --git a/doc/news/changes/minor/20180606SambitDas b/doc/news/changes/minor/20180606SambitDas new file mode 100644 index 0000000000..ac62bd9895 --- /dev/null +++ b/doc/news/changes/minor/20180606SambitDas @@ -0,0 +1,3 @@ +Improved: Extend ScaLAPACKMatrix::invert() to use pXtrtri for inversion of triangular matrices. +
+(Sambit Das 2018/06/06) diff --git a/include/deal.II/lac/scalapack.h b/include/deal.II/lac/scalapack.h index b5ae7283cd..fd8f0d40c2 100644 --- a/include/deal.II/lac/scalapack.h +++ b/include/deal.II/lac/scalapack.h @@ -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 pXpotri or - * pXgetri. + * pXgetri. If the matrix is triangular, the LU factorization + * step is skipped, and pXtrtri is used directly. * * If a Cholesky or LU factorization has been applied previously, * pXpotri or pXgetri are called directly. diff --git a/include/deal.II/lac/scalapack.templates.h b/include/deal.II/lac/scalapack.templates.h index 46d6748bf5..fa47c6ceee 100644 --- a/include/deal.II/lac/scalapack.templates.h +++ b/include/deal.II/lac/scalapack.templates.h @@ -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 +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 inline void diff --git a/source/lac/scalapack.cc b/source/lac/scalapack.cc index 24c5d765af..4e7badf55b 100644 --- a/source/lac/scalapack.cc +++ b/source/lac/scalapack.cc @@ -309,9 +309,9 @@ ScaLAPACKMatrix::copy_to(FullMatrix &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::copy_to(FullMatrix &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::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; diff --git a/tests/lapack/create_matrix.h b/tests/lapack/create_matrix.h index d51dca4ffa..c3fc9aeb11 100644 --- a/tests/lapack/create_matrix.h +++ b/tests/lapack/create_matrix.h @@ -41,7 +41,23 @@ create_spd(FullMatrix &A) } } - +// create random invertible lower triangular matrix +template +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(); + } +} template void diff --git a/tests/scalapack/scalapack_03_c.cc b/tests/scalapack/scalapack_03_c.cc new file mode 100644 index 0000000000..47da17eee5 --- /dev/null +++ b/tests/scalapack/scalapack_03_c.cc @@ -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 +#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 random lower triangular matrices of requested size: + FullMatrix full_in(size), inverse(size), full_out(size), + diff(size), prod1(size), prod2(size), one(size); + + std::shared_ptr grid = + std::make_shared( + mpi_communicator, size, size, block_size, block_size); + ScaLAPACKMatrix 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::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_c.mpirun=1.output b/tests/scalapack/scalapack_03_c.mpirun=1.output new file mode 100644 index 0000000000..21fc377823 --- /dev/null +++ b/tests/scalapack/scalapack_03_c.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_c.mpirun=10.output b/tests/scalapack/scalapack_03_c.mpirun=10.output new file mode 100644 index 0000000000..384c6b7acd --- /dev/null +++ b/tests/scalapack/scalapack_03_c.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_c.mpirun=4.output b/tests/scalapack/scalapack_03_c.mpirun=4.output new file mode 100644 index 0000000000..827afa216b --- /dev/null +++ b/tests/scalapack/scalapack_03_c.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 -- 2.39.5