17 #ifndef __deal2__sparse_matrix_h
18 #define __deal2__sparse_matrix_h
21 #include <deal.II/base/config.h>
22 #include <deal.II/base/subscriptor.h>
23 #include <deal.II/base/smartpointer.h>
24 #include <deal.II/lac/sparsity_pattern.h>
25 #include <deal.II/lac/identity_matrix.h>
26 #include <deal.II/lac/exceptions.h>
28 #include <deal.II/lac/vector.h>
32 template <
typename number>
class Vector;
37 #ifdef DEAL_II_WITH_TRILINOS
61 template <
typename number,
bool Constness>
74 template <
typename number,
bool Constness>
103 template <
typename number>
127 const std::size_t index_within_matrix);
142 number
value()
const;
164 template <
typename,
bool>
175 template <
typename number>
210 Reference (
const Accessor *accessor,
216 operator number ()
const;
221 const Reference &operator = (
const number n)
const;
226 const Reference &operator += (
const number n)
const;
231 const Reference &operator -= (
const number n)
const;
236 const Reference &operator *= (
const number n)
const;
241 const Reference &operator /= (
const number n)
const;
263 const size_type
index);
269 const std::size_t
index);
279 Reference
value()
const;
301 template <
typename,
bool>
342 template <
typename number,
bool Constness>
376 const std::size_t index_within_matrix);
492 template <
typename number>
658 virtual void clear ();
674 size_type
m ()
const;
680 size_type
n ()
const;
724 void compress (::VectorOperation::values);
736 void set (
const size_type i,
755 template <
typename number2>
756 void set (
const std::vector<size_type> &indices,
758 const bool elide_zero_values =
false);
765 template <
typename number2>
766 void set (
const std::vector<size_type> &row_indices,
767 const std::vector<size_type> &col_indices,
769 const bool elide_zero_values =
false);
781 template <
typename number2>
782 void set (
const size_type row,
783 const std::vector<size_type> &col_indices,
784 const std::vector<number2> &values,
785 const bool elide_zero_values =
false);
796 template <
typename number2>
797 void set (
const size_type row,
798 const size_type n_cols,
799 const size_type *col_indices,
800 const number2 *values,
801 const bool elide_zero_values =
false);
808 void add (
const size_type i,
826 template <
typename number2>
827 void add (
const std::vector<size_type> &indices,
829 const bool elide_zero_values =
true);
836 template <
typename number2>
837 void add (
const std::vector<size_type> &row_indices,
838 const std::vector<size_type> &col_indices,
840 const bool elide_zero_values =
true);
851 template <
typename number2>
852 void add (
const size_type row,
853 const std::vector<size_type> &col_indices,
854 const std::vector<number2> &values,
855 const bool elide_zero_values =
true);
866 template <
typename number2>
867 void add (
const size_type row,
868 const size_type n_cols,
869 const size_type *col_indices,
870 const number2 *values,
871 const bool elide_zero_values =
true,
872 const bool col_indices_are_sorted =
false);
911 template <
typename somenumber>
931 template <
typename ForwardIterator>
933 const ForwardIterator
end);
940 template <
typename somenumber>
943 #ifdef DEAL_II_WITH_TRILINOS
968 template <
typename somenumber>
969 void add (
const number factor,
992 const size_type j)
const;
1006 number
el (
const size_type i,
1007 const size_type j)
const;
1078 template <
class OutVector,
class InVector>
1079 void vmult (OutVector &dst,
1080 const InVector &src)
const;
1096 template <
class OutVector,
class InVector>
1097 void Tvmult (OutVector &dst,
1098 const InVector &src)
const;
1113 template <
class OutVector,
class InVector>
1115 const InVector &src)
const;
1131 template <
class OutVector,
class InVector>
1133 const InVector &src)
const;
1150 template <
typename somenumber>
1156 template <
typename somenumber>
1167 template <
typename somenumber>
1195 template <
typename numberB,
typename numberC>
1199 const bool rebuild_sparsity_pattern =
true)
const;
1225 template <
typename numberB,
typename numberC>
1229 const bool rebuild_sparsity_pattern =
true)
const;
1272 template <
typename somenumber>
1275 const number omega = 1.)
const;
1283 template <
typename somenumber>
1286 const number omega = 1.,
1287 const std::vector<std::size_t> &pos_right_of_diagonal=std::vector<std::size_t>())
const;
1292 template <
typename somenumber>
1295 const number om = 1.)
const;
1300 template <
typename somenumber>
1303 const number om = 1.)
const;
1310 template <
typename somenumber>
1312 const number omega = 1.)
const;
1318 template <
typename somenumber>
1320 const number om = 1.)
const;
1326 template <
typename somenumber>
1328 const number om = 1.)
const;
1340 template <
typename somenumber>
1342 const std::vector<size_type> &permutation,
1343 const std::vector<size_type> &inverse_permutation,
1344 const number om = 1.)
const;
1356 template <
typename somenumber>
1358 const std::vector<size_type> &permutation,
1359 const std::vector<size_type> &inverse_permutation,
1360 const number om = 1.)
const;
1367 template <
typename somenumber>
1370 const number om = 1.)
const;
1376 template <
typename somenumber>
1379 const number om = 1.)
const;
1385 template <
typename somenumber>
1388 const number om = 1.)
const;
1394 template <
typename somenumber>
1397 const number om = 1.)
const;
1498 template <
class STREAM>
1499 void print (STREAM &out,
1500 const bool across =
false,
1501 const bool diagonal_first =
true)
const;
1524 const unsigned int precision = 3,
1525 const bool scientific =
true,
1526 const unsigned int width = 0,
1527 const char *zero_string =
" ",
1528 const double denominator = 1.)
const;
1536 const double threshold = 0.)
const;
1576 <<
"The entry with index <" << arg1 <<
',' << arg2
1577 <<
"> does not exist.");
1583 <<
"The index " << arg1 <<
" is not in the allowed range.");
1593 <<
"The iterators denote a range of " << arg1
1594 <<
" elements, but the given number of rows was " << arg2);
1645 template <
typename somenumber>
friend class SparseMatrix;
1647 template <
typename>
friend class SparseILU;
1667 template <
typename number>
1676 template <
typename number>
1686 template <
typename number>
1695 const size_type
index = cols->operator()(i, j);
1703 ExcInvalidIndex(i, j));
1712 template <
typename number>
1713 template <
typename number2>
1718 const bool elide_zero_values)
1720 Assert (indices.size() == values.
m(),
1722 Assert (values.
m() == values.
n(), ExcNotQuadratic());
1724 for (size_type i=0; i<indices.size(); ++i)
1725 set (indices[i], indices.size(), &indices[0], &values(i,0),
1731 template <
typename number>
1732 template <
typename number2>
1736 const std::vector<size_type> &col_indices,
1738 const bool elide_zero_values)
1740 Assert (row_indices.size() == values.
m(),
1742 Assert (col_indices.size() == values.
n(),
1745 for (size_type i=0; i<row_indices.size(); ++i)
1746 set (row_indices[i], col_indices.size(), &col_indices[0], &values(i,0),
1752 template <
typename number>
1753 template <
typename number2>
1757 const std::vector<size_type> &col_indices,
1758 const std::vector<number2> &values,
1759 const bool elide_zero_values)
1761 Assert (col_indices.size() == values.size(),
1764 set (row, col_indices.size(), &col_indices[0], &values[0],
1770 template <
typename number>
1782 const size_type index = cols->operator()(i, j);
1790 ExcInvalidIndex(i, j));
1799 template <
typename number>
1800 template <
typename number2>
1805 const bool elide_zero_values)
1807 Assert (indices.size() == values.
m(),
1809 Assert (values.
m() == values.
n(), ExcNotQuadratic());
1811 for (size_type i=0; i<indices.size(); ++i)
1812 add (indices[i], indices.size(), &indices[0], &values(i,0),
1818 template <
typename number>
1819 template <
typename number2>
1823 const std::vector<size_type> &col_indices,
1825 const bool elide_zero_values)
1827 Assert (row_indices.size() == values.
m(),
1829 Assert (col_indices.size() == values.
n(),
1832 for (size_type i=0; i<row_indices.size(); ++i)
1833 add (row_indices[i], col_indices.size(), &col_indices[0], &values(i,0),
1839 template <
typename number>
1840 template <
typename number2>
1844 const std::vector<size_type> &col_indices,
1845 const std::vector<number2> &values,
1846 const bool elide_zero_values)
1848 Assert (col_indices.size() == values.size(),
1851 add (row, col_indices.size(), &col_indices[0], &values[0],
1857 template <
typename number>
1865 number *val_ptr = &val[0];
1866 const number *
const end_ptr = &val[cols->n_nonzero_elements()];
1868 while (val_ptr != end_ptr)
1869 *val_ptr++ *= factor;
1876 template <
typename number>
1885 const number factor_inv = 1. / factor;
1887 number *val_ptr = &val[0];
1888 const number *
const end_ptr = &val[cols->n_nonzero_elements()];
1890 while (val_ptr != end_ptr)
1891 *val_ptr++ *= factor_inv;
1898 template <
typename number>
1901 const size_type j)
const
1905 ExcInvalidIndex(i,j));
1906 return val[cols->operator()(i,j)];
1911 template <
typename number>
1914 const size_type j)
const
1917 const size_type index = cols->operator()(i,j);
1927 template <
typename number>
1932 Assert (m() == n(), ExcNotQuadratic());
1933 Assert (i<m(), ExcInvalidIndex1(i));
1937 return val[cols->rowstart[i]];
1942 template <
typename number>
1947 Assert (m() == n(), ExcNotQuadratic());
1948 Assert (i<m(), ExcInvalidIndex1(i));
1952 return val[cols->rowstart[i]];
1957 template <
typename number>
1961 const size_type index)
const
1964 Assert(index<cols->row_length(row),
1968 return val[cols->rowstart[
row]+
index];
1973 template <
typename number>
1978 Assert (j < cols->n_nonzero_elements(),
1986 template <
typename number>
1991 Assert (j < cols->n_nonzero_elements(),
1999 template <
typename number>
2000 template <
typename ForwardIterator>
2003 const ForwardIterator end)
2005 Assert (static_cast<size_type>(std::distance (begin, end)) == m(),
2006 ExcIteratorRange (std::distance (begin, end), m()));
2010 typedef typename std::iterator_traits<ForwardIterator>::value_type::const_iterator inner_iterator;
2012 for (ForwardIterator i=begin; i!=end; ++i, ++
row)
2014 const inner_iterator end_of_row = i->end();
2015 for (inner_iterator j=i->begin(); j!=end_of_row; ++j)
2017 set (row, j->first, j->second);
2027 template <
typename number>
2030 Accessor (
const MatrixType *matrix,
2031 const size_type row,
2032 const size_type index)
2041 template <
typename number>
2043 Accessor<number,true>::
2044 Accessor (
const MatrixType *matrix,
2045 const std::size_t index_within_matrix)
2048 index_within_matrix),
2054 template <
typename number>
2056 Accessor<number,true>::
2057 Accessor (
const MatrixType *matrix)
2065 template <
typename number>
2067 Accessor<number,true>::
2071 matrix (&a.get_matrix())
2076 template <
typename number>
2079 Accessor<number, true>::value ()
const
2082 return matrix->val[index_within_sparsity];
2087 template <
typename number>
2089 typename Accessor<number, true>::MatrixType &
2090 Accessor<number, true>::get_matrix ()
const
2097 template <
typename number>
2099 Accessor<number, false>::Reference::Reference (
2100 const Accessor *accessor,
2107 template <
typename number>
2109 Accessor<number, false>::Reference::operator number()
const
2111 AssertIndexRange(accessor->index_within_sparsity, accessor->matrix->n_nonzero_elements());
2112 return accessor->matrix->val[accessor->index_within_sparsity];
2117 template <
typename number>
2119 const typename Accessor<number, false>::Reference &
2120 Accessor<number, false>::Reference::operator = (
const number n)
const
2122 AssertIndexRange(accessor->index_within_sparsity, accessor->matrix->n_nonzero_elements());
2123 accessor->matrix->val[accessor->index_within_sparsity] = n;
2129 template <
typename number>
2131 const typename Accessor<number, false>::Reference &
2132 Accessor<number, false>::Reference::operator += (
const number n)
const
2134 AssertIndexRange(accessor->index_within_sparsity, accessor->matrix->n_nonzero_elements());
2135 accessor->matrix->val[accessor->index_within_sparsity] += n;
2141 template <
typename number>
2143 const typename Accessor<number, false>::Reference &
2144 Accessor<number, false>::Reference::operator -= (
const number n)
const
2146 AssertIndexRange(accessor->index_within_sparsity, accessor->matrix->n_nonzero_elements());
2147 accessor->matrix->val[accessor->index_within_sparsity] -= n;
2153 template <
typename number>
2155 const typename Accessor<number, false>::Reference &
2156 Accessor<number, false>::Reference::operator *= (
const number n)
const
2158 AssertIndexRange(accessor->index_within_sparsity, accessor->matrix->n_nonzero_elements());
2159 accessor->matrix->val[accessor->index_within_sparsity] *= n;
2165 template <
typename number>
2167 const typename Accessor<number, false>::Reference &
2168 Accessor<number, false>::Reference::operator /= (
const number n)
const
2170 AssertIndexRange(accessor->index_within_sparsity, accessor->matrix->n_nonzero_elements());
2171 accessor->matrix->val[accessor->index_within_sparsity] /= n;
2177 template <
typename number>
2179 Accessor<number,false>::
2180 Accessor (MatrixType *matrix,
2181 const size_type row,
2182 const size_type index)
2191 template <
typename number>
2193 Accessor<number,false>::
2194 Accessor (MatrixType *matrix,
2195 const std::size_t index)
2204 template <
typename number>
2206 Accessor<number,false>::
2207 Accessor (MatrixType *matrix)
2215 template <
typename number>
2217 typename Accessor<number, false>::Reference
2218 Accessor<number, false>::value()
const
2220 return Reference(
this,
true);
2226 template <
typename number>
2228 typename Accessor<number, false>::MatrixType &
2229 Accessor<number, false>::get_matrix ()
const
2236 template <
typename number,
bool Constness>
2238 Iterator<number, Constness>::
2239 Iterator (MatrixType *matrix,
2243 accessor(matrix, r, i)
2248 template <
typename number,
bool Constness>
2250 Iterator<number, Constness>::
2251 Iterator (MatrixType *matrix,
2252 const std::size_t index)
2254 accessor(matrix, index)
2259 template <
typename number,
bool Constness>
2261 Iterator<number, Constness>::
2262 Iterator (MatrixType *matrix)
2269 template <
typename number,
bool Constness>
2271 Iterator<number, Constness>::
2279 template <
typename number,
bool Constness>
2281 Iterator<number, Constness> &
2282 Iterator<number,Constness>::operator++ ()
2284 accessor.advance ();
2289 template <
typename number,
bool Constness>
2291 Iterator<number,Constness>
2292 Iterator<number,Constness>::operator++ (
int)
2294 const Iterator iter = *
this;
2295 accessor.advance ();
2300 template <
typename number,
bool Constness>
2302 const Accessor<number,Constness> &
2303 Iterator<number,Constness>::operator* ()
const
2309 template <
typename number,
bool Constness>
2311 const Accessor<number,Constness> *
2312 Iterator<number,Constness>::operator-> ()
const
2318 template <
typename number,
bool Constness>
2321 Iterator<number,Constness>::
2322 operator == (
const Iterator &other)
const
2324 return (accessor == other.accessor);
2328 template <
typename number,
bool Constness>
2331 Iterator<number,Constness>::
2332 operator != (
const Iterator &other)
const
2334 return ! (*
this == other);
2338 template <
typename number,
bool Constness>
2341 Iterator<number,Constness>::
2342 operator < (
const Iterator &other)
const
2344 Assert (&accessor.get_matrix() == &other.accessor.get_matrix(),
2347 return (accessor < other.accessor);
2351 template <
typename number,
bool Constness>
2354 Iterator<number,Constness>::
2355 operator > (
const Iterator &other)
const
2357 return (other < *
this);
2361 template <
typename number,
bool Constness>
2364 Iterator<number,Constness>::
2365 operator - (
const Iterator &other)
const
2367 Assert (&accessor.get_matrix() == &other.accessor.get_matrix(),
2370 return (*this)->index_within_sparsity - other->index_within_sparsity;
2375 template <
typename number,
bool Constness>
2377 Iterator<number,Constness>
2378 Iterator<number,Constness>::
2379 operator + (
const size_type n)
const
2382 for (size_type i=0; i<n; ++i)
2392 template <
typename number>
2397 return const_iterator(
this, 0);
2401 template <
typename number>
2406 return const_iterator(
this);
2410 template <
typename number>
2415 return iterator (
this, 0);
2419 template <
typename number>
2424 return iterator(
this, cols->rowstart[cols->rows]);
2428 template <
typename number>
2435 return const_iterator(
this, cols->rowstart[r]);
2440 template <
typename number>
2447 return const_iterator(
this, cols->rowstart[r+1]);
2452 template <
typename number>
2459 return iterator(
this, cols->rowstart[r]);
2464 template <
typename number>
2471 return iterator(
this, cols->rowstart[r+1]);
2476 template <
typename number>
2477 template <
class STREAM>
2481 const bool diagonal_first)
const
2486 bool hanging_diagonal =
false;
2487 number diagonal = number();
2489 for (size_type i=0; i<cols->rows; ++i)
2491 for (size_type j=cols->rowstart[i]; j<cols->rowstart[i+1]; ++j)
2493 if (!diagonal_first && i == cols->colnums[j])
2496 hanging_diagonal =
true;
2500 if (hanging_diagonal && cols->colnums[j]>i)
2503 out <<
' ' << i <<
',' << i <<
':' << diagonal;
2505 out <<
'(' << i <<
',' << i <<
") " << diagonal << std::endl;
2506 hanging_diagonal =
false;
2509 out <<
' ' << i <<
',' << cols->colnums[j] <<
':' << val[j];
2511 out <<
"(" << i <<
"," << cols->colnums[j] <<
") " << val[j] << std::endl;
2514 if (hanging_diagonal)
2517 out <<
' ' << i <<
',' << i <<
':' << diagonal;
2519 out <<
'(' << i <<
',' << i <<
") " << diagonal << std::endl;
2520 hanging_diagonal =
false;
2528 template <
typename number>
2539 template <
typename number>
2553 DEAL_II_NAMESPACE_CLOSE
bool operator>(const Iterator &) const
somenumber matrix_norm_square(const Vector< somenumber > &v) const
number raw_entry(const size_type row, const size_type index) const DEAL_II_DEPRECATED
SparseMatrix & operator/=(const number factor)
void precondition_Jacobi(Vector< somenumber > &dst, const Vector< somenumber > &src, const number omega=1.) const
bool operator!=(const Iterator &) const
void Tvmult(OutVector &dst, const InVector &src) const
void vmult(OutVector &dst, const InVector &src) const
void PSOR(Vector< somenumber > &v, const std::vector< size_type > &permutation, const std::vector< size_type > &inverse_permutation, const number om=1.) const
void Tvmult_add(OutVector &dst, const InVector &src) const
numbers::NumberTraits< number >::real_type real_type
void mmult(SparseMatrix< numberC > &C, const SparseMatrix< numberB > &B, const Vector< number > &V=Vector< number >(), const bool rebuild_sparsity_pattern=true) const
Accessor< number, Constness >::MatrixType MatrixType
number operator()(const size_type i, const size_type j) const
Iterator operator+(const size_type n) const
void SSOR(Vector< somenumber > &v, const number omega=1.) const
size_type get_row_length(const size_type row) const
void set(const size_type i, const size_type j, const number value)
DeclException1(ExcInvalidIndex1, int,<< "The index "<< arg1<< " is not in the allowed range.")
#define AssertIndexRange(index, range)
void block_read(std::istream &in)
DeclException2(ExcInvalidIndex, int, int,<< "The entry with index <"<< arg1<< ','<< arg2<< "> does not exist.")
SparseMatrixIterators::Iterator< number, false > iterator
size_type n_actually_nonzero_elements(const double threshold=0.) const
void Tmmult(SparseMatrix< numberC > &C, const SparseMatrix< numberB > &B, const Vector< number > &V=Vector< number >(), const bool rebuild_sparsity_pattern=true) const
void SOR(Vector< somenumber > &v, const number om=1.) const
bool is_finite(const double x)
const_iterator begin() const
number el(const size_type i, const size_type j) const
void print_pattern(std::ostream &out, const double threshold=0.) const
real_type l1_norm() const
void print(STREAM &out, const bool across=false, const bool diagonal_first=true) const
void precondition_TSOR(Vector< somenumber > &dst, const Vector< somenumber > &src, const number om=1.) const
SparseMatrix< number > & operator=(const SparseMatrix< number > &)
number diag_element(const size_type i) const
unsigned int global_dof_index
const SparseMatrix< number > & get_matrix() const
#define Assert(cond, exc)
std::size_t memory_consumption() const
types::global_dof_index size_type
const SparsityPattern & get_sparsity_pattern() 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
const Accessor< number, Constness > & operator*() const
const SparseMatrix< number > MatrixType
void precondition_SOR(Vector< somenumber > &dst, const Vector< somenumber > &src, const number om=1.) const
::ExceptionBase & ExcIndexRange(int arg1, int arg2, int arg3)
size_type n_nonzero_elements() const
const Accessor< number, Constness > * operator->() const
BlockCompressedSparsityPattern CompressedBlockSparsityPattern DEAL_II_DEPRECATED
somenumber residual(Vector< somenumber > &dst, const Vector< somenumber > &x, const Vector< somenumber > &b) const
somenumber matrix_scalar_product(const Vector< somenumber > &u, const Vector< somenumber > &v) const
void add(const size_type i, const size_type j, const number value)
void TPSOR(Vector< somenumber > &v, const std::vector< size_type > &permutation, const std::vector< size_type > &inverse_permutation, const number om=1.) const
const Accessor< number, Constness > & value_type
const_iterator end() const
Accessor< number, Constness > accessor
bool operator<(const Iterator &) const
SparseMatrix< number > & copy_from(const SparseMatrix< somenumber > &source)
int operator-(const Iterator &p) const
SparseMatrixIterators::Iterator< number, true > const_iterator
virtual void reinit(const SparsityPattern &sparsity)
void Jacobi_step(Vector< somenumber > &v, const Vector< somenumber > &b, const number om=1.) const
Iterator(MatrixType *matrix, const size_type row, const size_type index) DEAL_II_DEPRECATED
const Accessor * accessor
Accessor(const SparsityPattern *matrix, const size_type row, const size_type index) DEAL_II_DEPRECATED
void vmult_add(OutVector &dst, const InVector &src) const
types::global_dof_index size_type
::ExceptionBase & ExcNumberNotFinite()
void TSOR(Vector< somenumber > &v, const number om=1.) const
SparseMatrix< number > MatrixType
void SSOR_step(Vector< somenumber > &v, const Vector< somenumber > &b, const number om=1.) const
number global_entry(const size_type i) const DEAL_II_DEPRECATED
DeclException0(ExcDifferentSparsityPatterns)
SmartPointer< const SparsityPattern, SparseMatrix< number > > cols
void block_write(std::ostream &out) const
::ExceptionBase & ExcNotInitialized()
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
static const bool zero_addition_can_be_elided
::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
void SOR_step(Vector< somenumber > &v, const Vector< somenumber > &b, const number om=1.) const
::ExceptionBase & ExcInternalError()
::ExceptionBase & ExcDivideByZero()
bool operator==(const Iterator &) const
real_type frobenius_norm() const
void compress(::VectorOperation::values)
SparseMatrix & operator*=(const number factor)
real_type linfty_norm() const
static const size_type invalid_entry
void TSOR_step(Vector< somenumber > &v, const Vector< somenumber > &b, const number om=1.) const