15#ifndef dealii_chunk_sparse_matrix_h
16#define dealii_chunk_sparse_matrix_h
36template <
typename number>
38template <
typename number>
54 template <
typename number,
bool Constness>
67 template <
typename number,
bool Constness>
99 template <
typename number>
149 template <
typename,
bool>
160 template <
typename number>
198 operator number()
const;
280 template <
typename,
bool>
296 template <
typename number,
bool Constness>
415 template <
typename number,
bool Constness>
421 Iterator<number, Constness>::value_type;
423 Iterator<number, Constness>::difference_type;
447template <
typename number>
706 template <
typename number2>
758 template <
typename somenumber>
779 template <
typename ForwardIterator>
788 template <
typename somenumber>
803 template <
typename somenumber>
870 number *values)
const;
891 template <
class OutVector,
class InVector>
910 template <
class OutVector,
class InVector>
928 template <
class OutVector,
class InVector>
947 template <
class OutVector,
class InVector>
966 template <
typename somenumber>
973 template <
typename somenumber>
984 template <
typename somenumber>
1033 template <
typename somenumber>
1037 const number omega = 1.)
const;
1042 template <
typename somenumber>
1046 const number
om = 1.)
const;
1051 template <
typename somenumber>
1055 const number
om = 1.)
const;
1060 template <
typename somenumber>
1064 const number
om = 1.)
const;
1071 template <
typename somenumber>
1079 template <
typename somenumber>
1087 template <
typename somenumber>
1101 template <
typename somenumber>
1104 const std::vector<size_type> &permutation,
1105 const std::vector<size_type> &inverse_permutation,
1106 const number
om = 1.)
const;
1118 template <
typename somenumber>
1121 const std::vector<size_type> &permutation,
1122 const std::vector<size_type> &inverse_permutation,
1123 const number
om = 1.)
const;
1129 template <
typename somenumber>
1133 const number
om = 1.)
const;
1139 template <
typename somenumber>
1143 const number
om = 1.)
const;
1149 template <
typename somenumber>
1153 const number
om = 1.)
const;
1310 const unsigned int precision = 3,
1311 const bool scientific =
true,
1312 const unsigned int width = 0,
1315 const char *separator =
" ")
const;
1368 <<
"You are trying to access the matrix entry with index <"
1370 <<
">, but this entry does not exist in the sparsity pattern "
1373 "The most common cause for this problem is that you used "
1374 "a method to build the sparsity pattern that did not "
1375 "(completely) take into account all of the entries you "
1376 "will later try to write into. An example would be "
1377 "building a sparsity pattern that does not include "
1378 "the entries you will write into due to constraints "
1379 "on degrees of freedom such as hanging nodes or periodic "
1380 "boundary conditions. In such cases, building the "
1381 "sparsity pattern will succeed, but you will get errors "
1382 "such as the current one at one point or other when "
1383 "trying to write into the entries of the matrix.");
1394 <<
"The iterators denote a range of " <<
arg1
1395 <<
" elements, but the given number of rows was " <<
arg2);
1400 "You are attempting an operation on two vectors that "
1401 "are the same object, but the operation requires that the "
1402 "two objects are in fact different.");
1419 std::unique_ptr<number[]>
val;
1436 template <
typename somenumber>
1440 template <
typename,
bool>
1442 template <
typename,
bool>
1453template <
typename number>
1462template <
typename number>
1472template <
typename number>
1482template <
typename number>
1487 const size_type chunk_size = cols->get_chunk_size();
1489 cols->sparsity_pattern(i / chunk_size,
j / chunk_size);
1496 (i % chunk_size) * chunk_size + (
j % chunk_size));
1501template <
typename number>
1514 ExcInvalidIndex(i,
j));
1522template <
typename number>
1536 ExcInvalidIndex(i,
j));
1544template <
typename number>
1545template <
typename number2>
1555 for (size_type col = 0; col < n_cols; ++col)
1561template <
typename number>
1568 const size_type chunk_size = cols->get_chunk_size();
1577 cols->sparsity_pattern.n_nonzero_elements() * chunk_size * chunk_size;
1586template <
typename number>
1596 const size_type chunk_size = cols->get_chunk_size();
1605 cols->sparsity_pattern.n_nonzero_elements() * chunk_size * chunk_size;
1615template <
typename number>
1622 ExcInvalidIndex(i,
j));
1623 return val[compute_location(i,
j)];
1628template <
typename number>
1643template <
typename number>
1653 const size_type chunk_size = cols->get_chunk_size();
1654 return val[cols->sparsity_pattern.rowstart[i / chunk_size] * chunk_size *
1656 (i % chunk_size) * chunk_size + (i % chunk_size)];
1661template <
typename number>
1662template <
typename ForwardIterator>
1668 ExcIteratorRange(std::distance(begin, end), m()));
1673 typename std::iterator_traits<ForwardIterator>::value_type::const_iterator;
1680 set(row,
j->first,
j->second);
1691 template <
typename number>
1693 const unsigned int row)
1701 template <
typename number>
1702 inline Accessor<number, true>::Accessor(
const MatrixType *matrix)
1709 template <
typename number>
1710 inline Accessor<number, true>::Accessor(
1713 ,
matrix(&a.get_matrix())
1718 template <
typename number>
1720 Accessor<number, true>::value()
const
1722 const unsigned int chunk_size =
1723 matrix->get_sparsity_pattern().get_chunk_size();
1724 return matrix->val[reduced_index() * chunk_size * chunk_size +
1725 chunk_row * chunk_size + chunk_col];
1730 template <
typename number>
1731 inline const typename Accessor<number, true>::MatrixType &
1732 Accessor<number, true>::get_matrix()
const
1739 template <
typename number>
1740 inline Accessor<number, false>::Reference::Reference(
const Accessor *accessor,
1742 : accessor(accessor)
1746 template <
typename number>
1747 inline Accessor<number, false>::Reference::operator number()
const
1749 const unsigned int chunk_size =
1750 accessor->matrix->get_sparsity_pattern().get_chunk_size();
1751 return accessor->matrix
1752 ->val[accessor->reduced_index() * chunk_size * chunk_size +
1753 accessor->chunk_row * chunk_size + accessor->chunk_col];
1758 template <
typename number>
1759 inline const typename Accessor<number, false>::Reference &
1760 Accessor<number, false>::Reference::operator=(
const number n)
const
1762 const unsigned int chunk_size =
1763 accessor->matrix->get_sparsity_pattern().get_chunk_size();
1765 ->val[accessor->reduced_index() * chunk_size * chunk_size +
1766 accessor->chunk_row * chunk_size + accessor->chunk_col] = n;
1772 template <
typename number>
1773 inline const typename Accessor<number, false>::Reference &
1774 Accessor<number, false>::Reference::operator+=(
const number n)
const
1776 const unsigned int chunk_size =
1777 accessor->matrix->get_sparsity_pattern().get_chunk_size();
1779 ->val[accessor->reduced_index() * chunk_size * chunk_size +
1780 accessor->chunk_row * chunk_size + accessor->chunk_col] += n;
1786 template <
typename number>
1787 inline const typename Accessor<number, false>::Reference &
1788 Accessor<number, false>::Reference::operator-=(
const number n)
const
1790 const unsigned int chunk_size =
1791 accessor->matrix->get_sparsity_pattern().get_chunk_size();
1793 ->val[accessor->reduced_index() * chunk_size * chunk_size +
1794 accessor->chunk_row * chunk_size + accessor->chunk_col] -= n;
1800 template <
typename number>
1801 inline const typename Accessor<number, false>::Reference &
1802 Accessor<number, false>::Reference::operator*=(
const number n)
const
1804 const unsigned int chunk_size =
1805 accessor->matrix->get_sparsity_pattern().get_chunk_size();
1807 ->val[accessor->reduced_index() * chunk_size * chunk_size +
1808 accessor->chunk_row * chunk_size + accessor->chunk_col] *= n;
1814 template <
typename number>
1815 inline const typename Accessor<number, false>::Reference &
1816 Accessor<number, false>::Reference::operator/=(
const number n)
const
1818 const unsigned int chunk_size =
1819 accessor->matrix->get_sparsity_pattern().get_chunk_size();
1821 ->val[accessor->reduced_index() * chunk_size * chunk_size +
1822 accessor->chunk_row * chunk_size + accessor->chunk_col] /= n;
1828 template <
typename number>
1829 inline Accessor<number, false>::Accessor(MatrixType *matrix,
1830 const unsigned int row)
1838 template <
typename number>
1839 inline Accessor<number, false>::Accessor(MatrixType *matrix)
1846 template <
typename number>
1847 inline typename Accessor<number, false>::Reference
1848 Accessor<number, false>::value()
const
1850 return Reference(
this,
true);
1855 template <
typename number>
1856 inline typename Accessor<number, false>::MatrixType &
1857 Accessor<number, false>::get_matrix()
const
1864 template <
typename number,
bool Constness>
1865 inline Iterator<number, Constness>::Iterator(MatrixType *matrix,
1866 const unsigned int row)
1872 template <
typename number,
bool Constness>
1873 inline Iterator<number, Constness>::Iterator(MatrixType *matrix)
1879 template <
typename number,
bool Constness>
1880 inline Iterator<number, Constness>::Iterator(
1887 template <
typename number,
bool Constness>
1889 Iterator<number, Constness>::operator++()
1896 template <
typename number,
bool Constness>
1898 Iterator<number, Constness>::operator++(
int)
1900 const Iterator iter = *
this;
1906 template <
typename number,
bool Constness>
1908 Iterator<number, Constness>::operator*()
const
1914 template <
typename number,
bool Constness>
1916 Iterator<number, Constness>::operator->()
const
1922 template <
typename number,
bool Constness>
1924 Iterator<number, Constness>::operator==(
const Iterator &
other)
const
1926 return (accessor ==
other.accessor);
1930 template <
typename number,
bool Constness>
1932 Iterator<number, Constness>::operator!=(
const Iterator &
other)
const
1934 return !(*
this ==
other);
1938 template <
typename number,
bool Constness>
1940 Iterator<number, Constness>::operator<(
const Iterator &
other)
const
1942 Assert(&accessor.get_matrix() == &
other.accessor.get_matrix(),
1945 return (accessor <
other.accessor);
1949 template <
typename number,
bool Constness>
1951 Iterator<number, Constness>::operator>(
const Iterator &
other)
const
1953 return (
other < *
this);
1957 template <
typename number,
bool Constness>
1959 Iterator<number, Constness>::operator-(
const Iterator &
other)
const
1961 Assert(&accessor.get_matrix() == &
other.accessor.get_matrix(),
1968 Iterator
copy = *
this;
1969 while (copy !=
other)
1978 while (copy != *
this)
1989 template <
typename number,
bool Constness>
1991 Iterator<number, Constness>::operator+(
const unsigned int n)
const
1994 for (
unsigned int i = 0; i < n; ++i)
2004template <
typename number>
2008 return const_iterator(
this, 0);
2012template <
typename number>
2016 return const_iterator(
this);
2020template <
typename number>
2024 return iterator(
this, 0);
2028template <
typename number>
2032 return iterator(
this);
2036template <
typename number>
2041 return const_iterator(
this,
r);
2046template <
typename number>
2051 return const_iterator(
this,
r + 1);
2056template <
typename number>
2061 return iterator(
this,
r);
2066template <
typename number>
2071 return iterator(
this,
r + 1);
const Reference & operator=(const number n) const
const Reference & operator/=(const number n) const
const Reference & operator-=(const number n) const
const Accessor * accessor
const Reference & operator*=(const number n) const
const Reference & operator+=(const number n) const
Reference(const Accessor *accessor, const bool dummy)
MatrixType & get_matrix() const
Accessor(MatrixType *matrix)
Accessor(MatrixType *matrix, const unsigned int row)
Accessor(MatrixType *matrix, const unsigned int row)
const MatrixType & get_matrix() const
Accessor(const ChunkSparseMatrixIterators::Accessor< number, false > &a)
Accessor(MatrixType *matrix)
const ChunkSparseMatrix< number > & get_matrix() const
bool operator>(const Iterator &) const
Iterator(const ChunkSparseMatrixIterators::Iterator< number, false > &i)
const Accessor< number, Constness > & value_type
bool operator<(const Iterator &) const
Iterator(MatrixType *matrix, const unsigned int row)
const Accessor< number, Constness > * operator->() const
bool operator==(const Iterator &) const
Accessor< number, Constness > accessor
const Accessor< number, Constness > & operator*() const
int operator-(const Iterator &p) const
typename Accessor< number, Constness >::MatrixType MatrixType
bool operator!=(const Iterator &) const
Iterator operator+(const unsigned int n) const
Iterator(MatrixType *matrix)
somenumber matrix_scalar_product(const Vector< somenumber > &u, const Vector< somenumber > &v) const
void print_formatted(std::ostream &out, const unsigned int precision=3, const bool scientific=true, const unsigned int width=0, const char *zero_string=" ", const double denominator=1., const char *separator=" ") const
real_type linfty_norm() const
size_type compute_location(const size_type i, const size_type j) const
ChunkSparseMatrix(const ChunkSparsityPattern &sparsity)
void precondition_Jacobi(Vector< somenumber > &dst, const Vector< somenumber > &src, const number omega=1.) const
ChunkSparseMatrix< number > & operator=(const IdentityMatrix &id)
void Tvmult(OutVector &dst, const InVector &src) const
const_iterator begin() const
void add(const size_type i, const size_type j, const number value)
void set(const size_type i, const size_type j, const number value)
void precondition_SSOR(Vector< somenumber > &dst, const Vector< somenumber > &src, const number om=1.) const
const_iterator begin(const unsigned int r) const
void print(std::ostream &out) const
ChunkSparseMatrix & operator/=(const number factor)
ChunkSparseMatrix & operator=(const double d)
virtual void reinit(const ChunkSparsityPattern &sparsity)
void extract_row_copy(const size_type row, const size_type array_length, size_type &row_length, size_type *column_indices, number *values) const
iterator end(const unsigned int r)
void print_pattern(std::ostream &out, const double threshold=0.) const
void TPSOR(Vector< somenumber > &v, const std::vector< size_type > &permutation, const std::vector< size_type > &inverse_permutation, const number om=1.) const
void copy_from(const ForwardIterator begin, const ForwardIterator end)
real_type l1_norm() const
const_iterator end() const
void vmult_add(OutVector &dst, const InVector &src) const
void precondition_SOR(Vector< somenumber > &dst, const Vector< somenumber > &src, const number om=1.) const
ObserverPointer< const ChunkSparsityPattern, ChunkSparseMatrix< number > > cols
somenumber matrix_norm_square(const Vector< somenumber > &v) const
virtual ~ChunkSparseMatrix() override
std::unique_ptr< number[]> val
void block_write(std::ostream &out) const
void SSOR(Vector< somenumber > &v, const number omega=1.) const
void precondition_TSOR(Vector< somenumber > &dst, const Vector< somenumber > &src, const number om=1.) const
size_type n_actually_nonzero_elements() const
ChunkSparseMatrix< number > & operator=(const ChunkSparseMatrix< number > &)
ChunkSparseMatrix(const ChunkSparsityPattern &sparsity, const IdentityMatrix &id)
void vmult(OutVector &dst, const InVector &src) const
void TSOR(Vector< somenumber > &v, const number om=1.) const
void add(const number factor, const ChunkSparseMatrix< somenumber > &matrix)
size_type n_nonzero_elements() const
void copy_from(const FullMatrix< somenumber > &matrix)
ChunkSparseMatrix & operator*=(const number factor)
void block_read(std::istream &in)
number operator()(const size_type i, const size_type j) const
void Tvmult_add(OutVector &dst, const InVector &src) const
number el(const size_type i, const size_type j) const
ChunkSparseMatrix< number > & copy_from(const ChunkSparseMatrix< somenumber > &source)
void TSOR_step(Vector< somenumber > &v, const Vector< somenumber > &b, const number om=1.) const
typename numbers::NumberTraits< number >::real_type real_type
ChunkSparseMatrix(const ChunkSparseMatrix &)
number diag_element(const size_type i) const
void SOR(Vector< somenumber > &v, const number om=1.) const
iterator begin(const unsigned int r)
void SSOR_step(Vector< somenumber > &v, const Vector< somenumber > &b, const number om=1.) const
std::size_t memory_consumption() const
void add(const size_type row, const size_type n_cols, const size_type *col_indices, const number2 *values, const bool elide_zero_values=true, const bool col_indices_are_sorted=false)
real_type frobenius_norm() const
somenumber residual(Vector< somenumber > &dst, const Vector< somenumber > &x, const Vector< somenumber > &b) const
void PSOR(Vector< somenumber > &v, const std::vector< size_type > &permutation, const std::vector< size_type > &inverse_permutation, const number om=1.) const
void SOR_step(Vector< somenumber > &v, const Vector< somenumber > &b, const number om=1.) const
const_iterator end(const unsigned int r) const
const ChunkSparsityPattern & get_sparsity_pattern() const
Accessor(const ChunkSparsityPattern *matrix, const size_type row)
types::global_dof_index size_type
static const size_type invalid_entry
static constexpr size_type invalid_entry
std::conditional_t< rank_==1, std::array< Number, dim >, std::array< Tensor< rank_ - 1, dim, Number >, dim > > values
#define DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_CLOSE
#define DeclException0(Exception0)
static ::ExceptionBase & ExcDifferentChunkSparsityPatterns()
static ::ExceptionBase & ExcSourceEqualsDestination()
#define Assert(cond, exc)
#define AssertIsFinite(number)
#define DeclException2(Exception2, type1, type2, outsequence)
static ::ExceptionBase & ExcNeedsSparsityPattern()
#define AssertIndexRange(index, range)
#define DeclExceptionMsg(Exception, defaulttext)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcDivideByZero()
static ::ExceptionBase & ExcNotInitialized()
static ::ExceptionBase & ExcNotQuadratic()
static ::ExceptionBase & ExcInvalidIndex(int arg1, int arg2)
static ::ExceptionBase & ExcIteratorRange(int arg1, int arg2)
#define AssertThrow(cond, exc)
types::global_dof_index size_type
@ matrix
Contents is actually a matrix.
VectorType::value_type * end(VectorType &V)
void copy(const T *begin, const T *end, U *dest)
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)
unsigned int global_dof_index
static const bool zero_addition_can_be_elided
typename ::ChunkSparseMatrixIterators::Iterator< number, Constness >::difference_type difference_type
typename ::ChunkSparseMatrixIterators::Iterator< number, Constness >::value_type value_type