# include <deal.II/base/index_set.h>
# include <deal.II/lac/full_matrix.h>
# include <deal.II/lac/exceptions.h>
-# include <deal.II/lac/trilinos_vector_base.h>
-# include <deal.II/lac/parallel_vector.h>
+# include <deal.II/lac/trilinos_vector.h>
+# include <deal.II/lac/vector_view.h>
# include <vector>
# include <cmath>
namespace TrilinosWrappers
{
// forward declarations
- class VectorBase;
class SparseMatrix;
class SparsityPattern;
*
* Source and destination must not be the same vector.
*
+ * This function can be called with several different vector objects,
+ * namely TrilinosWrappers::Vector, TrilinosWrappers::MPI::Vector as well
+ * as deal.II's own vector classes Vector<double> and
+ * parallel::distributed::Vector<double>.
+ *
* Note that both vectors have to be distributed vectors generated using
* the same Map as was used for the matrix in case you work on a
* distributed memory architecture, using the interface in the
* running on one processor, since the matrix object is inherently
* distributed. Otherwise, and exception will be thrown.
*/
- void vmult (VectorBase &dst,
- const VectorBase &src) const;
-
- /**
- * Same as before, but working with deal.II's own distributed vector
- * class.
- */
- void vmult (parallel::distributed::Vector<TrilinosScalar> &dst,
- const parallel::distributed::Vector<TrilinosScalar> &src) const;
-
- /**
- * Same as before, but working with deal.II's own vector class.
- */
- void vmult (dealii::Vector<TrilinosScalar> &dst,
- const dealii::Vector<TrilinosScalar> &src) const;
+ template<typename VectorType>
+ void vmult (VectorType &dst,
+ const VectorType &src) const;
/**
* Matrix-vector multiplication: let <i>dst = M<sup>T</sup>*src</i> with
*
* Source and destination must not be the same vector.
*
+ * This function can be called with several different vector objects,
+ * namely TrilinosWrappers::Vector, TrilinosWrappers::MPI::Vector as well
+ * as deal.II's own vector classes Vector<double> and
+ * parallel::distributed::Vector<double>.
+ *
* Note that both vectors have to be distributed vectors generated using
* the same Map as was used for the matrix in case you work on a
* distributed memory architecture, using the interface in the
* running on one processor, since the matrix object is inherently
* distributed. Otherwise, and exception will be thrown.
*/
- void Tvmult (VectorBase &dst,
- const VectorBase &src) const;
-
- /**
- * Same as before, but working with deal.II's own distributed vector
- * class.
- */
- void Tvmult (parallel::distributed::Vector<TrilinosScalar> &dst,
- const parallel::distributed::Vector<TrilinosScalar> &src) const;
-
- /**
- * Same as before, but working with deal.II's own vector class.
- */
- void Tvmult (dealii::Vector<TrilinosScalar> &dst,
- const dealii::Vector<TrilinosScalar> &src) const;
+ template <typename VectorType>
+ void Tvmult (VectorType &dst,
+ const VectorType &src) const;
/**
* Adding matrix-vector multiplication. Add <i>M*src</i> on <i>dst</i>
- inline
- void
- SparseMatrix::vmult (VectorBase &dst,
- const VectorBase &src) const
+ namespace internal
{
- Assert (&src != &dst, ExcSourceEqualsDestination());
- Assert (matrix->Filled(), ExcMatrixNotCompressed());
-
- Assert (src.vector_partitioner().SameAs(matrix->DomainMap()) == true,
- ExcMessage ("Column map of matrix does not fit with vector map!"));
- Assert (dst.vector_partitioner().SameAs(matrix->RangeMap()) == true,
- ExcMessage ("Row map of matrix does not fit with vector map!"));
-
- const int ierr = matrix->Multiply (false, src.trilinos_vector(),
- dst.trilinos_vector());
- Assert (ierr == 0, ExcTrilinosError(ierr));
- (void)ierr; // removes -Wunused-variable in optimized mode
- }
-
-
-
- inline
- void
- SparseMatrix::vmult (parallel::distributed::Vector<TrilinosScalar> &dst,
- const parallel::distributed::Vector<TrilinosScalar> &src) const
- {
- Assert (&src != &dst, ExcSourceEqualsDestination());
- Assert (matrix->Filled(), ExcMatrixNotCompressed());
-
- AssertDimension (dst.local_size(), static_cast<unsigned int>(matrix->RangeMap().NumMyElements()));
- AssertDimension (src.local_size(), static_cast<unsigned int>(matrix->DomainMap().NumMyElements()));
-
- Epetra_Vector tril_dst (View, matrix->RangeMap(), dst.begin());
- Epetra_Vector tril_src (View, matrix->DomainMap(),
- const_cast<double *>(src.begin()));
+ namespace SparseMatrix
+ {
+ template <typename VectorType>
+ inline
+ void check_vector_map_equality(const Epetra_CrsMatrix &,
+ const VectorType &,
+ const VectorType &)
+ {
+ }
- const int ierr = matrix->Multiply (false, tril_src, tril_dst);
- Assert (ierr == 0, ExcTrilinosError(ierr));
- (void)ierr; // removes -Wunused-variable in optimized mode
+ inline
+ void check_vector_map_equality(const Epetra_CrsMatrix &m,
+ const TrilinosWrappers::MPI::Vector &in,
+ const TrilinosWrappers::MPI::Vector &out)
+ {
+ Assert (in.vector_partitioner().SameAs(m.DomainMap()) == true,
+ ExcMessage ("Column map of matrix does not fit with vector map!"));
+ Assert (out.vector_partitioner().SameAs(m.RangeMap()) == true,
+ ExcMessage ("Row map of matrix does not fit with vector map!"));
+ }
+ }
}
-
+ template <typename VectorType>
inline
void
- SparseMatrix::vmult (dealii::Vector<TrilinosScalar> &dst,
- const dealii::Vector<TrilinosScalar> &src) const
+ SparseMatrix::vmult (VectorType &dst,
+ const VectorType &src) const
{
Assert (&src != &dst, ExcSourceEqualsDestination());
Assert (matrix->Filled(), ExcMatrixNotCompressed());
- AssertDimension (static_cast<unsigned int>(matrix->DomainMap().NumMyElements()),
- static_cast<unsigned int>(matrix->DomainMap().NumGlobalElements()));
- AssertDimension (dst.size(), static_cast<unsigned int>(matrix->RangeMap().NumMyElements()));
- AssertDimension (src.size(), static_cast<unsigned int>(matrix->DomainMap().NumMyElements()));
+ internal::SparseMatrix::check_vector_map_equality(*matrix, src, dst);
+ const int dst_local_size = dst.end() - dst.begin();
+ AssertDimension (dst_local_size, matrix->RangeMap().NumMyElements());
+ const int src_local_size = src.end() - src.begin();
+ AssertDimension (src_local_size, matrix->DomainMap().NumMyElements());
Epetra_Vector tril_dst (View, matrix->RangeMap(), dst.begin());
Epetra_Vector tril_src (View, matrix->DomainMap(),
- const_cast<double *>(src.begin()));
+ const_cast<TrilinosScalar *>(src.begin()));
const int ierr = matrix->Multiply (false, tril_src, tril_dst);
Assert (ierr == 0, ExcTrilinosError(ierr));
}
-
- inline
- void
- SparseMatrix::Tvmult (VectorBase &dst,
- const VectorBase &src) const
- {
- Assert (&src != &dst, ExcSourceEqualsDestination());
- Assert (matrix->Filled(), ExcMatrixNotCompressed());
-
- Assert (src.vector_partitioner().SameAs(matrix->RangeMap()) == true,
- ExcMessage ("Column map of matrix does not fit with vector map!"));
- Assert (dst.vector_partitioner().SameAs(matrix->DomainMap()) == true,
- ExcMessage ("Row map of matrix does not fit with vector map!"));
-
- const int ierr = matrix->Multiply (true, src.trilinos_vector(),
- dst.trilinos_vector());
- Assert (ierr == 0, ExcTrilinosError(ierr));
- (void)ierr; // removes -Wunused-variable in optimized mode
- }
-
-
-
- inline
- void
- SparseMatrix::Tvmult (parallel::distributed::Vector<TrilinosScalar> &dst,
- const parallel::distributed::Vector<TrilinosScalar> &src) const
- {
- Assert (&src != &dst, ExcSourceEqualsDestination());
- Assert (matrix->Filled(), ExcMatrixNotCompressed());
-
- AssertDimension (dst.local_size(), static_cast<unsigned int>(matrix->DomainMap().NumMyElements()));
- AssertDimension (src.local_size(), static_cast<unsigned int>(matrix->RangeMap().NumMyElements()));
-
- Epetra_Vector tril_dst (View, matrix->DomainMap(), dst.begin());
- Epetra_Vector tril_src (View, matrix->RangeMap(),
- const_cast<double *>(src.begin()));
-
- const int ierr = matrix->Multiply (true, tril_src, tril_dst);
- Assert (ierr == 0, ExcTrilinosError(ierr));
- (void)ierr; // removes -Wunused-variable in optimized mode
- }
-
-
-
+
+ template <typename VectorType>
inline
void
- SparseMatrix::Tvmult (dealii::Vector<TrilinosScalar> &dst,
- const dealii::Vector<TrilinosScalar> &src) const
+ SparseMatrix::Tvmult (VectorType &dst,
+ const VectorType &src) const
{
Assert (&src != &dst, ExcSourceEqualsDestination());
Assert (matrix->Filled(), ExcMatrixNotCompressed());
- AssertDimension (static_cast<unsigned int>(matrix->DomainMap().NumMyElements()),
- static_cast<unsigned int>(matrix->DomainMap().NumGlobalElements()));
- AssertDimension (dst.size(), static_cast<unsigned int>(matrix->DomainMap().NumMyElements()));
- AssertDimension (src.size(), static_cast<unsigned int>(matrix->RangeMap().NumMyElements()));
+ internal::SparseMatrix::check_vector_map_equality(*matrix, dst, src);
+ const int dst_local_size = dst.end() - dst.begin();
+ AssertDimension (dst_local_size, matrix->DomainMap().NumMyElements());
+ const int src_local_size = src.end() - src.begin();
+ AssertDimension (src_local_size, matrix->RangeMap().NumMyElements());
Epetra_Vector tril_dst (View, matrix->DomainMap(), dst.begin());
Epetra_Vector tril_src (View, matrix->RangeMap(),
{
Assert (&src != &dst, ExcSourceEqualsDestination());
- // Choose to reinit the vector with fast argument set, which does not
- // overwrite the content -- this is what we needs since we're going to
- // overwrite that anyway in the vmult operation.
- VectorType temp_vector;
- temp_vector.reinit(dst, true);
-
- vmult (temp_vector, src);
- dst += temp_vector;
+ // Reinit a temporary vector with fast argument set, which does not
+ // overwrite the content (to save time). However, the
+ // TrilinosWrappers::Vector classes do not support this, so create a
+ // deal.II local vector that has this fast setting. It will be accepted in
+ // vmult because it only checks the local size.
+ dealii::Vector<TrilinosScalar> temp_vector;
+ temp_vector.reinit(dst.end()-dst.begin(), true);
+ dealii::VectorView<TrilinosScalar> src_view(src.end()-src.begin(), src.begin());
+ dealii::VectorView<TrilinosScalar> dst_view(dst.end()-dst.begin(), dst.begin());
+ vmult (temp_vector, static_cast<const dealii::Vector<TrilinosScalar>&>(src_view));
+ if (dst_view.size() > 0)
+ dst_view += temp_vector;
}
const VectorType &src) const
{
Assert (&src != &dst, ExcSourceEqualsDestination());
- VectorType temp_vector;
- temp_vector.reinit(dst, true);
- Tvmult (temp_vector, src);
- dst += temp_vector;
+ // Reinit a temporary vector with fast argument set, which does not
+ // overwrite the content (to save time). However, the
+ // TrilinosWrappers::Vector classes do not support this, so create a
+ // deal.II local vector that has this fast setting. It will be accepted in
+ // vmult because it only checks the local size.
+ dealii::Vector<TrilinosScalar> temp_vector;
+ temp_vector.reinit(dst.end()-dst.begin(), true);
+ dealii::VectorView<TrilinosScalar> src_view(src.end()-src.begin(), src.begin());
+ dealii::VectorView<TrilinosScalar> dst_view(dst.end()-dst.begin(), dst.begin());
+ Tvmult (temp_vector, static_cast<const dealii::Vector<TrilinosScalar>&>(src_view));
+ if (dst_view.size() > 0)
+ dst_view += temp_vector;
}
--- /dev/null
+//---------------------- trilinos_matvec_03.cc ---------------------------
+// $Id$
+// Version: $Name$
+//
+// Copyright (C) 2011, 2013 by the deal.II authors
+//
+// This file is subject to QPL and may not be distributed
+// without copyright and license information. Please refer
+// to the file deal.II/doc/license.html for the text and
+// further information on this license.
+//
+//---------------------- trilinos_matvec_03.cc ---------------------------
+
+
+// Test whether TrilinosWrappers::SparseMatrix::(T)vmult(_add) gives same
+// result with Trilinos vector and parallel distributed vector when one
+// processor does not have any rows
+
+#include "../tests.h"
+#include <deal.II/base/utilities.h>
+#include <deal.II/base/index_set.h>
+#include <deal.II/lac/trilinos_sparse_matrix.h>
+#include <deal.II/lac/trilinos_sparsity_pattern.h>
+#include <deal.II/lac/trilinos_vector.h>
+#include <deal.II/lac/parallel_vector.h>
+#include <fstream>
+#include <iostream>
+#include <vector>
+
+
+void test ()
+{
+ const unsigned int n_procs = Utilities::MPI::n_mpi_processes(MPI_COMM_WORLD);
+ const unsigned int my_id = Utilities::MPI::this_mpi_process(MPI_COMM_WORLD);
+
+ const unsigned int n_rows = 3;
+ const unsigned int n_cols = 4;
+
+ IndexSet row_partitioning (n_rows);
+ IndexSet col_partitioning (n_cols);
+
+ if (n_procs == 1)
+ {
+ row_partitioning.add_range(0, n_rows);
+ col_partitioning.add_range(0, n_cols);
+ }
+ else if (n_procs == 2)
+ {
+ // row_partitioning should be { [0, 2), [2, n_rows) }
+ // col_partitioning should be { [0, 2), [2, n_cols) }
+ // col_relevant_set should be { [0, 3), [1, n_cols) }
+ if (my_id == 0)
+ {
+ row_partitioning.add_range(0, n_rows);
+ col_partitioning.add_range(0, n_cols);
+ }
+ }
+ else
+ Assert (false, ExcNotImplemented());
+
+ TrilinosWrappers::SparsityPattern sp (row_partitioning,
+ col_partitioning, MPI_COMM_WORLD);
+ if (my_id == 0)
+ {
+ sp.add (0, 0);
+ sp.add (0, 2);
+ sp.add (2, 3);
+ }
+ sp.compress();
+
+ TrilinosWrappers::SparseMatrix A;
+ A.clear ();
+ A.reinit (sp);
+ if (my_id == 0)
+ {
+ A.add (0, 0, 1);
+ A.add (0, 2, 1);
+ A.add (2, 3, 2.0);
+ }
+ A.compress(VectorOperation::add);
+
+ TrilinosWrappers::MPI::Vector x, y;
+ x.reinit (col_partitioning, MPI_COMM_WORLD);
+ y.reinit (row_partitioning, MPI_COMM_WORLD);
+
+ parallel::distributed::Vector<double>
+ dx (col_partitioning, col_partitioning, MPI_COMM_WORLD),
+ dy (row_partitioning, row_partitioning, MPI_COMM_WORLD);
+
+ for (unsigned int i=0; i<col_partitioning.n_elements(); ++i)
+ {
+ const unsigned int global_index = col_partitioning.nth_index_in_set(i);
+ dx(global_index) = (double)rand()/RAND_MAX;
+ x(global_index) = dx(global_index);
+ }
+ dy = 1.;
+
+ A.vmult (y, x);
+ A.vmult (dy, dx);
+
+ // compare whether we got the same result
+ // (should be no roundoff difference)
+ for (unsigned int i=0; i<row_partitioning.n_elements(); ++i)
+ {
+ const unsigned int global_index = row_partitioning.nth_index_in_set(i);
+ Assert (dy(global_index) == y(global_index), ExcInternalError());
+ }
+
+ A.vmult_add (y, x);
+ A.vmult_add (dy, dx);
+
+ // compare whether we got the same result
+ // (should be no roundoff difference)
+ for (unsigned int i=0; i<row_partitioning.n_elements(); ++i)
+ {
+ const unsigned int global_index = row_partitioning.nth_index_in_set(i);
+ Assert (dy(global_index) == y(global_index), ExcInternalError());
+ }
+
+ A.Tvmult (x, y);
+ A.Tvmult (dx, dy);
+ for (unsigned int i=0; i<col_partitioning.n_elements(); ++i)
+ {
+ const unsigned int global_index = col_partitioning.nth_index_in_set(i);
+ Assert (dx(global_index) == x(global_index), ExcInternalError());
+ }
+
+ A.Tvmult_add (x, y);
+ A.Tvmult_add (dx, dy);
+ for (unsigned int i=0; i<col_partitioning.n_elements(); ++i)
+ {
+ const unsigned int global_index = col_partitioning.nth_index_in_set(i);
+ Assert (dx(global_index) == x(global_index), ExcInternalError());
+ }
+
+ if (my_id == 0) deallog << "OK" << std::endl;
+}
+
+
+
+int main (int argc, char **argv)
+{
+ Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv);
+
+ const unsigned int n_procs = Utilities::MPI::n_mpi_processes(MPI_COMM_WORLD);
+ unsigned int myid = Utilities::MPI::this_mpi_process (MPI_COMM_WORLD);
+ deallog.push(Utilities::int_to_string(myid));
+
+ if (myid == 0)
+ {
+ std::ofstream logfile(output_file_for_mpi("trilinos_matvec_03").c_str());
+ deallog.attach(logfile);
+ deallog << std::setprecision(4);
+ deallog.depth_console(0);
+ deallog.threshold_double(1.e-10);
+
+ test();
+ }
+ else
+ test();
+
+}
projects / dealii-svn.git / commitdiff