16 #ifndef dealii_chunk_sparse_matrix_h 17 # define dealii_chunk_sparse_matrix_h 20 # include <deal.II/base/config.h> 22 # include <deal.II/base/smartpointer.h> 23 # include <deal.II/base/subscriptor.h> 25 # include <deal.II/lac/chunk_sparsity_pattern.h> 26 # include <deal.II/lac/exceptions.h> 27 # include <deal.II/lac/identity_matrix.h> 32 DEAL_II_NAMESPACE_OPEN
34 template <
typename number>
36 template <
typename number>
50 template <
typename number,
bool Constness>
63 template <
typename number,
bool Constness>
95 template <
typename number>
147 template <
typename,
bool>
158 template <
typename number>
191 Reference(
const Accessor *accessor,
const bool dummy);
196 operator number()
const;
202 operator=(
const number n)
const;
208 operator+=(
const number n)
const;
214 operator-=(
const number n)
const;
220 operator*=(
const number n)
const;
226 operator/=(
const number n)
const;
280 template <
typename,
bool>
296 template <
typename number,
bool Constness>
423 template <
typename number>
475 static const bool zero_addition_can_be_elided =
true;
565 operator=(
const double d);
622 n_nonzero_elements()
const;
632 n_actually_nonzero_elements()
const;
643 get_sparsity_pattern()
const;
650 memory_consumption()
const;
682 template <
typename number2>
687 const number2 * values,
688 const bool elide_zero_values =
true,
689 const bool col_indices_are_sorted =
false);
695 operator*=(
const number factor);
701 operator/=(
const number factor);
734 template <
typename somenumber>
755 template <
typename ForwardIterator>
757 copy_from(
const ForwardIterator begin,
const ForwardIterator end);
764 template <
typename somenumber>
779 template <
typename somenumber>
846 number * values)
const;
867 template <
class OutVector,
class InVector>
869 vmult(OutVector &dst,
const InVector &src)
const;
886 template <
class OutVector,
class InVector>
888 Tvmult(OutVector &dst,
const InVector &src)
const;
904 template <
class OutVector,
class InVector>
906 vmult_add(OutVector &dst,
const InVector &src)
const;
923 template <
class OutVector,
class InVector>
925 Tvmult_add(OutVector &dst,
const InVector &src)
const;
942 template <
typename somenumber>
944 matrix_norm_square(
const Vector<somenumber> &v)
const;
949 template <
typename somenumber>
951 matrix_scalar_product(
const Vector<somenumber> &u,
952 const Vector<somenumber> &v)
const;
960 template <
typename somenumber>
962 residual(Vector<somenumber> & dst,
963 const Vector<somenumber> &x,
964 const Vector<somenumber> &b)
const;
997 frobenius_norm()
const;
1009 template <
typename somenumber>
1011 precondition_Jacobi(Vector<somenumber> & dst,
1012 const Vector<somenumber> &src,
1013 const number omega = 1.)
const;
1018 template <
typename somenumber>
1020 precondition_SSOR(Vector<somenumber> & dst,
1021 const Vector<somenumber> &src,
1022 const number om = 1.)
const;
1027 template <
typename somenumber>
1029 precondition_SOR(Vector<somenumber> & dst,
1030 const Vector<somenumber> &src,
1031 const number om = 1.)
const;
1036 template <
typename somenumber>
1038 precondition_TSOR(Vector<somenumber> & dst,
1039 const Vector<somenumber> &src,
1040 const number om = 1.)
const;
1047 template <
typename somenumber>
1049 SSOR(Vector<somenumber> &v,
const number omega = 1.)
const;
1055 template <
typename somenumber>
1057 SOR(Vector<somenumber> &v,
const number om = 1.)
const;
1063 template <
typename somenumber>
1065 TSOR(Vector<somenumber> &v,
const number om = 1.)
const;
1077 template <
typename somenumber>
1079 PSOR(Vector<somenumber> & v,
1080 const std::vector<size_type> &permutation,
1081 const std::vector<size_type> &inverse_permutation,
1082 const number om = 1.)
const;
1094 template <
typename somenumber>
1096 TPSOR(Vector<somenumber> & v,
1097 const std::vector<size_type> &permutation,
1098 const std::vector<size_type> &inverse_permutation,
1099 const number om = 1.)
const;
1105 template <
typename somenumber>
1107 SOR_step(Vector<somenumber> & v,
1108 const Vector<somenumber> &b,
1109 const number om = 1.)
const;
1115 template <
typename somenumber>
1117 TSOR_step(Vector<somenumber> & v,
1118 const Vector<somenumber> &b,
1119 const number om = 1.)
const;
1125 template <
typename somenumber>
1127 SSOR_step(Vector<somenumber> & v,
1128 const Vector<somenumber> &b,
1129 const number om = 1.)
const;
1197 begin(
const unsigned int r)
const;
1214 end(
const unsigned int r)
const;
1231 begin(
const unsigned int r);
1248 end(
const unsigned int r);
1260 print(std::ostream &out)
const;
1283 print_formatted(std::ostream & out,
1284 const unsigned int precision = 3,
1285 const bool scientific =
true,
1286 const unsigned int width = 0,
1287 const char * zero_string =
" ",
1288 const double denominator = 1.)
const;
1296 print_pattern(std::ostream &out,
const double threshold = 0.)
const;
1309 block_write(std::ostream &out)
const;
1328 block_read(std::istream &in);
1341 <<
"You are trying to access the matrix entry with index <" 1342 << arg1 <<
',' << arg2
1343 <<
">, but this entry does not exist in the sparsity pattern" 1346 "The most common cause for this problem is that you used " 1347 "a method to build the sparsity pattern that did not " 1348 "(completely) take into account all of the entries you " 1349 "will later try to write into. An example would be " 1350 "building a sparsity pattern that does not include " 1351 "the entries you will write into due to constraints " 1352 "on degrees of freedom such as hanging nodes or periodic " 1353 "boundary conditions. In such cases, building the " 1354 "sparsity pattern will succeed, but you will get errors " 1355 "such as the current one at one point or other when " 1356 "trying to write into the entries of the matrix.");
1367 <<
"The iterators denote a range of " << arg1
1368 <<
" elements, but the given number of rows was " << arg2);
1373 "You are attempting an operation on two matrices that " 1374 "are the same object, but the operation requires that the " 1375 "two objects are in fact different.");
1392 std::unique_ptr<number[]>
val;
1409 template <
typename somenumber>
1415 template <
typename,
bool>
1417 template <
typename,
bool>
1428 template <
typename number>
1437 template <
typename number>
1447 template <
typename number>
1457 template <
typename number>
1462 const size_type chunk_size = cols->get_chunk_size();
1464 cols->sparsity_pattern(i / chunk_size, j / chunk_size);
1470 return (chunk_index * chunk_size * chunk_size +
1471 (i % chunk_size) * chunk_size + (j % chunk_size));
1476 template <
typename number>
1487 const size_type index = compute_location(i, j);
1489 ExcInvalidIndex(i, j));
1497 template <
typename number>
1507 if (std::abs(value) != 0.)
1509 const size_type index = compute_location(i, j);
1511 ExcInvalidIndex(i, j));
1513 val[index] +=
value;
1519 template <
typename number>
1520 template <
typename number2>
1525 const number2 * values,
1530 for (
size_type col = 0; col < n_cols; ++col)
1531 add(row, col_indices[col], static_cast<number>(values[col]));
1536 template <
typename number>
1543 const size_type chunk_size = cols->get_chunk_size();
1549 number * val_ptr = val.get();
1550 const number *
const end_ptr =
1552 cols->sparsity_pattern.n_nonzero_elements() * chunk_size * chunk_size;
1553 while (val_ptr != end_ptr)
1554 *val_ptr++ *= factor;
1561 template <
typename number>
1569 const number factor_inv = 1. / factor;
1571 const size_type chunk_size = cols->get_chunk_size();
1577 number * val_ptr = val.get();
1578 const number *
const end_ptr =
1580 cols->sparsity_pattern.n_nonzero_elements() * chunk_size * chunk_size;
1582 while (val_ptr != end_ptr)
1583 *val_ptr++ *= factor_inv;
1590 template <
typename number>
1597 ExcInvalidIndex(i, j));
1598 return val[compute_location(i, j)];
1603 template <
typename number>
1608 const size_type index = compute_location(i, j);
1618 template <
typename number>
1628 const size_type chunk_size = cols->get_chunk_size();
1629 return val[cols->sparsity_pattern.rowstart[i / chunk_size] * chunk_size *
1631 (i % chunk_size) * chunk_size + (i % chunk_size)];
1636 template <
typename number>
1637 template <
typename ForwardIterator>
1640 const ForwardIterator end)
1642 Assert(static_cast<size_type>(std::distance(begin, end)) == m(),
1643 ExcIteratorRange(std::distance(begin, end), m()));
1647 using inner_iterator =
1648 typename std::iterator_traits<ForwardIterator>::value_type::const_iterator;
1650 for (ForwardIterator i = begin; i != end; ++i, ++
row)
1652 const inner_iterator end_of_row = i->end();
1653 for (inner_iterator j = i->begin(); j != end_of_row; ++j)
1655 set(row, j->first, j->second);
1666 template <
typename number>
1668 const unsigned int row)
1676 template <
typename number>
1677 inline Accessor<number, true>::Accessor(
const MatrixType *matrix)
1684 template <
typename number>
1685 inline Accessor<number, true>::Accessor(
1688 , matrix(&a.get_matrix())
1693 template <
typename number>
1695 Accessor<number, true>::value()
const 1697 const unsigned int chunk_size =
1698 matrix->get_sparsity_pattern().get_chunk_size();
1699 return matrix->val[reduced_index() * chunk_size * chunk_size +
1700 chunk_row * chunk_size + chunk_col];
1705 template <
typename number>
1706 inline const typename Accessor<number, true>::MatrixType &
1707 Accessor<number, true>::get_matrix()
const 1714 template <
typename number>
1715 inline Accessor<number, false>::Reference::Reference(
const Accessor *accessor,
1717 : accessor(accessor)
1721 template <
typename number>
1722 inline Accessor<number, false>::Reference::operator number()
const 1724 const unsigned int chunk_size =
1725 accessor->matrix->get_sparsity_pattern().get_chunk_size();
1726 return accessor->matrix
1727 ->val[accessor->reduced_index() * chunk_size * chunk_size +
1728 accessor->chunk_row * chunk_size + accessor->chunk_col];
1733 template <
typename number>
1734 inline const typename Accessor<number, false>::Reference &
1735 Accessor<number, false>::Reference::operator=(
const number
n)
const 1737 const unsigned int chunk_size =
1738 accessor->matrix->get_sparsity_pattern().get_chunk_size();
1740 ->val[accessor->reduced_index() * chunk_size * chunk_size +
1741 accessor->chunk_row * chunk_size + accessor->chunk_col] =
n;
1747 template <
typename number>
1748 inline const typename Accessor<number, false>::Reference &
1749 Accessor<number, false>::Reference::operator+=(
const number n)
const 1751 const unsigned int chunk_size =
1752 accessor->matrix->get_sparsity_pattern().get_chunk_size();
1754 ->val[accessor->reduced_index() * chunk_size * chunk_size +
1755 accessor->chunk_row * chunk_size + accessor->chunk_col] +=
n;
1761 template <
typename number>
1762 inline const typename Accessor<number, false>::Reference &
1763 Accessor<number, false>::Reference::operator-=(
const number n)
const 1765 const unsigned int chunk_size =
1766 accessor->matrix->get_sparsity_pattern().get_chunk_size();
1768 ->val[accessor->reduced_index() * chunk_size * chunk_size +
1769 accessor->chunk_row * chunk_size + accessor->chunk_col] -=
n;
1775 template <
typename number>
1776 inline const typename Accessor<number, false>::Reference &
1777 Accessor<number, false>::Reference::operator*=(
const number n)
const 1779 const unsigned int chunk_size =
1780 accessor->matrix->get_sparsity_pattern().get_chunk_size();
1782 ->val[accessor->reduced_index() * chunk_size * chunk_size +
1783 accessor->chunk_row * chunk_size + accessor->chunk_col] *=
n;
1789 template <
typename number>
1790 inline const typename Accessor<number, false>::Reference &
1791 Accessor<number, false>::Reference::operator/=(
const number n)
const 1793 const unsigned int chunk_size =
1794 accessor->matrix->get_sparsity_pattern().get_chunk_size();
1796 ->val[accessor->reduced_index() * chunk_size * chunk_size +
1797 accessor->chunk_row * chunk_size + accessor->chunk_col] /=
n;
1803 template <
typename number>
1804 inline Accessor<number, false>::Accessor(MatrixType * matrix,
1805 const unsigned int row)
1813 template <
typename number>
1814 inline Accessor<number, false>::Accessor(MatrixType *matrix)
1821 template <
typename number>
1822 inline typename Accessor<number, false>::Reference
1823 Accessor<number, false>::value()
const 1825 return Reference(
this,
true);
1830 template <
typename number>
1831 inline typename Accessor<number, false>::MatrixType &
1832 Accessor<number, false>::get_matrix()
const 1839 template <
typename number,
bool Constness>
1840 inline Iterator<number, Constness>::Iterator(MatrixType * matrix,
1841 const unsigned int row)
1842 : accessor(matrix, row)
1847 template <
typename number,
bool Constness>
1848 inline Iterator<number, Constness>::Iterator(MatrixType *matrix)
1854 template <
typename number,
bool Constness>
1855 inline Iterator<number, Constness>::Iterator(
1862 template <
typename number,
bool Constness>
1863 inline Iterator<number, Constness> &
1864 Iterator<number, Constness>::operator++()
1871 template <
typename number,
bool Constness>
1872 inline Iterator<number, Constness>
1873 Iterator<number, Constness>::operator++(
int)
1875 const Iterator iter = *
this;
1881 template <
typename number,
bool Constness>
1889 template <
typename number,
bool Constness>
1890 inline const Accessor<number, Constness> *Iterator<number, Constness>::
1897 template <
typename number,
bool Constness>
1899 Iterator<number, Constness>::operator==(
const Iterator &other)
const 1901 return (accessor == other.accessor);
1905 template <
typename number,
bool Constness>
1907 Iterator<number, Constness>::operator!=(
const Iterator &other)
const 1909 return !(*
this == other);
1913 template <
typename number,
bool Constness>
1915 Iterator<number, Constness>::operator<(
const Iterator &other)
const 1917 Assert(&accessor.get_matrix() == &other.accessor.get_matrix(),
1920 return (accessor < other.accessor);
1924 template <
typename number,
bool Constness>
1926 Iterator<number, Constness>::operator>(
const Iterator &other)
const 1928 return (other < *
this);
1932 template <
typename number,
bool Constness>
1936 Assert(&accessor.get_matrix() == &other.accessor.get_matrix(),
1943 Iterator copy = *
this;
1944 while (copy != other)
1952 Iterator copy = other;
1953 while (copy != *
this)
1964 template <
typename number,
bool Constness>
1965 inline Iterator<number, Constness>
1969 for (
unsigned int i = 0; i <
n; ++i)
1979 template <
typename number>
1987 template <
typename number>
1995 template <
typename number>
2003 template <
typename number>
2011 template <
typename number>
2021 template <
typename number>
2031 template <
typename number>
2041 template <
typename number>
2053 DEAL_II_NAMESPACE_CLOSE
typename numbers::NumberTraits< number >::real_type real_type
number el(const size_type i, const size_type j) const
#define DeclException2(Exception2, type1, type2, outsequence)
static const size_type invalid_entry
ChunkSparseMatrix< number > & copy_from(const ChunkSparseMatrix< somenumber > &source)
#define AssertIndexRange(index, range)
bool operator==(const Accessor &) const
const_iterator end() const
static::ExceptionBase & ExcNotInitialized()
#define AssertThrow(cond, exc)
SymmetricTensor< rank_, dim, typename ProductType< Number, OtherNumber >::type > operator+(const SymmetricTensor< rank_, dim, Number > &left, const SymmetricTensor< rank_, dim, OtherNumber > &right)
static::ExceptionBase & ExcIndexRange(int arg1, int arg2, int arg3)
SymmetricTensor< rank_, dim, Number > operator*(const SymmetricTensor< rank_, dim, Number > &t, const Number &factor)
std::unique_ptr< number[]> val
static::ExceptionBase & ExcDivideByZero()
unsigned long long int global_dof_index
ChunkSparseMatrixIterators::Iterator< number, false > iterator
Accessor< number, Constness > accessor
const Accessor * accessor
const ChunkSparsityPattern & get_sparsity_pattern() const
#define Assert(cond, exc)
number operator()(const size_type i, const size_type j) const
bool operator<(const Accessor &) const
ChunkSparseMatrix & operator*=(const number factor)
#define DeclExceptionMsg(Exception, defaulttext)
#define DeclException0(Exception0)
ChunkSparseMatrixIterators::Iterator< number, true > const_iterator
SymmetricTensor< 2, dim, Number > symmetrize(const Tensor< 2, dim, Number > &t)
SmartPointer< const ChunkSparsityPattern, ChunkSparseMatrix< number > > cols
size_type compute_location(const size_type i, const size_type j) const
Accessor(const ChunkSparsityPattern *matrix, const unsigned int row)
static::ExceptionBase & ExcNotQuadratic()
const ChunkSparseMatrix< number > & get_matrix() const
ChunkSparseMatrix & operator/=(const number factor)
const Accessor< number, Constness > & value_type
number diag_element(const size_type i) const
void set(const size_type i, const size_type j, const number value)
static const size_type invalid_entry
types::global_dof_index size_type
void add(const size_type i, const size_type j, const number value)
typename Accessor< number, Constness >::MatrixType MatrixType
const_iterator begin() const
SymmetricTensor< rank_, dim, typename ProductType< Number, OtherNumber >::type > operator-(const SymmetricTensor< rank_, dim, Number > &left, const SymmetricTensor< rank_, dim, OtherNumber > &right)
#define AssertIsFinite(number)
static::ExceptionBase & ExcInternalError()