17 #ifndef __deal2__block_matrix_base_h
18 #define __deal2__block_matrix_base_h
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>
35 DEAL_II_NAMESPACE_OPEN
51 namespace BlockMatrixIterators
58 template <
class BlockMatrix>
120 template <
class BlockMatrix,
bool ConstNess>
128 template <
class BlockMatrix>
213 bool operator == (
const Accessor &a)
const;
216 friend class Accessor<BlockMatrix, true>;
224 template <
class BlockMatrix>
309 bool operator == (
const Accessor &a)
const;
316 friend class ::MatrixIterator;
381 template <
typename MatrixType>
448 template <
class BlockMatrixType>
450 copy_from (
const BlockMatrixType &source);
457 block (
const unsigned int row,
458 const unsigned int column);
467 block (
const unsigned int row,
468 const unsigned int column)
const;
476 size_type
m ()
const;
484 size_type
n ()
const;
512 void set (
const size_type i,
543 template <
typename number>
544 void set (
const std::vector<size_type> &indices,
546 const bool elide_zero_values =
false);
555 template <
typename number>
556 void set (
const std::vector<size_type> &row_indices,
557 const std::vector<size_type> &col_indices,
559 const bool elide_zero_values =
false);
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);
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);
617 void add (
const size_type i,
648 template <
typename number>
649 void add (
const std::vector<size_type> &indices,
651 const bool elide_zero_values =
true);
660 template <
typename number>
661 void add (
const std::vector<size_type> &row_indices,
662 const std::vector<size_type> &col_indices,
664 const bool elide_zero_values =
true);
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);
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);
741 const size_type j)
const;
761 const size_type j)
const;
791 void compress (::VectorOperation::values operation);
816 template <class BlockVectorType>
818 const BlockVectorType &src) const;
830 template <class BlockVectorType>
832 const BlockVectorType &src) const;
858 template <class BlockVectorType>
866 template <class BlockVectorType>
869 const BlockVectorType &v) const;
876 template <class BlockVectorType>
878 const BlockVectorType &x,
879 const BlockVectorType &b) const;
891 void print (std::ostream &out,
892 const
bool alternative_output = false) const;
967 << "The blocks [" << arg1 << ',' << arg2 << "] and ["
968 << arg3 << ',' << arg4 << "] have differing row numbers.");
974 << "The blocks [" << arg1 << ',' << arg2 << "] and ["
975 << arg3 << ',' << arg4 << "] have differing column numbers.");
1063 template <class BlockVectorType>
1065 const BlockVectorType &src) const;
1085 template <class BlockVectorType,
1108 template <class BlockVectorType,
1111 const BlockVectorType &src) const;
1154 template <class BlockVectorType>
1156 const BlockVectorType &src) const;
1176 template <class BlockVectorType,
1199 template <class BlockVectorType,
1202 const BlockVectorType &src) const;
1324 template <
typename,
bool>
1338 namespace BlockMatrixIterators
1340 template <
class BlockMatrix>
1349 template <
class BlockMatrix>
1352 AccessorBase<BlockMatrix>::block_row()
const
1361 template <
class BlockMatrix>
1364 AccessorBase<BlockMatrix>::block_column()
const
1373 template <
class BlockMatrix>
1375 Accessor<BlockMatrix, true>::Accessor (
1376 const BlockMatrix *matrix,
1381 base_iterator(matrix->block(0,0).begin())
1387 if (row < matrix->
m())
1389 const std::pair<unsigned int,size_type> indices
1390 = matrix->row_block_indices.global_to_local(row);
1394 for (
unsigned int bc=0; bc<matrix->n_block_cols(); ++bc)
1397 = matrix->block(indices.first, bc).begin(indices.second);
1398 if (base_iterator !=
1399 matrix->block(indices.first, bc).end(indices.second))
1401 this->row_block = indices.first;
1402 this->col_block = bc;
1413 *
this = Accessor (matrix, row+1, 0);
1437 template <
class BlockMatrix>
1439 Accessor<BlockMatrix, true>::Accessor (
const Accessor<BlockMatrix, false> &other)
1441 matrix(other.matrix),
1442 base_iterator(other.base_iterator)
1444 this->row_block = other.row_block;
1445 this->col_block = other.col_block;
1449 template <
class BlockMatrix>
1451 typename Accessor<BlockMatrix, true>::size_type
1452 Accessor<BlockMatrix, true>::row()
const
1457 return (matrix->row_block_indices.local_to_global(this->row_block, 0) +
1458 base_iterator->row());
1462 template <
class BlockMatrix>
1464 typename Accessor<BlockMatrix, true>::size_type
1465 Accessor<BlockMatrix, true>::column()
const
1470 return (matrix->column_block_indices.local_to_global(this->col_block,0) +
1471 base_iterator->column());
1475 template <
class BlockMatrix>
1477 typename Accessor<BlockMatrix, true>::value_type
1478 Accessor<BlockMatrix, true>::value ()
const
1485 return base_iterator->value();
1490 template <
class BlockMatrix>
1493 Accessor<BlockMatrix, true>::advance ()
1501 size_type local_row = base_iterator->row();
1512 while (base_iterator ==
1513 matrix->block(this->row_block, this->col_block).end(local_row))
1518 if (this->col_block < matrix->n_block_cols()-1)
1522 = matrix->block(this->row_block, this->col_block).begin(local_row);
1528 this->col_block = 0;
1537 if (local_row == matrix->block(this->row_block, this->col_block).m())
1541 if (this->row_block == matrix->n_block_rows())
1550 = matrix->block(this->row_block, this->col_block).begin(local_row);
1556 template <
class BlockMatrix>
1559 Accessor<BlockMatrix, true>::operator == (
const Accessor &a)
const
1561 if (matrix != a.matrix)
1564 if (this->row_block == a.row_block
1565 && this->col_block == a.col_block)
1574 (base_iterator == a.base_iterator));
1582 template <
class BlockMatrix>
1584 Accessor<BlockMatrix, false>::Accessor (
1585 BlockMatrix *matrix,
1590 base_iterator(matrix->block(0,0).begin())
1595 if (row < matrix->
m())
1597 const std::pair<unsigned int,size_type> indices
1598 = matrix->row_block_indices.global_to_local(row);
1602 for (
size_type bc=0; bc<matrix->n_block_cols(); ++bc)
1605 = matrix->block(indices.first, bc).begin(indices.second);
1606 if (base_iterator !=
1607 matrix->block(indices.first, bc).end(indices.second))
1609 this->row_block = indices.first;
1610 this->col_block = bc;
1621 *
this = Accessor (matrix, row+1, 0);
1633 template <
class BlockMatrix>
1635 typename Accessor<BlockMatrix, false>::size_type
1636 Accessor<BlockMatrix, false>::row()
const
1641 return (matrix->row_block_indices.local_to_global(this->row_block, 0) +
1642 base_iterator->row());
1646 template <
class BlockMatrix>
1648 typename Accessor<BlockMatrix, false>::size_type
1649 Accessor<BlockMatrix, false>::column()
const
1654 return (matrix->column_block_indices.local_to_global(this->col_block,0) +
1655 base_iterator->column());
1659 template <
class BlockMatrix>
1661 typename Accessor<BlockMatrix, false>::value_type
1662 Accessor<BlockMatrix, false>::value ()
const
1669 return base_iterator->value();
1674 template <
class BlockMatrix>
1677 Accessor<BlockMatrix, false>::set_value (
typename Accessor<BlockMatrix, false>::value_type newval)
const
1684 base_iterator->value() = newval;
1689 template <
class BlockMatrix>
1692 Accessor<BlockMatrix, false>::advance ()
1700 size_type local_row = base_iterator->row();
1711 while (base_iterator ==
1712 matrix->block(this->row_block, this->col_block).end(local_row))
1717 if (this->col_block < matrix->n_block_cols()-1)
1721 = matrix->block(this->row_block, this->col_block).begin(local_row);
1727 this->col_block = 0;
1736 if (local_row == matrix->block(this->row_block, this->col_block).m())
1740 if (this->row_block == matrix->n_block_rows())
1749 = matrix->block(this->row_block, this->col_block).begin(local_row);
1756 template <
class BlockMatrix>
1759 Accessor<BlockMatrix, false>::operator == (
const Accessor &a)
const
1761 if (matrix != a.matrix)
1764 if (this->row_block == a.row_block
1765 && this->col_block == a.col_block)
1774 (base_iterator == a.base_iterator));
1784 template <
typename MatrixType>
1789 template <
typename MatrixType>
1797 template <
class MatrixType>
1798 template <
class BlockMatrixType>
1802 copy_from (
const BlockMatrixType &source)
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));
1812 template <
class MatrixType>
1823 sizeof(temporary_data.mutex);
1825 for (
unsigned int r=0; r<n_block_rows(); ++r)
1826 for (
unsigned int c=0; c<n_block_cols(); ++c)
1828 MatrixType *p = this->sub_objects[r][c];
1837 template <
class MatrixType>
1842 for (
unsigned int r=0; r<n_block_rows(); ++r)
1843 for (
unsigned int c=0; c<n_block_cols(); ++c)
1845 MatrixType *p = this->sub_objects[r][c];
1846 this->sub_objects[r][c] = 0;
1849 sub_objects.reinit (0,0);
1852 row_block_indices = column_block_indices =
BlockIndices ();
1857 template <
class MatrixType>
1861 const unsigned int column)
1863 Assert (row<n_block_rows(),
1865 Assert (column<n_block_cols(),
1868 return *sub_objects[row][column];
1873 template <
class MatrixType>
1877 const unsigned int column)
const
1879 Assert (row<n_block_rows(),
1881 Assert (column<n_block_cols(),
1884 return *sub_objects[row][column];
1888 template <
class MatrixType>
1890 typename BlockMatrixBase<MatrixType>::size_type
1893 return row_block_indices.total_size();
1898 template <
class MatrixType>
1900 typename BlockMatrixBase<MatrixType>::size_type
1903 return column_block_indices.total_size();
1908 template <
class MatrixType>
1913 return column_block_indices.size();
1918 template <
class MatrixType>
1923 return row_block_indices.size();
1931 template <
class MatrixType>
1936 const value_type value)
1938 prepare_set_operation();
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,
1952 template <
class MatrixType>
1953 template <
typename number>
1957 const std::vector<size_type> &col_indices,
1959 const bool elide_zero_values)
1961 Assert (row_indices.size() == values.
m(),
1963 Assert (col_indices.size() == values.
n(),
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),
1973 template <
class MatrixType>
1974 template <
typename number>
1979 const bool elide_zero_values)
1981 Assert (indices.size() == values.
m(),
1983 Assert (values.
n() == values.
m(), ExcNotQuadratic());
1985 for (
size_type i=0; i<indices.size(); ++i)
1986 set (indices[i], indices.size(), &indices[0], &values(i,0),
1992 template <
class MatrixType>
1993 template <
typename number>
1997 const std::vector<size_type> &col_indices,
1998 const std::vector<number> &values,
1999 const bool elide_zero_values)
2001 Assert (col_indices.size() == values.size(),
2004 set (row, col_indices.size(), &col_indices[0], &values[0],
2013 template <
class MatrixType>
2014 template <
typename number>
2020 const number *values,
2021 const bool elide_zero_values)
2023 prepare_set_operation();
2030 if (temporary_data.column_indices.size() < this->n_block_cols())
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());
2047 if (temporary_data.column_indices[0].size() < n_cols)
2049 for (
unsigned int i=0; i<this->n_block_cols(); ++i)
2051 temporary_data.column_indices[i].resize(n_cols);
2052 temporary_data.column_values[i].resize(n_cols);
2058 for (
unsigned int i=0; i<this->n_block_cols(); ++i)
2059 temporary_data.counter_within_block[i] = 0;
2072 double value = values[j];
2074 if (value == 0 && elide_zero_values ==
true)
2077 const std::pair<unsigned int, size_type>
2078 col_index = this->column_block_indices.global_to_local(col_indices[j]);
2080 const size_type local_index = temporary_data.counter_within_block[col_index.first]++;
2082 temporary_data.column_indices[col_index.first][local_index] = col_index.second;
2083 temporary_data.column_values[col_index.first][local_index] = value;
2090 for (
unsigned int i=0; i<this->n_block_cols(); ++i)
2091 length += temporary_data.counter_within_block[i];
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)
2104 if (temporary_data.counter_within_block[block_col] == 0)
2107 block(row_index.first, block_col).set
2109 temporary_data.counter_within_block[block_col],
2110 &temporary_data.column_indices[block_col][0],
2111 &temporary_data.column_values[block_col][0],
2118 template <
class MatrixType>
2123 const value_type value)
2128 prepare_add_operation();
2133 typedef typename MatrixType::Traits MatrixTraits;
2134 if ((MatrixTraits::zero_addition_can_be_elided ==
true)
2136 (value == value_type()))
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,
2149 template <
class MatrixType>
2150 template <
typename number>
2154 const std::vector<size_type> &col_indices,
2156 const bool elide_zero_values)
2158 Assert (row_indices.size() == values.
m(),
2160 Assert (col_indices.size() == values.
n(),
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),
2170 template <
class MatrixType>
2171 template <
typename number>
2176 const bool elide_zero_values)
2178 Assert (indices.size() == values.
m(),
2180 Assert (values.
n() == values.
m(), ExcNotQuadratic());
2182 for (
size_type i=0; i<indices.size(); ++i)
2183 add (indices[i], indices.size(), &indices[0], &values(i,0),
2189 template <
class MatrixType>
2190 template <
typename number>
2194 const std::vector<size_type> &col_indices,
2195 const std::vector<number> &values,
2196 const bool elide_zero_values)
2198 Assert (col_indices.size() == values.size(),
2201 add (row, col_indices.size(), &col_indices[0], &values[0],
2210 template <
class MatrixType>
2211 template <
typename number>
2217 const number *values,
2218 const bool elide_zero_values,
2219 const bool col_indices_are_sorted)
2221 prepare_add_operation();
2226 if (col_indices_are_sorted ==
true)
2233 if (col_indices[i] <= before)
2235 "indices appear to not be sorted."))
2237 before = col_indices[i];
2239 const std::pair<unsigned int,size_type>
2240 row_index = this->row_block_indices.global_to_local (row);
2242 if (this->n_block_cols() > 1)
2246 this->column_block_indices.block_start(1));
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,
2254 col_indices_are_sorted);
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,
2263 block(row_index.first, 0). add (row_index.second,
2268 col_indices_are_sorted);
2277 if (temporary_data.column_indices.size() < this->n_block_cols())
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());
2294 if (temporary_data.column_indices[0].size() < n_cols)
2296 for (
unsigned int i=0; i<this->n_block_cols(); ++i)
2298 temporary_data.column_indices[i].resize(n_cols);
2299 temporary_data.column_values[i].resize(n_cols);
2305 for (
unsigned int i=0; i<this->n_block_cols(); ++i)
2306 temporary_data.counter_within_block[i] = 0;
2319 double value = values[j];
2321 if (value == 0 && elide_zero_values ==
true)
2324 const std::pair<unsigned int, size_type>
2325 col_index = this->column_block_indices.global_to_local(col_indices[j]);
2327 const size_type local_index = temporary_data.counter_within_block[col_index.first]++;
2329 temporary_data.column_indices[col_index.first][local_index] = col_index.second;
2330 temporary_data.column_values[col_index.first][local_index] = value;
2337 for (
unsigned int i=0; i<this->n_block_cols(); ++i)
2338 length += temporary_data.counter_within_block[i];
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)
2351 if (temporary_data.counter_within_block[block_col] == 0)
2354 block(row_index.first, block_col).add
2356 temporary_data.counter_within_block[block_col],
2357 &temporary_data.column_indices[block_col][0],
2358 &temporary_data.column_values[block_col][0],
2360 col_indices_are_sorted);
2366 template <
class MatrixType>
2374 prepare_add_operation();
2379 typedef typename MatrixType::Traits MatrixTraits;
2380 if ((MatrixTraits::zero_addition_can_be_elided ==
true)
2385 for (
unsigned int row=0; row<n_block_rows(); ++row)
2386 for (
unsigned int col=0; col<n_block_cols(); ++col)
2389 block(row, col).add(factor, matrix.
block(row,col));
2394 template <
class MatrixType>
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,
2409 template <
class MatrixType>
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,
2424 template <
class MatrixType>
2429 Assert (n_block_rows() == n_block_cols(),
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);
2439 template <
class MatrixType>
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);
2449 template <
class MatrixType>
2454 compress(::VectorOperation::unknown);
2459 template <
class MatrixType>
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;
2476 template <
class MatrixType>
2485 const value_type factor_inv = 1. / factor;
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;
2496 template <
class MatrixType>
2500 return this->row_block_indices;
2505 template <
class MatrixType>
2509 return this->column_block_indices;
2514 template <
class MatrixType>
2515 template <
class BlockVectorType>
2519 const BlockVectorType &src)
const
2521 Assert (dst.n_blocks() == n_block_rows(),
2523 Assert (src.n_blocks() == n_block_cols(),
2526 for (
size_type row=0; row<n_block_rows(); ++row)
2528 block(row,0).vmult (dst.block(row),
2530 for (
size_type col=1; col<n_block_cols(); ++col)
2531 block(row,col).vmult_add (dst.block(row),
2538 template <
class MatrixType>
2539 template <
class BlockVectorType,
2544 const BlockVectorType &src)
const
2546 Assert (n_block_rows() == 1,
2548 Assert (src.n_blocks() == n_block_cols(),
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));
2558 template <
class MatrixType>
2559 template <
class BlockVectorType,
2566 Assert (dst.n_blocks() == n_block_rows(),
2568 Assert (1 == n_block_cols(),
2571 for (
size_type row=0; row<n_block_rows(); ++row)
2572 block(row,0).vmult (dst.block(row),
2578 template <
class MatrixType>
2579 template <
class VectorType>
2585 Assert (1 == n_block_rows(),
2587 Assert (1 == n_block_cols(),
2590 block(0,0).vmult (dst, src);
2595 template <
class MatrixType>
2596 template <
class BlockVectorType>
2599 const BlockVectorType &src)
const
2601 Assert (dst.n_blocks() == n_block_rows(),
2603 Assert (src.n_blocks() == n_block_cols(),
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),
2615 template <
class MatrixType>
2616 template <
class BlockVectorType>
2620 const BlockVectorType &src)
const
2622 Assert (dst.n_blocks() == n_block_cols(),
2624 Assert (src.n_blocks() == n_block_rows(),
2629 for (
unsigned int row=0; row<n_block_rows(); ++row)
2631 for (
unsigned int col=0; col<n_block_cols(); ++col)
2632 block(row,col).Tvmult_add (dst.block(col),
2639 template <
class MatrixType>
2640 template <
class BlockVectorType,
2647 Assert (dst.n_blocks() == n_block_cols(),
2649 Assert (1 == n_block_rows(),
2654 for (
unsigned int col=0; col<n_block_cols(); ++col)
2655 block(0,col).Tvmult_add (dst.block(col), src);
2660 template <
class MatrixType>
2661 template <
class BlockVectorType,
2666 const BlockVectorType &src)
const
2668 Assert (1 == n_block_cols(),
2670 Assert (src.n_blocks() == n_block_rows(),
2673 block(0,0).Tvmult (dst, src.block(0));
2675 for (
size_type row=1; row<n_block_rows(); ++row)
2676 block(row,0).Tvmult_add (dst, src.block(row));
2681 template <
class MatrixType>
2682 template <
class VectorType>
2688 Assert (1 == n_block_cols(),
2690 Assert (1 == n_block_rows(),
2693 block(0,0).Tvmult (dst, src);
2698 template <
class MatrixType>
2699 template <
class BlockVectorType>
2702 const BlockVectorType &src)
const
2704 Assert (dst.n_blocks() == n_block_cols(),
2706 Assert (src.n_blocks() == n_block_rows(),
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),
2717 template <
class MatrixType>
2718 template <
class BlockVectorType>
2722 Assert (n_block_rows() == n_block_cols(), ExcNotQuadratic());
2723 Assert (v.n_blocks() == n_block_rows(),
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)
2730 norm_sqr += block(row,col).matrix_norm_square (v.block(row));
2732 norm_sqr += block(row,col).matrix_scalar_product (v.block(row),
2739 template <
class MatrixType>
2740 template <
class BlockVectorType>
2744 const BlockVectorType &v)
const
2746 Assert (u.n_blocks() == n_block_rows(),
2748 Assert (v.n_blocks() == n_block_cols(),
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),
2761 template <
class MatrixType>
2762 template <
class BlockVectorType>
2766 const BlockVectorType &x,
2767 const BlockVectorType &b)
const
2769 Assert (dst.n_blocks() == n_block_rows(),
2771 Assert (b.n_blocks() == n_block_rows(),
2773 Assert (x.n_blocks() == n_block_cols(),
2789 for (
unsigned int row=0; row<n_block_rows(); ++row)
2791 block(row,0).residual (dst.block(row),
2795 for (
size_type i=0; i<dst.block(row).size(); ++i)
2796 dst.block(row)(i) = -dst.block(row)(i);
2798 for (
unsigned int col=1; col<n_block_cols(); ++col)
2799 block(row,col).vmult_add (dst.block(row),
2802 for (
size_type i=0; i<dst.block(row).size(); ++i)
2803 dst.block(row)(i) = -dst.block(row)(i);
2807 for (
size_type row=0; row<n_block_rows(); ++row)
2808 res += dst.block(row).norm_sqr ();
2809 return std::sqrt(res);
2814 template <
class MatrixType>
2818 const bool alternative_output)
const
2820 for (
unsigned int row=0; row<n_block_rows(); ++row)
2821 for (
unsigned int col=0; col<n_block_cols(); ++col)
2823 if (!alternative_output)
2824 out <<
"Block (" << row <<
", " << col <<
")" << std::endl;
2826 block(row, col).print(out, alternative_output);
2832 template <
class MatrixType>
2837 return const_iterator(
this, 0);
2842 template <
class MatrixType>
2847 return const_iterator(
this, m());
2852 template <
class MatrixType>
2858 return const_iterator(
this, r);
2863 template <
class MatrixType>
2869 return const_iterator(
this, r+1);
2874 template <
class MatrixType>
2879 return iterator(
this, 0);
2884 template <
class MatrixType>
2889 return iterator(
this, m());
2894 template <
class MatrixType>
2900 return iterator(
this, r);
2905 template <
class MatrixType>
2911 return iterator(
this, r+1);
2916 template <
class MatrixType>
2920 std::vector<size_type> row_sizes (this->n_block_rows());
2921 std::vector<size_type> col_sizes (this->n_block_cols());
2925 for (
unsigned int r=0; r<this->n_block_rows(); ++r)
2926 row_sizes[r] = sub_objects[r][0]->m();
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));
2937 this->row_block_indices.reinit (row_sizes);
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));
2950 this->column_block_indices.reinit (col_sizes);
2955 template <
class MatrixType>
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();
2966 template <
class MatrixType>
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();
2978 DEAL_II_NAMESPACE_CLOSE
2980 #endif // __deal2__block_matrix_base_h
Iterator lower_bound(Iterator first, Iterator last, const T &val)
const types::global_dof_index invalid_size_type
static const unsigned int invalid_unsigned_int
const BlockMatrix * matrix
const BlockMatrix MatrixType
void add(const size_type i, const size_type j, const value_type value)
void vmult_block_block(BlockVectorType &dst, const BlockVectorType &src) const
void Tvmult_add(BlockVectorType &dst, const BlockVectorType &src) const
unsigned int n_block_cols() const
BlockMatrix::BlockType::const_iterator base_iterator
unsigned int n_block_rows() const
Subscriptor & operator=(const Subscriptor &)
types::global_dof_index size_type
std::size_t memory_consumption() const
::ExceptionBase & ExcMessage(std::string arg1)
BlockIndices row_block_indices
Auxiliary class aiding in the handling of block structures like in BlockVector or FESystem...
unsigned int block_row() const
value_type diag_element(const size_type i) const
void compress() DEAL_II_DEPRECATED
void Tvmult_nonblock_nonblock(VectorType &dst, const VectorType &src) const
value_type matrix_scalar_product(const BlockVectorType &u, const BlockVectorType &v) const
void Tvmult_block_block(BlockVectorType &dst, const BlockVectorType &src) const
void vmult_nonblock_block(VectorType &dst, const BlockVectorType &src) const
std::vector< std::vector< double > > column_values
value_type matrix_norm_square(const BlockVectorType &v) const
bool is_finite(const double x)
void vmult_add(BlockVectorType &dst, const BlockVectorType &src) const
BlockMatrix::BlockType::iterator base_iterator
void Tvmult_nonblock_block(VectorType &dst, const BlockVectorType &src) const
void set(const size_type i, const size_type j, const value_type value)
value_type residual(BlockVectorType &dst, const BlockVectorType &x, const BlockVectorType &b) const
BlockType::value_type value_type
const BlockIndices & get_row_indices() const
unsigned int global_dof_index
BlockMatrix::value_type value_type
#define Assert(cond, exc)
std::size_t memory_consumption(const T &t)
DeclException4(ExcIncompatibleRowNumbers, int, int, int, int,<< "The blocks ["<< arg1<< ','<< arg2<< "] and ["<< arg3<< ','<< arg4<< "] have differing row numbers.")
types::global_dof_index size_type
void prepare_set_operation()
BlockMatrixBase & copy_from(const BlockMatrixType &source)
void vmult_nonblock_nonblock(VectorType &dst, const VectorType &src) const
void prepare_add_operation()
BlockMatrix::value_type value_type
::ExceptionBase & ExcIndexRange(int arg1, int arg2, int arg3)
Table< 2, SmartPointer< BlockType, BlockMatrixBase< MatrixType > > > sub_objects
void vmult_block_nonblock(BlockVectorType &dst, const VectorType &src) const
BlockCompressedSparsityPattern CompressedBlockSparsityPattern DEAL_II_DEPRECATED
value_type operator()(const size_type i, const size_type j) const
TemporaryData temporary_data
unsigned int block_column() const
const BlockIndices & get_column_indices() const
BlockMatrixBase & operator*=(const value_type factor)
::ExceptionBase & ExcIteratorPastEnd()
std::vector< size_type > counter_within_block
::ExceptionBase & ExcNumberNotFinite()
BlockMatrix::value_type value_type
void Tvmult_block_nonblock(BlockVectorType &dst, const VectorType &src) const
::ExceptionBase & ExcNotImplemented()
::ExceptionBase & ExcNotInitialized()
::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
void print(std::ostream &out, const bool alternative_output=false) const
types::global_dof_index size_type
BlockMatrixBase & operator/=(const value_type factor)
value_type el(const size_type i, const size_type j) const
std::vector< std::vector< size_type > > column_indices
::ExceptionBase & ExcInternalError()
::ExceptionBase & ExcDivideByZero()
BlockType & block(const unsigned int row, const unsigned int column)