class SparseMatrix;
# ifndef DOXYGEN
+template <typename MatrixType>
+class BlockMatrixBase;
+
namespace LinearAlgebra
{
namespace TpetraWrappers
*/
using size_type = dealii::types::global_dof_index;
+ /**
+ * A structure that describes some of the traits of this class in terms of
+ * its run-time behavior. Some other classes (such as the block matrix
+ * classes) that take one or other of the matrix classes as its template
+ * parameters can tune their behavior based on the variables in this
+ * class.
+ */
+ struct Traits
+ {
+ /**
+ * It is safe to elide additions of zeros to individual elements of this
+ * matrix.
+ */
+ static const bool zero_addition_can_be_elided = true;
+ };
+
/**
* Declare an alias for the type used to store matrix elements, in analogy
* to all the other container classes.
SparseMatrix &
operator/=(const Number factor);
+ /**
+ * Copy the given (Trilinos) matrix (sparsity pattern and entries).
+ */
+ void
+ copy_from(const SparseMatrix<Number, MemorySpace> &source);
+
/**
* Add @p value to the element (<i>i,j</i>).
* Just as the respective call in deal.II SparseMatrix<Number,
const bool elide_zero_values = true,
const bool col_indices_are_sorted = false);
+ /**
+ * Add <tt>matrix</tt> scaled by <tt>factor</tt> to this matrix, i.e. the
+ * matrix <tt>factor*matrix</tt> is added to <tt>this</tt>. If the
+ * sparsity pattern of the calling matrix does not contain all the
+ * elements in the sparsity pattern of the input matrix, this function
+ * will throw an exception.
+ */
+ void
+ add(const Number factor, const SparseMatrix<Number, MemorySpace> &matrix);
+
/**
* Set the element (<i>i,j</i>) to @p value.
*
const Number *values,
const bool elide_zero_values = false);
+ /**
+ * Release all memory and return to a state just like after having called
+ * the default constructor.
+ *
+ * This is a @ref GlossCollectiveOperation "collective operation" that needs to be called on all
+ * processors in order to avoid a dead lock.
+ */
+ void
+ clear();
+
/** @} */
/**
* @name Entry Access
void
Tvmult_add(Vector<Number, MemorySpace> &dst,
const Vector<Number, MemorySpace> &src) const;
+
+ /**
+ * Compute the residual of an equation <i>Mx=b</i>, where the residual is
+ * defined to be <i>r=b-Mx</i>. Write the residual into @p dst. The
+ * <i>l<sub>2</sub></i> norm of the residual vector is returned.
+ *
+ * Source <i>x</i> and destination <i>dst</i> must not be the same vector.
+ *
+ * The vectors @p dst and @p b have to be initialized with the same
+ * IndexSet that was used for the row indices of the matrix and the vector
+ * @p x has to be initialized with the same IndexSet that was used for the
+ * column indices of the matrix.
+ *
+ * In case of a localized Vector, this function will only work when
+ * running on one processor, since the matrix object is inherently
+ * distributed. Otherwise, an exception will be thrown.
+ */
+ Number
+ residual(Vector<Number, MemorySpace> &dst,
+ const Vector<Number, MemorySpace> &x,
+ const Vector<Number, MemorySpace> &b) const;
+
/** @} */
/**
*/
bool compressed;
+ /**
+ * For some matrix storage formats, in particular for the PETSc
+ * distributed blockmatrices, set and add operations on individual
+ * elements can not be freely mixed. Rather, one has to synchronize
+ * operations when one wants to switch from setting elements to adding to
+ * elements. BlockMatrixBase automatically synchronizes the access by
+ * calling this helper function for each block. This function ensures
+ * that the matrix is in a state that allows adding elements; if it
+ * previously already was in this state, the function does nothing.
+ *
+ * This function is called from BlockMatrixBase.
+ */
+ void
+ prepare_add();
+
+ /**
+ * Same as prepare_add() but prepare the matrix for setting elements if
+ * the representation of elements in this class requires such an
+ * operation.
+ *
+ * This function is called from BlockMatrixBase.
+ */
+ void
+ prepare_set();
+
+ // To allow calling protected prepare_add() and prepare_set().
+ friend class BlockMatrixBase<SparseMatrix<Number, MemorySpace>>;
}; // class SparseMatrix
/* ------------------------- Inline functions ---------------------- */
+ template <typename Number, typename MemorySpace>
+ inline void
+ SparseMatrix<Number, MemorySpace>::set(const size_type i,
+ const size_type j,
+ const Number value)
+ {
+ set(i, 1, &j, &value, false);
+ }
+
+
+
template <typename Number, typename MemorySpace>
inline void
SparseMatrix<Number, MemorySpace>::add(const size_type i,
+ template <typename Number, typename MemorySpace>
+ inline Number
+ SparseMatrix<Number, MemorySpace>::residual(
+ Vector<Number, MemorySpace> &dst,
+ const Vector<Number, MemorySpace> &x,
+ const Vector<Number, MemorySpace> &b) const
+ {
+ vmult(dst, x);
+ dst -= b;
+ dst *= -1.;
+
+ return dst.l2_norm();
+ }
+
+
+
template <typename Number, typename MemorySpace>
inline dealii::types::signed_global_dof_index
SparseMatrix<Number, MemorySpace>::m() const
+ template <typename Number, typename MemorySpace>
+ inline void
+ SparseMatrix<Number, MemorySpace>::prepare_add()
+ {
+ // nothing to do here
+ }
+
+
+
+ template <typename Number, typename MemorySpace>
+ inline void
+ SparseMatrix<Number, MemorySpace>::prepare_set()
+ {
+ // nothing to do here
+ }
+
+
+
template <typename Number, typename MemorySpace>
inline const Tpetra::CrsMatrix<
Number,
// 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 the matrix is marked as compressed, we need to
+ // call resumeFill() first.
if (compressed || matrix->isFillComplete())
{
matrix->resumeFill();
- template <typename Number, typename MemorySpace>
- inline void
- SparseMatrix<Number, MemorySpace>::set(const size_type i,
- const size_type j,
- const Number value)
- {
- AssertIsFinite(value);
-
- set(i, 1, &j, &value, false);
- }
-
-
-
template <typename Number, typename MemorySpace>
inline void
SparseMatrix<Number, MemorySpace>::set(
+ template <typename Number, typename MemorySpace>
+ void
+ SparseMatrix<Number, MemorySpace>::add(
+ const Number factor,
+ const SparseMatrix<Number, MemorySpace> &source)
+ {
+ AssertDimension(source.m(), m());
+ AssertDimension(source.n(), n());
+ AssertDimension(source.local_range().first, local_range().first);
+ AssertDimension(source.local_range().second, local_range().second);
+ Assert(matrix->getRowMap()->isSameAs(*source.matrix->getRowMap()),
+ ExcMessage(
+ "Can only add matrices with same distribution of rows"));
+ Assert(matrix->isFillComplete() && source.matrix->isFillComplete(),
+ ExcMessage("Addition of matrices only allowed if matrices are "
+ "filled, i.e., compress() has been called"));
+
+ matrix->add(factor,
+ *source.matrix,
+ 1.0,
+ matrix->getDomainMap(),
+ matrix->getRangeMap(),
+ Teuchos::null);
+ }
+
+
+
+ template <typename Number, typename MemorySpace>
+ void
+ SparseMatrix<Number, MemorySpace>::copy_from(
+ const SparseMatrix<Number, MemorySpace> &source)
+ {
+ if (this == &source)
+ return;
+
+ // release memory before reallocation
+ matrix.reset();
+ column_space_map.reset();
+
+ // TODO:
+ // If the source and the target matrix have the same structure, we do
+ // not need to perform a deep copy.
+
+ // Perform a deep copy
+ matrix =
+ Utilities::Trilinos::internal::make_rcp<MatrixType>(*source.matrix,
+ Teuchos::Copy);
+ column_space_map = Teuchos::rcp_const_cast<MapType>(matrix->getColMap());
+ compressed = source.compressed;
+ }
+
+
+
+ template <typename Number, typename MemorySpace>
+ void
+ SparseMatrix<Number, MemorySpace>::clear()
+ {
+ // When we clear the matrix, reset
+ // the pointer and generate an
+ // empty matrix.
+ column_space_map = Utilities::Trilinos::internal::make_rcp<MapType>(
+ 0, 0, Utilities::Trilinos::tpetra_comm_self());
+
+ // Prepare the graph
+ Teuchos::RCP<GraphType> graph =
+ Utilities::Trilinos::internal::make_rcp<GraphType>(column_space_map,
+ column_space_map,
+ 0);
+ graph->fillComplete();
+
+ // Create the matrix from the graph
+ matrix = Utilities::Trilinos::internal::make_rcp<MatrixType>(graph);
+
+ compressed = true;
+ }
+
+
+
// Multiplications
template <typename Number, typename MemorySpace>
}
+
template <typename Number, typename MemorySpace>
void
SparseMatrix<Number, MemorySpace>::print(
projects / dealii.git / commitdiff