Reference documentation for deal.II version 9.1.1
\(\newcommand{\dealcoloneq}{\mathrel{\vcenter{:}}=}\)
sparse_matrix_ez.h
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2002 - 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_sparse_matrix_ez_h
17 # define dealii_sparse_matrix_ez_h
18 
19 
20 # include <deal.II/base/config.h>
21 
22 # include <deal.II/base/smartpointer.h>
23 # include <deal.II/base/subscriptor.h>
24 
25 # include <deal.II/lac/exceptions.h>
26 
27 # include <vector>
28 
29 DEAL_II_NAMESPACE_OPEN
30 
31 template <typename number>
32 class Vector;
33 template <typename number>
34 class FullMatrix;
35 
103 template <typename number>
104 class SparseMatrixEZ : public Subscriptor
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  {
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 
238  friend class const_iterator;
239  };
240 
241  public:
246  const size_type row,
247  const unsigned short index);
248 
253  operator++();
254 
258  const Accessor &operator*() const;
259 
263  const Accessor *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 
296  };
297 
302  using value_type = number;
303 
311  SparseMatrixEZ();
312 
321 
328  explicit SparseMatrixEZ(const size_type n_rows,
329  const size_type n_columns,
330  const size_type default_row_length = 0,
331  const unsigned int default_increment = 1);
332 
336  ~SparseMatrixEZ() override = default;
337 
343 
353  operator=(const double d);
354 
362  void
363  reinit(const size_type n_rows,
364  const size_type n_columns,
365  size_type default_row_length = 0,
366  unsigned int default_increment = 1,
367  size_type reserve = 0);
368 
373  void
374  clear();
376 
384  bool
385  empty() const;
386 
391  size_type
392  m() const;
393 
398  size_type
399  n() const;
400 
404  size_type
405  get_row_length(const size_type row) const;
406 
410  size_type
411  n_nonzero_elements() const;
412 
417  std::size_t
418  memory_consumption() const;
419 
425  template <class StreamType>
426  void
427  print_statistics(StreamType &s, bool full = false);
428 
438  void
440  size_type & allocated,
441  size_type & reserved,
442  std::vector<size_type> &used_by_line,
443  const bool compute_by_line) const;
445 
465  void
466  set(const size_type i,
467  const size_type j,
468  const number value,
469  const bool elide_zero_values = true);
470 
476  void
477  add(const size_type i, const size_type j, const number value);
478 
493  template <typename number2>
494  void
495  add(const std::vector<size_type> &indices,
496  const FullMatrix<number2> & full_matrix,
497  const bool elide_zero_values = true);
498 
504  template <typename number2>
505  void
506  add(const std::vector<size_type> &row_indices,
507  const std::vector<size_type> &col_indices,
508  const FullMatrix<number2> & full_matrix,
509  const bool elide_zero_values = true);
510 
520  template <typename number2>
521  void
522  add(const size_type row,
523  const std::vector<size_type> &col_indices,
524  const std::vector<number2> & values,
525  const bool elide_zero_values = true);
526 
536  template <typename number2>
537  void
538  add(const size_type row,
539  const size_type n_cols,
540  const size_type *col_indices,
541  const number2 * values,
542  const bool elide_zero_values = true,
543  const bool col_indices_are_sorted = false);
544 
566  template <typename MatrixType>
568  copy_from(const MatrixType &source, const bool elide_zero_values = true);
569 
577  template <typename MatrixType>
578  void
579  add(const number factor, const MatrixType &matrix);
581 
594  number
595  operator()(const size_type i, const size_type j) const;
596 
601  number
602  el(const size_type i, const size_type j) const;
604 
612  template <typename somenumber>
613  void
614  vmult(Vector<somenumber> &dst, const Vector<somenumber> &src) const;
615 
621  template <typename somenumber>
622  void
623  Tvmult(Vector<somenumber> &dst, const Vector<somenumber> &src) const;
624 
629  template <typename somenumber>
630  void
631  vmult_add(Vector<somenumber> &dst, const Vector<somenumber> &src) const;
632 
638  template <typename somenumber>
639  void
640  Tvmult_add(Vector<somenumber> &dst, const Vector<somenumber> &src) const;
642 
649  number
650  l2_norm() const;
652 
661  template <typename somenumber>
662  void
663  precondition_Jacobi(Vector<somenumber> & dst,
664  const Vector<somenumber> &src,
665  const number omega = 1.) const;
666 
670  template <typename somenumber>
671  void
672  precondition_SSOR(Vector<somenumber> & dst,
673  const Vector<somenumber> & src,
674  const number om = 1.,
675  const std::vector<std::size_t> &pos_right_of_diagonal =
676  std::vector<std::size_t>()) const;
677 
682  template <typename somenumber>
683  void
684  precondition_SOR(Vector<somenumber> & dst,
685  const Vector<somenumber> &src,
686  const number om = 1.) const;
687 
692  template <typename somenumber>
693  void
694  precondition_TSOR(Vector<somenumber> & dst,
695  const Vector<somenumber> &src,
696  const number om = 1.) const;
697 
706  template <typename MatrixTypeA, typename MatrixTypeB>
707  void
708  conjugate_add(const MatrixTypeA &A,
709  const MatrixTypeB &B,
710  const bool transpose = false);
712 
720  begin() const;
721 
726  end() const;
727 
733  begin(const size_type r) const;
734 
740  end(const size_type r) const;
742 
750  void
751  print(std::ostream &out) const;
752 
773  void
774  print_formatted(std::ostream & out,
775  const unsigned int precision = 3,
776  const bool scientific = true,
777  const unsigned int width = 0,
778  const char * zero_string = " ",
779  const double denominator = 1.) const;
780 
786  void
787  block_write(std::ostream &out) const;
788 
799  void
800  block_read(std::istream &in);
802 
812 
817  int,
818  int,
819  << "The entry with index (" << arg1 << ',' << arg2
820  << ") does not exist.");
821 
823  int,
824  int,
825  << "An entry with index (" << arg1 << ',' << arg2
826  << ") cannot be allocated.");
828 private:
833  const Entry *
834  locate(const size_type row, const size_type col) const;
835 
840  Entry *
841  locate(const size_type row, const size_type col);
842 
846  Entry *
847  allocate(const size_type row, const size_type col);
848 
854  template <typename somenumber>
855  void
856  threaded_vmult(Vector<somenumber> & dst,
857  const Vector<somenumber> &src,
858  const size_type begin_row,
859  const size_type end_row) const;
860 
866  template <typename somenumber>
867  void
868  threaded_matrix_norm_square(const Vector<somenumber> &v,
869  const size_type begin_row,
870  const size_type end_row,
871  somenumber * partial_sum) const;
872 
878  template <typename somenumber>
879  void
880  threaded_matrix_scalar_product(const Vector<somenumber> &u,
881  const Vector<somenumber> &v,
882  const size_type begin_row,
883  const size_type end_row,
884  somenumber * partial_sum) const;
885 
890 
894  std::vector<RowInfo> row_info;
895 
899  std::vector<Entry> data;
900 
904  unsigned int increment;
905 
910 };
911 
915 /*---------------------- Inline functions -----------------------------------*/
916 
917 template <typename number>
919  const number & value)
920  : column(column)
921  , value(value)
922 {}
923 
924 
925 
926 template <typename number>
928  : column(invalid)
929  , value(0)
930 {}
931 
932 
933 template <typename number>
935  : start(start)
936  , length(0)
937  , diagonal(invalid_diagonal)
938 {}
939 
940 
941 //---------------------------------------------------------------------------
942 template <typename number>
944  const SparseMatrixEZ<number> *matrix,
945  const size_type r,
946  const unsigned short i)
947  : matrix(matrix)
948  , a_row(r)
949  , a_index(i)
950 {}
951 
952 
953 template <typename number>
954 inline typename SparseMatrixEZ<number>::size_type
956 {
957  return a_row;
958 }
959 
960 
961 template <typename number>
962 inline typename SparseMatrixEZ<number>::size_type
964 {
965  return matrix->data[matrix->row_info[a_row].start + a_index].column;
966 }
967 
968 
969 template <typename number>
970 inline unsigned short
972 {
973  return a_index;
974 }
975 
976 
977 
978 template <typename number>
979 inline number
981 {
982  return matrix->data[matrix->row_info[a_row].start + a_index].value;
983 }
984 
985 
986 template <typename number>
988  const SparseMatrixEZ<number> *matrix,
989  const size_type r,
990  const unsigned short i)
991  : accessor(matrix, r, i)
992 {
993  // Finish if this is the end()
994  if (r == accessor.matrix->m() && i == 0)
995  return;
996 
997  // Make sure we never construct an
998  // iterator pointing to a
999  // non-existing entry
1000 
1001  // If the index points beyond the
1002  // end of the row, try the next
1003  // row.
1004  if (accessor.a_index >= accessor.matrix->row_info[accessor.a_row].length)
1005  {
1006  do
1007  {
1008  ++accessor.a_row;
1009  }
1010  // Beware! If the next row is
1011  // empty, iterate until a
1012  // non-empty row is found or we
1013  // hit the end of the matrix.
1014  while (accessor.a_row < accessor.matrix->m() &&
1015  accessor.matrix->row_info[accessor.a_row].length == 0);
1016  }
1017 }
1018 
1019 
1020 template <typename number>
1023 {
1024  Assert(accessor.a_row < accessor.matrix->m(), ExcIteratorPastEnd());
1025 
1026  // Increment column index
1027  ++(accessor.a_index);
1028  // If index exceeds number of
1029  // entries in this row, proceed
1030  // with next row.
1031  if (accessor.a_index >= accessor.matrix->row_info[accessor.a_row].length)
1032  {
1033  accessor.a_index = 0;
1034  // Do this loop to avoid
1035  // elements in empty rows
1036  do
1037  {
1038  ++accessor.a_row;
1039  }
1040  while (accessor.a_row < accessor.matrix->m() &&
1041  accessor.matrix->row_info[accessor.a_row].length == 0);
1042  }
1043  return *this;
1044 }
1045 
1046 
1047 template <typename number>
1050 {
1051  return accessor;
1052 }
1053 
1054 
1055 template <typename number>
1058 {
1059  return &accessor;
1060 }
1061 
1062 
1063 template <typename number>
1064 inline bool
1066 operator==(const const_iterator &other) const
1067 {
1068  return (accessor.row() == other.accessor.row() &&
1069  accessor.index() == other.accessor.index());
1070 }
1071 
1072 
1073 template <typename number>
1074 inline bool
1076 operator!=(const const_iterator &other) const
1077 {
1078  return !(*this == other);
1079 }
1080 
1081 
1082 template <typename number>
1083 inline bool
1085 operator<(const const_iterator &other) const
1086 {
1087  return (accessor.row() < other.accessor.row() ||
1088  (accessor.row() == other.accessor.row() &&
1089  accessor.index() < other.accessor.index()));
1090 }
1091 
1092 
1093 //---------------------------------------------------------------------------
1094 template <typename number>
1095 inline typename SparseMatrixEZ<number>::size_type
1097 {
1098  return row_info.size();
1099 }
1100 
1101 
1102 template <typename number>
1103 inline typename SparseMatrixEZ<number>::size_type
1105 {
1106  return n_columns;
1107 }
1108 
1109 
1110 template <typename number>
1111 inline typename SparseMatrixEZ<number>::Entry *
1113 {
1114  Assert(row < m(), ExcIndexRange(row, 0, m()));
1115  Assert(col < n(), ExcIndexRange(col, 0, n()));
1116 
1117  const RowInfo & r = row_info[row];
1118  const size_type end = r.start + r.length;
1119  for (size_type i = r.start; i < end; ++i)
1120  {
1121  Entry *const entry = &data[i];
1122  if (entry->column == col)
1123  return entry;
1124  if (entry->column == Entry::invalid)
1125  return nullptr;
1126  }
1127  return nullptr;
1128 }
1129 
1130 
1131 
1132 template <typename number>
1133 inline const typename SparseMatrixEZ<number>::Entry *
1135 {
1136  SparseMatrixEZ<number> *t = const_cast<SparseMatrixEZ<number> *>(this);
1137  return t->locate(row, col);
1138 }
1139 
1140 
1141 template <typename number>
1142 inline typename SparseMatrixEZ<number>::Entry *
1144 {
1145  Assert(row < m(), ExcIndexRange(row, 0, m()));
1146  Assert(col < n(), ExcIndexRange(col, 0, n()));
1147 
1148  RowInfo & r = row_info[row];
1149  const size_type end = r.start + r.length;
1150 
1151  size_type i = r.start;
1152  // If diagonal exists and this
1153  // column is higher, start only
1154  // after diagonal.
1155  if (r.diagonal != RowInfo::invalid_diagonal && col >= row)
1156  i += r.diagonal;
1157  // Find position of entry
1158  while (i < end && data[i].column < col)
1159  ++i;
1160 
1161  // entry found
1162  if (i != end && data[i].column == col)
1163  return &data[i];
1164 
1165  // Now, we must insert the new
1166  // entry and move all successive
1167  // entries back.
1168 
1169  // If no more space is available
1170  // for this row, insert new
1171  // elements into the vector.
1172  // TODO:[GK] We should not extend this row if i<end
1173  if (row != row_info.size() - 1)
1174  {
1175  if (end >= row_info[row + 1].start)
1176  {
1177  // Failure if increment 0
1178  Assert(increment != 0, ExcEntryAllocationFailure(row, col));
1179 
1180  // Insert new entries
1181  data.insert(data.begin() + end, increment, Entry());
1182  // Update starts of
1183  // following rows
1184  for (size_type rn = row + 1; rn < row_info.size(); ++rn)
1185  row_info[rn].start += increment;
1186  }
1187  }
1188  else
1189  {
1190  if (end >= data.size())
1191  {
1192  // Here, appending a block
1193  // does not increase
1194  // performance.
1195  data.push_back(Entry());
1196  }
1197  }
1198 
1199  Entry *entry = &data[i];
1200  // Save original entry
1201  Entry temp = *entry;
1202  // Insert new entry here to
1203  // make sure all entries
1204  // are ordered by column
1205  // index
1206  entry->column = col;
1207  entry->value = 0;
1208  // Update row_info
1209  ++r.length;
1210  if (col == row)
1211  r.diagonal = i - r.start;
1212  else if (col < row && r.diagonal != RowInfo::invalid_diagonal)
1213  ++r.diagonal;
1214 
1215  if (i == end)
1216  return entry;
1217 
1218  // Move all entries in this
1219  // row up by one
1220  for (size_type j = i + 1; j < end; ++j)
1221  {
1222  // There should be no invalid
1223  // entry below end
1224  Assert(data[j].column != Entry::invalid, ExcInternalError());
1225 
1226  // TODO[GK]: This could be done more efficiently by moving starting at the
1227  // top rather than swapping starting at the bottom
1228  std::swap(data[j], temp);
1229  }
1231 
1232  data[end] = temp;
1233 
1234  return entry;
1235 }
1236 
1237 
1238 
1239 template <typename number>
1240 inline void
1242  const size_type j,
1243  const number value,
1244  const bool elide_zero_values)
1245 {
1246  AssertIsFinite(value);
1247 
1248  Assert(i < m(), ExcIndexRange(i, 0, m()));
1249  Assert(j < n(), ExcIndexRange(j, 0, n()));
1250 
1251  if (elide_zero_values && value == 0.)
1252  {
1253  Entry *entry = locate(i, j);
1254  if (entry != nullptr)
1255  entry->value = 0.;
1256  }
1257  else
1258  {
1259  Entry *entry = allocate(i, j);
1260  entry->value = value;
1261  }
1262 }
1263 
1264 
1265 
1266 template <typename number>
1267 inline void
1269  const size_type j,
1270  const number value)
1271 {
1272  AssertIsFinite(value);
1273 
1274  Assert(i < m(), ExcIndexRange(i, 0, m()));
1275  Assert(j < n(), ExcIndexRange(j, 0, n()));
1276 
1277  // ignore zero additions
1278  if (std::abs(value) == 0.)
1279  return;
1280 
1281  Entry *entry = allocate(i, j);
1282  entry->value += value;
1283 }
1284 
1285 
1286 template <typename number>
1287 template <typename number2>
1288 void
1289 SparseMatrixEZ<number>::add(const std::vector<size_type> &indices,
1290  const FullMatrix<number2> & full_matrix,
1291  const bool elide_zero_values)
1292 {
1293  // TODO: This function can surely be made more efficient
1294  for (size_type i = 0; i < indices.size(); ++i)
1295  for (size_type j = 0; j < indices.size(); ++j)
1296  if ((full_matrix(i, j) != 0) || (elide_zero_values == false))
1297  add(indices[i], indices[j], full_matrix(i, j));
1298 }
1299 
1300 
1301 
1302 template <typename number>
1303 template <typename number2>
1304 void
1305 SparseMatrixEZ<number>::add(const std::vector<size_type> &row_indices,
1306  const std::vector<size_type> &col_indices,
1307  const FullMatrix<number2> & full_matrix,
1308  const bool elide_zero_values)
1309 {
1310  // TODO: This function can surely be made more efficient
1311  for (size_type i = 0; i < row_indices.size(); ++i)
1312  for (size_type j = 0; j < col_indices.size(); ++j)
1313  if ((full_matrix(i, j) != 0) || (elide_zero_values == false))
1314  add(row_indices[i], col_indices[j], full_matrix(i, j));
1315 }
1316 
1317 
1318 
1319 template <typename number>
1320 template <typename number2>
1321 void
1323  const std::vector<size_type> &col_indices,
1324  const std::vector<number2> & values,
1325  const bool elide_zero_values)
1326 {
1327  // TODO: This function can surely be made more efficient
1328  for (size_type j = 0; j < col_indices.size(); ++j)
1329  if ((values[j] != 0) || (elide_zero_values == false))
1330  add(row, col_indices[j], values[j]);
1331 }
1332 
1333 
1334 
1335 template <typename number>
1336 template <typename number2>
1337 void
1339  const size_type n_cols,
1340  const size_type *col_indices,
1341  const number2 * values,
1342  const bool elide_zero_values,
1343  const bool /*col_indices_are_sorted*/)
1344 {
1345  // TODO: This function can surely be made more efficient
1346  for (size_type j = 0; j < n_cols; ++j)
1347  if ((std::abs(values[j]) != 0) || (elide_zero_values == false))
1348  add(row, col_indices[j], values[j]);
1349 }
1350 
1351 
1352 
1353 template <typename number>
1354 inline number
1356 {
1357  const Entry *entry = locate(i, j);
1358  if (entry)
1359  return entry->value;
1360  return 0.;
1361 }
1362 
1363 
1364 
1365 template <typename number>
1366 inline number
1368 {
1369  const Entry *entry = locate(i, j);
1370  if (entry)
1371  return entry->value;
1372  Assert(false, ExcInvalidEntry(i, j));
1373  return 0.;
1374 }
1375 
1376 
1377 template <typename number>
1380 {
1381  const_iterator result(this, 0, 0);
1382  return result;
1383 }
1384 
1385 template <typename number>
1388 {
1389  return const_iterator(this, m(), 0);
1390 }
1391 
1392 template <typename number>
1395 {
1396  Assert(r < m(), ExcIndexRange(r, 0, m()));
1397  const_iterator result(this, r, 0);
1398  return result;
1399 }
1400 
1401 template <typename number>
1404 {
1405  Assert(r < m(), ExcIndexRange(r, 0, m()));
1406  const_iterator result(this, r + 1, 0);
1407  return result;
1408 }
1409 
1410 template <typename number>
1411 template <typename MatrixType>
1412 inline SparseMatrixEZ<number> &
1414  const bool elide_zero_values)
1415 {
1416  reinit(M.m(), M.n(), this->saved_default_row_length, this->increment);
1417 
1418  // loop over the elements of the argument matrix row by row, as suggested
1419  // in the documentation of the sparse matrix iterator class, and
1420  // copy them into the current object
1421  for (size_type row = 0; row < M.m(); ++row)
1422  {
1423  const typename MatrixType::const_iterator end_row = M.end(row);
1424  for (typename MatrixType::const_iterator entry = M.begin(row);
1425  entry != end_row;
1426  ++entry)
1427  set(row, entry->column(), entry->value(), elide_zero_values);
1428  }
1429 
1430  return *this;
1431 }
1432 
1433 template <typename number>
1434 template <typename MatrixType>
1435 inline void
1436 SparseMatrixEZ<number>::add(const number factor, const MatrixType &M)
1437 {
1438  Assert(M.m() == m(), ExcDimensionMismatch(M.m(), m()));
1439  Assert(M.n() == n(), ExcDimensionMismatch(M.n(), n()));
1440 
1441  if (factor == 0.)
1442  return;
1443 
1444  // loop over the elements of the argument matrix row by row, as suggested
1445  // in the documentation of the sparse matrix iterator class, and
1446  // add them into the current object
1447  for (size_type row = 0; row < M.m(); ++row)
1448  {
1449  const typename MatrixType::const_iterator end_row = M.end(row);
1450  for (typename MatrixType::const_iterator entry = M.begin(row);
1451  entry != end_row;
1452  ++entry)
1453  if (entry->value() != 0)
1454  add(row, entry->column(), factor * entry->value());
1455  }
1456 }
1457 
1458 
1459 
1460 template <typename number>
1461 template <typename MatrixTypeA, typename MatrixTypeB>
1462 inline void
1464  const MatrixTypeB &B,
1465  const bool transpose)
1466 {
1467  // Compute the result
1468  // r_ij = \sum_kl b_ik b_jl a_kl
1469 
1470  // Assert (n() == B.m(), ExcDimensionMismatch(n(), B.m()));
1471  // Assert (m() == B.m(), ExcDimensionMismatch(m(), B.m()));
1472  // Assert (A.n() == B.n(), ExcDimensionMismatch(A.n(), B.n()));
1473  // Assert (A.m() == B.n(), ExcDimensionMismatch(A.m(), B.n()));
1474 
1475  // Somehow, we have to avoid making
1476  // this an operation of complexity
1477  // n^2. For the transpose case, we
1478  // can go through the non-zero
1479  // elements of A^-1 and use the
1480  // corresponding rows of B only.
1481  // For the non-transpose case, we
1482  // must find a trick.
1483  typename MatrixTypeB::const_iterator b1 = B.begin();
1484  const typename MatrixTypeB::const_iterator b_final = B.end();
1485  if (transpose)
1486  while (b1 != b_final)
1487  {
1488  const size_type i = b1->column();
1489  const size_type k = b1->row();
1490  typename MatrixTypeB::const_iterator b2 = B.begin();
1491  while (b2 != b_final)
1492  {
1493  const size_type j = b2->column();
1494  const size_type l = b2->row();
1495 
1496  const typename MatrixTypeA::value_type a = A.el(k, l);
1497 
1498  if (a != 0.)
1499  add(i, j, a * b1->value() * b2->value());
1500  ++b2;
1501  }
1502  ++b1;
1503  }
1504  else
1505  {
1506  // Determine minimal and
1507  // maximal row for a column in
1508  // advance.
1509 
1510  std::vector<size_type> minrow(B.n(), B.m());
1511  std::vector<size_type> maxrow(B.n(), 0);
1512  while (b1 != b_final)
1513  {
1514  const size_type r = b1->row();
1515  if (r < minrow[b1->column()])
1516  minrow[b1->column()] = r;
1517  if (r > maxrow[b1->column()])
1518  maxrow[b1->column()] = r;
1519  ++b1;
1520  }
1521 
1522  typename MatrixTypeA::const_iterator ai = A.begin();
1523  const typename MatrixTypeA::const_iterator ae = A.end();
1524 
1525  while (ai != ae)
1526  {
1527  const typename MatrixTypeA::value_type a = ai->value();
1528  // Don't do anything if
1529  // this entry is zero.
1530  if (a == 0.)
1531  continue;
1532 
1533  // Now, loop over all rows
1534  // having possibly a
1535  // nonzero entry in column
1536  // ai->row()
1537  b1 = B.begin(minrow[ai->row()]);
1538  const typename MatrixTypeB::const_iterator be1 =
1539  B.end(maxrow[ai->row()]);
1540  const typename MatrixTypeB::const_iterator be2 =
1541  B.end(maxrow[ai->column()]);
1542 
1543  while (b1 != be1)
1544  {
1545  const double b1v = b1->value();
1546  // We need the product
1547  // of both. If it is
1548  // zero, we can save
1549  // the work
1550  if (b1->column() == ai->row() && (b1v != 0.))
1551  {
1552  const size_type i = b1->row();
1553 
1554  typename MatrixTypeB::const_iterator b2 =
1555  B.begin(minrow[ai->column()]);
1556  while (b2 != be2)
1557  {
1558  if (b2->column() == ai->column())
1559  {
1560  const size_type j = b2->row();
1561  add(i, j, a * b1v * b2->value());
1562  }
1563  ++b2;
1564  }
1565  }
1566  ++b1;
1567  }
1568  ++ai;
1569  }
1570  }
1571 }
1572 
1573 
1574 template <typename number>
1575 template <class StreamType>
1576 inline void
1577 SparseMatrixEZ<number>::print_statistics(StreamType &out, bool full)
1578 {
1579  size_type used;
1580  size_type allocated;
1581  size_type reserved;
1582  std::vector<size_type> used_by_line;
1583 
1584  compute_statistics(used, allocated, reserved, used_by_line, full);
1585 
1586  out << "SparseMatrixEZ:used entries:" << used << std::endl
1587  << "SparseMatrixEZ:allocated entries:" << allocated << std::endl
1588  << "SparseMatrixEZ:reserved entries:" << reserved << std::endl;
1589 
1590  if (full)
1591  {
1592  for (size_type i = 0; i < used_by_line.size(); ++i)
1593  if (used_by_line[i] != 0)
1594  out << "SparseMatrixEZ:entries\t" << i << "\trows\t"
1595  << used_by_line[i] << std::endl;
1596  }
1597 }
1598 
1599 
1600 DEAL_II_NAMESPACE_CLOSE
1601 
1602 #endif
1603 /*---------------------------- sparse_matrix.h ---------------------------*/
size_type get_row_length(const size_type row) const
const types::global_dof_index invalid_size_type
Definition: types.h:182
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:541
void add(const size_type i, const size_type j, const number value)
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
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)
void conjugate_add(const MatrixTypeA &A, const MatrixTypeB &B, const bool transpose=false)
size_type n() const
void block_read(std::istream &in)
static ::ExceptionBase & ExcIndexRange(int arg1, int arg2, int arg3)
bool operator<(const const_iterator &) const
void vmult_add(Vector< somenumber > &dst, const Vector< somenumber > &src) const
LinearAlgebra::distributed::Vector< Number > Vector
void vmult(Vector< somenumber > &dst, const Vector< somenumber > &src) const
static const size_type invalid
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:1407
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
const_iterator(const SparseMatrixEZ< number > *matrix, const size_type row, const unsigned short index)
#define DeclException0(Exception0)
Definition: exceptions.h:473
const SparseMatrixEZ< number > * matrix
std::vector< Entry > data
Entry * allocate(const size_type row, const size_type col)
SymmetricTensor< rank_, dim, Number > transpose(const SymmetricTensor< rank_, dim, Number > &t)
RowInfo(const size_type start=Entry::invalid)
void print_statistics(StreamType &s, bool full=false)
static ::ExceptionBase & ExcIteratorPastEnd()
void swap(Vector< Number > &u, Vector< Number > &v)
Definition: vector.h:1376
unsigned int global_dof_index
Definition: types.h:89
void block_write(std::ostream &out) const
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
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
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
types::global_dof_index size_type
#define AssertIsFinite(number)
Definition: exceptions.h:1669
unsigned int saved_default_row_length
size_type m() const
std::vector< RowInfo > row_info
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