Reference documentation for deal.II version Git f0d1c24e5f 2021-10-18 08:09:39 -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 &
257  operator*() const;
258 
262  const Accessor *
263  operator->() const;
264 
268  bool
269  operator==(const const_iterator &) const;
273  bool
274  operator!=(const const_iterator &) const;
275 
280  bool
281  operator<(const const_iterator &) const;
282 
283  private:
288  };
289 
294  using value_type = number;
295 
303  SparseMatrixEZ();
304 
313 
320  explicit SparseMatrixEZ(const size_type n_rows,
321  const size_type n_columns,
322  const size_type default_row_length = 0,
323  const unsigned int default_increment = 1);
324 
328  ~SparseMatrixEZ() override = default;
329 
335 
345  operator=(const double d);
346 
354  void
355  reinit(const size_type n_rows,
356  const size_type n_columns,
357  size_type default_row_length = 0,
358  unsigned int default_increment = 1,
359  size_type reserve = 0);
360 
365  void
366  clear();
368 
376  bool
377  empty() const;
378 
383  size_type
384  m() const;
385 
390  size_type
391  n() const;
392 
396  size_type
397  get_row_length(const size_type row) const;
398 
402  size_type
403  n_nonzero_elements() const;
404 
409  std::size_t
410  memory_consumption() const;
411 
417  template <class StreamType>
418  void
419  print_statistics(StreamType &s, bool full = false);
420 
430  void
432  size_type & allocated,
433  size_type & reserved,
434  std::vector<size_type> &used_by_line,
435  const bool compute_by_line) const;
437 
458  void
459  set(const size_type i,
460  const size_type j,
461  const number value,
462  const bool elide_zero_values = true);
463 
474  void
475  add(const size_type i, const size_type j, const number value);
476 
491  template <typename number2>
492  void
493  add(const std::vector<size_type> &indices,
494  const FullMatrix<number2> & full_matrix,
495  const bool elide_zero_values = true);
496 
502  template <typename number2>
503  void
504  add(const std::vector<size_type> &row_indices,
505  const std::vector<size_type> &col_indices,
506  const FullMatrix<number2> & full_matrix,
507  const bool elide_zero_values = true);
508 
518  template <typename number2>
519  void
520  add(const size_type row,
521  const std::vector<size_type> &col_indices,
522  const std::vector<number2> & values,
523  const bool elide_zero_values = true);
524 
534  template <typename number2>
535  void
536  add(const size_type row,
537  const size_type n_cols,
538  const size_type *col_indices,
539  const number2 * values,
540  const bool elide_zero_values = true,
541  const bool col_indices_are_sorted = false);
542 
564  template <typename MatrixType>
566  copy_from(const MatrixType &source, const bool elide_zero_values = true);
567 
575  template <typename MatrixType>
576  void
577  add(const number factor, const MatrixType &matrix);
579 
592  number
593  operator()(const size_type i, const size_type j) const;
594 
599  number
600  el(const size_type i, const size_type j) const;
602 
610  template <typename somenumber>
611  void
612  vmult(Vector<somenumber> &dst, const Vector<somenumber> &src) const;
613 
619  template <typename somenumber>
620  void
621  Tvmult(Vector<somenumber> &dst, const Vector<somenumber> &src) const;
622 
627  template <typename somenumber>
628  void
629  vmult_add(Vector<somenumber> &dst, const Vector<somenumber> &src) const;
630 
636  template <typename somenumber>
637  void
638  Tvmult_add(Vector<somenumber> &dst, const Vector<somenumber> &src) const;
640 
647  number
648  l2_norm() const;
650 
659  template <typename somenumber>
660  void
661  precondition_Jacobi(Vector<somenumber> & dst,
662  const Vector<somenumber> &src,
663  const number omega = 1.) const;
664 
668  template <typename somenumber>
669  void
670  precondition_SSOR(Vector<somenumber> & dst,
671  const Vector<somenumber> & src,
672  const number om = 1.,
673  const std::vector<std::size_t> &pos_right_of_diagonal =
674  std::vector<std::size_t>()) const;
675 
680  template <typename somenumber>
681  void
682  precondition_SOR(Vector<somenumber> & dst,
683  const Vector<somenumber> &src,
684  const number om = 1.) const;
685 
690  template <typename somenumber>
691  void
692  precondition_TSOR(Vector<somenumber> & dst,
693  const Vector<somenumber> &src,
694  const number om = 1.) const;
695 
704  template <typename MatrixTypeA, typename MatrixTypeB>
705  void
706  conjugate_add(const MatrixTypeA &A,
707  const MatrixTypeB &B,
708  const bool transpose = false);
710 
718  begin() const;
719 
724  end() const;
725 
731  begin(const size_type r) const;
732 
738  end(const size_type r) const;
740 
748  void
749  print(std::ostream &out) const;
750 
771  void
772  print_formatted(std::ostream & out,
773  const unsigned int precision = 3,
774  const bool scientific = true,
775  const unsigned int width = 0,
776  const char * zero_string = " ",
777  const double denominator = 1.) const;
778 
784  void
785  block_write(std::ostream &out) const;
786 
797  void
798  block_read(std::istream &in);
800 
810 
815  int,
816  int,
817  << "The entry with index (" << arg1 << ',' << arg2
818  << ") does not exist.");
819 
821  int,
822  int,
823  << "An entry with index (" << arg1 << ',' << arg2
824  << ") cannot be allocated.");
826 private:
831  const Entry *
832  locate(const size_type row, const size_type col) const;
833 
838  Entry *
839  locate(const size_type row, const size_type col);
840 
844  Entry *
845  allocate(const size_type row, const size_type col);
846 
852  template <typename somenumber>
853  void
854  threaded_vmult(Vector<somenumber> & dst,
855  const Vector<somenumber> &src,
856  const size_type begin_row,
857  const size_type end_row) const;
858 
864  template <typename somenumber>
865  void
866  threaded_matrix_norm_square(const Vector<somenumber> &v,
867  const size_type begin_row,
868  const size_type end_row,
869  somenumber * partial_sum) const;
870 
876  template <typename somenumber>
877  void
878  threaded_matrix_scalar_product(const Vector<somenumber> &u,
879  const Vector<somenumber> &v,
880  const size_type begin_row,
881  const size_type end_row,
882  somenumber * partial_sum) const;
883 
888 
892  std::vector<RowInfo> row_info;
893 
897  std::vector<Entry> data;
898 
902  unsigned int increment;
903 
908 };
909 
913 /*---------------------- Inline functions -----------------------------------*/
914 
915 template <typename number>
917  const number & value)
918  : column(column)
919  , value(value)
920 {}
921 
922 
923 
924 template <typename number>
926  : column(invalid)
927  , value(0)
928 {}
929 
930 
931 template <typename number>
933  : start(start)
934  , length(0)
935  , diagonal(invalid_diagonal)
936 {}
937 
938 
939 //---------------------------------------------------------------------------
940 template <typename number>
943  const size_type r,
944  const unsigned short i)
945  : matrix(matrix)
946  , a_row(r)
947  , a_index(i)
948 {}
949 
950 
951 template <typename number>
952 inline typename SparseMatrixEZ<number>::size_type
954 {
955  return a_row;
956 }
957 
958 
959 template <typename number>
960 inline typename SparseMatrixEZ<number>::size_type
962 {
963  return matrix->data[matrix->row_info[a_row].start + a_index].column;
964 }
965 
966 
967 template <typename number>
968 inline unsigned short
970 {
971  return a_index;
972 }
973 
974 
975 
976 template <typename number>
977 inline number
979 {
980  return matrix->data[matrix->row_info[a_row].start + a_index].value;
981 }
982 
983 
984 template <typename number>
987  const size_type r,
988  const unsigned short i)
989  : accessor(matrix, r, i)
990 {
991  // Finish if this is the end()
992  if (r == accessor.matrix->m() && i == 0)
993  return;
994 
995  // Make sure we never construct an
996  // iterator pointing to a
997  // non-existing entry
998 
999  // If the index points beyond the
1000  // end of the row, try the next
1001  // row.
1002  if (accessor.a_index >= accessor.matrix->row_info[accessor.a_row].length)
1003  {
1004  do
1005  {
1006  ++accessor.a_row;
1007  }
1008  // Beware! If the next row is
1009  // empty, iterate until a
1010  // non-empty row is found or we
1011  // hit the end of the matrix.
1012  while (accessor.a_row < accessor.matrix->m() &&
1013  accessor.matrix->row_info[accessor.a_row].length == 0);
1014  }
1015 }
1016 
1017 
1018 template <typename number>
1021 {
1023 
1024  // Increment column index
1025  ++(accessor.a_index);
1026  // If index exceeds number of
1027  // entries in this row, proceed
1028  // with next row.
1029  if (accessor.a_index >= accessor.matrix->row_info[accessor.a_row].length)
1030  {
1031  accessor.a_index = 0;
1032  // Do this loop to avoid
1033  // elements in empty rows
1034  do
1035  {
1036  ++accessor.a_row;
1037  }
1038  while (accessor.a_row < accessor.matrix->m() &&
1039  accessor.matrix->row_info[accessor.a_row].length == 0);
1040  }
1041  return *this;
1042 }
1043 
1044 
1045 template <typename number>
1048 {
1049  return accessor;
1050 }
1051 
1052 
1053 template <typename number>
1056 {
1057  return &accessor;
1058 }
1059 
1060 
1061 template <typename number>
1062 inline bool
1064  const const_iterator &other) const
1065 {
1066  return (accessor.row() == other.accessor.row() &&
1067  accessor.index() == other.accessor.index());
1068 }
1069 
1070 
1071 template <typename number>
1072 inline bool
1074  const const_iterator &other) const
1075 {
1076  return !(*this == other);
1077 }
1078 
1079 
1080 template <typename number>
1081 inline bool
1083  const const_iterator &other) const
1084 {
1085  return (accessor.row() < other.accessor.row() ||
1086  (accessor.row() == other.accessor.row() &&
1087  accessor.index() < other.accessor.index()));
1088 }
1089 
1090 
1091 //---------------------------------------------------------------------------
1092 template <typename number>
1093 inline typename SparseMatrixEZ<number>::size_type
1095 {
1096  return row_info.size();
1097 }
1098 
1099 
1100 template <typename number>
1101 inline typename SparseMatrixEZ<number>::size_type
1103 {
1104  return n_columns;
1105 }
1106 
1107 
1108 template <typename number>
1109 inline typename SparseMatrixEZ<number>::Entry *
1111 {
1112  AssertIndexRange(row, m());
1113  AssertIndexRange(col, n());
1114 
1115  const RowInfo & r = row_info[row];
1116  const size_type end = r.start + r.length;
1117  for (size_type i = r.start; i < end; ++i)
1118  {
1119  Entry *const entry = &data[i];
1120  if (entry->column == col)
1121  return entry;
1122  if (entry->column == Entry::invalid)
1123  return nullptr;
1124  }
1125  return nullptr;
1126 }
1127 
1128 
1129 
1130 template <typename number>
1131 inline const typename SparseMatrixEZ<number>::Entry *
1133 {
1134  SparseMatrixEZ<number> *t = const_cast<SparseMatrixEZ<number> *>(this);
1135  return t->locate(row, col);
1136 }
1137 
1138 
1139 template <typename number>
1140 inline typename SparseMatrixEZ<number>::Entry *
1142 {
1143  AssertIndexRange(row, m());
1144  AssertIndexRange(col, n());
1145 
1146  RowInfo & r = row_info[row];
1147  const size_type end = r.start + r.length;
1148 
1149  size_type i = r.start;
1150  // If diagonal exists and this
1151  // column is higher, start only
1152  // after diagonal.
1153  if (r.diagonal != RowInfo::invalid_diagonal && col >= row)
1154  i += r.diagonal;
1155  // Find position of entry
1156  while (i < end && data[i].column < col)
1157  ++i;
1158 
1159  // entry found
1160  if (i != end && data[i].column == col)
1161  return &data[i];
1162 
1163  // Now, we must insert the new
1164  // entry and move all successive
1165  // entries back.
1166 
1167  // If no more space is available
1168  // for this row, insert new
1169  // elements into the vector.
1170  // TODO:[GK] We should not extend this row if i<end
1171  if (row != row_info.size() - 1)
1172  {
1173  if (end >= row_info[row + 1].start)
1174  {
1175  // Failure if increment 0
1176  Assert(increment != 0, ExcEntryAllocationFailure(row, col));
1177 
1178  // Insert new entries
1179  data.insert(data.begin() + end, increment, Entry());
1180  // Update starts of
1181  // following rows
1182  for (size_type rn = row + 1; rn < row_info.size(); ++rn)
1183  row_info[rn].start += increment;
1184  }
1185  }
1186  else
1187  {
1188  if (end >= data.size())
1189  {
1190  // Here, appending a block
1191  // does not increase
1192  // performance.
1193  data.push_back(Entry());
1194  }
1195  }
1196 
1197  Entry *entry = &data[i];
1198  // Save original entry
1199  Entry temp = *entry;
1200  // Insert new entry here to
1201  // make sure all entries
1202  // are ordered by column
1203  // index
1204  entry->column = col;
1205  entry->value = 0;
1206  // Update row_info
1207  ++r.length;
1208  if (col == row)
1209  r.diagonal = i - r.start;
1210  else if (col < row && r.diagonal != RowInfo::invalid_diagonal)
1211  ++r.diagonal;
1212 
1213  if (i == end)
1214  return entry;
1215 
1216  // Move all entries in this
1217  // row up by one
1218  for (size_type j = i + 1; j < end; ++j)
1219  {
1220  // There should be no invalid
1221  // entry below end
1222  Assert(data[j].column != Entry::invalid, ExcInternalError());
1223 
1224  // TODO[GK]: This could be done more efficiently by moving starting at the
1225  // top rather than swapping starting at the bottom
1226  std::swap(data[j], temp);
1227  }
1228  Assert(data[end].column == Entry::invalid, ExcInternalError());
1229 
1230  data[end] = temp;
1231 
1232  return entry;
1233 }
1234 
1235 
1236 
1237 template <typename number>
1238 inline void
1240  const size_type j,
1241  const number value,
1242  const bool elide_zero_values)
1243 {
1244  AssertIsFinite(value);
1245 
1246  AssertIndexRange(i, m());
1247  AssertIndexRange(j, n());
1248 
1249  if (elide_zero_values && value == 0.)
1250  {
1251  Entry *entry = locate(i, j);
1252  if (entry != nullptr)
1253  entry->value = 0.;
1254  }
1255  else
1256  {
1257  Entry *entry = allocate(i, j);
1258  entry->value = value;
1259  }
1260 }
1261 
1262 
1263 
1264 template <typename number>
1265 inline void
1267  const size_type j,
1268  const number value)
1269 {
1270  AssertIsFinite(value);
1271 
1272  AssertIndexRange(i, m());
1273  AssertIndexRange(j, n());
1274 
1275  // ignore zero additions
1276  if (std::abs(value) == 0.)
1277  return;
1278 
1279  Entry *entry = allocate(i, j);
1280  entry->value += value;
1281 }
1282 
1283 
1284 template <typename number>
1285 template <typename number2>
1286 void
1287 SparseMatrixEZ<number>::add(const std::vector<size_type> &indices,
1288  const FullMatrix<number2> & full_matrix,
1289  const bool elide_zero_values)
1290 {
1291  // TODO: This function can surely be made more efficient
1292  for (size_type i = 0; i < indices.size(); ++i)
1293  for (size_type j = 0; j < indices.size(); ++j)
1294  if ((full_matrix(i, j) != 0) || (elide_zero_values == false))
1295  add(indices[i], indices[j], full_matrix(i, j));
1296 }
1297 
1298 
1299 
1300 template <typename number>
1301 template <typename number2>
1302 void
1303 SparseMatrixEZ<number>::add(const std::vector<size_type> &row_indices,
1304  const std::vector<size_type> &col_indices,
1305  const FullMatrix<number2> & full_matrix,
1306  const bool elide_zero_values)
1307 {
1308  // TODO: This function can surely be made more efficient
1309  for (size_type i = 0; i < row_indices.size(); ++i)
1310  for (size_type j = 0; j < col_indices.size(); ++j)
1311  if ((full_matrix(i, j) != 0) || (elide_zero_values == false))
1312  add(row_indices[i], col_indices[j], full_matrix(i, j));
1313 }
1314 
1315 
1316 
1317 template <typename number>
1318 template <typename number2>
1319 void
1321  const std::vector<size_type> &col_indices,
1322  const std::vector<number2> & values,
1323  const bool elide_zero_values)
1324 {
1325  // TODO: This function can surely be made more efficient
1326  for (size_type j = 0; j < col_indices.size(); ++j)
1327  if ((values[j] != 0) || (elide_zero_values == false))
1328  add(row, col_indices[j], values[j]);
1329 }
1330 
1331 
1332 
1333 template <typename number>
1334 template <typename number2>
1335 void
1337  const size_type n_cols,
1338  const size_type *col_indices,
1339  const number2 * values,
1340  const bool elide_zero_values,
1341  const bool /*col_indices_are_sorted*/)
1342 {
1343  // TODO: This function can surely be made more efficient
1344  for (size_type j = 0; j < n_cols; ++j)
1345  if ((std::abs(values[j]) != 0) || (elide_zero_values == false))
1346  add(row, col_indices[j], values[j]);
1347 }
1348 
1349 
1350 
1351 template <typename number>
1352 inline number
1354 {
1355  const Entry *entry = locate(i, j);
1356  if (entry)
1357  return entry->value;
1358  return 0.;
1359 }
1360 
1361 
1362 
1363 template <typename number>
1364 inline number
1366 {
1367  const Entry *entry = locate(i, j);
1368  if (entry)
1369  return entry->value;
1370  Assert(false, ExcInvalidEntry(i, j));
1371  return 0.;
1372 }
1373 
1374 
1375 template <typename number>
1378 {
1379  const_iterator result(this, 0, 0);
1380  return result;
1381 }
1382 
1383 template <typename number>
1386 {
1387  return const_iterator(this, m(), 0);
1388 }
1389 
1390 template <typename number>
1393 {
1394  AssertIndexRange(r, m());
1395  const_iterator result(this, r, 0);
1396  return result;
1397 }
1398 
1399 template <typename number>
1402 {
1403  AssertIndexRange(r, m());
1404  const_iterator result(this, r + 1, 0);
1405  return result;
1406 }
1407 
1408 template <typename number>
1409 template <typename MatrixType>
1410 inline SparseMatrixEZ<number> &
1412  const bool elide_zero_values)
1413 {
1414  reinit(M.m(), M.n(), this->saved_default_row_length, this->increment);
1415 
1416  // loop over the elements of the argument matrix row by row, as suggested
1417  // in the documentation of the sparse matrix iterator class, and
1418  // copy them into the current object
1419  for (size_type row = 0; row < M.m(); ++row)
1420  {
1421  const typename MatrixType::const_iterator end_row = M.end(row);
1422  for (typename MatrixType::const_iterator entry = M.begin(row);
1423  entry != end_row;
1424  ++entry)
1425  set(row, entry->column(), entry->value(), elide_zero_values);
1426  }
1427 
1428  return *this;
1429 }
1430 
1431 template <typename number>
1432 template <typename MatrixType>
1433 inline void
1434 SparseMatrixEZ<number>::add(const number factor, const MatrixType &M)
1435 {
1436  Assert(M.m() == m(), ExcDimensionMismatch(M.m(), m()));
1437  Assert(M.n() == n(), ExcDimensionMismatch(M.n(), n()));
1438 
1439  if (factor == 0.)
1440  return;
1441 
1442  // loop over the elements of the argument matrix row by row, as suggested
1443  // in the documentation of the sparse matrix iterator class, and
1444  // add them into the current object
1445  for (size_type row = 0; row < M.m(); ++row)
1446  {
1447  const typename MatrixType::const_iterator end_row = M.end(row);
1448  for (typename MatrixType::const_iterator entry = M.begin(row);
1449  entry != end_row;
1450  ++entry)
1451  if (entry->value() != 0)
1452  add(row, entry->column(), factor * entry->value());
1453  }
1454 }
1455 
1456 
1457 
1458 template <typename number>
1459 template <typename MatrixTypeA, typename MatrixTypeB>
1460 inline void
1462  const MatrixTypeB &B,
1463  const bool transpose)
1464 {
1465  // Compute the result
1466  // r_ij = \sum_kl b_ik b_jl a_kl
1467 
1468  // Assert (n() == B.m(), ExcDimensionMismatch(n(), B.m()));
1469  // Assert (m() == B.m(), ExcDimensionMismatch(m(), B.m()));
1470  // Assert (A.n() == B.n(), ExcDimensionMismatch(A.n(), B.n()));
1471  // Assert (A.m() == B.n(), ExcDimensionMismatch(A.m(), B.n()));
1472 
1473  // Somehow, we have to avoid making
1474  // this an operation of complexity
1475  // n^2. For the transpose case, we
1476  // can go through the non-zero
1477  // elements of A^-1 and use the
1478  // corresponding rows of B only.
1479  // For the non-transpose case, we
1480  // must find a trick.
1481  typename MatrixTypeB::const_iterator b1 = B.begin();
1482  const typename MatrixTypeB::const_iterator b_final = B.end();
1483  if (transpose)
1484  while (b1 != b_final)
1485  {
1486  const size_type i = b1->column();
1487  const size_type k = b1->row();
1488  typename MatrixTypeB::const_iterator b2 = B.begin();
1489  while (b2 != b_final)
1490  {
1491  const size_type j = b2->column();
1492  const size_type l = b2->row();
1493 
1494  const typename MatrixTypeA::value_type a = A.el(k, l);
1495 
1496  if (a != 0.)
1497  add(i, j, a * b1->value() * b2->value());
1498  ++b2;
1499  }
1500  ++b1;
1501  }
1502  else
1503  {
1504  // Determine minimal and
1505  // maximal row for a column in
1506  // advance.
1507 
1508  std::vector<size_type> minrow(B.n(), B.m());
1509  std::vector<size_type> maxrow(B.n(), 0);
1510  while (b1 != b_final)
1511  {
1512  const size_type r = b1->row();
1513  if (r < minrow[b1->column()])
1514  minrow[b1->column()] = r;
1515  if (r > maxrow[b1->column()])
1516  maxrow[b1->column()] = r;
1517  ++b1;
1518  }
1519 
1520  typename MatrixTypeA::const_iterator ai = A.begin();
1521  const typename MatrixTypeA::const_iterator ae = A.end();
1522 
1523  while (ai != ae)
1524  {
1525  const typename MatrixTypeA::value_type a = ai->value();
1526  // Don't do anything if
1527  // this entry is zero.
1528  if (a == 0.)
1529  continue;
1530 
1531  // Now, loop over all rows
1532  // having possibly a
1533  // nonzero entry in column
1534  // ai->row()
1535  b1 = B.begin(minrow[ai->row()]);
1536  const typename MatrixTypeB::const_iterator be1 =
1537  B.end(maxrow[ai->row()]);
1538  const typename MatrixTypeB::const_iterator be2 =
1539  B.end(maxrow[ai->column()]);
1540 
1541  while (b1 != be1)
1542  {
1543  const double b1v = b1->value();
1544  // We need the product
1545  // of both. If it is
1546  // zero, we can save
1547  // the work
1548  if (b1->column() == ai->row() && (b1v != 0.))
1549  {
1550  const size_type i = b1->row();
1551 
1552  typename MatrixTypeB::const_iterator b2 =
1553  B.begin(minrow[ai->column()]);
1554  while (b2 != be2)
1555  {
1556  if (b2->column() == ai->column())
1557  {
1558  const size_type j = b2->row();
1559  add(i, j, a * b1v * b2->value());
1560  }
1561  ++b2;
1562  }
1563  }
1564  ++b1;
1565  }
1566  ++ai;
1567  }
1568  }
1569 }
1570 
1571 
1572 template <typename number>
1573 template <class StreamType>
1574 inline void
1575 SparseMatrixEZ<number>::print_statistics(StreamType &out, bool full)
1576 {
1577  size_type used;
1578  size_type allocated;
1579  size_type reserved;
1580  std::vector<size_type> used_by_line;
1581 
1582  compute_statistics(used, allocated, reserved, used_by_line, full);
1583 
1584  out << "SparseMatrixEZ:used entries:" << used << std::endl
1585  << "SparseMatrixEZ:allocated entries:" << allocated << std::endl
1586  << "SparseMatrixEZ:reserved entries:" << reserved << std::endl;
1587 
1588  if (full)
1589  {
1590  for (size_type i = 0; i < used_by_line.size(); ++i)
1591  if (used_by_line[i] != 0)
1592  out << "SparseMatrixEZ:entries\t" << i << "\trows\t"
1593  << used_by_line[i] << std::endl;
1594  }
1595 }
1596 
1597 
1599 
1600 #endif
1601 /*---------------------------- 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:532
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:1720
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:1461
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:464
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:401
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:400
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:1746
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