]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Follow-up of #5655: add PETScWrappers::SparseMatrix::mmult and PETScWrappers::MPI... 5732/head
authorJie Cheng <chengjiehust@gmail.com>
Sun, 14 Jan 2018 06:08:47 +0000 (01:08 -0500)
committerJie Cheng <chengjiehust@gmail.com>
Wed, 7 Mar 2018 00:46:25 +0000 (19:46 -0500)
14 files changed:
include/deal.II/lac/petsc_matrix_base.h
include/deal.II/lac/petsc_parallel_sparse_matrix.h
include/deal.II/lac/petsc_sparse_matrix.h
source/lac/petsc_matrix_base.cc
source/lac/petsc_parallel_sparse_matrix.cc
source/lac/petsc_sparse_matrix.cc
tests/petsc/sparse_matrix_matrix_01.cc [new file with mode: 0644]
tests/petsc/sparse_matrix_matrix_01.output [new file with mode: 0644]
tests/petsc/sparse_matrix_matrix_02.cc [new file with mode: 0644]
tests/petsc/sparse_matrix_matrix_02.output [new file with mode: 0644]
tests/petsc/sparse_matrix_matrix_03.cc [new file with mode: 0644]
tests/petsc/sparse_matrix_matrix_03.output [new file with mode: 0644]
tests/petsc/sparse_matrix_matrix_04.cc [new file with mode: 0644]
tests/petsc/sparse_matrix_matrix_04.output [new file with mode: 0644]

index 3e92e70dffc4c1a80b0b9a9a0c2c2cc2ae0ab3ff..2be755ce3c3116f9225f19251d4054f538798ab6 100644 (file)
@@ -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:
 
index 04c0849228371515868e72ecd2b1187f0f18d36f..13f99e5043555207365878e937c3578b2b936985 100644 (file)
@@ -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
+       * <tt>this</tt>, 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:
 
       /**
index 27c5cc715e36c552f612a3de455b6d1802248dcf..9f34cc8f1668234e58c51c01db6b778ad73312b0 100644 (file)
@@ -23,7 +23,7 @@
 
 #  include <deal.II/lac/exceptions.h>
 #  include <deal.II/lac/petsc_matrix_base.h>
-#  include <deal.II/lac/petsc_vector_base.h>
+#  include <deal.II/lac/petsc_parallel_vector.h>
 #  include <vector>
 
 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
+      * <tt>this</tt>, 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:
 
     /**
index 11f5b54f4182d21a18cbdb82022544d821e31dd3..65bc5e910768fb8ae221d0db9ecf6bde2ba87f31 100644 (file)
@@ -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);
index 65461a0e6b47bfbc135a3f7bcb95aff2ac78264d..9cee6d4b1ecaceb06320d2c7390e6577b107787e 100644 (file)
@@ -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);
+    }
+
   }
 }
 
index f1004246ccd5100672d70c4108afd1a155d99f60..b27ee6b6912527c7f7d642200220f0dcf6767a71 100644 (file)
@@ -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 (file)
index 0000000..018319a
--- /dev/null
@@ -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 <deal.II/lac/petsc_parallel_vector.h>
+#include <deal.II/lac/petsc_sparse_matrix.h>
+#include <iostream>
+#include <vector>
+
+
+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 (file)
index 0000000..0fd8fc1
--- /dev/null
@@ -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 (file)
index 0000000..3e0e860
--- /dev/null
@@ -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 <deal.II/lac/petsc_parallel_vector.h>
+#include <deal.II/lac/petsc_sparse_matrix.h>
+#include <iostream>
+#include <vector>
+
+
+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 (file)
index 0000000..0fd8fc1
--- /dev/null
@@ -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 (file)
index 0000000..df68d8e
--- /dev/null
@@ -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 <deal.II/lac/petsc_parallel_vector.h>
+#include <deal.II/lac/petsc_sparse_matrix.h>
+#include <iostream>
+#include <vector>
+
+
+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 (file)
index 0000000..0fd8fc1
--- /dev/null
@@ -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 (file)
index 0000000..09e8347
--- /dev/null
@@ -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 <deal.II/lac/petsc_parallel_vector.h>
+#include <deal.II/lac/petsc_sparse_matrix.h>
+#include <iostream>
+#include <vector>
+
+
+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 (file)
index 0000000..0fd8fc1
--- /dev/null
@@ -0,0 +1,2 @@
+
+DEAL::OK

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.