16 #ifndef dealii__sparse_matrix_h 17 #define dealii__sparse_matrix_h 20 #include <deal.II/base/config.h> 21 #include <deal.II/base/subscriptor.h> 22 #include <deal.II/base/smartpointer.h> 23 #include <deal.II/lac/sparsity_pattern.h> 24 #include <deal.II/lac/identity_matrix.h> 25 #include <deal.II/lac/exceptions.h> 26 #include <deal.II/lac/vector.h> 28 DEAL_II_NAMESPACE_OPEN
30 template <
typename number>
class Vector;
33 template <
typename number>
class SparseILU;
35 #ifdef DEAL_II_WITH_TRILINOS 59 template <
typename number,
bool Constness>
72 template <
typename number,
bool Constness>
101 template <
typename number>
115 const std::size_t index_within_matrix);
130 number
value()
const;
136 const MatrixType &get_matrix ()
const;
152 template <
typename,
bool>
163 template <
typename number>
198 Reference (
const Accessor *accessor,
204 operator number ()
const;
209 const Reference &operator = (
const number n)
const;
214 const Reference &operator += (
const number n)
const;
219 const Reference &operator -= (
const number n)
const;
224 const Reference &operator *= (
const number n)
const;
229 const Reference &operator /= (
const number n)
const;
250 const std::size_t
index);
260 Reference
value()
const;
266 MatrixType &get_matrix ()
const;
282 template <
typename,
bool>
317 template <
typename number,
bool Constness>
340 const std::size_t index_within_matrix);
376 bool operator == (
const Iterator &)
const;
381 bool operator != (
const Iterator &)
const;
390 bool operator < (
const Iterator &)
const;
396 bool operator > (
const Iterator &)
const;
404 int operator - (
const Iterator &p)
const;
409 Iterator operator + (
const size_type n)
const;
455 template <
typename number>
511 static const bool zero_addition_can_be_elided =
true;
540 #ifdef DEAL_II_WITH_CXX11 595 #ifdef DEAL_II_WITH_CXX11 648 virtual void clear ();
664 size_type m ()
const;
670 size_type n ()
const;
675 size_type get_row_length (
const size_type
row)
const;
682 std::size_t n_nonzero_elements ()
const;
693 std::size_t n_actually_nonzero_elements (
const double threshold = 0.)
const;
709 std::size_t memory_consumption ()
const;
726 void set (
const size_type i,
745 template <
typename number2>
746 void set (
const std::vector<size_type> &indices,
748 const bool elide_zero_values =
false);
755 template <
typename number2>
756 void set (
const std::vector<size_type> &row_indices,
757 const std::vector<size_type> &col_indices,
759 const bool elide_zero_values =
false);
771 template <
typename number2>
772 void set (
const size_type
row,
773 const std::vector<size_type> &col_indices,
774 const std::vector<number2> &values,
775 const bool elide_zero_values =
false);
786 template <
typename number2>
787 void set (
const size_type
row,
788 const size_type n_cols,
789 const size_type *col_indices,
790 const number2 *values,
791 const bool elide_zero_values =
false);
798 void add (
const size_type i,
816 template <
typename number2>
817 void add (
const std::vector<size_type> &indices,
819 const bool elide_zero_values =
true);
826 template <
typename number2>
827 void add (
const std::vector<size_type> &row_indices,
828 const std::vector<size_type> &col_indices,
830 const bool elide_zero_values =
true);
841 template <
typename number2>
842 void add (
const size_type row,
843 const std::vector<size_type> &col_indices,
844 const std::vector<number2> &values,
845 const bool elide_zero_values =
true);
856 template <
typename number2>
857 void add (
const size_type row,
858 const size_type n_cols,
859 const size_type *col_indices,
860 const number2 *values,
861 const bool elide_zero_values =
true,
862 const bool col_indices_are_sorted =
false);
904 template <
typename somenumber>
924 template <
typename ForwardIterator>
925 void copy_from (
const ForwardIterator begin,
926 const ForwardIterator end);
933 template <
typename somenumber>
936 #ifdef DEAL_II_WITH_TRILINOS 961 template <
typename somenumber>
962 void add (
const number factor,
984 number operator () (
const size_type i,
985 const size_type j)
const;
999 number el (
const size_type i,
1000 const size_type j)
const;
1011 number diag_element (
const size_type i)
const;
1017 number &diag_element (
const size_type i);
1040 template <
class OutVector,
class InVector>
1041 void vmult (OutVector &dst,
1042 const InVector &src)
const;
1059 template <
class OutVector,
class InVector>
1060 void Tvmult (OutVector &dst,
1061 const InVector &src)
const;
1079 template <
class OutVector,
class InVector>
1080 void vmult_add (OutVector &dst,
1081 const InVector &src)
const;
1098 template <
class OutVector,
class InVector>
1099 void Tvmult_add (OutVector &dst,
1100 const InVector &src)
const;
1119 template <
typename somenumber>
1127 template <
typename somenumber>
1140 template <
typename somenumber>
1180 template <
typename numberB,
typename numberC>
1184 const bool rebuild_sparsity_pattern =
true)
const;
1210 template <
typename numberB,
typename numberC>
1214 const bool rebuild_sparsity_pattern =
true)
const;
1229 real_type l1_norm ()
const;
1238 real_type linfty_norm ()
const;
1244 real_type frobenius_norm ()
const;
1256 template <
typename somenumber>
1259 const number omega = 1.)
const;
1267 template <
typename somenumber>
1270 const number omega = 1.,
1271 const std::vector<std::size_t> &pos_right_of_diagonal=std::vector<std::size_t>())
const;
1276 template <
typename somenumber>
1279 const number om = 1.)
const;
1284 template <
typename somenumber>
1287 const number om = 1.)
const;
1294 template <
typename somenumber>
1296 const number omega = 1.)
const;
1302 template <
typename somenumber>
1304 const number om = 1.)
const;
1310 template <
typename somenumber>
1312 const number om = 1.)
const;
1324 template <
typename somenumber>
1326 const std::vector<size_type> &permutation,
1327 const std::vector<size_type> &inverse_permutation,
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;
1351 template <
typename somenumber>
1354 const number om = 1.)
const;
1360 template <
typename somenumber>
1363 const number om = 1.)
const;
1369 template <
typename somenumber>
1372 const number om = 1.)
const;
1378 template <
typename somenumber>
1381 const number om = 1.)
const;
1425 iterator begin (
const size_type r);
1459 template <
class StreamType>
1460 void print (StreamType &out,
1461 const bool across =
false,
1462 const bool diagonal_first =
true)
const;
1484 void print_formatted (std::ostream &out,
1485 const unsigned int precision = 3,
1486 const bool scientific =
true,
1487 const unsigned int width = 0,
1488 const char *zero_string =
" ",
1489 const double denominator = 1.)
const;
1496 void print_pattern(std::ostream &out,
1497 const double threshold = 0.)
const;
1509 void block_write (std::ostream &out)
const;
1527 void block_read (std::istream &in);
1539 <<
"You are trying to access the matrix entry with index <" 1540 << arg1 <<
',' << arg2
1541 <<
">, but this entry does not exist in the sparsity pattern " 1544 "The most common cause for this problem is that you used " 1545 "a method to build the sparsity pattern that did not " 1546 "(completely) take into account all of the entries you " 1547 "will later try to write into. An example would be " 1548 "building a sparsity pattern that does not include " 1549 "the entries you will write into due to constraints " 1550 "on degrees of freedom such as hanging nodes or periodic " 1551 "boundary conditions. In such cases, building the " 1552 "sparsity pattern will succeed, but you will get errors " 1553 "such as the current one at one point or other when " 1554 "trying to write into the entries of the matrix.");
1559 "When copying one sparse matrix into another, " 1560 "or when adding one sparse matrix to another, " 1561 "both matrices need to refer to the same " 1562 "sparsity pattern.");
1568 <<
"The iterators denote a range of " << arg1
1569 <<
" elements, but the given number of rows was " << arg2);
1622 template <
typename>
friend class SparseILU;
1641 template <
typename number>
1650 template <
typename number>
1660 template <
typename number>
1669 const size_type
index = cols->operator()(i, j);
1676 (value == number()),
1677 ExcInvalidIndex(i, j));
1686 template <
typename number>
1687 template <
typename number2>
1692 const bool elide_zero_values)
1694 Assert (indices.size() == values.
m(),
1698 for (size_type i=0; i<indices.size(); ++i)
1699 set (indices[i], indices.size(), &indices[0], &values(i,0),
1705 template <
typename number>
1706 template <
typename number2>
1710 const std::vector<size_type> &col_indices,
1712 const bool elide_zero_values)
1714 Assert (row_indices.size() == values.
m(),
1716 Assert (col_indices.size() == values.
n(),
1719 for (size_type i=0; i<row_indices.size(); ++i)
1720 set (row_indices[i], col_indices.size(), &col_indices[0], &values(i,0),
1726 template <
typename number>
1727 template <
typename number2>
1731 const std::vector<size_type> &col_indices,
1732 const std::vector<number2> &values,
1733 const bool elide_zero_values)
1735 Assert (col_indices.size() == values.size(),
1738 set (
row, col_indices.size(), &col_indices[0], &values[0],
1744 template <
typename number>
1753 if (value == number())
1756 const size_type index = cols->operator()(i, j);
1763 (value == number()),
1764 ExcInvalidIndex(i, j));
1773 template <
typename number>
1774 template <
typename number2>
1779 const bool elide_zero_values)
1781 Assert (indices.size() == values.
m(),
1785 for (size_type i=0; i<indices.size(); ++i)
1786 add (indices[i], indices.size(), &indices[0], &values(i,0),
1792 template <
typename number>
1793 template <
typename number2>
1797 const std::vector<size_type> &col_indices,
1799 const bool elide_zero_values)
1801 Assert (row_indices.size() == values.
m(),
1803 Assert (col_indices.size() == values.
n(),
1806 for (size_type i=0; i<row_indices.size(); ++i)
1807 add (row_indices[i], col_indices.size(), &col_indices[0], &values(i,0),
1813 template <
typename number>
1814 template <
typename number2>
1818 const std::vector<size_type> &col_indices,
1819 const std::vector<number2> &values,
1820 const bool elide_zero_values)
1822 Assert (col_indices.size() == values.size(),
1825 add (row, col_indices.size(), &col_indices[0], &values[0],
1831 template <
typename number>
1839 number *val_ptr = &val[0];
1840 const number *
const end_ptr = &val[cols->n_nonzero_elements()];
1842 while (val_ptr != end_ptr)
1843 *val_ptr++ *= factor;
1850 template <
typename number>
1859 const number factor_inv = number(1.) / factor;
1861 number *val_ptr = &val[0];
1862 const number *
const end_ptr = &val[cols->n_nonzero_elements()];
1864 while (val_ptr != end_ptr)
1865 *val_ptr++ *= factor_inv;
1872 template <
typename number>
1875 const size_type j)
const 1879 ExcInvalidIndex(i,j));
1880 return val[cols->operator()(i,j)];
1885 template <
typename number>
1888 const size_type j)
const 1891 const size_type index = cols->operator()(i,j);
1901 template <
typename number>
1911 return val[cols->rowstart[i]];
1916 template <
typename number>
1926 return val[cols->rowstart[i]];
1931 template <
typename number>
1932 template <
typename ForwardIterator>
1935 const ForwardIterator end)
1937 Assert (static_cast<size_type>(std::distance (begin, end)) == m(),
1938 ExcIteratorRange (std::distance (begin, end), m()));
1942 typedef typename std::iterator_traits<ForwardIterator>::value_type::const_iterator inner_iterator;
1944 for (ForwardIterator i=begin; i!=end; ++i, ++
row)
1946 const inner_iterator end_of_row = i->end();
1947 for (inner_iterator j=i->begin(); j!=end_of_row; ++j)
1949 set (row, j->first, j->second);
1959 template <
typename number>
1961 Accessor<number,true>::
1962 Accessor (
const MatrixType *
matrix,
1963 const std::size_t index_within_matrix)
1966 index_within_matrix),
1972 template <
typename number>
1974 Accessor<number,true>::
1975 Accessor (
const MatrixType *matrix)
1983 template <
typename number>
1985 Accessor<number,true>::
1989 matrix (&a.get_matrix())
1994 template <
typename number>
1997 Accessor<number, true>::value ()
const 2000 return matrix->val[index_within_sparsity];
2005 template <
typename number>
2007 const typename Accessor<number, true>::MatrixType &
2008 Accessor<number, true>::get_matrix ()
const 2015 template <
typename number>
2017 Accessor<number, false>::Reference::Reference (
2018 const Accessor *accessor,
2025 template <
typename number>
2027 Accessor<number, false>::Reference::operator number()
const 2029 AssertIndexRange(accessor->index_within_sparsity, accessor->matrix->n_nonzero_elements());
2030 return accessor->matrix->val[accessor->index_within_sparsity];
2035 template <
typename number>
2037 const typename Accessor<number, false>::Reference &
2038 Accessor<number, false>::Reference::operator = (
const number
n)
const 2040 AssertIndexRange(accessor->index_within_sparsity, accessor->matrix->n_nonzero_elements());
2041 accessor->matrix->val[accessor->index_within_sparsity] =
n;
2047 template <
typename number>
2049 const typename Accessor<number, false>::Reference &
2050 Accessor<number, false>::Reference::operator += (
const number n)
const 2052 AssertIndexRange(accessor->index_within_sparsity, accessor->matrix->n_nonzero_elements());
2053 accessor->matrix->val[accessor->index_within_sparsity] +=
n;
2059 template <
typename number>
2061 const typename Accessor<number, false>::Reference &
2062 Accessor<number, false>::Reference::operator -= (
const number n)
const 2064 AssertIndexRange(accessor->index_within_sparsity, accessor->matrix->n_nonzero_elements());
2065 accessor->matrix->val[accessor->index_within_sparsity] -=
n;
2071 template <
typename number>
2073 const typename Accessor<number, false>::Reference &
2074 Accessor<number, false>::Reference::operator *= (
const number n)
const 2076 AssertIndexRange(accessor->index_within_sparsity, accessor->matrix->n_nonzero_elements());
2077 accessor->matrix->val[accessor->index_within_sparsity] *=
n;
2083 template <
typename number>
2085 const typename Accessor<number, false>::Reference &
2086 Accessor<number, false>::Reference::operator /= (
const number n)
const 2088 AssertIndexRange(accessor->index_within_sparsity, accessor->matrix->n_nonzero_elements());
2089 accessor->matrix->val[accessor->index_within_sparsity] /=
n;
2095 template <
typename number>
2097 Accessor<number,false>::
2098 Accessor (MatrixType *matrix,
2099 const std::size_t index)
2108 template <
typename number>
2110 Accessor<number,false>::
2111 Accessor (MatrixType *matrix)
2119 template <
typename number>
2121 typename Accessor<number, false>::Reference
2122 Accessor<number, false>::value()
const 2124 return Reference(
this,
true);
2130 template <
typename number>
2132 typename Accessor<number, false>::MatrixType &
2133 Accessor<number, false>::get_matrix ()
const 2140 template <
typename number,
bool Constness>
2142 Iterator<number, Constness>::
2143 Iterator (MatrixType *matrix,
2144 const std::size_t index)
2146 accessor(matrix, index)
2151 template <
typename number,
bool Constness>
2153 Iterator<number, Constness>::
2154 Iterator (MatrixType *matrix)
2161 template <
typename number,
bool Constness>
2163 Iterator<number, Constness>::
2171 template <
typename number,
bool Constness>
2173 Iterator<number, Constness> &
2174 Iterator<number,Constness>::operator++ ()
2176 accessor.advance ();
2181 template <
typename number,
bool Constness>
2183 Iterator<number,Constness>
2184 Iterator<number,Constness>::operator++ (
int)
2186 const Iterator iter = *
this;
2187 accessor.advance ();
2192 template <
typename number,
bool Constness>
2194 const Accessor<number,Constness> &
2195 Iterator<number,Constness>::operator* ()
const 2201 template <
typename number,
bool Constness>
2203 const Accessor<number,Constness> *
2204 Iterator<number,Constness>::operator-> ()
const 2210 template <
typename number,
bool Constness>
2213 Iterator<number,Constness>::
2214 operator == (
const Iterator &other)
const 2216 return (accessor == other.accessor);
2220 template <
typename number,
bool Constness>
2223 Iterator<number,Constness>::
2224 operator != (
const Iterator &other)
const 2226 return ! (*
this == other);
2230 template <
typename number,
bool Constness>
2233 Iterator<number,Constness>::
2234 operator < (
const Iterator &other)
const 2236 Assert (&accessor.get_matrix() == &other.accessor.get_matrix(),
2239 return (accessor < other.accessor);
2243 template <
typename number,
bool Constness>
2246 Iterator<number,Constness>::
2247 operator > (
const Iterator &other)
const 2249 return (other < *
this);
2253 template <
typename number,
bool Constness>
2256 Iterator<number,Constness>::
2257 operator - (
const Iterator &other)
const 2259 Assert (&accessor.get_matrix() == &other.accessor.get_matrix(),
2262 return (*this)->index_within_sparsity - other->index_within_sparsity;
2267 template <
typename number,
bool Constness>
2269 Iterator<number,Constness>
2270 Iterator<number,Constness>::
2271 operator + (
const size_type n)
const 2274 for (size_type i=0; i<
n; ++i)
2284 template <
typename number>
2293 template <
typename number>
2302 template <
typename number>
2311 template <
typename number>
2320 template <
typename number>
2332 template <
typename number>
2344 template <
typename number>
2356 template <
typename number>
2368 template <
typename number>
2369 template <
class StreamType>
2373 const bool diagonal_first)
const 2378 bool hanging_diagonal =
false;
2379 number diagonal = number();
2381 for (size_type i=0; i<
cols->
rows; ++i)
2383 for (size_type j=
cols->
rowstart[i]; j<cols->rowstart[i+1]; ++j)
2388 hanging_diagonal =
true;
2395 out <<
' ' << i <<
',' << i <<
':' << diagonal;
2397 out <<
'(' << i <<
',' << i <<
") " << diagonal << std::endl;
2398 hanging_diagonal =
false;
2403 out <<
"(" << i <<
"," <<
cols->
colnums[j] <<
") " <<
val[j] << std::endl;
2406 if (hanging_diagonal)
2409 out <<
' ' << i <<
',' << i <<
':' << diagonal;
2411 out <<
'(' << i <<
',' << i <<
") " << diagonal << std::endl;
2412 hanging_diagonal =
false;
2420 template <
typename number>
2431 template <
typename number>
2445 DEAL_II_NAMESPACE_CLOSE
#define DeclException2(Exception2, type1, type2, outsequence)
numbers::NumberTraits< number >::real_type real_type
Accessor< number, Constness >::MatrixType MatrixType
static const size_type invalid_entry
void print(StreamType &out, const bool across=false, const bool diagonal_first=true) const
#define AssertIndexRange(index, range)
TrilinosScalar diag_element(const size_type i) const
SparseMatrixIterators::Iterator< number, false > iterator
static::ExceptionBase & ExcNotInitialized()
static::ExceptionBase & ExcIndexRange(int arg1, int arg2, int arg3)
static::ExceptionBase & ExcDivideByZero()
const_iterator begin() const
SparseMatrix & operator*=(const TrilinosScalar factor)
TrilinosScalar value() const
TrilinosScalar el(const size_type i, const size_type j) const
unsigned int global_dof_index
void add(const size_type i, const size_type j, const TrilinosScalar value)
#define Assert(cond, exc)
static::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
types::global_dof_index size_type
#define DeclExceptionMsg(Exception, defaulttext)
#define DeclException0(Exception0)
const SparseMatrix< number > MatrixType
TrilinosScalar operator()(const size_type i, const size_type j) const
const Accessor< number, Constness > & value_type
const_iterator end() const
Accessor< number, Constness > accessor
SparseMatrixIterators::Iterator< number, true > const_iterator
static::ExceptionBase & ExcNotQuadratic()
const Accessor * accessor
types::global_dof_index size_type
const SparsityPattern & get_sparsity_pattern() const
SparseMatrix< number > MatrixType
void copy_from(const SparseMatrix &source)
SparseMatrix & operator/=(const TrilinosScalar factor)
SmartPointer< const SparsityPattern, SparseMatrix< number > > cols
void set(const size_type i, const size_type j, const TrilinosScalar value)
#define AssertIsFinite(number)
static::ExceptionBase & ExcInternalError()