16 #ifndef dealii_full_matrix_h
17 #define dealii_full_matrix_h
38 template <
typename number>
40 template <
typename number>
77 template <
typename number>
91 "The FullMatrix class only supports basic numeric types. In particular, it "
92 "does not support automatically differentiated numbers.");
185 template <
typename number2>
214 template <
typename number2>
224 template <
typename MatrixType>
233 template <
typename MatrixType>
247 const unsigned int src_r_i = 0,
248 const unsigned int src_r_j = dim - 1,
249 const unsigned int src_c_i = 0,
250 const unsigned int src_c_j = dim - 1,
268 const unsigned int dst_r = 0,
269 const unsigned int dst_c = 0)
const;
283 template <
typename MatrixType,
typename index_type>
286 const std::vector<index_type> &row_index_set,
287 const std::vector<index_type> &column_index_set);
301 template <
typename MatrixType,
typename index_type>
304 const std::vector<index_type> &column_index_set,
305 MatrixType &
matrix)
const;
317 template <
typename number2>
329 template <
typename number2>
344 template <
typename number2>
347 const std::vector<size_type> &p_rows,
348 const std::vector<size_type> &p_cols);
415 template <
typename number2>
428 template <
typename number2>
490 template <
typename StreamType>
493 const unsigned int width = 5,
494 const unsigned int precision = 2)
const;
524 const unsigned int precision = 3,
525 const bool scientific =
true,
526 const unsigned int width = 0,
527 const char *zero_string =
" ",
528 const double denominator = 1.,
529 const double threshold = 0.,
530 const char *separator =
" ")
const;
594 template <
typename number2>
605 template <
typename number2>
620 template <
typename number2>
640 template <
typename number2>
654 template <
typename number2>
669 template <
typename number2>
693 template <
typename number2,
typename index_type>
697 const index_type *col_indices,
699 const bool elide_zero_values =
true,
700 const bool col_indices_are_sorted =
false);
758 template <
typename number2>
765 template <
typename number2>
775 template <
typename number2>
816 template <
typename number2>
828 template <
typename number2>
836 template <
typename number2>
845 template <
typename number2>
854 template <
typename number2>
882 template <
typename number2>
886 const bool adding =
false)
const;
906 template <
typename number2>
910 const bool adding =
false)
const;
930 template <
typename number2>
934 const bool adding =
false)
const;
955 template <
typename number2>
959 const bool adding =
false)
const;
975 const bool transpose_B =
false,
976 const bool transpose_D =
false,
977 const number scaling = number(1.));
991 template <
typename number2>
995 const bool adding =
false)
const;
1002 template <
typename number2>
1019 template <
typename number2>
1023 const bool adding =
false)
const;
1031 template <
typename number2>
1040 template <
typename somenumber>
1044 const number omega = 1.)
const;
1052 template <
typename number2,
typename number3>
1068 template <
typename number2>
1079 template <
typename number2>
1101 <<
"The maximal pivot is " << arg1
1102 <<
", which is below the threshold. The matrix may be singular.");
1110 <<
"Target region not in matrix: size in this direction="
1111 << arg1 <<
", size of new matrix=" << arg2
1112 <<
", offset=" << arg3);
1117 "You are attempting an operation on two vectors that "
1118 "are the same object, but the operation requires that the "
1119 "two objects are in fact different.");
1134 template <
typename number>
1138 return this->n_rows();
1143 template <
typename number>
1147 return this->n_cols();
1152 template <
typename number>
1159 if (this->n_elements() != 0)
1160 this->reset_values();
1167 template <
typename number>
1168 template <
typename number2>
1177 template <
typename number>
1178 template <
typename MatrixType>
1182 this->
reinit(M.m(), M.n());
1187 for (
size_type row = 0; row < M.m(); ++row)
1189 const typename MatrixType::const_iterator end_row = M.end(row);
1190 for (
typename MatrixType::const_iterator entry = M.begin(row);
1193 this->el(row, entry->column()) = entry->value();
1199 template <
typename number>
1203 const unsigned int src_r_i,
1204 const unsigned int src_r_j,
1205 const unsigned int src_c_i,
1206 const unsigned int src_c_j,
1210 Assert(!this->empty(), ExcEmptyMatrix());
1218 for (
size_type i = 0; i < src_r_j - src_r_i + 1; ++i)
1219 for (
size_type j = 0; j < src_c_j - src_c_i + 1; ++j)
1221 const unsigned int src_r_index =
static_cast<unsigned int>(i + src_r_i);
1222 const unsigned int src_c_index =
static_cast<unsigned int>(j + src_c_i);
1223 (*this)(i + dst_r, j + dst_c) = number(
T[src_r_index][src_c_index]);
1229 template <
typename number>
1237 const unsigned int dst_r,
1238 const unsigned int dst_c)
const
1240 Assert(!this->empty(), ExcEmptyMatrix());
1248 for (
size_type i = 0; i < src_r_j - src_r_i + 1; ++i)
1249 for (
size_type j = 0; j < src_c_j - src_c_i + 1; ++j)
1251 const unsigned int dst_r_index =
static_cast<unsigned int>(i + dst_r);
1252 const unsigned int dst_c_index =
static_cast<unsigned int>(j + dst_c);
1253 T[dst_r_index][dst_c_index] =
double((*
this)(i + src_r_i, j + src_c_i));
1259 template <
typename number>
1260 template <
typename MatrixType>
1264 this->
reinit(M.n(), M.m());
1269 for (
size_type row = 0; row < M.m(); ++row)
1271 const typename MatrixType::const_iterator end_row = M.end(row);
1272 for (
typename MatrixType::const_iterator entry = M.begin(row);
1275 this->el(entry->column(), row) = entry->value();
1281 template <
typename number>
1282 template <
typename MatrixType,
typename index_type>
1285 const MatrixType &
matrix,
1286 const std::vector<index_type> &row_index_set,
1287 const std::vector<index_type> &column_index_set)
1292 const size_type n_rows_submatrix = row_index_set.size();
1293 const size_type n_cols_submatrix = column_index_set.size();
1295 for (
size_type sub_row = 0; sub_row < n_rows_submatrix; ++sub_row)
1296 for (
size_type sub_col = 0; sub_col < n_cols_submatrix; ++sub_col)
1297 (*
this)(sub_row, sub_col) =
1298 matrix.el(row_index_set[sub_row], column_index_set[sub_col]);
1303 template <
typename number>
1304 template <
typename MatrixType,
typename index_type>
1307 const std::vector<index_type> &row_index_set,
1308 const std::vector<index_type> &column_index_set,
1309 MatrixType &
matrix)
const
1314 const size_type n_rows_submatrix = row_index_set.size();
1315 const size_type n_cols_submatrix = column_index_set.size();
1317 for (
size_type sub_row = 0; sub_row < n_rows_submatrix; ++sub_row)
1318 for (
size_type sub_col = 0; sub_col < n_cols_submatrix; ++sub_col)
1319 matrix.set(row_index_set[sub_row],
1320 column_index_set[sub_col],
1321 (*
this)(sub_row, sub_col));
1325 template <
typename number>
1331 (*this)(i, j) = value;
1336 template <
typename number>
1337 template <
typename number2>
1346 template <
typename number>
1347 template <
typename number2>
1357 template <
typename number>
1362 return iterator(
this, r, 0);
1367 template <
typename number>
1372 return iterator(
this, r + 1, 0);
1377 template <
typename number>
1382 return const_iterator(
this, r, 0);
1387 template <
typename number>
1392 return const_iterator(
this, r + 1, 0);
1397 template <
typename number>
1404 this->operator()(r, c) += v;
1409 template <
typename number>
1410 template <
typename number2,
typename index_type>
1414 const index_type *col_indices,
1420 for (
size_type col = 0; col < n_cols; ++col)
1423 this->operator()(row, col_indices[col]) +=
values[col];
1428 template <
typename number>
1429 template <
typename StreamType>
1432 const unsigned int w,
1433 const unsigned int p)
const
1435 Assert(!this->empty(), ExcEmptyMatrix());
1438 const std::streamsize old_precision = s.precision(p);
1439 const std::streamsize old_width = s.width(
w);
1441 for (
size_type i = 0; i < this->m(); ++i)
1443 for (
size_type j = 0; j < this->n(); ++j)
1447 s << this->el(i, j);
1453 s.precision(old_precision);
typename numbers::NumberTraits< number >::real_type real_type
typename Table< 2, number >::const_iterator const_iterator
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 TmTmult(FullMatrix< number2 > &C, const FullMatrix< number2 > &B, const bool adding=false) const
FullMatrix(const size_type rows, const size_type cols)
number residual(Vector< number2 > &dst, const Vector< number2 > &x, const Vector< number3 > &b) const
std::size_t memory_consumption() const
void diagadd(const number s)
void fill_permutation(const FullMatrix< number2 > &src, const std::vector< size_type > &p_rows, const std::vector< size_type > &p_cols)
void add_row(const size_type i, const number s, const size_type j)
void mmult(FullMatrix< number2 > &C, const FullMatrix< number2 > &B, const bool adding=false) const
real_type relative_symmetry_norm2() const
void add(const size_type row, const size_type column, const number value)
void copy_from(const Tensor< 2, dim > &T, const unsigned int src_r_i=0, const unsigned int src_r_j=dim - 1, const unsigned int src_c_i=0, const unsigned int src_c_j=dim - 1, const size_type dst_r=0, const size_type dst_c=0)
void equ(const number a, const FullMatrix< number2 > &A, const number b, const FullMatrix< number2 > &B)
void right_invert(const FullMatrix< number2 > &M)
void set(const size_type i, const size_type j, const number value)
void add_row(const size_type i, const number s, const size_type j, const number t, const size_type k)
FullMatrix< number > & operator=(const FullMatrix< number2 > &)
typename Table< 2, number >::iterator iterator
void add(const number a, const FullMatrix< number2 > &A, const number b, const FullMatrix< number2 > &B, const number c, const FullMatrix< number2 > &C)
FullMatrix & operator/=(const number factor)
void Tadd(const number s, const FullMatrix< number2 > &B)
void vmult(Vector< number2 > &w, const Vector< number2 > &v, const bool adding=false) const
void add_col(const size_type i, const number s, const size_type j, const number t, const size_type k)
void scatter_matrix_to(const std::vector< index_type > &row_index_set, const std::vector< index_type > &column_index_set, MatrixType &matrix) const
void swap_row(const size_type i, const size_type j)
void add(const FullMatrix< number2 > &src, const number factor, 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 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 unsigned int dst_r=0, const unsigned int dst_c=0) const
number2 matrix_norm_square(const Vector< number2 > &v) const
void equ(const number a, const FullMatrix< number2 > &A)
const_iterator end(const size_type r) const
FullMatrix(const size_type rows, const size_type cols, const number *entries)
bool operator==(const FullMatrix< number > &) const
void Tvmult(Vector< number2 > &w, const Vector< number2 > &v, const bool adding=false) 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 char *separator=" ") const
void Tmmult(FullMatrix< number2 > &C, const FullMatrix< number2 > &B, const bool adding=false) const
number2 matrix_scalar_product(const Vector< number2 > &u, const Vector< number2 > &v) const
void add(const size_type row, const size_type n_cols, const index_type *col_indices, const number2 *values, const bool elide_zero_values=true, const bool col_indices_are_sorted=false)
void add_col(const size_type i, const number s, const size_type j)
void precondition_Jacobi(Vector< somenumber > &dst, const Vector< somenumber > &src, const number omega=1.) const
FullMatrix(const IdentityMatrix &id)
void vmult_add(Vector< number2 > &w, const Vector< number2 > &v) const
FullMatrix< number > & operator=(const IdentityMatrix &id)
void invert(const FullMatrix< number2 > &M)
void cholesky(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 add(const number a, const FullMatrix< number2 > &A)
FullMatrix & operator*=(const number factor)
FullMatrix< number > & operator=(const number d)
number determinant() const
void equ(const number a, const FullMatrix< number2 > &A, const number b, const FullMatrix< number2 > &B, const number c, const FullMatrix< number2 > &C)
void fill(const number2 *)
void copy_transposed(const MatrixType &)
void print(StreamType &s, const unsigned int width=5, const unsigned int precision=2) const
void Tvmult_add(Vector< number2 > &w, const Vector< number2 > &v) const
void mTmult(FullMatrix< number2 > &C, const FullMatrix< number2 > &B, const bool adding=false) const
void left_invert(const FullMatrix< number2 > &M)
iterator begin(const size_type r)
iterator end(const size_type r)
void add(const number a, const FullMatrix< number2 > &A, const number b, const FullMatrix< number2 > &B)
void outer_product(const Vector< number2 > &V, const Vector< number2 > &W)
real_type frobenius_norm() const
void forward(Vector< number2 > &dst, const Vector< number2 > &src) const
void Tadd(const FullMatrix< number2 > &src, const number factor, 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)
FullMatrix< number > & operator=(const LAPACKFullMatrix< number2 > &)
FullMatrix(const size_type n=0)
void swap_col(const size_type i, const size_type j)
void extract_submatrix_from(const MatrixType &matrix, const std::vector< index_type > &row_index_set, const std::vector< index_type > &column_index_set)
void backward(Vector< number2 > &dst, const Vector< number2 > &src) const
void copy_from(const MatrixType &)
const_iterator begin(const size_type r) const
real_type l1_norm() const
real_type linfty_norm() const
#define DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_CLOSE
#define DeclException0(Exception0)
static ::ExceptionBase & ExcSourceEqualsDestination()
static ::ExceptionBase & ExcScalarAssignmentOnlyForZeroValue()
#define Assert(cond, exc)
static ::ExceptionBase & ExcInvalidDestination(size_type arg1, size_type arg2, size_type arg3)
static ::ExceptionBase & ExcNotRegular(number arg1)
#define AssertDimension(dim1, dim2)
#define AssertIndexRange(index, range)
#define DeclExceptionMsg(Exception, defaulttext)
#define DeclException3(Exception3, type1, type2, type3, outsequence)
static ::ExceptionBase & ExcMatrixNotPositiveDefinite()
static ::ExceptionBase & ExcEmptyMatrix()
#define DeclException1(Exception1, type1, outsequence)
@ matrix
Contents is actually a matrix.
types::global_dof_index size_type
SymmetricTensor< 2, dim, Number > C(const Tensor< 2, dim, Number > &F)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
Tensor< 2, dim, Number > w(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
void reinit(MatrixBlock< MatrixType > &v, const BlockSparsityPattern &p)