17 #ifndef __deal2__full_matrix_h
18 #define __deal2__full_matrix_h
21 #include <deal.II/base/config.h>
22 #include <deal.II/base/numbers.h>
23 #include <deal.II/base/table.h>
24 #include <deal.II/lac/exceptions.h>
25 #include <deal.II/lac/identity_matrix.h>
26 #include <deal.II/base/tensor.h>
36 template <
typename number>
class Vector;
64 template <
typename number>
105 class const_iterator;
121 const size_type col);
128 size_type
row()
const;
140 number
value()
const;
176 const size_type col);
261 const size_type cols);
291 const size_type cols,
292 const number *entries);
327 template <
typename number2>
364 template <
typename number2>
377 template <
class MATRIX>
388 template <
class MATRIX>
407 const size_type src_r_i=0,
408 const size_type src_r_j=dim-1,
409 const size_type src_c_i=0,
410 const size_type src_c_j=dim-1,
411 const size_type dst_r=0,
412 const size_type dst_c=0);
430 const size_type src_r_i=0,
431 const size_type src_r_j=dim-1,
432 const size_type src_c_i=0,
433 const size_type src_c_j=dim-1,
434 const size_type dst_r=0,
435 const size_type dst_c=0)
const;
449 template <
typename MatrixType,
typename index_type>
451 const std::vector<index_type> &row_index_set,
452 const std::vector<index_type> &column_index_set);
466 template <
typename MatrixType,
typename index_type>
469 const std::vector<index_type> &column_index_set,
470 MatrixType &matrix)
const;
490 template <
typename number2>
492 const size_type dst_offset_i = 0,
493 const size_type dst_offset_j = 0,
494 const size_type src_offset_i = 0,
495 const size_type src_offset_j = 0);
502 template <
typename number2>
503 void fill (
const number2 *);
523 template <
typename number2>
525 const std::vector<size_type> &p_rows,
526 const std::vector<size_type> &p_cols);
538 void set (
const size_type i,
567 size_type
m ()
const;
574 size_type
n ()
const;
615 template <
typename number2>
633 template <
typename number2>
705 number
trace ()
const;
713 template <
class STREAM>
714 void print (STREAM &s,
715 const unsigned int width=5,
716 const unsigned int precision=2)
const;
760 const unsigned int precision=3,
761 const bool scientific =
true,
762 const unsigned int width = 0,
763 const char *zero_string =
" ",
764 const double denominator = 1.,
765 const double threshold = 0.)
const;
828 template <
typename number2>
829 void add (
const number a,
845 template <
typename number2>
846 void add (
const number a,
864 template <
typename number2>
865 void add (
const number a,
889 template <
typename number2>
892 const size_type dst_offset_i = 0,
893 const size_type dst_offset_j = 0,
894 const size_type src_offset_i = 0,
895 const size_type src_offset_j = 0);
904 template <
typename number2>
905 void Tadd (
const number s,
930 template <
typename number2>
933 const size_type dst_offset_i = 0,
934 const size_type dst_offset_j = 0,
935 const size_type src_offset_i = 0,
936 const size_type src_offset_j = 0);
942 void add (
const size_type row,
943 const size_type column,
962 template <
typename number2,
typename index_type>
963 void add (
const size_type row,
964 const unsigned int n_cols,
965 const index_type *col_indices,
967 const bool elide_zero_values =
true,
968 const bool col_indices_are_sorted =
false);
975 void add_row (
const size_type i,
984 void add_row (
const size_type i,
985 const number s,
const size_type j,
986 const number t,
const size_type k);
992 void add_col (
const size_type i,
1001 void add_col (
const size_type i,
1002 const number s,
const size_type j,
1003 const number t,
const size_type k);
1027 void diagadd (
const number s);
1033 template <
typename number2>
1034 void equ (
const number a,
1041 template <
typename number2>
1042 void equ (
const number a,
1051 template <
typename number2>
1052 void equ (
const number a,
1106 template <
typename number2>
1119 template <
typename number2>
1127 template <
typename number2>
1137 template <
typename number2>
1146 template <
typename number2>
1179 template <
typename number2>
1182 const bool adding=
false)
const;
1211 template <
typename number2>
1214 const bool adding=
false)
const;
1243 template <
typename number2>
1246 const bool adding=
false)
const;
1276 template <
typename number2>
1279 const bool adding=
false)
const;
1304 const bool transpose_B =
false,
1305 const bool transpose_D =
false,
1306 const number scaling = number(1.));
1325 template <
typename number2>
1328 const bool adding=
false)
const;
1337 template <
typename number2>
1360 template <
typename number2>
1363 const bool adding=
false)
const;
1373 template <
typename number2>
1387 template <
typename somenumber>
1390 const number omega = 1.)
const;
1401 template <
typename number2,
typename number3>
1423 template <
typename number2>
1437 template <
typename number2>
1456 <<
"The maximal pivot is " << arg1
1457 <<
", which is below the threshold. The matrix may be singular.");
1462 size_type, size_type, size_type,
1463 <<
"Target region not in matrix: size in this direction="
1464 << arg1 <<
", size of new matrix=" << arg2
1465 <<
", offset=" << arg3);
1487 template <
typename number>
1492 return this->n_rows();
1497 template <
typename number>
1502 return this->n_cols();
1507 template <
typename number>
1514 if (this->n_elements() != 0)
1515 std::memset (&this->values[0], 0, this->n_elements()*
sizeof(number));
1522 template <
typename number>
1523 template <
typename number2>
1532 template <
typename number>
1533 template <
class MATRIX>
1537 this->reinit (M.m(), M.n());
1542 for (size_type row = 0; row < M.m(); ++row)
1544 const typename MATRIX::const_iterator end_row = M.end(row);
1545 for (
typename MATRIX::const_iterator entry = M.begin(row);
1546 entry != end_row; ++entry)
1547 this->el(row, entry->column()) = entry->value();
1553 template <
typename number>
1554 template <
class MATRIX>
1558 this->reinit (M.n(), M.m());
1563 for (size_type row = 0; row < M.m(); ++row)
1565 const typename MATRIX::const_iterator end_row = M.end(row);
1566 for (
typename MATRIX::const_iterator entry = M.begin(row);
1567 entry != end_row; ++entry)
1568 this->el(entry->column(), row) = entry->value();
1574 template <
typename number>
1575 template <
typename MatrixType,
typename index_type>
1579 const std::vector<index_type> &row_index_set,
1580 const std::vector<index_type> &column_index_set)
1585 const size_type n_rows_submatrix = row_index_set.size();
1586 const size_type n_cols_submatrix = column_index_set.size();
1588 for (size_type sub_row = 0; sub_row < n_rows_submatrix; ++sub_row)
1589 for (size_type sub_col = 0; sub_col < n_cols_submatrix; ++sub_col)
1590 (*
this)(sub_row, sub_col) = matrix.el(row_index_set[sub_row], column_index_set[sub_col]);
1595 template <
typename number>
1596 template <
typename MatrixType,
typename index_type>
1600 const std::vector<index_type> &column_index_set,
1601 MatrixType &matrix)
const
1606 const size_type n_rows_submatrix = row_index_set.size();
1607 const size_type n_cols_submatrix = column_index_set.size();
1609 for (size_type sub_row = 0; sub_row < n_rows_submatrix; ++sub_row)
1610 for (size_type sub_col = 0; sub_col < n_cols_submatrix; ++sub_col)
1611 matrix.set(row_index_set[sub_row],
1612 column_index_set[sub_col],
1613 (*
this)(sub_row, sub_col));
1617 template <
typename number>
1624 (*this)(i,j) = value;
1629 template <
typename number>
1630 template <
typename number2>
1639 template <
typename number>
1640 template <
typename number2>
1652 template <
typename number>
1665 template <
typename number>
1674 template <
typename number>
1683 template <
typename number>
1689 return matrix->el(a_row, a_col);
1693 template <
typename number>
1700 accessor(matrix, r, c)
1704 template <
typename number>
1712 if (accessor.a_col >= accessor.matrix->n())
1721 template <
typename number>
1730 template <
typename number>
1739 template <
typename number>
1745 return (accessor.row() == other.accessor.row() &&
1746 accessor.column() == other.accessor.column());
1750 template <
typename number>
1756 return ! (*
this == other);
1760 template <
typename number>
1764 operator < (
const const_iterator &other)
const
1766 return (accessor.row() < other.accessor.row() ||
1767 (accessor.row() == other.accessor.row() &&
1768 accessor.column() < other.accessor.column()));
1772 template <
typename number>
1776 operator > (
const const_iterator &other)
const
1778 return (other < *
this);
1782 template <
typename number>
1787 return const_iterator(
this, 0, 0);
1791 template <
typename number>
1796 return const_iterator(
this,
m(), 0);
1800 template <
typename number>
1806 return const_iterator(
this, r, 0);
1811 template <
typename number>
1817 return const_iterator(
this, r+1, 0);
1822 template <
typename number>
1835 template <
typename number>
1836 template <
typename number2,
typename index_type>
1840 const unsigned int n_cols,
1841 const index_type *col_indices,
1847 for (size_type col=0; col<n_cols; ++col)
1850 this->
operator()(row,col_indices[col]) += values[col];
1855 template <
typename number>
1856 template <
class STREAM>
1860 const unsigned int w,
1861 const unsigned int p)
const
1866 const unsigned int old_precision = s.precision (p);
1867 const unsigned int old_width = s.width (w);
1869 for (size_type i=0; i<this->
m(); ++i)
1871 for (size_type j=0; j<this->
n(); ++j)
1881 s.precision(old_precision);
1888 DEAL_II_NAMESPACE_CLOSE
void vmult_add(Vector< number2 > &w, const Vector< number2 > &v) const
real_type l1_norm() const
void diagadd(const number s)
FullMatrix< number > & operator=(const FullMatrix< number > &)
#define AssertDimension(dim1, dim2)
void backward(Vector< number2 > &dst, const Vector< number2 > &src) const
const_iterator(const FullMatrix< number > *matrix, const size_type row, const size_type col)
const Accessor & operator*() const
FullMatrix & operator*=(const number factor)
bool operator<(const const_iterator &) const
bool operator==(const FullMatrix< number > &) const
void cholesky(const FullMatrix< number2 > &A)
void add_row(const size_type i, const number s, const size_type j)
real_type linfty_norm() const
number residual(Vector< number2 > &dst, const Vector< number2 > &x, const Vector< number3 > &b) const
void outer_product(const Vector< number2 > &V, const Vector< number2 > &W)
#define AssertIndexRange(index, range)
void copy_transposed(const MATRIX &)
void right_invert(const FullMatrix< number2 > &M)
void Tvmult_add(Vector< number2 > &w, const Vector< number2 > &v) const
void fill_permutation(const FullMatrix< number2 > &src, const std::vector< size_type > &p_rows, const std::vector< size_type > &p_cols)
void invert(const FullMatrix< number2 > &M)
void left_invert(const FullMatrix< number2 > &M)
void Tadd(const number s, const FullMatrix< number2 > &B)
bool operator!=(const const_iterator &) const
void vmult(Vector< number2 > &w, const Vector< number2 > &v, const bool adding=false) const
void scatter_matrix_to(const std::vector< index_type > &row_index_set, const std::vector< index_type > &column_index_set, MatrixType &matrix) const
bool is_finite(const double x)
const_iterator & operator++()
void set(const size_type i, const size_type j, const number value)
void triple_product(const FullMatrix< number > &A, const FullMatrix< number > &B, const FullMatrix< number > &D, const bool transpose_B=false, const bool transpose_D=false, const number scaling=number(1.))
void forward(Vector< number2 > &dst, const Vector< number2 > &src) const
real_type relative_symmetry_norm2() const
DeclException3(ExcInvalidDestination, size_type, size_type, size_type,<< "Target region not in matrix: size in this direction="<< arg1<< ", size of new matrix="<< arg2<< ", offset="<< arg3)
number2 matrix_scalar_product(const Vector< number2 > &u, const Vector< number2 > &v) const
const FullMatrix< number > * matrix
std::size_t memory_consumption() const
numbers::NumberTraits< number >::real_type real_type
Accessor(const FullMatrix< number > *matrix, const size_type row, const size_type col)
void swap_row(const size_type i, const size_type j)
number2 matrix_norm_square(const Vector< number2 > &v) const
std::vector< number >::reference operator()(const TableIndices< N > &indices)
#define Assert(cond, exc)
DeclException0(ExcEmptyMatrix)
real_type frobenius_norm() const
void extract_submatrix_from(const MatrixType &matrix, const std::vector< index_type > &row_index_set, const std::vector< index_type > &column_index_set)
::ExceptionBase & ExcIndexRange(int arg1, int arg2, int arg3)
void Tvmult(Vector< number2 > &w, const Vector< number2 > &v, const bool adding=false) const
void print(STREAM &s, const unsigned int width=5, const unsigned int precision=2) const
bool operator==(const const_iterator &) 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 double threshold=0.) const
const_iterator begin() const
void add_col(const size_type i, const number s, const size_type j)
::ExceptionBase & ExcIteratorPastEnd()
void mTmult(FullMatrix< number2 > &C, const FullMatrix< number2 > &B, const bool adding=false) const
FullMatrix & operator/=(const number factor)
void swap_col(const size_type i, const size_type j)
::ExceptionBase & ExcNumberNotFinite()
const Accessor * operator->() const
DeclException1(ExcNotRegular, number,<< "The maximal pivot is "<< arg1<< ", which is below the threshold. The matrix may be singular.")
void fill(const T2 *entries)
void add(const number a, const FullMatrix< number2 > &A)
void TmTmult(FullMatrix< number2 > &C, const FullMatrix< number2 > &B, const bool adding=false) const
void Tmmult(FullMatrix< number2 > &C, const FullMatrix< number2 > &B, const bool adding=false) const
FullMatrix(const size_type n=0)
std::vector< number >::reference el(const TableIndices< N > &indices)
number determinant() const
bool operator>(const const_iterator &) const
void mmult(FullMatrix< number2 > &C, const FullMatrix< number2 > &B, const bool adding=false) const
void copy_from(const MATRIX &)
::ExceptionBase & ExcScalarAssignmentOnlyForZeroValue()
const_iterator end() const
std::vector< number > values
void equ(const number a, const FullMatrix< number2 > &A)
void fill(const FullMatrix< number2 > &src, const size_type dst_offset_i=0, const size_type dst_offset_j=0, const size_type src_offset_i=0, const size_type src_offset_j=0)
void precondition_Jacobi(Vector< somenumber > &dst, const Vector< somenumber > &src, const number omega=1.) const
void copy_to(Tensor< 2, dim > &T, const size_type src_r_i=0, const size_type src_r_j=dim-1, const size_type src_c_i=0, const size_type src_c_j=dim-1, const size_type dst_r=0, const size_type dst_c=0) const