16#ifndef dealii_sparse_matrix_h
17# define dealii_sparse_matrix_h
29# ifdef DEAL_II_WITH_MPI
40template <
typename number>
42template <
typename number>
44template <
typename Matrix>
46template <
typename number>
48# ifdef DEAL_II_WITH_MPI
53 template <
typename Number>
60# ifdef DEAL_II_WITH_TRILINOS
85 template <
typename number,
bool Constness>
98 template <
typename number,
bool Constness>
130 template <
typename number>
180 template <
typename,
bool>
191 template <
typename number>
231 operator number()
const;
313 template <
typename,
bool>
348 template <
typename number,
bool Constness>
494template <
typename number>
781 template <
typename number2>
783 set(
const std::vector<size_type> &indices,
785 const bool elide_zero_values =
false);
792 template <
typename number2>
794 set(
const std::vector<size_type> &row_indices,
795 const std::vector<size_type> &col_indices,
797 const bool elide_zero_values =
false);
809 template <
typename number2>
812 const std::vector<size_type> &col_indices,
813 const std::vector<number2> &
values,
814 const bool elide_zero_values =
false);
825 template <
typename number2>
831 const bool elide_zero_values =
false);
855 template <
typename number2>
857 add(
const std::vector<size_type> &indices,
859 const bool elide_zero_values =
true);
866 template <
typename number2>
868 add(
const std::vector<size_type> &row_indices,
869 const std::vector<size_type> &col_indices,
871 const bool elide_zero_values =
true);
882 template <
typename number2>
885 const std::vector<size_type> &col_indices,
886 const std::vector<number2> &
values,
887 const bool elide_zero_values =
true);
898 template <
typename number2>
904 const bool elide_zero_values =
true,
905 const bool col_indices_are_sorted =
false);
950 template <
typename somenumber>
970 template <
typename ForwardIterator>
983 template <
typename somenumber>
987# ifdef DEAL_II_WITH_TRILINOS
1012 template <
typename somenumber>
1099 template <
class OutVector,
class InVector>
1101 vmult(OutVector &dst,
const InVector &src)
const;
1118 template <
class OutVector,
class InVector>
1120 Tvmult(OutVector &dst,
const InVector &src)
const;
1138 template <
class OutVector,
class InVector>
1157 template <
class OutVector,
class InVector>
1178 template <
typename somenumber>
1187 template <
typename somenumber>
1201 template <
typename somenumber>
1242 template <
typename numberB,
typename numberC>
1247 const bool rebuild_sparsity_pattern =
true)
const;
1273 template <
typename numberB,
typename numberC>
1278 const bool rebuild_sparsity_pattern =
true)
const;
1323 template <
typename somenumber>
1327 const number omega = 1.)
const;
1335 template <
typename somenumber>
1339 const number omega = 1.,
1340 const std::vector<std::size_t> &pos_right_of_diagonal =
1341 std::vector<std::size_t>())
const;
1346 template <
typename somenumber>
1350 const number om = 1.)
const;
1355 template <
typename somenumber>
1359 const number om = 1.)
const;
1366 template <
typename somenumber>
1374 template <
typename somenumber>
1382 template <
typename somenumber>
1396 template <
typename somenumber>
1399 const std::vector<size_type> &permutation,
1400 const std::vector<size_type> &inverse_permutation,
1401 const number om = 1.)
const;
1413 template <
typename somenumber>
1416 const std::vector<size_type> &permutation,
1417 const std::vector<size_type> &inverse_permutation,
1418 const number om = 1.)
const;
1425 template <
typename somenumber>
1429 const number om = 1.)
const;
1435 template <
typename somenumber>
1439 const number om = 1.)
const;
1445 template <
typename somenumber>
1449 const number om = 1.)
const;
1455 template <
typename somenumber>
1459 const number om = 1.)
const;
1545 template <
class StreamType>
1548 const bool across =
false,
1549 const bool diagonal_first =
true)
const;
1573 const unsigned int precision = 3,
1574 const bool scientific =
true,
1575 const unsigned int width = 0,
1576 const char * zero_string =
" ",
1577 const double denominator = 1.)
const;
1597 const unsigned int precision = 9)
const;
1642 <<
"You are trying to access the matrix entry with index <"
1643 << arg1 <<
',' << arg2
1644 <<
">, but this entry does not exist in the sparsity pattern "
1647 "The most common cause for this problem is that you used "
1648 "a method to build the sparsity pattern that did not "
1649 "(completely) take into account all of the entries you "
1650 "will later try to write into. An example would be "
1651 "building a sparsity pattern that does not include "
1652 "the entries you will write into due to constraints "
1653 "on degrees of freedom such as hanging nodes or periodic "
1654 "boundary conditions. In such cases, building the "
1655 "sparsity pattern will succeed, but you will get errors "
1656 "such as the current one at one point or other when "
1657 "trying to write into the entries of the matrix.");
1662 "When copying one sparse matrix into another, "
1663 "or when adding one sparse matrix to another, "
1664 "both matrices need to refer to the same "
1665 "sparsity pattern.");
1672 <<
"The iterators denote a range of " << arg1
1673 <<
" elements, but the given number of rows was " << arg2);
1678 "You are attempting an operation on two matrices that "
1679 "are the same object, but the operation requires that the "
1680 "two objects are in fact different.");
1719 std::unique_ptr<number[]>
val;
1730 template <
typename somenumber>
1732 template <
typename somenumber>
1742 template <
typename,
bool>
1744 template <
typename,
bool>
1747# ifdef DEAL_II_WITH_MPI
1749 template <
typename Number>
1762template <
typename number>
1771template <
typename number>
1781template <
typename number>
1796 ExcInvalidIndex(i, j));
1805template <
typename number>
1806template <
typename number2>
1810 const bool elide_zero_values)
1816 for (
size_type i = 0; i < indices.size(); ++i)
1826template <
typename number>
1827template <
typename number2>
1830 const std::vector<size_type> &col_indices,
1832 const bool elide_zero_values)
1839 for (
size_type i = 0; i < row_indices.size(); ++i)
1849template <
typename number>
1850template <
typename number2>
1853 const std::vector<size_type> &col_indices,
1854 const std::vector<number2> &
values,
1855 const bool elide_zero_values)
1869template <
typename number>
1877 if (
value == number())
1887 ExcInvalidIndex(i, j));
1896template <
typename number>
1897template <
typename number2>
1901 const bool elide_zero_values)
1907 for (
size_type i = 0; i < indices.size(); ++i)
1917template <
typename number>
1918template <
typename number2>
1921 const std::vector<size_type> &col_indices,
1923 const bool elide_zero_values)
1930 for (
size_type i = 0; i < row_indices.size(); ++i)
1940template <
typename number>
1941template <
typename number2>
1944 const std::vector<size_type> &col_indices,
1945 const std::vector<number2> &
values,
1946 const bool elide_zero_values)
1960template <
typename number>
1967 number * val_ptr = val.get();
1968 const number *
const end_ptr = val.get() + cols->n_nonzero_elements();
1970 while (val_ptr != end_ptr)
1971 *val_ptr++ *= factor;
1978template <
typename number>
1986 const number factor_inv = number(1.) / factor;
1988 number * val_ptr = val.get();
1989 const number *
const end_ptr = val.get() + cols->n_nonzero_elements();
1991 while (val_ptr != end_ptr)
1992 *val_ptr++ *= factor_inv;
1999template <
typename number>
2000inline const number &
2005 ExcInvalidIndex(i, j));
2006 return val[cols->operator()(i, j)];
2011template <
typename number>
2017 ExcInvalidIndex(i, j));
2018 return val[cols->operator()(i, j)];
2023template <
typename number>
2038template <
typename number>
2048 return val[cols->rowstart[i]];
2053template <
typename number>
2063 return val[cols->rowstart[i]];
2068template <
typename number>
2069template <
typename ForwardIterator>
2072 const ForwardIterator
end)
2075 ExcIteratorRange(std::distance(
begin,
end), m()));
2079 using inner_iterator =
2080 typename std::iterator_traits<ForwardIterator>::value_type::const_iterator;
2082 for (ForwardIterator i =
begin; i !=
end; ++i, ++
row)
2084 const inner_iterator end_of_row = i->end();
2085 for (inner_iterator j = i->begin(); j != end_of_row; ++j)
2087 set(
row, j->first, j->second);
2097 template <
typename number>
2099 const std::size_t index_within_matrix)
2101 index_within_matrix)
2107 template <
typename number>
2108 inline Accessor<number, true>::Accessor(
const MatrixType *
matrix)
2115 template <
typename number>
2116 inline Accessor<number, true>::Accessor(
2119 ,
matrix(&a.get_matrix())
2124 template <
typename number>
2126 Accessor<number, true>::value()
const
2129 return matrix->val[linear_index];
2134 template <
typename number>
2135 inline const typename Accessor<number, true>::MatrixType &
2136 Accessor<number, true>::get_matrix()
const
2143 template <
typename number>
2144 inline Accessor<number, false>::Reference::Reference(
const Accessor *accessor,
2146 : accessor(accessor)
2150 template <
typename number>
2151 inline Accessor<number, false>::Reference::operator number()
const
2154 accessor->matrix->n_nonzero_elements());
2155 return accessor->matrix->val[accessor->linear_index];
2160 template <
typename number>
2161 inline const typename Accessor<number, false>::Reference &
2162 Accessor<number, false>::Reference::operator=(
const number n)
const
2165 accessor->matrix->n_nonzero_elements());
2166 accessor->matrix->val[accessor->linear_index] = n;
2172 template <
typename number>
2173 inline const typename Accessor<number, false>::Reference &
2174 Accessor<number, false>::Reference::operator+=(
const number n)
const
2177 accessor->matrix->n_nonzero_elements());
2178 accessor->matrix->val[accessor->linear_index] += n;
2184 template <
typename number>
2185 inline const typename Accessor<number, false>::Reference &
2186 Accessor<number, false>::Reference::operator-=(
const number n)
const
2189 accessor->matrix->n_nonzero_elements());
2190 accessor->matrix->val[accessor->linear_index] -= n;
2196 template <
typename number>
2197 inline const typename Accessor<number, false>::Reference &
2198 Accessor<number, false>::Reference::operator*=(
const number n)
const
2201 accessor->matrix->n_nonzero_elements());
2202 accessor->matrix->val[accessor->linear_index] *= n;
2208 template <
typename number>
2209 inline const typename Accessor<number, false>::Reference &
2210 Accessor<number, false>::Reference::operator/=(
const number n)
const
2213 accessor->matrix->n_nonzero_elements());
2214 accessor->matrix->val[accessor->linear_index] /= n;
2220 template <
typename number>
2221 inline Accessor<number, false>::Accessor(MatrixType *
matrix,
2222 const std::size_t index)
2229 template <
typename number>
2230 inline Accessor<number, false>::Accessor(MatrixType *
matrix)
2237 template <
typename number>
2238 inline typename Accessor<number, false>::Reference
2239 Accessor<number, false>::value()
const
2241 return Reference(
this,
true);
2246 template <
typename number>
2247 inline typename Accessor<number, false>::MatrixType &
2248 Accessor<number, false>::get_matrix()
const
2255 template <
typename number,
bool Constness>
2256 inline Iterator<number, Constness>::Iterator(MatrixType *
matrix,
2257 const std::size_t index)
2258 : accessor(
matrix, index)
2263 template <
typename number,
bool Constness>
2264 inline Iterator<number, Constness>::Iterator(MatrixType *
matrix)
2270 template <
typename number,
bool Constness>
2271 inline Iterator<number, Constness>::Iterator(
2278 template <
typename number,
bool Constness>
2279 inline const Iterator<number, Constness> &
2280 Iterator<number, Constness>::
2289 template <
typename number,
bool Constness>
2290 inline Iterator<number, Constness> &
2298 template <
typename number,
bool Constness>
2299 inline Iterator<number, Constness>
2302 const Iterator iter = *
this;
2308 template <
typename number,
bool Constness>
2316 template <
typename number,
bool Constness>
2317 inline const Accessor<number, Constness> *Iterator<number, Constness>::
2324 template <
typename number,
bool Constness>
2328 return (accessor == other.accessor);
2332 template <
typename number,
bool Constness>
2336 return !(*
this == other);
2340 template <
typename number,
bool Constness>
2344 Assert(&accessor.get_matrix() == &other.accessor.get_matrix(),
2347 return (accessor < other.accessor);
2351 template <
typename number,
bool Constness>
2355 return (other < *
this);
2359 template <
typename number,
bool Constness>
2363 Assert(&accessor.get_matrix() == &other.accessor.get_matrix(),
2366 return (*this)->linear_index - other->linear_index;
2371 template <
typename number,
bool Constness>
2372 inline Iterator<number, Constness>
2386template <
typename number>
2390 return const_iterator(
this, 0);
2394template <
typename number>
2398 return const_iterator(
this);
2402template <
typename number>
2406 return iterator(
this, 0);
2410template <
typename number>
2414 return iterator(
this, cols->rowstart[cols->rows]);
2418template <
typename number>
2424 return const_iterator(
this, cols->rowstart[r]);
2429template <
typename number>
2435 return const_iterator(
this, cols->rowstart[r + 1]);
2440template <
typename number>
2446 return iterator(
this, cols->rowstart[r]);
2451template <
typename number>
2457 return iterator(
this, cols->rowstart[r + 1]);
2462template <
typename number>
2463template <
class StreamType>
2467 const bool diagonal_first)
const
2472 bool hanging_diagonal =
false;
2475 for (
size_type i = 0; i < cols->rows; ++i)
2477 for (
size_type j = cols->rowstart[i]; j < cols->rowstart[i + 1]; ++j)
2479 if (!diagonal_first && i == cols->colnums[j])
2482 hanging_diagonal =
true;
2486 if (hanging_diagonal && cols->colnums[j] > i)
2489 out <<
' ' << i <<
',' << i <<
':' <<
diagonal;
2491 out <<
'(' << i <<
',' << i <<
") " <<
diagonal
2493 hanging_diagonal =
false;
2496 out <<
' ' << i <<
',' << cols->colnums[j] <<
':' << val[j];
2498 out <<
"(" << i <<
"," << cols->colnums[j] <<
") " << val[j]
2502 if (hanging_diagonal)
2505 out <<
' ' << i <<
',' << i <<
':' <<
diagonal;
2507 out <<
'(' << i <<
',' << i <<
") " <<
diagonal << std::endl;
2508 hanging_diagonal =
false;
2516template <
typename number>
2525template <
typename number>
bool operator!=(const AlignedVector< T > &lhs, const AlignedVector< T > &rhs)
bool operator==(const AlignedVector< T > &lhs, const AlignedVector< T > &rhs)
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 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
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
void compress(::VectorOperation::values)
SparseMatrix< number > & operator=(const SparseMatrix< number > &)
const_iterator end() const
void SSOR_step(Vector< somenumber > &v, const Vector< somenumber > &b, const number om=1.) const
void set(const std::vector< size_type > &indices, const FullMatrix< number2 > &full_matrix, const bool elide_zero_values=false)
void TSOR_step(Vector< somenumber > &v, const Vector< somenumber > &b, const number om=1.) const
void mmult(SparseMatrix< numberC > &C, const SparseMatrix< numberB > &B, const Vector< number > &V=Vector< number >(), const bool rebuild_sparsity_pattern=true) const
const SparsityPattern & get_sparsity_pattern() const
void set(const size_type i, const size_type j, const number value)
const_iterator begin() const
SparseMatrix< number > & copy_from(const SparseMatrix< somenumber > &source)
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
SparseMatrix< number > & copy_from(const TrilinosWrappers::SparseMatrix &matrix)
SparseMatrix(const SparsityPattern &sparsity, const IdentityMatrix &id)
void SOR_step(Vector< somenumber > &v, const Vector< somenumber > &b, const number om=1.) const
void SOR(Vector< somenumber > &v, const number om=1.) const
void vmult(OutVector &dst, const InVector &src) const
std::size_t memory_consumption() const
void copy_from(const ForwardIterator begin, const ForwardIterator end)
void PSOR(Vector< somenumber > &v, const std::vector< size_type > &permutation, const std::vector< size_type > &inverse_permutation, const number om=1.) const
void precondition_TSOR(Vector< somenumber > &dst, const Vector< somenumber > &src, const number om=1.) const
const number & operator()(const size_type i, const size_type j) const
number & operator()(const size_type i, const size_type j)
void Jacobi_step(Vector< somenumber > &v, const Vector< somenumber > &b, const number om=1.) const
number el(const size_type i, const size_type j) const
void TPSOR(Vector< somenumber > &v, const std::vector< size_type > &permutation, const std::vector< size_type > &inverse_permutation, const number om=1.) 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
const_iterator end(const size_type r) const
SparseMatrix & operator/=(const number factor)
void TSOR(Vector< somenumber > &v, const number om=1.) const
iterator begin(const size_type r)
typename numbers::NumberTraits< number >::real_type real_type
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)
void precondition_SOR(Vector< somenumber > &dst, const Vector< somenumber > &src, const number om=1.) const
SparseMatrix & operator*=(const number factor)
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)
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)
std::enable_if< std::is_floating_point< T >::value &&std::is_floating_point< U >::value, typenameProductType< std::complex< T >, std::complex< U > >::type >::type operator*(const std::complex< T > &left, const std::complex< U > &right)
#define DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_CLOSE
__global__ void set(Number *val, const Number s, const size_type N)
static ::ExceptionBase & ExcDifferentSparsityPatterns()
void print_as_numpy_arrays(std::ostream &out, const unsigned int precision=9) const
static ::ExceptionBase & ExcInvalidIndex(int arg1, int arg2)
SmartPointer< const SparsityPattern, SparseMatrix< number > > cols
static const bool zero_addition_can_be_elided
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
#define Assert(cond, exc)
static ::ExceptionBase & ExcSourceEqualsDestination()
void print_pattern(std::ostream &out, const double threshold=0.) const
#define AssertIsFinite(number)
void block_write(std::ostream &out) const
#define DeclException2(Exception2, type1, type2, outsequence)
void block_read(std::istream &in)
void print(StreamType &out, const bool across=false, const bool diagonal_first=true) const
static ::ExceptionBase & ExcNeedsSparsityPattern()
#define AssertIndexRange(index, range)
#define DeclExceptionMsg(Exception, defaulttext)
Expression operator>(const Expression &lhs, const Expression &rhs)
static ::ExceptionBase & ExcInternalError()
std::unique_ptr< number[]> val
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()
static const size_type invalid_entry
@ matrix
Contents is actually a matrix.
@ diagonal
Matrix is diagonal.
SymmetricTensor< 2, dim, Number > C(const Tensor< 2, dim, Number > &F)
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
types::global_dof_index size_type
VectorType::value_type * end(VectorType &V)
VectorType::value_type * begin(VectorType &V)
T sum(const T &t, const MPI_Comm &mpi_communicator)
unsigned int global_dof_index
BarycentricPolynomial< dim, Number1 > operator-(const Number2 &a, const BarycentricPolynomial< dim, Number1 > &bp)
BarycentricPolynomial< dim, Number1 > operator+(const Number2 &a, const BarycentricPolynomial< dim, Number1 > &bp)
SynchronousIterators< Iterators > operator++(SynchronousIterators< Iterators > &a)
bool operator<(const SynchronousIterators< Iterators > &a, const SynchronousIterators< Iterators > &b)