From 81ac29687c22b1abb2a6b4a4e9ce72f8ce6a0854 Mon Sep 17 00:00:00 2001 From: Daniel Arndt Date: Wed, 28 Feb 2024 09:14:25 -0500 Subject: [PATCH] TpetraWrappers: Implement iterators for SparseMatrix --- .../lac/trilinos_tpetra_sparse_matrix.h | 1061 ++++++++++++++++- .../lac/trilinos_tpetra_sparsity_pattern.cc | 14 +- .../trilinos_sparse_matrix_iterator_01.cc | 89 ++ .../trilinos_sparse_matrix_iterator_01.output | 5 + .../trilinos_sparse_matrix_iterator_02.cc | 96 ++ .../trilinos_sparse_matrix_iterator_02.output | 15 + .../trilinos_sparse_matrix_iterator_03.cc | 95 ++ .../trilinos_sparse_matrix_iterator_03.output | 15 + .../trilinos_sparse_matrix_iterator_04.cc | 93 ++ .../trilinos_sparse_matrix_iterator_04.output | 15 + .../trilinos_sparse_matrix_iterator_05.cc | 99 ++ .../trilinos_sparse_matrix_iterator_05.output | 15 + .../trilinos_sparse_matrix_iterator_06.cc | 99 ++ .../trilinos_sparse_matrix_iterator_06.output | 15 + .../trilinos_sparse_matrix_iterator_07.cc | 99 ++ .../trilinos_sparse_matrix_iterator_07.output | 15 + .../trilinos_sparse_matrix_iterator_08.cc | 99 ++ .../trilinos_sparse_matrix_iterator_08.output | 15 + .../trilinos_sparse_matrix_iterator_09.cc | 90 ++ .../trilinos_sparse_matrix_iterator_09.output | 3 + .../trilinos_sparse_matrix_iterator_10.cc | 92 ++ .../trilinos_sparse_matrix_iterator_10.output | 3 + .../trilinos_sparse_matrix_iterator_11.cc | 100 ++ .../trilinos_sparse_matrix_iterator_11.output | 2 + .../trilinos_sparse_matrix_iterator_12.cc | 100 ++ .../trilinos_sparse_matrix_iterator_12.output | 2 + 26 files changed, 2334 insertions(+), 12 deletions(-) create mode 100644 tests/trilinos_tpetra/trilinos_sparse_matrix_iterator_01.cc create mode 100644 tests/trilinos_tpetra/trilinos_sparse_matrix_iterator_01.output create mode 100644 tests/trilinos_tpetra/trilinos_sparse_matrix_iterator_02.cc create mode 100644 tests/trilinos_tpetra/trilinos_sparse_matrix_iterator_02.output create mode 100644 tests/trilinos_tpetra/trilinos_sparse_matrix_iterator_03.cc create mode 100644 tests/trilinos_tpetra/trilinos_sparse_matrix_iterator_03.output create mode 100644 tests/trilinos_tpetra/trilinos_sparse_matrix_iterator_04.cc create mode 100644 tests/trilinos_tpetra/trilinos_sparse_matrix_iterator_04.output create mode 100644 tests/trilinos_tpetra/trilinos_sparse_matrix_iterator_05.cc create mode 100644 tests/trilinos_tpetra/trilinos_sparse_matrix_iterator_05.output create mode 100644 tests/trilinos_tpetra/trilinos_sparse_matrix_iterator_06.cc create mode 100644 tests/trilinos_tpetra/trilinos_sparse_matrix_iterator_06.output create mode 100644 tests/trilinos_tpetra/trilinos_sparse_matrix_iterator_07.cc create mode 100644 tests/trilinos_tpetra/trilinos_sparse_matrix_iterator_07.output create mode 100644 tests/trilinos_tpetra/trilinos_sparse_matrix_iterator_08.cc create mode 100644 tests/trilinos_tpetra/trilinos_sparse_matrix_iterator_08.output create mode 100644 tests/trilinos_tpetra/trilinos_sparse_matrix_iterator_09.cc create mode 100644 tests/trilinos_tpetra/trilinos_sparse_matrix_iterator_09.output create mode 100644 tests/trilinos_tpetra/trilinos_sparse_matrix_iterator_10.cc create mode 100644 tests/trilinos_tpetra/trilinos_sparse_matrix_iterator_10.output create mode 100644 tests/trilinos_tpetra/trilinos_sparse_matrix_iterator_11.cc create mode 100644 tests/trilinos_tpetra/trilinos_sparse_matrix_iterator_11.output create mode 100644 tests/trilinos_tpetra/trilinos_sparse_matrix_iterator_12.cc create mode 100644 tests/trilinos_tpetra/trilinos_sparse_matrix_iterator_12.output diff --git a/include/deal.II/lac/trilinos_tpetra_sparse_matrix.h b/include/deal.II/lac/trilinos_tpetra_sparse_matrix.h index f29536400e..044a85d8db 100644 --- a/include/deal.II/lac/trilinos_tpetra_sparse_matrix.h +++ b/include/deal.II/lac/trilinos_tpetra_sparse_matrix.h @@ -50,6 +50,12 @@ namespace LinearAlgebra { template class SparsityPattern; + + namespace SparseMatrixIterators + { + template + class Iterator; + } } // namespace TpetraWrappers } // namespace LinearAlgebra # endif @@ -106,6 +112,16 @@ namespace LinearAlgebra */ using size_type = dealii::types::global_dof_index; + /** + * Exception + */ + DeclException1(ExcAccessToNonlocalRow, + std::size_t, + << "You tried to access row " << arg1 + << " of a non-contiguous locally owned row set." + << " The row " << arg1 + << " is not stored locally and can't be accessed."); + /** * 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 @@ -122,6 +138,18 @@ namespace LinearAlgebra static const bool zero_addition_can_be_elided = true; }; + /** + * Declare an alias for the iterator class. + */ + using iterator = + SparseMatrixIterators::Iterator; + + /** + * Declare an alias for the const iterator class. + */ + using const_iterator = + SparseMatrixIterators::Iterator; + /** * Declare an alias for the type used to store matrix elements, in analogy * to all the other container classes. @@ -437,6 +465,13 @@ namespace LinearAlgebra std::pair local_range() const; + /** + * Return whether @p index is in the local range or not, see also + * local_range(). + */ + bool + in_local_range(const size_type index) const; + /** * Return the total number of nonzero elements of this matrix (summed * over all MPI processes). @@ -444,6 +479,12 @@ namespace LinearAlgebra size_t n_nonzero_elements() const; + /** + * Number of entries in a specific row. + */ + unsigned int + row_length(const size_type row) const; + /** * Return the state of the matrix, i.e., whether compress() needs to be * called after an operation requiring data exchange. A call to compress() @@ -1043,6 +1084,107 @@ namespace LinearAlgebra /** @} */ + /** + * @name Iterators + */ + /** @{ */ + + /** + * Return an iterator pointing to the first element of the matrix. + * + * The elements accessed by iterators within each row are ordered in the + * way in which Trilinos stores them, though the implementation guarantees + * that all elements of one row are accessed before the elements of the + * next row. If your algorithm relies on visiting elements within one row, + * you will need to consult with the Trilinos documentation on the order + * in which it stores data. It is, however, generally not a good and + * long-term stable idea to rely on the order in which receive elements if + * you iterate over them. + * + * When you iterate over the elements of a parallel matrix, you will only + * be able to access the locally owned rows. (You can access the other + * rows as well, but they will look empty.) In that case, you probably + * want to call the begin() function that takes the row as an argument to + * limit the range of elements to loop over. + */ + const_iterator + begin() const; + + /** + * Like the function above, but for non-const matrices. + */ + iterator + begin(); + + /** + * Return an iterator pointing the element past the last one of this + * matrix. + */ + const_iterator + end() const; + + /** + * Like the function above, but for non-const matrices. + */ + iterator + end(); + + /** + * Return an iterator pointing to the first element of row @p r. + * + * Note that if the given row is empty, i.e. does not contain any nonzero + * entries, then the iterator returned by this function equals + * end(r). The returned iterator may not be dereferenceable in + * that case if neither row @p r nor any of the following rows contain any + * nonzero entries. + * + * The elements accessed by iterators within each row are ordered in the + * way in which Trilinos stores them, though the implementation guarantees + * that all elements of one row are accessed before the elements of the + * next row. If your algorithm relies on visiting elements within one row, + * you will need to consult with the Trilinos documentation on the order + * in which it stores data. It is, however, generally not a good and + * long-term stable idea to rely on the order in which receive elements if + * you iterate over them. + * + * @note When you access the elements of a parallel matrix, you can only + * access the elements of rows that are actually stored locally. (You can + * access the other rows as well, but they will look empty.) Even then, if + * another processor has since written into, or added to, an element of + * the matrix that is stored on the current processor, then you will still + * see the old value of this entry unless you have called compress() + * between modifying the matrix element on the remote processor and + * accessing it on the current processor. See the documentation of the + * compress() function for more information. + */ + const_iterator + begin(const size_type r) const; + + /** + * Like the function above, but for non-const matrices. + */ + iterator + begin(const size_type r); + + /** + * Return an iterator pointing the element past the last one of row @p r , + * or past the end of the entire sparsity pattern if none of the rows + * after @p r contain any entries at all. + * + * Note that the end iterator is not necessarily dereferenceable. This is + * in particular the case if it is the end iterator for the last row of a + * matrix. + */ + const_iterator + end(const size_type r) const; + + /** + * Like the function above, but for non-const matrices. + */ + iterator + end(const size_type r); + + /** @} */ /** * @addtogroup Exceptions */ @@ -1167,9 +1309,428 @@ namespace LinearAlgebra friend class BlockMatrixBase>; }; // class SparseMatrix + /** + * Iterators for Trilinos matrices + */ + namespace SparseMatrixIterators + { + /** + * Exception + */ + DeclException0(ExcBeyondEndOfMatrix); + + /** + * Handling of indices for both constant and non constant Accessor objects + * + * For a regular dealii::SparseMatrix, we would use an accessor for the + * sparsity pattern. For Trilinos matrices, this does not seem so simple, + * therefore, we write a little base class here. + */ + template + class AccessorBase + { + public: + /** + * Declare the type for container size. + */ + using size_type = dealii::types::global_dof_index; + + /** + * Constructor. + */ + AccessorBase(SparseMatrix *matrix, + const size_type row, + const size_type index); + + /** + * Row number of the element represented by this object. + */ + size_type + row() const; + + /** + * Index in row of the element represented by this object. + */ + size_type + index() const; + + /** + * Column number of the element represented by this object. + */ + size_type + column() const; + + protected: + /** + * Pointer to the matrix object. This object should be handled as a + * const pointer or non-const by the appropriate derived classes. In + * order to be able to implement both, it is not const here, so handle + * with care! + */ + mutable SparseMatrix *matrix; + /** + * Current row number. + */ + size_type a_row; + + /** + * Current index in row. + */ + size_type a_index; + + /** + * Discard the old row caches (they may still be used by other + * accessors) and generate new ones for the row pointed to presently by + * this accessor. + */ + void + visit_present_row(); + + /** + * Cache where we store the column indices of the present row. This is + * necessary, since Trilinos makes access to the elements of its + * matrices rather hard, and it is much more efficient to copy all + * column entries of a row once when we enter it than repeatedly asking + * Trilinos for individual ones. This also makes some sense since it is + * likely that we will access them sequentially anyway. + * + * In order to make copying of iterators/accessor of acceptable + * performance, we keep a shared pointer to these entries so that more + * than one accessor can access this data if necessary. + */ + std::shared_ptr> + colnum_cache; + + /** + * Cache for the values of this row. + */ + std::shared_ptr> value_cache; + + private: + friend class Iterator; + friend class Iterator; + }; + + /** + * General template for sparse matrix accessors. The first template + * argument denotes the underlying numeric type, the second the constness + * of the matrix. + * + * The general template is not implemented, only the specializations for + * the two possible values of the second template argument. Therefore, the + * interface listed here only serves as a template provided since doxygen + * does not link the specializations. + */ + template + class Accessor : public AccessorBase + { + /** + * Value of this matrix entry. + */ + Number + value() const; + + /** + * Value of this matrix entry. + */ + Number & + value(); + }; + + /** + * The specialization for a const Accessor. + */ + template + class Accessor + : public AccessorBase + { + public: + /** + * Typedef for the type (including constness) of the matrix to be used + * here. + */ + using MatrixType = const SparseMatrix; + + /** + * Typedef for the size type of the matrix to be used here. + */ + using size_type = typename AccessorBase::size_type; + + /** + * Constructor. Since we use accessors only for read access, a const + * matrix pointer is sufficient. + */ + Accessor(MatrixType *matrix, + const size_type row, + const size_type index); + + /** + * Copy constructor to get from a const or non-const accessor to a const + * accessor. + */ + template + Accessor(const Accessor &a); + + /** + * Value of this matrix entry. + */ + Number + value() const; + }; + + /** + * The specialization for a mutable Accessor. + */ + template + class Accessor + : public AccessorBase + { + class Reference + { + public: + /** + * Constructor. + */ + Reference(const Accessor &accessor); + + /** + * Conversion operator to the data type of the matrix. + */ + operator Number() const; + + /** + * Set the element of the matrix we presently point to to @p n. + */ + const Reference & + operator=(const Number n) const; + + /** + * Add @p n to the element of the matrix we presently point to. + */ + const Reference & + operator+=(const Number n) const; + + /** + * Subtract @p n from the element of the matrix we presently point to. + */ + const Reference & + operator-=(const Number n) const; + + /** + * Multiply the element of the matrix we presently point to by @p n. + */ + const Reference & + operator*=(const Number n) const; + + /** + * Divide the element of the matrix we presently point to by @p n. + */ + const Reference & + operator/=(const Number n) const; + + private: + /** + * Pointer to the accessor that denotes which element we presently + * point to. + */ + Accessor &accessor; + }; + + public: + /** + * Typedef for the type (including constness) of the matrix to be used + * here. + */ + using MatrixType = SparseMatrix; + + /** + * Typedef for the size type of the matrix to be used here. + */ + using size_type = typename AccessorBase::size_type; + + /** + * Constructor. Since we use accessors only for read access, a const + * matrix pointer is sufficient. + */ + Accessor(MatrixType *matrix, + const size_type row, + const size_type index); + + /** + * Value of this matrix entry. + */ + Reference + value() const; + + private: + // Make Reference object a friend. + friend class Reference; + }; + + /** + * This class acts as an iterator walking over the elements of Trilinos + * matrices. The implementation of this class is similar to the one for + * PETSc matrices. + * + * Note that Trilinos stores the elements within each row in ascending + * order. This is opposed to the deal.II sparse matrix style where the + * diagonal element (if it exists) is stored before all other values, and + * the PETSc sparse matrices, where one can't guarantee a certain order of + * the elements. + * + * @ingroup TpetraWrappers + */ + template + class Iterator + { + public: + /** + * Declare type for container size. + */ + using size_type = dealii::types::global_dof_index; + + /** + * A type that denotes what data types is used to express the difference + * between two iterators. + */ + using difference_type = dealii::types::global_dof_index; + + /** + * An alias for the type you get when you dereference an iterator of the + * current kind. + */ + using value_type = Number; + + /** + * Typedef for the matrix type (including constness) we are to operate + * on. + */ + using MatrixType = + typename Accessor::MatrixType; + + /** + * Constructor. Create an iterator into the matrix @p matrix for the + * given row and the index within it. + */ + Iterator(MatrixType *matrix, + const size_type row, + const size_type index); + + /** + * Copy constructor with optional change of constness. + */ + template + Iterator(const Iterator &other); + + /** + * Prefix increment. + */ + Iterator & + operator++(); + + /** + * Postfix increment. + */ + Iterator + operator++(int); + + /** + * Dereferencing operator. + */ + const Accessor & + operator*() const; + + /** + * Dereferencing operator. + */ + const Accessor * + operator->() const; + + /** + * Comparison. True, if both iterators point to the same matrix + * position. + */ + template + bool + operator==(const Iterator &) const; + + /** + * Inverse of ==. + */ + template + bool + operator!=(const Iterator &) const; + + /** + * Comparison operator. Result is true if either the first row number is + * smaller or if the row numbers are equal and the first index is + * smaller. + */ + template + bool + operator<(const Iterator &) const; + + /** + * Comparison operator. The opposite of the previous operator + */ + template + bool + operator>(const Iterator &) const; + + /** + * Exception + */ + DeclException2(ExcInvalidIndexWithinRow, + size_type, + size_type, + << "Attempt to access element " << arg2 << " of row " + << arg1 << " which doesn't have that many elements."); + + private: + /** + * Store an object of the accessor class. + */ + Accessor accessor; + + friend class Iterator; + friend class Iterator; + }; + } // namespace SparseMatrixIterators + + } // namespace TpetraWrappers + +} // namespace LinearAlgebra + - /* ------------------------- Inline functions ---------------------- */ +DEAL_II_NAMESPACE_CLOSE +namespace std +{ + template + struct iterator_traits< + dealii::LinearAlgebra::TpetraWrappers::SparseMatrixIterators:: + Iterator> + { + using iterator_category = forward_iterator_tag; + using value_type = + typename dealii::LinearAlgebra::TpetraWrappers::SparseMatrixIterators:: + Iterator::value_type; + using difference_type = + typename dealii::LinearAlgebra::TpetraWrappers::SparseMatrixIterators:: + Iterator::difference_type; + }; +} // namespace std + +/* ------------------------- Inline functions ---------------------- */ + + +DEAL_II_NAMESPACE_OPEN + +namespace LinearAlgebra +{ + + namespace TpetraWrappers + { template inline void SparseMatrix::set(const size_type i, @@ -1239,6 +1800,134 @@ namespace LinearAlgebra + template + inline typename SparseMatrix::const_iterator + SparseMatrix::begin() const + { + return begin(0); + } + + + + template + inline typename SparseMatrix::const_iterator + SparseMatrix::end() const + { + return const_iterator(this, m(), 0); + } + + + + template + inline typename SparseMatrix::const_iterator + SparseMatrix::begin(const size_type r) const + { + AssertIndexRange(r, m()); + if (in_local_range(r) && (row_length(r) > 0)) + return const_iterator(this, r, 0); + else + return end(r); + } + + + template + inline typename SparseMatrix::const_iterator + SparseMatrix::end(const size_type r) const + { + AssertIndexRange(r, m()); + + // place the iterator on the first entry + // past this line, or at the end of the + // matrix + for (size_type i = r + 1; i < m(); ++i) + if (in_local_range(i) && (row_length(i) > 0)) + return const_iterator(this, i, 0); + + // if there is no such line, then take the + // end iterator of the matrix + return end(); + } + + + + template + inline typename SparseMatrix::iterator + SparseMatrix::begin() + { + return begin(0); + } + + + + template + inline typename SparseMatrix::iterator + SparseMatrix::end() + { + return iterator(this, m(), 0); + } + + + + template + inline typename SparseMatrix::iterator + SparseMatrix::begin(const size_type r) + { + AssertIndexRange(r, m()); + if (in_local_range(r) && (row_length(r) > 0)) + return iterator(this, r, 0); + else + return end(r); + } + + + + template + inline typename SparseMatrix::iterator + SparseMatrix::end(const size_type r) + { + AssertIndexRange(r, m()); + + // place the iterator on the first entry + // past this line, or at the end of the + // matrix + for (size_type i = r + 1; i < m(); ++i) + if (in_local_range(i) && (row_length(i) > 0)) + return iterator(this, i, 0); + + // if there is no such line, then take the + // end iterator of the matrix + return end(); + } + + + + template + inline bool + SparseMatrix::in_local_range( + const size_type index) const + { + auto begin = matrix->getRowMap()->getMinGlobalIndex(); + auto end = matrix->getRowMap()->getMaxGlobalIndex() + 1; + + return ((index >= begin) && (index < end)); + } + + + + template + unsigned int + SparseMatrix::row_length(const size_type row) const + { + auto n_entries = matrix->getNumEntriesInGlobalRow(row); + Assert(n_entries != + Teuchos::OrdinalTraits::invalid(), + ExcAccessToNonlocalRow(row)); + + return n_entries; + } + + + template inline void SparseMatrix::prepare_add() @@ -1325,7 +2014,375 @@ namespace LinearAlgebra return IndexSet(matrix->getRangeMap()); } - } // namespace TpetraWrappers + + namespace SparseMatrixIterators + { + template + inline AccessorBase::AccessorBase( + SparseMatrix *matrix, + size_type row, + size_type index) + : matrix(matrix) + , a_row(row) + , a_index(index) + { + visit_present_row(); + } + + + template + inline typename AccessorBase::size_type + AccessorBase::row() const + { + Assert(a_row < matrix->m(), ExcBeyondEndOfMatrix()); + return a_row; + } + + + template + inline typename AccessorBase::size_type + AccessorBase::column() const + { + Assert(a_row < matrix->m(), ExcBeyondEndOfMatrix()); + return (*colnum_cache)[a_index]; + } + + + template + inline typename AccessorBase::size_type + AccessorBase::index() const + { + Assert(a_row < matrix->m(), ExcBeyondEndOfMatrix()); + return a_index; + } + + + + template + void + AccessorBase::visit_present_row() + { + // if we are asked to visit the past-the-end line, then simply + // release all our caches and go on with life. + // + // do the same if the row we're supposed to visit is not locally + // owned. this is simply going to make non-locally owned rows + // look like they're empty + if ((this->a_row == matrix->m()) || + (matrix->in_local_range(this->a_row) == false)) + { + colnum_cache.reset(); + value_cache.reset(); + + return; + } + + // get a representation of the present row + size_t ncols; + TrilinosWrappers::types::int_type colnums = + matrix->row_length(this->a_row); + if (value_cache.get() == nullptr) + { + value_cache = std::make_shared>(colnums); + colnum_cache = std::make_shared< + std::vector>(colnums); + } + else + { + value_cache->resize(colnums); + colnum_cache->resize(colnums); + } + +# if DEAL_II_TRILINOS_VERSION_GTE(13, 2, 0) + typename SparseMatrix::MatrixType:: + nonconst_global_inds_host_view_type col_indices(colnum_cache->data(), + colnums); + typename SparseMatrix::MatrixType:: + nonconst_values_host_view_type values(value_cache->data(), colnums); +# else + Teuchos::ArrayView col_indices( + *colnum_cache); + Teuchos::ArrayView values(*value_cache); +# endif + matrix->trilinos_matrix().getGlobalRowCopy(this->a_row, + col_indices, + values, + ncols); + + AssertDimension(ncols, colnums); + + // copy it into our caches if the + // line isn't empty. if it is, then + // we've done something wrong, since + // we shouldn't have initialized an + // iterator for an empty line (what + // would it point to?) + } + + + + template + inline Accessor::Accessor( + MatrixType *matrix, + const size_type row, + const size_type index) + : AccessorBase( + const_cast *>(matrix), + row, + index) + {} + + + + template + template + inline Accessor::Accessor( + const Accessor &other) + : AccessorBase(other) + {} + + + + template + inline Number + Accessor::value() const + { + Assert((AccessorBase::a_row < + AccessorBase::matrix->m()), + ExcBeyondEndOfMatrix()); + return (*AccessorBase::value_cache) + [AccessorBase::a_index]; + } + + + + template + inline Accessor::Reference::Reference( + const Accessor &acc) + : accessor(const_cast &>(acc)) + {} + + + + template + inline Accessor::Reference::operator Number() + const + { + return (*accessor.value_cache)[accessor.a_index]; + } + + + + template + inline const typename Accessor::Reference & + Accessor::Reference::operator=( + const Number n) const + { + (*accessor.value_cache)[accessor.a_index] = n; + accessor.matrix->set(accessor.row(), + accessor.column(), + static_cast(*this)); + return *this; + } + + + + template + inline const typename Accessor::Reference & + Accessor::Reference::operator+=( + const Number n) const + { + (*accessor.value_cache)[accessor.a_index] += n; + accessor.matrix->set(accessor.row(), + accessor.column(), + static_cast(*this)); + return *this; + } + + + + template + inline const typename Accessor::Reference & + Accessor::Reference::operator-=( + const Number n) const + { + (*accessor.value_cache)[accessor.a_index] -= n; + accessor.matrix->set(accessor.row(), + accessor.column(), + static_cast(*this)); + return *this; + } + + + template + inline const typename Accessor::Reference & + Accessor::Reference::operator*=( + const Number n) const + { + (*accessor.value_cache)[accessor.a_index] *= n; + accessor.matrix->set(accessor.row(), + accessor.column(), + static_cast(*this)); + return *this; + } + + + template + inline const typename Accessor::Reference & + Accessor::Reference::operator/=( + const Number n) const + { + (*accessor.value_cache)[accessor.a_index] /= n; + accessor.matrix->set(accessor.row(), + accessor.column(), + static_cast(*this)); + return *this; + } + + + template + inline Accessor::Accessor( + MatrixType *matrix, + const size_type row, + const size_type index) + : AccessorBase(matrix, row, index) + {} + + + template + inline typename Accessor::Reference + Accessor::value() const + { + Assert((AccessorBase::a_row < + AccessorBase::matrix->m()), + ExcBeyondEndOfMatrix()); + return {*this}; + } + + + + template + inline Iterator::Iterator( + MatrixType *matrix, + const size_type row, + const size_type index) + : accessor(matrix, row, index) + {} + + + template + template + inline Iterator::Iterator( + const Iterator &other) + : accessor(other.accessor) + {} + + + + template + inline Iterator & + Iterator::operator++() + { + Assert(accessor.a_row < accessor.matrix->m(), ExcIteratorPastEnd()); + + ++accessor.a_index; + + // If at end of line: do one + // step, then cycle until we + // find a row with a nonzero + // number of entries. + if (accessor.a_index >= accessor.colnum_cache->size()) + { + accessor.a_index = 0; + ++accessor.a_row; + + while ( + (accessor.a_row < accessor.matrix->m()) && + ((accessor.matrix->in_local_range(accessor.a_row) == false) || + (accessor.matrix->row_length(accessor.a_row) == 0))) + ++accessor.a_row; + + accessor.visit_present_row(); + } + return *this; + } + + + template + inline Iterator + Iterator::operator++(int) + { + const Iterator old_state = *this; + ++(*this); + return old_state; + } + + + + template + inline const Accessor & + Iterator::operator*() const + { + return accessor; + } + + + + template + inline const Accessor * + Iterator::operator->() const + { + return &accessor; + } + + + + template + template + inline bool + Iterator::operator==( + const Iterator &other) const + { + return (accessor.a_row == other.accessor.a_row && + accessor.a_index == other.accessor.a_index); + } + + + + template + template + inline bool + Iterator::operator!=( + const Iterator &other) const + { + return !(*this == other); + } + + + + template + template + inline bool + Iterator::operator<( + const Iterator &other) const + { + return (accessor.row() < other.accessor.row() || + (accessor.row() == other.accessor.row() && + accessor.index() < other.accessor.index())); + } + + + template + template + inline bool + Iterator::operator>( + const Iterator &other) const + { + return (other < *this); + } + + } // namespace SparseMatrixIterators + } // namespace TpetraWrappers } // namespace LinearAlgebra diff --git a/source/lac/trilinos_tpetra_sparsity_pattern.cc b/source/lac/trilinos_tpetra_sparsity_pattern.cc index a3d4cb41c9..488d757267 100644 --- a/source/lac/trilinos_tpetra_sparsity_pattern.cc +++ b/source/lac/trilinos_tpetra_sparsity_pattern.cc @@ -778,21 +778,15 @@ namespace LinearAlgebra ExcInternalError()); # endif - nonlocal_graph->fillComplete(column_space_map, graph->getRangeMap()); - graph->fillComplete(column_space_map, graph->getRangeMap()); - } - else - { - graph->globalAssemble(); + nonlocal_graph->fillComplete(column_space_map, graph->getRowMap()); } + graph->fillComplete(column_space_map, graph->getRowMap()); // Check consistency between the sizes set at the beginning and what // Trilinos stores: using namespace deal_II_exceptions::internals; - Assert(compare_for_equality(n_rows(), graph->getGlobalNumEntries()), - ExcInternalError()); - Assert(compare_for_equality(n_cols(), graph->getGlobalNumEntries()), - ExcInternalError()); + AssertDimension(n_rows(), graph->getGlobalNumRows()); + AssertDimension(n_cols(), graph->getGlobalNumCols()); } diff --git a/tests/trilinos_tpetra/trilinos_sparse_matrix_iterator_01.cc b/tests/trilinos_tpetra/trilinos_sparse_matrix_iterator_01.cc new file mode 100644 index 0000000000..0487caaa6c --- /dev/null +++ b/tests/trilinos_tpetra/trilinos_sparse_matrix_iterator_01.cc @@ -0,0 +1,89 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2004 - 2018 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.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + + + +// LinearAlgebra::TpetraWrappers::SparseMatrix::const_iterator::operator++(int) +// was declared but not defined for some time, leading to a linker error + +#include +#include + +#include "../tests.h" + + +void +test() +{ + LinearAlgebra::TpetraWrappers::SparsityPattern sp(5, 5, 3); + for (unsigned int i = 0; i < 5; ++i) + for (unsigned int j = 0; j < 5; ++j) + if ((i + 2 * j + 1) % 3 == 0) + sp.add(i, j); + sp.compress(); + + LinearAlgebra::TpetraWrappers::SparseMatrix m(sp); + LinearAlgebra::TpetraWrappers::SparseMatrix::const_iterator i = + m.begin(); + deallog << i->value() << std::endl; + ++i; + deallog << i->value() << std::endl; + i++; + deallog << i->value() << std::endl; + + deallog << "OK" << std::endl; +} + + + +int +main(int argc, char **argv) +{ + initlog(); + + Utilities::MPI::MPI_InitFinalize mpi_initialization( + argc, argv, testing_max_num_threads()); + + try + { + test(); + } + catch (const std::exception &exc) + { + deallog << std::endl + << std::endl + << "----------------------------------------------------" + << std::endl; + deallog << "Exception on processing: " << std::endl + << exc.what() << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + + return 1; + } + catch (...) + { + deallog << std::endl + << std::endl + << "----------------------------------------------------" + << std::endl; + deallog << "Unknown exception!" << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + return 1; + }; +} diff --git a/tests/trilinos_tpetra/trilinos_sparse_matrix_iterator_01.output b/tests/trilinos_tpetra/trilinos_sparse_matrix_iterator_01.output new file mode 100644 index 0000000000..df15178e0e --- /dev/null +++ b/tests/trilinos_tpetra/trilinos_sparse_matrix_iterator_01.output @@ -0,0 +1,5 @@ + +DEAL::0 +DEAL::0 +DEAL::0 +DEAL::OK diff --git a/tests/trilinos_tpetra/trilinos_sparse_matrix_iterator_02.cc b/tests/trilinos_tpetra/trilinos_sparse_matrix_iterator_02.cc new file mode 100644 index 0000000000..e08ba5e997 --- /dev/null +++ b/tests/trilinos_tpetra/trilinos_sparse_matrix_iterator_02.cc @@ -0,0 +1,96 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2004 - 2018 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.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + + + +// test setting some elements and reading them back from a const matrix +// iterator + +#include +#include + +#include "../tests.h" + + +void +test() +{ + LinearAlgebra::TpetraWrappers::SparsityPattern sp(5, 5, 3); + for (unsigned int i = 0; i < 5; ++i) + for (unsigned int j = 0; j < 5; ++j) + if (((i + 2 * j + 1) % 3 == 0) || (i == j)) + sp.add(i, j); + sp.compress(); + + LinearAlgebra::TpetraWrappers::SparseMatrix m(sp); + for (unsigned int i = 0; i < 5; ++i) + for (unsigned int j = 0; j < 5; ++j) + if (((i + 2 * j + 1) % 3 == 0) || (i == j)) + m.set(i, j, i * j); + + LinearAlgebra::TpetraWrappers::SparseMatrix::const_iterator i = + m.begin(); + for (; i != m.end(); ++i) + { + deallog << i->row() << ' ' << i->column() << ' ' << i->value() + << std::endl; + Assert(std::fabs(i->value() - i->row() * i->column()) < 1e-14, + ExcInternalError()); + } + + deallog << "OK" << std::endl; +} + + + +int +main(int argc, char **argv) +{ + initlog(); + + Utilities::MPI::MPI_InitFinalize mpi_initialization( + argc, argv, testing_max_num_threads()); + + try + { + test(); + } + catch (const std::exception &exc) + { + deallog << std::endl + << std::endl + << "----------------------------------------------------" + << std::endl; + deallog << "Exception on processing: " << std::endl + << exc.what() << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + + return 1; + } + catch (...) + { + deallog << std::endl + << std::endl + << "----------------------------------------------------" + << std::endl; + deallog << "Unknown exception!" << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + return 1; + }; +} diff --git a/tests/trilinos_tpetra/trilinos_sparse_matrix_iterator_02.output b/tests/trilinos_tpetra/trilinos_sparse_matrix_iterator_02.output new file mode 100644 index 0000000000..3262f31a6d --- /dev/null +++ b/tests/trilinos_tpetra/trilinos_sparse_matrix_iterator_02.output @@ -0,0 +1,15 @@ + +DEAL::0 0 0 +DEAL::0 1 0 +DEAL::0 4 0 +DEAL::1 1 1.00000 +DEAL::1 2 2.00000 +DEAL::2 0 0 +DEAL::2 2 4.00000 +DEAL::2 3 6.00000 +DEAL::3 1 3.00000 +DEAL::3 3 9.00000 +DEAL::3 4 12.0000 +DEAL::4 2 8.00000 +DEAL::4 4 16.0000 +DEAL::OK diff --git a/tests/trilinos_tpetra/trilinos_sparse_matrix_iterator_03.cc b/tests/trilinos_tpetra/trilinos_sparse_matrix_iterator_03.cc new file mode 100644 index 0000000000..3646f465e8 --- /dev/null +++ b/tests/trilinos_tpetra/trilinos_sparse_matrix_iterator_03.cc @@ -0,0 +1,95 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2004 - 2018 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.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + + + +// test setting some elements and reading them back from a non-const matrix +// iterator + +#include +#include + +#include "../tests.h" + + +void +test() +{ + LinearAlgebra::TpetraWrappers::SparsityPattern sp(5, 5, 3); + for (unsigned int i = 0; i < 5; ++i) + for (unsigned int j = 0; j < 5; ++j) + if (((i + 2 * j + 1) % 3 == 0) || (i == j)) + sp.add(i, j); + sp.compress(); + + LinearAlgebra::TpetraWrappers::SparseMatrix m(sp); + for (unsigned int i = 0; i < 5; ++i) + for (unsigned int j = 0; j < 5; ++j) + if (((i + 2 * j + 1) % 3 == 0) || (i == j)) + m.set(i, j, i * j); + + LinearAlgebra::TpetraWrappers::SparseMatrix::iterator i = m.begin(); + for (; i != m.end(); ++i) + { + deallog << i->row() << ' ' << i->column() << ' ' << i->value() + << std::endl; + Assert(std::fabs(i->value() - i->row() * i->column()) < 1e-14, + ExcInternalError()); + } + + deallog << "OK" << std::endl; +} + + + +int +main(int argc, char **argv) +{ + initlog(); + + Utilities::MPI::MPI_InitFinalize mpi_initialization( + argc, argv, testing_max_num_threads()); + + try + { + test(); + } + catch (const std::exception &exc) + { + deallog << std::endl + << std::endl + << "----------------------------------------------------" + << std::endl; + deallog << "Exception on processing: " << std::endl + << exc.what() << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + + return 1; + } + catch (...) + { + deallog << std::endl + << std::endl + << "----------------------------------------------------" + << std::endl; + deallog << "Unknown exception!" << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + return 1; + }; +} diff --git a/tests/trilinos_tpetra/trilinos_sparse_matrix_iterator_03.output b/tests/trilinos_tpetra/trilinos_sparse_matrix_iterator_03.output new file mode 100644 index 0000000000..50fad50cd5 --- /dev/null +++ b/tests/trilinos_tpetra/trilinos_sparse_matrix_iterator_03.output @@ -0,0 +1,15 @@ + +DEAL::0 0 0.00000 +DEAL::0 1 0.00000 +DEAL::0 4 0.00000 +DEAL::1 1 1.00000 +DEAL::1 2 2.00000 +DEAL::2 0 0.00000 +DEAL::2 2 4.00000 +DEAL::2 3 6.00000 +DEAL::3 1 3.00000 +DEAL::3 3 9.00000 +DEAL::3 4 12.0000 +DEAL::4 2 8.00000 +DEAL::4 4 16.0000 +DEAL::OK diff --git a/tests/trilinos_tpetra/trilinos_sparse_matrix_iterator_04.cc b/tests/trilinos_tpetra/trilinos_sparse_matrix_iterator_04.cc new file mode 100644 index 0000000000..a9e2d582a2 --- /dev/null +++ b/tests/trilinos_tpetra/trilinos_sparse_matrix_iterator_04.cc @@ -0,0 +1,93 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2004 - 2018 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.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + + +// test setting some elements using a non-const matrix iterator and operator=, +// and reading them back through the matrix itself + +#include +#include + +#include "../tests.h" + + +void +test() +{ + LinearAlgebra::TpetraWrappers::SparsityPattern sp(5, 5, 3); + for (unsigned int i = 0; i < 5; ++i) + for (unsigned int j = 0; j < 5; ++j) + if (((i + 2 * j + 1) % 3 == 0) || (i == j)) + sp.add(i, j); + sp.compress(); + LinearAlgebra::TpetraWrappers::SparseMatrix m(sp); + + LinearAlgebra::TpetraWrappers::SparseMatrix::iterator i = m.begin(); + for (; i != m.end(); ++i) + i->value() = i->row() * i->column(); + m.compress(VectorOperation::insert); + + for (unsigned int i = 0; i < 5; ++i) + for (unsigned int j = 0; j < 5; ++j) + if (((i + 2 * j + 1) % 3 == 0) || (i == j)) + { + deallog << i << ' ' << j << ' ' << m.el(i, j) << std::endl; + Assert(std::fabs(m.el(i, j) - i * j) < 1e-14, ExcInternalError()); + } + + deallog << "OK" << std::endl; +} + + + +int +main(int argc, char **argv) +{ + initlog(); + + Utilities::MPI::MPI_InitFinalize mpi_initialization( + argc, argv, testing_max_num_threads()); + + try + { + test(); + } + catch (const std::exception &exc) + { + deallog << std::endl + << std::endl + << "----------------------------------------------------" + << std::endl; + deallog << "Exception on processing: " << std::endl + << exc.what() << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + + return 1; + } + catch (...) + { + deallog << std::endl + << std::endl + << "----------------------------------------------------" + << std::endl; + deallog << "Unknown exception!" << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + return 1; + }; +} diff --git a/tests/trilinos_tpetra/trilinos_sparse_matrix_iterator_04.output b/tests/trilinos_tpetra/trilinos_sparse_matrix_iterator_04.output new file mode 100644 index 0000000000..3262f31a6d --- /dev/null +++ b/tests/trilinos_tpetra/trilinos_sparse_matrix_iterator_04.output @@ -0,0 +1,15 @@ + +DEAL::0 0 0 +DEAL::0 1 0 +DEAL::0 4 0 +DEAL::1 1 1.00000 +DEAL::1 2 2.00000 +DEAL::2 0 0 +DEAL::2 2 4.00000 +DEAL::2 3 6.00000 +DEAL::3 1 3.00000 +DEAL::3 3 9.00000 +DEAL::3 4 12.0000 +DEAL::4 2 8.00000 +DEAL::4 4 16.0000 +DEAL::OK diff --git a/tests/trilinos_tpetra/trilinos_sparse_matrix_iterator_05.cc b/tests/trilinos_tpetra/trilinos_sparse_matrix_iterator_05.cc new file mode 100644 index 0000000000..a91dc2e4c0 --- /dev/null +++ b/tests/trilinos_tpetra/trilinos_sparse_matrix_iterator_05.cc @@ -0,0 +1,99 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2004 - 2018 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.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + + +// test setting some elements using a non-const matrix iterator and operator+=, +// and reading them back through the matrix itself + +#include +#include + +#include "../tests.h" + + +void +test() +{ + LinearAlgebra::TpetraWrappers::SparsityPattern sp(5, 5, 3); + for (unsigned int i = 0; i < 5; ++i) + for (unsigned int j = 0; j < 5; ++j) + if (((i + 2 * j + 1) % 3 == 0) || (i == j)) + sp.add(i, j); + sp.compress(); + + LinearAlgebra::TpetraWrappers::SparseMatrix m(sp); + for (unsigned int i = 0; i < 5; ++i) + for (unsigned int j = 0; j < 5; ++j) + if (((i + 2 * j + 1) % 3 == 0) || (i == j)) + m.set(i, j, 1.); + + LinearAlgebra::TpetraWrappers::SparseMatrix::iterator i = m.begin(); + for (; i != m.end(); ++i) + i->value() += i->row() * i->column(); + m.compress(VectorOperation::insert); + + for (unsigned int i = 0; i < 5; ++i) + for (unsigned int j = 0; j < 5; ++j) + if (((i + 2 * j + 1) % 3 == 0) || (i == j)) + { + deallog << i << ' ' << j << ' ' << m.el(i, j) << std::endl; + Assert(std::fabs(m.el(i, j) - (1. + i * j)) < 1e-14, + ExcInternalError()); + } + + deallog << "OK" << std::endl; +} + + + +int +main(int argc, char **argv) +{ + initlog(); + + Utilities::MPI::MPI_InitFinalize mpi_initialization( + argc, argv, testing_max_num_threads()); + + try + { + test(); + } + catch (const std::exception &exc) + { + deallog << std::endl + << std::endl + << "----------------------------------------------------" + << std::endl; + deallog << "Exception on processing: " << std::endl + << exc.what() << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + + return 1; + } + catch (...) + { + deallog << std::endl + << std::endl + << "----------------------------------------------------" + << std::endl; + deallog << "Unknown exception!" << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + return 1; + }; +} diff --git a/tests/trilinos_tpetra/trilinos_sparse_matrix_iterator_05.output b/tests/trilinos_tpetra/trilinos_sparse_matrix_iterator_05.output new file mode 100644 index 0000000000..f07ba8bff9 --- /dev/null +++ b/tests/trilinos_tpetra/trilinos_sparse_matrix_iterator_05.output @@ -0,0 +1,15 @@ + +DEAL::0 0 1.00000 +DEAL::0 1 1.00000 +DEAL::0 4 1.00000 +DEAL::1 1 2.00000 +DEAL::1 2 3.00000 +DEAL::2 0 1.00000 +DEAL::2 2 5.00000 +DEAL::2 3 7.00000 +DEAL::3 1 4.00000 +DEAL::3 3 10.0000 +DEAL::3 4 13.0000 +DEAL::4 2 9.00000 +DEAL::4 4 17.0000 +DEAL::OK diff --git a/tests/trilinos_tpetra/trilinos_sparse_matrix_iterator_06.cc b/tests/trilinos_tpetra/trilinos_sparse_matrix_iterator_06.cc new file mode 100644 index 0000000000..3fb1a0c420 --- /dev/null +++ b/tests/trilinos_tpetra/trilinos_sparse_matrix_iterator_06.cc @@ -0,0 +1,99 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2004 - 2018 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.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + + +// test setting some elements using a non-const matrix iterator and operator-=, +// and reading them back through the matrix itself + +#include +#include + +#include "../tests.h" + + +void +test() +{ + LinearAlgebra::TpetraWrappers::SparsityPattern sp(5, 5, 3); + for (unsigned int i = 0; i < 5; ++i) + for (unsigned int j = 0; j < 5; ++j) + if (((i + 2 * j + 1) % 3 == 0) || (i == j)) + sp.add(i, j); + sp.compress(); + + LinearAlgebra::TpetraWrappers::SparseMatrix m(sp); + for (unsigned int i = 0; i < 5; ++i) + for (unsigned int j = 0; j < 5; ++j) + if (((i + 2 * j + 1) % 3 == 0) || (i == j)) + m.set(i, j, 1.); + + LinearAlgebra::TpetraWrappers::SparseMatrix::iterator i = m.begin(); + for (; i != m.end(); ++i) + i->value() -= i->row() * i->column(); + m.compress(VectorOperation::insert); + + for (unsigned int i = 0; i < 5; ++i) + for (unsigned int j = 0; j < 5; ++j) + if (((i + 2 * j + 1) % 3 == 0) || (i == j)) + { + deallog << i << ' ' << j << ' ' << m.el(i, j) << std::endl; + Assert(std::fabs(m.el(i, j) - (1. - i * j)) < 1e-14, + ExcInternalError()); + } + + deallog << "OK" << std::endl; +} + + + +int +main(int argc, char **argv) +{ + initlog(); + + Utilities::MPI::MPI_InitFinalize mpi_initialization( + argc, argv, testing_max_num_threads()); + + try + { + test(); + } + catch (const std::exception &exc) + { + deallog << std::endl + << std::endl + << "----------------------------------------------------" + << std::endl; + deallog << "Exception on processing: " << std::endl + << exc.what() << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + + return 1; + } + catch (...) + { + deallog << std::endl + << std::endl + << "----------------------------------------------------" + << std::endl; + deallog << "Unknown exception!" << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + return 1; + }; +} diff --git a/tests/trilinos_tpetra/trilinos_sparse_matrix_iterator_06.output b/tests/trilinos_tpetra/trilinos_sparse_matrix_iterator_06.output new file mode 100644 index 0000000000..1c38f7a4b7 --- /dev/null +++ b/tests/trilinos_tpetra/trilinos_sparse_matrix_iterator_06.output @@ -0,0 +1,15 @@ + +DEAL::0 0 1.00000 +DEAL::0 1 1.00000 +DEAL::0 4 1.00000 +DEAL::1 1 0 +DEAL::1 2 -1.00000 +DEAL::2 0 1.00000 +DEAL::2 2 -3.00000 +DEAL::2 3 -5.00000 +DEAL::3 1 -2.00000 +DEAL::3 3 -8.00000 +DEAL::3 4 -11.0000 +DEAL::4 2 -7.00000 +DEAL::4 4 -15.0000 +DEAL::OK diff --git a/tests/trilinos_tpetra/trilinos_sparse_matrix_iterator_07.cc b/tests/trilinos_tpetra/trilinos_sparse_matrix_iterator_07.cc new file mode 100644 index 0000000000..5f1d8d6a5a --- /dev/null +++ b/tests/trilinos_tpetra/trilinos_sparse_matrix_iterator_07.cc @@ -0,0 +1,99 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2004 - 2018 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.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + + +// test setting some elements using a non-const matrix iterator and operator*=, +// and reading them back through the matrix itself + +#include +#include + +#include "../tests.h" + + +void +test() +{ + LinearAlgebra::TpetraWrappers::SparsityPattern sp(5, 5, 3); + for (unsigned int i = 0; i < 5; ++i) + for (unsigned int j = 0; j < 5; ++j) + if (((i + 2 * j + 1) % 3 == 0) || (i == j)) + sp.add(i, j); + sp.compress(); + + LinearAlgebra::TpetraWrappers::SparseMatrix m(sp); + for (unsigned int i = 0; i < 5; ++i) + for (unsigned int j = 0; j < 5; ++j) + if (((i + 2 * j + 1) % 3 == 0) || (i == j)) + m.set(i, j, i * j); + + LinearAlgebra::TpetraWrappers::SparseMatrix::iterator i = m.begin(); + for (; i != m.end(); ++i) + i->value() *= 2; + m.compress(VectorOperation::insert); + + for (unsigned int i = 0; i < 5; ++i) + for (unsigned int j = 0; j < 5; ++j) + if (((i + 2 * j + 1) % 3 == 0) || (i == j)) + { + deallog << i << ' ' << j << ' ' << m.el(i, j) << std::endl; + Assert(std::fabs(m.el(i, j) - (2. * i * j)) < 1e-14, + ExcInternalError()); + } + + deallog << "OK" << std::endl; +} + + + +int +main(int argc, char **argv) +{ + initlog(); + + Utilities::MPI::MPI_InitFinalize mpi_initialization( + argc, argv, testing_max_num_threads()); + + try + { + test(); + } + catch (const std::exception &exc) + { + deallog << std::endl + << std::endl + << "----------------------------------------------------" + << std::endl; + deallog << "Exception on processing: " << std::endl + << exc.what() << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + + return 1; + } + catch (...) + { + deallog << std::endl + << std::endl + << "----------------------------------------------------" + << std::endl; + deallog << "Unknown exception!" << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + return 1; + }; +} diff --git a/tests/trilinos_tpetra/trilinos_sparse_matrix_iterator_07.output b/tests/trilinos_tpetra/trilinos_sparse_matrix_iterator_07.output new file mode 100644 index 0000000000..081fbcc5c6 --- /dev/null +++ b/tests/trilinos_tpetra/trilinos_sparse_matrix_iterator_07.output @@ -0,0 +1,15 @@ + +DEAL::0 0 0 +DEAL::0 1 0 +DEAL::0 4 0 +DEAL::1 1 2.00000 +DEAL::1 2 4.00000 +DEAL::2 0 0 +DEAL::2 2 8.00000 +DEAL::2 3 12.0000 +DEAL::3 1 6.00000 +DEAL::3 3 18.0000 +DEAL::3 4 24.0000 +DEAL::4 2 16.0000 +DEAL::4 4 32.0000 +DEAL::OK diff --git a/tests/trilinos_tpetra/trilinos_sparse_matrix_iterator_08.cc b/tests/trilinos_tpetra/trilinos_sparse_matrix_iterator_08.cc new file mode 100644 index 0000000000..9b8f907ca6 --- /dev/null +++ b/tests/trilinos_tpetra/trilinos_sparse_matrix_iterator_08.cc @@ -0,0 +1,99 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2004 - 2018 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.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + + +// test setting some elements using a non-const matrix iterator and operator/=, +// and reading them back through the matrix itself + +#include +#include + +#include "../tests.h" + + +void +test() +{ + LinearAlgebra::TpetraWrappers::SparsityPattern sp(5, 5, 3); + for (unsigned int i = 0; i < 5; ++i) + for (unsigned int j = 0; j < 5; ++j) + if (((i + 2 * j + 1) % 3 == 0) || (i == j)) + sp.add(i, j); + sp.compress(); + + LinearAlgebra::TpetraWrappers::SparseMatrix m(sp); + for (unsigned int i = 0; i < 5; ++i) + for (unsigned int j = 0; j < 5; ++j) + if (((i + 2 * j + 1) % 3 == 0) || (i == j)) + m.set(i, j, i * j); + + LinearAlgebra::TpetraWrappers::SparseMatrix::iterator i = m.begin(); + for (; i != m.end(); ++i) + i->value() /= 2; + m.compress(VectorOperation::insert); + + for (unsigned int i = 0; i < 5; ++i) + for (unsigned int j = 0; j < 5; ++j) + if (((i + 2 * j + 1) % 3 == 0) || (i == j)) + { + deallog << i << ' ' << j << ' ' << m.el(i, j) << std::endl; + Assert(std::fabs(m.el(i, j) - (i * j / 2.)) < 1e-14, + ExcInternalError()); + } + + deallog << "OK" << std::endl; +} + + + +int +main(int argc, char **argv) +{ + initlog(); + + Utilities::MPI::MPI_InitFinalize mpi_initialization( + argc, argv, testing_max_num_threads()); + + try + { + test(); + } + catch (const std::exception &exc) + { + deallog << std::endl + << std::endl + << "----------------------------------------------------" + << std::endl; + deallog << "Exception on processing: " << std::endl + << exc.what() << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + + return 1; + } + catch (...) + { + deallog << std::endl + << std::endl + << "----------------------------------------------------" + << std::endl; + deallog << "Unknown exception!" << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + return 1; + }; +} diff --git a/tests/trilinos_tpetra/trilinos_sparse_matrix_iterator_08.output b/tests/trilinos_tpetra/trilinos_sparse_matrix_iterator_08.output new file mode 100644 index 0000000000..9290c35eba --- /dev/null +++ b/tests/trilinos_tpetra/trilinos_sparse_matrix_iterator_08.output @@ -0,0 +1,15 @@ + +DEAL::0 0 0 +DEAL::0 1 0 +DEAL::0 4 0 +DEAL::1 1 0.500000 +DEAL::1 2 1.00000 +DEAL::2 0 0 +DEAL::2 2 2.00000 +DEAL::2 3 3.00000 +DEAL::3 1 1.50000 +DEAL::3 3 4.50000 +DEAL::3 4 6.00000 +DEAL::4 2 4.00000 +DEAL::4 4 8.00000 +DEAL::OK diff --git a/tests/trilinos_tpetra/trilinos_sparse_matrix_iterator_09.cc b/tests/trilinos_tpetra/trilinos_sparse_matrix_iterator_09.cc new file mode 100644 index 0000000000..cc2b3f998c --- /dev/null +++ b/tests/trilinos_tpetra/trilinos_sparse_matrix_iterator_09.cc @@ -0,0 +1,90 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2004 - 2018 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.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + + + +// this test, extracted from dof_constraints_09 and +// trilinos_sparse_matrix_iterator_09, used to fail with aborts + +#include +#include +#include + +#include "../tests.h" + + +void +test() +{ + // create a sparsity pattern with totally + // empty lines (not even diagonals, since + // not quadratic) + LinearAlgebra::TpetraWrappers::SparsityPattern sparsity(4, 5, 1); + sparsity.add(1, 1); + sparsity.add(3, 1); + sparsity.compress(); + + // attach a sparse matrix to it + LinearAlgebra::TpetraWrappers::SparseMatrix A(sparsity); + + // and loop over the elements of it + for (LinearAlgebra::TpetraWrappers::SparseMatrix::const_iterator k = + A.begin(); + k != A.end(); + ++k) + deallog << k->row() << ' ' << k->column() << ' ' << k->value() << std::endl; +} + + + +int +main(int argc, char **argv) +{ + initlog(); + + Utilities::MPI::MPI_InitFinalize mpi_initialization( + argc, argv, testing_max_num_threads()); + + try + { + test(); + } + catch (const std::exception &exc) + { + deallog << std::endl + << std::endl + << "----------------------------------------------------" + << std::endl; + deallog << "Exception on processing: " << std::endl + << exc.what() << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + + return 1; + } + catch (...) + { + deallog << std::endl + << std::endl + << "----------------------------------------------------" + << std::endl; + deallog << "Unknown exception!" << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + return 1; + }; +} diff --git a/tests/trilinos_tpetra/trilinos_sparse_matrix_iterator_09.output b/tests/trilinos_tpetra/trilinos_sparse_matrix_iterator_09.output new file mode 100644 index 0000000000..62d67f6c93 --- /dev/null +++ b/tests/trilinos_tpetra/trilinos_sparse_matrix_iterator_09.output @@ -0,0 +1,3 @@ + +DEAL::1 1 0 +DEAL::3 1 0 diff --git a/tests/trilinos_tpetra/trilinos_sparse_matrix_iterator_10.cc b/tests/trilinos_tpetra/trilinos_sparse_matrix_iterator_10.cc new file mode 100644 index 0000000000..5ff1074836 --- /dev/null +++ b/tests/trilinos_tpetra/trilinos_sparse_matrix_iterator_10.cc @@ -0,0 +1,92 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2004 - 2018 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.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + + + +// this test, extracted from dof_constraints_09 and +// trilinos_sparse_matrix_iterator_10, used to fail with aborts. this test is +// equivalent to trilinos_sparse_matrix_iterator_09, except that we use the +// postfix operator++ instead of the prefix form + +#include +#include +#include + +#include "../tests.h" + + +void +test() +{ + // create a sparsity pattern with totally + // empty lines (not even diagonals, since + // not quadratic) + LinearAlgebra::TpetraWrappers::SparsityPattern sparsity(4, 5, 1); + sparsity.add(1, 1); + sparsity.add(3, 1); + sparsity.compress(); + + // attach a sparse matrix to it + LinearAlgebra::TpetraWrappers::SparseMatrix A(sparsity); + + // and loop over the elements of it + for (LinearAlgebra::TpetraWrappers::SparseMatrix::const_iterator k = + A.begin(); + k != A.end(); + k++) + deallog << k->row() << ' ' << k->column() << ' ' << k->value() << std::endl; +} + + + +int +main(int argc, char **argv) +{ + initlog(); + + Utilities::MPI::MPI_InitFinalize mpi_initialization( + argc, argv, testing_max_num_threads()); + + try + { + test(); + } + catch (const std::exception &exc) + { + deallog << std::endl + << std::endl + << "----------------------------------------------------" + << std::endl; + deallog << "Exception on processing: " << std::endl + << exc.what() << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + + return 1; + } + catch (...) + { + deallog << std::endl + << std::endl + << "----------------------------------------------------" + << std::endl; + deallog << "Unknown exception!" << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + return 1; + }; +} diff --git a/tests/trilinos_tpetra/trilinos_sparse_matrix_iterator_10.output b/tests/trilinos_tpetra/trilinos_sparse_matrix_iterator_10.output new file mode 100644 index 0000000000..62d67f6c93 --- /dev/null +++ b/tests/trilinos_tpetra/trilinos_sparse_matrix_iterator_10.output @@ -0,0 +1,3 @@ + +DEAL::1 1 0 +DEAL::3 1 0 diff --git a/tests/trilinos_tpetra/trilinos_sparse_matrix_iterator_11.cc b/tests/trilinos_tpetra/trilinos_sparse_matrix_iterator_11.cc new file mode 100644 index 0000000000..b29ba58ad0 --- /dev/null +++ b/tests/trilinos_tpetra/trilinos_sparse_matrix_iterator_11.cc @@ -0,0 +1,100 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2004 - 2021 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.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + + + +// certain comparisons between sparse matrix iterators didn't compile. + +#include +#include +#include + +#include "../tests.h" + + +void +test() +{ + // create a sparsity pattern with totally + // empty lines (not even diagonals, since + // not quadratic) + LinearAlgebra::TpetraWrappers::SparsityPattern sparsity(4, 5, 1); + sparsity.add(1, 1); + sparsity.add(3, 1); + sparsity.compress(); + + // attach a sparse matrix to it + LinearAlgebra::TpetraWrappers::SparseMatrix A(sparsity); + + LinearAlgebra::TpetraWrappers::SparseMatrix::iterator k = A.begin(), + j = std::next( + A.begin()); + + AssertThrow(k < j, ExcInternalError()); + AssertThrow(j > k, ExcInternalError()); + + AssertThrow(!(j < k), ExcInternalError()); + AssertThrow(!(k > j), ExcInternalError()); + + AssertThrow(k != j, ExcInternalError()); + AssertThrow(!(k == j), ExcInternalError()); + + AssertThrow(k == k, ExcInternalError()); + AssertThrow(!(k != k), ExcInternalError()); + + deallog << "OK" << std::endl; +} + + + +int +main(int argc, char **argv) +{ + initlog(); + + Utilities::MPI::MPI_InitFinalize mpi_initialization( + argc, argv, testing_max_num_threads()); + + try + { + test(); + } + catch (const std::exception &exc) + { + deallog << std::endl + << std::endl + << "----------------------------------------------------" + << std::endl; + deallog << "Exception on processing: " << std::endl + << exc.what() << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + + return 1; + } + catch (...) + { + deallog << std::endl + << std::endl + << "----------------------------------------------------" + << std::endl; + deallog << "Unknown exception!" << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + return 1; + }; +} diff --git a/tests/trilinos_tpetra/trilinos_sparse_matrix_iterator_11.output b/tests/trilinos_tpetra/trilinos_sparse_matrix_iterator_11.output new file mode 100644 index 0000000000..0fd8fc12f0 --- /dev/null +++ b/tests/trilinos_tpetra/trilinos_sparse_matrix_iterator_11.output @@ -0,0 +1,2 @@ + +DEAL::OK diff --git a/tests/trilinos_tpetra/trilinos_sparse_matrix_iterator_12.cc b/tests/trilinos_tpetra/trilinos_sparse_matrix_iterator_12.cc new file mode 100644 index 0000000000..a5e1a997c0 --- /dev/null +++ b/tests/trilinos_tpetra/trilinos_sparse_matrix_iterator_12.cc @@ -0,0 +1,100 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2004 - 2021 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.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + + + +// like trilinos_sparse_matrix_iterator_11, but for const_iterators + +#include +#include +#include + +#include "../tests.h" + + +void +test() +{ + // create a sparsity pattern with totally + // empty lines (not even diagonals, since + // not quadratic) + LinearAlgebra::TpetraWrappers::SparsityPattern sparsity(4, 5, 1); + sparsity.add(1, 1); + sparsity.add(3, 1); + sparsity.compress(); + + // attach a sparse matrix to it + LinearAlgebra::TpetraWrappers::SparseMatrix A(sparsity); + + LinearAlgebra::TpetraWrappers::SparseMatrix::const_iterator + k = A.begin(), + j = std::next(A.begin()); + + AssertThrow(k < j, ExcInternalError()); + AssertThrow(j > k, ExcInternalError()); + + AssertThrow(!(j < k), ExcInternalError()); + AssertThrow(!(k > j), ExcInternalError()); + + AssertThrow(k != j, ExcInternalError()); + AssertThrow(!(k == j), ExcInternalError()); + + AssertThrow(k == k, ExcInternalError()); + AssertThrow(!(k != k), ExcInternalError()); + + deallog << "OK" << std::endl; +} + + + +int +main(int argc, char **argv) +{ + initlog(); + + Utilities::MPI::MPI_InitFinalize mpi_initialization( + argc, argv, testing_max_num_threads()); + + try + { + test(); + } + catch (const std::exception &exc) + { + deallog << std::endl + << std::endl + << "----------------------------------------------------" + << std::endl; + deallog << "Exception on processing: " << std::endl + << exc.what() << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + + return 1; + } + catch (...) + { + deallog << std::endl + << std::endl + << "----------------------------------------------------" + << std::endl; + deallog << "Unknown exception!" << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + return 1; + }; +} diff --git a/tests/trilinos_tpetra/trilinos_sparse_matrix_iterator_12.output b/tests/trilinos_tpetra/trilinos_sparse_matrix_iterator_12.output new file mode 100644 index 0000000000..0fd8fc12f0 --- /dev/null +++ b/tests/trilinos_tpetra/trilinos_sparse_matrix_iterator_12.output @@ -0,0 +1,2 @@ + +DEAL::OK -- 2.39.5