*/
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
*/
void Tmmult (MatrixBase &C,
const MatrixBase &B,
- const VectorBase &V = VectorBase()) const;
+ const VectorBase &V) const;
private:
*/
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:
/**
# 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
*/
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:
/**
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);
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);
+ }
+
}
}
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
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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;
+ };
+}
--- /dev/null
+
+DEAL::OK
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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;
+ };
+}
--- /dev/null
+
+DEAL::OK
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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;
+ };
+}
--- /dev/null
+
+DEAL::OK
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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;
+ };
+}
--- /dev/null
+
+DEAL::OK
projects / dealii.git / commitdiff