17 #ifndef __deal2__chunk_sparse_matrix_h
18 #define __deal2__chunk_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/chunk_sparsity_pattern.h>
25 #include <deal.II/lac/identity_matrix.h>
26 #include <deal.II/lac/exceptions.h>
28 DEAL_II_NAMESPACE_OPEN
30 template<
typename number>
class Vector;
41 namespace ChunkSparseMatrixIterators
44 template <
typename number,
bool Constness>
57 template <
typename number,
bool Constness>
86 template <
typename number>
100 const unsigned int row);
115 number
value()
const;
137 template <
typename,
bool>
148 template <
typename number>
181 Reference (
const Accessor *accessor,
187 operator number ()
const;
192 const Reference &operator = (
const number n)
const;
197 const Reference &operator += (
const number n)
const;
202 const Reference &operator -= (
const number n)
const;
207 const Reference &operator *= (
const number n)
const;
212 const Reference &operator /= (
const number n)
const;
233 const unsigned int row);
243 Reference
value()
const;
265 template <
typename,
bool>
286 template <
typename number,
bool Constness>
309 const unsigned int row);
408 template <
typename number>
575 virtual void clear ();
663 template <
typename number2>
667 const number2 *values,
668 const bool elide_zero_values =
true,
669 const bool col_indices_are_sorted =
false);
708 template <
typename somenumber>
729 template <
typename ForwardIterator>
731 const ForwardIterator
end);
738 template <
typename somenumber>
752 template <
typename somenumber>
753 void add (
const number factor,
831 template <
class OutVector,
class InVector>
832 void vmult (OutVector &dst,
833 const InVector &src)
const;
849 template <
class OutVector,
class InVector>
850 void Tvmult (OutVector &dst,
851 const InVector &src)
const;
866 template <
class OutVector,
class InVector>
868 const InVector &src)
const;
884 template <
class OutVector,
class InVector>
886 const InVector &src)
const;
903 template <
typename somenumber>
909 template <
typename somenumber>
919 template <
typename somenumber>
964 template <
typename somenumber>
967 const number omega = 1.)
const;
972 template <
typename somenumber>
975 const number om = 1.)
const;
980 template <
typename somenumber>
983 const number om = 1.)
const;
988 template <
typename somenumber>
991 const number om = 1.)
const;
998 template <
typename somenumber>
1000 const number omega = 1.)
const;
1006 template <
typename somenumber>
1008 const number om = 1.)
const;
1014 template <
typename somenumber>
1016 const number om = 1.)
const;
1028 template <
typename somenumber>
1030 const std::vector<size_type> &permutation,
1031 const std::vector<size_type> &inverse_permutation,
1032 const number om = 1.)
const;
1044 template <
typename somenumber>
1046 const std::vector<size_type> &permutation,
1047 const std::vector<size_type> &inverse_permutation,
1048 const number om = 1.)
const;
1054 template <
typename somenumber>
1057 const number om = 1.)
const;
1063 template <
typename somenumber>
1066 const number om = 1.)
const;
1072 template <
typename somenumber>
1075 const number om = 1.)
const;
1197 void print (std::ostream &out)
const;
1220 const unsigned int precision = 3,
1221 const bool scientific =
true,
1222 const unsigned int width = 0,
1223 const char *zero_string =
" ",
1224 const double denominator = 1.)
const;
1232 const double threshold = 0.)
const;
1272 <<
"The entry with index <" << arg1 <<
',' << arg2
1273 <<
"> does not exist.");
1279 <<
"The index " << arg1 <<
" is not in the allowed range.");
1289 <<
"The iterators denote a range of " << arg1
1290 <<
" elements, but the given number of rows was " << arg2);
1344 template <
typename number>
1354 template <
typename number>
1365 template <
typename number>
1376 template <
typename number>
1382 const size_type chunk_size = cols->get_chunk_size();
1384 = cols->sparsity_pattern(i/chunk_size, j/chunk_size);
1390 return (chunk_index * chunk_size * chunk_size
1392 (i % chunk_size) * chunk_size
1399 template <
typename number>
1411 const size_type index = compute_location(i,j);
1414 ExcInvalidIndex(i,j));
1422 template <
typename number>
1435 const size_type index = compute_location(i,j);
1437 ExcInvalidIndex(i,j));
1439 val[index] +=
value;
1445 template <
typename number>
1446 template <
typename number2>
1451 const number2 *values,
1456 for (
size_type col=0; col<n_cols; ++col)
1457 add(row, col_indices[col], static_cast<number>(values[col]));
1462 template <
typename number>
1470 const size_type chunk_size = cols->get_chunk_size();
1476 number *val_ptr = val;
1477 const number *
const end_ptr = val +
1478 cols->sparsity_pattern.n_nonzero_elements()
1480 chunk_size * chunk_size;
1481 while (val_ptr != end_ptr)
1482 *val_ptr++ *= factor;
1489 template <
typename number>
1498 const number factor_inv = 1. / factor;
1500 const size_type chunk_size = cols->get_chunk_size();
1506 number *val_ptr = val;
1507 const number *
const end_ptr = val +
1508 cols->sparsity_pattern.n_nonzero_elements()
1510 chunk_size * chunk_size;
1512 while (val_ptr != end_ptr)
1513 *val_ptr++ *= factor_inv;
1520 template <
typename number>
1527 ExcInvalidIndex(i,j));
1528 return val[compute_location(i,j)];
1533 template <
typename number>
1539 const size_type index = compute_location(i,j);
1549 template <
typename number>
1554 Assert (m() == n(), ExcNotQuadratic());
1555 Assert (i<m(), ExcInvalidIndex1(i));
1559 const size_type chunk_size = cols->get_chunk_size();
1560 return val[cols->sparsity_pattern.rowstart[i/chunk_size]
1562 chunk_size * chunk_size
1564 (i % chunk_size) * chunk_size
1571 template <
typename number>
1572 template <
typename ForwardIterator>
1576 const ForwardIterator end)
1578 Assert (static_cast<size_type >(std::distance (begin, end)) == m(),
1579 ExcIteratorRange (std::distance (begin, end), m()));
1583 typedef typename std::iterator_traits<ForwardIterator>::value_type::const_iterator inner_iterator;
1585 for (ForwardIterator i=begin; i!=end; ++i, ++
row)
1587 const inner_iterator end_of_row = i->end();
1588 for (inner_iterator j=i->begin(); j!=end_of_row; ++j)
1590 set (row, j->first, j->second);
1599 namespace ChunkSparseMatrixIterators
1601 template <
typename number>
1604 Accessor (
const MatrixType *matrix,
1605 const unsigned int row)
1607 ChunkSparsityPatternIterators::Accessor (&matrix->get_sparsity_pattern(),
1614 template <
typename number>
1616 Accessor<number,true>::
1617 Accessor (
const MatrixType *matrix)
1619 ChunkSparsityPatternIterators::Accessor (&matrix->get_sparsity_pattern()),
1625 template <
typename number>
1627 Accessor<number,true>::
1630 ChunkSparsityPatternIterators::Accessor (a),
1631 matrix (&a.get_matrix())
1636 template <
typename number>
1639 Accessor<number, true>::value ()
const
1641 const unsigned int chunk_size = matrix->get_sparsity_pattern().get_chunk_size();
1642 return matrix->val[reduced_index() * chunk_size * chunk_size
1644 chunk_row * chunk_size
1651 template <
typename number>
1653 typename Accessor<number, true>::MatrixType &
1654 Accessor<number, true>::get_matrix ()
const
1661 template <
typename number>
1663 Accessor<number, false>::Reference::Reference (
1664 const Accessor *accessor,
1671 template <
typename number>
1673 Accessor<number, false>::Reference::operator number()
const
1675 const unsigned int chunk_size = accessor->matrix->get_sparsity_pattern().get_chunk_size();
1676 return accessor->matrix->val[accessor->reduced_index() * chunk_size * chunk_size
1678 accessor->chunk_row * chunk_size
1680 accessor->chunk_col];
1685 template <
typename number>
1687 const typename Accessor<number, false>::Reference &
1688 Accessor<number, false>::Reference::operator = (
const number n)
const
1690 const unsigned int chunk_size = accessor->matrix->get_sparsity_pattern().get_chunk_size();
1691 accessor->matrix->val[accessor->reduced_index() * chunk_size * chunk_size
1693 accessor->chunk_row * chunk_size
1695 accessor->chunk_col] = n;
1701 template <
typename number>
1703 const typename Accessor<number, false>::Reference &
1704 Accessor<number, false>::Reference::operator += (
const number n)
const
1706 const unsigned int chunk_size = accessor->matrix->get_sparsity_pattern().get_chunk_size();
1707 accessor->matrix->val[accessor->reduced_index() * chunk_size * chunk_size
1709 accessor->chunk_row * chunk_size
1711 accessor->chunk_col] += n;
1717 template <
typename number>
1719 const typename Accessor<number, false>::Reference &
1720 Accessor<number, false>::Reference::operator -= (
const number n)
const
1722 const unsigned int chunk_size = accessor->matrix->get_sparsity_pattern().get_chunk_size();
1723 accessor->matrix->val[accessor->reduced_index() * chunk_size * chunk_size
1725 accessor->chunk_row * chunk_size
1727 accessor->chunk_col] -= n;
1733 template <
typename number>
1735 const typename Accessor<number, false>::Reference &
1736 Accessor<number, false>::Reference::operator *= (
const number n)
const
1738 const unsigned int chunk_size = accessor->matrix->get_sparsity_pattern().get_chunk_size();
1739 accessor->matrix->val[accessor->reduced_index() * chunk_size * chunk_size
1741 accessor->chunk_row * chunk_size
1743 accessor->chunk_col] *= n;
1749 template <
typename number>
1751 const typename Accessor<number, false>::Reference &
1752 Accessor<number, false>::Reference::operator /= (
const number n)
const
1754 const unsigned int chunk_size = accessor->matrix->get_sparsity_pattern().get_chunk_size();
1755 accessor->matrix->val[accessor->reduced_index() * chunk_size * chunk_size
1757 accessor->chunk_row * chunk_size
1759 accessor->chunk_col] /= n;
1765 template <
typename number>
1767 Accessor<number,false>::
1768 Accessor (MatrixType *matrix,
1769 const unsigned int row)
1771 ChunkSparsityPatternIterators::Accessor (&matrix->get_sparsity_pattern(),
1778 template <
typename number>
1780 Accessor<number,false>::
1781 Accessor (MatrixType *matrix)
1783 ChunkSparsityPatternIterators::Accessor (&matrix->get_sparsity_pattern()),
1789 template <
typename number>
1791 typename Accessor<number, false>::Reference
1792 Accessor<number, false>::value()
const
1794 return Reference(
this,
true);
1800 template <
typename number>
1802 typename Accessor<number, false>::MatrixType &
1803 Accessor<number, false>::get_matrix ()
const
1810 template <
typename number,
bool Constness>
1812 Iterator<number, Constness>::
1813 Iterator (MatrixType *matrix,
1814 const unsigned int row)
1816 accessor(matrix, row)
1821 template <
typename number,
bool Constness>
1823 Iterator<number, Constness>::
1824 Iterator (MatrixType *matrix)
1831 template <
typename number,
bool Constness>
1833 Iterator<number, Constness>::
1841 template <
typename number,
bool Constness>
1843 Iterator<number, Constness> &
1844 Iterator<number,Constness>::operator++ ()
1846 accessor.advance ();
1851 template <
typename number,
bool Constness>
1853 Iterator<number,Constness>
1854 Iterator<number,Constness>::operator++ (
int)
1856 const Iterator iter = *
this;
1857 accessor.advance ();
1862 template <
typename number,
bool Constness>
1864 const Accessor<number,Constness> &
1865 Iterator<number,Constness>::operator* ()
const
1871 template <
typename number,
bool Constness>
1873 const Accessor<number,Constness> *
1874 Iterator<number,Constness>::operator-> ()
const
1880 template <
typename number,
bool Constness>
1883 Iterator<number,Constness>::
1884 operator == (
const Iterator &other)
const
1886 return (accessor == other.accessor);
1890 template <
typename number,
bool Constness>
1893 Iterator<number,Constness>::
1894 operator != (
const Iterator &other)
const
1896 return ! (*
this == other);
1900 template <
typename number,
bool Constness>
1903 Iterator<number,Constness>::
1904 operator < (
const Iterator &other)
const
1906 Assert (&accessor.get_matrix() == &other.accessor.get_matrix(),
1909 return (accessor < other.accessor);
1913 template <
typename number,
bool Constness>
1916 Iterator<number,Constness>::
1917 operator > (
const Iterator &other)
const
1919 return (other < *
this);
1923 template <
typename number,
bool Constness>
1926 Iterator<number,Constness>::
1927 operator - (
const Iterator &other)
const
1929 Assert (&accessor.get_matrix() == &other.accessor.get_matrix(),
1936 Iterator copy = *
this;
1937 while (copy != other)
1945 Iterator copy = other;
1946 while (copy != *
this)
1957 template <
typename number,
bool Constness>
1959 Iterator<number,Constness>
1960 Iterator<number,Constness>::
1961 operator + (
const unsigned int n)
const
1964 for (
unsigned int i=0; i<n; ++i)
1974 template <
typename number>
1979 return const_iterator(
this, 0);
1983 template <
typename number>
1988 return const_iterator(
this);
1992 template <
typename number>
1997 return iterator(
this, 0);
2001 template <
typename number>
2006 return iterator(
this);
2010 template <
typename number>
2016 return const_iterator(
this, r);
2021 template <
typename number>
2027 return const_iterator(
this, r+1);
2032 template <
typename number>
2038 return iterator(
this, r);
2043 template <
typename number>
2049 return iterator(
this, r+1);
2060 DEAL_II_NAMESPACE_CLOSE
number el(const size_type i, const size_type j) const
Iterator operator+(const unsigned int n) const
bool operator>(const Iterator &) const
ChunkSparseMatrixIterators::Iterator< number, false > iterator
DeclException0(ExcDifferentChunkSparsityPatterns)
void TSOR_step(Vector< somenumber > &v, const Vector< somenumber > &b, const number om=1.) const
void TSOR(Vector< somenumber > &v, const number om=1.) const
void SSOR(Vector< somenumber > &v, const number omega=1.) const
void print(std::ostream &out) const
ChunkSparseMatrixIterators::Iterator< number, true > const_iterator
void print_pattern(std::ostream &out, const double threshold=0.) const
const_iterator end() const
types::global_dof_index size_type
Iterator(MatrixType *matrix, const unsigned int row)
#define AssertThrow(cond, exc)
void SOR(Vector< somenumber > &v, const number om=1.) const
size_type n_actually_nonzero_elements() const
DeclException1(ExcInvalidIndex1, int,<< "The index "<< arg1<< " is not in the allowed range.")
bool is_finite(const double x)
void vmult_add(OutVector &dst, const InVector &src) const
size_type n_nonzero_elements() const
void SSOR_step(Vector< somenumber > &v, const Vector< somenumber > &b, const number om=1.) const
void vmult(OutVector &dst, const InVector &src) const
void SOR_step(Vector< somenumber > &v, const Vector< somenumber > &b, const number om=1.) const
numbers::NumberTraits< number >::real_type real_type
const Accessor * accessor
bool operator!=(const Iterator &) const
const ChunkSparsityPattern & get_sparsity_pattern() const
virtual ~ChunkSparseMatrix()
unsigned int global_dof_index
#define Assert(cond, exc)
bool operator==(const Iterator &) const
real_type frobenius_norm() const
number operator()(const size_type i, const size_type j) const
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
void block_write(std::ostream &out) const
somenumber matrix_scalar_product(const Vector< somenumber > &u, const Vector< somenumber > &v) const
ChunkSparseMatrix & operator*=(const number factor)
void PSOR(Vector< somenumber > &v, const std::vector< size_type > &permutation, const std::vector< size_type > &inverse_permutation, const number om=1.) const
void precondition_Jacobi(Vector< somenumber > &dst, const Vector< somenumber > &src, const number omega=1.) const
ChunkSparseMatrix< number > & operator=(const ChunkSparseMatrix< number > &)
::ExceptionBase & ExcIndexRange(int arg1, int arg2, int arg3)
ChunkSparseMatrix< number > & copy_from(const ChunkSparseMatrix< somenumber > &source)
void precondition_TSOR(Vector< somenumber > &dst, const Vector< somenumber > &src, const number om=1.) const
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)
DeclException2(ExcInvalidIndex, int, int,<< "The entry with index <"<< arg1<< ','<< arg2<< "> does not exist.")
ChunkSparseMatrix< number > MatrixType
::ExceptionBase & ExcNumberNotFinite()
Accessor< number, Constness > accessor
static const bool zero_addition_can_be_elided
somenumber residual(Vector< somenumber > &dst, const Vector< somenumber > &x, const Vector< somenumber > &b) const
void precondition_SOR(Vector< somenumber > &dst, const Vector< somenumber > &src, const number om=1.) const
const ChunkSparseMatrix< number > & get_matrix() const
Accessor< number, Constness >::MatrixType MatrixType
void precondition_SSOR(Vector< somenumber > &dst, const Vector< somenumber > &src, const number om=1.) const
void block_read(std::istream &in)
ChunkSparseMatrix & operator/=(const number factor)
::ExceptionBase & ExcNotInitialized()
bool operator<(const Iterator &) const
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
std::size_t memory_consumption() const
virtual void reinit(const ChunkSparsityPattern &sparsity)
void add(const size_type i, const size_type j, const number value)
real_type l1_norm() const
SmartPointer< const ChunkSparsityPattern, ChunkSparseMatrix< number > > cols
const ChunkSparseMatrix< number > MatrixType
somenumber matrix_norm_square(const Vector< somenumber > &v) const
const_iterator begin() const
::ExceptionBase & ExcInternalError()
::ExceptionBase & ExcDivideByZero()
int operator-(const Iterator &p) const
const Accessor< number, Constness > & operator*() const
void Tvmult(OutVector &dst, const InVector &src) const
const Accessor< number, Constness > * operator->() const
real_type linfty_norm() const
void TPSOR(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
static const size_type invalid_entry