Reference documentation for deal.II version 9.4.1
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
Loading...
Searching...
No Matches
full_matrix.h
Go to the documentation of this file.
1// ---------------------------------------------------------------------
2//
3// Copyright (C) 1999 - 2021 by the deal.II authors
4//
5// This file is part of the deal.II library.
6//
7// The deal.II library is free software; you can use it, redistribute
8// it, and/or modify it under the terms of the GNU Lesser General
9// Public License as published by the Free Software Foundation; either
10// version 2.1 of the License, or (at your option) any later version.
11// The full text of the license can be found in the file LICENSE.md at
12// the top level directory of deal.II.
13//
14// ---------------------------------------------------------------------
15
16#ifndef dealii_full_matrix_h
17#define dealii_full_matrix_h
18
19
20#include <deal.II/base/config.h>
21
23#include <deal.II/base/table.h>
24#include <deal.II/base/tensor.h>
25
27
30
31#include <cstring>
32#include <iomanip>
33#include <vector>
34
36
37
38// forward declarations
39#ifndef DOXYGEN
40template <typename number>
41class Vector;
42template <typename number>
44#endif
45
68template <typename number>
69class FullMatrix : public Table<2, number>
70{
71public:
72 // The assertion in full_matrix.templates.h for whether or not a number is
73 // finite is not compatible for AD number types.
74 static_assert(
76 "The FullMatrix class does not support auto-differentiable numbers.");
77
81 using size_type = std::size_t;
82
87 using value_type = number;
88
93
98
102 using Table<2, number>::begin;
103
107 using Table<2, number>::end;
108
119
124
134 explicit FullMatrix(const size_type n = 0);
135
139 FullMatrix(const size_type rows, const size_type cols);
140
145 FullMatrix(const size_type rows, const size_type cols, const number *entries);
146
170 template <typename number2>
173
182 operator=(const number d);
183
194
199 template <typename number2>
202
203
209 template <typename MatrixType>
210 void
211 copy_from(const MatrixType &);
212
218 template <typename MatrixType>
219 void
220 copy_transposed(const MatrixType &);
221
229 template <int dim>
230 void
232 const unsigned int src_r_i = 0,
233 const unsigned int src_r_j = dim - 1,
234 const unsigned int src_c_i = 0,
235 const unsigned int src_c_j = dim - 1,
236 const size_type dst_r = 0,
237 const size_type dst_c = 0);
238
246 template <int dim>
247 void
249 const size_type src_r_i = 0,
250 const size_type src_r_j = dim - 1,
251 const size_type src_c_i = 0,
252 const size_type src_c_j = dim - 1,
253 const unsigned int dst_r = 0,
254 const unsigned int dst_c = 0) const;
255
268 template <typename MatrixType, typename index_type>
269 void
270 extract_submatrix_from(const MatrixType & matrix,
271 const std::vector<index_type> &row_index_set,
272 const std::vector<index_type> &column_index_set);
273
286 template <typename MatrixType, typename index_type>
287 void
288 scatter_matrix_to(const std::vector<index_type> &row_index_set,
289 const std::vector<index_type> &column_index_set,
290 MatrixType & matrix) const;
291
302 template <typename number2>
303 void
305 const size_type dst_offset_i = 0,
306 const size_type dst_offset_j = 0,
307 const size_type src_offset_i = 0,
308 const size_type src_offset_j = 0);
309
310
314 template <typename number2>
315 void
316 fill(const number2 *);
317
329 template <typename number2>
330 void
332 const std::vector<size_type> &p_rows,
333 const std::vector<size_type> &p_cols);
334
345 void
346 set(const size_type i, const size_type j, const number value);
362 bool
364
370 m() const;
371
377 n() const;
378
384 bool
385 all_zero() const;
386
402 template <typename number2>
403 number2
405
415 template <typename number2>
416 number2
418 const Vector<number2> &v) const;
419
425 l1_norm() const;
426
432 linfty_norm() const;
433
443
454
460 number
461 determinant() const;
462
468 number
469 trace() const;
470
477 template <class StreamType>
478 void
479 print(StreamType & s,
480 const unsigned int width = 5,
481 const unsigned int precision = 2) const;
482
505 void
506 print_formatted(std::ostream & out,
507 const unsigned int precision = 3,
508 const bool scientific = true,
509 const unsigned int width = 0,
510 const char * zero_string = " ",
511 const double denominator = 1.,
512 const double threshold = 0.) const;
513
518 std::size_t
520
522
524
529 begin(const size_type r);
530
535 end(const size_type r);
536
541 begin(const size_type r) const;
542
547 end(const size_type r) const;
548
550
552
556 FullMatrix &
557 operator*=(const number factor);
558
562 FullMatrix &
563 operator/=(const number factor);
564
572 template <typename number2>
573 void
574 add(const number a, const FullMatrix<number2> &A);
575
583 template <typename number2>
584 void
585 add(const number a,
586 const FullMatrix<number2> &A,
587 const number b,
588 const FullMatrix<number2> &B);
589
598 template <typename number2>
599 void
600 add(const number a,
601 const FullMatrix<number2> &A,
602 const number b,
603 const FullMatrix<number2> &B,
604 const number c,
605 const FullMatrix<number2> &C);
606
618 template <typename number2>
619 void
621 const number factor,
622 const size_type dst_offset_i = 0,
623 const size_type dst_offset_j = 0,
624 const size_type src_offset_i = 0,
625 const size_type src_offset_j = 0);
626
632 template <typename number2>
633 void
634 Tadd(const number s, const FullMatrix<number2> &B);
635
647 template <typename number2>
648 void
650 const number factor,
651 const size_type dst_offset_i = 0,
652 const size_type dst_offset_j = 0,
653 const size_type src_offset_i = 0,
654 const size_type src_offset_j = 0);
655
659 void
660 add(const size_type row, const size_type column, const number value);
661
671 template <typename number2, typename index_type>
672 void
673 add(const size_type row,
674 const size_type n_cols,
675 const index_type *col_indices,
676 const number2 * values,
677 const bool elide_zero_values = true,
678 const bool col_indices_are_sorted = false);
679
683 void
684 add_row(const size_type i, const number s, const size_type j);
685
690 void
692 const number s,
693 const size_type j,
694 const number t,
695 const size_type k);
696
700 void
701 add_col(const size_type i, const number s, const size_type j);
702
707 void
709 const number s,
710 const size_type j,
711 const number t,
712 const size_type k);
713
717 void
718 swap_row(const size_type i, const size_type j);
719
723 void
724 swap_col(const size_type i, const size_type j);
725
730 void
731 diagadd(const number s);
732
736 template <typename number2>
737 void
738 equ(const number a, const FullMatrix<number2> &A);
739
743 template <typename number2>
744 void
745 equ(const number a,
746 const FullMatrix<number2> &A,
747 const number b,
748 const FullMatrix<number2> &B);
749
753 template <typename number2>
754 void
755 equ(const number a,
756 const FullMatrix<number2> &A,
757 const number b,
758 const FullMatrix<number2> &B,
759 const number c,
760 const FullMatrix<number2> &C);
761
768 void
770
785 void
787
794 template <typename number2>
795 void
797
806 template <typename number2>
807 void
809
814 template <typename number2>
815 void
817
823 template <typename number2>
824 void
826
832 template <typename number2>
833 void
835
837
839
858 template <typename number2>
859 void
861 const FullMatrix<number2> &B,
862 const bool adding = false) const;
863
882 template <typename number2>
883 void
885 const FullMatrix<number2> &B,
886 const bool adding = false) const;
887
906 template <typename number2>
907 void
909 const FullMatrix<number2> &B,
910 const bool adding = false) const;
911
931 template <typename number2>
932 void
934 const FullMatrix<number2> &B,
935 const bool adding = false) const;
936
947 void
949 const FullMatrix<number> &B,
950 const FullMatrix<number> &D,
951 const bool transpose_B = false,
952 const bool transpose_D = false,
953 const number scaling = number(1.));
954
967 template <typename number2>
968 void
970 const Vector<number2> &v,
971 const bool adding = false) const;
972
978 template <typename number2>
979 void
981
995 template <typename number2>
996 void
998 const Vector<number2> &v,
999 const bool adding = false) const;
1000
1007 template <typename number2>
1008 void
1010
1016 template <typename somenumber>
1017 void
1019 const Vector<somenumber> &src,
1020 const number omega = 1.) const;
1021
1028 template <typename number2, typename number3>
1029 number
1031 const Vector<number2> &x,
1032 const Vector<number3> &b) const;
1033
1044 template <typename number2>
1045 void
1046 forward(Vector<number2> &dst, const Vector<number2> &src) const;
1047
1055 template <typename number2>
1056 void
1058
1060
1070
1076 number,
1077 << "The maximal pivot is " << arg1
1078 << ", which is below the threshold. The matrix may be singular.");
1083 size_type,
1084 size_type,
1085 size_type,
1086 << "Target region not in matrix: size in this direction="
1087 << arg1 << ", size of new matrix=" << arg2
1088 << ", offset=" << arg3);
1093 "You are attempting an operation on two matrices that "
1094 "are the same object, but the operation requires that the "
1095 "two objects are in fact different.");
1101};
1102
1105#ifndef DOXYGEN
1106/*-------------------------Inline functions -------------------------------*/
1107
1108
1109
1110template <typename number>
1111inline typename FullMatrix<number>::size_type
1113{
1114 return this->n_rows();
1115}
1116
1117
1118
1119template <typename number>
1120inline typename FullMatrix<number>::size_type
1122{
1123 return this->n_cols();
1124}
1125
1126
1127
1128template <typename number>
1130FullMatrix<number>::operator=(const number d)
1131{
1133 (void)d; // removes -Wunused-parameter warning in optimized mode
1134
1135 if (this->n_elements() != 0)
1136 this->reset_values();
1137
1138 return *this;
1139}
1140
1141
1142
1143template <typename number>
1144template <typename number2>
1145inline void
1146FullMatrix<number>::fill(const number2 *src)
1147{
1149}
1150
1151
1152
1153template <typename number>
1154template <typename MatrixType>
1155void
1156FullMatrix<number>::copy_from(const MatrixType &M)
1157{
1158 this->reinit(M.m(), M.n());
1159
1160 // loop over the elements of the argument matrix row by row, as suggested
1161 // in the documentation of the sparse matrix iterator class, and
1162 // copy them into the current object
1163 for (size_type row = 0; row < M.m(); ++row)
1164 {
1165 const typename MatrixType::const_iterator end_row = M.end(row);
1166 for (typename MatrixType::const_iterator entry = M.begin(row);
1167 entry != end_row;
1168 ++entry)
1169 this->el(row, entry->column()) = entry->value();
1170 }
1171}
1172
1173
1174
1175template <typename number>
1176template <typename MatrixType>
1177void
1178FullMatrix<number>::copy_transposed(const MatrixType &M)
1179{
1180 this->reinit(M.n(), M.m());
1181
1182 // loop over the elements of the argument matrix row by row, as suggested
1183 // in the documentation of the sparse matrix iterator class, and
1184 // copy them into the current object
1185 for (size_type row = 0; row < M.m(); ++row)
1186 {
1187 const typename MatrixType::const_iterator end_row = M.end(row);
1188 for (typename MatrixType::const_iterator entry = M.begin(row);
1189 entry != end_row;
1190 ++entry)
1191 this->el(entry->column(), row) = entry->value();
1192 }
1193}
1194
1195
1196
1197template <typename number>
1198template <typename MatrixType, typename index_type>
1199inline void
1201 const MatrixType & matrix,
1202 const std::vector<index_type> &row_index_set,
1203 const std::vector<index_type> &column_index_set)
1204{
1205 AssertDimension(row_index_set.size(), this->n_rows());
1206 AssertDimension(column_index_set.size(), this->n_cols());
1207
1208 const size_type n_rows_submatrix = row_index_set.size();
1209 const size_type n_cols_submatrix = column_index_set.size();
1210
1211 for (size_type sub_row = 0; sub_row < n_rows_submatrix; ++sub_row)
1212 for (size_type sub_col = 0; sub_col < n_cols_submatrix; ++sub_col)
1213 (*this)(sub_row, sub_col) =
1214 matrix.el(row_index_set[sub_row], column_index_set[sub_col]);
1215}
1216
1217
1218
1219template <typename number>
1220template <typename MatrixType, typename index_type>
1221inline void
1223 const std::vector<index_type> &row_index_set,
1224 const std::vector<index_type> &column_index_set,
1225 MatrixType & matrix) const
1226{
1227 AssertDimension(row_index_set.size(), this->n_rows());
1228 AssertDimension(column_index_set.size(), this->n_cols());
1229
1230 const size_type n_rows_submatrix = row_index_set.size();
1231 const size_type n_cols_submatrix = column_index_set.size();
1232
1233 for (size_type sub_row = 0; sub_row < n_rows_submatrix; ++sub_row)
1234 for (size_type sub_col = 0; sub_col < n_cols_submatrix; ++sub_col)
1235 matrix.set(row_index_set[sub_row],
1236 column_index_set[sub_col],
1237 (*this)(sub_row, sub_col));
1238}
1239
1240
1241template <typename number>
1242inline void
1243FullMatrix<number>::set(const size_type i,
1244 const size_type j,
1245 const number value)
1246{
1247 (*this)(i, j) = value;
1248}
1249
1250
1251
1252template <typename number>
1253template <typename number2>
1254void
1256 const Vector<number2> &v) const
1257{
1258 vmult(w, v, true);
1259}
1260
1261
1262template <typename number>
1263template <typename number2>
1264void
1266 const Vector<number2> &v) const
1267{
1268 Tvmult(w, v, true);
1269}
1270
1271
1272//---------------------------------------------------------------------------
1273template <typename number>
1274inline typename FullMatrix<number>::iterator
1275FullMatrix<number>::begin(const size_type r)
1276{
1277 AssertIndexRange(r, m());
1278 return iterator(this, r, 0);
1279}
1280
1281
1282
1283template <typename number>
1284inline typename FullMatrix<number>::iterator
1285FullMatrix<number>::end(const size_type r)
1286{
1287 AssertIndexRange(r, m());
1288 return iterator(this, r + 1, 0);
1289}
1290
1291
1292
1293template <typename number>
1295FullMatrix<number>::begin(const size_type r) const
1296{
1297 AssertIndexRange(r, m());
1298 return const_iterator(this, r, 0);
1299}
1300
1301
1302
1303template <typename number>
1305FullMatrix<number>::end(const size_type r) const
1306{
1307 AssertIndexRange(r, m());
1308 return const_iterator(this, r + 1, 0);
1309}
1310
1311
1312
1313template <typename number>
1314inline void
1315FullMatrix<number>::add(const size_type r, const size_type c, const number v)
1316{
1317 AssertIndexRange(r, this->m());
1318 AssertIndexRange(c, this->n());
1319
1320 this->operator()(r, c) += v;
1321}
1322
1323
1324
1325template <typename number>
1326template <typename number2, typename index_type>
1327inline void
1328FullMatrix<number>::add(const size_type row,
1329 const size_type n_cols,
1330 const index_type *col_indices,
1331 const number2 * values,
1332 const bool,
1333 const bool)
1334{
1335 AssertIndexRange(row, this->m());
1336 for (size_type col = 0; col < n_cols; ++col)
1337 {
1338 AssertIndexRange(col_indices[col], this->n());
1339 this->operator()(row, col_indices[col]) += values[col];
1340 }
1341}
1342
1343
1344template <typename number>
1345template <class StreamType>
1346inline void
1347FullMatrix<number>::print(StreamType & s,
1348 const unsigned int w,
1349 const unsigned int p) const
1350{
1351 Assert(!this->empty(), ExcEmptyMatrix());
1352
1353 // save the state of out stream
1354 const std::streamsize old_precision = s.precision(p);
1355 const std::streamsize old_width = s.width(w);
1356
1357 for (size_type i = 0; i < this->m(); ++i)
1358 {
1359 for (size_type j = 0; j < this->n(); ++j)
1360 {
1361 s.width(w);
1362 s.precision(p);
1363 s << this->el(i, j);
1364 }
1365 s << std::endl;
1366 }
1367
1368 // reset output format
1369 s.precision(old_precision);
1370 s.width(old_width);
1371}
1372
1373
1374#endif // DOXYGEN
1375
1377
1378#endif
typename numbers::NumberTraits< number >::real_type real_type
Definition: full_matrix.h:118
typename Table< 2, number >::const_iterator const_iterator
Definition: full_matrix.h:97
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.))
FullMatrix< number > & operator=(const number d)
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 symmetrize()
void equ(const number a, const FullMatrix< number2 > &A, const number b, const FullMatrix< number2 > &B)
number trace() const
void right_invert(const FullMatrix< number2 > &M)
FullMatrix< number > & operator=(const FullMatrix< number2 > &)
FullMatrix & operator/=(const number factor)
std::size_t size_type
Definition: full_matrix.h:81
void set(const size_type i, const size_type j, const number value)
size_type n() const
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 IdentityMatrix &id)
number value_type
Definition: full_matrix.h:87
typename Table< 2, number >::iterator iterator
Definition: full_matrix.h:92
void add(const number a, const FullMatrix< number2 > &A, const number b, const FullMatrix< number2 > &B, const number c, const FullMatrix< number2 > &C)
void Tadd(const number s, const FullMatrix< number2 > &B)
void vmult(Vector< number2 > &w, const Vector< number2 > &v, const bool adding=false) const
bool all_zero() 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 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 gauss_jordan()
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
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)
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
FullMatrix< number > & operator=(const LAPACKFullMatrix< number2 > &)
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
FullMatrix & operator*=(const number factor)
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(const size_type n=0)
size_type m() const
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
Definition: tensor.h:503
Definition: vector.h:109
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:442
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:443
#define DeclException0(Exception0)
Definition: exceptions.h:464
static ::ExceptionBase & ExcEmptyMatrix()
static ::ExceptionBase & ExcScalarAssignmentOnlyForZeroValue()
#define Assert(cond, exc)
Definition: exceptions.h:1473
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1667
static ::ExceptionBase & ExcNotRegular(number arg1)
#define AssertIndexRange(index, range)
Definition: exceptions.h:1732
#define DeclExceptionMsg(Exception, defaulttext)
Definition: exceptions.h:487
static ::ExceptionBase & ExcMatrixNotPositiveDefinite()
static ::ExceptionBase & ExcSourceEqualsDestination()
#define DeclException3(Exception3, type1, type2, type3, outsequence)
Definition: exceptions.h:555
static ::ExceptionBase & ExcInvalidDestination(size_type arg1, size_type arg2, size_type arg3)
#define DeclException1(Exception1, type1, outsequence)
Definition: exceptions.h:509
@ matrix
Contents is actually a matrix.
types::global_dof_index size_type
Definition: cuda_kernels.h:45
void reinit(MatrixBlock< MatrixType > &v, const BlockSparsityPattern &p)
Definition: matrix_block.h:618