From 93f3c810548e7f58b30715d833f74ab70dce0a8f Mon Sep 17 00:00:00 2001 From: Benjamin Brands Date: Tue, 6 Feb 2018 11:26:23 +0100 Subject: [PATCH] fix bug in computation of norms for ScaLAPACKMatrix --- include/deal.II/lac/scalapack.h | 8 +- include/deal.II/lac/scalapack.templates.h | 68 ++++++------ source/lac/scalapack.cc | 61 +++++++++-- tests/scalapack/scalapack_04_b.cc | 100 ++++++++++++++++++ .../scalapack/scalapack_04_b.mpirun=1.output | 36 +++++++ .../scalapack/scalapack_04_b.mpirun=11.output | 36 +++++++ .../scalapack/scalapack_04_b.mpirun=4.output | 36 +++++++ 7 files changed, 304 insertions(+), 41 deletions(-) create mode 100644 tests/scalapack/scalapack_04_b.cc create mode 100644 tests/scalapack/scalapack_04_b.mpirun=1.output create mode 100644 tests/scalapack/scalapack_04_b.mpirun=11.output create mode 100644 tests/scalapack/scalapack_04_b.mpirun=4.output diff --git a/include/deal.II/lac/scalapack.h b/include/deal.II/lac/scalapack.h index 42468c81fc..27e5a2b8ee 100644 --- a/include/deal.II/lac/scalapack.h +++ b/include/deal.II/lac/scalapack.h @@ -391,11 +391,17 @@ public: private: + /** + * Calculate the norm of a distributed symmetric dense matrix using ScaLAPACK's + * internal function. + */ + NumberType norm_symmetric(const char type) const; + /** * Calculate the norm of a distributed dense matrix using ScaLAPACK's * internal function. */ - NumberType norm(const char type) const; + NumberType norm_general(const char type) const; /** * Computing selected eigenvalues and, optionally, the eigenvectors. diff --git a/include/deal.II/lac/scalapack.templates.h b/include/deal.II/lac/scalapack.templates.h index 6ca40dfbe0..7479394c8e 100644 --- a/include/deal.II/lac/scalapack.templates.h +++ b/include/deal.II/lac/scalapack.templates.h @@ -403,20 +403,20 @@ extern "C" * or the element of largest absolute value of a distributed matrix */ double pdlange_(char const *norm, - int const &m, - int const &n, - double *A, - int const &ia, - int const &ja, - int *desca, + const int *m, + const int *n, + const double *A, + const int *ia, + const int *ja, + const int *desca, double *work); - float pslange_(char const *norm, - int const &m, - int const &n, - float *A, - int const &ia, - int const &ja, - int *desca, + float pslange_(const char *norm, + const int *m, + const int *n, + const float *A, + const int *ia, + const int *ja, + const int *desca, float *work); /** @@ -1043,37 +1043,37 @@ inline void pgemm(const char *transa, template -inline number plange(char const * /*norm*/, - int const &/*m*/, - int const &/*n*/, - number * /*A*/, - int const &/*ia*/, - int const &/*ja*/, - int * /*desca*/, +inline number plange(const char * /*norm*/, + const int * /*m*/, + const int * /*n*/, + const number * /*A*/, + const int * /*ia*/, + const int * /*ja*/, + const int * /*desca*/, number * /*work*/) { Assert (false, dealii::ExcNotImplemented()); } -inline double plange(char const *norm, - int const &m, - int const &n, - double *A, - int const &ia, - int const &ja, - int *desca, +inline double plange(const char *norm, + const int *m, + const int *n, + const double *A, + const int *ia, + const int *ja, + const int *desca, double *work) { return pdlange_(norm, m, n, A, ia, ja,desca, work); } -inline float plange(char const *norm, - int const &m, - int const &n, - float *A, - int const &ia, - int const &ja, - int *desca, +inline float plange(const char *norm, + const int *m, + const int *n, + const float *A, + const int *ia, + const int *ja, + const int *desca, float *work) { return pslange_(norm, m, n, A, ia, ja,desca, work); diff --git a/source/lac/scalapack.cc b/source/lac/scalapack.cc index 55d592b31a..236123895b 100644 --- a/source/lac/scalapack.cc +++ b/source/lac/scalapack.cc @@ -851,7 +851,11 @@ template NumberType ScaLAPACKMatrix::l1_norm() const { const char type('O'); - return norm(type); + + if (property == LAPACKSupport::symmetric) + return norm_symmetric(type); + else + return norm_general(type); } @@ -860,7 +864,11 @@ template NumberType ScaLAPACKMatrix::linfty_norm() const { const char type('I'); - return norm(type); + + if (property == LAPACKSupport::symmetric) + return norm_symmetric(type); + else + return norm_general(type); } @@ -869,17 +877,58 @@ template NumberType ScaLAPACKMatrix::frobenius_norm() const { const char type('F'); - return norm(type); + + if (property == LAPACKSupport::symmetric) + return norm_symmetric(type); + else + return norm_general(type); +} + + + +template +NumberType ScaLAPACKMatrix::norm_general(const char type) const +{ + Assert (state == LAPACKSupport::matrix || + state == LAPACKSupport::inverse_matrix, + ExcMessage("norms can be called in matrix state only.")); + Threads::Mutex::ScopedLock lock (mutex); + NumberType res = 0.; + + if (grid->mpi_process_is_active) + { + const int iarow = indxg2p_(&submatrix_row, &row_block_size, &(grid->this_process_row), &first_process_row, &(grid->n_process_rows)); + const int iacol = indxg2p_(&submatrix_column, &column_block_size, &(grid->this_process_column), &first_process_column, &(grid->n_process_columns)); + const int mp0 = numroc_(&n_rows, &row_block_size, &(grid->this_process_row), &iarow, &(grid->n_process_rows)); + const int nq0 = numroc_(&n_columns, &column_block_size, &(grid->this_process_column), &iacol, &(grid->n_process_columns)); + + // type='M': compute largest absolute value + // type='F' || type='E': compute Frobenius norm + // type='0' || type='1': compute infinity norm + int lwork=0; // for type == 'M' || type == 'F' || type == 'E' + if (type=='O' || type=='1') + lwork = nq0; + else if (type=='I') + lwork = mp0; + + work.resize(lwork); + const NumberType *A_loc = this->values.begin(); + res = plange(&type, &n_rows, &n_columns, A_loc, &submatrix_row, &submatrix_column, descriptor, work.data()); + } + grid->send_to_inactive(&res); + return res; } template -NumberType ScaLAPACKMatrix::norm(const char type) const +NumberType ScaLAPACKMatrix::norm_symmetric(const char type) const { Assert (state == LAPACKSupport::matrix || state == LAPACKSupport::inverse_matrix, ExcMessage("norms can be called in matrix state only.")); + Assert (property == LAPACKSupport::symmetric, + ExcMessage("Matrix has to be symmetric for this operation.")); Threads::Mutex::ScopedLock lock (mutex); NumberType res = 0.; @@ -904,8 +953,8 @@ NumberType ScaLAPACKMatrix::norm(const char type) const 0 : 2*Nq0+Np0+ldw; work.resize(lwork); - const NumberType *A_loc = &this->values[0]; - res = plansy(&type, &uplo, &n_columns, A_loc, &submatrix_row, &submatrix_column, descriptor, &work[0]); + const NumberType *A_loc = this->values.begin(); + res = plansy(&type, &uplo, &n_columns, A_loc, &submatrix_row, &submatrix_column, descriptor, work.data()); } grid->send_to_inactive(&res); return res; diff --git a/tests/scalapack/scalapack_04_b.cc b/tests/scalapack/scalapack_04_b.cc new file mode 100644 index 0000000000..f472085f24 --- /dev/null +++ b/tests/scalapack/scalapack_04_b.cc @@ -0,0 +1,100 @@ +// --------------------------------------------------------------------- +// +// 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 norms of ScaLAPACK vs FullMatrix for general matrices + +#include +#include +#include +#include +#include + +#include + +#include + +#include +#include + +template +void test(const unsigned int size, const unsigned int block_size, const NumberType tol) +{ + 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(size); + + std::shared_ptr grid = std::make_shared(mpi_communicator,size,size,block_size,block_size); + ScaLAPACKMatrix scalapack (size, size, grid, block_size, block_size, + LAPACKSupport::Property::general); + + pcout << size << " " << block_size << " " << grid->get_process_grid_rows() << " " << grid->get_process_grid_columns() << std::endl; + + create_random (full); + scalapack = full; + + const NumberType l1 = full.l1_norm(); + const NumberType linfty = full.linfty_norm(); + const NumberType frobenius = full.frobenius_norm(); + + // local result on this core: + const NumberType s_l1 = scalapack.l1_norm(); + const NumberType s_linfty = scalapack.linfty_norm(); + const NumberType s_frobenius = scalapack.frobenius_norm(); + + // make sure we have the same result on all cores, do average: + const NumberType as_l1 = dealii::Utilities::MPI::sum(s_l1, mpi_communicator)/n_mpi_processes; + const NumberType as_linfty = dealii::Utilities::MPI::sum(s_linfty, mpi_communicator)/n_mpi_processes; + const NumberType as_frobenius = dealii::Utilities::MPI::sum(s_frobenius, mpi_communicator)/n_mpi_processes; + + pcout << l1 << " " << s_l1 << " " << as_l1 << std::endl + << linfty << " " << s_linfty << " " << as_linfty << std::endl + << frobenius << " " << s_frobenius << " " << as_frobenius << std::endl; + + AssertThrow (std::abs(l1 -as_l1) < tol*std::abs(l1), dealii::ExcInternalError()); + AssertThrow (std::abs(linfty -as_linfty) < tol*std::abs(linfty), dealii::ExcInternalError()); + AssertThrow (std::abs(frobenius -as_frobenius) < tol*std::abs(frobenius), dealii::ExcInternalError()); +} + + + +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}}; + + const double tol_double = 1e-10; + const float tol_float = 1e-5; + + /*for (const auto &s : sizes) + for (const auto &b : blocks) + if (b <= s) + test(s,b,tol_float);*/ + + for (const auto &s : sizes) + for (const auto &b : blocks) + if (b <= s) + test(s,b,tol_double); +} diff --git a/tests/scalapack/scalapack_04_b.mpirun=1.output b/tests/scalapack/scalapack_04_b.mpirun=1.output new file mode 100644 index 0000000000..820f498747 --- /dev/null +++ b/tests/scalapack/scalapack_04_b.mpirun=1.output @@ -0,0 +1,36 @@ +32 32 1 1 +19.2741 19.2741 19.2741 +20.6796 20.6796 20.6796 +18.5785 18.5785 18.5785 +64 32 1 1 +38.5552 38.5552 38.5552 +38.2739 38.2739 38.2739 +36.8389 36.8389 36.8389 +64 64 1 1 +37.8671 37.8671 37.8671 +37.8657 37.8657 37.8657 +36.7171 36.7171 36.7171 +120 32 1 1 +69.0663 69.0663 69.0663 +68.0208 68.0208 68.0208 +69.1394 69.1394 69.1394 +120 64 1 1 +68.2078 68.2078 68.2078 +66.8542 66.8542 66.8542 +69.0112 69.0112 69.0112 +320 32 1 1 +174.387 174.387 174.387 +172.758 172.758 172.758 +184.876 184.876 184.876 +320 64 1 1 +182.948 182.948 182.948 +174.945 174.945 174.945 +184.521 184.521 184.521 +640 32 1 1 +338.935 338.935 338.935 +342.898 342.898 342.898 +369.65 369.65 369.65 +640 64 1 1 +342.421 342.421 342.421 +341.617 341.617 341.617 +369.497 369.497 369.497 diff --git a/tests/scalapack/scalapack_04_b.mpirun=11.output b/tests/scalapack/scalapack_04_b.mpirun=11.output new file mode 100644 index 0000000000..9f198a0d7d --- /dev/null +++ b/tests/scalapack/scalapack_04_b.mpirun=11.output @@ -0,0 +1,36 @@ +32 32 1 1 +19.2741 19.2741 19.2741 +20.6796 20.6796 20.6796 +18.5785 18.5785 18.5785 +64 32 2 2 +38.5552 38.5552 38.5552 +38.2739 38.2739 38.2739 +36.8389 36.8389 36.8389 +64 64 1 1 +37.8671 37.8671 37.8671 +37.8657 37.8657 37.8657 +36.7171 36.7171 36.7171 +120 32 3 3 +69.0663 69.0663 69.0663 +68.0208 68.0208 68.0208 +69.1394 69.1394 69.1394 +120 64 2 2 +68.2078 68.2078 68.2078 +66.8542 66.8542 66.8542 +69.0112 69.0112 69.0112 +320 32 3 3 +174.387 174.387 174.387 +172.758 172.758 172.758 +184.876 184.876 184.876 +320 64 3 3 +182.948 182.948 182.948 +174.945 174.945 174.945 +184.521 184.521 184.521 +640 32 3 3 +338.935 338.935 338.935 +342.898 342.898 342.898 +369.65 369.65 369.65 +640 64 3 3 +342.421 342.421 342.421 +341.617 341.617 341.617 +369.497 369.497 369.497 diff --git a/tests/scalapack/scalapack_04_b.mpirun=4.output b/tests/scalapack/scalapack_04_b.mpirun=4.output new file mode 100644 index 0000000000..603d0a4cb7 --- /dev/null +++ b/tests/scalapack/scalapack_04_b.mpirun=4.output @@ -0,0 +1,36 @@ +32 32 1 1 +19.2741 19.2741 19.2741 +20.6796 20.6796 20.6796 +18.5785 18.5785 18.5785 +64 32 2 2 +38.5552 38.5552 38.5552 +38.2739 38.2739 38.2739 +36.8389 36.8389 36.8389 +64 64 1 1 +37.8671 37.8671 37.8671 +37.8657 37.8657 37.8657 +36.7171 36.7171 36.7171 +120 32 2 2 +69.0663 69.0663 69.0663 +68.0208 68.0208 68.0208 +69.1394 69.1394 69.1394 +120 64 2 2 +68.2078 68.2078 68.2078 +66.8542 66.8542 66.8542 +69.0112 69.0112 69.0112 +320 32 2 2 +174.387 174.387 174.387 +172.758 172.758 172.758 +184.876 184.876 184.876 +320 64 2 2 +182.948 182.948 182.948 +174.945 174.945 174.945 +184.521 184.521 184.521 +640 32 2 2 +338.935 338.935 338.935 +342.898 342.898 342.898 +369.65 369.65 369.65 +640 64 2 2 +342.421 342.421 342.421 +341.617 341.617 341.617 +369.497 369.497 369.497 -- 2.39.5