Reference documentation for deal.II version 8.1.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Modules Pages
block_matrix_base.h
1 // ---------------------------------------------------------------------
2 // @f$Id: block_matrix_base.h 31932 2013-12-08 02:15:54Z heister @f$
3 //
4 // Copyright (C) 2004 - 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__block_matrix_base_h
18 #define __deal2__block_matrix_base_h
19 
20 
21 #include <deal.II/base/config.h>
22 #include <deal.II/base/table.h>
23 #include <deal.II/base/thread_management.h>
24 #include <deal.II/base/utilities.h>
25 #include <deal.II/base/smartpointer.h>
26 #include <deal.II/base/memory_consumption.h>
27 #include <deal.II/lac/block_indices.h>
28 #include <deal.II/lac/exceptions.h>
29 #include <deal.II/lac/full_matrix.h>
30 #include <deal.II/lac/matrix_iterator.h>
31 #include <deal.II/lac/vector.h>
32 
33 #include <cmath>
34 
36 
37 
38 template <typename> class MatrixIterator;
39 
40 
41 
52 {
58  template <class BlockMatrix>
60  {
61  public:
66 
71  typedef typename BlockMatrix::value_type value_type;
72 
77  AccessorBase ();
78 
84  unsigned int block_row() const;
85 
91  unsigned int block_column() const;
92 
93  protected:
98  unsigned int row_block;
99 
104  unsigned int col_block;
105 
110  template <typename>
111  friend class MatrixIterator;
112  };
113 
114 
115 
120  template <class BlockMatrix, bool ConstNess>
121  class Accessor;
122 
123 
128  template <class BlockMatrix>
129  class Accessor<BlockMatrix, false>
130  :
131  public AccessorBase<BlockMatrix>
132  {
133  public:
138 
143  typedef BlockMatrix MatrixType;
144 
149  typedef typename BlockMatrix::value_type value_type;
150 
163  Accessor (BlockMatrix *m,
164  const size_type row,
165  const size_type col);
166 
172  size_type row() const;
173 
179  size_type column() const;
180 
185  value_type value() const;
186 
190  void set_value(value_type newval) const;
191 
192  protected:
196  BlockMatrix *matrix;
197 
202  typename BlockMatrix::BlockType::iterator base_iterator;
203 
207  void advance ();
208 
213  bool operator == (const Accessor &a) const;
214 
215  template <typename> friend class MatrixIterator;
216  friend class Accessor<BlockMatrix, true>;
217  };
218 
224  template <class BlockMatrix>
225  class Accessor<BlockMatrix, true>
226  :
227  public AccessorBase<BlockMatrix>
228  {
229  public:
234 
239  typedef const BlockMatrix MatrixType;
240 
245  typedef typename BlockMatrix::value_type value_type;
246 
259  Accessor (const BlockMatrix *m,
260  const size_type row,
261  const size_type col);
262 
268 
274  size_type row() const;
275 
281  size_type column() const;
282 
287  value_type value() const;
288  protected:
292  const BlockMatrix *matrix;
293 
298  typename BlockMatrix::BlockType::const_iterator base_iterator;
299 
303  void advance ();
304 
309  bool operator == (const Accessor &a) const;
310 
315  template <typename>
316  friend class ::MatrixIterator;
317  };
318 }
319 
320 
321 
381 template <typename MatrixType>
382 class BlockMatrixBase : public Subscriptor
383 {
384 public:
389  typedef MatrixType BlockType;
390 
396  typedef value_type *pointer;
397  typedef const value_type *const_pointer;
398  typedef value_type &reference;
399  typedef const value_type &const_reference;
400  typedef types::global_dof_index size_type;
401 
402  typedef
404  iterator;
405 
406  typedef
409 
410 
414  BlockMatrixBase ();
415 
419  ~BlockMatrixBase ();
420 
448  template <class BlockMatrixType>
450  copy_from (const BlockMatrixType &source);
451 
456  BlockType &
457  block (const unsigned int row,
458  const unsigned int column);
459 
460 
466  const BlockType &
467  block (const unsigned int row,
468  const unsigned int column) const;
469 
476  size_type m () const;
477 
484  size_type n () const;
485 
486 
493  unsigned int n_block_rows () const;
494 
501  unsigned int n_block_cols () const;
502 
512  void set (const size_type i,
513  const size_type j,
514  const value_type value);
515 
543  template <typename number>
544  void set (const std::vector<size_type> &indices,
545  const FullMatrix<number> &full_matrix,
546  const bool elide_zero_values = false);
547 
555  template <typename number>
556  void set (const std::vector<size_type> &row_indices,
557  const std::vector<size_type> &col_indices,
558  const FullMatrix<number> &full_matrix,
559  const bool elide_zero_values = false);
560 
579  template <typename number>
580  void set (const size_type row,
581  const std::vector<size_type> &col_indices,
582  const std::vector<number> &values,
583  const bool elide_zero_values = false);
584 
601  template <typename number>
602  void set (const size_type row,
603  const size_type n_cols,
604  const size_type *col_indices,
605  const number *values,
606  const bool elide_zero_values = false);
607 
617  void add (const size_type i,
618  const size_type j,
619  const value_type value);
620 
648  template <typename number>
649  void add (const std::vector<size_type> &indices,
650  const FullMatrix<number> &full_matrix,
651  const bool elide_zero_values = true);
652 
660  template <typename number>
661  void add (const std::vector<size_type> &row_indices,
662  const std::vector<size_type> &col_indices,
663  const FullMatrix<number> &full_matrix,
664  const bool elide_zero_values = true);
665 
683  template <typename number>
684  void add (const size_type row,
685  const std::vector<size_type> &col_indices,
686  const std::vector<number> &values,
687  const bool elide_zero_values = true);
688 
706  template <typename number>
707  void add (const size_type row,
708  const size_type n_cols,
709  const size_type *col_indices,
710  const number *values,
711  const bool elide_zero_values = true,
712  const bool col_indices_are_sorted = false);
713 
726  void add (const value_type factor,
727  const BlockMatrixBase<MatrixType> &matrix);
728 
740  value_type operator () (const size_type i,
741  const size_type j) const;
742 
760  value_type el (const size_type i,
761  const size_type j) const;
762 
780  value_type diag_element (const size_type i) const;
781 
791  void compress (::VectorOperation::values operation);
792 
797 
802  BlockMatrixBase &operator *= (const value_type factor);
803 
808  BlockMatrixBase &operator /= (const value_type factor);
809 
816  template <class BlockVectorType>
817  void vmult_add (BlockVectorType &dst,
818  const BlockVectorType &src) const;
819 
830  template <class BlockVectorType>
831  void Tvmult_add (BlockVectorType &dst,
832  const BlockVectorType &src) const;
833 
858  template <class BlockVectorType>
859  value_type
860  matrix_norm_square (const BlockVectorType &v) const;
861 
866  template <class BlockVectorType>
867  value_type
868  matrix_scalar_product (const BlockVectorType &u,
869  const BlockVectorType &v) const;
870 
876  template <class BlockVectorType>
877  value_type residual (BlockVectorType &dst,
878  const BlockVectorType &x,
879  const BlockVectorType &b) const;
880 
891  void print (std::ostream &out,
892  const bool alternative_output = false) const;
893 
898  iterator begin ();
899 
903  iterator end ();
904 
909  iterator begin (const size_type r);
910 
914  iterator end (const size_type r);
919  const_iterator begin () const;
920 
924  const_iterator end () const;
925 
930  const_iterator begin (const size_type r) const;
931 
935  const_iterator end (const size_type r) const;
936 
941  const BlockIndices &get_row_indices () const;
942 
947  const BlockIndices &get_column_indices () const;
948 
957  std::size_t memory_consumption () const;
958 
965  DeclException4 (ExcIncompatibleRowNumbers,
966  int, int, int, int,
967  << "The blocks [" << arg1 << ',' << arg2 << "] and ["
968  << arg3 << ',' << arg4 << "] have differing row numbers.");
972  DeclException4 (ExcIncompatibleColNumbers,
973  int, int, int, int,
974  << "The blocks [" << arg1 << ',' << arg2 << "] and ["
975  << arg3 << ',' << arg4 << "] have differing column numbers.");
977 protected:
997  void clear ();
998 
1003  BlockIndices column_block_indices;
1004 
1008  Table<2,SmartPointer<BlockType, BlockMatrixBase<MatrixType> > > sub_objects;
1009 
1045  void collect_sizes ();
1046 
1063  template <class BlockVectorType>
1064  void vmult_block_block (BlockVectorType &dst,
1065  const BlockVectorType &src) const;
1066 
1085  template <class BlockVectorType,
1086  class VectorType>
1087  void vmult_block_nonblock (BlockVectorType &dst,
1088  const VectorType &src) const;
1089 
1108  template <class BlockVectorType,
1109  class VectorType>
1110  void vmult_nonblock_block (VectorType &dst,
1111  const BlockVectorType &src) const;
1112 
1131  template <class VectorType>
1132  void vmult_nonblock_nonblock (VectorType &dst,
1133  const VectorType &src) const;
1134 
1154  template <class BlockVectorType>
1155  void Tvmult_block_block (BlockVectorType &dst,
1156  const BlockVectorType &src) const;
1157 
1176  template <class BlockVectorType,
1177  class VectorType>
1178  void Tvmult_block_nonblock (BlockVectorType &dst,
1179  const VectorType &src) const;
1180 
1199  template <class BlockVectorType,
1200  class VectorType>
1201  void Tvmult_nonblock_block (VectorType &dst,
1202  const BlockVectorType &src) const;
1203 
1222  template <class VectorType>
1223  void Tvmult_nonblock_nonblock (VectorType &dst,
1224  const VectorType &src) const;
1225 
1226 
1227 protected:
1228 
1239  void prepare_add_operation();
1240 
1246  void prepare_set_operation();
1247 
1248 
1249 private:
1250 
1262  {
1269  std::vector<size_type> counter_within_block;
1270 
1277  std::vector<std::vector<size_type> > column_indices;
1278 
1285  std::vector<std::vector<double> > column_values;
1286 
1292 
1304  TemporaryData &operator = (const TemporaryData &)
1305  {
1306  return *this;
1307  }
1308  };
1309 
1317  TemporaryData temporary_data;
1318 
1324  template <typename, bool>
1326 
1327  template <typename>
1328  friend class MatrixIterator;
1329 };
1330 
1331 
1334 #ifndef DOXYGEN
1335 /* ------------------------- Template functions ---------------------- */
1336 
1337 
1338 namespace BlockMatrixIterators
1339 {
1340  template <class BlockMatrix>
1341  inline
1343  :
1344  row_block(0),
1345  col_block(0)
1346  {}
1347 
1348 
1349  template <class BlockMatrix>
1350  inline
1351  unsigned int
1352  AccessorBase<BlockMatrix>::block_row() const
1353  {
1354  Assert (row_block != numbers::invalid_unsigned_int,
1355  ExcIteratorPastEnd());
1356 
1357  return row_block;
1358  }
1359 
1360 
1361  template <class BlockMatrix>
1362  inline
1363  unsigned int
1364  AccessorBase<BlockMatrix>::block_column() const
1365  {
1366  Assert (col_block != numbers::invalid_unsigned_int,
1367  ExcIteratorPastEnd());
1368 
1369  return col_block;
1370  }
1371 
1372 
1373  template <class BlockMatrix>
1374  inline
1375  Accessor<BlockMatrix, true>::Accessor (
1376  const BlockMatrix *matrix,
1377  const size_type row,
1378  const size_type col)
1379  :
1380  matrix(matrix),
1381  base_iterator(matrix->block(0,0).begin())
1382  {
1383  Assert(col==0, ExcNotImplemented());
1384 
1385  // check if this is a regular row or
1386  // the end of the matrix
1387  if (row < matrix->m())
1388  {
1389  const std::pair<unsigned int,size_type> indices
1390  = matrix->row_block_indices.global_to_local(row);
1391 
1392  // find the first block that does
1393  // have an entry in this row
1394  for (unsigned int bc=0; bc<matrix->n_block_cols(); ++bc)
1395  {
1396  base_iterator
1397  = matrix->block(indices.first, bc).begin(indices.second);
1398  if (base_iterator !=
1399  matrix->block(indices.first, bc).end(indices.second))
1400  {
1401  this->row_block = indices.first;
1402  this->col_block = bc;
1403  return;
1404  }
1405  }
1406 
1407  // hm, there is no block that has
1408  // an entry in this column. we need
1409  // to take the next entry then,
1410  // which may be the first entry of
1411  // the next row, or recursively the
1412  // next row, or so on
1413  *this = Accessor (matrix, row+1, 0);
1414  }
1415  else
1416  {
1417  // we were asked to create the end
1418  // iterator for this matrix
1419  this->row_block = numbers::invalid_unsigned_int;
1420  this->col_block = numbers::invalid_unsigned_int;
1421  }
1422  }
1423 
1424 
1425 // template <class BlockMatrix>
1426 // inline
1427 // Accessor<BlockMatrix, true>::Accessor (const Accessor<BlockMatrix, true>& other)
1428 // :
1429 // matrix(other.matrix),
1430 // base_iterator(other.base_iterator)
1431 // {
1432 // this->row_block = other.row_block;
1433 // this->col_block = other.col_block;
1434 // }
1435 
1436 
1437  template <class BlockMatrix>
1438  inline
1439  Accessor<BlockMatrix, true>::Accessor (const Accessor<BlockMatrix, false> &other)
1440  :
1441  matrix(other.matrix),
1442  base_iterator(other.base_iterator)
1443  {
1444  this->row_block = other.row_block;
1445  this->col_block = other.col_block;
1446  }
1447 
1448 
1449  template <class BlockMatrix>
1450  inline
1451  typename Accessor<BlockMatrix, true>::size_type
1452  Accessor<BlockMatrix, true>::row() const
1453  {
1454  Assert (this->row_block != numbers::invalid_unsigned_int,
1455  ExcIteratorPastEnd());
1456 
1457  return (matrix->row_block_indices.local_to_global(this->row_block, 0) +
1458  base_iterator->row());
1459  }
1460 
1461 
1462  template <class BlockMatrix>
1463  inline
1464  typename Accessor<BlockMatrix, true>::size_type
1465  Accessor<BlockMatrix, true>::column() const
1466  {
1467  Assert (this->col_block != numbers::invalid_unsigned_int,
1468  ExcIteratorPastEnd());
1469 
1470  return (matrix->column_block_indices.local_to_global(this->col_block,0) +
1471  base_iterator->column());
1472  }
1473 
1474 
1475  template <class BlockMatrix>
1476  inline
1477  typename Accessor<BlockMatrix, true>::value_type
1478  Accessor<BlockMatrix, true>::value () const
1479  {
1480  Assert (this->row_block != numbers::invalid_unsigned_int,
1481  ExcIteratorPastEnd());
1482  Assert (this->col_block != numbers::invalid_unsigned_int,
1483  ExcIteratorPastEnd());
1484 
1485  return base_iterator->value();
1486  }
1487 
1488 
1489 
1490  template <class BlockMatrix>
1491  inline
1492  void
1493  Accessor<BlockMatrix, true>::advance ()
1494  {
1495  Assert (this->row_block != numbers::invalid_unsigned_int,
1496  ExcIteratorPastEnd());
1497  Assert (this->col_block != numbers::invalid_unsigned_int,
1498  ExcIteratorPastEnd());
1499 
1500  // Remember current row inside block
1501  size_type local_row = base_iterator->row();
1502 
1503  // Advance one element inside the
1504  // current block
1505  ++base_iterator;
1506 
1507  // while we hit the end of the row of a
1508  // block (which may happen multiple
1509  // times if rows inside a block are
1510  // empty), we have to jump to the next
1511  // block and take the
1512  while (base_iterator ==
1513  matrix->block(this->row_block, this->col_block).end(local_row))
1514  {
1515  // jump to next block in this block
1516  // row, if possible, otherwise go
1517  // to next row
1518  if (this->col_block < matrix->n_block_cols()-1)
1519  {
1520  ++this->col_block;
1521  base_iterator
1522  = matrix->block(this->row_block, this->col_block).begin(local_row);
1523  }
1524  else
1525  {
1526  // jump back to next row in
1527  // first block column
1528  this->col_block = 0;
1529  ++local_row;
1530 
1531  // see if this has brought us
1532  // past the number of rows in
1533  // this block. if so see
1534  // whether we've just fallen
1535  // off the end of the whole
1536  // matrix
1537  if (local_row == matrix->block(this->row_block, this->col_block).m())
1538  {
1539  local_row = 0;
1540  ++this->row_block;
1541  if (this->row_block == matrix->n_block_rows())
1542  {
1543  this->row_block = numbers::invalid_unsigned_int;
1544  this->col_block = numbers::invalid_unsigned_int;
1545  return;
1546  }
1547  }
1548 
1549  base_iterator
1550  = matrix->block(this->row_block, this->col_block).begin(local_row);
1551  }
1552  }
1553  }
1554 
1555 
1556  template <class BlockMatrix>
1557  inline
1558  bool
1559  Accessor<BlockMatrix, true>::operator == (const Accessor &a) const
1560  {
1561  if (matrix != a.matrix)
1562  return false;
1563 
1564  if (this->row_block == a.row_block
1565  && this->col_block == a.col_block)
1566  // end iterators do not necessarily
1567  // have to have the same
1568  // base_iterator representation, but
1569  // valid iterators have to
1570  return (((this->row_block == numbers::invalid_unsigned_int)
1571  &&
1572  (this->col_block == numbers::invalid_unsigned_int))
1573  ||
1574  (base_iterator == a.base_iterator));
1575 
1576  return false;
1577  }
1578 
1579 //----------------------------------------------------------------------//
1580 
1581 
1582  template <class BlockMatrix>
1583  inline
1584  Accessor<BlockMatrix, false>::Accessor (
1585  BlockMatrix *matrix,
1586  const size_type row,
1587  const size_type col)
1588  :
1589  matrix(matrix),
1590  base_iterator(matrix->block(0,0).begin())
1591  {
1592  Assert(col==0, ExcNotImplemented());
1593  // check if this is a regular row or
1594  // the end of the matrix
1595  if (row < matrix->m())
1596  {
1597  const std::pair<unsigned int,size_type> indices
1598  = matrix->row_block_indices.global_to_local(row);
1599 
1600  // find the first block that does
1601  // have an entry in this row
1602  for (size_type bc=0; bc<matrix->n_block_cols(); ++bc)
1603  {
1604  base_iterator
1605  = matrix->block(indices.first, bc).begin(indices.second);
1606  if (base_iterator !=
1607  matrix->block(indices.first, bc).end(indices.second))
1608  {
1609  this->row_block = indices.first;
1610  this->col_block = bc;
1611  return;
1612  }
1613  }
1614 
1615  // hm, there is no block that has
1616  // an entry in this column. we need
1617  // to take the next entry then,
1618  // which may be the first entry of
1619  // the next row, or recursively the
1620  // next row, or so on
1621  *this = Accessor (matrix, row+1, 0);
1622  }
1623  else
1624  {
1625  // we were asked to create the end
1626  // iterator for this matrix
1627  this->row_block = numbers::invalid_size_type;
1628  this->col_block = numbers::invalid_size_type;
1629  }
1630  }
1631 
1632 
1633  template <class BlockMatrix>
1634  inline
1635  typename Accessor<BlockMatrix, false>::size_type
1636  Accessor<BlockMatrix, false>::row() const
1637  {
1638  Assert (this->row_block != numbers::invalid_size_type,
1639  ExcIteratorPastEnd());
1640 
1641  return (matrix->row_block_indices.local_to_global(this->row_block, 0) +
1642  base_iterator->row());
1643  }
1644 
1645 
1646  template <class BlockMatrix>
1647  inline
1648  typename Accessor<BlockMatrix, false>::size_type
1649  Accessor<BlockMatrix, false>::column() const
1650  {
1651  Assert (this->col_block != numbers::invalid_size_type,
1652  ExcIteratorPastEnd());
1653 
1654  return (matrix->column_block_indices.local_to_global(this->col_block,0) +
1655  base_iterator->column());
1656  }
1657 
1658 
1659  template <class BlockMatrix>
1660  inline
1661  typename Accessor<BlockMatrix, false>::value_type
1662  Accessor<BlockMatrix, false>::value () const
1663  {
1664  Assert (this->row_block != numbers::invalid_size_type,
1665  ExcIteratorPastEnd());
1666  Assert (this->col_block != numbers::invalid_size_type,
1667  ExcIteratorPastEnd());
1668 
1669  return base_iterator->value();
1670  }
1671 
1672 
1673 
1674  template <class BlockMatrix>
1675  inline
1676  void
1677  Accessor<BlockMatrix, false>::set_value (typename Accessor<BlockMatrix, false>::value_type newval) const
1678  {
1679  Assert (this->row_block != numbers::invalid_size_type,
1680  ExcIteratorPastEnd());
1681  Assert (this->col_block != numbers::invalid_size_type,
1682  ExcIteratorPastEnd());
1683 
1684  base_iterator->value() = newval;
1685  }
1686 
1687 
1688 
1689  template <class BlockMatrix>
1690  inline
1691  void
1692  Accessor<BlockMatrix, false>::advance ()
1693  {
1694  Assert (this->row_block != numbers::invalid_size_type,
1695  ExcIteratorPastEnd());
1696  Assert (this->col_block != numbers::invalid_size_type,
1697  ExcIteratorPastEnd());
1698 
1699  // Remember current row inside block
1700  size_type local_row = base_iterator->row();
1701 
1702  // Advance one element inside the
1703  // current block
1704  ++base_iterator;
1705 
1706  // while we hit the end of the row of a
1707  // block (which may happen multiple
1708  // times if rows inside a block are
1709  // empty), we have to jump to the next
1710  // block and take the
1711  while (base_iterator ==
1712  matrix->block(this->row_block, this->col_block).end(local_row))
1713  {
1714  // jump to next block in this block
1715  // row, if possible, otherwise go
1716  // to next row
1717  if (this->col_block < matrix->n_block_cols()-1)
1718  {
1719  ++this->col_block;
1720  base_iterator
1721  = matrix->block(this->row_block, this->col_block).begin(local_row);
1722  }
1723  else
1724  {
1725  // jump back to next row in
1726  // first block column
1727  this->col_block = 0;
1728  ++local_row;
1729 
1730  // see if this has brought us
1731  // past the number of rows in
1732  // this block. if so see
1733  // whether we've just fallen
1734  // off the end of the whole
1735  // matrix
1736  if (local_row == matrix->block(this->row_block, this->col_block).m())
1737  {
1738  local_row = 0;
1739  ++this->row_block;
1740  if (this->row_block == matrix->n_block_rows())
1741  {
1742  this->row_block = numbers::invalid_size_type;
1743  this->col_block = numbers::invalid_size_type;
1744  return;
1745  }
1746  }
1747 
1748  base_iterator
1749  = matrix->block(this->row_block, this->col_block).begin(local_row);
1750  }
1751  }
1752  }
1753 
1754 
1755 
1756  template <class BlockMatrix>
1757  inline
1758  bool
1759  Accessor<BlockMatrix, false>::operator == (const Accessor &a) const
1760  {
1761  if (matrix != a.matrix)
1762  return false;
1763 
1764  if (this->row_block == a.row_block
1765  && this->col_block == a.col_block)
1766  // end iterators do not necessarily
1767  // have to have the same
1768  // base_iterator representation, but
1769  // valid iterators have to
1770  return (((this->row_block == numbers::invalid_size_type)
1771  &&
1772  (this->col_block == numbers::invalid_size_type))
1773  ||
1774  (base_iterator == a.base_iterator));
1775 
1776  return false;
1777  }
1778 }
1779 
1780 
1781 //---------------------------------------------------------------------------
1782 
1783 
1784 template <typename MatrixType>
1785 inline
1787 {}
1788 
1789 template <typename MatrixType>
1790 inline
1792 {
1793  clear ();
1794 }
1795 
1796 
1797 template <class MatrixType>
1798 template <class BlockMatrixType>
1799 inline
1802 copy_from (const BlockMatrixType &source)
1803 {
1804  for (unsigned int r=0; r<n_block_rows(); ++r)
1805  for (unsigned int c=0; c<n_block_cols(); ++c)
1806  block(r,c).copy_from (source.block(r,c));
1807 
1808  return *this;
1809 }
1810 
1811 
1812 template <class MatrixType>
1813 std::size_t
1815 {
1816  std::size_t mem =
1817  MemoryConsumption::memory_consumption(row_block_indices)+
1818  MemoryConsumption::memory_consumption(column_block_indices)+
1820  MemoryConsumption::memory_consumption(temporary_data.counter_within_block)+
1821  MemoryConsumption::memory_consumption(temporary_data.column_indices)+
1822  MemoryConsumption::memory_consumption(temporary_data.column_values)+
1823  sizeof(temporary_data.mutex);
1824 
1825  for (unsigned int r=0; r<n_block_rows(); ++r)
1826  for (unsigned int c=0; c<n_block_cols(); ++c)
1827  {
1828  MatrixType *p = this->sub_objects[r][c];
1830  }
1831 
1832  return mem;
1833 }
1834 
1835 
1836 
1837 template <class MatrixType>
1838 inline
1839 void
1841 {
1842  for (unsigned int r=0; r<n_block_rows(); ++r)
1843  for (unsigned int c=0; c<n_block_cols(); ++c)
1844  {
1845  MatrixType *p = this->sub_objects[r][c];
1846  this->sub_objects[r][c] = 0;
1847  delete p;
1848  }
1849  sub_objects.reinit (0,0);
1850 
1851  // reset block indices to empty
1852  row_block_indices = column_block_indices = BlockIndices ();
1853 }
1854 
1855 
1856 
1857 template <class MatrixType>
1858 inline
1860 BlockMatrixBase<MatrixType>::block (const unsigned int row,
1861  const unsigned int column)
1862 {
1863  Assert (row<n_block_rows(),
1864  ExcIndexRange (row, 0, n_block_rows()));
1865  Assert (column<n_block_cols(),
1866  ExcIndexRange (column, 0, n_block_cols()));
1867 
1868  return *sub_objects[row][column];
1869 }
1870 
1871 
1872 
1873 template <class MatrixType>
1874 inline
1876 BlockMatrixBase<MatrixType>::block (const unsigned int row,
1877  const unsigned int column) const
1878 {
1879  Assert (row<n_block_rows(),
1880  ExcIndexRange (row, 0, n_block_rows()));
1881  Assert (column<n_block_cols(),
1882  ExcIndexRange (column, 0, n_block_cols()));
1883 
1884  return *sub_objects[row][column];
1885 }
1886 
1887 
1888 template <class MatrixType>
1889 inline
1890 typename BlockMatrixBase<MatrixType>::size_type
1892 {
1893  return row_block_indices.total_size();
1894 }
1895 
1896 
1897 
1898 template <class MatrixType>
1899 inline
1900 typename BlockMatrixBase<MatrixType>::size_type
1902 {
1903  return column_block_indices.total_size();
1904 }
1905 
1906 
1907 
1908 template <class MatrixType>
1909 inline
1910 unsigned int
1912 {
1913  return column_block_indices.size();
1914 }
1915 
1916 
1917 
1918 template <class MatrixType>
1919 inline
1920 unsigned int
1922 {
1923  return row_block_indices.size();
1924 }
1925 
1926 
1927 
1928 // Write the single set manually,
1929 // since the other function has a lot
1930 // of overhead in that case.
1931 template <class MatrixType>
1932 inline
1933 void
1935  const size_type j,
1936  const value_type value)
1937 {
1938  prepare_set_operation();
1939 
1941 
1942  const std::pair<unsigned int,size_type>
1943  row_index = row_block_indices.global_to_local (i),
1944  col_index = column_block_indices.global_to_local (j);
1945  block(row_index.first,col_index.first).set (row_index.second,
1946  col_index.second,
1947  value);
1948 }
1949 
1950 
1951 
1952 template <class MatrixType>
1953 template <typename number>
1954 inline
1955 void
1956 BlockMatrixBase<MatrixType>::set (const std::vector<size_type> &row_indices,
1957  const std::vector<size_type> &col_indices,
1958  const FullMatrix<number> &values,
1959  const bool elide_zero_values)
1960 {
1961  Assert (row_indices.size() == values.m(),
1962  ExcDimensionMismatch(row_indices.size(), values.m()));
1963  Assert (col_indices.size() == values.n(),
1964  ExcDimensionMismatch(col_indices.size(), values.n()));
1965 
1966  for (size_type i=0; i<row_indices.size(); ++i)
1967  set (row_indices[i], col_indices.size(), &col_indices[0], &values(i,0),
1968  elide_zero_values);
1969 }
1970 
1971 
1972 
1973 template <class MatrixType>
1974 template <typename number>
1975 inline
1976 void
1977 BlockMatrixBase<MatrixType>::set (const std::vector<size_type> &indices,
1978  const FullMatrix<number> &values,
1979  const bool elide_zero_values)
1980 {
1981  Assert (indices.size() == values.m(),
1982  ExcDimensionMismatch(indices.size(), values.m()));
1983  Assert (values.n() == values.m(), ExcNotQuadratic());
1984 
1985  for (size_type i=0; i<indices.size(); ++i)
1986  set (indices[i], indices.size(), &indices[0], &values(i,0),
1987  elide_zero_values);
1988 }
1989 
1990 
1991 
1992 template <class MatrixType>
1993 template <typename number>
1994 inline
1995 void
1997  const std::vector<size_type> &col_indices,
1998  const std::vector<number> &values,
1999  const bool elide_zero_values)
2000 {
2001  Assert (col_indices.size() == values.size(),
2002  ExcDimensionMismatch(col_indices.size(), values.size()));
2003 
2004  set (row, col_indices.size(), &col_indices[0], &values[0],
2005  elide_zero_values);
2006 }
2007 
2008 
2009 
2010 // This is a very messy function, since
2011 // we need to calculate to each position
2012 // the location in the global array.
2013 template <class MatrixType>
2014 template <typename number>
2015 inline
2016 void
2018  const size_type n_cols,
2019  const size_type *col_indices,
2020  const number *values,
2021  const bool elide_zero_values)
2022 {
2023  prepare_set_operation();
2024 
2025  // lock access to the temporary data structure to
2026  // allow multiple threads to call this function concurrently
2027  Threads::Mutex::ScopedLock lock (temporary_data.mutex);
2028 
2029  // Resize scratch arrays
2030  if (temporary_data.column_indices.size() < this->n_block_cols())
2031  {
2032  temporary_data.column_indices.resize (this->n_block_cols());
2033  temporary_data.column_values.resize (this->n_block_cols());
2034  temporary_data.counter_within_block.resize (this->n_block_cols());
2035  }
2036 
2037  // Resize sub-arrays to n_cols. This
2038  // is a bit wasteful, but we resize
2039  // only a few times (then the maximum
2040  // row length won't increase that
2041  // much any more). At least we know
2042  // that all arrays are going to be of
2043  // the same size, so we can check
2044  // whether the size of one is large
2045  // enough before actually going
2046  // through all of them.
2047  if (temporary_data.column_indices[0].size() < n_cols)
2048  {
2049  for (unsigned int i=0; i<this->n_block_cols(); ++i)
2050  {
2051  temporary_data.column_indices[i].resize(n_cols);
2052  temporary_data.column_values[i].resize(n_cols);
2053  }
2054  }
2055 
2056  // Reset the number of added elements
2057  // in each block to zero.
2058  for (unsigned int i=0; i<this->n_block_cols(); ++i)
2059  temporary_data.counter_within_block[i] = 0;
2060 
2061  // Go through the column indices to
2062  // find out which portions of the
2063  // values should be set in which
2064  // block of the matrix. We need to
2065  // touch all the data, since we can't
2066  // be sure that the data of one block
2067  // is stored contiguously (in fact,
2068  // indices will be intermixed when it
2069  // comes from an element matrix).
2070  for (size_type j=0; j<n_cols; ++j)
2071  {
2072  double value = values[j];
2073 
2074  if (value == 0 && elide_zero_values == true)
2075  continue;
2076 
2077  const std::pair<unsigned int, size_type>
2078  col_index = this->column_block_indices.global_to_local(col_indices[j]);
2079 
2080  const size_type local_index = temporary_data.counter_within_block[col_index.first]++;
2081 
2082  temporary_data.column_indices[col_index.first][local_index] = col_index.second;
2083  temporary_data.column_values[col_index.first][local_index] = value;
2084  }
2085 
2086 #ifdef DEBUG
2087  // If in debug mode, do a check whether
2088  // the right length has been obtained.
2089  size_type length = 0;
2090  for (unsigned int i=0; i<this->n_block_cols(); ++i)
2091  length += temporary_data.counter_within_block[i];
2092  Assert (length <= n_cols, ExcInternalError());
2093 #endif
2094 
2095  // Now we found out about where the
2096  // individual columns should start and
2097  // where we should start reading out
2098  // data. Now let's write the data into
2099  // the individual blocks!
2100  const std::pair<unsigned int,size_type>
2101  row_index = this->row_block_indices.global_to_local (row);
2102  for (unsigned int block_col=0; block_col<n_block_cols(); ++block_col)
2103  {
2104  if (temporary_data.counter_within_block[block_col] == 0)
2105  continue;
2106 
2107  block(row_index.first, block_col).set
2108  (row_index.second,
2109  temporary_data.counter_within_block[block_col],
2110  &temporary_data.column_indices[block_col][0],
2111  &temporary_data.column_values[block_col][0],
2112  false);
2113  }
2114 }
2115 
2116 
2117 
2118 template <class MatrixType>
2119 inline
2120 void
2122  const size_type j,
2123  const value_type value)
2124 {
2125 
2127 
2128  prepare_add_operation();
2129 
2130  // save some cycles for zero additions, but
2131  // only if it is safe for the matrix we are
2132  // working with
2133  typedef typename MatrixType::Traits MatrixTraits;
2134  if ((MatrixTraits::zero_addition_can_be_elided == true)
2135  &&
2136  (value == value_type()))
2137  return;
2138 
2139  const std::pair<unsigned int,size_type>
2140  row_index = row_block_indices.global_to_local (i),
2141  col_index = column_block_indices.global_to_local (j);
2142  block(row_index.first,col_index.first).add (row_index.second,
2143  col_index.second,
2144  value);
2145 }
2146 
2147 
2148 
2149 template <class MatrixType>
2150 template <typename number>
2151 inline
2152 void
2153 BlockMatrixBase<MatrixType>::add (const std::vector<size_type> &row_indices,
2154  const std::vector<size_type> &col_indices,
2155  const FullMatrix<number> &values,
2156  const bool elide_zero_values)
2157 {
2158  Assert (row_indices.size() == values.m(),
2159  ExcDimensionMismatch(row_indices.size(), values.m()));
2160  Assert (col_indices.size() == values.n(),
2161  ExcDimensionMismatch(col_indices.size(), values.n()));
2162 
2163  for (size_type i=0; i<row_indices.size(); ++i)
2164  add (row_indices[i], col_indices.size(), &col_indices[0], &values(i,0),
2165  elide_zero_values);
2166 }
2167 
2168 
2169 
2170 template <class MatrixType>
2171 template <typename number>
2172 inline
2173 void
2174 BlockMatrixBase<MatrixType>::add (const std::vector<size_type> &indices,
2175  const FullMatrix<number> &values,
2176  const bool elide_zero_values)
2177 {
2178  Assert (indices.size() == values.m(),
2179  ExcDimensionMismatch(indices.size(), values.m()));
2180  Assert (values.n() == values.m(), ExcNotQuadratic());
2181 
2182  for (size_type i=0; i<indices.size(); ++i)
2183  add (indices[i], indices.size(), &indices[0], &values(i,0),
2184  elide_zero_values);
2185 }
2186 
2187 
2188 
2189 template <class MatrixType>
2190 template <typename number>
2191 inline
2192 void
2194  const std::vector<size_type> &col_indices,
2195  const std::vector<number> &values,
2196  const bool elide_zero_values)
2197 {
2198  Assert (col_indices.size() == values.size(),
2199  ExcDimensionMismatch(col_indices.size(), values.size()));
2200 
2201  add (row, col_indices.size(), &col_indices[0], &values[0],
2202  elide_zero_values);
2203 }
2204 
2205 
2206 
2207 // This is a very messy function, since
2208 // we need to calculate to each position
2209 // the location in the global array.
2210 template <class MatrixType>
2211 template <typename number>
2212 inline
2213 void
2215  const size_type n_cols,
2216  const size_type *col_indices,
2217  const number *values,
2218  const bool elide_zero_values,
2219  const bool col_indices_are_sorted)
2220 {
2221  prepare_add_operation();
2222 
2223  // TODO: Look over this to find out
2224  // whether we can do that more
2225  // efficiently.
2226  if (col_indices_are_sorted == true)
2227  {
2228 #ifdef DEBUG
2229  // check whether indices really are
2230  // sorted.
2231  size_type before = col_indices[0];
2232  for (size_type i=1; i<n_cols; ++i)
2233  if (col_indices[i] <= before)
2234  Assert (false, ExcMessage ("Flag col_indices_are_sorted is set, but "
2235  "indices appear to not be sorted."))
2236  else
2237  before = col_indices[i];
2238 #endif
2239  const std::pair<unsigned int,size_type>
2240  row_index = this->row_block_indices.global_to_local (row);
2241 
2242  if (this->n_block_cols() > 1)
2243  {
2244  const size_type *first_block = Utilities::lower_bound (col_indices,
2245  col_indices+n_cols,
2246  this->column_block_indices.block_start(1));
2247 
2248  const size_type n_zero_block_indices = first_block - col_indices;
2249  block(row_index.first, 0).add (row_index.second,
2250  n_zero_block_indices,
2251  col_indices,
2252  values,
2253  elide_zero_values,
2254  col_indices_are_sorted);
2255 
2256  if (n_zero_block_indices < n_cols)
2257  this->add(row, n_cols - n_zero_block_indices, first_block,
2258  values + n_zero_block_indices, elide_zero_values,
2259  false);
2260  }
2261  else
2262  {
2263  block(row_index.first, 0). add (row_index.second,
2264  n_cols,
2265  col_indices,
2266  values,
2267  elide_zero_values,
2268  col_indices_are_sorted);
2269  }
2270 
2271  return;
2272  }
2273 
2274  // Lock scratch arrays, then resize them
2275  Threads::Mutex::ScopedLock lock (temporary_data.mutex);
2276 
2277  if (temporary_data.column_indices.size() < this->n_block_cols())
2278  {
2279  temporary_data.column_indices.resize (this->n_block_cols());
2280  temporary_data.column_values.resize (this->n_block_cols());
2281  temporary_data.counter_within_block.resize (this->n_block_cols());
2282  }
2283 
2284  // Resize sub-arrays to n_cols. This
2285  // is a bit wasteful, but we resize
2286  // only a few times (then the maximum
2287  // row length won't increase that
2288  // much any more). At least we know
2289  // that all arrays are going to be of
2290  // the same size, so we can check
2291  // whether the size of one is large
2292  // enough before actually going
2293  // through all of them.
2294  if (temporary_data.column_indices[0].size() < n_cols)
2295  {
2296  for (unsigned int i=0; i<this->n_block_cols(); ++i)
2297  {
2298  temporary_data.column_indices[i].resize(n_cols);
2299  temporary_data.column_values[i].resize(n_cols);
2300  }
2301  }
2302 
2303  // Reset the number of added elements
2304  // in each block to zero.
2305  for (unsigned int i=0; i<this->n_block_cols(); ++i)
2306  temporary_data.counter_within_block[i] = 0;
2307 
2308  // Go through the column indices to
2309  // find out which portions of the
2310  // values should be written into
2311  // which block of the matrix. We need
2312  // to touch all the data, since we
2313  // can't be sure that the data of one
2314  // block is stored contiguously (in
2315  // fact, data will be intermixed when
2316  // it comes from an element matrix).
2317  for (size_type j=0; j<n_cols; ++j)
2318  {
2319  double value = values[j];
2320 
2321  if (value == 0 && elide_zero_values == true)
2322  continue;
2323 
2324  const std::pair<unsigned int, size_type>
2325  col_index = this->column_block_indices.global_to_local(col_indices[j]);
2326 
2327  const size_type local_index = temporary_data.counter_within_block[col_index.first]++;
2328 
2329  temporary_data.column_indices[col_index.first][local_index] = col_index.second;
2330  temporary_data.column_values[col_index.first][local_index] = value;
2331  }
2332 
2333 #ifdef DEBUG
2334  // If in debug mode, do a check whether
2335  // the right length has been obtained.
2336  size_type length = 0;
2337  for (unsigned int i=0; i<this->n_block_cols(); ++i)
2338  length += temporary_data.counter_within_block[i];
2339  Assert (length <= n_cols, ExcInternalError());
2340 #endif
2341 
2342  // Now we found out about where the
2343  // individual columns should start and
2344  // where we should start reading out
2345  // data. Now let's write the data into
2346  // the individual blocks!
2347  const std::pair<unsigned int,size_type>
2348  row_index = this->row_block_indices.global_to_local (row);
2349  for (unsigned int block_col=0; block_col<n_block_cols(); ++block_col)
2350  {
2351  if (temporary_data.counter_within_block[block_col] == 0)
2352  continue;
2353 
2354  block(row_index.first, block_col).add
2355  (row_index.second,
2356  temporary_data.counter_within_block[block_col],
2357  &temporary_data.column_indices[block_col][0],
2358  &temporary_data.column_values[block_col][0],
2359  false,
2360  col_indices_are_sorted);
2361  }
2362 }
2363 
2364 
2365 
2366 template <class MatrixType>
2367 inline
2368 void
2369 BlockMatrixBase<MatrixType>::add (const value_type factor,
2370  const BlockMatrixBase<MatrixType> &matrix)
2371 {
2373 
2374  prepare_add_operation();
2375 
2376  // save some cycles for zero additions, but
2377  // only if it is safe for the matrix we are
2378  // working with
2379  typedef typename MatrixType::Traits MatrixTraits;
2380  if ((MatrixTraits::zero_addition_can_be_elided == true)
2381  &&
2382  (factor == 0))
2383  return;
2384 
2385  for (unsigned int row=0; row<n_block_rows(); ++row)
2386  for (unsigned int col=0; col<n_block_cols(); ++col)
2387  // This function should throw if the sparsity
2388  // patterns of the two blocks differ
2389  block(row, col).add(factor, matrix.block(row,col));
2390 }
2391 
2392 
2393 
2394 template <class MatrixType>
2395 inline
2398  const size_type j) const
2399 {
2400  const std::pair<unsigned int,size_type>
2401  row_index = row_block_indices.global_to_local (i),
2402  col_index = column_block_indices.global_to_local (j);
2403  return block(row_index.first,col_index.first) (row_index.second,
2404  col_index.second);
2405 }
2406 
2407 
2408 
2409 template <class MatrixType>
2410 inline
2413  const size_type j) const
2414 {
2415  const std::pair<unsigned int,size_type>
2416  row_index = row_block_indices.global_to_local (i),
2417  col_index = column_block_indices.global_to_local (j);
2418  return block(row_index.first,col_index.first).el (row_index.second,
2419  col_index.second);
2420 }
2421 
2422 
2423 
2424 template <class MatrixType>
2425 inline
2428 {
2429  Assert (n_block_rows() == n_block_cols(),
2430  ExcNotQuadratic());
2431 
2432  const std::pair<unsigned int,size_type>
2433  index = row_block_indices.global_to_local (i);
2434  return block(index.first,index.first).diag_element(index.second);
2435 }
2436 
2437 
2438 
2439 template <class MatrixType>
2440 inline
2441 void
2442 BlockMatrixBase<MatrixType>::compress (::VectorOperation::values operation)
2443 {
2444  for (unsigned int r=0; r<n_block_rows(); ++r)
2445  for (unsigned int c=0; c<n_block_cols(); ++c)
2446  block(r,c).compress (operation);
2447 }
2448 
2449 template <class MatrixType>
2450 inline
2451 void
2453 {
2454  compress(::VectorOperation::unknown);
2455 }
2456 
2457 
2458 
2459 template <class MatrixType>
2460 inline
2462 BlockMatrixBase<MatrixType>::operator *= (const value_type factor)
2463 {
2464  Assert (n_block_cols() != 0, ExcNotInitialized());
2465  Assert (n_block_rows() != 0, ExcNotInitialized());
2466 
2467  for (unsigned int r=0; r<n_block_rows(); ++r)
2468  for (unsigned int c=0; c<n_block_cols(); ++c)
2469  block(r,c) *= factor;
2470 
2471  return *this;
2472 }
2473 
2474 
2475 
2476 template <class MatrixType>
2477 inline
2479 BlockMatrixBase<MatrixType>::operator /= (const value_type factor)
2480 {
2481  Assert (n_block_cols() != 0, ExcNotInitialized());
2482  Assert (n_block_rows() != 0, ExcNotInitialized());
2483  Assert (factor !=0, ExcDivideByZero());
2484 
2485  const value_type factor_inv = 1. / factor;
2486 
2487  for (unsigned int r=0; r<n_block_rows(); ++r)
2488  for (unsigned int c=0; c<n_block_cols(); ++c)
2489  block(r,c) *= factor_inv;
2490 
2491  return *this;
2492 }
2493 
2494 
2495 
2496 template <class MatrixType>
2497 const BlockIndices &
2499 {
2500  return this->row_block_indices;
2501 }
2502 
2503 
2504 
2505 template <class MatrixType>
2506 const BlockIndices &
2508 {
2509  return this->column_block_indices;
2510 }
2511 
2512 
2513 
2514 template <class MatrixType>
2515 template <class BlockVectorType>
2516 void
2518 vmult_block_block (BlockVectorType &dst,
2519  const BlockVectorType &src) const
2520 {
2521  Assert (dst.n_blocks() == n_block_rows(),
2522  ExcDimensionMismatch(dst.n_blocks(), n_block_rows()));
2523  Assert (src.n_blocks() == n_block_cols(),
2524  ExcDimensionMismatch(src.n_blocks(), n_block_cols()));
2525 
2526  for (size_type row=0; row<n_block_rows(); ++row)
2527  {
2528  block(row,0).vmult (dst.block(row),
2529  src.block(0));
2530  for (size_type col=1; col<n_block_cols(); ++col)
2531  block(row,col).vmult_add (dst.block(row),
2532  src.block(col));
2533  };
2534 }
2535 
2536 
2537 
2538 template <class MatrixType>
2539 template <class BlockVectorType,
2540  class VectorType>
2541 void
2543 vmult_nonblock_block (VectorType &dst,
2544  const BlockVectorType &src) const
2545 {
2546  Assert (n_block_rows() == 1,
2547  ExcDimensionMismatch(1, n_block_rows()));
2548  Assert (src.n_blocks() == n_block_cols(),
2549  ExcDimensionMismatch(src.n_blocks(), n_block_cols()));
2550 
2551  block(0,0).vmult (dst, src.block(0));
2552  for (size_type col=1; col<n_block_cols(); ++col)
2553  block(0,col).vmult_add (dst, src.block(col));
2554 }
2555 
2556 
2557 
2558 template <class MatrixType>
2559 template <class BlockVectorType,
2560  class VectorType>
2561 void
2563 vmult_block_nonblock (BlockVectorType &dst,
2564  const VectorType &src) const
2565 {
2566  Assert (dst.n_blocks() == n_block_rows(),
2567  ExcDimensionMismatch(dst.n_blocks(), n_block_rows()));
2568  Assert (1 == n_block_cols(),
2569  ExcDimensionMismatch(1, n_block_cols()));
2570 
2571  for (size_type row=0; row<n_block_rows(); ++row)
2572  block(row,0).vmult (dst.block(row),
2573  src);
2574 }
2575 
2576 
2577 
2578 template <class MatrixType>
2579 template <class VectorType>
2580 void
2582 vmult_nonblock_nonblock (VectorType &dst,
2583  const VectorType &src) const
2584 {
2585  Assert (1 == n_block_rows(),
2586  ExcDimensionMismatch(1, n_block_rows()));
2587  Assert (1 == n_block_cols(),
2588  ExcDimensionMismatch(1, n_block_cols()));
2589 
2590  block(0,0).vmult (dst, src);
2591 }
2592 
2593 
2594 
2595 template <class MatrixType>
2596 template <class BlockVectorType>
2597 void
2598 BlockMatrixBase<MatrixType>::vmult_add (BlockVectorType &dst,
2599  const BlockVectorType &src) const
2600 {
2601  Assert (dst.n_blocks() == n_block_rows(),
2602  ExcDimensionMismatch(dst.n_blocks(), n_block_rows()));
2603  Assert (src.n_blocks() == n_block_cols(),
2604  ExcDimensionMismatch(src.n_blocks(), n_block_cols()));
2605 
2606  for (unsigned int row=0; row<n_block_rows(); ++row)
2607  for (unsigned int col=0; col<n_block_cols(); ++col)
2608  block(row,col).vmult_add (dst.block(row),
2609  src.block(col));
2610 }
2611 
2612 
2613 
2614 
2615 template <class MatrixType>
2616 template <class BlockVectorType>
2617 void
2619 Tvmult_block_block (BlockVectorType &dst,
2620  const BlockVectorType &src) const
2621 {
2622  Assert (dst.n_blocks() == n_block_cols(),
2623  ExcDimensionMismatch(dst.n_blocks(), n_block_cols()));
2624  Assert (src.n_blocks() == n_block_rows(),
2625  ExcDimensionMismatch(src.n_blocks(), n_block_rows()));
2626 
2627  dst = 0.;
2628 
2629  for (unsigned int row=0; row<n_block_rows(); ++row)
2630  {
2631  for (unsigned int col=0; col<n_block_cols(); ++col)
2632  block(row,col).Tvmult_add (dst.block(col),
2633  src.block(row));
2634  };
2635 }
2636 
2637 
2638 
2639 template <class MatrixType>
2640 template <class BlockVectorType,
2641  class VectorType>
2642 void
2644 Tvmult_block_nonblock (BlockVectorType &dst,
2645  const VectorType &src) const
2646 {
2647  Assert (dst.n_blocks() == n_block_cols(),
2648  ExcDimensionMismatch(dst.n_blocks(), n_block_cols()));
2649  Assert (1 == n_block_rows(),
2650  ExcDimensionMismatch(1, n_block_rows()));
2651 
2652  dst = 0.;
2653 
2654  for (unsigned int col=0; col<n_block_cols(); ++col)
2655  block(0,col).Tvmult_add (dst.block(col), src);
2656 }
2657 
2658 
2659 
2660 template <class MatrixType>
2661 template <class BlockVectorType,
2662  class VectorType>
2663 void
2665 Tvmult_nonblock_block (VectorType &dst,
2666  const BlockVectorType &src) const
2667 {
2668  Assert (1 == n_block_cols(),
2669  ExcDimensionMismatch(1, n_block_cols()));
2670  Assert (src.n_blocks() == n_block_rows(),
2671  ExcDimensionMismatch(src.n_blocks(), n_block_rows()));
2672 
2673  block(0,0).Tvmult (dst, src.block(0));
2674 
2675  for (size_type row=1; row<n_block_rows(); ++row)
2676  block(row,0).Tvmult_add (dst, src.block(row));
2677 }
2678 
2679 
2680 
2681 template <class MatrixType>
2682 template <class VectorType>
2683 void
2685 Tvmult_nonblock_nonblock (VectorType &dst,
2686  const VectorType &src) const
2687 {
2688  Assert (1 == n_block_cols(),
2689  ExcDimensionMismatch(1, n_block_cols()));
2690  Assert (1 == n_block_rows(),
2691  ExcDimensionMismatch(1, n_block_rows()));
2692 
2693  block(0,0).Tvmult (dst, src);
2694 }
2695 
2696 
2697 
2698 template <class MatrixType>
2699 template <class BlockVectorType>
2700 void
2701 BlockMatrixBase<MatrixType>::Tvmult_add (BlockVectorType &dst,
2702  const BlockVectorType &src) const
2703 {
2704  Assert (dst.n_blocks() == n_block_cols(),
2705  ExcDimensionMismatch(dst.n_blocks(), n_block_cols()));
2706  Assert (src.n_blocks() == n_block_rows(),
2707  ExcDimensionMismatch(src.n_blocks(), n_block_rows()));
2708 
2709  for (unsigned int row=0; row<n_block_rows(); ++row)
2710  for (unsigned int col=0; col<n_block_cols(); ++col)
2711  block(row,col).Tvmult_add (dst.block(col),
2712  src.block(row));
2713 }
2714 
2715 
2716 
2717 template <class MatrixType>
2718 template <class BlockVectorType>
2720 BlockMatrixBase<MatrixType>::matrix_norm_square (const BlockVectorType &v) const
2721 {
2722  Assert (n_block_rows() == n_block_cols(), ExcNotQuadratic());
2723  Assert (v.n_blocks() == n_block_rows(),
2724  ExcDimensionMismatch(v.n_blocks(), n_block_rows()));
2725 
2726  value_type norm_sqr = 0;
2727  for (unsigned int row=0; row<n_block_rows(); ++row)
2728  for (unsigned int col=0; col<n_block_cols(); ++col)
2729  if (row==col)
2730  norm_sqr += block(row,col).matrix_norm_square (v.block(row));
2731  else
2732  norm_sqr += block(row,col).matrix_scalar_product (v.block(row),
2733  v.block(col));
2734  return norm_sqr;
2735 }
2736 
2737 
2738 
2739 template <class MatrixType>
2740 template <class BlockVectorType>
2743 matrix_scalar_product (const BlockVectorType &u,
2744  const BlockVectorType &v) const
2745 {
2746  Assert (u.n_blocks() == n_block_rows(),
2747  ExcDimensionMismatch(u.n_blocks(), n_block_rows()));
2748  Assert (v.n_blocks() == n_block_cols(),
2749  ExcDimensionMismatch(v.n_blocks(), n_block_cols()));
2750 
2751  value_type result = 0;
2752  for (unsigned int row=0; row<n_block_rows(); ++row)
2753  for (unsigned int col=0; col<n_block_cols(); ++col)
2754  result += block(row,col).matrix_scalar_product (u.block(row),
2755  v.block(col));
2756  return result;
2757 }
2758 
2759 
2760 
2761 template <class MatrixType>
2762 template <class BlockVectorType>
2765 residual (BlockVectorType &dst,
2766  const BlockVectorType &x,
2767  const BlockVectorType &b) const
2768 {
2769  Assert (dst.n_blocks() == n_block_rows(),
2770  ExcDimensionMismatch(dst.n_blocks(), n_block_rows()));
2771  Assert (b.n_blocks() == n_block_rows(),
2772  ExcDimensionMismatch(b.n_blocks(), n_block_rows()));
2773  Assert (x.n_blocks() == n_block_cols(),
2774  ExcDimensionMismatch(x.n_blocks(), n_block_cols()));
2775  // in block notation, the residual is
2776  // r_i = b_i - \sum_j A_ij x_j.
2777  // this can be written as
2778  // r_i = b_i - A_i0 x_0 - \sum_{j>0} A_ij x_j.
2779  //
2780  // for the first two terms, we can
2781  // call the residual function of
2782  // A_i0. for the other terms, we
2783  // use vmult_add. however, we want
2784  // to subtract, so in order to
2785  // avoid a temporary vector, we
2786  // perform a sign change of the
2787  // first two term before, and after
2788  // adding up
2789  for (unsigned int row=0; row<n_block_rows(); ++row)
2790  {
2791  block(row,0).residual (dst.block(row),
2792  x.block(0),
2793  b.block(row));
2794 
2795  for (size_type i=0; i<dst.block(row).size(); ++i)
2796  dst.block(row)(i) = -dst.block(row)(i);
2797 
2798  for (unsigned int col=1; col<n_block_cols(); ++col)
2799  block(row,col).vmult_add (dst.block(row),
2800  x.block(col));
2801 
2802  for (size_type i=0; i<dst.block(row).size(); ++i)
2803  dst.block(row)(i) = -dst.block(row)(i);
2804  };
2805 
2806  value_type res = 0;
2807  for (size_type row=0; row<n_block_rows(); ++row)
2808  res += dst.block(row).norm_sqr ();
2809  return std::sqrt(res);
2810 }
2811 
2812 
2813 
2814 template <class MatrixType>
2815 inline
2816 void
2817 BlockMatrixBase<MatrixType>::print (std::ostream &out,
2818  const bool alternative_output) const
2819 {
2820  for (unsigned int row=0; row<n_block_rows(); ++row)
2821  for (unsigned int col=0; col<n_block_cols(); ++col)
2822  {
2823  if (!alternative_output)
2824  out << "Block (" << row << ", " << col << ")" << std::endl;
2825 
2826  block(row, col).print(out, alternative_output);
2827  }
2828 }
2829 
2830 
2831 
2832 template <class MatrixType>
2833 inline
2836 {
2837  return const_iterator(this, 0);
2838 }
2839 
2840 
2841 
2842 template <class MatrixType>
2843 inline
2846 {
2847  return const_iterator(this, m());
2848 }
2849 
2850 
2851 
2852 template <class MatrixType>
2853 inline
2856 {
2857  Assert (r<m(), ExcIndexRange(r,0,m()));
2858  return const_iterator(this, r);
2859 }
2860 
2861 
2862 
2863 template <class MatrixType>
2864 inline
2867 {
2868  Assert (r<m(), ExcIndexRange(r,0,m()));
2869  return const_iterator(this, r+1);
2870 }
2871 
2872 
2873 
2874 template <class MatrixType>
2875 inline
2878 {
2879  return iterator(this, 0);
2880 }
2881 
2882 
2883 
2884 template <class MatrixType>
2885 inline
2888 {
2889  return iterator(this, m());
2890 }
2891 
2892 
2893 
2894 template <class MatrixType>
2895 inline
2898 {
2899  Assert (r<m(), ExcIndexRange(r,0,m()));
2900  return iterator(this, r);
2901 }
2902 
2903 
2904 
2905 template <class MatrixType>
2906 inline
2909 {
2910  Assert (r<m(), ExcIndexRange(r,0,m()));
2911  return iterator(this, r+1);
2912 }
2913 
2914 
2915 
2916 template <class MatrixType>
2917 void
2919 {
2920  std::vector<size_type> row_sizes (this->n_block_rows());
2921  std::vector<size_type> col_sizes (this->n_block_cols());
2922 
2923  // first find out the row sizes
2924  // from the first block column
2925  for (unsigned int r=0; r<this->n_block_rows(); ++r)
2926  row_sizes[r] = sub_objects[r][0]->m();
2927  // then check that the following
2928  // block columns have the same
2929  // sizes
2930  for (unsigned int c=1; c<this->n_block_cols(); ++c)
2931  for (unsigned int r=0; r<this->n_block_rows(); ++r)
2932  Assert (row_sizes[r] == sub_objects[r][c]->m(),
2933  ExcIncompatibleRowNumbers (r,0,r,c));
2934 
2935  // finally initialize the row
2936  // indices with this array
2937  this->row_block_indices.reinit (row_sizes);
2938 
2939 
2940  // then do the same with the columns
2941  for (unsigned int c=0; c<this->n_block_cols(); ++c)
2942  col_sizes[c] = sub_objects[0][c]->n();
2943  for (unsigned int r=1; r<this->n_block_rows(); ++r)
2944  for (unsigned int c=0; c<this->n_block_cols(); ++c)
2945  Assert (col_sizes[c] == sub_objects[r][c]->n(),
2946  ExcIncompatibleRowNumbers (0,c,r,c));
2947 
2948  // finally initialize the row
2949  // indices with this array
2950  this->column_block_indices.reinit (col_sizes);
2951 }
2952 
2953 
2954 
2955 template <class MatrixType>
2956 void
2958 {
2959  for (unsigned int row=0; row<n_block_rows(); ++row)
2960  for (unsigned int col=0; col<n_block_cols(); ++col)
2961  block(row, col).prepare_add();
2962 }
2963 
2964 
2965 
2966 template <class MatrixType>
2967 void
2969 {
2970  for (unsigned int row=0; row<n_block_rows(); ++row)
2971  for (unsigned int col=0; col<n_block_cols(); ++col)
2972  block(row, col).prepare_set();
2973 }
2974 
2975 #endif // DOXYGEN
2976 
2977 
2978 DEAL_II_NAMESPACE_CLOSE
2979 
2980 #endif // __deal2__block_matrix_base_h
Iterator lower_bound(Iterator first, Iterator last, const T &val)
Definition: utilities.h:731
const types::global_dof_index invalid_size_type
Definition: types.h:211
static const unsigned int invalid_unsigned_int
Definition: types.h:191
void set(const size_type i, const size_type j, const value_type value)
value_type matrix_scalar_product(const BlockVectorType &u, const BlockVectorType &v) const
void vmult_block_nonblock(BlockVectorType &dst, const VectorType &src) const
std::vector< size_type > counter_within_block
number value_type
BlockMatrix::BlockType::const_iterator base_iterator
BlockIndices row_block_indices
std::vector< std::vector< size_type > > column_indices
size_type m() const
Subscriptor & operator=(const Subscriptor &)
::ExceptionBase & ExcMessage(std::string arg1)
Auxiliary class aiding in the handling of block structures like in BlockVector or FESystem...
Definition: block_indices.h:55
const BlockIndices & get_column_indices() const
unsigned int block_row() const
Table< 2, SmartPointer< BlockType, BlockMatrixBase< MatrixType > > > sub_objects
STL namespace.
void vmult_nonblock_nonblock(VectorType &dst, const VectorType &src) const
void vmult_block_block(BlockVectorType &dst, const BlockVectorType &src) const
std::size_t memory_consumption() const
size_type n() const
value_type operator()(const size_type i, const size_type j) const
void prepare_add_operation()
bool is_finite(const double x)
value_type diag_element(const size_type i) const
void Tvmult_block_block(BlockVectorType &dst, const BlockVectorType &src) const
void prepare_set_operation()
void print(std::ostream &out, const bool alternative_output=false) const
void add(const size_type i, const size_type j, const value_type value)
unsigned int n_block_rows() const
void vmult_add(BlockVectorType &dst, const BlockVectorType &src) const
BlockMatrixBase & operator/=(const value_type factor)
unsigned int global_dof_index
Definition: types.h:100
#define Assert(cond, exc)
Definition: exceptions.h:299
unsigned int n_block_cols() const
std::size_t memory_consumption(const T &t)
BlockMatrix::value_type value_type
::ExceptionBase & ExcIndexRange(int arg1, int arg2, int arg3)
BlockType & block(const unsigned int row, const unsigned int column)
value_type residual(BlockVectorType &dst, const BlockVectorType &x, const BlockVectorType &b) const
BlockCompressedSparsityPattern CompressedBlockSparsityPattern DEAL_II_DEPRECATED
BlockType::value_type value_type
value_type el(const size_type i, const size_type j) const
void Tvmult_add(BlockVectorType &dst, const BlockVectorType &src) const
unsigned int block_column() const
DeclException4(ExcIncompatibleRowNumbers, int, int, int, int,<< "The blocks ["<< arg1<< ','<< arg2<< "] and ["<< arg3<< ','<< arg4<< "] have differing row numbers.")
void compress() DEAL_II_DEPRECATED
::ExceptionBase & ExcIteratorPastEnd()
void vmult_nonblock_block(VectorType &dst, const BlockVectorType &src) const
void Tvmult_nonblock_block(VectorType &dst, const BlockVectorType &src) const
BlockMatrixBase & operator*=(const value_type factor)
::ExceptionBase & ExcNumberNotFinite()
size_type m() const
iterator end()
TemporaryData temporary_data
void Tvmult_block_nonblock(BlockVectorType &dst, const VectorType &src) const
std::vector< std::vector< double > > column_values
void collect_sizes()
size_type n() const
::ExceptionBase & ExcNotImplemented()
::ExceptionBase & ExcNotInitialized()
::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
Definition: table.h:33
types::global_dof_index size_type
::ExceptionBase & ExcInternalError()
::ExceptionBase & ExcDivideByZero()
void Tvmult_nonblock_nonblock(VectorType &dst, const VectorType &src) const
const BlockIndices & get_row_indices() const
value_type matrix_norm_square(const BlockVectorType &v) const
iterator begin()
BlockMatrixBase & copy_from(const BlockMatrixType &source)