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;
35 template <
typename number>
class SparseILU;
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;
148 MatrixType &get_matrix ()
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;
285 MatrixType &get_matrix ()
const;
301 template <
typename,
bool>
342 template <
typename number,
bool Constness>
376 const std::size_t index_within_matrix);
412 bool operator == (
const Iterator &)
const;
417 bool operator != (
const Iterator &)
const;
426 bool operator < (
const Iterator &)
const;
432 bool operator > (
const Iterator &)
const;
440 int operator - (
const Iterator &p)
const;
445 Iterator operator + (
const size_type n)
const;
492 template <
typename number>
548 static const bool zero_addition_can_be_elided =
true;
658 virtual void clear ();
674 size_type m ()
const;
680 size_type n ()
const;
685 size_type get_row_length (
const size_type
row)
const;
692 size_type n_nonzero_elements ()
const;
703 size_type n_actually_nonzero_elements (
const double threshold = 0.)
const;
719 std::size_t memory_consumption ()
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>
932 void copy_from (
const ForwardIterator begin,
933 const ForwardIterator end);
940 template <
typename somenumber>
943 #ifdef DEAL_II_WITH_TRILINOS 968 template <
typename somenumber>
969 void add (
const number factor,
991 number operator () (
const size_type i,
992 const size_type j)
const;
1006 number el (
const size_type i,
1007 const size_type j)
const;
1019 number diag_element (
const size_type i)
const;
1025 number &diag_element (
const size_type i);
1035 number raw_entry (
const size_type row,
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>
1114 void vmult_add (OutVector &dst,
1115 const InVector &src)
const;
1131 template <
class OutVector,
class InVector>
1132 void Tvmult_add (OutVector &dst,
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;
1244 real_type l1_norm ()
const;
1254 real_type linfty_norm ()
const;
1260 real_type frobenius_norm ()
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;
1469 iterator begin (
const size_type r);
1498 template <
class STREAM>
1499 void print (STREAM &out,
1500 const bool across =
false,
1501 const bool diagonal_first =
true)
const;
1523 void print_formatted (std::ostream &out,
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;
1535 void print_pattern(std::ostream &out,
1536 const double threshold = 0.)
const;
1548 void block_write (std::ostream &out)
const;
1566 void block_read (std::istream &in);
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);
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>
2029 Accessor<number,true>::
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>
2401 template <
typename number>
2410 template <
typename number>
2419 template <
typename number>
2428 template <
typename number>
2440 template <
typename number>
2452 template <
typename number>
2464 template <
typename number>
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)
2496 hanging_diagonal =
true;
2503 out <<
' ' << i <<
',' << i <<
':' << diagonal;
2505 out <<
'(' << i <<
',' << i <<
") " << diagonal << std::endl;
2506 hanging_diagonal =
false;
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
number raw_entry(const size_type row, const size_type index) const DEAL_II_DEPRECATED
#define DeclException2(Exception2, type1, type2, outsequence)
numbers::NumberTraits< number >::real_type real_type
Accessor< number, Constness >::MatrixType MatrixType
#define AssertIndexRange(index, range)
SparseMatrixIterators::Iterator< number, false > iterator
DeclException0(ExcBeyondEndOfMatrix)
bool is_finite(const double x)
const_iterator begin() const
void print(STREAM &out, const bool across=false, const bool diagonal_first=true) const
TrilinosScalar value() const
TrilinosScalar operator()(const size_type i, const size_type j) const
#define DeclException1(Exception1, type1, outsequence)
unsigned int global_dof_index
void add(const size_type i, const size_type j, const TrilinosScalar value)
#define Assert(cond, exc)
SparseMatrix & operator/=(const TrilinosScalar factor)
types::global_dof_index size_type
const SparsityPattern & get_sparsity_pattern() const
const SparseMatrix< number > MatrixType
::ExceptionBase & ExcIndexRange(int arg1, int arg2, int arg3)
BlockCompressedSparsityPattern CompressedBlockSparsityPattern DEAL_II_DEPRECATED
const Accessor< number, Constness > & value_type
const_iterator end() const
Accessor< number, Constness > accessor
void copy_from(const SparseMatrix &source)
SparseMatrixIterators::Iterator< number, true > const_iterator
SparseMatrix & operator*=(const TrilinosScalar factor)
const Accessor * accessor
TrilinosScalar el(const size_type i, const size_type j) const
types::global_dof_index size_type
::ExceptionBase & ExcNumberNotFinite()
SparseMatrix< number > MatrixType
number global_entry(const size_type i) const DEAL_II_DEPRECATED
TrilinosScalar diag_element(const size_type i) const
SmartPointer< const SparsityPattern, SparseMatrix< number > > cols
::ExceptionBase & ExcNotInitialized()
::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
::ExceptionBase & ExcInternalError()
::ExceptionBase & ExcDivideByZero()
void set(const size_type i, const size_type j, const TrilinosScalar value)
static const size_type invalid_entry