16 #ifndef dealii_sparse_matrix_h 17 # define dealii_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/exceptions.h> 26 # include <deal.II/lac/identity_matrix.h> 27 # include <deal.II/lac/sparsity_pattern.h> 28 # include <deal.II/lac/vector_operation.h> 29 # ifdef DEAL_II_WITH_MPI 36 DEAL_II_NAMESPACE_OPEN
38 template <
typename number>
40 template <
typename number>
42 template <
typename Matrix>
44 template <
typename number>
46 # ifdef DEAL_II_WITH_MPI 51 template <
typename Number>
58 # ifdef DEAL_II_WITH_TRILINOS 82 template <
typename number,
bool Constness>
95 template <
typename number,
bool Constness>
127 template <
typename number>
179 template <
typename,
bool>
190 template <
typename number>
225 Reference(
const Accessor *accessor,
const bool dummy);
230 operator number()
const;
236 operator=(
const number n)
const;
242 operator+=(
const number n)
const;
248 operator-=(
const number n)
const;
254 operator*=(
const number n)
const;
260 operator/=(
const number n)
const;
314 template <
typename,
bool>
349 template <
typename number,
bool Constness>
491 template <
typename number>
543 static const bool zero_addition_can_be_elided =
true;
650 operator=(const
double d);
705 get_row_length(const
size_type row) const;
713 n_nonzero_elements() const;
725 n_actually_nonzero_elements(const
double threshold = 0.) const;
736 get_sparsity_pattern() const;
743 memory_consumption() const;
778 template <typename number2>
782 const
bool elide_zero_values = false);
789 template <typename number2>
794 const
bool elide_zero_values = false);
806 template <typename number2>
810 const
std::vector<number2> & values,
811 const
bool elide_zero_values = false);
822 template <typename number2>
827 const number2 * values,
828 const
bool elide_zero_values = false);
852 template <typename number2>
856 const
bool elide_zero_values = true);
863 template <typename number2>
868 const
bool elide_zero_values = true);
879 template <typename number2>
883 const
std::vector<number2> & values,
884 const
bool elide_zero_values = true);
895 template <typename number2>
900 const number2 * values,
901 const
bool elide_zero_values = true,
902 const
bool col_indices_are_sorted = false);
908 operator*=(const number factor);
914 operator/=(const number factor);
947 template <typename somenumber>
967 template <typename ForwardIterator>
969 copy_from(const ForwardIterator begin, const ForwardIterator end);
976 template <typename somenumber>
978 copy_from(const
FullMatrix<somenumber> &matrix);
980 # ifdef DEAL_II_WITH_TRILINOS 1005 template <
typename somenumber>
1092 template <
class OutVector,
class InVector>
1094 vmult(OutVector &dst,
const InVector &src)
const;
1111 template <
class OutVector,
class InVector>
1113 Tvmult(OutVector &dst,
const InVector &src)
const;
1131 template <
class OutVector,
class InVector>
1133 vmult_add(OutVector &dst,
const InVector &src)
const;
1150 template <
class OutVector,
class InVector>
1152 Tvmult_add(OutVector &dst,
const InVector &src)
const;
1171 template <
typename somenumber>
1173 matrix_norm_square(
const Vector<somenumber> &v)
const;
1180 template <
typename somenumber>
1182 matrix_scalar_product(
const Vector<somenumber> &u,
1183 const Vector<somenumber> &v)
const;
1194 template <
typename somenumber>
1196 residual(Vector<somenumber> & dst,
1197 const Vector<somenumber> &x,
1198 const Vector<somenumber> &b)
const;
1235 template <
typename numberB,
typename numberC>
1240 const bool rebuild_sparsity_pattern =
true)
const;
1266 template <
typename numberB,
typename numberC>
1271 const bool rebuild_sparsity_pattern =
true)
const;
1297 linfty_norm()
const;
1304 frobenius_norm()
const;
1316 template <
typename somenumber>
1318 precondition_Jacobi(Vector<somenumber> & dst,
1319 const Vector<somenumber> &src,
1320 const number omega = 1.)
const;
1328 template <
typename somenumber>
1330 precondition_SSOR(Vector<somenumber> & dst,
1331 const Vector<somenumber> & src,
1332 const number omega = 1.,
1333 const std::vector<std::size_t> &pos_right_of_diagonal =
1334 std::vector<std::size_t>())
const;
1339 template <
typename somenumber>
1341 precondition_SOR(Vector<somenumber> & dst,
1342 const Vector<somenumber> &src,
1343 const number om = 1.)
const;
1348 template <
typename somenumber>
1350 precondition_TSOR(Vector<somenumber> & dst,
1351 const Vector<somenumber> &src,
1352 const number om = 1.)
const;
1359 template <
typename somenumber>
1361 SSOR(Vector<somenumber> &v,
const number omega = 1.)
const;
1367 template <
typename somenumber>
1369 SOR(Vector<somenumber> &v,
const number om = 1.)
const;
1375 template <
typename somenumber>
1377 TSOR(Vector<somenumber> &v,
const number om = 1.)
const;
1389 template <
typename somenumber>
1391 PSOR(Vector<somenumber> & v,
1392 const std::vector<size_type> &permutation,
1393 const std::vector<size_type> &inverse_permutation,
1394 const number om = 1.)
const;
1406 template <
typename somenumber>
1408 TPSOR(Vector<somenumber> & v,
1409 const std::vector<size_type> &permutation,
1410 const std::vector<size_type> &inverse_permutation,
1411 const number om = 1.)
const;
1418 template <
typename somenumber>
1420 Jacobi_step(Vector<somenumber> & v,
1421 const Vector<somenumber> &b,
1422 const number om = 1.)
const;
1428 template <
typename somenumber>
1430 SOR_step(Vector<somenumber> & v,
1431 const Vector<somenumber> &b,
1432 const number om = 1.)
const;
1438 template <
typename somenumber>
1440 TSOR_step(Vector<somenumber> & v,
1441 const Vector<somenumber> &b,
1442 const number om = 1.)
const;
1448 template <
typename somenumber>
1450 SSOR_step(Vector<somenumber> & v,
1451 const Vector<somenumber> &b,
1452 const number om = 1.)
const;
1538 template <
class StreamType>
1540 print(StreamType &out,
1541 const bool across =
false,
1542 const bool diagonal_first =
true)
const;
1565 print_formatted(std::ostream & out,
1566 const unsigned int precision = 3,
1567 const bool scientific =
true,
1568 const unsigned int width = 0,
1569 const char * zero_string =
" ",
1570 const double denominator = 1.)
const;
1578 print_pattern(std::ostream &out,
const double threshold = 0.)
const;
1591 block_write(std::ostream &out)
const;
1610 block_read(std::istream &in);
1623 <<
"You are trying to access the matrix entry with index <" 1624 << arg1 <<
',' << arg2
1625 <<
">, but this entry does not exist in the sparsity pattern " 1628 "The most common cause for this problem is that you used " 1629 "a method to build the sparsity pattern that did not " 1630 "(completely) take into account all of the entries you " 1631 "will later try to write into. An example would be " 1632 "building a sparsity pattern that does not include " 1633 "the entries you will write into due to constraints " 1634 "on degrees of freedom such as hanging nodes or periodic " 1635 "boundary conditions. In such cases, building the " 1636 "sparsity pattern will succeed, but you will get errors " 1637 "such as the current one at one point or other when " 1638 "trying to write into the entries of the matrix.");
1643 "When copying one sparse matrix into another, " 1644 "or when adding one sparse matrix to another, " 1645 "both matrices need to refer to the same " 1646 "sparsity pattern.");
1653 <<
"The iterators denote a range of " << arg1
1654 <<
" elements, but the given number of rows was " << arg2);
1659 "You are attempting an operation on two matrices that " 1660 "are the same object, but the operation requires that the " 1661 "two objects are in fact different.");
1700 std::unique_ptr<number[]>
val;
1711 template <
typename somenumber>
1713 template <
typename somenumber>
1727 template <
typename,
bool>
1729 template <
typename,
bool>
1732 # ifdef DEAL_II_WITH_MPI 1736 template <
typename Number>
1749 template <
typename number>
1758 template <
typename number>
1768 template <
typename number>
1783 ExcInvalidIndex(i, j));
1792 template <
typename number>
1793 template <
typename number2>
1797 const bool elide_zero_values)
1799 Assert(indices.size() == values.
m(),
1803 for (
size_type i = 0; i < indices.size(); ++i)
1813 template <
typename number>
1814 template <
typename number2>
1817 const std::vector<size_type> &col_indices,
1819 const bool elide_zero_values)
1821 Assert(row_indices.size() == values.
m(),
1823 Assert(col_indices.size() == values.
n(),
1826 for (
size_type i = 0; i < row_indices.size(); ++i)
1836 template <
typename number>
1837 template <
typename number2>
1840 const std::vector<size_type> &col_indices,
1841 const std::vector<number2> & values,
1842 const bool elide_zero_values)
1844 Assert(col_indices.size() == values.size(),
1856 template <
typename number>
1864 if (value == number())
1867 const size_type index = cols->operator()(i, j);
1874 ExcInvalidIndex(i, j));
1883 template <
typename number>
1884 template <
typename number2>
1888 const bool elide_zero_values)
1890 Assert(indices.size() == values.
m(),
1894 for (
size_type i = 0; i < indices.size(); ++i)
1904 template <
typename number>
1905 template <
typename number2>
1908 const std::vector<size_type> &col_indices,
1910 const bool elide_zero_values)
1912 Assert(row_indices.size() == values.
m(),
1914 Assert(col_indices.size() == values.
n(),
1917 for (
size_type i = 0; i < row_indices.size(); ++i)
1927 template <
typename number>
1928 template <
typename number2>
1931 const std::vector<size_type> &col_indices,
1932 const std::vector<number2> & values,
1933 const bool elide_zero_values)
1935 Assert(col_indices.size() == values.size(),
1947 template <
typename number>
1954 number * val_ptr = val.get();
1955 const number *
const end_ptr = val.get() + cols->n_nonzero_elements();
1957 while (val_ptr != end_ptr)
1958 *val_ptr++ *= factor;
1965 template <
typename number>
1973 const number factor_inv = number(1.) / factor;
1975 number * val_ptr = val.get();
1976 const number *
const end_ptr = val.get() + cols->n_nonzero_elements();
1978 while (val_ptr != end_ptr)
1979 *val_ptr++ *= factor_inv;
1986 template <
typename number>
1987 inline const number &
1992 ExcInvalidIndex(i, j));
1993 return val[cols->operator()(i, j)];
1998 template <
typename number>
2004 ExcInvalidIndex(i, j));
2005 return val[cols->operator()(i, j)];
2010 template <
typename number>
2015 const size_type index = cols->operator()(i, j);
2025 template <
typename number>
2035 return val[cols->rowstart[i]];
2040 template <
typename number>
2050 return val[cols->rowstart[i]];
2055 template <
typename number>
2056 template <
typename ForwardIterator>
2059 const ForwardIterator end)
2061 Assert(static_cast<size_type>(std::distance(begin, end)) == m(),
2062 ExcIteratorRange(std::distance(begin, end), m()));
2066 using inner_iterator =
2067 typename std::iterator_traits<ForwardIterator>::value_type::const_iterator;
2069 for (ForwardIterator i = begin; i != end; ++i, ++
row)
2071 const inner_iterator end_of_row = i->end();
2072 for (inner_iterator j = i->begin(); j != end_of_row; ++j)
2074 set(row, j->first, j->second);
2084 template <
typename number>
2085 inline Accessor<number, true>::Accessor(
const MatrixType *matrix,
2086 const std::size_t index_within_matrix)
2088 index_within_matrix)
2094 template <
typename number>
2095 inline Accessor<number, true>::Accessor(
const MatrixType *matrix)
2102 template <
typename number>
2103 inline Accessor<number, true>::Accessor(
2106 , matrix(&a.get_matrix())
2111 template <
typename number>
2113 Accessor<number, true>::value()
const 2116 return matrix->val[index_within_sparsity];
2121 template <
typename number>
2122 inline const typename Accessor<number, true>::MatrixType &
2123 Accessor<number, true>::get_matrix()
const 2130 template <
typename number>
2131 inline Accessor<number, false>::Reference::Reference(
const Accessor *accessor,
2133 : accessor(accessor)
2137 template <
typename number>
2138 inline Accessor<number, false>::Reference::operator number()
const 2141 accessor->matrix->n_nonzero_elements());
2142 return accessor->matrix->val[accessor->index_within_sparsity];
2147 template <
typename number>
2148 inline const typename Accessor<number, false>::Reference &
2149 Accessor<number, false>::Reference::operator=(
const number
n)
const 2152 accessor->matrix->n_nonzero_elements());
2153 accessor->matrix->val[accessor->index_within_sparsity] =
n;
2159 template <
typename number>
2160 inline const typename Accessor<number, false>::Reference &
2161 Accessor<number, false>::Reference::operator+=(
const number n)
const 2164 accessor->matrix->n_nonzero_elements());
2165 accessor->matrix->val[accessor->index_within_sparsity] +=
n;
2171 template <
typename number>
2172 inline const typename Accessor<number, false>::Reference &
2173 Accessor<number, false>::Reference::operator-=(
const number n)
const 2176 accessor->matrix->n_nonzero_elements());
2177 accessor->matrix->val[accessor->index_within_sparsity] -=
n;
2183 template <
typename number>
2184 inline const typename Accessor<number, false>::Reference &
2185 Accessor<number, false>::Reference::operator*=(
const number n)
const 2188 accessor->matrix->n_nonzero_elements());
2189 accessor->matrix->val[accessor->index_within_sparsity] *=
n;
2195 template <
typename number>
2196 inline const typename Accessor<number, false>::Reference &
2197 Accessor<number, false>::Reference::operator/=(
const number n)
const 2200 accessor->matrix->n_nonzero_elements());
2201 accessor->matrix->val[accessor->index_within_sparsity] /=
n;
2207 template <
typename number>
2208 inline Accessor<number, false>::Accessor(MatrixType * matrix,
2209 const std::size_t index)
2216 template <
typename number>
2217 inline Accessor<number, false>::Accessor(MatrixType *matrix)
2224 template <
typename number>
2225 inline typename Accessor<number, false>::Reference
2226 Accessor<number, false>::value()
const 2228 return Reference(
this,
true);
2233 template <
typename number>
2234 inline typename Accessor<number, false>::MatrixType &
2235 Accessor<number, false>::get_matrix()
const 2242 template <
typename number,
bool Constness>
2243 inline Iterator<number, Constness>::Iterator(MatrixType * matrix,
2244 const std::size_t index)
2245 : accessor(matrix, index)
2250 template <
typename number,
bool Constness>
2251 inline Iterator<number, Constness>::Iterator(MatrixType *matrix)
2257 template <
typename number,
bool Constness>
2258 inline Iterator<number, Constness>::Iterator(
2265 template <
typename number,
bool Constness>
2266 inline Iterator<number, Constness> &
2267 Iterator<number, Constness>::operator++()
2274 template <
typename number,
bool Constness>
2275 inline Iterator<number, Constness>
2276 Iterator<number, Constness>::operator++(
int)
2278 const Iterator iter = *
this;
2284 template <
typename number,
bool Constness>
2292 template <
typename number,
bool Constness>
2293 inline const Accessor<number, Constness> *Iterator<number, Constness>::
2300 template <
typename number,
bool Constness>
2302 Iterator<number, Constness>::operator==(
const Iterator &other)
const 2304 return (accessor == other.accessor);
2308 template <
typename number,
bool Constness>
2310 Iterator<number, Constness>::operator!=(
const Iterator &other)
const 2312 return !(*
this == other);
2316 template <
typename number,
bool Constness>
2318 Iterator<number, Constness>::operator<(
const Iterator &other)
const 2320 Assert(&accessor.get_matrix() == &other.accessor.get_matrix(),
2323 return (accessor < other.accessor);
2327 template <
typename number,
bool Constness>
2329 Iterator<number, Constness>::operator>(
const Iterator &other)
const 2331 return (other < *
this);
2335 template <
typename number,
bool Constness>
2339 Assert(&accessor.get_matrix() == &other.accessor.get_matrix(),
2342 return (*this)->index_within_sparsity - other->index_within_sparsity;
2347 template <
typename number,
bool Constness>
2348 inline Iterator<number, Constness>
2362 template <
typename number>
2370 template <
typename number>
2378 template <
typename number>
2386 template <
typename number>
2394 template <
typename number>
2405 template <
typename number>
2416 template <
typename number>
2427 template <
typename number>
2438 template <
typename number>
2439 template <
class StreamType>
2443 const bool diagonal_first)
const 2448 bool hanging_diagonal =
false;
2449 number diagonal = number();
2458 hanging_diagonal =
true;
2465 out <<
' ' << i <<
',' << i <<
':' << diagonal;
2467 out <<
'(' << i <<
',' << i <<
") " << diagonal
2469 hanging_diagonal =
false;
2478 if (hanging_diagonal)
2481 out <<
' ' << i <<
',' << i <<
':' << diagonal;
2483 out <<
'(' << i <<
',' << i <<
") " << diagonal << std::endl;
2484 hanging_diagonal =
false;
2492 template <
typename number>
2501 template <
typename number>
2513 DEAL_II_NAMESPACE_CLOSE
typename Accessor< number, Constness >::MatrixType MatrixType
SparseMatrix & operator/=(const number factor)
#define DeclException2(Exception2, type1, type2, outsequence)
types::global_dof_index size_type
typename numbers::NumberTraits< number >::real_type real_type
static const size_type invalid_entry
std::unique_ptr< number[]> val
void print(StreamType &out, const bool across=false, const bool diagonal_first=true) const
void set(const size_type i, const size_type j, const number value)
#define AssertIndexRange(index, range)
static::ExceptionBase & ExcNotInitialized()
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)
static::ExceptionBase & ExcDivideByZero()
unsigned long long int global_dof_index
const_iterator begin() const
number el(const size_type i, const size_type j) const
std::unique_ptr< size_type[]> colnums
number diag_element(const size_type i) const
T sum(const T &t, const MPI_Comm &mpi_communicator)
#define Assert(cond, exc)
static::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
#define DeclExceptionMsg(Exception, defaulttext)
SparseMatrixIterators::Iterator< number, false > iterator
SmartPointer< const SparsityPattern, SparseMatrix< number > > cols
Accessor< number, Constness > accessor
SymmetricTensor< 2, dim, Number > symmetrize(const Tensor< 2, dim, Number > &t)
void add(const size_type i, const size_type j, const number value)
const_iterator end() const
types::global_dof_index size_type
static::ExceptionBase & ExcNotQuadratic()
const Accessor * accessor
const SparsityPattern & get_sparsity_pattern() const
std::unique_ptr< std::size_t[]> rowstart
SparseMatrix< number > & copy_from(const SparseMatrix< somenumber > &source)
const number & operator()(const size_type i, const size_type j) const
SparseMatrixIterators::Iterator< number, true > const_iterator
SymmetricTensor< rank_, dim, typename ProductType< Number, OtherNumber >::type > operator-(const SymmetricTensor< rank_, dim, Number > &left, const SymmetricTensor< rank_, dim, OtherNumber > &right)
#define AssertIsFinite(number)
const Accessor< number, Constness > & value_type
SparseMatrix & operator*=(const number factor)
static::ExceptionBase & ExcInternalError()