Reference documentation for deal.II version Git 7026f387cc 2019-10-15 14:19:01 -0400
\(\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 // Forward declarations
32 # ifndef DOXYGEN
33 template <typename number>
34 class Vector;
35 template <typename number>
36 class FullMatrix;
37 # endif
38 
106 template <typename number>
108 {
109 public:
114 
119  struct Entry
120  {
124  Entry();
125 
129  Entry(const size_type column, const number &value);
130 
135 
139  number value;
140 
145  };
146 
151  struct RowInfo
152  {
156  RowInfo(const size_type start = Entry::invalid);
157 
165  unsigned short length;
169  unsigned short diagonal;
173  static const unsigned short invalid_diagonal =
174  static_cast<unsigned short>(-1);
175  };
176 
177 public:
182  {
183  private:
187  class Accessor
188  {
189  public:
194  Accessor(const SparseMatrixEZ<number> *matrix,
195  const size_type row,
196  const unsigned short index);
197 
201  size_type
202  row() const;
203 
207  unsigned short
208  index() const;
209 
213  size_type
214  column() const;
215 
219  number
220  value() const;
221 
222  protected:
227 
232 
236  unsigned short a_index;
237 
238  // Make enclosing class a friend.
239  friend class const_iterator;
240  };
241 
242  public:
247  const size_type row,
248  const unsigned short index);
249 
254  operator++();
255 
259  const Accessor &operator*() const;
260 
264  const Accessor *operator->() const;
265 
269  bool
270  operator==(const const_iterator &) const;
274  bool
275  operator!=(const const_iterator &) const;
276 
281  bool
282  operator<(const const_iterator &) const;
283 
284  private:
289  };
290 
295  using value_type = number;
296 
304  SparseMatrixEZ();
305 
314 
321  explicit SparseMatrixEZ(const size_type n_rows,
322  const size_type n_columns,
323  const size_type default_row_length = 0,
324  const unsigned int default_increment = 1);
325 
329  ~SparseMatrixEZ() override = default;
330 
336 
346  operator=(const double d);
347 
355  void
356  reinit(const size_type n_rows,
357  const size_type n_columns,
358  size_type default_row_length = 0,
359  unsigned int default_increment = 1,
360  size_type reserve = 0);
361 
366  void
367  clear();
369 
377  bool
378  empty() const;
379 
384  size_type
385  m() const;
386 
391  size_type
392  n() const;
393 
397  size_type
398  get_row_length(const size_type row) const;
399 
403  size_type
404  n_nonzero_elements() const;
405 
410  std::size_t
411  memory_consumption() const;
412 
418  template <class StreamType>
419  void
420  print_statistics(StreamType &s, bool full = false);
421 
431  void
433  size_type & allocated,
434  size_type & reserved,
435  std::vector<size_type> &used_by_line,
436  const bool compute_by_line) const;
438 
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 
469  void
470  add(const size_type i, const size_type j, const number value);
471 
486  template <typename number2>
487  void
488  add(const std::vector<size_type> &indices,
489  const FullMatrix<number2> & full_matrix,
490  const bool elide_zero_values = true);
491 
497  template <typename number2>
498  void
499  add(const std::vector<size_type> &row_indices,
500  const std::vector<size_type> &col_indices,
501  const FullMatrix<number2> & full_matrix,
502  const bool elide_zero_values = true);
503 
513  template <typename number2>
514  void
515  add(const size_type row,
516  const std::vector<size_type> &col_indices,
517  const std::vector<number2> & values,
518  const bool elide_zero_values = true);
519 
529  template <typename number2>
530  void
531  add(const size_type row,
532  const size_type n_cols,
533  const size_type *col_indices,
534  const number2 * values,
535  const bool elide_zero_values = true,
536  const bool col_indices_are_sorted = false);
537 
559  template <typename MatrixType>
561  copy_from(const MatrixType &source, const bool elide_zero_values = true);
562 
570  template <typename MatrixType>
571  void
572  add(const number factor, const MatrixType &matrix);
574 
587  number
588  operator()(const size_type i, const size_type j) const;
589 
594  number
595  el(const size_type i, const size_type j) const;
597 
605  template <typename somenumber>
606  void
607  vmult(Vector<somenumber> &dst, const Vector<somenumber> &src) const;
608 
614  template <typename somenumber>
615  void
616  Tvmult(Vector<somenumber> &dst, const Vector<somenumber> &src) const;
617 
622  template <typename somenumber>
623  void
624  vmult_add(Vector<somenumber> &dst, const Vector<somenumber> &src) const;
625 
631  template <typename somenumber>
632  void
633  Tvmult_add(Vector<somenumber> &dst, const Vector<somenumber> &src) const;
635 
642  number
643  l2_norm() const;
645 
654  template <typename somenumber>
655  void
656  precondition_Jacobi(Vector<somenumber> & dst,
657  const Vector<somenumber> &src,
658  const number omega = 1.) const;
659 
663  template <typename somenumber>
664  void
665  precondition_SSOR(Vector<somenumber> & dst,
666  const Vector<somenumber> & src,
667  const number om = 1.,
668  const std::vector<std::size_t> &pos_right_of_diagonal =
669  std::vector<std::size_t>()) const;
670 
675  template <typename somenumber>
676  void
677  precondition_SOR(Vector<somenumber> & dst,
678  const Vector<somenumber> &src,
679  const number om = 1.) const;
680 
685  template <typename somenumber>
686  void
687  precondition_TSOR(Vector<somenumber> & dst,
688  const Vector<somenumber> &src,
689  const number om = 1.) const;
690 
699  template <typename MatrixTypeA, typename MatrixTypeB>
700  void
701  conjugate_add(const MatrixTypeA &A,
702  const MatrixTypeB &B,
703  const bool transpose = false);
705 
712  const_iterator
713  begin() const;
714 
718  const_iterator
719  end() const;
720 
725  const_iterator
726  begin(const size_type r) const;
727 
732  const_iterator
733  end(const size_type r) const;
735 
743  void
744  print(std::ostream &out) const;
745 
766  void
767  print_formatted(std::ostream & out,
768  const unsigned int precision = 3,
769  const bool scientific = true,
770  const unsigned int width = 0,
771  const char * zero_string = " ",
772  const double denominator = 1.) const;
773 
779  void
780  block_write(std::ostream &out) const;
781 
792  void
793  block_read(std::istream &in);
795 
805 
810  int,
811  int,
812  << "The entry with index (" << arg1 << ',' << arg2
813  << ") does not exist.");
814 
816  int,
817  int,
818  << "An entry with index (" << arg1 << ',' << arg2
819  << ") cannot be allocated.");
821 private:
826  const Entry *
827  locate(const size_type row, const size_type col) const;
828 
833  Entry *
834  locate(const size_type row, const size_type col);
835 
839  Entry *
840  allocate(const size_type row, const size_type col);
841 
847  template <typename somenumber>
848  void
849  threaded_vmult(Vector<somenumber> & dst,
850  const Vector<somenumber> &src,
851  const size_type begin_row,
852  const size_type end_row) const;
853 
859  template <typename somenumber>
860  void
861  threaded_matrix_norm_square(const Vector<somenumber> &v,
862  const size_type begin_row,
863  const size_type end_row,
864  somenumber * partial_sum) const;
865 
871  template <typename somenumber>
872  void
873  threaded_matrix_scalar_product(const Vector<somenumber> &u,
874  const Vector<somenumber> &v,
875  const size_type begin_row,
876  const size_type end_row,
877  somenumber * partial_sum) const;
878 
883 
887  std::vector<RowInfo> row_info;
888 
892  std::vector<Entry> data;
893 
897  unsigned int increment;
898 
903 };
904 
908 /*---------------------- Inline functions -----------------------------------*/
909 
910 template <typename number>
912  const number & value)
913  : column(column)
914  , value(value)
915 {}
916 
917 
918 
919 template <typename number>
921  : column(invalid)
922  , value(0)
923 {}
924 
925 
926 template <typename number>
928  : start(start)
929  , length(0)
930  , diagonal(invalid_diagonal)
931 {}
932 
933 
934 //---------------------------------------------------------------------------
935 template <typename number>
937  const SparseMatrixEZ<number> *matrix,
938  const size_type r,
939  const unsigned short i)
940  : matrix(matrix)
941  , a_row(r)
942  , a_index(i)
943 {}
944 
945 
946 template <typename number>
947 inline typename SparseMatrixEZ<number>::size_type
949 {
950  return a_row;
951 }
952 
953 
954 template <typename number>
955 inline typename SparseMatrixEZ<number>::size_type
957 {
958  return matrix->data[matrix->row_info[a_row].start + a_index].column;
959 }
960 
961 
962 template <typename number>
963 inline unsigned short
965 {
966  return a_index;
967 }
968 
969 
970 
971 template <typename number>
972 inline number
974 {
975  return matrix->data[matrix->row_info[a_row].start + a_index].value;
976 }
977 
978 
979 template <typename number>
982  const size_type r,
983  const unsigned short i)
984  : accessor(matrix, r, i)
985 {
986  // Finish if this is the end()
987  if (r == accessor.matrix->m() && i == 0)
988  return;
989 
990  // Make sure we never construct an
991  // iterator pointing to a
992  // non-existing entry
993 
994  // If the index points beyond the
995  // end of the row, try the next
996  // row.
997  if (accessor.a_index >= accessor.matrix->row_info[accessor.a_row].length)
998  {
999  do
1000  {
1001  ++accessor.a_row;
1002  }
1003  // Beware! If the next row is
1004  // empty, iterate until a
1005  // non-empty row is found or we
1006  // hit the end of the matrix.
1007  while (accessor.a_row < accessor.matrix->m() &&
1008  accessor.matrix->row_info[accessor.a_row].length == 0);
1009  }
1010 }
1011 
1012 
1013 template <typename number>
1016 {
1018 
1019  // Increment column index
1020  ++(accessor.a_index);
1021  // If index exceeds number of
1022  // entries in this row, proceed
1023  // with next row.
1024  if (accessor.a_index >= accessor.matrix->row_info[accessor.a_row].length)
1025  {
1026  accessor.a_index = 0;
1027  // Do this loop to avoid
1028  // elements in empty rows
1029  do
1030  {
1031  ++accessor.a_row;
1032  }
1033  while (accessor.a_row < accessor.matrix->m() &&
1034  accessor.matrix->row_info[accessor.a_row].length == 0);
1035  }
1036  return *this;
1037 }
1038 
1039 
1040 template <typename number>
1043 {
1044  return accessor;
1045 }
1046 
1047 
1048 template <typename number>
1051 {
1052  return &accessor;
1053 }
1054 
1055 
1056 template <typename number>
1057 inline bool
1059 operator==(const const_iterator &other) const
1060 {
1061  return (accessor.row() == other.accessor.row() &&
1062  accessor.index() == other.accessor.index());
1063 }
1064 
1065 
1066 template <typename number>
1067 inline bool
1069 operator!=(const const_iterator &other) const
1070 {
1071  return !(*this == other);
1072 }
1073 
1074 
1075 template <typename number>
1076 inline bool
1078 operator<(const const_iterator &other) const
1079 {
1080  return (accessor.row() < other.accessor.row() ||
1081  (accessor.row() == other.accessor.row() &&
1082  accessor.index() < other.accessor.index()));
1083 }
1084 
1085 
1086 //---------------------------------------------------------------------------
1087 template <typename number>
1088 inline typename SparseMatrixEZ<number>::size_type
1090 {
1091  return row_info.size();
1092 }
1093 
1094 
1095 template <typename number>
1096 inline typename SparseMatrixEZ<number>::size_type
1098 {
1099  return n_columns;
1100 }
1101 
1102 
1103 template <typename number>
1104 inline typename SparseMatrixEZ<number>::Entry *
1106 {
1107  Assert(row < m(), ExcIndexRange(row, 0, m()));
1108  Assert(col < n(), ExcIndexRange(col, 0, n()));
1109 
1110  const RowInfo & r = row_info[row];
1111  const size_type end = r.start + r.length;
1112  for (size_type i = r.start; i < end; ++i)
1113  {
1114  Entry *const entry = &data[i];
1115  if (entry->column == col)
1116  return entry;
1117  if (entry->column == Entry::invalid)
1118  return nullptr;
1119  }
1120  return nullptr;
1121 }
1122 
1123 
1124 
1125 template <typename number>
1126 inline const typename SparseMatrixEZ<number>::Entry *
1128 {
1129  SparseMatrixEZ<number> *t = const_cast<SparseMatrixEZ<number> *>(this);
1130  return t->locate(row, col);
1131 }
1132 
1133 
1134 template <typename number>
1135 inline typename SparseMatrixEZ<number>::Entry *
1137 {
1138  Assert(row < m(), ExcIndexRange(row, 0, m()));
1139  Assert(col < n(), ExcIndexRange(col, 0, n()));
1140 
1141  RowInfo & r = row_info[row];
1142  const size_type end = r.start + r.length;
1143 
1144  size_type i = r.start;
1145  // If diagonal exists and this
1146  // column is higher, start only
1147  // after diagonal.
1148  if (r.diagonal != RowInfo::invalid_diagonal && col >= row)
1149  i += r.diagonal;
1150  // Find position of entry
1151  while (i < end && data[i].column < col)
1152  ++i;
1153 
1154  // entry found
1155  if (i != end && data[i].column == col)
1156  return &data[i];
1157 
1158  // Now, we must insert the new
1159  // entry and move all successive
1160  // entries back.
1161 
1162  // If no more space is available
1163  // for this row, insert new
1164  // elements into the vector.
1165  // TODO:[GK] We should not extend this row if i<end
1166  if (row != row_info.size() - 1)
1167  {
1168  if (end >= row_info[row + 1].start)
1169  {
1170  // Failure if increment 0
1171  Assert(increment != 0, ExcEntryAllocationFailure(row, col));
1172 
1173  // Insert new entries
1174  data.insert(data.begin() + end, increment, Entry());
1175  // Update starts of
1176  // following rows
1177  for (size_type rn = row + 1; rn < row_info.size(); ++rn)
1178  row_info[rn].start += increment;
1179  }
1180  }
1181  else
1182  {
1183  if (end >= data.size())
1184  {
1185  // Here, appending a block
1186  // does not increase
1187  // performance.
1188  data.push_back(Entry());
1189  }
1190  }
1191 
1192  Entry *entry = &data[i];
1193  // Save original entry
1194  Entry temp = *entry;
1195  // Insert new entry here to
1196  // make sure all entries
1197  // are ordered by column
1198  // index
1199  entry->column = col;
1200  entry->value = 0;
1201  // Update row_info
1202  ++r.length;
1203  if (col == row)
1204  r.diagonal = i - r.start;
1205  else if (col < row && r.diagonal != RowInfo::invalid_diagonal)
1206  ++r.diagonal;
1207 
1208  if (i == end)
1209  return entry;
1210 
1211  // Move all entries in this
1212  // row up by one
1213  for (size_type j = i + 1; j < end; ++j)
1214  {
1215  // There should be no invalid
1216  // entry below end
1217  Assert(data[j].column != Entry::invalid, ExcInternalError());
1218 
1219  // TODO[GK]: This could be done more efficiently by moving starting at the
1220  // top rather than swapping starting at the bottom
1221  std::swap(data[j], temp);
1222  }
1223  Assert(data[end].column == Entry::invalid, ExcInternalError());
1224 
1225  data[end] = temp;
1226 
1227  return entry;
1228 }
1229 
1230 
1231 
1232 template <typename number>
1233 inline void
1235  const size_type j,
1236  const number value,
1237  const bool elide_zero_values)
1238 {
1239  AssertIsFinite(value);
1240 
1241  Assert(i < m(), ExcIndexRange(i, 0, m()));
1242  Assert(j < n(), ExcIndexRange(j, 0, n()));
1243 
1244  if (elide_zero_values && value == 0.)
1245  {
1246  Entry *entry = locate(i, j);
1247  if (entry != nullptr)
1248  entry->value = 0.;
1249  }
1250  else
1251  {
1252  Entry *entry = allocate(i, j);
1253  entry->value = value;
1254  }
1255 }
1256 
1257 
1258 
1259 template <typename number>
1260 inline void
1262  const size_type j,
1263  const number value)
1264 {
1265  AssertIsFinite(value);
1266 
1267  Assert(i < m(), ExcIndexRange(i, 0, m()));
1268  Assert(j < n(), ExcIndexRange(j, 0, n()));
1269 
1270  // ignore zero additions
1271  if (std::abs(value) == 0.)
1272  return;
1273 
1274  Entry *entry = allocate(i, j);
1275  entry->value += value;
1276 }
1277 
1278 
1279 template <typename number>
1280 template <typename number2>
1281 void
1282 SparseMatrixEZ<number>::add(const std::vector<size_type> &indices,
1283  const FullMatrix<number2> & full_matrix,
1284  const bool elide_zero_values)
1285 {
1286  // TODO: This function can surely be made more efficient
1287  for (size_type i = 0; i < indices.size(); ++i)
1288  for (size_type j = 0; j < indices.size(); ++j)
1289  if ((full_matrix(i, j) != 0) || (elide_zero_values == false))
1290  add(indices[i], indices[j], full_matrix(i, j));
1291 }
1292 
1293 
1294 
1295 template <typename number>
1296 template <typename number2>
1297 void
1298 SparseMatrixEZ<number>::add(const std::vector<size_type> &row_indices,
1299  const std::vector<size_type> &col_indices,
1300  const FullMatrix<number2> & full_matrix,
1301  const bool elide_zero_values)
1302 {
1303  // TODO: This function can surely be made more efficient
1304  for (size_type i = 0; i < row_indices.size(); ++i)
1305  for (size_type j = 0; j < col_indices.size(); ++j)
1306  if ((full_matrix(i, j) != 0) || (elide_zero_values == false))
1307  add(row_indices[i], col_indices[j], full_matrix(i, j));
1308 }
1309 
1310 
1311 
1312 template <typename number>
1313 template <typename number2>
1314 void
1316  const std::vector<size_type> &col_indices,
1317  const std::vector<number2> & values,
1318  const bool elide_zero_values)
1319 {
1320  // TODO: This function can surely be made more efficient
1321  for (size_type j = 0; j < col_indices.size(); ++j)
1322  if ((values[j] != 0) || (elide_zero_values == false))
1323  add(row, col_indices[j], values[j]);
1324 }
1325 
1326 
1327 
1328 template <typename number>
1329 template <typename number2>
1330 void
1332  const size_type n_cols,
1333  const size_type *col_indices,
1334  const number2 * values,
1335  const bool elide_zero_values,
1336  const bool /*col_indices_are_sorted*/)
1337 {
1338  // TODO: This function can surely be made more efficient
1339  for (size_type j = 0; j < n_cols; ++j)
1340  if ((std::abs(values[j]) != 0) || (elide_zero_values == false))
1341  add(row, col_indices[j], values[j]);
1342 }
1343 
1344 
1345 
1346 template <typename number>
1347 inline number
1349 {
1350  const Entry *entry = locate(i, j);
1351  if (entry)
1352  return entry->value;
1353  return 0.;
1354 }
1355 
1356 
1357 
1358 template <typename number>
1359 inline number
1361 {
1362  const Entry *entry = locate(i, j);
1363  if (entry)
1364  return entry->value;
1365  Assert(false, ExcInvalidEntry(i, j));
1366  return 0.;
1367 }
1368 
1369 
1370 template <typename number>
1373 {
1374  const_iterator result(this, 0, 0);
1375  return result;
1376 }
1377 
1378 template <typename number>
1381 {
1382  return const_iterator(this, m(), 0);
1383 }
1384 
1385 template <typename number>
1388 {
1389  Assert(r < m(), ExcIndexRange(r, 0, m()));
1390  const_iterator result(this, r, 0);
1391  return result;
1392 }
1393 
1394 template <typename number>
1397 {
1398  Assert(r < m(), ExcIndexRange(r, 0, m()));
1399  const_iterator result(this, r + 1, 0);
1400  return result;
1401 }
1402 
1403 template <typename number>
1404 template <typename MatrixType>
1405 inline SparseMatrixEZ<number> &
1407  const bool elide_zero_values)
1408 {
1409  reinit(M.m(), M.n(), this->saved_default_row_length, this->increment);
1410 
1411  // loop over the elements of the argument matrix row by row, as suggested
1412  // in the documentation of the sparse matrix iterator class, and
1413  // copy them into the current object
1414  for (size_type row = 0; row < M.m(); ++row)
1415  {
1416  const typename MatrixType::const_iterator end_row = M.end(row);
1417  for (typename MatrixType::const_iterator entry = M.begin(row);
1418  entry != end_row;
1419  ++entry)
1420  set(row, entry->column(), entry->value(), elide_zero_values);
1421  }
1422 
1423  return *this;
1424 }
1425 
1426 template <typename number>
1427 template <typename MatrixType>
1428 inline void
1429 SparseMatrixEZ<number>::add(const number factor, const MatrixType &M)
1430 {
1431  Assert(M.m() == m(), ExcDimensionMismatch(M.m(), m()));
1432  Assert(M.n() == n(), ExcDimensionMismatch(M.n(), n()));
1433 
1434  if (factor == 0.)
1435  return;
1436 
1437  // loop over the elements of the argument matrix row by row, as suggested
1438  // in the documentation of the sparse matrix iterator class, and
1439  // add them into the current object
1440  for (size_type row = 0; row < M.m(); ++row)
1441  {
1442  const typename MatrixType::const_iterator end_row = M.end(row);
1443  for (typename MatrixType::const_iterator entry = M.begin(row);
1444  entry != end_row;
1445  ++entry)
1446  if (entry->value() != 0)
1447  add(row, entry->column(), factor * entry->value());
1448  }
1449 }
1450 
1451 
1452 
1453 template <typename number>
1454 template <typename MatrixTypeA, typename MatrixTypeB>
1455 inline void
1457  const MatrixTypeB &B,
1458  const bool transpose)
1459 {
1460  // Compute the result
1461  // r_ij = \sum_kl b_ik b_jl a_kl
1462 
1463  // Assert (n() == B.m(), ExcDimensionMismatch(n(), B.m()));
1464  // Assert (m() == B.m(), ExcDimensionMismatch(m(), B.m()));
1465  // Assert (A.n() == B.n(), ExcDimensionMismatch(A.n(), B.n()));
1466  // Assert (A.m() == B.n(), ExcDimensionMismatch(A.m(), B.n()));
1467 
1468  // Somehow, we have to avoid making
1469  // this an operation of complexity
1470  // n^2. For the transpose case, we
1471  // can go through the non-zero
1472  // elements of A^-1 and use the
1473  // corresponding rows of B only.
1474  // For the non-transpose case, we
1475  // must find a trick.
1476  typename MatrixTypeB::const_iterator b1 = B.begin();
1477  const typename MatrixTypeB::const_iterator b_final = B.end();
1478  if (transpose)
1479  while (b1 != b_final)
1480  {
1481  const size_type i = b1->column();
1482  const size_type k = b1->row();
1483  typename MatrixTypeB::const_iterator b2 = B.begin();
1484  while (b2 != b_final)
1485  {
1486  const size_type j = b2->column();
1487  const size_type l = b2->row();
1488 
1489  const typename MatrixTypeA::value_type a = A.el(k, l);
1490 
1491  if (a != 0.)
1492  add(i, j, a * b1->value() * b2->value());
1493  ++b2;
1494  }
1495  ++b1;
1496  }
1497  else
1498  {
1499  // Determine minimal and
1500  // maximal row for a column in
1501  // advance.
1502 
1503  std::vector<size_type> minrow(B.n(), B.m());
1504  std::vector<size_type> maxrow(B.n(), 0);
1505  while (b1 != b_final)
1506  {
1507  const size_type r = b1->row();
1508  if (r < minrow[b1->column()])
1509  minrow[b1->column()] = r;
1510  if (r > maxrow[b1->column()])
1511  maxrow[b1->column()] = r;
1512  ++b1;
1513  }
1514 
1515  typename MatrixTypeA::const_iterator ai = A.begin();
1516  const typename MatrixTypeA::const_iterator ae = A.end();
1517 
1518  while (ai != ae)
1519  {
1520  const typename MatrixTypeA::value_type a = ai->value();
1521  // Don't do anything if
1522  // this entry is zero.
1523  if (a == 0.)
1524  continue;
1525 
1526  // Now, loop over all rows
1527  // having possibly a
1528  // nonzero entry in column
1529  // ai->row()
1530  b1 = B.begin(minrow[ai->row()]);
1531  const typename MatrixTypeB::const_iterator be1 =
1532  B.end(maxrow[ai->row()]);
1533  const typename MatrixTypeB::const_iterator be2 =
1534  B.end(maxrow[ai->column()]);
1535 
1536  while (b1 != be1)
1537  {
1538  const double b1v = b1->value();
1539  // We need the product
1540  // of both. If it is
1541  // zero, we can save
1542  // the work
1543  if (b1->column() == ai->row() && (b1v != 0.))
1544  {
1545  const size_type i = b1->row();
1546 
1547  typename MatrixTypeB::const_iterator b2 =
1548  B.begin(minrow[ai->column()]);
1549  while (b2 != be2)
1550  {
1551  if (b2->column() == ai->column())
1552  {
1553  const size_type j = b2->row();
1554  add(i, j, a * b1v * b2->value());
1555  }
1556  ++b2;
1557  }
1558  }
1559  ++b1;
1560  }
1561  ++ai;
1562  }
1563  }
1564 }
1565 
1566 
1567 template <typename number>
1568 template <class StreamType>
1569 inline void
1570 SparseMatrixEZ<number>::print_statistics(StreamType &out, bool full)
1571 {
1572  size_type used;
1573  size_type allocated;
1574  size_type reserved;
1575  std::vector<size_type> used_by_line;
1576 
1577  compute_statistics(used, allocated, reserved, used_by_line, full);
1578 
1579  out << "SparseMatrixEZ:used entries:" << used << std::endl
1580  << "SparseMatrixEZ:allocated entries:" << allocated << std::endl
1581  << "SparseMatrixEZ:reserved entries:" << reserved << std::endl;
1582 
1583  if (full)
1584  {
1585  for (size_type i = 0; i < used_by_line.size(); ++i)
1586  if (used_by_line[i] != 0)
1587  out << "SparseMatrixEZ:entries\t" << i << "\trows\t"
1588  << used_by_line[i] << std::endl;
1589  }
1590 }
1591 
1592 
1593 DEAL_II_NAMESPACE_CLOSE
1594 
1595 #endif
1596 /*---------------------------- sparse_matrix.h ---------------------------*/
size_type get_row_length(const size_type row) const
const types::global_dof_index invalid_size_type
Definition: types.h:196
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)
DEAL_II_CONSTEXPR SymmetricTensor< rank_, dim, Number > operator*(const SymmetricTensor< rank_, dim, Number > &t, const Number &factor)
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)
constexpr SymmetricTensor< rank_, dim, Number > transpose(const SymmetricTensor< rank_, dim, Number > &t)
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)
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