--- /dev/null
+Improved: Extend ScaLAPACKMatrix::invert() to use pXtrtri for inversion of triangular matrices.
+<br>
+(Sambit Das 2018/06/06)
* 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.
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
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
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);
// 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);
+ }
}
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;
}
}
-
+// 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
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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);
+}
--- /dev/null
+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
--- /dev/null
+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
--- /dev/null
+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
projects / dealii.git / commitdiff