16 #ifndef dealii__chunk_sparse_matrix_h 17 #define dealii__chunk_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/chunk_sparsity_pattern.h> 24 #include <deal.II/lac/identity_matrix.h> 25 #include <deal.II/lac/exceptions.h> 27 DEAL_II_NAMESPACE_OPEN
29 template<
typename number>
class Vector;
43 template <
typename number,
bool Constness>
56 template <
typename number,
bool Constness>
85 template <
typename number>
99 const unsigned int row);
114 number
value()
const;
136 template <
typename,
bool>
147 template <
typename number>
180 Reference (
const Accessor *accessor,
186 operator number ()
const;
191 const Reference &operator = (
const number n)
const;
196 const Reference &operator += (
const number n)
const;
201 const Reference &operator -= (
const number n)
const;
206 const Reference &operator *= (
const number n)
const;
211 const Reference &operator /= (
const number n)
const;
232 const unsigned int row);
242 Reference
value()
const;
264 template <
typename,
bool>
280 template <
typename number,
bool Constness>
303 const unsigned int row);
344 bool operator != (
const Iterator &)
const;
359 bool operator > (
const Iterator &)
const;
367 int operator - (
const Iterator &p)
const;
372 Iterator operator + (
const unsigned int n)
const;
403 template <
typename number>
459 static const bool zero_addition_can_be_elided =
true;
569 virtual void clear ();
585 size_type m ()
const;
591 size_type n ()
const;
598 size_type n_nonzero_elements ()
const;
607 size_type n_actually_nonzero_elements ()
const;
623 std::size_t memory_consumption ()
const;
635 void set (
const size_type i,
644 void add (
const size_type i,
657 template <
typename number2>
658 void add (
const size_type
row,
659 const size_type n_cols,
660 const size_type *col_indices,
661 const number2 *values,
662 const bool elide_zero_values =
true,
663 const bool col_indices_are_sorted =
false);
705 template <
typename somenumber>
726 template <
typename ForwardIterator>
727 void copy_from (
const ForwardIterator begin,
728 const ForwardIterator end);
735 template <
typename somenumber>
749 template <
typename somenumber>
750 void add (
const number factor,
772 number operator () (
const size_type i,
773 const size_type j)
const;
787 number el (
const size_type i,
788 const size_type j)
const;
799 number diag_element (
const size_type i)
const;
805 number &diag_element (
const size_type i);
816 void extract_row_copy (
const size_type
row,
817 const size_type array_length,
818 size_type &row_length,
819 size_type *column_indices,
820 number *values)
const;
841 template <
class OutVector,
class InVector>
842 void vmult (OutVector &dst,
843 const InVector &src)
const;
860 template <
class OutVector,
class InVector>
861 void Tvmult (OutVector &dst,
862 const InVector &src)
const;
878 template <
class OutVector,
class InVector>
879 void vmult_add (OutVector &dst,
880 const InVector &src)
const;
897 template <
class OutVector,
class InVector>
898 void Tvmult_add (OutVector &dst,
899 const InVector &src)
const;
916 template <
typename somenumber>
922 template <
typename somenumber>
932 template <
typename somenumber>
950 real_type l1_norm ()
const;
959 real_type linfty_norm ()
const;
965 real_type frobenius_norm ()
const;
977 template <
typename somenumber>
980 const number omega = 1.)
const;
985 template <
typename somenumber>
988 const number om = 1.)
const;
993 template <
typename somenumber>
996 const number om = 1.)
const;
1001 template <
typename somenumber>
1004 const number om = 1.)
const;
1011 template <
typename somenumber>
1013 const number omega = 1.)
const;
1019 template <
typename somenumber>
1021 const number om = 1.)
const;
1027 template <
typename somenumber>
1029 const number om = 1.)
const;
1041 template <
typename somenumber>
1043 const std::vector<size_type> &permutation,
1044 const std::vector<size_type> &inverse_permutation,
1045 const number om = 1.)
const;
1057 template <
typename somenumber>
1059 const std::vector<size_type> &permutation,
1060 const std::vector<size_type> &inverse_permutation,
1061 const number om = 1.)
const;
1067 template <
typename somenumber>
1070 const number om = 1.)
const;
1076 template <
typename somenumber>
1079 const number om = 1.)
const;
1085 template <
typename somenumber>
1088 const number om = 1.)
const;
1183 iterator begin (
const unsigned int r);
1199 iterator end (
const unsigned int r);
1210 void print (std::ostream &out)
const;
1232 void print_formatted (std::ostream &out,
1233 const unsigned int precision = 3,
1234 const bool scientific =
true,
1235 const unsigned int width = 0,
1236 const char *zero_string =
" ",
1237 const double denominator = 1.)
const;
1244 void print_pattern(std::ostream &out,
1245 const double threshold = 0.)
const;
1257 void block_write (std::ostream &out)
const;
1275 void block_read (std::istream &in);
1287 <<
"You are trying to access the matrix entry with index <" 1288 << arg1 <<
',' << arg2
1289 <<
">, but this entry does not exist in the sparsity pattern" 1292 "The most common cause for this problem is that you used " 1293 "a method to build the sparsity pattern that did not " 1294 "(completely) take into account all of the entries you " 1295 "will later try to write into. An example would be " 1296 "building a sparsity pattern that does not include " 1297 "the entries you will write into due to constraints " 1298 "on degrees of freedom such as hanging nodes or periodic " 1299 "boundary conditions. In such cases, building the " 1300 "sparsity pattern will succeed, but you will get errors " 1301 "such as the current one at one point or other when " 1302 "trying to write into the entries of the matrix.");
1312 <<
"The iterators denote a range of " << arg1
1313 <<
" elements, but the given number of rows was " << arg2);
1346 size_type compute_location (
const size_type i,
1347 const size_type j)
const;
1366 template <
typename number>
1376 template <
typename number>
1387 template <
typename number>
1398 template <
typename number>
1404 const size_type chunk_size = cols->get_chunk_size();
1406 = cols->sparsity_pattern(i/chunk_size, j/chunk_size);
1412 return (chunk_index * chunk_size * chunk_size
1414 (i % chunk_size) * chunk_size
1421 template <
typename number>
1433 const size_type index = compute_location(i,j);
1436 ExcInvalidIndex(i,j));
1444 template <
typename number>
1457 const size_type index = compute_location(i,j);
1459 ExcInvalidIndex(i,j));
1461 val[index] +=
value;
1467 template <
typename number>
1468 template <
typename number2>
1473 const number2 *values,
1478 for (
size_type col=0; col<n_cols; ++col)
1479 add(row, col_indices[col], static_cast<number>(values[col]));
1484 template <
typename number>
1492 const size_type chunk_size = cols->get_chunk_size();
1498 number *val_ptr = val;
1499 const number *
const end_ptr = val +
1500 cols->sparsity_pattern.n_nonzero_elements()
1502 chunk_size * chunk_size;
1503 while (val_ptr != end_ptr)
1504 *val_ptr++ *= factor;
1511 template <
typename number>
1520 const number factor_inv = 1. / factor;
1522 const size_type chunk_size = cols->get_chunk_size();
1528 number *val_ptr = val;
1529 const number *
const end_ptr = val +
1530 cols->sparsity_pattern.n_nonzero_elements()
1532 chunk_size * chunk_size;
1534 while (val_ptr != end_ptr)
1535 *val_ptr++ *= factor_inv;
1542 template <
typename number>
1549 ExcInvalidIndex(i,j));
1550 return val[compute_location(i,j)];
1555 template <
typename number>
1561 const size_type index = compute_location(i,j);
1571 template <
typename number>
1581 const size_type chunk_size = cols->get_chunk_size();
1582 return val[cols->sparsity_pattern.rowstart[i/chunk_size]
1584 chunk_size * chunk_size
1586 (i % chunk_size) * chunk_size
1593 template <
typename number>
1594 template <
typename ForwardIterator>
1598 const ForwardIterator end)
1600 Assert (static_cast<size_type >(std::distance (begin, end)) == m(),
1601 ExcIteratorRange (std::distance (begin, end), m()));
1605 typedef typename std::iterator_traits<ForwardIterator>::value_type::const_iterator inner_iterator;
1607 for (ForwardIterator i=begin; i!=end; ++i, ++
row)
1609 const inner_iterator end_of_row = i->end();
1610 for (inner_iterator j=i->begin(); j!=end_of_row; ++j)
1612 set (row, j->first, j->second);
1623 template <
typename number>
1626 Accessor (
const MatrixType *matrix,
1627 const unsigned int row)
1636 template <
typename number>
1638 Accessor<number,true>::
1639 Accessor (
const MatrixType *matrix)
1647 template <
typename number>
1649 Accessor<number,true>::
1653 matrix (&a.get_matrix())
1658 template <
typename number>
1661 Accessor<number, true>::value ()
const 1663 const unsigned int chunk_size = matrix->get_sparsity_pattern().get_chunk_size();
1664 return matrix->val[reduced_index() * chunk_size * chunk_size
1666 chunk_row * chunk_size
1673 template <
typename number>
1675 const typename Accessor<number, true>::MatrixType &
1676 Accessor<number, true>::get_matrix ()
const 1683 template <
typename number>
1685 Accessor<number, false>::Reference::Reference (
1686 const Accessor *accessor,
1693 template <
typename number>
1695 Accessor<number, false>::Reference::operator number()
const 1697 const unsigned int chunk_size = accessor->matrix->get_sparsity_pattern().get_chunk_size();
1698 return accessor->matrix->val[accessor->reduced_index() * chunk_size * chunk_size
1700 accessor->chunk_row * chunk_size
1702 accessor->chunk_col];
1707 template <
typename number>
1709 const typename Accessor<number, false>::Reference &
1710 Accessor<number, false>::Reference::operator = (
const number
n)
const 1712 const unsigned int chunk_size = accessor->matrix->get_sparsity_pattern().get_chunk_size();
1713 accessor->matrix->val[accessor->reduced_index() * chunk_size * chunk_size
1715 accessor->chunk_row * chunk_size
1717 accessor->chunk_col] =
n;
1723 template <
typename number>
1725 const typename Accessor<number, false>::Reference &
1726 Accessor<number, false>::Reference::operator += (
const number n)
const 1728 const unsigned int chunk_size = accessor->matrix->get_sparsity_pattern().get_chunk_size();
1729 accessor->matrix->val[accessor->reduced_index() * chunk_size * chunk_size
1731 accessor->chunk_row * chunk_size
1733 accessor->chunk_col] +=
n;
1739 template <
typename number>
1741 const typename Accessor<number, false>::Reference &
1742 Accessor<number, false>::Reference::operator -= (
const number n)
const 1744 const unsigned int chunk_size = accessor->matrix->get_sparsity_pattern().get_chunk_size();
1745 accessor->matrix->val[accessor->reduced_index() * chunk_size * chunk_size
1747 accessor->chunk_row * chunk_size
1749 accessor->chunk_col] -=
n;
1755 template <
typename number>
1757 const typename Accessor<number, false>::Reference &
1758 Accessor<number, false>::Reference::operator *= (
const number n)
const 1760 const unsigned int chunk_size = accessor->matrix->get_sparsity_pattern().get_chunk_size();
1761 accessor->matrix->val[accessor->reduced_index() * chunk_size * chunk_size
1763 accessor->chunk_row * chunk_size
1765 accessor->chunk_col] *=
n;
1771 template <
typename number>
1773 const typename Accessor<number, false>::Reference &
1774 Accessor<number, false>::Reference::operator /= (
const number n)
const 1776 const unsigned int chunk_size = accessor->matrix->get_sparsity_pattern().get_chunk_size();
1777 accessor->matrix->val[accessor->reduced_index() * chunk_size * chunk_size
1779 accessor->chunk_row * chunk_size
1781 accessor->chunk_col] /=
n;
1787 template <
typename number>
1789 Accessor<number,false>::
1790 Accessor (MatrixType *matrix,
1791 const unsigned int row)
1800 template <
typename number>
1802 Accessor<number,false>::
1803 Accessor (MatrixType *matrix)
1811 template <
typename number>
1813 typename Accessor<number, false>::Reference
1814 Accessor<number, false>::value()
const 1816 return Reference(
this,
true);
1822 template <
typename number>
1824 typename Accessor<number, false>::MatrixType &
1825 Accessor<number, false>::get_matrix ()
const 1832 template <
typename number,
bool Constness>
1834 Iterator<number, Constness>::
1835 Iterator (MatrixType *matrix,
1836 const unsigned int row)
1838 accessor(matrix, row)
1843 template <
typename number,
bool Constness>
1845 Iterator<number, Constness>::
1846 Iterator (MatrixType *matrix)
1853 template <
typename number,
bool Constness>
1855 Iterator<number, Constness>::
1863 template <
typename number,
bool Constness>
1865 Iterator<number, Constness> &
1866 Iterator<number,Constness>::operator++ ()
1868 accessor.advance ();
1873 template <
typename number,
bool Constness>
1875 Iterator<number,Constness>
1876 Iterator<number,Constness>::operator++ (
int)
1878 const Iterator iter = *
this;
1879 accessor.advance ();
1884 template <
typename number,
bool Constness>
1886 const Accessor<number,Constness> &
1887 Iterator<number,Constness>::operator* ()
const 1893 template <
typename number,
bool Constness>
1895 const Accessor<number,Constness> *
1896 Iterator<number,Constness>::operator-> ()
const 1902 template <
typename number,
bool Constness>
1905 Iterator<number,Constness>::
1906 operator == (
const Iterator &other)
const 1908 return (accessor == other.accessor);
1912 template <
typename number,
bool Constness>
1915 Iterator<number,Constness>::
1916 operator != (
const Iterator &other)
const 1918 return ! (*
this == other);
1922 template <
typename number,
bool Constness>
1925 Iterator<number,Constness>::
1926 operator < (
const Iterator &other)
const 1928 Assert (&accessor.get_matrix() == &other.accessor.get_matrix(),
1931 return (accessor < other.accessor);
1935 template <
typename number,
bool Constness>
1938 Iterator<number,Constness>::
1939 operator > (
const Iterator &other)
const 1941 return (other < *
this);
1945 template <
typename number,
bool Constness>
1948 Iterator<number,Constness>::
1949 operator - (
const Iterator &other)
const 1951 Assert (&accessor.get_matrix() == &other.accessor.get_matrix(),
1958 Iterator copy = *
this;
1959 while (copy != other)
1967 Iterator copy = other;
1968 while (copy != *
this)
1979 template <
typename number,
bool Constness>
1981 Iterator<number,Constness>
1982 Iterator<number,Constness>::
1983 operator + (
const unsigned int n)
const 1986 for (
unsigned int i=0; i<
n; ++i)
1996 template <
typename number>
2005 template <
typename number>
2014 template <
typename number>
2023 template <
typename number>
2032 template <
typename number>
2043 template <
typename number>
2054 template <
typename number>
2065 template <
typename number>
2082 DEAL_II_NAMESPACE_CLOSE
number el(const size_type i, const size_type j) const
ChunkSparseMatrixIterators::Iterator< number, false > iterator
#define DeclException2(Exception2, type1, type2, outsequence)
static const size_type invalid_entry
ChunkSparseMatrix< number > & copy_from(const ChunkSparseMatrix< somenumber > &source)
#define AssertIndexRange(index, range)
ChunkSparseMatrixIterators::Iterator< number, true > const_iterator
bool operator==(const Accessor &) const
const_iterator end() const
types::global_dof_index size_type
static::ExceptionBase & ExcNotInitialized()
#define AssertThrow(cond, exc)
static::ExceptionBase & ExcIndexRange(int arg1, int arg2, int arg3)
static::ExceptionBase & ExcDivideByZero()
numbers::NumberTraits< number >::real_type real_type
const Accessor * accessor
const ChunkSparsityPattern & get_sparsity_pattern() const
unsigned int global_dof_index
#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 DeclException0(Exception0)
size_type compute_location(const size_type i, const size_type j) const
const Accessor< number, Constness > & value_type
Accessor(const ChunkSparsityPattern *matrix, const unsigned int row)
static::ExceptionBase & ExcNotQuadratic()
ChunkSparseMatrix< number > MatrixType
Accessor< number, Constness > accessor
const ChunkSparseMatrix< number > & get_matrix() const
Accessor< number, Constness >::MatrixType MatrixType
ChunkSparseMatrix & operator/=(const number factor)
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
void add(const size_type i, const size_type j, const number value)
SmartPointer< const ChunkSparsityPattern, ChunkSparseMatrix< number > > cols
const ChunkSparseMatrix< number > MatrixType
const_iterator begin() const
#define AssertIsFinite(number)
static::ExceptionBase & ExcInternalError()