From 8a26b20dbd1b3ad0d10234127bfb3c0d6e535b26 Mon Sep 17 00:00:00 2001 From: Jie Cheng Date: Sun, 14 Jan 2018 01:08:47 -0500 Subject: [PATCH] Follow-up of #5655: add PETScWrappers::SparseMatrix::mmult and PETScWrappers::MPI::SparseMatrix::mmult --- include/deal.II/lac/petsc_matrix_base.h | 4 +- .../lac/petsc_parallel_sparse_matrix.h | 20 ++++ include/deal.II/lac/petsc_sparse_matrix.h | 22 +++- source/lac/petsc_matrix_base.cc | 4 + source/lac/petsc_parallel_sparse_matrix.cc | 22 ++++ source/lac/petsc_sparse_matrix.cc | 22 ++++ tests/petsc/sparse_matrix_matrix_01.cc | 99 ++++++++++++++++ tests/petsc/sparse_matrix_matrix_01.output | 2 + tests/petsc/sparse_matrix_matrix_02.cc | 99 ++++++++++++++++ tests/petsc/sparse_matrix_matrix_02.output | 2 + tests/petsc/sparse_matrix_matrix_03.cc | 107 ++++++++++++++++++ tests/petsc/sparse_matrix_matrix_03.output | 2 + tests/petsc/sparse_matrix_matrix_04.cc | 107 ++++++++++++++++++ tests/petsc/sparse_matrix_matrix_04.output | 2 + 14 files changed, 511 insertions(+), 3 deletions(-) create mode 100644 tests/petsc/sparse_matrix_matrix_01.cc create mode 100644 tests/petsc/sparse_matrix_matrix_01.output create mode 100644 tests/petsc/sparse_matrix_matrix_02.cc create mode 100644 tests/petsc/sparse_matrix_matrix_02.output create mode 100644 tests/petsc/sparse_matrix_matrix_03.cc create mode 100644 tests/petsc/sparse_matrix_matrix_03.output create mode 100644 tests/petsc/sparse_matrix_matrix_04.cc create mode 100644 tests/petsc/sparse_matrix_matrix_04.output diff --git a/include/deal.II/lac/petsc_matrix_base.h b/include/deal.II/lac/petsc_matrix_base.h index 3e92e70dff..2be755ce3c 100644 --- a/include/deal.II/lac/petsc_matrix_base.h +++ b/include/deal.II/lac/petsc_matrix_base.h @@ -975,7 +975,7 @@ namespace PETScWrappers */ void mmult (MatrixBase &C, const MatrixBase &B, - const VectorBase &V = VectorBase()) const; + const VectorBase &V) const; /** * Base function to perform the matrix-matrix multiplication with @@ -995,7 +995,7 @@ namespace PETScWrappers */ void Tmmult (MatrixBase &C, const MatrixBase &B, - const VectorBase &V = VectorBase()) const; + const VectorBase &V) const; private: diff --git a/include/deal.II/lac/petsc_parallel_sparse_matrix.h b/include/deal.II/lac/petsc_parallel_sparse_matrix.h index 04c0849228..13f99e5043 100644 --- a/include/deal.II/lac/petsc_parallel_sparse_matrix.h +++ b/include/deal.II/lac/petsc_parallel_sparse_matrix.h @@ -408,6 +408,26 @@ namespace PETScWrappers */ IndexSet locally_owned_range_indices() const; + /** + * Perform the matrix-matrix multiplication $C = AB$, or, + * $C = A \text{diag}(V) B$ given a compatible vector $V$. + * + * This function calls MatrixBase::mmult() to do the actual work. + */ + void mmult (SparseMatrix &C, + const SparseMatrix &B, + const MPI::Vector &V = MPI::Vector()) const; + + /** + * Perform the matrix-matrix multiplication with the transpose of + * this, i.e., $C = A^T B$, or, + * $C = A^T \text{diag}(V) B$ given a compatible vector $V$. + * + * This function calls MatrixBase::Tmmult() to do the actual work. + */ + void Tmmult (SparseMatrix &C, + const SparseMatrix &B, + const MPI::Vector &V = MPI::Vector()) const; private: /** diff --git a/include/deal.II/lac/petsc_sparse_matrix.h b/include/deal.II/lac/petsc_sparse_matrix.h index 27c5cc715e..9f34cc8f16 100644 --- a/include/deal.II/lac/petsc_sparse_matrix.h +++ b/include/deal.II/lac/petsc_sparse_matrix.h @@ -23,7 +23,7 @@ # include # include -# include +# include # include DEAL_II_NAMESPACE_OPEN @@ -209,6 +209,26 @@ namespace PETScWrappers */ size_t n() const; + /** + * Perform the matrix-matrix multiplication $C = AB$, or, + * $C = A \text{diag}(V) B$ given a compatible vector $V$. + * + * This function calls MatrixBase::mmult() to do the actual work. + */ + void mmult (SparseMatrix &C, + const SparseMatrix &B, + const MPI::Vector &V = MPI::Vector()) const; + + /** + * Perform the matrix-matrix multiplication with the transpose of + * this, i.e., $C = A^T B$, or, + * $C = A^T \text{diag}(V) B$ given a compatible vector $V$. + * + * This function calls MatrixBase::Tmmult() to do the actual work. + */ + void Tmmult (SparseMatrix &C, + const SparseMatrix &B, + const MPI::Vector &V = MPI::Vector()) const; private: /** diff --git a/source/lac/petsc_matrix_base.cc b/source/lac/petsc_matrix_base.cc index 11f5b54f41..65bc5e9107 100644 --- a/source/lac/petsc_matrix_base.cc +++ b/source/lac/petsc_matrix_base.cc @@ -541,7 +541,11 @@ namespace PETScWrappers AssertThrow (ierr == 0, ExcPETScError(ierr)); if (transpose_left) { +#if DEAL_II_PETSC_VERSION_LT(3,8,0) ierr = MatTranspose(tmp, MAT_REUSE_MATRIX, &tmp); +#else + ierr = MatTranspose(tmp, MAT_INPLACE_MATRIX, &tmp); +#endif AssertThrow (ierr == 0, ExcPETScError(ierr)); } ierr = MatDiagonalScale (tmp, nullptr, V); diff --git a/source/lac/petsc_parallel_sparse_matrix.cc b/source/lac/petsc_parallel_sparse_matrix.cc index 65461a0e6b..9cee6d4b1e 100644 --- a/source/lac/petsc_parallel_sparse_matrix.cc +++ b/source/lac/petsc_parallel_sparse_matrix.cc @@ -707,6 +707,28 @@ namespace PETScWrappers return indices; } + void + SparseMatrix::mmult (SparseMatrix &C, + const SparseMatrix &B, + const MPI::Vector &V) const + { + // Simply forward to the protected member function of the base class + // that takes abstract matrix and vector arguments (to which the compiler + // automatically casts the arguments). + MatrixBase::mmult (C, B, V); + } + + void + SparseMatrix::Tmmult (SparseMatrix &C, + const SparseMatrix &B, + const MPI::Vector &V) const + { + // Simply forward to the protected member function of the base class + // that takes abstract matrix and vector arguments (to which the compiler + // automatically casts the arguments). + MatrixBase::Tmmult (C, B, V); + } + } } diff --git a/source/lac/petsc_sparse_matrix.cc b/source/lac/petsc_sparse_matrix.cc index f1004246cc..b27ee6b691 100644 --- a/source/lac/petsc_sparse_matrix.cc +++ b/source/lac/petsc_sparse_matrix.cc @@ -262,6 +262,28 @@ namespace PETScWrappers return n; } + void + SparseMatrix::mmult (SparseMatrix &C, + const SparseMatrix &B, + const MPI::Vector &V) const + { + // Simply forward to the protected member function of the base class + // that takes abstract matrix and vector arguments (to which the compiler + // automatically casts the arguments). + MatrixBase::mmult (C, B, V); + } + + void + SparseMatrix::Tmmult (SparseMatrix &C, + const SparseMatrix &B, + const MPI::Vector &V) const + { + // Simply forward to the protected member function of the base class + // that takes abstract matrix and vector arguments (to which the compiler + // automatically casts the arguments). + MatrixBase::Tmmult (C, B, V); + } + // Explicit instantiations // template diff --git a/tests/petsc/sparse_matrix_matrix_01.cc b/tests/petsc/sparse_matrix_matrix_01.cc new file mode 100644 index 0000000000..018319ae29 --- /dev/null +++ b/tests/petsc/sparse_matrix_matrix_01.cc @@ -0,0 +1,99 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2004 - 2015 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. +// +// --------------------------------------------------------------------- + + + +// check SparseMatrix::mmult + +#include "../tests.h" +#include +#include +#include +#include + + +void test () +{ + // A = [1, 0, 0; 4, 1, 0; 7, 2, 1] + PETScWrappers::SparseMatrix A(3, 3, 3); + A.set (0, 0, 1); + A.set (1, 0, 4); + A.set (1, 1, 1); + A.set (2, 0, 7); + A.set (2, 1, 2); + A.set (2, 2, 1); + A.compress (VectorOperation::insert); + + // B = [1, 2, 3; 0, -3, -6; 0, 0, 0] + PETScWrappers::SparseMatrix B(3, 3, 3); + B.set (0, 0, 1); + B.set (0, 1, 2); + B.set (0, 2, 3); + B.set (1, 1, -3); + B.set (1, 2, -6); + B.compress (VectorOperation::insert); + + PETScWrappers::SparseMatrix C(3, 3, 3); + + // C := AB = [1, 2, 3; 4, 5, 6; 7, 8, 9] + A.mmult(C, B); + + // make sure we get the expected result + for (unsigned int i = 0; i < C.m(); ++i) + for (unsigned int j = 0; j < C.n(); ++j) + AssertThrow (C(i, j) == 3 * i + j + 1, ExcInternalError()); + + deallog << "OK" << std::endl; +} + + + +int main (int argc, char **argv) +{ + initlog(); + + try + { + Utilities::MPI::MPI_InitFinalize mpi_initialization (argc, argv, 1); + { + test (); + } + + } + catch (std::exception &exc) + { + std::cerr << std::endl << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Exception on processing: " << std::endl + << exc.what() << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + + return 1; + } + catch (...) + { + std::cerr << std::endl << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Unknown exception!" << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + return 1; + }; +} diff --git a/tests/petsc/sparse_matrix_matrix_01.output b/tests/petsc/sparse_matrix_matrix_01.output new file mode 100644 index 0000000000..0fd8fc12f0 --- /dev/null +++ b/tests/petsc/sparse_matrix_matrix_01.output @@ -0,0 +1,2 @@ + +DEAL::OK diff --git a/tests/petsc/sparse_matrix_matrix_02.cc b/tests/petsc/sparse_matrix_matrix_02.cc new file mode 100644 index 0000000000..3e0e8602ab --- /dev/null +++ b/tests/petsc/sparse_matrix_matrix_02.cc @@ -0,0 +1,99 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2004 - 2015 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. +// +// --------------------------------------------------------------------- + + + +// check SparseMatrix::Tmmult + +#include "../tests.h" +#include +#include +#include +#include + + +void test () +{ + // A = [1, 4, 7; 0, 1, 2; 0, 0, 1] + PETScWrappers::SparseMatrix A(3, 3, 3); + A.set (0, 0, 1); + A.set (0, 1, 4); + A.set (0, 2, 7); + A.set (1, 1, 1); + A.set (1, 2, 2); + A.set (2, 2, 1); + A.compress (VectorOperation::insert); + + // B = [1, 2, 3; 0, -3, -6; 0, 0, 0] + PETScWrappers::SparseMatrix B(3, 3, 3); + B.set (0, 0, 1); + B.set (0, 1, 2); + B.set (0, 2, 3); + B.set (1, 1, -3); + B.set (1, 2, -6); + B.compress (VectorOperation::insert); + + PETScWrappers::SparseMatrix C(3, 3, 3); + + // C := A^T B = [1, 2, 3; 4, 5, 6; 7, 8, 9] + A.Tmmult(C, B); + + // make sure we get the expected result + for (unsigned int i = 0; i < C.m(); ++i) + for (unsigned int j = 0; j < C.n(); ++j) + AssertThrow (C(i, j) == 3 * i + j + 1, ExcInternalError()); + + deallog << "OK" << std::endl; +} + + + +int main (int argc, char **argv) +{ + initlog(); + + try + { + Utilities::MPI::MPI_InitFinalize mpi_initialization (argc, argv, 1); + { + test (); + } + + } + catch (std::exception &exc) + { + std::cerr << std::endl << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Exception on processing: " << std::endl + << exc.what() << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + + return 1; + } + catch (...) + { + std::cerr << std::endl << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Unknown exception!" << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + return 1; + }; +} diff --git a/tests/petsc/sparse_matrix_matrix_02.output b/tests/petsc/sparse_matrix_matrix_02.output new file mode 100644 index 0000000000..0fd8fc12f0 --- /dev/null +++ b/tests/petsc/sparse_matrix_matrix_02.output @@ -0,0 +1,2 @@ + +DEAL::OK diff --git a/tests/petsc/sparse_matrix_matrix_03.cc b/tests/petsc/sparse_matrix_matrix_03.cc new file mode 100644 index 0000000000..df68d8ec96 --- /dev/null +++ b/tests/petsc/sparse_matrix_matrix_03.cc @@ -0,0 +1,107 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2004 - 2015 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. +// +// --------------------------------------------------------------------- + + + +// check SparseMatrix::mmult + +#include "../tests.h" +#include +#include +#include +#include + + +void test () +{ + // A = [1, 0, 0; 4, 1, 0; 7, 2, 1] + PETScWrappers::SparseMatrix A(3, 3, 3); + A.set (0, 0, 1); + A.set (1, 0, 4); + A.set (1, 1, 1); + A.set (2, 0, 7); + A.set (2, 1, 2); + A.set (2, 2, 1); + A.compress (VectorOperation::insert); + + // B = [1, 2, 3; 0, -3, -6; 0, 0, 0] + PETScWrappers::SparseMatrix B(3, 3, 3); + B.set (0, 0, 1); + B.set (0, 1, 2); + B.set (0, 2, 3); + B.set (1, 1, -3); + B.set (1, 2, -6); + B.compress (VectorOperation::insert); + + // v = [2, 2, 2] + IndexSet indices(3); + indices.add_range(0, 3); + PETScWrappers::MPI::Vector v (indices, MPI_COMM_WORLD); + for (unsigned int i = 0; i < v.size(); ++i) + v(i) = 2; + v.compress (VectorOperation::insert); + + PETScWrappers::SparseMatrix C(3, 3, 3); + + // C := A diag(v) B = [2, 4, 6; 8, 10, 12; 14, 16, 18] + A.mmult(C, B, v); + + // make sure we get the expected result + for (unsigned int i = 0; i < C.m(); ++i) + for (unsigned int j = 0; j < C.n(); ++j) + AssertThrow (C(i, j) == 6 * i + 2 * j + 2, ExcInternalError()); + + deallog << "OK" << std::endl; +} + + + +int main (int argc, char **argv) +{ + initlog(); + + try + { + Utilities::MPI::MPI_InitFinalize mpi_initialization (argc, argv, 1); + { + test (); + } + + } + catch (std::exception &exc) + { + std::cerr << std::endl << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Exception on processing: " << std::endl + << exc.what() << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + + return 1; + } + catch (...) + { + std::cerr << std::endl << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Unknown exception!" << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + return 1; + }; +} diff --git a/tests/petsc/sparse_matrix_matrix_03.output b/tests/petsc/sparse_matrix_matrix_03.output new file mode 100644 index 0000000000..0fd8fc12f0 --- /dev/null +++ b/tests/petsc/sparse_matrix_matrix_03.output @@ -0,0 +1,2 @@ + +DEAL::OK diff --git a/tests/petsc/sparse_matrix_matrix_04.cc b/tests/petsc/sparse_matrix_matrix_04.cc new file mode 100644 index 0000000000..09e834766f --- /dev/null +++ b/tests/petsc/sparse_matrix_matrix_04.cc @@ -0,0 +1,107 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2004 - 2015 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. +// +// --------------------------------------------------------------------- + + + +// check SparseMatrix::Tmmult + +#include "../tests.h" +#include +#include +#include +#include + + +void test () +{ + // A = [1, 4, 7; 0, 1, 2; 0, 0, 1] + PETScWrappers::SparseMatrix A(3, 3, 3); + A.set (0, 0, 1); + A.set (0, 1, 4); + A.set (0, 2, 7); + A.set (1, 1, 1); + A.set (1, 2, 2); + A.set (2, 2, 1); + A.compress (VectorOperation::insert); + + // B = [1, 2, 3; 0, -3, -6; 0, 0, 0] + PETScWrappers::SparseMatrix B(3, 3, 3); + B.set (0, 0, 1); + B.set (0, 1, 2); + B.set (0, 2, 3); + B.set (1, 1, -3); + B.set (1, 2, -6); + B.compress (VectorOperation::insert); + + // v = [2, 2, 2] + IndexSet indices(3); + indices.add_range(0, 3); + PETScWrappers::MPI::Vector v (indices, MPI_COMM_WORLD); + for (unsigned int i = 0; i < v.size(); ++i) + v(i) = 2; + v.compress (VectorOperation::insert); + + PETScWrappers::SparseMatrix C(3, 3, 3); + + // C := A^T diag(v) B = [2, 4, 6; 8, 10, 12; 14, 16, 18] + A.Tmmult(C, B, v); + + // make sure we get the expected result + for (unsigned int i = 0; i < C.m(); ++i) + for (unsigned int j = 0; j < C.n(); ++j) + AssertThrow (C(i, j) == 6 * i + 2 * j + 2, ExcInternalError()); + + deallog << "OK" << std::endl; +} + + + +int main (int argc, char **argv) +{ + initlog(); + + try + { + Utilities::MPI::MPI_InitFinalize mpi_initialization (argc, argv, 1); + { + test (); + } + + } + catch (std::exception &exc) + { + std::cerr << std::endl << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Exception on processing: " << std::endl + << exc.what() << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + + return 1; + } + catch (...) + { + std::cerr << std::endl << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Unknown exception!" << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + return 1; + }; +} diff --git a/tests/petsc/sparse_matrix_matrix_04.output b/tests/petsc/sparse_matrix_matrix_04.output new file mode 100644 index 0000000000..0fd8fc12f0 --- /dev/null +++ b/tests/petsc/sparse_matrix_matrix_04.output @@ -0,0 +1,2 @@ + +DEAL::OK -- 2.39.5