Reference documentation for deal.II version 8.4.1
block_matrix_base.h
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2004 - 2016 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 at
12 // the top level of the deal.II distribution.
13 //
14 // ---------------------------------------------------------------------
15 
16 #ifndef dealii__block_matrix_base_h
17 #define dealii__block_matrix_base_h
18 
19 
20 #include <deal.II/base/config.h>
21 #include <deal.II/base/table.h>
22 #include <deal.II/base/thread_management.h>
23 #include <deal.II/base/utilities.h>
24 #include <deal.II/base/smartpointer.h>
25 #include <deal.II/base/memory_consumption.h>
26 #include <deal.II/lac/block_indices.h>
27 #include <deal.II/lac/exceptions.h>
28 #include <deal.II/lac/full_matrix.h>
29 #include <deal.II/lac/matrix_iterator.h>
30 #include <deal.II/lac/vector.h>
31 
32 #include <cmath>
33 
34 DEAL_II_NAMESPACE_OPEN
35 
36 
37 template <typename> class MatrixIterator;
38 
39 
40 
51 {
56  template <class BlockMatrixType>
58  {
59  public:
64 
68  typedef typename BlockMatrixType::value_type value_type;
69 
73  AccessorBase ();
74 
78  unsigned int block_row() const;
79 
83  unsigned int block_column() const;
84 
85  protected:
89  unsigned int row_block;
90 
94  unsigned int col_block;
95 
99  template <typename>
100  friend class MatrixIterator;
101  };
102 
103 
104 
108  template <class BlockMatrixType, bool Constness>
109  class Accessor;
110 
111 
115  template <class BlockMatrixType>
116  class Accessor<BlockMatrixType, false>
117  :
118  public AccessorBase<BlockMatrixType>
119  {
120  public:
125 
129  typedef BlockMatrixType MatrixType;
130 
134  typedef typename BlockMatrixType::value_type value_type;
135 
144  Accessor (BlockMatrixType *m,
145  const size_type row,
146  const size_type col);
147 
151  size_type row() const;
152 
156  size_type column() const;
157 
161  value_type value() const;
162 
166  void set_value(value_type newval) const;
167 
168  protected:
172  BlockMatrixType *matrix;
173 
177  typename BlockMatrixType::BlockType::iterator base_iterator;
178 
182  void advance ();
183 
187  bool operator == (const Accessor &a) const;
188 
189  template <typename> friend class MatrixIterator;
190  friend class Accessor<BlockMatrixType, true>;
191  };
192 
197  template <class BlockMatrixType>
198  class Accessor<BlockMatrixType, true>
199  :
200  public AccessorBase<BlockMatrixType>
201  {
202  public:
207 
211  typedef const BlockMatrixType MatrixType;
212 
216  typedef typename BlockMatrixType::value_type value_type;
217 
226  Accessor (const BlockMatrixType *m,
227  const size_type row,
228  const size_type col);
229 
234 
238  size_type row() const;
239 
243  size_type column() const;
244 
248  value_type value() const;
249  protected:
253  const BlockMatrixType *matrix;
254 
258  typename BlockMatrixType::BlockType::const_iterator base_iterator;
259 
263  void advance ();
264 
268  bool operator == (const Accessor &a) const;
269 
273  template <typename>
274  friend class ::MatrixIterator;
275  };
276 }
277 
278 
279 
340 template <typename MatrixType>
342 {
343 public:
347  typedef MatrixType BlockType;
348 
354  typedef value_type *pointer;
355  typedef const value_type *const_pointer;
356  typedef value_type &reference;
357  typedef const value_type &const_reference;
358  typedef types::global_dof_index size_type;
359 
360  typedef
362  iterator;
363 
364  typedef
367 
368 
372  BlockMatrixBase ();
373 
377  ~BlockMatrixBase ();
378 
395  template <class BlockMatrixType>
397  copy_from (const BlockMatrixType &source);
398 
402  BlockType &
403  block (const unsigned int row,
404  const unsigned int column);
405 
406 
411  const BlockType &
412  block (const unsigned int row,
413  const unsigned int column) const;
414 
419  size_type m () const;
420 
425  size_type n () const;
426 
427 
432  unsigned int n_block_rows () const;
433 
438  unsigned int n_block_cols () const;
439 
445  void set (const size_type i,
446  const size_type j,
447  const value_type value);
448 
464  template <typename number>
465  void set (const std::vector<size_type> &indices,
466  const FullMatrix<number> &full_matrix,
467  const bool elide_zero_values = false);
468 
474  template <typename number>
475  void set (const std::vector<size_type> &row_indices,
476  const std::vector<size_type> &col_indices,
477  const FullMatrix<number> &full_matrix,
478  const bool elide_zero_values = false);
479 
490  template <typename number>
491  void set (const size_type row,
492  const std::vector<size_type> &col_indices,
493  const std::vector<number> &values,
494  const bool elide_zero_values = false);
495 
505  template <typename number>
506  void set (const size_type row,
507  const size_type n_cols,
508  const size_type *col_indices,
509  const number *values,
510  const bool elide_zero_values = false);
511 
517  void add (const size_type i,
518  const size_type j,
519  const value_type value);
520 
535  template <typename number>
536  void add (const std::vector<size_type> &indices,
537  const FullMatrix<number> &full_matrix,
538  const bool elide_zero_values = true);
539 
545  template <typename number>
546  void add (const std::vector<size_type> &row_indices,
547  const std::vector<size_type> &col_indices,
548  const FullMatrix<number> &full_matrix,
549  const bool elide_zero_values = true);
550 
560  template <typename number>
561  void add (const size_type row,
562  const std::vector<size_type> &col_indices,
563  const std::vector<number> &values,
564  const bool elide_zero_values = true);
565 
575  template <typename number>
576  void add (const size_type row,
577  const size_type n_cols,
578  const size_type *col_indices,
579  const number *values,
580  const bool elide_zero_values = true,
581  const bool col_indices_are_sorted = false);
582 
594  void add (const value_type factor,
595  const BlockMatrixBase<MatrixType> &matrix);
596 
603  value_type operator () (const size_type i,
604  const size_type j) const;
605 
614  value_type el (const size_type i,
615  const size_type j) const;
616 
627  value_type diag_element (const size_type i) const;
628 
637  void compress (::VectorOperation::values operation);
638 
642  BlockMatrixBase &operator *= (const value_type factor);
643 
647  BlockMatrixBase &operator /= (const value_type factor);
648 
653  template <class BlockVectorType>
654  void vmult_add (BlockVectorType &dst,
655  const BlockVectorType &src) const;
656 
662  template <class BlockVectorType>
663  void Tvmult_add (BlockVectorType &dst,
664  const BlockVectorType &src) const;
665 
678  template <class BlockVectorType>
679  value_type
680  matrix_norm_square (const BlockVectorType &v) const;
681 
685  template <class BlockVectorType>
686  value_type
687  matrix_scalar_product (const BlockVectorType &u,
688  const BlockVectorType &v) const;
689 
693  template <class BlockVectorType>
694  value_type residual (BlockVectorType &dst,
695  const BlockVectorType &x,
696  const BlockVectorType &b) const;
697 
704  void print (std::ostream &out,
705  const bool alternative_output = false) const;
706 
710  iterator begin ();
711 
715  iterator end ();
716 
720  iterator begin (const size_type r);
721 
725  iterator end (const size_type r);
729  const_iterator begin () const;
730 
734  const_iterator end () const;
735 
739  const_iterator begin (const size_type r) const;
740 
744  const_iterator end (const size_type r) const;
745 
749  const BlockIndices &get_row_indices () const;
750 
754  const BlockIndices &get_column_indices () const;
755 
761  std::size_t memory_consumption () const;
762 
771  DeclException4 (ExcIncompatibleRowNumbers,
772  int, int, int, int,
773  << "The blocks [" << arg1 << ',' << arg2 << "] and ["
774  << arg3 << ',' << arg4 << "] have differing row numbers.");
778  DeclException4 (ExcIncompatibleColNumbers,
779  int, int, int, int,
780  << "The blocks [" << arg1 << ',' << arg2 << "] and ["
781  << arg3 << ',' << arg4 << "] have differing column numbers.");
783 protected:
796  void clear ();
797 
802  BlockIndices column_block_indices;
803 
808 
827  void collect_sizes ();
828 
839  template <class BlockVectorType>
840  void vmult_block_block (BlockVectorType &dst,
841  const BlockVectorType &src) const;
842 
853  template <class BlockVectorType,
854  class VectorType>
855  void vmult_block_nonblock (BlockVectorType &dst,
856  const VectorType &src) const;
857 
868  template <class BlockVectorType,
869  class VectorType>
870  void vmult_nonblock_block (VectorType &dst,
871  const BlockVectorType &src) const;
872 
883  template <class VectorType>
884  void vmult_nonblock_nonblock (VectorType &dst,
885  const VectorType &src) const;
886 
898  template <class BlockVectorType>
899  void Tvmult_block_block (BlockVectorType &dst,
900  const BlockVectorType &src) const;
901 
912  template <class BlockVectorType,
913  class VectorType>
914  void Tvmult_block_nonblock (BlockVectorType &dst,
915  const VectorType &src) const;
916 
927  template <class BlockVectorType,
928  class VectorType>
929  void Tvmult_nonblock_block (VectorType &dst,
930  const BlockVectorType &src) const;
931 
942  template <class VectorType>
943  void Tvmult_nonblock_nonblock (VectorType &dst,
944  const VectorType &src) const;
945 
946 
947 protected:
948 
955  void prepare_add_operation();
956 
961  void prepare_set_operation();
962 
963 
964 private:
965 
975  {
980  std::vector<size_type> counter_within_block;
981 
986  std::vector<std::vector<size_type> > column_indices;
987 
992  std::vector<std::vector<value_type> > column_values;
993 
999 
1010  {
1011  return *this;
1012  }
1013  };
1014 
1021  TemporaryData temporary_data;
1022 
1027  template <typename, bool>
1029 
1030  template <typename>
1031  friend class MatrixIterator;
1032 };
1033 
1034 
1037 #ifndef DOXYGEN
1038 /* ------------------------- Template functions ---------------------- */
1039 
1040 
1041 namespace BlockMatrixIterators
1042 {
1043  template <class BlockMatrixType>
1044  inline
1046  :
1047  row_block(0),
1048  col_block(0)
1049  {}
1050 
1051 
1052  template <class BlockMatrixType>
1053  inline
1054  unsigned int
1055  AccessorBase<BlockMatrixType>::block_row() const
1056  {
1057  Assert (row_block != numbers::invalid_unsigned_int,
1058  ExcIteratorPastEnd());
1059 
1060  return row_block;
1061  }
1062 
1063 
1064  template <class BlockMatrixType>
1065  inline
1066  unsigned int
1067  AccessorBase<BlockMatrixType>::block_column() const
1068  {
1069  Assert (col_block != numbers::invalid_unsigned_int,
1070  ExcIteratorPastEnd());
1071 
1072  return col_block;
1073  }
1074 
1075 
1076  template <class BlockMatrixType>
1077  inline
1078  Accessor<BlockMatrixType, true>::Accessor (
1079  const BlockMatrixType *matrix,
1080  const size_type row,
1081  const size_type col)
1082  :
1083  matrix(matrix),
1084  base_iterator(matrix->block(0,0).begin())
1085  {
1086  (void)col;
1087  Assert(col==0, ExcNotImplemented());
1088 
1089  // check if this is a regular row or
1090  // the end of the matrix
1091  if (row < matrix->m())
1092  {
1093  const std::pair<unsigned int,size_type> indices
1094  = matrix->row_block_indices.global_to_local(row);
1095 
1096  // find the first block that does
1097  // have an entry in this row
1098  for (unsigned int bc=0; bc<matrix->n_block_cols(); ++bc)
1099  {
1100  base_iterator
1101  = matrix->block(indices.first, bc).begin(indices.second);
1102  if (base_iterator !=
1103  matrix->block(indices.first, bc).end(indices.second))
1104  {
1105  this->row_block = indices.first;
1106  this->col_block = bc;
1107  return;
1108  }
1109  }
1110 
1111  // hm, there is no block that has
1112  // an entry in this column. we need
1113  // to take the next entry then,
1114  // which may be the first entry of
1115  // the next row, or recursively the
1116  // next row, or so on
1117  *this = Accessor (matrix, row+1, 0);
1118  }
1119  else
1120  {
1121  // we were asked to create the end
1122  // iterator for this matrix
1123  this->row_block = numbers::invalid_unsigned_int;
1124  this->col_block = numbers::invalid_unsigned_int;
1125  }
1126  }
1127 
1128 
1129 // template <class BlockMatrixType>
1130 // inline
1131 // Accessor<BlockMatrixType, true>::Accessor (const Accessor<BlockMatrixType, true>& other)
1132 // :
1133 // matrix(other.matrix),
1134 // base_iterator(other.base_iterator)
1135 // {
1136 // this->row_block = other.row_block;
1137 // this->col_block = other.col_block;
1138 // }
1139 
1140 
1141  template <class BlockMatrixType>
1142  inline
1143  Accessor<BlockMatrixType, true>::Accessor (const Accessor<BlockMatrixType, false> &other)
1144  :
1145  matrix(other.matrix),
1146  base_iterator(other.base_iterator)
1147  {
1148  this->row_block = other.row_block;
1149  this->col_block = other.col_block;
1150  }
1151 
1152 
1153  template <class BlockMatrixType>
1154  inline
1155  typename Accessor<BlockMatrixType, true>::size_type
1156  Accessor<BlockMatrixType, true>::row() const
1157  {
1158  Assert (this->row_block != numbers::invalid_unsigned_int,
1159  ExcIteratorPastEnd());
1160 
1161  return (matrix->row_block_indices.local_to_global(this->row_block, 0) +
1162  base_iterator->row());
1163  }
1164 
1165 
1166  template <class BlockMatrixType>
1167  inline
1168  typename Accessor<BlockMatrixType, true>::size_type
1169  Accessor<BlockMatrixType, true>::column() const
1170  {
1171  Assert (this->col_block != numbers::invalid_unsigned_int,
1172  ExcIteratorPastEnd());
1173 
1174  return (matrix->column_block_indices.local_to_global(this->col_block,0) +
1175  base_iterator->column());
1176  }
1177 
1178 
1179  template <class BlockMatrixType>
1180  inline
1181  typename Accessor<BlockMatrixType, true>::value_type
1182  Accessor<BlockMatrixType, true>::value () const
1183  {
1184  Assert (this->row_block != numbers::invalid_unsigned_int,
1185  ExcIteratorPastEnd());
1186  Assert (this->col_block != numbers::invalid_unsigned_int,
1187  ExcIteratorPastEnd());
1188 
1189  return base_iterator->value();
1190  }
1191 
1192 
1193 
1194  template <class BlockMatrixType>
1195  inline
1196  void
1197  Accessor<BlockMatrixType, true>::advance ()
1198  {
1199  Assert (this->row_block != numbers::invalid_unsigned_int,
1200  ExcIteratorPastEnd());
1201  Assert (this->col_block != numbers::invalid_unsigned_int,
1202  ExcIteratorPastEnd());
1203 
1204  // Remember current row inside block
1205  size_type local_row = base_iterator->row();
1206 
1207  // Advance one element inside the
1208  // current block
1209  ++base_iterator;
1210 
1211  // while we hit the end of the row of a
1212  // block (which may happen multiple
1213  // times if rows inside a block are
1214  // empty), we have to jump to the next
1215  // block and take the
1216  while (base_iterator ==
1217  matrix->block(this->row_block, this->col_block).end(local_row))
1218  {
1219  // jump to next block in this block
1220  // row, if possible, otherwise go
1221  // to next row
1222  if (this->col_block < matrix->n_block_cols()-1)
1223  {
1224  ++this->col_block;
1225  base_iterator
1226  = matrix->block(this->row_block, this->col_block).begin(local_row);
1227  }
1228  else
1229  {
1230  // jump back to next row in
1231  // first block column
1232  this->col_block = 0;
1233  ++local_row;
1234 
1235  // see if this has brought us
1236  // past the number of rows in
1237  // this block. if so see
1238  // whether we've just fallen
1239  // off the end of the whole
1240  // matrix
1241  if (local_row == matrix->block(this->row_block, this->col_block).m())
1242  {
1243  local_row = 0;
1244  ++this->row_block;
1245  if (this->row_block == matrix->n_block_rows())
1246  {
1247  this->row_block = numbers::invalid_unsigned_int;
1248  this->col_block = numbers::invalid_unsigned_int;
1249  return;
1250  }
1251  }
1252 
1253  base_iterator
1254  = matrix->block(this->row_block, this->col_block).begin(local_row);
1255  }
1256  }
1257  }
1258 
1259 
1260  template <class BlockMatrixType>
1261  inline
1262  bool
1263  Accessor<BlockMatrixType, true>::operator == (const Accessor &a) const
1264  {
1265  if (matrix != a.matrix)
1266  return false;
1267 
1268  if (this->row_block == a.row_block
1269  && this->col_block == a.col_block)
1270  // end iterators do not necessarily
1271  // have to have the same
1272  // base_iterator representation, but
1273  // valid iterators have to
1274  return (((this->row_block == numbers::invalid_unsigned_int)
1275  &&
1276  (this->col_block == numbers::invalid_unsigned_int))
1277  ||
1278  (base_iterator == a.base_iterator));
1279 
1280  return false;
1281  }
1282 
1283 //----------------------------------------------------------------------//
1284 
1285 
1286  template <class BlockMatrixType>
1287  inline
1288  Accessor<BlockMatrixType, false>::Accessor (
1289  BlockMatrixType *matrix,
1290  const size_type row,
1291  const size_type col)
1292  :
1293  matrix(matrix),
1294  base_iterator(matrix->block(0,0).begin())
1295  {
1296  (void)col;
1297  Assert(col==0, ExcNotImplemented());
1298  // check if this is a regular row or
1299  // the end of the matrix
1300  if (row < matrix->m())
1301  {
1302  const std::pair<unsigned int,size_type> indices
1303  = matrix->row_block_indices.global_to_local(row);
1304 
1305  // find the first block that does
1306  // have an entry in this row
1307  for (size_type bc=0; bc<matrix->n_block_cols(); ++bc)
1308  {
1309  base_iterator
1310  = matrix->block(indices.first, bc).begin(indices.second);
1311  if (base_iterator !=
1312  matrix->block(indices.first, bc).end(indices.second))
1313  {
1314  this->row_block = indices.first;
1315  this->col_block = bc;
1316  return;
1317  }
1318  }
1319 
1320  // hm, there is no block that has
1321  // an entry in this column. we need
1322  // to take the next entry then,
1323  // which may be the first entry of
1324  // the next row, or recursively the
1325  // next row, or so on
1326  *this = Accessor (matrix, row+1, 0);
1327  }
1328  else
1329  {
1330  // we were asked to create the end
1331  // iterator for this matrix
1332  this->row_block = numbers::invalid_size_type;
1333  this->col_block = numbers::invalid_size_type;
1334  }
1335  }
1336 
1337 
1338  template <class BlockMatrixType>
1339  inline
1340  typename Accessor<BlockMatrixType, false>::size_type
1341  Accessor<BlockMatrixType, false>::row() const
1342  {
1343  Assert (this->row_block != numbers::invalid_size_type,
1344  ExcIteratorPastEnd());
1345 
1346  return (matrix->row_block_indices.local_to_global(this->row_block, 0) +
1347  base_iterator->row());
1348  }
1349 
1350 
1351  template <class BlockMatrixType>
1352  inline
1353  typename Accessor<BlockMatrixType, false>::size_type
1354  Accessor<BlockMatrixType, false>::column() const
1355  {
1356  Assert (this->col_block != numbers::invalid_size_type,
1357  ExcIteratorPastEnd());
1358 
1359  return (matrix->column_block_indices.local_to_global(this->col_block,0) +
1360  base_iterator->column());
1361  }
1362 
1363 
1364  template <class BlockMatrixType>
1365  inline
1366  typename Accessor<BlockMatrixType, false>::value_type
1367  Accessor<BlockMatrixType, false>::value () const
1368  {
1369  Assert (this->row_block != numbers::invalid_size_type,
1370  ExcIteratorPastEnd());
1371  Assert (this->col_block != numbers::invalid_size_type,
1372  ExcIteratorPastEnd());
1373 
1374  return base_iterator->value();
1375  }
1376 
1377 
1378 
1379  template <class BlockMatrixType>
1380  inline
1381  void
1382  Accessor<BlockMatrixType, false>::set_value (typename Accessor<BlockMatrixType, false>::value_type newval) const
1383  {
1384  Assert (this->row_block != numbers::invalid_size_type,
1385  ExcIteratorPastEnd());
1386  Assert (this->col_block != numbers::invalid_size_type,
1387  ExcIteratorPastEnd());
1388 
1389  base_iterator->value() = newval;
1390  }
1391 
1392 
1393 
1394  template <class BlockMatrixType>
1395  inline
1396  void
1397  Accessor<BlockMatrixType, false>::advance ()
1398  {
1399  Assert (this->row_block != numbers::invalid_size_type,
1400  ExcIteratorPastEnd());
1401  Assert (this->col_block != numbers::invalid_size_type,
1402  ExcIteratorPastEnd());
1403 
1404  // Remember current row inside block
1405  size_type local_row = base_iterator->row();
1406 
1407  // Advance one element inside the
1408  // current block
1409  ++base_iterator;
1410 
1411  // while we hit the end of the row of a
1412  // block (which may happen multiple
1413  // times if rows inside a block are
1414  // empty), we have to jump to the next
1415  // block and take the
1416  while (base_iterator ==
1417  matrix->block(this->row_block, this->col_block).end(local_row))
1418  {
1419  // jump to next block in this block
1420  // row, if possible, otherwise go
1421  // to next row
1422  if (this->col_block < matrix->n_block_cols()-1)
1423  {
1424  ++this->col_block;
1425  base_iterator
1426  = matrix->block(this->row_block, this->col_block).begin(local_row);
1427  }
1428  else
1429  {
1430  // jump back to next row in
1431  // first block column
1432  this->col_block = 0;
1433  ++local_row;
1434 
1435  // see if this has brought us
1436  // past the number of rows in
1437  // this block. if so see
1438  // whether we've just fallen
1439  // off the end of the whole
1440  // matrix
1441  if (local_row == matrix->block(this->row_block, this->col_block).m())
1442  {
1443  local_row = 0;
1444  ++this->row_block;
1445  if (this->row_block == matrix->n_block_rows())
1446  {
1447  this->row_block = numbers::invalid_size_type;
1448  this->col_block = numbers::invalid_size_type;
1449  return;
1450  }
1451  }
1452 
1453  base_iterator
1454  = matrix->block(this->row_block, this->col_block).begin(local_row);
1455  }
1456  }
1457  }
1458 
1459 
1460 
1461  template <class BlockMatrixType>
1462  inline
1463  bool
1464  Accessor<BlockMatrixType, false>::operator == (const Accessor &a) const
1465  {
1466  if (matrix != a.matrix)
1467  return false;
1468 
1469  if (this->row_block == a.row_block
1470  && this->col_block == a.col_block)
1471  // end iterators do not necessarily
1472  // have to have the same
1473  // base_iterator representation, but
1474  // valid iterators have to
1475  return (((this->row_block == numbers::invalid_size_type)
1476  &&
1477  (this->col_block == numbers::invalid_size_type))
1478  ||
1479  (base_iterator == a.base_iterator));
1480 
1481  return false;
1482  }
1483 }
1484 
1485 
1486 //---------------------------------------------------------------------------
1487 
1488 
1489 template <typename MatrixType>
1490 inline
1492 {}
1493 
1494 template <typename MatrixType>
1495 inline
1497 {
1498  clear ();
1499 }
1500 
1501 
1502 template <class MatrixType>
1503 template <class BlockMatrixType>
1504 inline
1507 copy_from (const BlockMatrixType &source)
1508 {
1509  for (unsigned int r=0; r<n_block_rows(); ++r)
1510  for (unsigned int c=0; c<n_block_cols(); ++c)
1511  block(r,c).copy_from (source.block(r,c));
1512 
1513  return *this;
1514 }
1515 
1516 
1517 template <class MatrixType>
1518 std::size_t
1520 {
1521  std::size_t mem =
1522  MemoryConsumption::memory_consumption(row_block_indices)+
1523  MemoryConsumption::memory_consumption(column_block_indices)+
1525  MemoryConsumption::memory_consumption(temporary_data.counter_within_block)+
1526  MemoryConsumption::memory_consumption(temporary_data.column_indices)+
1527  MemoryConsumption::memory_consumption(temporary_data.column_values)+
1528  sizeof(temporary_data.mutex);
1529 
1530  for (unsigned int r=0; r<n_block_rows(); ++r)
1531  for (unsigned int c=0; c<n_block_cols(); ++c)
1532  {
1533  MatrixType *p = this->sub_objects[r][c];
1535  }
1536 
1537  return mem;
1538 }
1539 
1540 
1541 
1542 template <class MatrixType>
1543 inline
1544 void
1546 {
1547  for (unsigned int r=0; r<n_block_rows(); ++r)
1548  for (unsigned int c=0; c<n_block_cols(); ++c)
1549  {
1550  MatrixType *p = this->sub_objects[r][c];
1551  this->sub_objects[r][c] = 0;
1552  delete p;
1553  }
1554  sub_objects.reinit (0,0);
1555 
1556  // reset block indices to empty
1557  row_block_indices = column_block_indices = BlockIndices ();
1558 }
1559 
1560 
1561 
1562 template <class MatrixType>
1563 inline
1565 BlockMatrixBase<MatrixType>::block (const unsigned int row,
1566  const unsigned int column)
1567 {
1568  Assert (row<n_block_rows(),
1569  ExcIndexRange (row, 0, n_block_rows()));
1570  Assert (column<n_block_cols(),
1571  ExcIndexRange (column, 0, n_block_cols()));
1572 
1573  return *sub_objects[row][column];
1574 }
1575 
1576 
1577 
1578 template <class MatrixType>
1579 inline
1581 BlockMatrixBase<MatrixType>::block (const unsigned int row,
1582  const unsigned int column) const
1583 {
1584  Assert (row<n_block_rows(),
1585  ExcIndexRange (row, 0, n_block_rows()));
1586  Assert (column<n_block_cols(),
1587  ExcIndexRange (column, 0, n_block_cols()));
1588 
1589  return *sub_objects[row][column];
1590 }
1591 
1592 
1593 template <class MatrixType>
1594 inline
1595 typename BlockMatrixBase<MatrixType>::size_type
1597 {
1598  return row_block_indices.total_size();
1599 }
1600 
1601 
1602 
1603 template <class MatrixType>
1604 inline
1605 typename BlockMatrixBase<MatrixType>::size_type
1607 {
1608  return column_block_indices.total_size();
1609 }
1610 
1611 
1612 
1613 template <class MatrixType>
1614 inline
1615 unsigned int
1617 {
1618  return column_block_indices.size();
1619 }
1620 
1621 
1622 
1623 template <class MatrixType>
1624 inline
1625 unsigned int
1627 {
1628  return row_block_indices.size();
1629 }
1630 
1631 
1632 
1633 // Write the single set manually,
1634 // since the other function has a lot
1635 // of overhead in that case.
1636 template <class MatrixType>
1637 inline
1638 void
1639 BlockMatrixBase<MatrixType>::set (const size_type i,
1640  const size_type j,
1641  const value_type value)
1642 {
1643  prepare_set_operation();
1644 
1645  AssertIsFinite(value);
1646 
1647  const std::pair<unsigned int,size_type>
1648  row_index = row_block_indices.global_to_local (i),
1649  col_index = column_block_indices.global_to_local (j);
1650  block(row_index.first,col_index.first).set (row_index.second,
1651  col_index.second,
1652  value);
1653 }
1654 
1655 
1656 
1657 template <class MatrixType>
1658 template <typename number>
1659 inline
1660 void
1661 BlockMatrixBase<MatrixType>::set (const std::vector<size_type> &row_indices,
1662  const std::vector<size_type> &col_indices,
1663  const FullMatrix<number> &values,
1664  const bool elide_zero_values)
1665 {
1666  Assert (row_indices.size() == values.m(),
1667  ExcDimensionMismatch(row_indices.size(), values.m()));
1668  Assert (col_indices.size() == values.n(),
1669  ExcDimensionMismatch(col_indices.size(), values.n()));
1670 
1671  for (size_type i=0; i<row_indices.size(); ++i)
1672  set (row_indices[i], col_indices.size(), &col_indices[0], &values(i,0),
1673  elide_zero_values);
1674 }
1675 
1676 
1677 
1678 template <class MatrixType>
1679 template <typename number>
1680 inline
1681 void
1682 BlockMatrixBase<MatrixType>::set (const std::vector<size_type> &indices,
1683  const FullMatrix<number> &values,
1684  const bool elide_zero_values)
1685 {
1686  Assert (indices.size() == values.m(),
1687  ExcDimensionMismatch(indices.size(), values.m()));
1688  Assert (values.n() == values.m(), ExcNotQuadratic());
1689 
1690  for (size_type i=0; i<indices.size(); ++i)
1691  set (indices[i], indices.size(), &indices[0], &values(i,0),
1692  elide_zero_values);
1693 }
1694 
1695 
1696 
1697 template <class MatrixType>
1698 template <typename number>
1699 inline
1700 void
1701 BlockMatrixBase<MatrixType>::set (const size_type row,
1702  const std::vector<size_type> &col_indices,
1703  const std::vector<number> &values,
1704  const bool elide_zero_values)
1705 {
1706  Assert (col_indices.size() == values.size(),
1707  ExcDimensionMismatch(col_indices.size(), values.size()));
1708 
1709  set (row, col_indices.size(), &col_indices[0], &values[0],
1710  elide_zero_values);
1711 }
1712 
1713 
1714 
1715 // This is a very messy function, since
1716 // we need to calculate to each position
1717 // the location in the global array.
1718 template <class MatrixType>
1719 template <typename number>
1720 inline
1721 void
1722 BlockMatrixBase<MatrixType>::set (const size_type row,
1723  const size_type n_cols,
1724  const size_type *col_indices,
1725  const number *values,
1726  const bool elide_zero_values)
1727 {
1728  prepare_set_operation();
1729 
1730  // lock access to the temporary data structure to
1731  // allow multiple threads to call this function concurrently
1732  Threads::Mutex::ScopedLock lock (temporary_data.mutex);
1733 
1734  // Resize scratch arrays
1735  if (temporary_data.column_indices.size() < this->n_block_cols())
1736  {
1737  temporary_data.column_indices.resize (this->n_block_cols());
1738  temporary_data.column_values.resize (this->n_block_cols());
1739  temporary_data.counter_within_block.resize (this->n_block_cols());
1740  }
1741 
1742  // Resize sub-arrays to n_cols. This
1743  // is a bit wasteful, but we resize
1744  // only a few times (then the maximum
1745  // row length won't increase that
1746  // much any more). At least we know
1747  // that all arrays are going to be of
1748  // the same size, so we can check
1749  // whether the size of one is large
1750  // enough before actually going
1751  // through all of them.
1752  if (temporary_data.column_indices[0].size() < n_cols)
1753  {
1754  for (unsigned int i=0; i<this->n_block_cols(); ++i)
1755  {
1756  temporary_data.column_indices[i].resize(n_cols);
1757  temporary_data.column_values[i].resize(n_cols);
1758  }
1759  }
1760 
1761  // Reset the number of added elements
1762  // in each block to zero.
1763  for (unsigned int i=0; i<this->n_block_cols(); ++i)
1764  temporary_data.counter_within_block[i] = 0;
1765 
1766  // Go through the column indices to
1767  // find out which portions of the
1768  // values should be set in which
1769  // block of the matrix. We need to
1770  // touch all the data, since we can't
1771  // be sure that the data of one block
1772  // is stored contiguously (in fact,
1773  // indices will be intermixed when it
1774  // comes from an element matrix).
1775  for (size_type j=0; j<n_cols; ++j)
1776  {
1777  number value = values[j];
1778 
1779  if (value == number() && elide_zero_values == true)
1780  continue;
1781 
1782  const std::pair<unsigned int, size_type>
1783  col_index = this->column_block_indices.global_to_local(col_indices[j]);
1784 
1785  const size_type local_index = temporary_data.counter_within_block[col_index.first]++;
1786 
1787  temporary_data.column_indices[col_index.first][local_index] = col_index.second;
1788  temporary_data.column_values[col_index.first][local_index] = value;
1789  }
1790 
1791 #ifdef DEBUG
1792  // If in debug mode, do a check whether
1793  // the right length has been obtained.
1794  size_type length = 0;
1795  for (unsigned int i=0; i<this->n_block_cols(); ++i)
1796  length += temporary_data.counter_within_block[i];
1797  Assert (length <= n_cols, ExcInternalError());
1798 #endif
1799 
1800  // Now we found out about where the
1801  // individual columns should start and
1802  // where we should start reading out
1803  // data. Now let's write the data into
1804  // the individual blocks!
1805  const std::pair<unsigned int,size_type>
1806  row_index = this->row_block_indices.global_to_local (row);
1807  for (unsigned int block_col=0; block_col<n_block_cols(); ++block_col)
1808  {
1809  if (temporary_data.counter_within_block[block_col] == 0)
1810  continue;
1811 
1812  block(row_index.first, block_col).set
1813  (row_index.second,
1814  temporary_data.counter_within_block[block_col],
1815  &temporary_data.column_indices[block_col][0],
1816  &temporary_data.column_values[block_col][0],
1817  false);
1818  }
1819 }
1820 
1821 
1822 
1823 template <class MatrixType>
1824 inline
1825 void
1826 BlockMatrixBase<MatrixType>::add (const size_type i,
1827  const size_type j,
1828  const value_type value)
1829 {
1830 
1831  AssertIsFinite(value);
1832 
1833  prepare_add_operation();
1834 
1835  // save some cycles for zero additions, but
1836  // only if it is safe for the matrix we are
1837  // working with
1838  typedef typename MatrixType::Traits MatrixTraits;
1839  if ((MatrixTraits::zero_addition_can_be_elided == true)
1840  &&
1841  (value == value_type()))
1842  return;
1843 
1844  const std::pair<unsigned int,size_type>
1845  row_index = row_block_indices.global_to_local (i),
1846  col_index = column_block_indices.global_to_local (j);
1847  block(row_index.first,col_index.first).add (row_index.second,
1848  col_index.second,
1849  value);
1850 }
1851 
1852 
1853 
1854 template <class MatrixType>
1855 template <typename number>
1856 inline
1857 void
1858 BlockMatrixBase<MatrixType>::add (const std::vector<size_type> &row_indices,
1859  const std::vector<size_type> &col_indices,
1860  const FullMatrix<number> &values,
1861  const bool elide_zero_values)
1862 {
1863  Assert (row_indices.size() == values.m(),
1864  ExcDimensionMismatch(row_indices.size(), values.m()));
1865  Assert (col_indices.size() == values.n(),
1866  ExcDimensionMismatch(col_indices.size(), values.n()));
1867 
1868  for (size_type i=0; i<row_indices.size(); ++i)
1869  add (row_indices[i], col_indices.size(), &col_indices[0], &values(i,0),
1870  elide_zero_values);
1871 }
1872 
1873 
1874 
1875 template <class MatrixType>
1876 template <typename number>
1877 inline
1878 void
1879 BlockMatrixBase<MatrixType>::add (const std::vector<size_type> &indices,
1880  const FullMatrix<number> &values,
1881  const bool elide_zero_values)
1882 {
1883  Assert (indices.size() == values.m(),
1884  ExcDimensionMismatch(indices.size(), values.m()));
1885  Assert (values.n() == values.m(), ExcNotQuadratic());
1886 
1887  for (size_type i=0; i<indices.size(); ++i)
1888  add (indices[i], indices.size(), &indices[0], &values(i,0),
1889  elide_zero_values);
1890 }
1891 
1892 
1893 
1894 template <class MatrixType>
1895 template <typename number>
1896 inline
1897 void
1898 BlockMatrixBase<MatrixType>::add (const size_type row,
1899  const std::vector<size_type> &col_indices,
1900  const std::vector<number> &values,
1901  const bool elide_zero_values)
1902 {
1903  Assert (col_indices.size() == values.size(),
1904  ExcDimensionMismatch(col_indices.size(), values.size()));
1905 
1906  add (row, col_indices.size(), &col_indices[0], &values[0],
1907  elide_zero_values);
1908 }
1909 
1910 
1911 
1912 // This is a very messy function, since
1913 // we need to calculate to each position
1914 // the location in the global array.
1915 template <class MatrixType>
1916 template <typename number>
1917 inline
1918 void
1919 BlockMatrixBase<MatrixType>::add (const size_type row,
1920  const size_type n_cols,
1921  const size_type *col_indices,
1922  const number *values,
1923  const bool elide_zero_values,
1924  const bool col_indices_are_sorted)
1925 {
1926  prepare_add_operation();
1927 
1928  // TODO: Look over this to find out
1929  // whether we can do that more
1930  // efficiently.
1931  if (col_indices_are_sorted == true)
1932  {
1933 #ifdef DEBUG
1934  // check whether indices really are
1935  // sorted.
1936  size_type before = col_indices[0];
1937  for (size_type i=1; i<n_cols; ++i)
1938  if (col_indices[i] <= before)
1939  Assert (false, ExcMessage ("Flag col_indices_are_sorted is set, but "
1940  "indices appear to not be sorted."))
1941  else
1942  before = col_indices[i];
1943 #endif
1944  const std::pair<unsigned int,size_type>
1945  row_index = this->row_block_indices.global_to_local (row);
1946 
1947  if (this->n_block_cols() > 1)
1948  {
1949  const size_type *first_block = Utilities::lower_bound (col_indices,
1950  col_indices+n_cols,
1951  this->column_block_indices.block_start(1));
1952 
1953  const size_type n_zero_block_indices = first_block - col_indices;
1954  block(row_index.first, 0).add (row_index.second,
1955  n_zero_block_indices,
1956  col_indices,
1957  values,
1958  elide_zero_values,
1959  col_indices_are_sorted);
1960 
1961  if (n_zero_block_indices < n_cols)
1962  this->add(row, n_cols - n_zero_block_indices, first_block,
1963  values + n_zero_block_indices, elide_zero_values,
1964  false);
1965  }
1966  else
1967  {
1968  block(row_index.first, 0). add (row_index.second,
1969  n_cols,
1970  col_indices,
1971  values,
1972  elide_zero_values,
1973  col_indices_are_sorted);
1974  }
1975 
1976  return;
1977  }
1978 
1979  // Lock scratch arrays, then resize them
1980  Threads::Mutex::ScopedLock lock (temporary_data.mutex);
1981 
1982  if (temporary_data.column_indices.size() < this->n_block_cols())
1983  {
1984  temporary_data.column_indices.resize (this->n_block_cols());
1985  temporary_data.column_values.resize (this->n_block_cols());
1986  temporary_data.counter_within_block.resize (this->n_block_cols());
1987  }
1988 
1989  // Resize sub-arrays to n_cols. This
1990  // is a bit wasteful, but we resize
1991  // only a few times (then the maximum
1992  // row length won't increase that
1993  // much any more). At least we know
1994  // that all arrays are going to be of
1995  // the same size, so we can check
1996  // whether the size of one is large
1997  // enough before actually going
1998  // through all of them.
1999  if (temporary_data.column_indices[0].size() < n_cols)
2000  {
2001  for (unsigned int i=0; i<this->n_block_cols(); ++i)
2002  {
2003  temporary_data.column_indices[i].resize(n_cols);
2004  temporary_data.column_values[i].resize(n_cols);
2005  }
2006  }
2007 
2008  // Reset the number of added elements
2009  // in each block to zero.
2010  for (unsigned int i=0; i<this->n_block_cols(); ++i)
2011  temporary_data.counter_within_block[i] = 0;
2012 
2013  // Go through the column indices to
2014  // find out which portions of the
2015  // values should be written into
2016  // which block of the matrix. We need
2017  // to touch all the data, since we
2018  // can't be sure that the data of one
2019  // block is stored contiguously (in
2020  // fact, data will be intermixed when
2021  // it comes from an element matrix).
2022  for (size_type j=0; j<n_cols; ++j)
2023  {
2024  number value = values[j];
2025 
2026  if (value == number() && elide_zero_values == true)
2027  continue;
2028 
2029  const std::pair<unsigned int, size_type>
2030  col_index = this->column_block_indices.global_to_local(col_indices[j]);
2031 
2032  const size_type local_index = temporary_data.counter_within_block[col_index.first]++;
2033 
2034  temporary_data.column_indices[col_index.first][local_index] = col_index.second;
2035  temporary_data.column_values[col_index.first][local_index] = value;
2036  }
2037 
2038 #ifdef DEBUG
2039  // If in debug mode, do a check whether
2040  // the right length has been obtained.
2041  size_type length = 0;
2042  for (unsigned int i=0; i<this->n_block_cols(); ++i)
2043  length += temporary_data.counter_within_block[i];
2044  Assert (length <= n_cols, ExcInternalError());
2045 #endif
2046 
2047  // Now we found out about where the
2048  // individual columns should start and
2049  // where we should start reading out
2050  // data. Now let's write the data into
2051  // the individual blocks!
2052  const std::pair<unsigned int,size_type>
2053  row_index = this->row_block_indices.global_to_local (row);
2054  for (unsigned int block_col=0; block_col<n_block_cols(); ++block_col)
2055  {
2056  if (temporary_data.counter_within_block[block_col] == 0)
2057  continue;
2058 
2059  block(row_index.first, block_col).add
2060  (row_index.second,
2061  temporary_data.counter_within_block[block_col],
2062  &temporary_data.column_indices[block_col][0],
2063  &temporary_data.column_values[block_col][0],
2064  false,
2065  col_indices_are_sorted);
2066  }
2067 }
2068 
2069 
2070 
2071 template <class MatrixType>
2072 inline
2073 void
2074 BlockMatrixBase<MatrixType>::add (const value_type factor,
2075  const BlockMatrixBase<MatrixType> &matrix)
2076 {
2077  AssertIsFinite(factor);
2078 
2079  prepare_add_operation();
2080 
2081  // save some cycles for zero additions, but
2082  // only if it is safe for the matrix we are
2083  // working with
2084  typedef typename MatrixType::Traits MatrixTraits;
2085  if ((MatrixTraits::zero_addition_can_be_elided == true)
2086  &&
2087  (factor == 0))
2088  return;
2089 
2090  for (unsigned int row=0; row<n_block_rows(); ++row)
2091  for (unsigned int col=0; col<n_block_cols(); ++col)
2092  // This function should throw if the sparsity
2093  // patterns of the two blocks differ
2094  block(row, col).add(factor, matrix.block(row,col));
2095 }
2096 
2097 
2098 
2099 template <class MatrixType>
2100 inline
2103  const size_type j) const
2104 {
2105  const std::pair<unsigned int,size_type>
2106  row_index = row_block_indices.global_to_local (i),
2107  col_index = column_block_indices.global_to_local (j);
2108  return block(row_index.first,col_index.first) (row_index.second,
2109  col_index.second);
2110 }
2111 
2112 
2113 
2114 template <class MatrixType>
2115 inline
2117 BlockMatrixBase<MatrixType>::el (const size_type i,
2118  const size_type j) const
2119 {
2120  const std::pair<unsigned int,size_type>
2121  row_index = row_block_indices.global_to_local (i),
2122  col_index = column_block_indices.global_to_local (j);
2123  return block(row_index.first,col_index.first).el (row_index.second,
2124  col_index.second);
2125 }
2126 
2127 
2128 
2129 template <class MatrixType>
2130 inline
2132 BlockMatrixBase<MatrixType>::diag_element (const size_type i) const
2133 {
2134  Assert (n_block_rows() == n_block_cols(),
2135  ExcNotQuadratic());
2136 
2137  const std::pair<unsigned int,size_type>
2138  index = row_block_indices.global_to_local (i);
2139  return block(index.first,index.first).diag_element(index.second);
2140 }
2141 
2142 
2143 
2144 template <class MatrixType>
2145 inline
2146 void
2147 BlockMatrixBase<MatrixType>::compress (::VectorOperation::values operation)
2148 {
2149  for (unsigned int r=0; r<n_block_rows(); ++r)
2150  for (unsigned int c=0; c<n_block_cols(); ++c)
2151  block(r,c).compress (operation);
2152 }
2153 
2154 
2155 
2156 template <class MatrixType>
2157 inline
2159 BlockMatrixBase<MatrixType>::operator *= (const value_type factor)
2160 {
2161  Assert (n_block_cols() != 0, ExcNotInitialized());
2162  Assert (n_block_rows() != 0, ExcNotInitialized());
2163 
2164  for (unsigned int r=0; r<n_block_rows(); ++r)
2165  for (unsigned int c=0; c<n_block_cols(); ++c)
2166  block(r,c) *= factor;
2167 
2168  return *this;
2169 }
2170 
2171 
2172 
2173 template <class MatrixType>
2174 inline
2176 BlockMatrixBase<MatrixType>::operator /= (const value_type factor)
2177 {
2178  Assert (n_block_cols() != 0, ExcNotInitialized());
2179  Assert (n_block_rows() != 0, ExcNotInitialized());
2180  Assert (factor !=0, ExcDivideByZero());
2181 
2182  const value_type factor_inv = 1. / factor;
2183 
2184  for (unsigned int r=0; r<n_block_rows(); ++r)
2185  for (unsigned int c=0; c<n_block_cols(); ++c)
2186  block(r,c) *= factor_inv;
2187 
2188  return *this;
2189 }
2190 
2191 
2192 
2193 template <class MatrixType>
2194 const BlockIndices &
2196 {
2197  return this->row_block_indices;
2198 }
2199 
2200 
2201 
2202 template <class MatrixType>
2203 const BlockIndices &
2205 {
2206  return this->column_block_indices;
2207 }
2208 
2209 
2210 
2211 template <class MatrixType>
2212 template <class BlockVectorType>
2213 void
2215 vmult_block_block (BlockVectorType &dst,
2216  const BlockVectorType &src) const
2217 {
2218  Assert (dst.n_blocks() == n_block_rows(),
2219  ExcDimensionMismatch(dst.n_blocks(), n_block_rows()));
2220  Assert (src.n_blocks() == n_block_cols(),
2221  ExcDimensionMismatch(src.n_blocks(), n_block_cols()));
2222 
2223  for (size_type row=0; row<n_block_rows(); ++row)
2224  {
2225  block(row,0).vmult (dst.block(row),
2226  src.block(0));
2227  for (size_type col=1; col<n_block_cols(); ++col)
2228  block(row,col).vmult_add (dst.block(row),
2229  src.block(col));
2230  };
2231 }
2232 
2233 
2234 
2235 template <class MatrixType>
2236 template <class BlockVectorType,
2237  class VectorType>
2238 void
2240 vmult_nonblock_block (VectorType &dst,
2241  const BlockVectorType &src) const
2242 {
2243  Assert (n_block_rows() == 1,
2244  ExcDimensionMismatch(1, n_block_rows()));
2245  Assert (src.n_blocks() == n_block_cols(),
2246  ExcDimensionMismatch(src.n_blocks(), n_block_cols()));
2247 
2248  block(0,0).vmult (dst, src.block(0));
2249  for (size_type col=1; col<n_block_cols(); ++col)
2250  block(0,col).vmult_add (dst, src.block(col));
2251 }
2252 
2253 
2254 
2255 template <class MatrixType>
2256 template <class BlockVectorType,
2257  class VectorType>
2258 void
2260 vmult_block_nonblock (BlockVectorType &dst,
2261  const VectorType &src) const
2262 {
2263  Assert (dst.n_blocks() == n_block_rows(),
2264  ExcDimensionMismatch(dst.n_blocks(), n_block_rows()));
2265  Assert (1 == n_block_cols(),
2266  ExcDimensionMismatch(1, n_block_cols()));
2267 
2268  for (size_type row=0; row<n_block_rows(); ++row)
2269  block(row,0).vmult (dst.block(row),
2270  src);
2271 }
2272 
2273 
2274 
2275 template <class MatrixType>
2276 template <class VectorType>
2277 void
2279 vmult_nonblock_nonblock (VectorType &dst,
2280  const VectorType &src) const
2281 {
2282  Assert (1 == n_block_rows(),
2283  ExcDimensionMismatch(1, n_block_rows()));
2284  Assert (1 == n_block_cols(),
2285  ExcDimensionMismatch(1, n_block_cols()));
2286 
2287  block(0,0).vmult (dst, src);
2288 }
2289 
2290 
2291 
2292 template <class MatrixType>
2293 template <class BlockVectorType>
2294 void
2295 BlockMatrixBase<MatrixType>::vmult_add (BlockVectorType &dst,
2296  const BlockVectorType &src) const
2297 {
2298  Assert (dst.n_blocks() == n_block_rows(),
2299  ExcDimensionMismatch(dst.n_blocks(), n_block_rows()));
2300  Assert (src.n_blocks() == n_block_cols(),
2301  ExcDimensionMismatch(src.n_blocks(), n_block_cols()));
2302 
2303  for (unsigned int row=0; row<n_block_rows(); ++row)
2304  for (unsigned int col=0; col<n_block_cols(); ++col)
2305  block(row,col).vmult_add (dst.block(row),
2306  src.block(col));
2307 }
2308 
2309 
2310 
2311 
2312 template <class MatrixType>
2313 template <class BlockVectorType>
2314 void
2316 Tvmult_block_block (BlockVectorType &dst,
2317  const BlockVectorType &src) const
2318 {
2319  Assert (dst.n_blocks() == n_block_cols(),
2320  ExcDimensionMismatch(dst.n_blocks(), n_block_cols()));
2321  Assert (src.n_blocks() == n_block_rows(),
2322  ExcDimensionMismatch(src.n_blocks(), n_block_rows()));
2323 
2324  dst = 0.;
2325 
2326  for (unsigned int row=0; row<n_block_rows(); ++row)
2327  {
2328  for (unsigned int col=0; col<n_block_cols(); ++col)
2329  block(row,col).Tvmult_add (dst.block(col),
2330  src.block(row));
2331  };
2332 }
2333 
2334 
2335 
2336 template <class MatrixType>
2337 template <class BlockVectorType,
2338  class VectorType>
2339 void
2341 Tvmult_block_nonblock (BlockVectorType &dst,
2342  const VectorType &src) const
2343 {
2344  Assert (dst.n_blocks() == n_block_cols(),
2345  ExcDimensionMismatch(dst.n_blocks(), n_block_cols()));
2346  Assert (1 == n_block_rows(),
2347  ExcDimensionMismatch(1, n_block_rows()));
2348 
2349  dst = 0.;
2350 
2351  for (unsigned int col=0; col<n_block_cols(); ++col)
2352  block(0,col).Tvmult_add (dst.block(col), src);
2353 }
2354 
2355 
2356 
2357 template <class MatrixType>
2358 template <class BlockVectorType,
2359  class VectorType>
2360 void
2362 Tvmult_nonblock_block (VectorType &dst,
2363  const BlockVectorType &src) const
2364 {
2365  Assert (1 == n_block_cols(),
2366  ExcDimensionMismatch(1, n_block_cols()));
2367  Assert (src.n_blocks() == n_block_rows(),
2368  ExcDimensionMismatch(src.n_blocks(), n_block_rows()));
2369 
2370  block(0,0).Tvmult (dst, src.block(0));
2371 
2372  for (size_type row=1; row<n_block_rows(); ++row)
2373  block(row,0).Tvmult_add (dst, src.block(row));
2374 }
2375 
2376 
2377 
2378 template <class MatrixType>
2379 template <class VectorType>
2380 void
2382 Tvmult_nonblock_nonblock (VectorType &dst,
2383  const VectorType &src) const
2384 {
2385  Assert (1 == n_block_cols(),
2386  ExcDimensionMismatch(1, n_block_cols()));
2387  Assert (1 == n_block_rows(),
2388  ExcDimensionMismatch(1, n_block_rows()));
2389 
2390  block(0,0).Tvmult (dst, src);
2391 }
2392 
2393 
2394 
2395 template <class MatrixType>
2396 template <class BlockVectorType>
2397 void
2398 BlockMatrixBase<MatrixType>::Tvmult_add (BlockVectorType &dst,
2399  const BlockVectorType &src) const
2400 {
2401  Assert (dst.n_blocks() == n_block_cols(),
2402  ExcDimensionMismatch(dst.n_blocks(), n_block_cols()));
2403  Assert (src.n_blocks() == n_block_rows(),
2404  ExcDimensionMismatch(src.n_blocks(), n_block_rows()));
2405 
2406  for (unsigned int row=0; row<n_block_rows(); ++row)
2407  for (unsigned int col=0; col<n_block_cols(); ++col)
2408  block(row,col).Tvmult_add (dst.block(col),
2409  src.block(row));
2410 }
2411 
2412 
2413 
2414 template <class MatrixType>
2415 template <class BlockVectorType>
2417 BlockMatrixBase<MatrixType>::matrix_norm_square (const BlockVectorType &v) const
2418 {
2419  Assert (n_block_rows() == n_block_cols(), ExcNotQuadratic());
2420  Assert (v.n_blocks() == n_block_rows(),
2421  ExcDimensionMismatch(v.n_blocks(), n_block_rows()));
2422 
2423  value_type norm_sqr = 0;
2424  for (unsigned int row=0; row<n_block_rows(); ++row)
2425  for (unsigned int col=0; col<n_block_cols(); ++col)
2426  if (row==col)
2427  norm_sqr += block(row,col).matrix_norm_square (v.block(row));
2428  else
2429  norm_sqr += block(row,col).matrix_scalar_product (v.block(row),
2430  v.block(col));
2431  return norm_sqr;
2432 }
2433 
2434 
2435 
2436 template <class MatrixType>
2437 template <class BlockVectorType>
2440 matrix_scalar_product (const BlockVectorType &u,
2441  const BlockVectorType &v) const
2442 {
2443  Assert (u.n_blocks() == n_block_rows(),
2444  ExcDimensionMismatch(u.n_blocks(), n_block_rows()));
2445  Assert (v.n_blocks() == n_block_cols(),
2446  ExcDimensionMismatch(v.n_blocks(), n_block_cols()));
2447 
2448  value_type result = 0;
2449  for (unsigned int row=0; row<n_block_rows(); ++row)
2450  for (unsigned int col=0; col<n_block_cols(); ++col)
2451  result += block(row,col).matrix_scalar_product (u.block(row),
2452  v.block(col));
2453  return result;
2454 }
2455 
2456 
2457 
2458 template <class MatrixType>
2459 template <class BlockVectorType>
2462 residual (BlockVectorType &dst,
2463  const BlockVectorType &x,
2464  const BlockVectorType &b) const
2465 {
2466  Assert (dst.n_blocks() == n_block_rows(),
2467  ExcDimensionMismatch(dst.n_blocks(), n_block_rows()));
2468  Assert (b.n_blocks() == n_block_rows(),
2469  ExcDimensionMismatch(b.n_blocks(), n_block_rows()));
2470  Assert (x.n_blocks() == n_block_cols(),
2471  ExcDimensionMismatch(x.n_blocks(), n_block_cols()));
2472  // in block notation, the residual is
2473  // r_i = b_i - \sum_j A_ij x_j.
2474  // this can be written as
2475  // r_i = b_i - A_i0 x_0 - \sum_{j>0} A_ij x_j.
2476  //
2477  // for the first two terms, we can
2478  // call the residual function of
2479  // A_i0. for the other terms, we
2480  // use vmult_add. however, we want
2481  // to subtract, so in order to
2482  // avoid a temporary vector, we
2483  // perform a sign change of the
2484  // first two term before, and after
2485  // adding up
2486  for (unsigned int row=0; row<n_block_rows(); ++row)
2487  {
2488  block(row,0).residual (dst.block(row),
2489  x.block(0),
2490  b.block(row));
2491 
2492  for (size_type i=0; i<dst.block(row).size(); ++i)
2493  dst.block(row)(i) = -dst.block(row)(i);
2494 
2495  for (unsigned int col=1; col<n_block_cols(); ++col)
2496  block(row,col).vmult_add (dst.block(row),
2497  x.block(col));
2498 
2499  for (size_type i=0; i<dst.block(row).size(); ++i)
2500  dst.block(row)(i) = -dst.block(row)(i);
2501  };
2502 
2503  value_type res = 0;
2504  for (size_type row=0; row<n_block_rows(); ++row)
2505  res += dst.block(row).norm_sqr ();
2506  return std::sqrt(res);
2507 }
2508 
2509 
2510 
2511 template <class MatrixType>
2512 inline
2513 void
2514 BlockMatrixBase<MatrixType>::print (std::ostream &out,
2515  const bool alternative_output) const
2516 {
2517  for (unsigned int row=0; row<n_block_rows(); ++row)
2518  for (unsigned int col=0; col<n_block_cols(); ++col)
2519  {
2520  if (!alternative_output)
2521  out << "Block (" << row << ", " << col << ")" << std::endl;
2522 
2523  block(row, col).print(out, alternative_output);
2524  }
2525 }
2526 
2527 
2528 
2529 template <class MatrixType>
2530 inline
2533 {
2534  return const_iterator(this, 0);
2535 }
2536 
2537 
2538 
2539 template <class MatrixType>
2540 inline
2543 {
2544  return const_iterator(this, m());
2545 }
2546 
2547 
2548 
2549 template <class MatrixType>
2550 inline
2552 BlockMatrixBase<MatrixType>::begin (const size_type r) const
2553 {
2554  Assert (r<m(), ExcIndexRange(r,0,m()));
2555  return const_iterator(this, r);
2556 }
2557 
2558 
2559 
2560 template <class MatrixType>
2561 inline
2563 BlockMatrixBase<MatrixType>::end (const size_type r) const
2564 {
2565  Assert (r<m(), ExcIndexRange(r,0,m()));
2566  return const_iterator(this, r+1);
2567 }
2568 
2569 
2570 
2571 template <class MatrixType>
2572 inline
2575 {
2576  return iterator(this, 0);
2577 }
2578 
2579 
2580 
2581 template <class MatrixType>
2582 inline
2585 {
2586  return iterator(this, m());
2587 }
2588 
2589 
2590 
2591 template <class MatrixType>
2592 inline
2594 BlockMatrixBase<MatrixType>::begin (const size_type r)
2595 {
2596  Assert (r<m(), ExcIndexRange(r,0,m()));
2597  return iterator(this, r);
2598 }
2599 
2600 
2601 
2602 template <class MatrixType>
2603 inline
2605 BlockMatrixBase<MatrixType>::end (const size_type r)
2606 {
2607  Assert (r<m(), ExcIndexRange(r,0,m()));
2608  return iterator(this, r+1);
2609 }
2610 
2611 
2612 
2613 template <class MatrixType>
2614 void
2616 {
2617  std::vector<size_type> row_sizes (this->n_block_rows());
2618  std::vector<size_type> col_sizes (this->n_block_cols());
2619 
2620  // first find out the row sizes
2621  // from the first block column
2622  for (unsigned int r=0; r<this->n_block_rows(); ++r)
2623  row_sizes[r] = sub_objects[r][0]->m();
2624  // then check that the following
2625  // block columns have the same
2626  // sizes
2627  for (unsigned int c=1; c<this->n_block_cols(); ++c)
2628  for (unsigned int r=0; r<this->n_block_rows(); ++r)
2629  Assert (row_sizes[r] == sub_objects[r][c]->m(),
2630  ExcIncompatibleRowNumbers (r,0,r,c));
2631 
2632  // finally initialize the row
2633  // indices with this array
2634  this->row_block_indices.reinit (row_sizes);
2635 
2636 
2637  // then do the same with the columns
2638  for (unsigned int c=0; c<this->n_block_cols(); ++c)
2639  col_sizes[c] = sub_objects[0][c]->n();
2640  for (unsigned int r=1; r<this->n_block_rows(); ++r)
2641  for (unsigned int c=0; c<this->n_block_cols(); ++c)
2642  Assert (col_sizes[c] == sub_objects[r][c]->n(),
2643  ExcIncompatibleRowNumbers (0,c,r,c));
2644 
2645  // finally initialize the row
2646  // indices with this array
2647  this->column_block_indices.reinit (col_sizes);
2648 }
2649 
2650 
2651 
2652 template <class MatrixType>
2653 void
2655 {
2656  for (unsigned int row=0; row<n_block_rows(); ++row)
2657  for (unsigned int col=0; col<n_block_cols(); ++col)
2658  block(row, col).prepare_add();
2659 }
2660 
2661 
2662 
2663 template <class MatrixType>
2664 void
2666 {
2667  for (unsigned int row=0; row<n_block_rows(); ++row)
2668  for (unsigned int col=0; col<n_block_cols(); ++col)
2669  block(row, col).prepare_set();
2670 }
2671 
2672 #endif // DOXYGEN
2673 
2674 
2675 DEAL_II_NAMESPACE_CLOSE
2676 
2677 #endif // dealii__block_matrix_base_h
Iterator lower_bound(Iterator first, Iterator last, const T &val)
Definition: utilities.h:604
const types::global_dof_index invalid_size_type
Definition: types.h:173
value_type matrix_scalar_product(const BlockVectorType &u, const BlockVectorType &v) const
static const unsigned int invalid_unsigned_int
Definition: types.h:164
void collect_sizes()
BlockMatrixBase & copy_from(const BlockMatrixType &source)
void Tvmult_nonblock_nonblock(VectorType &dst, const VectorType &src) const
void vmult_block_block(BlockVectorType &dst, const BlockVectorType &src) const
TemporaryData & operator=(const TemporaryData &)
std::vector< size_type > counter_within_block
size_type n() const
void vmult_nonblock_nonblock(VectorType &dst, const VectorType &src) const
std::vector< std::vector< size_type > > column_indices
void set(const size_type i, const size_type j, const value_type value)
size_type m() const
void print(std::ostream &out, const bool alternative_output=false) const
::ExceptionBase & ExcMessage(std::string arg1)
value_type residual(BlockVectorType &dst, const BlockVectorType &x, const BlockVectorType &b) const
Auxiliary class aiding in the handling of block structures like in BlockVector or FESystem...
Definition: block_indices.h:54
BlockMatrixType::value_type value_type
void Tvmult_block_nonblock(BlockVectorType &dst, const VectorType &src) const
void add(const size_type i, const size_type j, const value_type value)
unsigned int n_block_cols() const
void vmult_block_nonblock(BlockVectorType &dst, const VectorType &src) const
TemporaryData temporary_data
types::global_dof_index size_type
value_type el(const size_type i, const size_type j) const
BlockMatrixBase & operator*=(const value_type factor)
std::vector< std::vector< value_type > > column_values
BlockMatrixBase & operator/=(const value_type factor)
size_type n() const
void vmult_nonblock_block(VectorType &dst, const BlockVectorType &src) const
const BlockIndices & get_row_indices() const
unsigned int global_dof_index
Definition: types.h:88
#define Assert(cond, exc)
Definition: exceptions.h:294
void Tvmult_add(BlockVectorType &dst, const BlockVectorType &src) const
number value_type
const BlockIndices & get_column_indices() const
void prepare_add_operation()
iterator end()
BlockType::value_type value_type
size_type m() const
value_type operator()(const size_type i, const size_type j) const
value_type diag_element(const size_type i) const
unsigned int block_row() const
value_type matrix_norm_square(const BlockVectorType &v) const
void vmult_add(BlockVectorType &dst, const BlockVectorType &src) const
BlockMatrixType::BlockType::const_iterator base_iterator
DeclException4(ExcIncompatibleRowNumbers, int, int, int, int,<< "The blocks ["<< arg1<< ','<< arg2<< "] and ["<< arg3<< ','<< arg4<< "] have differing row numbers.")
std_cxx11::enable_if< std_cxx11::is_fundamental< T >::value, std::size_t >::type memory_consumption(const T &t)
std::size_t memory_consumption() const
void Tvmult_nonblock_block(VectorType &dst, const BlockVectorType &src) const
unsigned int block_column() const
void Tvmult_block_block(BlockVectorType &dst, const BlockVectorType &src) const
BlockType & block(const unsigned int row, const unsigned int column)
Table< 2, SmartPointer< BlockType, BlockMatrixBase< MatrixType > > > sub_objects
BlockIndices row_block_indices
Definition: table.h:33
#define AssertIsFinite(number)
Definition: exceptions.h:1096
void compress(::VectorOperation::values operation)
void prepare_set_operation()
unsigned int n_block_rows() const
iterator begin()