Reference documentation for deal.II version Git 6b6ce30 2018-10-20 08:13:43 -0400
\(\newcommand{\dealcoloneq}{\mathrel{\vcenter{:}}=}\)
full_matrix.h
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 1999 - 2018 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 
22 #include <deal.II/base/numbers.h>
23 #include <deal.II/base/table.h>
24 #include <deal.II/base/tensor.h>
25 
26 #include <deal.II/differentiation/ad/ad_number_traits.h>
27 
28 #include <deal.II/lac/exceptions.h>
29 #include <deal.II/lac/identity_matrix.h>
30 
31 #include <cstring>
32 #include <iomanip>
33 #include <vector>
34 
35 DEAL_II_NAMESPACE_OPEN
36 
37 
38 // forward declarations
39 template <typename number>
40 class Vector;
41 template <typename number>
43 
44 
69 template <typename number>
70 class FullMatrix : public Table<2, number>
71 {
72 public:
73  // The assertion in full_matrix.templates.h for whether or not a number is
74  // finite is not compatible for AD number types.
75  static_assert(
77  "The FullMatrix class does not support auto-differentiable numbers.");
78 
82  using size_type = std::size_t;
83 
88  using value_type = number;
89 
94 
99 
104 
108  using Table<2, number>::end;
109 
120 
125 
135  explicit FullMatrix(const size_type n = 0);
136 
140  FullMatrix(const size_type rows, const size_type cols);
141 
146  FullMatrix(const size_type rows, const size_type cols, const number *entries);
147 
156  FullMatrix(const IdentityMatrix &id);
171  template <typename number2>
174 
183  operator=(const number d);
184 
194  operator=(const IdentityMatrix &id);
195 
200  template <typename number2>
203 
204 
210  template <typename MatrixType>
211  void
212  copy_from(const MatrixType &);
213 
219  template <typename MatrixType>
220  void
221  copy_transposed(const MatrixType &);
222 
230  template <int dim>
231  void
232  copy_from(const Tensor<2, dim> &T,
233  const unsigned int src_r_i = 0,
234  const unsigned int src_r_j = dim - 1,
235  const unsigned int src_c_i = 0,
236  const unsigned int src_c_j = dim - 1,
237  const size_type dst_r = 0,
238  const size_type dst_c = 0);
239 
247  template <int dim>
248  void copy_to(Tensor<2, dim> & T,
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
304  fill(const FullMatrix<number2> &src,
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
363  operator==(const FullMatrix<number> &) const;
364 
369  size_type
370  m() const;
371 
376  size_type
377  n() const;
378 
384  bool
385  all_zero() const;
386 
402  template <typename number2>
403  number2
404  matrix_norm_square(const Vector<number2> &v) const;
405 
415  template <typename number2>
416  number2
417  matrix_scalar_product(const Vector<number2> &u,
418  const Vector<number2> &v) const;
419 
424  real_type
425  l1_norm() const;
426 
431  real_type
432  linfty_norm() const;
433 
441  real_type
442  frobenius_norm() const;
443 
452  real_type
453  relative_symmetry_norm2() const;
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
519  memory_consumption() const;
520 
522 
524 
528  iterator
529  begin(const size_type r);
530 
534  iterator
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
620  add(const FullMatrix<number2> &src,
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
649  Tadd(const FullMatrix<number2> &src,
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
691  add_row(const size_type i,
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
708  add_col(const size_type i,
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
769  symmetrize();
770 
785  void
786  gauss_jordan();
787 
794  template <typename number2>
795  void
796  invert(const FullMatrix<number2> &M);
797 
806  template <typename number2>
807  void
808  cholesky(const FullMatrix<number2> &A);
809 
814  template <typename number2>
815  void
816  outer_product(const Vector<number2> &V, const Vector<number2> &W);
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
969  vmult(Vector<number2> & w,
970  const Vector<number2> &v,
971  const bool adding = false) const;
972 
978  template <typename number2>
979  void
980  vmult_add(Vector<number2> &w, const Vector<number2> &v) const;
981 
995  template <typename number2>
996  void
997  Tvmult(Vector<number2> & w,
998  const Vector<number2> &v,
999  const bool adding = false) const;
1000 
1007  template <typename number2>
1008  void
1009  Tvmult_add(Vector<number2> &w, const Vector<number2> &v) const;
1010 
1016  template <typename somenumber>
1017  void
1018  precondition_Jacobi(Vector<somenumber> & dst,
1019  const Vector<somenumber> &src,
1020  const number omega = 1.) const;
1021 
1028  template <typename number2, typename number3>
1029  number
1030  residual(Vector<number2> & dst,
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
1057  backward(Vector<number2> &dst, const Vector<number2> &src) const;
1058 
1060 
1070 
1075  ExcNotRegular,
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 
1110 template <typename number>
1111 inline typename FullMatrix<number>::size_type
1112 FullMatrix<number>::m() const
1113 {
1114  return this->n_rows();
1115 }
1116 
1117 
1118 
1119 template <typename number>
1120 inline typename FullMatrix<number>::size_type
1121 FullMatrix<number>::n() const
1122 {
1123  return this->n_cols();
1124 }
1125 
1126 
1127 
1128 template <typename number>
1130 FullMatrix<number>::operator=(const number d)
1131 {
1132  Assert(d == number(0), ExcScalarAssignmentOnlyForZeroValue());
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 
1143 template <typename number>
1144 template <typename number2>
1145 inline void
1146 FullMatrix<number>::fill(const number2 *src)
1147 {
1149 }
1150 
1151 
1152 
1153 template <typename number>
1154 template <typename MatrixType>
1155 void
1156 FullMatrix<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 
1175 template <typename number>
1176 template <typename MatrixType>
1177 void
1178 FullMatrix<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 
1197 template <typename number>
1198 template <typename MatrixType, typename index_type>
1199 inline 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 
1219 template <typename number>
1220 template <typename MatrixType, typename index_type>
1221 inline 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 
1241 template <typename number>
1242 inline void
1243 FullMatrix<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 
1252 template <typename number>
1253 template <typename number2>
1254 void
1255 FullMatrix<number>::vmult_add(Vector<number2> & w,
1256  const Vector<number2> &v) const
1257 {
1258  vmult(w, v, true);
1259 }
1260 
1261 
1262 template <typename number>
1263 template <typename number2>
1264 void
1265 FullMatrix<number>::Tvmult_add(Vector<number2> & w,
1266  const Vector<number2> &v) const
1267 {
1268  Tvmult(w, v, true);
1269 }
1270 
1271 
1272 //---------------------------------------------------------------------------
1273 template <typename number>
1274 inline typename FullMatrix<number>::iterator
1275 FullMatrix<number>::begin(const size_type r)
1276 {
1277  AssertIndexRange(r, m());
1278  return iterator(this, r, 0);
1279 }
1280 
1281 
1282 
1283 template <typename number>
1284 inline typename FullMatrix<number>::iterator
1285 FullMatrix<number>::end(const size_type r)
1286 {
1287  AssertIndexRange(r, m());
1288  return iterator(this, r + 1, 0);
1289 }
1290 
1291 
1292 
1293 template <typename number>
1294 inline typename FullMatrix<number>::const_iterator
1295 FullMatrix<number>::begin(const size_type r) const
1296 {
1297  AssertIndexRange(r, m());
1298  return const_iterator(this, r, 0);
1299 }
1300 
1301 
1302 
1303 template <typename number>
1304 inline typename FullMatrix<number>::const_iterator
1305 FullMatrix<number>::end(const size_type r) const
1306 {
1307  AssertIndexRange(r, m());
1308  return const_iterator(this, r + 1, 0);
1309 }
1310 
1311 
1312 
1313 template <typename number>
1314 inline void
1315 FullMatrix<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 
1325 template <typename number>
1326 template <typename number2, typename index_type>
1327 inline void
1328 FullMatrix<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 
1344 template <typename number>
1345 template <class StreamType>
1346 inline void
1347 FullMatrix<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 
1376 DEAL_II_NAMESPACE_CLOSE
1377 
1378 #endif
void vmult_add(Vector< number2 > &w, const Vector< number2 > &v) const
void diagadd(const number s)
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
typename Table< 2, std::complex< double > >::const_iterator const_iterator
Definition: full_matrix.h:98
static::ExceptionBase & ExcEmptyMatrix()
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1366
void TmTmult(FullMatrix< number2 > &C, const FullMatrix< number2 > &B, const bool adding=false) const
real_type linfty_norm() const
typename numbers::NumberTraits< std::complex< double > >::real_type real_type
Definition: full_matrix.h:119
void backward(Vector< number2 > &dst, const Vector< number2 > &src) const
FullMatrix(const size_type n=0)
FullMatrix & operator/=(const number factor)
std::size_t size_type
Definition: full_matrix.h:82
void gauss_jordan()
static::ExceptionBase & ExcNotRegular(number arg1)
static::ExceptionBase & ExcScalarAssignmentOnlyForZeroValue()
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)
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)
Definition: exceptions.h:1407
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)
number trace() const
void vmult(Vector< number2 > &w, const Vector< number2 > &v, const bool adding=false) const
void left_invert(const FullMatrix< number2 > &M)
void Tadd(const number s, const FullMatrix< number2 > &B)
void copy_transposed(const MatrixType &)
void scatter_matrix_to(const std::vector< index_type > &row_index_set, const std::vector< index_type > &column_index_set, MatrixType &matrix) const
real_type relative_symmetry_norm2() const
FullMatrix & operator*=(const number factor)
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
number2 matrix_scalar_product(const Vector< number2 > &u, const Vector< number2 > &v) const
std::size_t memory_consumption() const
size_type n() const
void swap_row(const size_type i, const size_type j)
#define DeclException1(Exception1, type1, outsequence)
Definition: exceptions.h:408
number2 matrix_norm_square(const Vector< number2 > &v) const
void invert(const FullMatrix< number2 > &M)
#define Assert(cond, exc)
Definition: exceptions.h:1239
#define DeclExceptionMsg(Exception, defaulttext)
Definition: exceptions.h:397
#define DeclException0(Exception0)
Definition: exceptions.h:385
void mTmult(FullMatrix< number2 > &C, const FullMatrix< number2 > &B, const bool adding=false) const
void extract_submatrix_from(const MatrixType &matrix, const std::vector< index_type > &row_index_set, const std::vector< index_type > &column_index_set)
void Tvmult(Vector< number2 > &w, const Vector< number2 > &v, const bool adding=false) const
iterator end(const size_type r)
void mmult(FullMatrix< number2 > &C, const FullMatrix< number2 > &B, 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
size_type m() const
FullMatrix< number > & operator=(const FullMatrix< number2 > &)
void Tmmult(FullMatrix< number2 > &C, const FullMatrix< number2 > &B, const bool adding=false) const
void add_col(const size_type i, const number s, const size_type j)
static::ExceptionBase & ExcSourceEqualsDestination()
void swap_col(const size_type i, const size_type j)
void copy_from(const MatrixType &)
typename AlignedVector< T >::size_type size_type
Definition: table.h:429
real_type l1_norm() const
void add(const number a, const FullMatrix< number2 > &A)
static::ExceptionBase & ExcInvalidDestination(size_type arg1, size_type arg2, size_type arg3)
real_type frobenius_norm() const
Definition: table.h:37
void symmetrize()
#define DeclException3(Exception3, type1, type2, type3, outsequence)
Definition: exceptions.h:432
number determinant() const
static::ExceptionBase & ExcMatrixNotPositiveDefinite()
iterator begin(const size_type r)
AlignedVector< T > values
Definition: table.h:681
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 fill(InputIterator entries, const bool C_style_indexing=true)
void print(StreamType &s, const unsigned int width=5, const unsigned int precision=2) const
bool all_zero() const
void precondition_Jacobi(Vector< somenumber > &dst, const Vector< somenumber > &src, const number omega=1.) const
typename Table< 2, std::complex< double > >::iterator iterator
Definition: full_matrix.h:93