15#ifndef dealii_sparse_matrix_h
16#define dealii_sparse_matrix_h
37template <
typename number>
39template <
typename number>
41template <
typename Matrix>
43template <
typename number>
45# ifdef DEAL_II_WITH_MPI
50 template <
typename Number>
57# ifdef DEAL_II_WITH_TRILINOS
82 template <
typename number,
bool Constness>
95 template <
typename number,
bool Constness>
127 template <
typename number>
177 template <
typename,
bool>
188 template <
typename number>
228 operator number()
const;
310 template <
typename,
bool>
345 template <
typename number,
bool Constness>
470 template <
typename number,
bool Constness>
476 typename ::SparseMatrixIterators::Iterator<number,
477 Constness>::value_type;
479 Iterator<number, Constness>::difference_type;
518template <
typename number>
701 template <
typename number2>
815 template <
typename number2>
817 set(
const std::vector<size_type> &indices,
819 const bool elide_zero_values =
false);
826 template <
typename number2>
828 set(
const std::vector<size_type> &row_indices,
829 const std::vector<size_type> &col_indices,
831 const bool elide_zero_values =
false);
843 template <
typename number2>
846 const std::vector<size_type> &col_indices,
847 const std::vector<number2> &values,
848 const bool elide_zero_values =
false);
859 template <
typename number2>
864 const number2 *values,
865 const bool elide_zero_values =
false);
889 template <
typename number2>
891 add(
const std::vector<size_type> &indices,
893 const bool elide_zero_values =
true);
900 template <
typename number2>
902 add(
const std::vector<size_type> &row_indices,
903 const std::vector<size_type> &col_indices,
905 const bool elide_zero_values =
true);
916 template <
typename number2>
919 const std::vector<size_type> &col_indices,
920 const std::vector<number2> &values,
921 const bool elide_zero_values =
true);
932 template <
typename number2>
937 const number2 *values,
938 const bool elide_zero_values =
true,
939 const bool col_indices_are_sorted =
false);
984 template <
typename somenumber>
1004 template <
typename ForwardIterator>
1017 template <
typename somenumber>
1021#ifdef DEAL_II_WITH_TRILINOS
1046 template <
typename somenumber>
1133 template <
class OutVector,
class InVector>
1135 vmult(OutVector &dst,
const InVector &src)
const;
1152 template <
class OutVector,
class InVector>
1154 Tvmult(OutVector &dst,
const InVector &src)
const;
1172 template <
class OutVector,
class InVector>
1191 template <
class OutVector,
class InVector>
1212 template <
typename somenumber>
1221 template <
typename somenumber>
1235 template <
typename somenumber>
1276 template <
typename numberB,
typename numberC>
1281 const bool rebuild_sparsity_pattern =
true)
const;
1307 template <
typename numberB,
typename numberC>
1312 const bool rebuild_sparsity_pattern =
true)
const;
1357 template <
typename somenumber>
1361 const number omega = 1.)
const;
1369 template <
typename somenumber>
1373 const number omega = 1.,
1374 const std::vector<std::size_t> &pos_right_of_diagonal =
1375 std::vector<std::size_t>())
const;
1380 template <
typename somenumber>
1384 const number omega = 1.)
const;
1389 template <
typename somenumber>
1393 const number omega = 1.)
const;
1400 template <
typename somenumber>
1408 template <
typename somenumber>
1416 template <
typename somenumber>
1430 template <
typename somenumber>
1433 const std::vector<size_type> &permutation,
1434 const std::vector<size_type> &inverse_permutation,
1435 const number omega = 1.)
const;
1447 template <
typename somenumber>
1450 const std::vector<size_type> &permutation,
1451 const std::vector<size_type> &inverse_permutation,
1452 const number omega = 1.)
const;
1459 template <
typename somenumber>
1463 const number omega = 1.)
const;
1469 template <
typename somenumber>
1473 const number omega = 1.)
const;
1479 template <
typename somenumber>
1483 const number omega = 1.)
const;
1489 template <
typename somenumber>
1493 const number omega = 1.)
const;
1579 template <
typename StreamType>
1582 const bool across =
false,
1583 const bool diagonal_first =
true)
const;
1609 const unsigned int precision = 3,
1610 const bool scientific =
true,
1611 const unsigned int width = 0,
1612 const char *zero_string =
" ",
1613 const double denominator = 1.,
1614 const char *separator =
" ")
const;
1634 const unsigned int precision = 9)
const;
1679 <<
"You are trying to access the matrix entry with index <"
1680 << arg1 <<
',' << arg2
1681 <<
">, but this entry does not exist in the sparsity pattern "
1684 "The most common cause for this problem is that you used "
1685 "a method to build the sparsity pattern that did not "
1686 "(completely) take into account all of the entries you "
1687 "will later try to write into. An example would be "
1688 "building a sparsity pattern that does not include "
1689 "the entries you will write into due to constraints "
1690 "on degrees of freedom such as hanging nodes or periodic "
1691 "boundary conditions. In such cases, building the "
1692 "sparsity pattern will succeed, but you will get errors "
1693 "such as the current one at one point or other when "
1694 "trying to write into the entries of the matrix.");
1699 "When copying one sparse matrix into another, "
1700 "or when adding one sparse matrix to another, "
1701 "both matrices need to refer to the same "
1702 "sparsity pattern.");
1709 <<
"The iterators denote a range of " << arg1
1710 <<
" elements, but the given number of rows was " << arg2);
1715 "You are attempting an operation on two vectors that "
1716 "are the same object, but the operation requires that the "
1717 "two objects are in fact different.");
1756 std::unique_ptr<number[]>
val;
1767 template <
typename somenumber>
1769 template <
typename somenumber>
1779 template <
typename,
bool>
1781 template <
typename,
bool>
1784#ifdef DEAL_II_WITH_MPI
1786 template <
typename Number>
1799template <
typename number>
1800template <
typename number2>
1809template <
typename number>
1819template <
typename number>
1829template <
typename number>
1840template <
typename number>
1855 ExcInvalidIndex(i, j));
1864template <
typename number>
1865template <
typename number2>
1869 const bool elide_zero_values)
1875 for (
size_type i = 0; i < indices.size(); ++i)
1885template <
typename number>
1886template <
typename number2>
1889 const std::vector<size_type> &col_indices,
1891 const bool elide_zero_values)
1898 for (
size_type i = 0; i < row_indices.size(); ++i)
1908template <
typename number>
1909template <
typename number2>
1912 const std::vector<size_type> &col_indices,
1913 const std::vector<number2> &values,
1914 const bool elide_zero_values)
1928template <
typename number>
1936 if (
value == number())
1946 ExcInvalidIndex(i, j));
1955template <
typename number>
1956template <
typename number2>
1960 const bool elide_zero_values)
1966 for (
size_type i = 0; i < indices.size(); ++i)
1976template <
typename number>
1977template <
typename number2>
1980 const std::vector<size_type> &col_indices,
1982 const bool elide_zero_values)
1989 for (
size_type i = 0; i < row_indices.size(); ++i)
1999template <
typename number>
2000template <
typename number2>
2003 const std::vector<size_type> &col_indices,
2004 const std::vector<number2> &values,
2005 const bool elide_zero_values)
2015 std::is_sorted(col_indices.begin(), col_indices.end()));
2020template <
typename number>
2027 number *val_ptr = val.get();
2028 const number *
const end_ptr = val.get() + cols->n_nonzero_elements();
2030 while (val_ptr != end_ptr)
2031 *val_ptr++ *= factor;
2038template <
typename number>
2046 const number factor_inv = number(1.) / factor;
2048 number *val_ptr = val.get();
2049 const number *
const end_ptr = val.get() + cols->n_nonzero_elements();
2051 while (val_ptr != end_ptr)
2052 *val_ptr++ *= factor_inv;
2059template <
typename number>
2060inline const number &
2065 ExcInvalidIndex(i, j));
2066 return val[cols->operator()(i, j)];
2071template <
typename number>
2077 ExcInvalidIndex(i, j));
2078 return val[cols->operator()(i, j)];
2083template <
typename number>
2098template <
typename number>
2108 return val[cols->rowstart[i]];
2113template <
typename number>
2123 return val[cols->rowstart[i]];
2128template <
typename number>
2129template <
typename ForwardIterator>
2132 const ForwardIterator end)
2135 ExcIteratorRange(std::distance(begin, end), m()));
2139 using inner_iterator =
2140 typename std::iterator_traits<ForwardIterator>::value_type::const_iterator;
2142 for (ForwardIterator i = begin; i !=
end; ++i, ++
row)
2144 const inner_iterator end_of_row = i->end();
2145 for (inner_iterator j = i->begin(); j != end_of_row; ++j)
2147 set(
row, j->first, j->second);
2157 template <
typename number>
2159 const std::size_t index_within_matrix)
2161 index_within_matrix)
2167 template <
typename number>
2168 inline Accessor<number, true>::Accessor(
const MatrixType *matrix)
2175 template <
typename number>
2176 inline Accessor<number, true>::Accessor(
2179 ,
matrix(&a.get_matrix())
2184 template <
typename number>
2186 Accessor<number, true>::value()
const
2189 return matrix->val[linear_index];
2194 template <
typename number>
2195 inline const typename Accessor<number, true>::MatrixType &
2196 Accessor<number, true>::get_matrix()
const
2203 template <
typename number>
2204 inline Accessor<number, false>::Reference::Reference(
const Accessor *accessor,
2206 : accessor(accessor)
2210 template <
typename number>
2211 inline Accessor<number, false>::Reference::operator number()
const
2214 accessor->matrix->n_nonzero_elements());
2215 return accessor->matrix->val[accessor->linear_index];
2220 template <
typename number>
2221 inline const typename Accessor<number, false>::Reference &
2222 Accessor<number, false>::Reference::operator=(
const number n)
const
2225 accessor->matrix->n_nonzero_elements());
2226 accessor->matrix->val[accessor->linear_index] = n;
2232 template <
typename number>
2233 inline const typename Accessor<number, false>::Reference &
2234 Accessor<number, false>::Reference::operator+=(
const number n)
const
2237 accessor->matrix->n_nonzero_elements());
2238 accessor->matrix->val[accessor->linear_index] += n;
2244 template <
typename number>
2245 inline const typename Accessor<number, false>::Reference &
2246 Accessor<number, false>::Reference::operator-=(
const number n)
const
2249 accessor->matrix->n_nonzero_elements());
2250 accessor->matrix->val[accessor->linear_index] -= n;
2256 template <
typename number>
2257 inline const typename Accessor<number, false>::Reference &
2258 Accessor<number, false>::Reference::operator*=(
const number n)
const
2261 accessor->matrix->n_nonzero_elements());
2262 accessor->matrix->val[accessor->linear_index] *= n;
2268 template <
typename number>
2269 inline const typename Accessor<number, false>::Reference &
2270 Accessor<number, false>::Reference::operator/=(
const number n)
const
2273 accessor->matrix->n_nonzero_elements());
2274 accessor->matrix->val[accessor->linear_index] /= n;
2280 template <
typename number>
2281 inline Accessor<number, false>::Accessor(MatrixType *matrix,
2282 const std::size_t index)
2289 template <
typename number>
2290 inline Accessor<number, false>::Accessor(MatrixType *matrix)
2297 template <
typename number>
2298 inline typename Accessor<number, false>::Reference
2299 Accessor<number, false>::value()
const
2301 return Reference(
this,
true);
2306 template <
typename number>
2307 inline typename Accessor<number, false>::MatrixType &
2308 Accessor<number, false>::get_matrix()
const
2315 template <
typename number,
bool Constness>
2316 inline Iterator<number, Constness>::Iterator(MatrixType *matrix,
2317 const std::size_t index)
2323 template <
typename number,
bool Constness>
2324 inline Iterator<number, Constness>::Iterator(MatrixType *matrix)
2330 template <
typename number,
bool Constness>
2331 inline Iterator<number, Constness>::Iterator(
2338 template <
typename number,
bool Constness>
2339 inline const Iterator<number, Constness> &
2340 Iterator<number, Constness>::operator=(
2349 template <
typename number,
bool Constness>
2350 inline Iterator<number, Constness> &
2351 Iterator<number, Constness>::operator++()
2358 template <
typename number,
bool Constness>
2359 inline Iterator<number, Constness>
2360 Iterator<number, Constness>::operator++(
int)
2362 const Iterator iter = *
this;
2368 template <
typename number,
bool Constness>
2369 inline const Accessor<number, Constness> &
2370 Iterator<number, Constness>::operator*()
const
2376 template <
typename number,
bool Constness>
2377 inline const Accessor<number, Constness> *
2378 Iterator<number, Constness>::operator->()
const
2384 template <
typename number,
bool Constness>
2386 Iterator<number, Constness>::operator==(
const Iterator &other)
const
2388 return (accessor == other.accessor);
2392 template <
typename number,
bool Constness>
2394 Iterator<number, Constness>::operator!=(
const Iterator &other)
const
2396 return !(*
this == other);
2400 template <
typename number,
bool Constness>
2402 Iterator<number, Constness>::operator<(
const Iterator &other)
const
2404 Assert(&accessor.get_matrix() == &other.accessor.get_matrix(),
2407 return (accessor < other.accessor);
2411 template <
typename number,
bool Constness>
2413 Iterator<number, Constness>::operator>(
const Iterator &other)
const
2415 return (other < *
this);
2419 template <
typename number,
bool Constness>
2421 Iterator<number, Constness>::operator-(
const Iterator &other)
const
2423 Assert(&accessor.get_matrix() == &other.accessor.get_matrix(),
2426 return (*this)->linear_index - other->linear_index;
2431 template <
typename number,
bool Constness>
2432 inline Iterator<number, Constness>
2433 Iterator<number, Constness>::operator+(
const size_type n)
const
2436 for (size_type i = 0; i < n; ++i)
2446template <
typename number>
2450 return const_iterator(
this, 0);
2454template <
typename number>
2458 return const_iterator(
this);
2462template <
typename number>
2466 return iterator(
this, 0);
2470template <
typename number>
2474 return iterator(
this, cols->rowstart[cols->rows]);
2478template <
typename number>
2484 return const_iterator(
this, cols->rowstart[r]);
2489template <
typename number>
2495 return const_iterator(
this, cols->rowstart[r + 1]);
2500template <
typename number>
2506 return iterator(
this, cols->rowstart[r]);
2511template <
typename number>
2517 return iterator(
this, cols->rowstart[r + 1]);
2522template <
typename number>
2523template <
typename StreamType>
2527 const bool diagonal_first)
const
2532 bool hanging_diagonal =
false;
2535 for (size_type i = 0; i < cols->rows; ++i)
2537 for (size_type j = cols->rowstart[i]; j < cols->rowstart[i + 1]; ++j)
2539 if (!diagonal_first && i == cols->colnums[j])
2542 hanging_diagonal =
true;
2546 if (hanging_diagonal && cols->colnums[j] > i)
2549 out <<
' ' << i <<
',' << i <<
':' <<
diagonal;
2551 out <<
'(' << i <<
',' << i <<
") " <<
diagonal
2553 hanging_diagonal =
false;
2556 out <<
' ' << i <<
',' << cols->colnums[j] <<
':' << val[j];
2558 out <<
'(' << i <<
',' << cols->colnums[j] <<
") " << val[j]
2562 if (hanging_diagonal)
2565 out <<
' ' << i <<
',' << i <<
':' <<
diagonal;
2567 out <<
'(' << i <<
',' << i <<
") " <<
diagonal << std::endl;
2568 hanging_diagonal =
false;
2576template <
typename number>
2585template <
typename number>
const Reference & operator-=(const number n) const
Reference(const Accessor *accessor, const bool dummy)
const Reference & operator/=(const number n) const
const Accessor * accessor
const Reference & operator+=(const number n) const
const Reference & operator*=(const number n) const
const Reference & operator=(const number n) const
SparseMatrix< number > MatrixType
MatrixType & get_matrix() const
Accessor(MatrixType *matrix, const std::size_t index)
Accessor(MatrixType *matrix)
Accessor(MatrixType *matrix)
Accessor(const SparseMatrixIterators::Accessor< number, false > &a)
Accessor(MatrixType *matrix, const std::size_t index_within_matrix)
const MatrixType & get_matrix() const
const SparseMatrix< number > MatrixType
const SparseMatrix< number > & get_matrix() const
const Iterator< number, Constness > & operator=(const SparseMatrixIterators::Iterator< number, false > &i)
bool operator>(const Iterator &) const
bool operator==(const Iterator &) const
int operator-(const Iterator &p) const
bool operator<(const Iterator &) const
Iterator(MatrixType *matrix)
const Accessor< number, Constness > & value_type
Iterator operator+(const size_type n) const
const Accessor< number, Constness > & operator*() const
Accessor< number, Constness > accessor
Iterator(const SparseMatrixIterators::Iterator< number, false > &i)
Iterator(MatrixType *matrix, const std::size_t index_within_matrix)
const Accessor< number, Constness > * operator->() const
typename Accessor< number, Constness >::MatrixType MatrixType
bool operator!=(const Iterator &) const
void set(const size_type row, const std::vector< size_type > &col_indices, const std::vector< number2 > &values, const bool elide_zero_values=false)
somenumber matrix_scalar_product(const Vector< somenumber > &u, const Vector< somenumber > &v) const
void TSOR_step(Vector< somenumber > &v, const Vector< somenumber > &b, const number omega=1.) const
void precondition_Jacobi(Vector< somenumber > &dst, const Vector< somenumber > &src, const number omega=1.) const
SparseMatrix(const SparsityPattern &sparsity)
void add(const number factor, const SparseMatrix< somenumber > &matrix)
std::size_t n_nonzero_elements() const
size_type get_row_length(const size_type row) const
void SOR_step(Vector< somenumber > &v, const Vector< somenumber > &b, const number omega=1.) const
somenumber residual(Vector< somenumber > &dst, const Vector< somenumber > &x, const Vector< somenumber > &b) const
number & diag_element(const size_type i)
void Tmmult(SparseMatrix< numberC > &C, const SparseMatrix< numberB > &B, const Vector< number > &V=Vector< number >(), const bool rebuild_sparsity_pattern=true) const
void Tvmult(OutVector &dst, const InVector &src) const
const_iterator begin(const size_type r) const
SparseMatrix< number > & operator=(const SparseMatrix< number > &)
const_iterator end() const
void SSOR_step(Vector< somenumber > &v, const Vector< somenumber > &b, const number omega=1.) const
void set(const std::vector< size_type > &indices, const FullMatrix< number2 > &full_matrix, const bool elide_zero_values=false)
void print_as_numpy_arrays(std::ostream &out, const unsigned int precision=9) const
void mmult(SparseMatrix< numberC > &C, const SparseMatrix< numberB > &B, const Vector< number > &V=Vector< number >(), const bool rebuild_sparsity_pattern=true) const
ObserverPointer< const SparsityPattern, SparseMatrix< number > > cols
const SparsityPattern & get_sparsity_pattern() const
void set(const size_type i, const size_type j, const number value)
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
const_iterator begin() const
SparseMatrix< number > & copy_from(const SparseMatrix< somenumber > &source)
void precondition_TSOR(Vector< somenumber > &dst, const Vector< somenumber > &src, const number omega=1.) const
virtual ~SparseMatrix() override
void vmult_add(OutVector &dst, const InVector &src) const
number diag_element(const size_type i) const
void add(const std::vector< size_type > &row_indices, const std::vector< size_type > &col_indices, const FullMatrix< number2 > &full_matrix, const bool elide_zero_values=true)
void SSOR(Vector< somenumber > &v, const number omega=1.) const
somenumber matrix_norm_square(const Vector< somenumber > &v) const
SparseMatrix(SparseMatrix< number > &&m) noexcept
void TSOR(Vector< somenumber > &v, const number omega=1.) const
SparseMatrix< number > & copy_from(const TrilinosWrappers::SparseMatrix &matrix)
SparseMatrix(const SparsityPattern &sparsity, const IdentityMatrix &id)
void print_pattern(std::ostream &out, const double threshold=0.) const
void vmult(OutVector &dst, const InVector &src) const
std::size_t memory_consumption() const
void PSOR(Vector< somenumber > &v, const std::vector< size_type > &permutation, const std::vector< size_type > &inverse_permutation, const number omega=1.) const
void block_write(std::ostream &out) const
void copy_from(const ForwardIterator begin, const ForwardIterator end)
const number & operator()(const size_type i, const size_type j) const
number & operator()(const size_type i, const size_type j)
void block_read(std::istream &in)
number el(const size_type i, const size_type j) const
SparseMatrix< number > & operator=(SparseMatrix< number > &&m) noexcept
void precondition_SSOR(Vector< somenumber > &dst, const Vector< somenumber > &src, const number omega=1., const std::vector< std::size_t > &pos_right_of_diagonal=std::vector< std::size_t >()) const
void print(StreamType &out, const bool across=false, const bool diagonal_first=true) const
const_iterator end(const size_type r) const
void SOR(Vector< somenumber > &v, const number omega=1.) const
SparseMatrix & operator/=(const number factor)
std::unique_ptr< number[]> val
iterator begin(const size_type r)
void precondition_SOR(Vector< somenumber > &dst, const Vector< somenumber > &src, const number omega=1.) const
typename numbers::NumberTraits< number >::real_type real_type
void reinit(const SparseMatrix< number2 > &sparse_matrix)
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)
void add(const size_type i, const size_type j, const number value)
void copy_from(const FullMatrix< somenumber > &matrix)
iterator end(const size_type r)
SparseMatrix & operator*=(const number factor)
void compress(VectorOperation::values)
SparseMatrix(const SparseMatrix &)
SparseMatrix< number > & operator=(const IdentityMatrix &id)
SparseMatrix & operator=(const double d)
real_type frobenius_norm() const
void Tvmult_add(OutVector &dst, const InVector &src) const
real_type linfty_norm() const
void add(const size_type row, const std::vector< size_type > &col_indices, const std::vector< number2 > &values, const bool elide_zero_values=true)
void set(const std::vector< size_type > &row_indices, const std::vector< size_type > &col_indices, const FullMatrix< number2 > &full_matrix, const bool elide_zero_values=false)
real_type l1_norm() const
std::size_t n_actually_nonzero_elements(const double threshold=0.) const
void add(const std::vector< size_type > &indices, const FullMatrix< number2 > &full_matrix, const bool elide_zero_values=true)
void TPSOR(Vector< somenumber > &v, const std::vector< size_type > &permutation, const std::vector< size_type > &inverse_permutation, const number omega=1.) const
virtual void reinit(const SparsityPattern &sparsity)
void set(const size_type row, const size_type n_cols, const size_type *col_indices, const number2 *values, const bool elide_zero_values=false)
void Jacobi_step(Vector< somenumber > &v, const Vector< somenumber > &b, const number omega=1.) const
friend class ChunkSparsityPatternIterators::Accessor
SparsityPatternIterators::size_type size_type
static constexpr size_type invalid_entry
#define DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_CLOSE
static ::ExceptionBase & ExcDifferentSparsityPatterns()
static ::ExceptionBase & ExcInvalidIndex(int arg1, int arg2)
#define Assert(cond, exc)
static ::ExceptionBase & ExcSourceEqualsDestination()
#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 & ExcIteratorRange(int arg1, int arg2)
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
static ::ExceptionBase & ExcNotInitialized()
static ::ExceptionBase & ExcNotQuadratic()
@ matrix
Contents is actually a matrix.
@ diagonal
Matrix is diagonal.
types::global_dof_index size_type
VectorType::value_type * end(VectorType &V)
T sum(const T &t, const MPI_Comm mpi_communicator)
unsigned int global_dof_index
static const bool zero_addition_can_be_elided
typename ::SparseMatrixIterators::Iterator< number, Constness >::difference_type difference_type
forward_iterator_tag iterator_category
typename ::SparseMatrixIterators::Iterator< number, Constness >::value_type value_type