Reference documentation for deal.II version 8.1.0
sparse_matrix_ez.h
1 // ---------------------------------------------------------------------
2 // @f$Id: sparse_matrix_ez.h 30040 2013-07-18 17:06:48Z maier @f$
3 //
4 // Copyright (C) 2002 - 2013 by the deal.II authors
5 //
6 // This file is part of the deal.II library.
7 //
8 // The deal.II library is free software; you can use it, redistribute
9 // it, and/or modify it under the terms of the GNU Lesser General
10 // Public License as published by the Free Software Foundation; either
11 // version 2.1 of the License, or (at your option) any later version.
12 // The full text of the license can be found in the file LICENSE at
13 // the top level of the deal.II distribution.
14 //
15 // ---------------------------------------------------------------------
16 
17 #ifndef __deal2__sparse_matrix_ez_h
18 #define __deal2__sparse_matrix_ez_h
19 
20 
21 #include <deal.II/base/config.h>
22 #include <deal.II/base/subscriptor.h>
23 #include <deal.II/base/smartpointer.h>
24 #include <deal.II/lac/exceptions.h>
25 
26 #include <vector>
27 
29 
30 template<typename number> class Vector;
31 template<typename number> class FullMatrix;
32 
101 template <typename number>
102 class SparseMatrixEZ : public Subscriptor
103 {
104 public:
109 
115  struct Entry
116  {
122  Entry();
123 
128  Entry (const size_type column,
129  const number &value);
130 
134  size_type column;
135 
139  number value;
140 
144  static const size_type invalid = numbers::invalid_size_type;
145  };
146 
153  struct RowInfo
154  {
158  RowInfo (size_type start = Entry::invalid);
159 
164  size_type start;
169  unsigned short length;
175  unsigned short diagonal;
179  static const unsigned short
180  invalid_diagonal = static_cast<unsigned short>(-1);
181  };
182 
183 public:
184 
189  {
190  private:
194  class Accessor
195  {
196  public:
204  const size_type row,
205  const unsigned short index);
206 
212  size_type row() const;
213 
219  unsigned short index() const;
220 
226  size_type column() const;
227 
231  number value() const;
232 
233  protected:
238 
242  size_type a_row;
243 
247  unsigned short a_index;
248 
253  friend class const_iterator;
254  };
255 
256  public:
261  const size_type row,
262  const unsigned short index);
263 
270 
277 
281  const Accessor &operator* () const;
282 
286  const Accessor *operator-> () const;
287 
294  bool operator == (const const_iterator &) const;
298  bool operator != (const const_iterator &) const;
299 
309  bool operator < (const const_iterator &) const;
310 
311  private:
317 
333  };
334 
339  typedef number value_type;
340 
350  SparseMatrixEZ ();
351 
362  SparseMatrixEZ (const SparseMatrixEZ &);
363 
377  explicit SparseMatrixEZ (const size_type n_rows,
378  const size_type n_columns,
379  const size_type default_row_length = 0,
380  const unsigned int default_increment = 1);
381 
387  ~SparseMatrixEZ ();
388 
394 
411  SparseMatrixEZ<number> &operator = (const double d);
412 
428  void reinit (const size_type n_rows,
429  const size_type n_columns,
430  size_type default_row_length = 0,
431  unsigned int default_increment = 1,
432  size_type reserve = 0);
433 
441  void clear ();
443 
452  bool empty () const;
453 
460  size_type m () const;
461 
468  size_type n () const;
469 
474  size_type get_row_length (const size_type row) const;
475 
480  size_type n_nonzero_elements () const;
481 
487  std::size_t memory_consumption () const;
488 
498  template <class STREAM>
499  void print_statistics (STREAM &s, bool full = false);
500 
514  void compute_statistics (size_type &used,
515  size_type &allocated,
516  size_type &reserved,
517  std::vector<size_type> &used_by_line,
518  const bool compute_by_line) const;
520 
533  void set (const size_type i, const size_type j,
534  const number value);
535 
545  void add (const size_type i,
546  const size_type j,
547  const number value);
548 
576  template <typename number2>
577  void add (const std::vector<size_type> &indices,
578  const FullMatrix<number2> &full_matrix,
579  const bool elide_zero_values = true);
580 
588  template <typename number2>
589  void add (const std::vector<size_type> &row_indices,
590  const std::vector<size_type> &col_indices,
591  const FullMatrix<number2> &full_matrix,
592  const bool elide_zero_values = true);
593 
611  template <typename number2>
612  void add (const size_type row,
613  const std::vector<size_type> &col_indices,
614  const std::vector<number2> &values,
615  const bool elide_zero_values = true);
616 
634  template <typename number2>
635  void add (const size_type row,
636  const size_type n_cols,
637  const size_type *col_indices,
638  const number2 *values,
639  const bool elide_zero_values = true,
640  const bool col_indices_are_sorted = false);
641 
668  template <class MATRIX>
670  copy_from (const MATRIX &source);
671 
683  template <class MATRIX>
684  void add (const number factor,
685  const MATRIX &matrix);
687 
709  number operator () (const size_type i,
710  const size_type j) const;
711 
717  number el (const size_type i,
718  const size_type j) const;
720 
729  template <typename somenumber>
730  void vmult (Vector<somenumber> &dst,
731  const Vector<somenumber> &src) const;
732 
741  template <typename somenumber>
742  void Tvmult (Vector<somenumber> &dst,
743  const Vector<somenumber> &src) const;
744 
751  template <typename somenumber>
752  void vmult_add (Vector<somenumber> &dst,
753  const Vector<somenumber> &src) const;
754 
763  template <typename somenumber>
764  void Tvmult_add (Vector<somenumber> &dst,
765  const Vector<somenumber> &src) const;
767 
774  number l2_norm () const;
776 
790  template <typename somenumber>
792  const Vector<somenumber> &src,
793  const number omega = 1.) const;
794 
799  template <typename somenumber>
801  const Vector<somenumber> &src,
802  const number om = 1.,
803  const std::vector<std::size_t> &pos_right_of_diagonal = std::vector<std::size_t>()) const;
804 
810  template <typename somenumber>
812  const Vector<somenumber> &src,
813  const number om = 1.) const;
814 
820  template <typename somenumber>
822  const Vector<somenumber> &src,
823  const number om = 1.) const;
824 
839  template <class MATRIXA, class MATRIXB>
840  void conjugate_add (const MATRIXA &A,
841  const MATRIXB &B,
842  const bool transpose = false);
844 
852  const_iterator begin () const;
853 
857  const_iterator end () const;
858 
866  const_iterator begin (const size_type r) const;
867 
873  const_iterator end (const size_type r) const;
875 
886  void print (std::ostream &out) const;
887 
927  void print_formatted (std::ostream &out,
928  const unsigned int precision = 3,
929  const bool scientific = true,
930  const unsigned int width = 0,
931  const char *zero_string = " ",
932  const double denominator = 1.) const;
933 
941  void block_write (std::ostream &out) const;
942 
960  void block_read (std::istream &in);
962 
969  DeclException0(ExcNoDiagonal);
970 
974  DeclException2 (ExcInvalidEntry,
975  int, int,
976  << "The entry with index (" << arg1 << ',' << arg2
977  << ") does not exist.");
978 
979  DeclException2(ExcEntryAllocationFailure,
980  int, int,
981  << "An entry with index (" << arg1 << ',' << arg2
982  << ") cannot be allocated.");
984 private:
991  const Entry *locate (const size_type row,
992  const size_type col) const;
993 
1000  Entry *locate (const size_type row,
1001  const size_type col);
1002 
1006  Entry *allocate (const size_type row,
1007  const size_type col);
1008 
1018  template <typename somenumber>
1020  const Vector<somenumber> &src,
1021  const size_type begin_row,
1022  const size_type end_row) const;
1023 
1035  template <typename somenumber>
1037  const size_type begin_row,
1038  const size_type end_row,
1039  somenumber *partial_sum) const;
1040 
1052  template <typename somenumber>
1054  const Vector<somenumber> &v,
1055  const size_type begin_row,
1056  const size_type end_row,
1057  somenumber *partial_sum) const;
1058 
1064  size_type n_columns;
1065 
1069  std::vector<RowInfo> row_info;
1070 
1074  std::vector<Entry> data;
1075 
1079  unsigned int increment;
1080 
1088 };
1089 
1093 /*---------------------- Inline functions -----------------------------------*/
1094 
1095 template <typename number>
1096 inline
1098  const number &value)
1099  :
1100  column(column),
1101  value(value)
1102 {}
1103 
1104 
1105 
1106 template <typename number>
1107 inline
1109  :
1110  column(invalid),
1111  value(0)
1112 {}
1113 
1114 
1115 template <typename number>
1116 inline
1118  :
1119  start(start),
1120  length(0),
1121  diagonal(invalid_diagonal)
1122 {}
1123 
1124 
1125 //---------------------------------------------------------------------------
1126 template <typename number>
1127 inline
1130  const size_type r,
1131  const unsigned short i)
1132  :
1133  matrix(matrix),
1134  a_row(r),
1135  a_index(i)
1136 {}
1137 
1138 
1139 template <typename number>
1140 inline
1143 {
1144  return a_row;
1145 }
1146 
1147 
1148 template <typename number>
1149 inline
1152 {
1153  return matrix->data[matrix->row_info[a_row].start+a_index].column;
1154 }
1155 
1156 
1157 template <typename number>
1158 inline
1159 unsigned short
1161 {
1162  return a_index;
1163 }
1164 
1165 
1166 
1167 template <typename number>
1168 inline
1169 number
1171 {
1172  return matrix->data[matrix->row_info[a_row].start+a_index].value;
1173 }
1174 
1175 
1176 template <typename number>
1177 inline
1180  const size_type r,
1181  const unsigned short i)
1182  :
1183  accessor(matrix, r, i)
1184 {
1185  // Finish if this is the end()
1186  if (r==accessor.matrix->m() && i==0) return;
1187 
1188  // Make sure we never construct an
1189  // iterator pointing to a
1190  // non-existing entry
1191 
1192  // If the index points beyond the
1193  // end of the row, try the next
1194  // row.
1195  if (accessor.a_index >= accessor.matrix->row_info[accessor.a_row].length)
1196  {
1197  do
1198  {
1199  ++accessor.a_row;
1200  }
1201  // Beware! If the next row is
1202  // empty, iterate until a
1203  // non-empty row is found or we
1204  // hit the end of the matrix.
1205  while (accessor.a_row < accessor.matrix->m()
1206  && accessor.matrix->row_info[accessor.a_row].length == 0);
1207  }
1208 }
1209 
1210 
1211 template <typename number>
1212 inline
1215 {
1216  Assert (accessor.a_row < accessor.matrix->m(), ExcIteratorPastEnd());
1217 
1218  // Increment column index
1219  ++(accessor.a_index);
1220  // If index exceeds number of
1221  // entries in this row, proceed
1222  // with next row.
1223  if (accessor.a_index >= accessor.matrix->row_info[accessor.a_row].length)
1224  {
1225  accessor.a_index = 0;
1226  // Do this loop to avoid
1227  // elements in empty rows
1228  do
1229  {
1230  ++accessor.a_row;
1231  }
1232  while (accessor.a_row < accessor.matrix->m()
1233  && accessor.matrix->row_info[accessor.a_row].length == 0);
1234  }
1235  return *this;
1236 }
1237 
1238 
1239 template <typename number>
1240 inline
1243 {
1244  return accessor;
1245 }
1246 
1247 
1248 template <typename number>
1249 inline
1252 {
1253  return &accessor;
1254 }
1255 
1256 
1257 template <typename number>
1258 inline
1259 bool
1261  const const_iterator &other) const
1262 {
1263  return (accessor.row() == other.accessor.row() &&
1264  accessor.index() == other.accessor.index());
1265 }
1266 
1267 
1268 template <typename number>
1269 inline
1270 bool
1272 operator != (const const_iterator &other) const
1273 {
1274  return ! (*this == other);
1275 }
1276 
1277 
1278 template <typename number>
1279 inline
1280 bool
1282 operator < (const const_iterator &other) const
1283 {
1284  return (accessor.row() < other.accessor.row() ||
1285  (accessor.row() == other.accessor.row() &&
1286  accessor.index() < other.accessor.index()));
1287 }
1288 
1289 
1290 //---------------------------------------------------------------------------
1291 template <typename number>
1292 inline
1294 {
1295  return row_info.size();
1296 }
1297 
1298 
1299 template <typename number>
1300 inline
1302 {
1303  return n_columns;
1304 }
1305 
1306 
1307 template <typename number>
1308 inline
1310 SparseMatrixEZ<number>::locate (const size_type row,
1311  const size_type col)
1312 {
1313  Assert (row<m(), ExcIndexRange(row,0,m()));
1314  Assert (col<n(), ExcIndexRange(col,0,n()));
1315 
1316  const RowInfo &r = row_info[row];
1317  const size_type end = r.start + r.length;
1318  for (size_type i=r.start; i<end; ++i)
1319  {
1320  Entry *const entry = &data[i];
1321  if (entry->column == col)
1322  return entry;
1323  if (entry->column == Entry::invalid)
1324  return 0;
1325  }
1326  return 0;
1327 }
1328 
1329 
1330 
1331 template <typename number>
1332 inline
1333 const typename SparseMatrixEZ<number>::Entry *
1334 SparseMatrixEZ<number>::locate (const size_type row,
1335  const size_type col) const
1336 {
1337  SparseMatrixEZ<number> *t = const_cast<SparseMatrixEZ<number>*> (this);
1338  return t->locate(row,col);
1339 }
1340 
1341 
1342 template <typename number>
1343 inline
1346  const size_type col)
1347 {
1348  Assert (row<m(), ExcIndexRange(row,0,m()));
1349  Assert (col<n(), ExcIndexRange(col,0,n()));
1350 
1351  RowInfo &r = row_info[row];
1352  const size_type end = r.start + r.length;
1353 
1354  size_type i = r.start;
1355  // If diagonal exists and this
1356  // column is higher, start only
1357  // after diagonal.
1358  if (r.diagonal != RowInfo::invalid_diagonal && col >= row)
1359  i += r.diagonal;
1360  // Find position of entry
1361  while (i<end && data[i].column < col) ++i;
1362 
1363  // entry found
1364  if (i != end && data[i].column == col)
1365  return &data[i];
1366 
1367  // Now, we must insert the new
1368  // entry and move all successive
1369  // entries back.
1370 
1371  // If no more space is available
1372  // for this row, insert new
1373  // elements into the vector.
1374 //TODO:[GK] We should not extend this row if i<end
1375  if (row != row_info.size()-1)
1376  {
1377  if (end >= row_info[row+1].start)
1378  {
1379  // Failure if increment 0
1380  Assert(increment!=0,ExcEntryAllocationFailure(row,col));
1381 
1382  // Insert new entries
1383  data.insert(data.begin()+end, increment, Entry());
1384  // Update starts of
1385  // following rows
1386  for (size_type rn=row+1; rn<row_info.size(); ++rn)
1387  row_info[rn].start += increment;
1388  }
1389  }
1390  else
1391  {
1392  if (end >= data.size())
1393  {
1394  // Here, appending a block
1395  // does not increase
1396  // performance.
1397  data.push_back(Entry());
1398  }
1399  }
1400 
1401  Entry *entry = &data[i];
1402  // Save original entry
1403  Entry temp = *entry;
1404  // Insert new entry here to
1405  // make sure all entries
1406  // are ordered by column
1407  // index
1408  entry->column = col;
1409  entry->value = 0;
1410  // Update row_info
1411  ++r.length;
1412  if (col == row)
1413  r.diagonal = i - r.start;
1414  else if (col<row && r.diagonal!= RowInfo::invalid_diagonal)
1415  ++r.diagonal;
1416 
1417  if (i == end)
1418  return entry;
1419 
1420  // Move all entries in this
1421  // row up by one
1422  for (size_type j = i+1; j < end; ++j)
1423  {
1424  // There should be no invalid
1425  // entry below end
1426  Assert (data[j].column != Entry::invalid, ExcInternalError());
1427 
1428 //TODO[GK]: This could be done more efficiently by moving starting at the top rather than swapping starting at the bottom
1429  std::swap (data[j], temp);
1430  }
1431  Assert (data[end].column == Entry::invalid, ExcInternalError());
1432 
1433  data[end] = temp;
1434 
1435  return entry;
1436 }
1437 
1438 
1439 
1440 template <typename number>
1441 inline
1442 void SparseMatrixEZ<number>::set (const size_type i,
1443  const size_type j,
1444  const number value)
1445 {
1446 
1448 
1449  Assert (i<m(), ExcIndexRange(i,0,m()));
1450  Assert (j<n(), ExcIndexRange(j,0,n()));
1451 
1452  if (value == 0.)
1453  {
1454  Entry *entry = locate(i,j);
1455  if (entry != 0)
1456  entry->value = 0.;
1457  }
1458  else
1459  {
1460  Entry *entry = allocate(i,j);
1461  entry->value = value;
1462  }
1463 }
1464 
1465 
1466 
1467 template <typename number>
1468 inline
1469 void SparseMatrixEZ<number>::add (const size_type i,
1470  const size_type j,
1471  const number value)
1472 {
1473 
1475 
1476  Assert (i<m(), ExcIndexRange(i,0,m()));
1477  Assert (j<n(), ExcIndexRange(j,0,n()));
1478 
1479  // ignore zero additions
1480  if (value == 0)
1481  return;
1482 
1483  Entry *entry = allocate(i,j);
1484  entry->value += value;
1485 }
1486 
1487 
1488 template <typename number>
1489 template <typename number2>
1490 void SparseMatrixEZ<number>::add (const std::vector<size_type> &indices,
1491  const FullMatrix<number2> &full_matrix,
1492  const bool elide_zero_values)
1493 {
1494 //TODO: This function can surely be made more efficient
1495  for (size_type i=0; i<indices.size(); ++i)
1496  for (size_type j=0; j<indices.size(); ++j)
1497  if ((full_matrix(i,j) != 0) || (elide_zero_values == false))
1498  add (indices[i], indices[j], full_matrix(i,j));
1499 }
1500 
1501 
1502 
1503 template <typename number>
1504 template <typename number2>
1505 void SparseMatrixEZ<number>::add (const std::vector<size_type> &row_indices,
1506  const std::vector<size_type> &col_indices,
1507  const FullMatrix<number2> &full_matrix,
1508  const bool elide_zero_values)
1509 {
1510 //TODO: This function can surely be made more efficient
1511  for (size_type i=0; i<row_indices.size(); ++i)
1512  for (size_type j=0; j<col_indices.size(); ++j)
1513  if ((full_matrix(i,j) != 0) || (elide_zero_values == false))
1514  add (row_indices[i], col_indices[j], full_matrix(i,j));
1515 }
1516 
1517 
1518 
1519 
1520 template <typename number>
1521 template <typename number2>
1522 void SparseMatrixEZ<number>::add (const size_type row,
1523  const std::vector<size_type> &col_indices,
1524  const std::vector<number2> &values,
1525  const bool elide_zero_values)
1526 {
1527 //TODO: This function can surely be made more efficient
1528  for (size_type j=0; j<col_indices.size(); ++j)
1529  if ((values[j] != 0) || (elide_zero_values == false))
1530  add (row, col_indices[j], values[j]);
1531 }
1532 
1533 
1534 
1535 template <typename number>
1536 template <typename number2>
1537 void SparseMatrixEZ<number>::add (const size_type row,
1538  const size_type n_cols,
1539  const size_type *col_indices,
1540  const number2 *values,
1541  const bool elide_zero_values,
1542  const bool /*col_indices_are_sorted*/)
1543 {
1544 //TODO: This function can surely be made more efficient
1545  for (size_type j=0; j<n_cols; ++j)
1546  if ((values[j] != 0) || (elide_zero_values == false))
1547  add (row, col_indices[j], values[j]);
1548 }
1549 
1550 
1551 
1552 
1553 template <typename number>
1554 inline
1555 number SparseMatrixEZ<number>::el (const size_type i,
1556  const size_type j) const
1557 {
1558  const Entry *entry = locate(i,j);
1559  if (entry)
1560  return entry->value;
1561  return 0.;
1562 }
1563 
1564 
1565 
1566 template <typename number>
1567 inline
1568 number SparseMatrixEZ<number>::operator() (const size_type i,
1569  const size_type j) const
1570 {
1571  const Entry *entry = locate(i,j);
1572  if (entry)
1573  return entry->value;
1574  Assert(false, ExcInvalidEntry(i,j));
1575  return 0.;
1576 }
1577 
1578 
1579 template <typename number>
1580 inline
1583 {
1584  const_iterator result(this, 0, 0);
1585  return result;
1586 }
1587 
1588 template <typename number>
1589 inline
1592 {
1593  return const_iterator(this, m(), 0);
1594 }
1595 
1596 template <typename number>
1597 inline
1599 SparseMatrixEZ<number>::begin (const size_type r) const
1600 {
1601  Assert (r<m(), ExcIndexRange(r,0,m()));
1602  const_iterator result (this, r, 0);
1603  return result;
1604 }
1605 
1606 template <typename number>
1607 inline
1609 SparseMatrixEZ<number>::end (const size_type r) const
1610 {
1611  Assert (r<m(), ExcIndexRange(r,0,m()));
1612  const_iterator result(this, r+1, 0);
1613  return result;
1614 }
1615 
1616 template<typename number>
1617 template <class MATRIX>
1618 inline
1621 {
1622  reinit(M.m(), M.n());
1623 
1624  // loop over the elements of the argument matrix row by row, as suggested
1625  // in the documentation of the sparse matrix iterator class, and
1626  // copy them into the current object
1627  for (size_type row = 0; row < M.m(); ++row)
1628  {
1629  const typename MATRIX::const_iterator end_row = M.end(row);
1630  for (typename MATRIX::const_iterator entry = M.begin(row);
1631  entry != end_row; ++entry)
1632  if (entry->value() != 0)
1633  set(row, entry->column(), entry->value());
1634  }
1635 
1636  return *this;
1637 }
1638 
1639 template<typename number>
1640 template <class MATRIX>
1641 inline
1642 void
1643 SparseMatrixEZ<number>::add (const number factor,
1644  const MATRIX &M)
1645 {
1646  Assert (M.m() == m(), ExcDimensionMismatch(M.m(), m()));
1647  Assert (M.n() == n(), ExcDimensionMismatch(M.n(), n()));
1648 
1649  if (factor == 0.)
1650  return;
1651 
1652  // loop over the elements of the argument matrix row by row, as suggested
1653  // in the documentation of the sparse matrix iterator class, and
1654  // add them into the current object
1655  for (size_type row = 0; row < M.m(); ++row)
1656  {
1657  const typename MATRIX::const_iterator end_row = M.end(row);
1658  for (typename MATRIX::const_iterator entry = M.begin(row);
1659  entry != end_row; ++entry)
1660  if (entry->value() != 0)
1661  add(row, entry->column(), factor * entry->value());
1662  }
1663 }
1664 
1665 
1666 
1667 template<typename number>
1668 template <class MATRIXA, class MATRIXB>
1669 inline void
1671  const MATRIXB &B,
1672  const bool transpose)
1673 {
1674 // Compute the result
1675 // r_ij = \sum_kl b_ik b_jl a_kl
1676 
1677 // Assert (n() == B.m(), ExcDimensionMismatch(n(), B.m()));
1678 // Assert (m() == B.m(), ExcDimensionMismatch(m(), B.m()));
1679 // Assert (A.n() == B.n(), ExcDimensionMismatch(A.n(), B.n()));
1680 // Assert (A.m() == B.n(), ExcDimensionMismatch(A.m(), B.n()));
1681 
1682  // Somehow, we have to avoid making
1683  // this an operation of complexity
1684  // n^2. For the transpose case, we
1685  // can go through the non-zero
1686  // elements of A^-1 and use the
1687  // corresponding rows of B only.
1688  // For the non-transpose case, we
1689  // must find a trick.
1690  typename MATRIXB::const_iterator b1 = B.begin();
1691  const typename MATRIXB::const_iterator b_final = B.end();
1692  if (transpose)
1693  while (b1 != b_final)
1694  {
1695  const size_type i = b1->column();
1696  const size_type k = b1->row();
1697  typename MATRIXB::const_iterator b2 = B.begin();
1698  while (b2 != b_final)
1699  {
1700  const size_type j = b2->column();
1701  const size_type l = b2->row();
1702 
1703  const typename MATRIXA::value_type a = A.el(k,l);
1704 
1705  if (a != 0.)
1706  add (i, j, a * b1->value() * b2->value());
1707  ++b2;
1708  }
1709  ++b1;
1710  }
1711  else
1712  {
1713  // Determine minimal and
1714  // maximal row for a column in
1715  // advance.
1716 
1717  std::vector<size_type> minrow(B.n(), B.m());
1718  std::vector<size_type> maxrow(B.n(), 0);
1719  while (b1 != b_final)
1720  {
1721  const size_type r = b1->row();
1722  if (r < minrow[b1->column()])
1723  minrow[b1->column()] = r;
1724  if (r > maxrow[b1->column()])
1725  maxrow[b1->column()] = r;
1726  ++b1;
1727  }
1728 
1729  typename MATRIXA::const_iterator ai = A.begin();
1730  const typename MATRIXA::const_iterator ae = A.end();
1731 
1732  while (ai != ae)
1733  {
1734  const typename MATRIXA::value_type a = ai->value();
1735  // Don't do anything if
1736  // this entry is zero.
1737  if (a == 0.) continue;
1738 
1739  // Now, loop over all rows
1740  // having possibly a
1741  // nonzero entry in column
1742  // ai->row()
1743  b1 = B.begin(minrow[ai->row()]);
1744  const typename MATRIXB::const_iterator
1745  be1 = B.end(maxrow[ai->row()]);
1746  const typename MATRIXB::const_iterator
1747  be2 = B.end(maxrow[ai->column()]);
1748 
1749  while (b1 != be1)
1750  {
1751  const double b1v = b1->value();
1752  // We need the product
1753  // of both. If it is
1754  // zero, we can save
1755  // the work
1756  if (b1->column() == ai->row() && (b1v != 0.))
1757  {
1758  const size_type i = b1->row();
1759 
1760  typename MATRIXB::const_iterator
1761  b2 = B.begin(minrow[ai->column()]);
1762  while (b2 != be2)
1763  {
1764  if (b2->column() == ai->column())
1765  {
1766  const size_type j = b2->row();
1767  add (i, j, a * b1v * b2->value());
1768  }
1769  ++b2;
1770  }
1771  }
1772  ++b1;
1773  }
1774  ++ai;
1775  }
1776  }
1777 }
1778 
1779 
1780 template <typename number>
1781 template <class STREAM>
1782 inline
1783 void
1785 {
1786  size_type used;
1787  size_type allocated;
1788  size_type reserved;
1789  std::vector<size_type> used_by_line;
1790 
1791  compute_statistics (used, allocated, reserved, used_by_line, full);
1792 
1793  out << "SparseMatrixEZ:used entries:" << used << std::endl
1794  << "SparseMatrixEZ:allocated entries:" << allocated << std::endl
1795  << "SparseMatrixEZ:reserved entries:" << reserved << std::endl;
1796 
1797  if (full)
1798  {
1799  for (size_type i=0; i< used_by_line.size(); ++i)
1800  if (used_by_line[i] != 0)
1801  out << "SparseMatrixEZ:entries\t" << i
1802  << "\trows\t" << used_by_line[i]
1803  << std::endl;
1804 
1805  }
1806 }
1807 
1808 
1809 DEAL_II_NAMESPACE_CLOSE
1810 
1811 #endif
1812 /*---------------------------- sparse_matrix.h ---------------------------*/
const types::global_dof_index invalid_size_type
Definition: types.h:211
void precondition_Jacobi(Vector< somenumber > &dst, const Vector< somenumber > &src, const number omega=1.) const
std::size_t memory_consumption() const
void add(const size_type i, const size_type j, const number value)
SparseMatrixEZ< number > & operator=(const SparseMatrixEZ< number > &)
unsigned int increment
const Accessor & operator*() const
void vmult(Vector< somenumber > &dst, const Vector< somenumber > &src) const
void threaded_vmult(Vector< somenumber > &dst, const Vector< somenumber > &src, const size_type begin_row, const size_type end_row) const
bool operator<(const const_iterator &) const
const Entry * locate(const size_type row, const size_type col) const
void vmult_add(Vector< somenumber > &dst, const Vector< somenumber > &src) const
void block_read(std::istream &in)
void precondition_TSOR(Vector< somenumber > &dst, const Vector< somenumber > &src, const number om=1.) const
bool is_finite(const double x)
void conjugate_add(const MATRIXA &A, const MATRIXB &B, const bool transpose=false)
static const size_type invalid
Iterator is invalid, probably due to an error.
DeclException0(ExcNoDiagonal)
Accessor(const SparseMatrixEZ< number > *matrix, const size_type row, const unsigned short index)
void precondition_SOR(Vector< somenumber > &dst, const Vector< somenumber > &src, const number om=1.) const
unsigned int global_dof_index
Definition: types.h:100
#define Assert(cond, exc)
Definition: exceptions.h:299
const Accessor * operator->() const
size_type n() const
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)
const_iterator(const SparseMatrixEZ< number > *matrix, const size_type row, const unsigned short index)
number el(const size_type i, const size_type j) const
const_iterator end() const
number operator()(const size_type i, const size_type j) const
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
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
::ExceptionBase & ExcIndexRange(int arg1, int arg2, int arg3)
const SparseMatrixEZ< number > * matrix
std::vector< Entry > data
Entry * allocate(const size_type row, const size_type col)
SparseMatrixEZ< number > & copy_from(const MATRIX &source)
void Tvmult(Vector< somenumber > &dst, const Vector< somenumber > &src) const
const_iterator begin() const
void set(const size_type i, const size_type j, const number value)
RowInfo(size_type start=Entry::invalid)
::ExceptionBase & ExcIteratorPastEnd()
types::global_dof_index size_type
size_type get_row_length(const size_type row) const
::ExceptionBase & ExcNumberNotFinite()
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
void print_statistics(STREAM &s, bool full=false)
::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
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
static const unsigned short invalid_diagonal
void Tvmult_add(Vector< somenumber > &dst, const Vector< somenumber > &src) const
void threaded_matrix_norm_square(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 block_write(std::ostream &out) const
size_type m() const
::ExceptionBase & ExcInternalError()
void print(std::ostream &out) const
bool operator==(const const_iterator &) const
std::vector< RowInfo > row_info
DeclException2(ExcInvalidEntry, int, int,<< "The entry with index ("<< arg1<< ','<< arg2<< ") does not exist.")