<ol>
+<li> New: The classes TrilinosWrappers::SparseMatrix and
+TrilinosWrappers::BlockSparseMatrix now fully implement vmult and Tvmult with
+deal.II's own vector classes Vector<double> and
+parallel::distributed::Vector<double>.
+<br>
+(Martin Kronbichler, 2013/05/08)
+</li>
+
<li> Improved: The matrix-vector product ChunkSparseMatrix::vmult now runs in
parallel in shared memory.
<br>
#ifdef DEAL_II_WITH_TRILINOS
# include <deal.II/base/table.h>
+# include <deal.II/base/template_constraints.h>
# include <deal.II/lac/block_matrix_base.h>
# include <deal.II/lac/trilinos_sparse_matrix.h>
# include <deal.II/lac/trilinos_block_vector.h>
unsigned int n_nonzero_elements () const;
/**
- * Matrix-vector multiplication:
- * let $dst = M*src$ with $M$
- * being this matrix.
+ * Matrix-vector multiplication: let $dst = M*src$ with $M$ being this
+ * matrix. The vector types can be block vectors or non-block vectors
+ * (only if the matrix has only one row or column, respectively), and need
+ * to define TrilinosWrappers::SparseMatrix::vmult.
*/
- void vmult (MPI::BlockVector &dst,
- const MPI::BlockVector &src) const;
-
-
- /**
- * Matrix-vector multiplication:
- * let $dst = M*src$ with $M$
- * being this matrix, now applied
- * to localized block vectors
- * (works only when run on one
- * processor).
- */
- void vmult (BlockVector &dst,
- const BlockVector &src) const;
-
- /**
- * Matrix-vector
- * multiplication. Just like the
- * previous function, but only
- * applicable if the matrix has
- * only one block column.
- */
- void vmult (MPI::BlockVector &dst,
- const MPI::Vector &src) const;
-
- /**
- * Matrix-vector
- * multiplication. Just like the
- * previous function, but only
- * applicable if the matrix has
- * only one block column, now
- * applied to localized vectors
- * (works only when run on one
- * processor).
- */
- void vmult (BlockVector &dst,
- const Vector &src) const;
-
- /**
- * Matrix-vector
- * multiplication. Just like the
- * previous function, but only
- * applicable if the matrix has
- * only one block row.
- */
- void vmult (MPI::Vector &dst,
- const MPI::BlockVector &src) const;
-
- /**
- * Matrix-vector
- * multiplication. Just like the
- * previous function, but only
- * applicable if the matrix has
- * only one block row, now
- * applied to localized vectors
- * (works only when run on one
- * processor).
- */
- void vmult (Vector &dst,
- const BlockVector &src) const;
-
- /**
- * Matrix-vector
- * multiplication. Just like the
- * previous function, but only
- * applicable if the matrix has
- * only one block.
- */
- void vmult (VectorBase &dst,
- const VectorBase &src) const;
-
- /**
- * Matrix-vector multiplication:
- * let $dst = M^T*src$ with $M$
- * being this matrix. This
- * function does the same as
- * vmult() but takes the
- * transposed matrix.
- */
- void Tvmult (MPI::BlockVector &dst,
- const MPI::BlockVector &src) const;
-
- /**
- * Matrix-vector multiplication:
- * let $dst = M^T*src$ with $M$
- * being this matrix. This
- * function does the same as
- * vmult() but takes the
- * transposed matrix, now applied
- * to localized Trilinos vectors
- * (works only when run on one
- * processor).
- */
- void Tvmult (BlockVector &dst,
- const BlockVector &src) const;
-
- /**
- * Matrix-vector
- * multiplication. Just like the
- * previous function, but only
- * applicable if the matrix has
- * only one block row.
- */
- void Tvmult (MPI::BlockVector &dst,
- const MPI::Vector &src) const;
+ template <typename VectorType1, typename VectorType2>
+ void vmult (VectorType1 &dst,
+ const VectorType2 &src) const;
/**
- * Matrix-vector
- * multiplication. Just like the
- * previous function, but only
- * applicable if the matrix has
- * only one block row, now
- * applied to localized Trilinos
- * vectors (works only when run
- * on one processor).
- */
- void Tvmult (BlockVector &dst,
- const Vector &src) const;
-
- /**
- * Matrix-vector
- * multiplication. Just like the
- * previous function, but only
- * applicable if the matrix has
- * only one block column.
- */
- void Tvmult (MPI::Vector &dst,
- const MPI::BlockVector &src) const;
-
- /**
- * Matrix-vector
- * multiplication. Just like the
- * previous function, but only
- * applicable if the matrix has
- * only one block column, now
- * applied to localized Trilinos
- * vectors (works only when run
- * on one processor).
- */
- void Tvmult (Vector &dst,
- const BlockVector &src) const;
-
- /**
- * Matrix-vector
- * multiplication. Just like the
- * previous function, but only
- * applicable if the matrix has
- * only one block.
+ * Matrix-vector multiplication: let $dst = M^T*src$ with $M$ being this
+ * matrix. This function does the same as vmult() but takes the transposed
+ * matrix.
*/
- void Tvmult (VectorBase &dst,
- const VectorBase &src) const;
+ template <typename VectorType1, typename VectorType2>
+ void Tvmult (VectorType1 &dst,
+ const VectorType2 &src) const;
/**
* Compute the residual of an
<< "The blocks [" << arg1 << ',' << arg2 << "] and ["
<< arg3 << ',' << arg4 << "] have differing column numbers.");
///@}
+
+ private:
+ /**
+ * Internal version of (T)vmult with two block vectors
+ */
+ template <typename VectorType1, typename VectorType2>
+ void vmult (VectorType1 &dst,
+ const VectorType2 &src,
+ const bool transpose,
+ const dealii::internal::bool2type<true>,
+ const dealii::internal::bool2type<true>) const;
+
+ /**
+ * Internal version of (T)vmult where the source vector is a block vector
+ * but the destination vector is a non-block vector
+ */
+ template <typename VectorType1, typename VectorType2>
+ void vmult (VectorType1 &dst,
+ const VectorType2 &src,
+ const bool transpose,
+ const dealii::internal::bool2type<false>,
+ const dealii::internal::bool2type<true>) const;
+
+ /**
+ * Internal version of (T)vmult where the source vector is a non-block
+ * vector but the destination vector is a block vector
+ */
+ template <typename VectorType1, typename VectorType2>
+ void vmult (VectorType1 &dst,
+ const VectorType2 &src,
+ const bool transpose,
+ const dealii::internal::bool2type<true>,
+ const dealii::internal::bool2type<false>) const;
+
+ /**
+ * Internal version of (T)vmult where both source vector and the
+ * destination vector are non-block vectors (only defined if the matrix
+ * consists of only one block)
+ */
+ template <typename VectorType1, typename VectorType2>
+ void vmult (VectorType1 &dst,
+ const VectorType2 &src,
+ const bool transpose,
+ const dealii::internal::bool2type<false>,
+ const dealii::internal::bool2type<false>) const;
};
+ template <typename VectorType1, typename VectorType2>
inline
void
- BlockSparseMatrix::vmult (MPI::BlockVector &dst,
- const MPI::BlockVector &src) const
+ BlockSparseMatrix::vmult (VectorType1 &dst,
+ const VectorType2 &src) const
{
- BaseClass::vmult_block_block (dst, src);
+ vmult(dst, src, false,
+ dealii::internal::bool2type<IsBlockVector<VectorType1>::value>(),
+ dealii::internal::bool2type<IsBlockVector<VectorType2>::value>());
}
+ template <typename VectorType1, typename VectorType2>
inline
void
- BlockSparseMatrix::vmult (BlockVector &dst,
- const BlockVector &src) const
+ BlockSparseMatrix::Tvmult (VectorType1 &dst,
+ const VectorType2 &src) const
{
- BaseClass::vmult_block_block (dst, src);
+ vmult(dst, src, true,
+ dealii::internal::bool2type<IsBlockVector<VectorType1>::value>(),
+ dealii::internal::bool2type<IsBlockVector<VectorType2>::value>());
}
+ template <typename VectorType1, typename VectorType2>
inline
void
- BlockSparseMatrix::vmult (MPI::BlockVector &dst,
- const MPI::Vector &src) const
+ BlockSparseMatrix::vmult (VectorType1 &dst,
+ const VectorType2 &src,
+ const bool transpose,
+ dealii::internal::bool2type<true>,
+ dealii::internal::bool2type<true>) const
{
- BaseClass::vmult_block_nonblock (dst, src);
+ if (transpose == true)
+ BaseClass::Tvmult_block_block (dst, src);
+ else
+ BaseClass::vmult_block_block (dst, src);
}
- inline
- void
- BlockSparseMatrix::vmult (BlockVector &dst,
- const Vector &src) const
- {
- BaseClass::vmult_block_nonblock (dst, src);
- }
-
-
-
- inline
- void
- BlockSparseMatrix::vmult (MPI::Vector &dst,
- const MPI::BlockVector &src) const
- {
- BaseClass::vmult_nonblock_block (dst, src);
- }
-
-
-
- inline
- void
- BlockSparseMatrix::vmult (Vector &dst,
- const BlockVector &src) const
- {
- BaseClass::vmult_nonblock_block (dst, src);
- }
-
-
-
- inline
- void
- BlockSparseMatrix::vmult (VectorBase &dst,
- const VectorBase &src) const
- {
- BaseClass::vmult_nonblock_nonblock (dst, src);
- }
-
-
-
- inline
- void
- BlockSparseMatrix::Tvmult (MPI::BlockVector &dst,
- const MPI::BlockVector &src) const
- {
- BaseClass::Tvmult_block_block (dst, src);
- }
-
-
-
- inline
- void
- BlockSparseMatrix::Tvmult (BlockVector &dst,
- const BlockVector &src) const
- {
- BaseClass::Tvmult_block_block (dst, src);
- }
-
-
-
- inline
- void
- BlockSparseMatrix::Tvmult (MPI::BlockVector &dst,
- const MPI::Vector &src) const
- {
- BaseClass::Tvmult_block_nonblock (dst, src);
- }
-
-
-
- inline
- void
- BlockSparseMatrix::Tvmult (BlockVector &dst,
- const Vector &src) const
- {
- BaseClass::Tvmult_block_nonblock (dst, src);
- }
-
-
+ template <typename VectorType1, typename VectorType2>
inline
void
- BlockSparseMatrix::Tvmult (MPI::Vector &dst,
- const MPI::BlockVector &src) const
+ BlockSparseMatrix::vmult (VectorType1 &dst,
+ const VectorType2 &src,
+ const bool transpose,
+ dealii::internal::bool2type<false>,
+ dealii::internal::bool2type<true>) const
{
- BaseClass::Tvmult_nonblock_block (dst, src);
+ if (transpose == true)
+ BaseClass::Tvmult_nonblock_block (dst, src);
+ else
+ BaseClass::vmult_nonblock_block (dst, src);
}
+ template <typename VectorType1, typename VectorType2>
inline
void
- BlockSparseMatrix::Tvmult (Vector &dst,
- const BlockVector &src) const
+ BlockSparseMatrix::vmult (VectorType1 &dst,
+ const VectorType2 &src,
+ const bool transpose,
+ dealii::internal::bool2type<true>,
+ dealii::internal::bool2type<false>) const
{
- BaseClass::Tvmult_nonblock_block (dst, src);
+ if (transpose == true)
+ BaseClass::Tvmult_block_nonblock (dst, src);
+ else
+ BaseClass::vmult_block_nonblock (dst, src);
}
+ template <typename VectorType1, typename VectorType2>
inline
void
- BlockSparseMatrix::Tvmult (VectorBase &dst,
- const VectorBase &src) const
+ BlockSparseMatrix::vmult (VectorType1 &dst,
+ const VectorType2 &src,
+ const bool transpose,
+ dealii::internal::bool2type<false>,
+ dealii::internal::bool2type<false>) const
{
- BaseClass::Tvmult_nonblock_nonblock (dst, src);
+ if (transpose == true)
+ BaseClass::Tvmult_nonblock_nonblock (dst, src);
+ else
+ BaseClass::vmult_nonblock_nonblock (dst, src);
}
}
//@{
/**
- * Matrix-vector multiplication:
- * let <i>dst = M*src</i> with
- * <i>M</i> being this matrix.
+ * Matrix-vector multiplication: let <i>dst = M*src</i> with <i>M</i>
+ * being this matrix.
*
- * Source and destination must
- * not be the same vector.
+ * Source and destination must not be the same vector.
*
- * 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
- * TrilinosWrappers::VectorBase
- * class (or one of the two
- * derived classes Vector and
- * MPI::Vector).
+ * 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
+ * TrilinosWrappers::VectorBase class (or one of the two derived classes
+ * Vector and MPI::Vector).
*
- * In case of a localized Vector,
- * this function will only work
- * when running on one processor,
- * since the matrix object is
- * inherently
- * distributed. Otherwise, and
- * exception will be thrown.
+ * In case of a localized Vector, this function will only work when
+ * 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
+ * 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.
+ * Same as before, but working with deal.II's own vector class.
*/
void vmult (dealii::Vector<TrilinosScalar> &dst,
const dealii::Vector<TrilinosScalar> &src) const;
/**
- * Matrix-vector multiplication:
- * let <i>dst =
- * M<sup>T</sup>*src</i> with
- * <i>M</i> being this
- * matrix. This function does the
- * same as vmult() but takes the
- * transposed matrix.
+ * Matrix-vector multiplication: let <i>dst = M<sup>T</sup>*src</i> with
+ * <i>M</i> being this matrix. This function does the same as vmult() but
+ * takes the transposed matrix.
*
- * Source and destination must
- * not be the same vector.
+ * Source and destination must not be the same vector.
*
- * 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
- * TrilinosWrappers::VectorBase
- * class (or one of the two
- * derived classes Vector and
- * MPI::Vector).
+ * 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
+ * TrilinosWrappers::VectorBase class (or one of the two derived classes
+ * Vector and MPI::Vector).
*
- * In case of a localized Vector,
- * this function will only work
- * when running on one processor,
- * since the matrix object is
- * inherently
- * distributed. Otherwise, and
- * exception will be thrown.
+ * In case of a localized Vector, this function will only work when
+ * 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
+ * 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.
+ * Same as before, but working with deal.II's own vector class.
*/
void Tvmult (dealii::Vector<TrilinosScalar> &dst,
const dealii::Vector<TrilinosScalar> &src) const;
/**
- * Adding Matrix-vector
- * multiplication. Add
- * <i>M*src</i> on <i>dst</i>
- * with <i>M</i> being this
- * matrix.
+ * Adding matrix-vector multiplication. Add <i>M*src</i> on <i>dst</i>
+ * with <i>M</i> being this matrix.
*
- * Source and destination must
- * not be the same vector.
+ * Source and destination must not be the same vector.
*
- * 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
- * TrilinosWrappers::VectorBase
- * class (or one of the two
- * derived classes Vector and
- * MPI::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>.
+ *
+ * When using a vector of type TrilinosWrappers::MPI::Vector, both vectors
+ * have to be distributed vectors generated using the same Map as was used
+ * for the matrix rows and columns in case you work on a distributed
+ * memory architecture, using the interface in the
+ * TrilinosWrappers::VectorBase class.
+ *
+ * In case of a localized Vector (i.e., TrilinosWrappers::Vector or
+ * Vector<double>), this function will only work when running on one
+ * processor, since the matrix object is inherently
+ * distributed. Otherwise, and exception will be thrown.
*
- * In case of a localized Vector,
- * this function will only work
- * when running on one processor,
- * since the matrix object is
- * inherently
- * distributed. Otherwise, and
- * exception will be thrown.
*/
- void vmult_add (VectorBase &dst,
- const VectorBase &src) const;
+ template<typename VectorType>
+ void vmult_add (VectorType &dst,
+ const VectorType &src) const;
/**
- * Adding Matrix-vector
- * multiplication. Add
- * <i>M<sup>T</sup>*src</i> to
- * <i>dst</i> with <i>M</i> being
- * this matrix. This function
- * does the same as vmult_add()
- * but takes the transposed
- * matrix.
+ * Adding matrix-vector multiplication. Add <i>M<sup>T</sup>*src</i> to
+ * <i>dst</i> with <i>M</i> being this matrix. This function does the same
+ * as vmult_add() but takes the transposed matrix.
*
- * Source and destination must
- * not be the same vector.
+ * Source and destination must not be the same vector.
*
- * 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
- * TrilinosWrappers::VectorBase
- * class (or one of the two
- * derived classes Vector and
- * MPI::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>.
*
- * In case of a localized Vector,
- * this function will only work
- * when running on one processor,
- * since the matrix object is
- * inherently
- * distributed. Otherwise, and
- * exception will be thrown.
+ * When using a vector of type TrilinosWrappers::MPI::Vector, both vectors
+ * have to be distributed vectors generated using the same Map as was used
+ * for the matrix rows and columns in case you work on a distributed
+ * memory architecture, using the interface in the
+ * TrilinosWrappers::VectorBase class.
+ *
+ * In case of a localized Vector (i.e., TrilinosWrappers::Vector or
+ * Vector<double>), this function will only work when running on one
+ * processor, since the matrix object is inherently
+ * distributed. Otherwise, and exception will be thrown.
*/
- void Tvmult_add (VectorBase &dst,
- const VectorBase &src) const;
+ template <typename VectorType>
+ void Tvmult_add (VectorType &dst,
+ const VectorType &src) const;
/**
* Return the square of the norm
*/
bool compressed;
- /**
- * An internal Trilinos vector that
- * is used for accelerating vmult_add
- * functions (in order not to need to
- * recreate temporary vectors every
- * time that function is called).
- */
- mutable VectorBase temp_vector;
-
- /**
- * An internal array of integer
- * values that is used to store the
- * column indices when
- * adding/inserting local data into
- * the (large) sparse matrix.
- */
- std::vector<unsigned int> column_indices;
-
- /**
- * An internal array of double values
- * that is used to store the column
- * indices when adding/inserting
- * local data into the (large) sparse
- * matrix.
- */
- std::vector<TrilinosScalar> column_values;
-
/**
* To allow calling protected
* prepare_add() and
last_action = Insert;
int *col_index_ptr;
- TrilinosScalar const *col_value_ptr;
+ TrilinosScalar *col_value_ptr;
int n_columns;
- // If we don't elide zeros, the pointers
- // are already available...
+ TrilinosScalar short_val_array[100];
+ int short_index_array[100];
+ std::vector<TrilinosScalar> long_val_array;
+ std::vector<int> long_index_array;
+
+
+ // If we don't elide zeros, the pointers are already available... need to
+ // cast to non-const pointers as that is the format taken by Trilinos (but
+ // we will not modify const data)
if (elide_zero_values == false)
{
col_index_ptr = (int *)col_indices;
- col_value_ptr = values;
+ col_value_ptr = const_cast<TrilinosScalar*>(values);
n_columns = n_cols;
}
else
{
- // Otherwise, extract nonzero values in
- // each row and get the respective
- // indices.
- if (column_indices.size() < n_cols)
+ // Otherwise, extract nonzero values in each row and get the
+ // respective indices.
+ if (n_cols > 100)
{
- column_indices.resize(n_cols);
- column_values.resize(n_cols);
+ long_val_array.resize(n_cols);
+ long_index_array.resize(n_cols);
+ col_index_ptr = &long_index_array[0];
+ col_value_ptr = &long_val_array[0];
+ }
+ else
+ {
+ col_index_ptr = &short_index_array[0];
+ col_value_ptr = &short_val_array[0];
}
n_columns = 0;
Assert (numbers::is_finite(value), ExcNumberNotFinite());
if (value != 0)
{
- column_indices[n_columns] = col_indices[j];
- column_values[n_columns] = value;
+ col_index_ptr[n_columns] = col_indices[j];
+ col_value_ptr[n_columns] = value;
n_columns++;
}
}
Assert(n_columns <= (int)n_cols, ExcInternalError());
-
- col_index_ptr = (int *)&column_indices[0];
- col_value_ptr = &column_values[0];
}
- // If the calling matrix owns the row to
- // which we want to insert values, we
- // can directly call the Epetra_CrsMatrix
- // input function, which is much faster
- // than the Epetra_FECrsMatrix
- // function. We distinguish between two
- // cases: the first one is when the matrix
- // is not filled (i.e., it is possible to
- // add new elements to the sparsity pattern),
- // and the second one is when the pattern is
- // already fixed. In the former case, we
- // add the possibility to insert new values,
- // and in the second we just replace
+ // If the calling matrix owns the row to which we want to insert values,
+ // we can directly call the Epetra_CrsMatrix input function, which is much
+ // faster than the Epetra_FECrsMatrix function. We distinguish between two
+ // cases: the first one is when the matrix is not filled (i.e., it is
+ // possible to add new elements to the sparsity pattern), and the second
+ // one is when the pattern is already fixed. In the former case, we add
+ // the possibility to insert new values, and in the second we just replace
// data.
if (row_partitioner().MyGID(static_cast<int>(row)) == true)
{
if (matrix->Filled() == false)
{
ierr = matrix->Epetra_CrsMatrix::InsertGlobalValues(row, n_columns,
- const_cast<double *>(col_value_ptr),
+ col_value_ptr,
col_index_ptr);
- // When inserting elements, we do
- // not want to create exceptions in
- // the case when inserting non-local
- // data (since that's what we want
- // to do right now).
+ // When inserting elements, we do not want to create exceptions in
+ // the case when inserting non-local data (since that's what we
+ // want to do right now).
if (ierr > 0)
ierr = 0;
}
else
ierr = matrix->Epetra_CrsMatrix::ReplaceGlobalValues(row, n_columns,
- const_cast<double *>(col_value_ptr),
+ col_value_ptr,
col_index_ptr);
}
else
{
- // When we're at off-processor data, we
- // have to stick with the standard
- // Insert/ReplaceGlobalValues
- // function. Nevertheless, the way we
- // call it is the fastest one (any other
- // will lead to repeated allocation and
- // deallocation of memory in order to
- // call the function we already use,
- // which is very unefficient if writing
- // one element at a time).
+ // When we're at off-processor data, we have to stick with the
+ // standard Insert/ReplaceGlobalValues function. Nevertheless, the way
+ // we call it is the fastest one (any other will lead to repeated
+ // allocation and deallocation of memory in order to call the function
+ // we already use, which is very unefficient if writing one element at
+ // a time).
compressed = false;
if (matrix->Filled() == false)
if (value == 0)
{
- // we have to do checkings on Insert/Add
- // in any case
- // to be consistent with the MPI
- // communication model (see the comments
- // in the documentation of
- // TrilinosWrappers::Vector), but we can
- // save some work if the addend is
- // zero. However, these actions are done
- // in case we pass on to the other
- // function.
+ // we have to do checkings on Insert/Add in any case to be consistent
+ // with the MPI communication model (see the comments in the
+ // documentation of TrilinosWrappers::Vector), but we can save some
+ // work if the addend is zero. However, these actions are done in case
+ // we pass on to the other function.
+
// TODO: fix this (do not run compress here, but fail)
if (last_action == Insert)
{
last_action = Add;
int *col_index_ptr;
- TrilinosScalar const *col_value_ptr;
+ TrilinosScalar *col_value_ptr;
int n_columns;
- // If we don't elide zeros, the pointers
- // are already available...
+ double short_val_array[100];
+ int short_index_array[100];
+ std::vector<TrilinosScalar> long_val_array;
+ std::vector<int> long_index_array;
+
+ // If we don't elide zeros, the pointers are already available... need to
+ // cast to non-const pointers as that is the format taken by Trilinos (but
+ // we will not modify const data)
if (elide_zero_values == false)
{
col_index_ptr = (int *)col_indices;
- col_value_ptr = values;
+ col_value_ptr = const_cast<TrilinosScalar*>(values);
n_columns = n_cols;
#ifdef DEBUG
for (unsigned int j=0; j<n_cols; ++j)
}
else
{
- // Otherwise, extract nonzero values in
- // each row and the corresponding index.
- if (column_indices.size() < n_cols)
+ // Otherwise, extract nonzero values in each row and the corresponding
+ // index.
+ if (n_cols > 100)
+ {
+ long_val_array.resize(n_cols);
+ long_index_array.resize(n_cols);
+ col_index_ptr = &long_index_array[0];
+ col_value_ptr = &long_val_array[0];
+ }
+ else
{
- column_indices.resize(n_cols);
- column_values.resize(n_cols);
+ col_index_ptr = &short_index_array[0];
+ col_value_ptr = &short_val_array[0];
}
n_columns = 0;
for (unsigned int j=0; j<n_cols; ++j)
{
const double value = values[j];
+
Assert (numbers::is_finite(value), ExcNumberNotFinite());
if (value != 0)
{
- column_indices[n_columns] = col_indices[j];
- column_values[n_columns] = value;
+ col_index_ptr[n_columns] = col_indices[j];
+ col_value_ptr[n_columns] = value;
n_columns++;
}
}
Assert(n_columns <= (int)n_cols, ExcInternalError());
- col_index_ptr = (int *)&column_indices[0];
- col_value_ptr = &column_values[0];
}
- // If the calling matrix owns the row to
- // which we want to add values, we
- // can directly call the Epetra_CrsMatrix
- // input function, which is much faster
- // than the Epetra_FECrsMatrix function.
+ // If the calling matrix owns the row to which we want to add values, we
+ // can directly call the Epetra_CrsMatrix input function, which is much
+ // faster than the Epetra_FECrsMatrix function.
if (row_partitioner().MyGID(static_cast<int>(row)) == true)
{
ierr = matrix->Epetra_CrsMatrix::SumIntoGlobalValues(row, n_columns,
- const_cast<double *>(col_value_ptr),
+ col_value_ptr,
col_index_ptr);
}
else
{
- // When we're at off-processor data, we
- // have to stick with the standard
- // SumIntoGlobalValues
- // function. Nevertheless, the way we
- // call it is the fastest one (any other
- // will lead to repeated allocation and
- // deallocation of memory in order to
- // call the function we already use,
- // which is very unefficient if writing
- // one element at a time).
+ // When we're at off-processor data, we have to stick with the
+ // standard SumIntoGlobalValues function. Nevertheless, the way we
+ // call it is the fastest one (any other will lead to repeated
+ // allocation and deallocation of memory in order to call the function
+ // we already use, which is very inefficient if writing one element at
+ // a time).
compressed = false;
ierr = matrix->SumIntoGlobalValues (1, (int *)&row, n_columns,
+ template <typename VectorType>
inline
void
- SparseMatrix::vmult_add (VectorBase &dst,
- const VectorBase &src) const
+ SparseMatrix::vmult_add (VectorType &dst,
+ const VectorType &src) const
{
Assert (&src != &dst, ExcSourceEqualsDestination());
- // Choose to reinit the vector with fast
- // argument set, which does not overwrite
- // the content -- this is what we need
- // since we're going to overwrite that
- // anyway in the vmult operation.
+ // 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);
+ template <typename VectorType>
inline
void
- SparseMatrix::Tvmult_add (VectorBase &dst,
- const VectorBase &src) const
+ SparseMatrix::Tvmult_add (VectorType &dst,
+ const VectorType &src) const
{
Assert (&src != &dst, ExcSourceEqualsDestination());
-
+ VectorType temp_vector;
temp_vector.reinit(dst, true);
Tvmult (temp_vector, src);
Assert (row_partitioner().SameAs(domain_partitioner()),
ExcNotQuadratic());
- temp_vector.reinit(v);
+ VectorBase temp_vector;
+ temp_vector.reinit(v, true);
vmult (temp_vector, v);
return temp_vector*v;
Assert (row_partitioner().SameAs(domain_partitioner()),
ExcNotQuadratic());
- temp_vector.reinit(v);
+ VectorBase temp_vector;
+ temp_vector.reinit(v, true);
vmult (temp_vector, v);
return u*temp_vector;
--- /dev/null
+//---------------------------- trilinos_block_sparse_matrix_vector_01.cc ---------------------------
+// $Id$
+// Version: $Name$
+//
+// Copyright (C) 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_block_sparse_matrix_vector_01.cc ---------------------------
+
+
+// check BlockSparseMatrix::vmult, Tvmult with deal.II block vector
+
+#include "../tests.h"
+#include <deal.II/base/utilities.h>
+#include <deal.II/lac/block_vector.h>
+#include <deal.II/lac/trilinos_block_sparse_matrix.h>
+#include <fstream>
+#include <iostream>
+#include <vector>
+
+
+void test (const unsigned int size_1, const unsigned int size_2)
+{
+ TrilinosWrappers::SparseMatrix m(size_1, size_2, size_2);
+ for (unsigned int i=0; i<m.m(); ++i)
+ for (unsigned int j=0; j<m.n(); ++j)
+ m.set (i,j, i+2*j);
+ m.compress (VectorOperation::insert);
+
+ TrilinosWrappers::BlockSparseMatrix m_block;
+ m_block.reinit(2,2);
+ m_block.block(0,0).copy_from(m);
+ m_block.block(1,0).copy_from(m);
+ m_block.block(0,1).copy_from(m);
+ m_block.block(1,1).copy_from(m);
+ m_block.collect_sizes();
+
+ BlockVector<double> v(2), w(2);
+ v.block(0).reinit(size_2);
+ v.block(1).reinit(size_2);
+ v.collect_sizes();
+ w.block(0).reinit(size_1);
+ w.block(1).reinit(size_1);
+ w.collect_sizes();
+ for (unsigned int i=0; i<v.size(); ++i)
+ v(i) = i;
+
+ // w:=Mv
+ m_block.vmult (w,v);
+
+ // make sure we get the expected result
+ for (unsigned int i=0; i<m.m(); ++i)
+ {
+ double result = 0;
+ for (unsigned int j=0; j<m.n(); ++j)
+ result += (i+2*j)*j;
+ for (unsigned int j=0; j<m.n(); ++j)
+ result += (i+2*j)*(j+m.n());
+ Assert (w(i) == result, ExcInternalError());
+ Assert (w(i+m.m()) == result, ExcInternalError());
+ }
+
+ for (unsigned int i=0; i<w.size(); ++i)
+ w(i) = i;
+
+ m_block.Tvmult (v, w);
+ // make sure we get the expected result
+ for (unsigned int i=0; i<m.n(); ++i)
+ {
+ double result = 0;
+ for (unsigned int j=0; j<m.m(); ++j)
+ result += (j+2*i)*j;
+ for (unsigned int j=0; j<m.m(); ++j)
+ result += (j+2*i)*(j+m.m());
+ Assert (v(i) == result, ExcInternalError());
+ Assert (v(i+m.n()) == result, ExcInternalError());
+ }
+
+ deallog << "OK" << std::endl;
+}
+
+
+
+int main (int argc, char **argv)
+{
+ std::ofstream logfile("block_sparse_matrix_vector_01/output");
+ deallog.attach(logfile);
+ deallog.depth_console(0);
+ deallog.threshold_double(1.e-10);
+
+ Utilities::MPI::MPI_InitFinalize mpi_initialization (argc, argv);
+
+
+ try
+ {
+ test (50,47);
+ }
+ 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
+//---------------------------- trilinos_block_sparse_matrix_vector_02.cc ---------------------------
+// $Id$
+// Version: $Name$
+//
+// Copyright (C) 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_block_sparse_matrix_vector_02.cc ---------------------------
+
+
+// check BlockSparseMatrix::vmult, Tvmult with deal.II block vector/vector
+// combination
+
+#include "../tests.h"
+#include <deal.II/base/utilities.h>
+#include <deal.II/lac/block_vector.h>
+#include <deal.II/lac/trilinos_block_sparse_matrix.h>
+#include <fstream>
+#include <iostream>
+#include <vector>
+
+
+void test (const unsigned int size_1, const unsigned int size_2)
+{
+ TrilinosWrappers::SparseMatrix m(size_1, size_2, size_2);
+ for (unsigned int i=0; i<m.m(); ++i)
+ for (unsigned int j=0; j<m.n(); ++j)
+ m.set (i,j, i+2*j);
+ m.compress (VectorOperation::insert);
+
+ TrilinosWrappers::BlockSparseMatrix m_block;
+ m_block.reinit(2,1);
+ m_block.block(0,0).copy_from(m);
+ m_block.block(1,0).copy_from(m);
+ m_block.collect_sizes();
+
+ BlockVector<double> w(2);
+ w.block(0).reinit(size_1);
+ w.block(1).reinit(size_1);
+ w.collect_sizes();
+ Vector<double>v(size_2);
+ for (unsigned int i=0; i<v.size(); ++i)
+ v(i) = i;
+
+ // w:=Mv
+ m_block.vmult (w,v);
+
+ // make sure we get the expected result
+ for (unsigned int i=0; i<m.m(); ++i)
+ {
+ double result = 0;
+ for (unsigned int j=0; j<m.n(); ++j)
+ result += (i+2*j)*j;
+ Assert (w(i) == result, ExcInternalError());
+ Assert (w(i+m.m()) == result, ExcInternalError());
+ }
+
+ for (unsigned int i=0; i<w.size(); ++i)
+ w(i) = i;
+
+ m_block.Tvmult (v, w);
+ // make sure we get the expected result
+ for (unsigned int i=0; i<m.n(); ++i)
+ {
+ double result = 0;
+ for (unsigned int j=0; j<m.m(); ++j)
+ result += (j+2*i)*j;
+ for (unsigned int j=0; j<m.m(); ++j)
+ result += (j+2*i)*(j+m.m());
+ Assert (v(i) == result, ExcInternalError());
+ }
+
+ deallog << "OK" << std::endl;
+}
+
+
+
+int main (int argc, char **argv)
+{
+ std::ofstream logfile("block_sparse_matrix_vector_02/output");
+ deallog.attach(logfile);
+ deallog.depth_console(0);
+ deallog.threshold_double(1.e-10);
+
+ Utilities::MPI::MPI_InitFinalize mpi_initialization (argc, argv);
+
+
+ try
+ {
+ test (50,47);
+ }
+ 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
+//---------------------------- trilinos_block_sparse_matrix_vector_03.cc ---------------------------
+// $Id$
+// Version: $Name$
+//
+// Copyright (C) 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_block_sparse_matrix_vector_03.cc ---------------------------
+
+
+// check BlockSparseMatrix::vmult, Tvmult with deal.II vector/block vector
+// combination
+
+#include "../tests.h"
+#include <deal.II/base/utilities.h>
+#include <deal.II/lac/block_vector.h>
+#include <deal.II/lac/trilinos_block_sparse_matrix.h>
+#include <fstream>
+#include <iostream>
+#include <vector>
+
+
+void test (const unsigned int size_1, const unsigned int size_2)
+{
+ TrilinosWrappers::SparseMatrix m(size_1, size_2, size_2);
+ for (unsigned int i=0; i<m.m(); ++i)
+ for (unsigned int j=0; j<m.n(); ++j)
+ m.set (i,j, i+2*j);
+ m.compress (VectorOperation::insert);
+
+ TrilinosWrappers::BlockSparseMatrix m_block;
+ m_block.reinit(1,2);
+ m_block.block(0,0).copy_from(m);
+ m_block.block(0,1).copy_from(m);
+ m_block.collect_sizes();
+
+ BlockVector<double> v(2);
+ v.block(0).reinit(size_2);
+ v.block(1).reinit(size_2);
+ v.collect_sizes();
+ Vector<double>w(size_1);
+ for (unsigned int i=0; i<v.size(); ++i)
+ v(i) = i;
+
+ // w:=Mv
+ m_block.vmult (w,v);
+
+ // make sure we get the expected result
+ for (unsigned int i=0; i<m.m(); ++i)
+ {
+ double result = 0;
+ for (unsigned int j=0; j<m.n(); ++j)
+ result += (i+2*j)*j;
+ for (unsigned int j=0; j<m.n(); ++j)
+ result += (i+2*j)*(j+m.n());
+ Assert (w(i) == result, ExcInternalError());
+ }
+
+ for (unsigned int i=0; i<w.size(); ++i)
+ w(i) = i;
+
+ m_block.Tvmult (v, w);
+ // make sure we get the expected result
+ for (unsigned int i=0; i<m.n(); ++i)
+ {
+ double result = 0;
+ for (unsigned int j=0; j<m.m(); ++j)
+ result += (j+2*i)*j;
+ Assert (v(i) == result, ExcInternalError());
+ Assert (v(i+m.n()) == result, ExcInternalError());
+ }
+
+ deallog << "OK" << std::endl;
+}
+
+
+
+int main (int argc, char **argv)
+{
+ std::ofstream logfile("block_sparse_matrix_vector_03/output");
+ deallog.attach(logfile);
+ deallog.depth_console(0);
+ deallog.threshold_double(1.e-10);
+
+ Utilities::MPI::MPI_InitFinalize mpi_initialization (argc, argv);
+
+
+ try
+ {
+ test (50,47);
+ }
+ 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
+//---------------------------- trilinos_block_sparse_matrix_vector_04.cc ---------------------------
+// $Id$
+// Version: $Name$
+//
+// Copyright (C) 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_block_sparse_matrix_vector_04.cc ---------------------------
+
+
+// check BlockSparseMatrix::vmult, Tvmult with deal.II vector/vector
+// combination
+
+#include "../tests.h"
+#include <deal.II/base/utilities.h>
+#include <deal.II/lac/block_vector.h>
+#include <deal.II/lac/trilinos_block_sparse_matrix.h>
+#include <fstream>
+#include <iostream>
+#include <vector>
+
+
+void test (const unsigned int size_1, const unsigned int size_2)
+{
+ TrilinosWrappers::SparseMatrix m(size_1, size_2, size_2);
+ for (unsigned int i=0; i<m.m(); ++i)
+ for (unsigned int j=0; j<m.n(); ++j)
+ m.set (i,j, i+2*j);
+ m.compress (VectorOperation::insert);
+
+ TrilinosWrappers::BlockSparseMatrix m_block;
+ m_block.reinit(1,1);
+ m_block.block(0,0).copy_from(m);
+ m_block.collect_sizes();
+
+ Vector<double> w(size_1), v(size_2);
+ for (unsigned int i=0; i<v.size(); ++i)
+ v(i) = i;
+
+ // w:=Mv
+ m_block.vmult (w,v);
+
+ // make sure we get the expected result
+ for (unsigned int i=0; i<m.m(); ++i)
+ {
+ double result = 0;
+ for (unsigned int j=0; j<m.n(); ++j)
+ result += (i+2*j)*j;
+ Assert (w(i) == result, ExcInternalError());
+ }
+
+ for (unsigned int i=0; i<w.size(); ++i)
+ w(i) = i;
+
+ m_block.Tvmult (v, w);
+ // make sure we get the expected result
+ for (unsigned int i=0; i<m.n(); ++i)
+ {
+ double result = 0;
+ for (unsigned int j=0; j<m.m(); ++j)
+ result += (j+2*i)*j;
+ Assert (v(i) == result, ExcInternalError());
+ }
+
+ deallog << "OK" << std::endl;
+}
+
+
+
+int main (int argc, char **argv)
+{
+ std::ofstream logfile("block_sparse_matrix_vector_04/output");
+ deallog.attach(logfile);
+ deallog.depth_console(0);
+ deallog.threshold_double(1.e-10);
+
+ Utilities::MPI::MPI_InitFinalize mpi_initialization (argc, argv);
+
+
+ try
+ {
+ test (50,47);
+ }
+ 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
+//---------------------------- trilinos_sparse_matrix_vector_08.cc ---------------------------
+// $Id$
+// Version: $Name$
+//
+// Copyright (C) 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_sparse_matrix_vector_08.cc ---------------------------
+
+
+// check SparseMatrix::vmult, vmult_add with deal.II vector
+
+#include "../tests.h"
+#include <deal.II/base/utilities.h>
+#include <deal.II/lac/vector.h>
+#include <deal.II/lac/trilinos_sparse_matrix.h>
+#include <fstream>
+#include <iostream>
+#include <vector>
+
+
+void test (Vector<double> &v,
+ Vector<double> &w)
+{
+ TrilinosWrappers::SparseMatrix m(w.size(),v.size(),v.size());
+ for (unsigned int i=0; i<m.m(); ++i)
+ for (unsigned int j=0; j<m.n(); ++j)
+ m.set (i,j, i+2*j);
+
+ for (unsigned int i=0; i<v.size(); ++i)
+ v(i) = i;
+
+ m.compress (VectorOperation::insert);
+
+ // w:=Mv
+ m.vmult (w,v);
+
+ // make sure we get the expected result
+ for (unsigned int i=0; i<m.m(); ++i)
+ {
+ double result = 0;
+ for (unsigned int j=0; j<m.n(); ++j)
+ result += (i+2*j)*j;
+ Assert (w(i) == result, ExcInternalError());
+ }
+
+ m.vmult_add (w, v);
+ // make sure we get the expected result
+ for (unsigned int i=0; i<m.m(); ++i)
+ {
+ double result = 0;
+ for (unsigned int j=0; j<m.n(); ++j)
+ result += (i+2*j)*j;
+ Assert (w(i) == result+result, ExcInternalError());
+ }
+
+ deallog << "OK" << std::endl;
+}
+
+
+
+int main (int argc, char **argv)
+{
+ std::ofstream logfile("sparse_matrix_vector_08/output");
+ deallog.attach(logfile);
+ deallog.depth_console(0);
+ deallog.threshold_double(1.e-10);
+
+ Utilities::MPI::MPI_InitFinalize mpi_initialization (argc, argv);
+
+
+ try
+ {
+ {
+ Vector<double> v (100);
+ Vector<double> w (95);
+ test (v,w);
+ }
+ }
+ 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
+//---------------------------- trilinos_sparse_matrix_vector_09.cc ---------------------------
+// $Id$
+// Version: $Name$
+//
+// Copyright (C) 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_sparse_matrix_vector_09.cc ---------------------------
+
+
+// check SparseMatrix::Tvmult, Tvmult_add with deal.II vector
+
+#include "../tests.h"
+#include <deal.II/base/utilities.h>
+#include <deal.II/lac/vector.h>
+#include <deal.II/lac/trilinos_sparse_matrix.h>
+#include <fstream>
+#include <iostream>
+#include <vector>
+
+
+void test (Vector<double> &v,
+ Vector<double> &w)
+{
+ TrilinosWrappers::SparseMatrix m(v.size(),w.size(),w.size());
+ for (unsigned int i=0; i<m.m(); ++i)
+ for (unsigned int j=0; j<m.n(); ++j)
+ m.set (i,j, i+2*j);
+
+ for (unsigned int i=0; i<v.size(); ++i)
+ v(i) = i;
+
+ m.compress (VectorOperation::insert);
+
+ // w:=Mv
+ m.Tvmult (w,v);
+
+ // make sure we get the expected result
+ for (unsigned int i=0; i<m.n(); ++i)
+ {
+ double result = 0;
+ for (unsigned int j=0; j<m.m(); ++j)
+ result += (j+2*i)*j;
+ Assert (w(i) == result, ExcInternalError());
+ }
+
+ m.Tvmult_add (w, v);
+ // make sure we get the expected result
+ for (unsigned int i=0; i<m.n(); ++i)
+ {
+ double result = 0;
+ for (unsigned int j=0; j<m.m(); ++j)
+ result += (j+2*i)*j;
+ Assert (w(i) == result+result, ExcInternalError());
+ }
+
+ deallog << "OK" << std::endl;
+}
+
+
+
+int main (int argc, char **argv)
+{
+ std::ofstream logfile("sparse_matrix_vector_09/output");
+ deallog.attach(logfile);
+ deallog.depth_console(0);
+ deallog.threshold_double(1.e-10);
+
+ Utilities::MPI::MPI_InitFinalize mpi_initialization (argc, argv);
+
+
+ try
+ {
+ {
+ Vector<double> v (95);
+ Vector<double> w (100);
+ test (v,w);
+ }
+ }
+ 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
+//---------------------------- trilinos_sparse_matrix_vector_10.cc ---------------------------
+// $Id$
+// Version: $Name$
+//
+// Copyright (C) 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_sparse_matrix_vector_10.cc ---------------------------
+
+
+// check SparseMatrix::vmult, vmult_add with distributed deal.II vector (but
+// without using distributed things)
+
+#include "../tests.h"
+#include <deal.II/base/utilities.h>
+#include <deal.II/lac/parallel_vector.h>
+#include <deal.II/lac/trilinos_sparse_matrix.h>
+#include <fstream>
+#include <iostream>
+#include <vector>
+
+
+void test (parallel::distributed::Vector<double> &v,
+ parallel::distributed::Vector<double> &w)
+{
+ TrilinosWrappers::SparseMatrix m(w.size(),v.size(),v.size());
+ for (unsigned int i=0; i<m.m(); ++i)
+ for (unsigned int j=0; j<m.n(); ++j)
+ m.set (i,j, i+2*j);
+
+ for (unsigned int i=0; i<v.size(); ++i)
+ v(i) = i;
+
+ m.compress (VectorOperation::insert);
+
+ // w:=Mv
+ m.vmult (w,v);
+
+ // make sure we get the expected result
+ for (unsigned int i=0; i<m.m(); ++i)
+ {
+ double result = 0;
+ for (unsigned int j=0; j<m.n(); ++j)
+ result += (i+2*j)*j;
+ Assert (w(i) == result, ExcInternalError());
+ }
+
+ m.vmult_add (w, v);
+ // make sure we get the expected result
+ for (unsigned int i=0; i<m.m(); ++i)
+ {
+ double result = 0;
+ for (unsigned int j=0; j<m.n(); ++j)
+ result += (i+2*j)*j;
+ Assert (w(i) == result+result, ExcInternalError());
+ }
+
+ deallog << "OK" << std::endl;
+}
+
+
+
+int main (int argc, char **argv)
+{
+ std::ofstream logfile("sparse_matrix_vector_10/output");
+ deallog.attach(logfile);
+ deallog.depth_console(0);
+ deallog.threshold_double(1.e-10);
+
+ Utilities::MPI::MPI_InitFinalize mpi_initialization (argc, argv);
+
+
+ try
+ {
+ {
+ parallel::distributed::Vector<double> v (100);
+ parallel::distributed::Vector<double> w (95);
+ test (v,w);
+ }
+ }
+ 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
+//---------------------------- trilinos_sparse_matrix_vector_11.cc ---------------------------
+// $Id$
+// Version: $Name$
+//
+// Copyright (C) 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_sparse_matrix_vector_11.cc ---------------------------
+
+
+// check SparseMatrix::Tvmult, Tvmult_add with distributed deal.II vector
+
+#include "../tests.h"
+#include <deal.II/base/utilities.h>
+#include <deal.II/lac/parallel_vector.h>
+#include <deal.II/lac/trilinos_sparse_matrix.h>
+#include <fstream>
+#include <iostream>
+#include <vector>
+
+
+void test (parallel::distributed::Vector<double> &v,
+ parallel::distributed::Vector<double> &w)
+{
+ TrilinosWrappers::SparseMatrix m(v.size(),w.size(),w.size());
+ for (unsigned int i=0; i<m.m(); ++i)
+ for (unsigned int j=0; j<m.n(); ++j)
+ m.set (i,j, i+2*j);
+
+ for (unsigned int i=0; i<v.size(); ++i)
+ v(i) = i;
+
+ m.compress (VectorOperation::insert);
+
+ // w:=Mv
+ m.Tvmult (w,v);
+
+ // make sure we get the expected result
+ for (unsigned int i=0; i<m.n(); ++i)
+ {
+ double result = 0;
+ for (unsigned int j=0; j<m.m(); ++j)
+ result += (j+2*i)*j;
+ Assert (w(i) == result, ExcInternalError());
+ }
+
+ m.Tvmult_add (w, v);
+ // make sure we get the expected result
+ for (unsigned int i=0; i<m.n(); ++i)
+ {
+ double result = 0;
+ for (unsigned int j=0; j<m.m(); ++j)
+ result += (j+2*i)*j;
+ Assert (w(i) == result+result, ExcInternalError());
+ }
+
+ deallog << "OK" << std::endl;
+}
+
+
+
+int main (int argc, char **argv)
+{
+ std::ofstream logfile("sparse_matrix_vector_11/output");
+ deallog.attach(logfile);
+ deallog.depth_console(0);
+ deallog.threshold_double(1.e-10);
+
+ Utilities::MPI::MPI_InitFinalize mpi_initialization (argc, argv);
+
+
+ try
+ {
+ {
+ parallel::distributed::Vector<double> v (95);
+ parallel::distributed::Vector<double> w (100);
+ test (v,w);
+ }
+ }
+ 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-svn.git / commitdiff