|
Reference documentation for deal.II version 9.2.0
|
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\)
\(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\)
\(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\)
\(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
Go to the documentation of this file.
16 #ifndef dealii_sparse_matrix_h
17 # define dealii_sparse_matrix_h
29 # ifdef DEAL_II_WITH_MPI
40 template <
typename number>
42 template <
typename number>
44 template <
typename Matrix>
46 template <
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;
237 operator=(
const number n)
const;
243 operator+=(
const number n)
const;
249 operator-=(
const number n)
const;
255 operator*=(
const number n)
const;
261 operator/=(
const number n)
const;
313 template <
typename,
bool>
348 template <
typename number,
bool Constness>
496 template <
typename number>
783 template <
typename number2>
785 set(
const std::vector<size_type> &indices,
787 const bool elide_zero_values =
false);
794 template <
typename number2>
796 set(
const std::vector<size_type> &row_indices,
797 const std::vector<size_type> &col_indices,
799 const bool elide_zero_values =
false);
811 template <
typename number2>
814 const std::vector<size_type> &col_indices,
815 const std::vector<number2> & values,
816 const bool elide_zero_values =
false);
827 template <
typename number2>
832 const number2 * values,
833 const bool elide_zero_values =
false);
857 template <
typename number2>
859 add(
const std::vector<size_type> &indices,
861 const bool elide_zero_values =
true);
868 template <
typename number2>
870 add(
const std::vector<size_type> &row_indices,
871 const std::vector<size_type> &col_indices,
873 const bool elide_zero_values =
true);
884 template <
typename number2>
887 const std::vector<size_type> &col_indices,
888 const std::vector<number2> & values,
889 const bool elide_zero_values =
true);
900 template <
typename number2>
905 const number2 * values,
906 const bool elide_zero_values =
true,
907 const bool col_indices_are_sorted =
false);
952 template <
typename somenumber>
972 template <
typename ForwardIterator>
981 template <
typename somenumber>
985 # ifdef DEAL_II_WITH_TRILINOS
1010 template <
typename somenumber>
1097 template <
class OutVector,
class InVector>
1099 vmult(OutVector &dst,
const InVector &src)
const;
1116 template <
class OutVector,
class InVector>
1118 Tvmult(OutVector &dst,
const InVector &src)
const;
1136 template <
class OutVector,
class InVector>
1138 vmult_add(OutVector &dst,
const InVector &src)
const;
1155 template <
class OutVector,
class InVector>
1157 Tvmult_add(OutVector &dst,
const InVector &src)
const;
1176 template <
typename somenumber>
1185 template <
typename somenumber>
1188 const Vector<somenumber> &v)
const;
1199 template <
typename somenumber>
1202 const Vector<somenumber> &x,
1203 const Vector<somenumber> &
b)
const;
1240 template <
typename numberB,
typename numberC>
1245 const bool rebuild_sparsity_pattern =
true)
const;
1271 template <
typename numberB,
typename numberC>
1276 const bool rebuild_sparsity_pattern =
true)
const;
1321 template <
typename somenumber>
1324 const Vector<somenumber> &src,
1325 const number omega = 1.)
const;
1333 template <
typename somenumber>
1336 const Vector<somenumber> & src,
1337 const number omega = 1.,
1338 const std::vector<std::size_t> &pos_right_of_diagonal =
1339 std::vector<std::size_t>())
const;
1344 template <
typename somenumber>
1347 const Vector<somenumber> &src,
1348 const number om = 1.)
const;
1353 template <
typename somenumber>
1356 const Vector<somenumber> &src,
1357 const number om = 1.)
const;
1364 template <
typename somenumber>
1366 SSOR(Vector<somenumber> &v,
const number omega = 1.)
const;
1372 template <
typename somenumber>
1374 SOR(Vector<somenumber> &v,
const number om = 1.)
const;
1380 template <
typename somenumber>
1382 TSOR(Vector<somenumber> &v,
const number om = 1.)
const;
1394 template <
typename somenumber>
1396 PSOR(Vector<somenumber> & v,
1397 const std::vector<size_type> &permutation,
1398 const std::vector<size_type> &inverse_permutation,
1399 const number om = 1.)
const;
1411 template <
typename somenumber>
1413 TPSOR(Vector<somenumber> & v,
1414 const std::vector<size_type> &permutation,
1415 const std::vector<size_type> &inverse_permutation,
1416 const number om = 1.)
const;
1423 template <
typename somenumber>
1426 const Vector<somenumber> &
b,
1427 const number om = 1.)
const;
1433 template <
typename somenumber>
1436 const Vector<somenumber> &
b,
1437 const number om = 1.)
const;
1443 template <
typename somenumber>
1446 const Vector<somenumber> &
b,
1447 const number om = 1.)
const;
1453 template <
typename somenumber>
1456 const Vector<somenumber> &
b,
1457 const number om = 1.)
const;
1543 template <
class StreamType>
1545 print(StreamType &out,
1546 const bool across =
false,
1547 const bool diagonal_first =
true)
const;
1571 const unsigned int precision = 3,
1572 const bool scientific =
true,
1573 const unsigned int width = 0,
1574 const char * zero_string =
" ",
1575 const double denominator = 1.)
const;
1583 print_pattern(std::ostream &out,
const double threshold = 0.)
const;
1595 const unsigned int precision = 9)
const;
1640 <<
"You are trying to access the matrix entry with index <"
1641 << arg1 <<
',' << arg2
1642 <<
">, but this entry does not exist in the sparsity pattern "
1645 "The most common cause for this problem is that you used "
1646 "a method to build the sparsity pattern that did not "
1647 "(completely) take into account all of the entries you "
1648 "will later try to write into. An example would be "
1649 "building a sparsity pattern that does not include "
1650 "the entries you will write into due to constraints "
1651 "on degrees of freedom such as hanging nodes or periodic "
1652 "boundary conditions. In such cases, building the "
1653 "sparsity pattern will succeed, but you will get errors "
1654 "such as the current one at one point or other when "
1655 "trying to write into the entries of the matrix.");
1660 "When copying one sparse matrix into another, "
1661 "or when adding one sparse matrix to another, "
1662 "both matrices need to refer to the same "
1663 "sparsity pattern.");
1670 <<
"The iterators denote a range of " << arg1
1671 <<
" elements, but the given number of rows was " << arg2);
1676 "You are attempting an operation on two matrices that "
1677 "are the same object, but the operation requires that the "
1678 "two objects are in fact different.");
1717 std::unique_ptr<number[]>
val;
1728 template <
typename somenumber>
1730 template <
typename somenumber>
1740 template <
typename,
bool>
1742 template <
typename,
bool>
1745 # ifdef DEAL_II_WITH_MPI
1747 template <
typename Number>
1760 template <
typename number>
1769 template <
typename number>
1779 template <
typename number>
1794 ExcInvalidIndex(i, j));
1803 template <
typename number>
1804 template <
typename number2>
1808 const bool elide_zero_values)
1810 Assert(indices.size() == values.
m(),
1814 for (
size_type i = 0; i < indices.size(); ++i)
1824 template <
typename number>
1825 template <
typename number2>
1828 const std::vector<size_type> &col_indices,
1830 const bool elide_zero_values)
1832 Assert(row_indices.size() == values.
m(),
1834 Assert(col_indices.size() == values.
n(),
1837 for (
size_type i = 0; i < row_indices.size(); ++i)
1847 template <
typename number>
1848 template <
typename number2>
1851 const std::vector<size_type> &col_indices,
1852 const std::vector<number2> & values,
1853 const bool elide_zero_values)
1855 Assert(col_indices.size() == values.size(),
1867 template <
typename number>
1875 if (
value == number())
1885 ExcInvalidIndex(i, j));
1894 template <
typename number>
1895 template <
typename number2>
1899 const bool elide_zero_values)
1901 Assert(indices.size() == values.
m(),
1905 for (
size_type i = 0; i < indices.size(); ++i)
1915 template <
typename number>
1916 template <
typename number2>
1919 const std::vector<size_type> &col_indices,
1921 const bool elide_zero_values)
1923 Assert(row_indices.size() == values.
m(),
1925 Assert(col_indices.size() == values.
n(),
1928 for (
size_type i = 0; i < row_indices.size(); ++i)
1938 template <
typename number>
1939 template <
typename number2>
1942 const std::vector<size_type> &col_indices,
1943 const std::vector<number2> & values,
1944 const bool elide_zero_values)
1946 Assert(col_indices.size() == values.size(),
1958 template <
typename number>
1965 number * val_ptr = val.get();
1966 const number *
const end_ptr = val.get() + cols->n_nonzero_elements();
1968 while (val_ptr != end_ptr)
1969 *val_ptr++ *= factor;
1976 template <
typename number>
1984 const number factor_inv = number(1.) / factor;
1986 number * val_ptr = val.get();
1987 const number *
const end_ptr = val.get() + cols->n_nonzero_elements();
1989 while (val_ptr != end_ptr)
1990 *val_ptr++ *= factor_inv;
1997 template <
typename number>
1998 inline const number &
2003 ExcInvalidIndex(i, j));
2004 return val[cols->operator()(i, j)];
2009 template <
typename number>
2015 ExcInvalidIndex(i, j));
2016 return val[cols->operator()(i, j)];
2021 template <
typename number>
2036 template <
typename number>
2046 return val[cols->rowstart[i]];
2051 template <
typename number>
2061 return val[cols->rowstart[i]];
2066 template <
typename number>
2067 template <
typename ForwardIterator>
2070 const ForwardIterator
end)
2073 ExcIteratorRange(std::distance(
begin,
end), m()));
2077 using inner_iterator =
2078 typename std::iterator_traits<ForwardIterator>::value_type::const_iterator;
2080 for (ForwardIterator i =
begin; i !=
end; ++i, ++
row)
2082 const inner_iterator end_of_row = i->end();
2083 for (inner_iterator j = i->begin(); j != end_of_row; ++j)
2085 set(
row, j->first, j->second);
2095 template <
typename number>
2097 const std::size_t index_within_matrix)
2099 index_within_matrix)
2105 template <
typename number>
2113 template <
typename number>
2117 ,
matrix(&a.get_matrix())
2122 template <
typename number>
2127 return matrix->val[linear_index];
2132 template <
typename number>
2133 inline const typename Accessor<number, true>::MatrixType &
2134 Accessor<number, true>::get_matrix()
const
2141 template <
typename number>
2142 inline Accessor<number, false>::Reference::Reference(
const Accessor *accessor,
2144 : accessor(accessor)
2148 template <
typename number>
2149 inline Accessor<number, false>::Reference::operator number()
const
2152 accessor->matrix->n_nonzero_elements());
2153 return accessor->matrix->val[accessor->linear_index];
2158 template <
typename number>
2159 inline const typename Accessor<number, false>::Reference &
2160 Accessor<number, false>::Reference::operator=(
const number n)
const
2163 accessor->matrix->n_nonzero_elements());
2164 accessor->matrix->val[accessor->linear_index] = n;
2170 template <
typename number>
2171 inline const typename Accessor<number, false>::Reference &
2172 Accessor<number, false>::Reference::operator+=(
const number n)
const
2175 accessor->matrix->n_nonzero_elements());
2176 accessor->matrix->val[accessor->linear_index] += n;
2182 template <
typename number>
2183 inline const typename Accessor<number, false>::Reference &
2184 Accessor<number, false>::Reference::operator-=(
const number n)
const
2187 accessor->matrix->n_nonzero_elements());
2188 accessor->matrix->val[accessor->linear_index] -= n;
2194 template <
typename number>
2195 inline const typename Accessor<number, false>::Reference &
2196 Accessor<number, false>::Reference::operator*=(
const number n)
const
2199 accessor->matrix->n_nonzero_elements());
2200 accessor->matrix->val[accessor->linear_index] *= n;
2206 template <
typename number>
2207 inline const typename Accessor<number, false>::Reference &
2208 Accessor<number, false>::Reference::operator/=(
const number n)
const
2211 accessor->matrix->n_nonzero_elements());
2212 accessor->matrix->val[accessor->linear_index] /= n;
2218 template <
typename number>
2220 const std::size_t index)
2227 template <
typename number>
2235 template <
typename number>
2236 inline typename Accessor<number, false>::Reference
2239 return Reference(
this,
true);
2244 template <
typename number>
2245 inline typename Accessor<number, false>::MatrixType &
2246 Accessor<number, false>::get_matrix()
const
2253 template <
typename number,
bool Constness>
2255 const std::size_t index)
2256 : accessor(
matrix, index)
2261 template <
typename number,
bool Constness>
2268 template <
typename number,
bool Constness>
2276 template <
typename number,
bool Constness>
2277 inline const Iterator<number, Constness> &
2278 Iterator<number, Constness>::
2287 template <
typename number,
bool Constness>
2288 inline Iterator<number, Constness> &
2296 template <
typename number,
bool Constness>
2297 inline Iterator<number, Constness>
2306 template <
typename number,
bool Constness>
2314 template <
typename number,
bool Constness>
2315 inline const Accessor<number, Constness> *Iterator<number, Constness>::
2322 template <
typename number,
bool Constness>
2326 return (accessor == other.accessor);
2330 template <
typename number,
bool Constness>
2334 return !(*
this == other);
2338 template <
typename number,
bool Constness>
2342 Assert(&accessor.get_matrix() == &other.accessor.get_matrix(),
2345 return (accessor < other.accessor);
2349 template <
typename number,
bool Constness>
2353 return (other < *
this);
2357 template <
typename number,
bool Constness>
2361 Assert(&accessor.get_matrix() == &other.accessor.get_matrix(),
2364 return (*this)->linear_index - other->linear_index;
2369 template <
typename number,
bool Constness>
2370 inline Iterator<number, Constness>
2384 template <
typename number>
2388 return const_iterator(
this, 0);
2392 template <
typename number>
2396 return const_iterator(
this);
2400 template <
typename number>
2404 return iterator(
this, 0);
2408 template <
typename number>
2412 return iterator(
this, cols->rowstart[cols->rows]);
2416 template <
typename number>
2422 return const_iterator(
this, cols->rowstart[r]);
2427 template <
typename number>
2433 return const_iterator(
this, cols->rowstart[r + 1]);
2438 template <
typename number>
2444 return iterator(
this, cols->rowstart[r]);
2449 template <
typename number>
2455 return iterator(
this, cols->rowstart[r + 1]);
2460 template <
typename number>
2461 template <
class StreamType>
2465 const bool diagonal_first)
const
2470 bool hanging_diagonal =
false;
2473 for (
size_type i = 0; i < cols->rows; ++i)
2475 for (
size_type j = cols->rowstart[i]; j < cols->rowstart[i + 1]; ++j)
2477 if (!diagonal_first && i == cols->colnums[j])
2480 hanging_diagonal =
true;
2484 if (hanging_diagonal && cols->colnums[j] > i)
2487 out <<
' ' << i <<
',' << i <<
':' <<
diagonal;
2489 out <<
'(' << i <<
',' << i <<
") " <<
diagonal
2491 hanging_diagonal =
false;
2494 out <<
' ' << i <<
',' << cols->colnums[j] <<
':' << val[j];
2496 out <<
"(" << i <<
"," << cols->colnums[j] <<
") " << val[j]
2500 if (hanging_diagonal)
2503 out <<
' ' << i <<
',' << i <<
':' <<
diagonal;
2505 out <<
'(' << i <<
',' << i <<
") " <<
diagonal << std::endl;
2506 hanging_diagonal =
false;
2514 template <
typename number>
2523 template <
typename number>
#define DeclExceptionMsg(Exception, defaulttext)
void add(const size_type i, const size_type j, const number value)
bool operator!=(const AlignedVector< T > &lhs, const AlignedVector< T > &rhs)
SparseMatrix< number > & copy_from(const SparseMatrix< somenumber > &source)
@ diagonal
Matrix is diagonal.
void print_as_numpy_arrays(std::ostream &out, const unsigned int precision=9) const
static const bool zero_addition_can_be_elided
SynchronousIterators< Iterators > operator++(SynchronousIterators< Iterators > &a)
typename Accessor< number, Constness >::MatrixType MatrixType
Iterator operator+(const size_type n) const
bool operator<(const SynchronousIterators< Iterators > &a, const SynchronousIterators< Iterators > &b)
void block_read(std::istream &in)
void precondition_SOR(Vector< somenumber > &dst, const Vector< somenumber > &src, const number om=1.) const
void SOR(Vector< somenumber > &v, const number om=1.) const
static ::ExceptionBase & ExcInvalidIndex(int arg1, int arg2)
const Iterator< number, Constness > & operator=(const SparseMatrixIterators::Iterator< number, false > &i)
number diag_element(const size_type i) const
bool operator==(const AlignedVector< T > &lhs, const AlignedVector< T > &rhs)
void print(StreamType &out, const bool across=false, const bool diagonal_first=true) const
MatrixTableIterators::Accessor< TransposeTable< T >, Constness, MatrixTableIterators::Storage::column_major > Accessor
const Accessor< number, Constness > & value_type
SparseMatrix & operator/=(const number factor)
const SparsityPattern & get_sparsity_pattern() const
real_type frobenius_norm() const
#define AssertIndexRange(index, range)
void block_write(std::ostream &out) const
SymmetricTensor< 2, dim, Number > C(const Tensor< 2, dim, Number > &F)
const_iterator end() const
typename numbers::NumberTraits< double >::real_type real_type
static ::ExceptionBase & ExcDivideByZero()
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
somenumber residual(Vector< somenumber > &dst, const Vector< somenumber > &x, const Vector< somenumber > &b) const
std::unique_ptr< number[]> val
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
friend class SparseMatrix
int operator-(const Iterator &p) const
__global__ void set(Number *val, const Number s, const size_type N)
void print_pattern(std::ostream &out, const double threshold=0.) const
somenumber matrix_norm_square(const Vector< somenumber > &v) const
static ::ExceptionBase & ExcNotQuadratic()
void vmult_add(OutVector &dst, const InVector &src) const
bool operator>(const Iterator &) const
void precondition_Jacobi(Vector< somenumber > &dst, const Vector< somenumber > &src, const number omega=1.) const
const Accessor< number, Constness > * operator->() const
const Accessor * accessor
void compress(::VectorOperation::values)
bool operator==(const Iterator &) const
Expression operator>(const Expression &lhs, const Expression &rhs)
static ::ExceptionBase & ExcIteratorRange(int arg1, int arg2)
void PSOR(Vector< somenumber > &v, const std::vector< size_type > &permutation, const std::vector< size_type > &inverse_permutation, 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
MatrixTableIterators::Iterator< TransposeTable< T >, Constness, MatrixTableIterators::Storage::column_major > Iterator
static const size_type invalid_entry
#define AssertIsFinite(number)
void Tvmult(OutVector &dst, const InVector &src) const
void SSOR(Vector< somenumber > &v, const number omega=1.) const
unsigned int global_dof_index
VectorType::value_type * begin(VectorType &V)
virtual void reinit(const SparsityPattern &sparsity)
const number & operator()(const size_type i, const size_type j) const
void SOR_step(Vector< somenumber > &v, const Vector< somenumber > &b, const number om=1.) const
#define DEAL_II_NAMESPACE_OPEN
std::size_t n_nonzero_elements() const
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
@ matrix
Contents is actually a matrix.
std::size_t n_actually_nonzero_elements(const double threshold=0.) const
const SparseMatrix< number > MatrixType
const_iterator begin() const
VectorType::value_type * end(VectorType &V)
T sum(const T &t, const MPI_Comm &mpi_communicator)
static ::ExceptionBase & ExcNotInitialized()
types::global_dof_index size_type
size_type get_row_length(const size_type row) const
SparseMatrix< number > MatrixType
constexpr SymmetricTensor< rank_, dim, typename ProductType< Number, OtherNumber >::type > operator-(const SymmetricTensor< rank_, dim, Number > &left, const SymmetricTensor< rank_, dim, OtherNumber > &right)
static ::ExceptionBase & ExcDifferentSparsityPatterns()
void set(const size_type i, const size_type j, const number value)
real_type linfty_norm() const
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
static ::ExceptionBase & ExcSourceEqualsDestination()
const SparseMatrix< number > & get_matrix() const
const types::global_dof_index * dummy()
virtual ~SparseMatrix() override
static ::ExceptionBase & ExcInternalError()
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 om=1.) const
#define Assert(cond, exc)
SmartPointer< const SparsityPattern, SparseMatrix< number > > cols
void Tmmult(SparseMatrix< numberC > &C, const SparseMatrix< numberB > &B, const Vector< number > &V=Vector< number >(), const bool rebuild_sparsity_pattern=true) const
bool operator<(const Iterator &) const
Accessor< number, Constness > accessor
void TPSOR(Vector< somenumber > &v, const std::vector< size_type > &permutation, const std::vector< size_type > &inverse_permutation, const number om=1.) const
Iterator(MatrixType *matrix, const std::size_t index_within_matrix)
void Jacobi_step(Vector< somenumber > &v, const Vector< somenumber > &b, const number om=1.) const
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
void precondition_TSOR(Vector< somenumber > &dst, const Vector< somenumber > &src, const number om=1.) const
bool operator!=(const Iterator &) const
real_type l1_norm() const
SparseMatrix & operator*=(const number factor)
const Accessor< number, Constness > & operator*() const
#define DEAL_II_NAMESPACE_CLOSE
number el(const size_type i, const size_type j) const
std::size_t memory_consumption() const
constexpr SymmetricTensor< rank_, dim, typename ProductType< Number, OtherNumber >::type > operator+(const SymmetricTensor< rank_, dim, Number > &left, const SymmetricTensor< rank_, dim, OtherNumber > &right)
#define DeclException2(Exception2, type1, type2, outsequence)
void vmult(OutVector &dst, const InVector &src) const
void Tvmult_add(OutVector &dst, const InVector &src) const
SparseMatrix< number > & operator=(const SparseMatrix< number > &)
void TSOR(Vector< somenumber > &v, const number om=1.) const
void SSOR_step(Vector< somenumber > &v, const Vector< somenumber > &b, const number om=1.) const
std::enable_if< std::is_floating_point< T >::value &&std::is_floating_point< U >::value, typename ProductType< std::complex< T >, std::complex< U > >::type >::type operator*(const std::complex< T > &left, const std::complex< U > &right)