Reference documentation for deal.II version Git d2d482dcec 2020-10-26 10:29:03 -0400
\(\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\}}\)
sparse_matrix_ez.h
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2002 - 2020 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_sparse_matrix_ez_h
17 # define dealii_sparse_matrix_ez_h
18 
19 
20 # include <deal.II/base/config.h>
21 
24 
25 # include <deal.II/lac/exceptions.h>
26 
27 # include <vector>
28 
30 
31 // Forward declarations
32 # ifndef DOXYGEN
33 template <typename number>
34 class Vector;
35 template <typename number>
36 class FullMatrix;
37 # endif
38 
103 template <typename number>
105 {
106 public:
111 
116  struct Entry
117  {
121  Entry();
122 
126  Entry(const size_type column, const number &value);
127 
132 
136  number value;
137 
142  };
143 
148  struct RowInfo
149  {
153  RowInfo(const size_type start = Entry::invalid);
154 
162  unsigned short length;
166  unsigned short diagonal;
170  static const unsigned short invalid_diagonal =
171  static_cast<unsigned short>(-1);
172  };
173 
174 public:
179  {
180  private:
184  class Accessor
185  {
186  public:
192  const size_type row,
193  const unsigned short index);
194 
198  size_type
199  row() const;
200 
204  unsigned short
205  index() const;
206 
210  size_type
211  column() const;
212 
216  number
217  value() const;
218 
219  protected:
224 
229 
233  unsigned short a_index;
234 
235  // Make enclosing class a friend.
236  friend class const_iterator;
237  };
238 
239  public:
244  const size_type row,
245  const unsigned short index);
246 
251  operator++();
252 
256  const Accessor &operator*() const;
257 
261  const Accessor *operator->() const;
262 
266  bool
267  operator==(const const_iterator &) const;
271  bool
272  operator!=(const const_iterator &) const;
273 
278  bool
279  operator<(const const_iterator &) const;
280 
281  private:
286  };
287 
292  using value_type = number;
293 
301  SparseMatrixEZ();
302 
311 
318  explicit SparseMatrixEZ(const size_type n_rows,
319  const size_type n_columns,
320  const size_type default_row_length = 0,
321  const unsigned int default_increment = 1);
322 
326  ~SparseMatrixEZ() override = default;
327 
333 
343  operator=(const double d);
344 
352  void
353  reinit(const size_type n_rows,
354  const size_type n_columns,
355  size_type default_row_length = 0,
356  unsigned int default_increment = 1,
357  size_type reserve = 0);
358 
363  void
364  clear();
366 
374  bool
375  empty() const;
376 
381  size_type
382  m() const;
383 
388  size_type
389  n() const;
390 
394  size_type
395  get_row_length(const size_type row) const;
396 
400  size_type
401  n_nonzero_elements() const;
402 
407  std::size_t
408  memory_consumption() const;
409 
415  template <class StreamType>
416  void
417  print_statistics(StreamType &s, bool full = false);
418 
428  void
430  size_type & allocated,
431  size_type & reserved,
432  std::vector<size_type> &used_by_line,
433  const bool compute_by_line) const;
435 
456  void
457  set(const size_type i,
458  const size_type j,
459  const number value,
460  const bool elide_zero_values = true);
461 
472  void
473  add(const size_type i, const size_type j, const number value);
474 
489  template <typename number2>
490  void
491  add(const std::vector<size_type> &indices,
492  const FullMatrix<number2> & full_matrix,
493  const bool elide_zero_values = true);
494 
500  template <typename number2>
501  void
502  add(const std::vector<size_type> &row_indices,
503  const std::vector<size_type> &col_indices,
504  const FullMatrix<number2> & full_matrix,
505  const bool elide_zero_values = true);
506 
516  template <typename number2>
517  void
518  add(const size_type row,
519  const std::vector<size_type> &col_indices,
520  const std::vector<number2> & values,
521  const bool elide_zero_values = true);
522 
532  template <typename number2>
533  void
534  add(const size_type row,
535  const size_type n_cols,
536  const size_type *col_indices,
537  const number2 * values,
538  const bool elide_zero_values = true,
539  const bool col_indices_are_sorted = false);
540 
562  template <typename MatrixType>
564  copy_from(const MatrixType &source, const bool elide_zero_values = true);
565 
573  template <typename MatrixType>
574  void
575  add(const number factor, const MatrixType &matrix);
577 
590  number
591  operator()(const size_type i, const size_type j) const;
592 
597  number
598  el(const size_type i, const size_type j) const;
600 
608  template <typename somenumber>
609  void
610  vmult(Vector<somenumber> &dst, const Vector<somenumber> &src) const;
611 
617  template <typename somenumber>
618  void
619  Tvmult(Vector<somenumber> &dst, const Vector<somenumber> &src) const;
620 
625  template <typename somenumber>
626  void
627  vmult_add(Vector<somenumber> &dst, const Vector<somenumber> &src) const;
628 
634  template <typename somenumber>
635  void
636  Tvmult_add(Vector<somenumber> &dst, const Vector<somenumber> &src) const;
638 
645  number
646  l2_norm() const;
648 
657  template <typename somenumber>
658  void
659  precondition_Jacobi(Vector<somenumber> & dst,
660  const Vector<somenumber> &src,
661  const number omega = 1.) const;
662 
666  template <typename somenumber>
667  void
668  precondition_SSOR(Vector<somenumber> & dst,
669  const Vector<somenumber> & src,
670  const number om = 1.,
671  const std::vector<std::size_t> &pos_right_of_diagonal =
672  std::vector<std::size_t>()) const;
673 
678  template <typename somenumber>
679  void
680  precondition_SOR(Vector<somenumber> & dst,
681  const Vector<somenumber> &src,
682  const number om = 1.) const;
683 
688  template <typename somenumber>
689  void
690  precondition_TSOR(Vector<somenumber> & dst,
691  const Vector<somenumber> &src,
692  const number om = 1.) const;
693 
702  template <typename MatrixTypeA, typename MatrixTypeB>
703  void
704  conjugate_add(const MatrixTypeA &A,
705  const MatrixTypeB &B,
706  const bool transpose = false);
708 
715  const_iterator
716  begin() const;
717 
721  const_iterator
722  end() const;
723 
728  const_iterator
729  begin(const size_type r) const;
730 
735  const_iterator
736  end(const size_type r) const;
738 
746  void
747  print(std::ostream &out) const;
748 
769  void
770  print_formatted(std::ostream & out,
771  const unsigned int precision = 3,
772  const bool scientific = true,
773  const unsigned int width = 0,
774  const char * zero_string = " ",
775  const double denominator = 1.) const;
776 
782  void
783  block_write(std::ostream &out) const;
784 
795  void
796  block_read(std::istream &in);
798 
808 
813  int,
814  int,
815  << "The entry with index (" << arg1 << ',' << arg2
816  << ") does not exist.");
817 
819  int,
820  int,
821  << "An entry with index (" << arg1 << ',' << arg2
822  << ") cannot be allocated.");
824 private:
829  const Entry *
830  locate(const size_type row, const size_type col) const;
831 
836  Entry *
837  locate(const size_type row, const size_type col);
838 
842  Entry *
843  allocate(const size_type row, const size_type col);
844 
850  template <typename somenumber>
851  void
852  threaded_vmult(Vector<somenumber> & dst,
853  const Vector<somenumber> &src,
854  const size_type begin_row,
855  const size_type end_row) const;
856 
862  template <typename somenumber>
863  void
864  threaded_matrix_norm_square(const Vector<somenumber> &v,
865  const size_type begin_row,
866  const size_type end_row,
867  somenumber * partial_sum) const;
868 
874  template <typename somenumber>
875  void
876  threaded_matrix_scalar_product(const Vector<somenumber> &u,
877  const Vector<somenumber> &v,
878  const size_type begin_row,
879  const size_type end_row,
880  somenumber * partial_sum) const;
881 
886 
890  std::vector<RowInfo> row_info;
891 
895  std::vector<Entry> data;
896 
900  unsigned int increment;
901 
906 };
907 
911 /*---------------------- Inline functions -----------------------------------*/
912 
913 template <typename number>
915  const number & value)
916  : column(column)
917  , value(value)
918 {}
919 
920 
921 
922 template <typename number>
924  : column(invalid)
925  , value(0)
926 {}
927 
928 
929 template <typename number>
931  : start(start)
932  , length(0)
933  , diagonal(invalid_diagonal)
934 {}
935 
936 
937 //---------------------------------------------------------------------------
938 template <typename number>
941  const size_type r,
942  const unsigned short i)
943  : matrix(matrix)
944  , a_row(r)
945  , a_index(i)
946 {}
947 
948 
949 template <typename number>
950 inline typename SparseMatrixEZ<number>::size_type
952 {
953  return a_row;
954 }
955 
956 
957 template <typename number>
958 inline typename SparseMatrixEZ<number>::size_type
960 {
961  return matrix->data[matrix->row_info[a_row].start + a_index].column;
962 }
963 
964 
965 template <typename number>
966 inline unsigned short
968 {
969  return a_index;
970 }
971 
972 
973 
974 template <typename number>
975 inline number
977 {
978  return matrix->data[matrix->row_info[a_row].start + a_index].value;
979 }
980 
981 
982 template <typename number>
985  const size_type r,
986  const unsigned short i)
987  : accessor(matrix, r, i)
988 {
989  // Finish if this is the end()
990  if (r == accessor.matrix->m() && i == 0)
991  return;
992 
993  // Make sure we never construct an
994  // iterator pointing to a
995  // non-existing entry
996 
997  // If the index points beyond the
998  // end of the row, try the next
999  // row.
1000  if (accessor.a_index >= accessor.matrix->row_info[accessor.a_row].length)
1001  {
1002  do
1003  {
1004  ++accessor.a_row;
1005  }
1006  // Beware! If the next row is
1007  // empty, iterate until a
1008  // non-empty row is found or we
1009  // hit the end of the matrix.
1010  while (accessor.a_row < accessor.matrix->m() &&
1011  accessor.matrix->row_info[accessor.a_row].length == 0);
1012  }
1013 }
1014 
1015 
1016 template <typename number>
1019 {
1021 
1022  // Increment column index
1023  ++(accessor.a_index);
1024  // If index exceeds number of
1025  // entries in this row, proceed
1026  // with next row.
1027  if (accessor.a_index >= accessor.matrix->row_info[accessor.a_row].length)
1028  {
1029  accessor.a_index = 0;
1030  // Do this loop to avoid
1031  // elements in empty rows
1032  do
1033  {
1034  ++accessor.a_row;
1035  }
1036  while (accessor.a_row < accessor.matrix->m() &&
1037  accessor.matrix->row_info[accessor.a_row].length == 0);
1038  }
1039  return *this;
1040 }
1041 
1042 
1043 template <typename number>
1046 {
1047  return accessor;
1048 }
1049 
1050 
1051 template <typename number>
1054 {
1055  return &accessor;
1056 }
1057 
1058 
1059 template <typename number>
1060 inline bool
1062 operator==(const const_iterator &other) const
1063 {
1064  return (accessor.row() == other.accessor.row() &&
1065  accessor.index() == other.accessor.index());
1066 }
1067 
1068 
1069 template <typename number>
1070 inline bool
1072 operator!=(const const_iterator &other) const
1073 {
1074  return !(*this == other);
1075 }
1076 
1077 
1078 template <typename number>
1079 inline bool
1081 operator<(const const_iterator &other) const
1082 {
1083  return (accessor.row() < other.accessor.row() ||
1084  (accessor.row() == other.accessor.row() &&
1085  accessor.index() < other.accessor.index()));
1086 }
1087 
1088 
1089 //---------------------------------------------------------------------------
1090 template <typename number>
1091 inline typename SparseMatrixEZ<number>::size_type
1093 {
1094  return row_info.size();
1095 }
1096 
1097 
1098 template <typename number>
1099 inline typename SparseMatrixEZ<number>::size_type
1101 {
1102  return n_columns;
1103 }
1104 
1105 
1106 template <typename number>
1107 inline typename SparseMatrixEZ<number>::Entry *
1109 {
1110  AssertIndexRange(row, m());
1111  AssertIndexRange(col, n());
1112 
1113  const RowInfo & r = row_info[row];
1114  const size_type end = r.start + r.length;
1115  for (size_type i = r.start; i < end; ++i)
1116  {
1117  Entry *const entry = &data[i];
1118  if (entry->column == col)
1119  return entry;
1120  if (entry->column == Entry::invalid)
1121  return nullptr;
1122  }
1123  return nullptr;
1124 }
1125 
1126 
1127 
1128 template <typename number>
1129 inline const typename SparseMatrixEZ<number>::Entry *
1131 {
1132  SparseMatrixEZ<number> *t = const_cast<SparseMatrixEZ<number> *>(this);
1133  return t->locate(row, col);
1134 }
1135 
1136 
1137 template <typename number>
1138 inline typename SparseMatrixEZ<number>::Entry *
1140 {
1141  AssertIndexRange(row, m());
1142  AssertIndexRange(col, n());
1143 
1144  RowInfo & r = row_info[row];
1145  const size_type end = r.start + r.length;
1146 
1147  size_type i = r.start;
1148  // If diagonal exists and this
1149  // column is higher, start only
1150  // after diagonal.
1151  if (r.diagonal != RowInfo::invalid_diagonal && col >= row)
1152  i += r.diagonal;
1153  // Find position of entry
1154  while (i < end && data[i].column < col)
1155  ++i;
1156 
1157  // entry found
1158  if (i != end && data[i].column == col)
1159  return &data[i];
1160 
1161  // Now, we must insert the new
1162  // entry and move all successive
1163  // entries back.
1164 
1165  // If no more space is available
1166  // for this row, insert new
1167  // elements into the vector.
1168  // TODO:[GK] We should not extend this row if i<end
1169  if (row != row_info.size() - 1)
1170  {
1171  if (end >= row_info[row + 1].start)
1172  {
1173  // Failure if increment 0
1174  Assert(increment != 0, ExcEntryAllocationFailure(row, col));
1175 
1176  // Insert new entries
1177  data.insert(data.begin() + end, increment, Entry());
1178  // Update starts of
1179  // following rows
1180  for (size_type rn = row + 1; rn < row_info.size(); ++rn)
1181  row_info[rn].start += increment;
1182  }
1183  }
1184  else
1185  {
1186  if (end >= data.size())
1187  {
1188  // Here, appending a block
1189  // does not increase
1190  // performance.
1191  data.push_back(Entry());
1192  }
1193  }
1194 
1195  Entry *entry = &data[i];
1196  // Save original entry
1197  Entry temp = *entry;
1198  // Insert new entry here to
1199  // make sure all entries
1200  // are ordered by column
1201  // index
1202  entry->column = col;
1203  entry->value = 0;
1204  // Update row_info
1205  ++r.length;
1206  if (col == row)
1207  r.diagonal = i - r.start;
1208  else if (col < row && r.diagonal != RowInfo::invalid_diagonal)
1209  ++r.diagonal;
1210 
1211  if (i == end)
1212  return entry;
1213 
1214  // Move all entries in this
1215  // row up by one
1216  for (size_type j = i + 1; j < end; ++j)
1217  {
1218  // There should be no invalid
1219  // entry below end
1220  Assert(data[j].column != Entry::invalid, ExcInternalError());
1221 
1222  // TODO[GK]: This could be done more efficiently by moving starting at the
1223  // top rather than swapping starting at the bottom
1224  std::swap(data[j], temp);
1225  }
1226  Assert(data[end].column == Entry::invalid, ExcInternalError());
1227 
1228  data[end] = temp;
1229 
1230  return entry;
1231 }
1232 
1233 
1234 
1235 template <typename number>
1236 inline void
1238  const size_type j,
1239  const number value,
1240  const bool elide_zero_values)
1241 {
1242  AssertIsFinite(value);
1243 
1244  AssertIndexRange(i, m());
1245  AssertIndexRange(j, n());
1246 
1247  if (elide_zero_values && value == 0.)
1248  {
1249  Entry *entry = locate(i, j);
1250  if (entry != nullptr)
1251  entry->value = 0.;
1252  }
1253  else
1254  {
1255  Entry *entry = allocate(i, j);
1256  entry->value = value;
1257  }
1258 }
1259 
1260 
1261 
1262 template <typename number>
1263 inline void
1265  const size_type j,
1266  const number value)
1267 {
1268  AssertIsFinite(value);
1269 
1270  AssertIndexRange(i, m());
1271  AssertIndexRange(j, n());
1272 
1273  // ignore zero additions
1274  if (std::abs(value) == 0.)
1275  return;
1276 
1277  Entry *entry = allocate(i, j);
1278  entry->value += value;
1279 }
1280 
1281 
1282 template <typename number>
1283 template <typename number2>
1284 void
1285 SparseMatrixEZ<number>::add(const std::vector<size_type> &indices,
1286  const FullMatrix<number2> & full_matrix,
1287  const bool elide_zero_values)
1288 {
1289  // TODO: This function can surely be made more efficient
1290  for (size_type i = 0; i < indices.size(); ++i)
1291  for (size_type j = 0; j < indices.size(); ++j)
1292  if ((full_matrix(i, j) != 0) || (elide_zero_values == false))
1293  add(indices[i], indices[j], full_matrix(i, j));
1294 }
1295 
1296 
1297 
1298 template <typename number>
1299 template <typename number2>
1300 void
1301 SparseMatrixEZ<number>::add(const std::vector<size_type> &row_indices,
1302  const std::vector<size_type> &col_indices,
1303  const FullMatrix<number2> & full_matrix,
1304  const bool elide_zero_values)
1305 {
1306  // TODO: This function can surely be made more efficient
1307  for (size_type i = 0; i < row_indices.size(); ++i)
1308  for (size_type j = 0; j < col_indices.size(); ++j)
1309  if ((full_matrix(i, j) != 0) || (elide_zero_values == false))
1310  add(row_indices[i], col_indices[j], full_matrix(i, j));
1311 }
1312 
1313 
1314 
1315 template <typename number>
1316 template <typename number2>
1317 void
1319  const std::vector<size_type> &col_indices,
1320  const std::vector<number2> & values,
1321  const bool elide_zero_values)
1322 {
1323  // TODO: This function can surely be made more efficient
1324  for (size_type j = 0; j < col_indices.size(); ++j)
1325  if ((values[j] != 0) || (elide_zero_values == false))
1326  add(row, col_indices[j], values[j]);
1327 }
1328 
1329 
1330 
1331 template <typename number>
1332 template <typename number2>
1333 void
1335  const size_type n_cols,
1336  const size_type *col_indices,
1337  const number2 * values,
1338  const bool elide_zero_values,
1339  const bool /*col_indices_are_sorted*/)
1340 {
1341  // TODO: This function can surely be made more efficient
1342  for (size_type j = 0; j < n_cols; ++j)
1343  if ((std::abs(values[j]) != 0) || (elide_zero_values == false))
1344  add(row, col_indices[j], values[j]);
1345 }
1346 
1347 
1348 
1349 template <typename number>
1350 inline number
1352 {
1353  const Entry *entry = locate(i, j);
1354  if (entry)
1355  return entry->value;
1356  return 0.;
1357 }
1358 
1359 
1360 
1361 template <typename number>
1362 inline number
1364 {
1365  const Entry *entry = locate(i, j);
1366  if (entry)
1367  return entry->value;
1368  Assert(false, ExcInvalidEntry(i, j));
1369  return 0.;
1370 }
1371 
1372 
1373 template <typename number>
1376 {
1377  const_iterator result(this, 0, 0);
1378  return result;
1379 }
1380 
1381 template <typename number>
1384 {
1385  return const_iterator(this, m(), 0);
1386 }
1387 
1388 template <typename number>
1391 {
1392  AssertIndexRange(r, m());
1393  const_iterator result(this, r, 0);
1394  return result;
1395 }
1396 
1397 template <typename number>
1400 {
1401  AssertIndexRange(r, m());
1402  const_iterator result(this, r + 1, 0);
1403  return result;
1404 }
1405 
1406 template <typename number>
1407 template <typename MatrixType>
1408 inline SparseMatrixEZ<number> &
1410  const bool elide_zero_values)
1411 {
1412  reinit(M.m(), M.n(), this->saved_default_row_length, this->increment);
1413 
1414  // loop over the elements of the argument matrix row by row, as suggested
1415  // in the documentation of the sparse matrix iterator class, and
1416  // copy them into the current object
1417  for (size_type row = 0; row < M.m(); ++row)
1418  {
1419  const typename MatrixType::const_iterator end_row = M.end(row);
1420  for (typename MatrixType::const_iterator entry = M.begin(row);
1421  entry != end_row;
1422  ++entry)
1423  set(row, entry->column(), entry->value(), elide_zero_values);
1424  }
1425 
1426  return *this;
1427 }
1428 
1429 template <typename number>
1430 template <typename MatrixType>
1431 inline void
1432 SparseMatrixEZ<number>::add(const number factor, const MatrixType &M)
1433 {
1434  Assert(M.m() == m(), ExcDimensionMismatch(M.m(), m()));
1435  Assert(M.n() == n(), ExcDimensionMismatch(M.n(), n()));
1436 
1437  if (factor == 0.)
1438  return;
1439 
1440  // loop over the elements of the argument matrix row by row, as suggested
1441  // in the documentation of the sparse matrix iterator class, and
1442  // add them into the current object
1443  for (size_type row = 0; row < M.m(); ++row)
1444  {
1445  const typename MatrixType::const_iterator end_row = M.end(row);
1446  for (typename MatrixType::const_iterator entry = M.begin(row);
1447  entry != end_row;
1448  ++entry)
1449  if (entry->value() != 0)
1450  add(row, entry->column(), factor * entry->value());
1451  }
1452 }
1453 
1454 
1455 
1456 template <typename number>
1457 template <typename MatrixTypeA, typename MatrixTypeB>
1458 inline void
1460  const MatrixTypeB &B,
1461  const bool transpose)
1462 {
1463  // Compute the result
1464  // r_ij = \sum_kl b_ik b_jl a_kl
1465 
1466  // Assert (n() == B.m(), ExcDimensionMismatch(n(), B.m()));
1467  // Assert (m() == B.m(), ExcDimensionMismatch(m(), B.m()));
1468  // Assert (A.n() == B.n(), ExcDimensionMismatch(A.n(), B.n()));
1469  // Assert (A.m() == B.n(), ExcDimensionMismatch(A.m(), B.n()));
1470 
1471  // Somehow, we have to avoid making
1472  // this an operation of complexity
1473  // n^2. For the transpose case, we
1474  // can go through the non-zero
1475  // elements of A^-1 and use the
1476  // corresponding rows of B only.
1477  // For the non-transpose case, we
1478  // must find a trick.
1479  typename MatrixTypeB::const_iterator b1 = B.begin();
1480  const typename MatrixTypeB::const_iterator b_final = B.end();
1481  if (transpose)
1482  while (b1 != b_final)
1483  {
1484  const size_type i = b1->column();
1485  const size_type k = b1->row();
1486  typename MatrixTypeB::const_iterator b2 = B.begin();
1487  while (b2 != b_final)
1488  {
1489  const size_type j = b2->column();
1490  const size_type l = b2->row();
1491 
1492  const typename MatrixTypeA::value_type a = A.el(k, l);
1493 
1494  if (a != 0.)
1495  add(i, j, a * b1->value() * b2->value());
1496  ++b2;
1497  }
1498  ++b1;
1499  }
1500  else
1501  {
1502  // Determine minimal and
1503  // maximal row for a column in
1504  // advance.
1505 
1506  std::vector<size_type> minrow(B.n(), B.m());
1507  std::vector<size_type> maxrow(B.n(), 0);
1508  while (b1 != b_final)
1509  {
1510  const size_type r = b1->row();
1511  if (r < minrow[b1->column()])
1512  minrow[b1->column()] = r;
1513  if (r > maxrow[b1->column()])
1514  maxrow[b1->column()] = r;
1515  ++b1;
1516  }
1517 
1518  typename MatrixTypeA::const_iterator ai = A.begin();
1519  const typename MatrixTypeA::const_iterator ae = A.end();
1520 
1521  while (ai != ae)
1522  {
1523  const typename MatrixTypeA::value_type a = ai->value();
1524  // Don't do anything if
1525  // this entry is zero.
1526  if (a == 0.)
1527  continue;
1528 
1529  // Now, loop over all rows
1530  // having possibly a
1531  // nonzero entry in column
1532  // ai->row()
1533  b1 = B.begin(minrow[ai->row()]);
1534  const typename MatrixTypeB::const_iterator be1 =
1535  B.end(maxrow[ai->row()]);
1536  const typename MatrixTypeB::const_iterator be2 =
1537  B.end(maxrow[ai->column()]);
1538 
1539  while (b1 != be1)
1540  {
1541  const double b1v = b1->value();
1542  // We need the product
1543  // of both. If it is
1544  // zero, we can save
1545  // the work
1546  if (b1->column() == ai->row() && (b1v != 0.))
1547  {
1548  const size_type i = b1->row();
1549 
1550  typename MatrixTypeB::const_iterator b2 =
1551  B.begin(minrow[ai->column()]);
1552  while (b2 != be2)
1553  {
1554  if (b2->column() == ai->column())
1555  {
1556  const size_type j = b2->row();
1557  add(i, j, a * b1v * b2->value());
1558  }
1559  ++b2;
1560  }
1561  }
1562  ++b1;
1563  }
1564  ++ai;
1565  }
1566  }
1567 }
1568 
1569 
1570 template <typename number>
1571 template <class StreamType>
1572 inline void
1573 SparseMatrixEZ<number>::print_statistics(StreamType &out, bool full)
1574 {
1575  size_type used;
1576  size_type allocated;
1577  size_type reserved;
1578  std::vector<size_type> used_by_line;
1579 
1580  compute_statistics(used, allocated, reserved, used_by_line, full);
1581 
1582  out << "SparseMatrixEZ:used entries:" << used << std::endl
1583  << "SparseMatrixEZ:allocated entries:" << allocated << std::endl
1584  << "SparseMatrixEZ:reserved entries:" << reserved << std::endl;
1585 
1586  if (full)
1587  {
1588  for (size_type i = 0; i < used_by_line.size(); ++i)
1589  if (used_by_line[i] != 0)
1590  out << "SparseMatrixEZ:entries\t" << i << "\trows\t"
1591  << used_by_line[i] << std::endl;
1592  }
1593 }
1594 
1595 
1597 
1598 #endif
1599 /*---------------------------- sparse_matrix.h ---------------------------*/
size_type get_row_length(const size_type row) const
const types::global_dof_index invalid_size_type
Definition: types.h:205
void Tvmult_add(Vector< somenumber > &dst, const Vector< somenumber > &src) const
static ::ExceptionBase & ExcEntryAllocationFailure(int arg1, int arg2)
#define DeclException2(Exception2, type1, type2, outsequence)
Definition: exceptions.h:538
void add(const size_type i, const size_type j, const number value)
Contents is actually a matrix.
bool operator!=(const AlignedVector< T > &lhs, const AlignedVector< T > &rhs)
unsigned int increment
void compute_statistics(size_type &used, size_type &allocated, size_type &reserved, std::vector< size_type > &used_by_line, const bool compute_by_line) const
#define AssertIndexRange(index, range)
Definition: exceptions.h:1648
number operator()(const size_type i, const size_type j) const
SparseMatrixEZ< number > & copy_from(const MatrixType &source, const bool elide_zero_values=true)
void set(const size_type i, const size_type j, const number value, const bool elide_zero_values=true)
bool operator<(const SynchronousIterators< Iterators > &a, const SynchronousIterators< Iterators > &b)
void conjugate_add(const MatrixTypeA &A, const MatrixTypeB &B, const bool transpose=false)
size_type n() const
void block_read(std::istream &in)
bool operator<(const const_iterator &) const
void vmult_add(Vector< somenumber > &dst, const Vector< somenumber > &src) const
bool operator==(const AlignedVector< T > &lhs, const AlignedVector< T > &rhs)
void vmult(Vector< somenumber > &dst, const Vector< somenumber > &src) const
static const size_type invalid
Matrix is diagonal.
const Accessor & operator*() const
std::size_t memory_consumption() const
void threaded_matrix_norm_square(const Vector< somenumber > &v, const size_type begin_row, const size_type end_row, somenumber *partial_sum) const
void threaded_matrix_scalar_product(const Vector< somenumber > &u, const Vector< somenumber > &v, const size_type begin_row, const size_type end_row, somenumber *partial_sum) const
bool operator==(const const_iterator &) const
void print(std::ostream &out) const
static ::ExceptionBase & ExcNoDiagonal()
void precondition_SOR(Vector< somenumber > &dst, const Vector< somenumber > &src, const number om=1.) const
SparseMatrixEZ< number > & operator=(const SparseMatrixEZ< number > &)
Accessor(const SparseMatrixEZ< number > *matrix, const size_type row, const unsigned short index)
#define Assert(cond, exc)
Definition: exceptions.h:1423
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
void reinit(const size_type n_rows, const size_type n_columns, size_type default_row_length=0, unsigned int default_increment=1, size_type reserve=0)
bool operator!=(const const_iterator &) const
DerivativeForm< 1, spacedim, dim, Number > transpose(const DerivativeForm< 1, dim, spacedim, Number > &DF)
const_iterator(const SparseMatrixEZ< number > *matrix, const size_type row, const unsigned short index)
#define DeclException0(Exception0)
Definition: exceptions.h:470
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:369
const SparseMatrixEZ< number > * matrix
std::vector< Entry > data
SynchronousIterators< Iterators > operator++(SynchronousIterators< Iterators > &a)
Entry * allocate(const size_type row, const size_type col)
RowInfo(const size_type start=Entry::invalid)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
void print_statistics(StreamType &s, bool full=false)
static ::ExceptionBase & ExcIteratorPastEnd()
static const char A
unsigned int global_dof_index
Definition: types.h:76
void block_write(std::ostream &out) const
void swap(MemorySpaceData< Number, MemorySpace > &, MemorySpaceData< Number, MemorySpace > &)
const Accessor * operator->() const
void precondition_TSOR(Vector< somenumber > &dst, const Vector< somenumber > &src, const number om=1.) const
number el(const size_type i, const size_type j) const
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:368
const Entry * locate(const size_type row, const size_type col) const
const_iterator begin() const
static ::ExceptionBase & ExcInvalidEntry(int arg1, int arg2)
void precondition_Jacobi(Vector< somenumber > &dst, const Vector< somenumber > &src, const number omega=1.) const
void threaded_vmult(Vector< somenumber > &dst, const Vector< somenumber > &src, const size_type begin_row, const size_type end_row) const
~SparseMatrixEZ() override=default
void Tvmult(Vector< somenumber > &dst, const Vector< somenumber > &src) const
bool empty() const
std::enable_if< std::is_floating_point< T >::value &&std::is_floating_point< U >::value, typename ProductType< std::complex< T >, std::complex< U > >::type >::type operator*(const std::complex< T > &left, const std::complex< U > &right)
static const unsigned short invalid_diagonal
number l2_norm() const
void precondition_SSOR(Vector< somenumber > &dst, const Vector< somenumber > &src, const number om=1., const std::vector< std::size_t > &pos_right_of_diagonal=std::vector< std::size_t >()) const
const_iterator end() const
#define AssertIsFinite(number)
Definition: exceptions.h:1679
unsigned int saved_default_row_length
size_type m() const
std::vector< RowInfo > row_info
Tensor< 2, dim, Number > l(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
static ::ExceptionBase & ExcInternalError()
size_type n_nonzero_elements() 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