15#ifndef dealii_block_matrix_base_h
16#define dealii_block_matrix_base_h
61 template <
typename BlockMatrixType>
109 template <
typename BlockMatrixType,
bool Constness>
116 template <
typename BlockMatrixType>
193 friend class ::MatrixIterator;
195 friend class Accessor<BlockMatrixType, true>;
204 template <
typename BlockMatrixType>
283 friend class ::MatrixIterator;
348template <
typename MatrixType>
402 template <
typename BlockMatrixType>
410 block(
const unsigned int row,
const unsigned int column);
418 block(
const unsigned int row,
const unsigned int column)
const;
472 template <
typename number>
474 set(
const std::vector<size_type> &indices,
476 const bool elide_zero_values =
false);
483 template <
typename number>
485 set(
const std::vector<size_type> &row_indices,
486 const std::vector<size_type> &col_indices,
488 const bool elide_zero_values =
false);
500 template <
typename number>
503 const std::vector<size_type> &col_indices,
504 const std::vector<number> &values,
505 const bool elide_zero_values =
false);
516 template <
typename number>
521 const number *values,
522 const bool elide_zero_values =
false);
546 template <
typename number>
548 add(
const std::vector<size_type> &indices,
550 const bool elide_zero_values =
true);
557 template <
typename number>
559 add(
const std::vector<size_type> &row_indices,
560 const std::vector<size_type> &col_indices,
562 const bool elide_zero_values =
true);
573 template <
typename number>
576 const std::vector<size_type> &col_indices,
577 const std::vector<number> &values,
578 const bool elide_zero_values =
true);
589 template <
typename number>
594 const number *values,
595 const bool elide_zero_values =
true,
596 const bool col_indices_are_sorted =
false);
672 template <
typename BlockVectorType>
674 vmult_add(BlockVectorType &dst,
const BlockVectorType &src)
const;
681 template <
typename BlockVectorType>
683 Tvmult_add(BlockVectorType &dst,
const BlockVectorType &src)
const;
697 template <
typename BlockVectorType>
711 template <
typename BlockVectorType>
714 const BlockVectorType &v)
const;
719 template <
typename BlockVectorType>
722 const BlockVectorType &x,
723 const BlockVectorType &b)
const;
732 print(std::ostream &out,
const bool alternative_output =
false)
const;
814 <<
"The blocks [" << arg1 <<
',' << arg2 <<
"] and [" << arg3
815 <<
',' << arg4 <<
"] have differing row numbers.");
824 <<
"The blocks [" << arg1 <<
',' << arg2 <<
"] and [" << arg3
825 <<
',' << arg4 <<
"] have differing column numbers.");
885 template <
typename BlockVectorType>
899 template <
typename BlockVectorType,
typename VectorType>
913 template <
typename BlockVectorType,
typename VectorType>
927 template <
typename VectorType>
942 template <
typename BlockVectorType>
956 template <
typename BlockVectorType,
typename VectorType>
970 template <
typename BlockVectorType,
typename VectorType>
984 template <
typename VectorType>
1068 template <
typename,
bool>
1084 template <
typename BlockMatrixType>
1091 template <
typename BlockMatrixType>
1093 AccessorBase<BlockMatrixType>::block_row()
const
1101 template <
typename BlockMatrixType>
1103 AccessorBase<BlockMatrixType>::block_column()
const
1111 template <
typename BlockMatrixType>
1112 inline Accessor<BlockMatrixType, true>::Accessor(
1113 const BlockMatrixType *matrix,
1114 const size_type row,
1115 const size_type col)
1124 if (row < matrix->
m())
1126 const std::pair<unsigned int, size_type> indices =
1127 matrix->row_block_indices.global_to_local(row);
1131 for (
unsigned int bc = 0; bc <
matrix->n_block_cols(); ++bc)
1134 matrix->block(indices.first, bc).begin(indices.second);
1135 if (base_iterator !=
1136 matrix->block(indices.first, bc).end(indices.second))
1138 this->row_block = indices.first;
1139 this->col_block = bc;
1150 *
this =
Accessor(matrix, row + 1, 0);
1175 template <
typename BlockMatrixType>
1176 inline Accessor<BlockMatrixType, true>::Accessor(
1177 const Accessor<BlockMatrixType, false> &other)
1179 , base_iterator(other.base_iterator)
1181 this->row_block = other.row_block;
1182 this->col_block = other.col_block;
1186 template <
typename BlockMatrixType>
1187 inline typename Accessor<BlockMatrixType, true>::size_type
1188 Accessor<BlockMatrixType, true>::row()
const
1193 return (
matrix->row_block_indices.local_to_global(this->row_block, 0) +
1194 base_iterator->row());
1198 template <
typename BlockMatrixType>
1199 inline typename Accessor<BlockMatrixType, true>::size_type
1200 Accessor<BlockMatrixType, true>::column()
const
1205 return (
matrix->column_block_indices.local_to_global(this->col_block, 0) +
1206 base_iterator->column());
1210 template <
typename BlockMatrixType>
1211 inline typename Accessor<BlockMatrixType, true>::value_type
1212 Accessor<BlockMatrixType, true>::value()
const
1219 return base_iterator->value();
1224 template <
typename BlockMatrixType>
1226 Accessor<BlockMatrixType, true>::advance()
1234 size_type local_row = base_iterator->row();
1245 while (base_iterator ==
1246 matrix->block(this->row_block, this->col_block).end(local_row))
1251 if (this->col_block < matrix->n_block_cols() - 1)
1255 matrix->block(this->row_block, this->col_block).begin(local_row);
1261 this->col_block = 0;
1271 matrix->block(this->row_block, this->col_block).m())
1275 if (this->row_block ==
matrix->n_block_rows())
1284 matrix->block(this->row_block, this->col_block).begin(local_row);
1290 template <
typename BlockMatrixType>
1292 Accessor<BlockMatrixType, true>::operator==(
const Accessor &a)
const
1294 if (matrix != a.matrix)
1297 if (this->row_block == a.row_block && this->col_block == a.col_block)
1304 (base_iterator == a.base_iterator));
1312 template <
typename BlockMatrixType>
1313 inline Accessor<BlockMatrixType, false>::Accessor(BlockMatrixType *matrix,
1314 const size_type row,
1315 const size_type col)
1323 if (row < matrix->
m())
1325 const std::pair<unsigned int, size_type> indices =
1326 matrix->row_block_indices.global_to_local(row);
1333 matrix->block(indices.first, bc).begin(indices.second);
1334 if (base_iterator !=
1335 matrix->block(indices.first, bc).end(indices.second))
1337 this->row_block = indices.first;
1338 this->col_block = bc;
1349 *
this =
Accessor(matrix, row + 1, 0);
1361 template <
typename BlockMatrixType>
1362 inline typename Accessor<BlockMatrixType, false>::size_type
1363 Accessor<BlockMatrixType, false>::row()
const
1367 return (
matrix->row_block_indices.local_to_global(this->row_block, 0) +
1368 base_iterator->row());
1372 template <
typename BlockMatrixType>
1373 inline typename Accessor<BlockMatrixType, false>::size_type
1374 Accessor<BlockMatrixType, false>::column()
const
1378 return (
matrix->column_block_indices.local_to_global(this->col_block, 0) +
1379 base_iterator->column());
1383 template <
typename BlockMatrixType>
1384 inline typename Accessor<BlockMatrixType, false>::value_type
1385 Accessor<BlockMatrixType, false>::value()
const
1390 return base_iterator->value();
1395 template <
typename BlockMatrixType>
1397 Accessor<BlockMatrixType, false>::set_value(
1398 typename Accessor<BlockMatrixType, false>::value_type newval)
const
1403 base_iterator->value() = newval;
1408 template <
typename BlockMatrixType>
1410 Accessor<BlockMatrixType, false>::advance()
1416 size_type local_row = base_iterator->row();
1427 while (base_iterator ==
1428 matrix->block(this->row_block, this->col_block).end(local_row))
1433 if (this->col_block < matrix->n_block_cols() - 1)
1437 matrix->block(this->row_block, this->col_block).begin(local_row);
1443 this->col_block = 0;
1453 matrix->block(this->row_block, this->col_block).m())
1457 if (this->row_block ==
matrix->n_block_rows())
1466 matrix->block(this->row_block, this->col_block).begin(local_row);
1473 template <
typename BlockMatrixType>
1475 Accessor<BlockMatrixType, false>::operator==(
const Accessor &a)
const
1477 if (matrix != a.matrix)
1480 if (this->row_block == a.row_block && this->col_block == a.col_block)
1487 (base_iterator == a.base_iterator));
1496template <
typename MatrixType>
1508template <
typename MatrixType>
1509template <
typename BlockMatrixType>
1513 for (
unsigned int r = 0; r < n_block_rows(); ++r)
1514 for (
unsigned int c = 0; c < n_block_cols(); ++c)
1515 block(r, c).
copy_from(source.block(r, c));
1521template <
typename MatrixType>
1532 sizeof(temporary_data.mutex);
1534 for (
unsigned int r = 0; r < n_block_rows(); ++r)
1535 for (
unsigned int c = 0; c < n_block_cols(); ++c)
1546template <
typename MatrixType>
1550 for (
unsigned int r = 0; r < n_block_rows(); ++r)
1551 for (
unsigned int c = 0; c < n_block_cols(); ++c)
1554 this->sub_objects[r][c] =
nullptr;
1557 sub_objects.reinit(0, 0);
1560 row_block_indices = column_block_indices =
BlockIndices();
1565template <
typename MatrixType>
1568 const unsigned int column)
1573 return *sub_objects[row][column];
1578template <
typename MatrixType>
1581 const unsigned int column)
const
1586 return *sub_objects[row][column];
1590template <
typename MatrixType>
1594 return row_block_indices.total_size();
1599template <
typename MatrixType>
1603 return column_block_indices.total_size();
1608template <
typename MatrixType>
1612 return column_block_indices.size();
1617template <
typename MatrixType>
1621 return row_block_indices.size();
1629template <
typename MatrixType>
1633 const value_type value)
1635 prepare_set_operation();
1639 const std::pair<unsigned int, size_type>
1640 row_index = row_block_indices.global_to_local(i),
1641 col_index = column_block_indices.global_to_local(j);
1642 block(row_index.first, col_index.first)
1643 .set(row_index.second, col_index.second, value);
1648template <
typename MatrixType>
1649template <
typename number>
1652 const std::vector<size_type> &col_indices,
1654 const bool elide_zero_values)
1661 for (size_type i = 0; i < row_indices.size(); ++i)
1671template <
typename MatrixType>
1672template <
typename number>
1676 const bool elide_zero_values)
1682 for (size_type i = 0; i < indices.size(); ++i)
1692template <
typename MatrixType>
1693template <
typename number>
1696 const std::vector<size_type> &col_indices,
1697 const std::vector<number> &values,
1698 const bool elide_zero_values)
1715template <
typename MatrixType>
1716template <
typename number>
1719 const size_type n_cols,
1720 const size_type *col_indices,
1721 const number *values,
1722 const bool elide_zero_values)
1724 prepare_set_operation();
1728 std::lock_guard<std::mutex> lock(temporary_data.mutex);
1731 if (temporary_data.column_indices.size() < this->n_block_cols())
1733 temporary_data.column_indices.resize(this->n_block_cols());
1734 temporary_data.column_values.resize(this->n_block_cols());
1735 temporary_data.counter_within_block.resize(this->n_block_cols());
1748 if (temporary_data.column_indices[0].size() < n_cols)
1750 for (
unsigned int i = 0; i < this->n_block_cols(); ++i)
1752 temporary_data.column_indices[i].resize(n_cols);
1753 temporary_data.column_values[i].resize(n_cols);
1759 for (
unsigned int i = 0; i < this->n_block_cols(); ++i)
1760 temporary_data.counter_within_block[i] = 0;
1771 for (size_type j = 0; j < n_cols; ++j)
1775 if (value == number() && elide_zero_values ==
true)
1778 const std::pair<unsigned int, size_type> col_index =
1779 this->column_block_indices.global_to_local(col_indices[j]);
1782 temporary_data.counter_within_block[col_index.first]++;
1784 temporary_data.column_indices[col_index.first][local_index] =
1786 temporary_data.column_values[col_index.first][local_index] =
value;
1794 for (
unsigned int i = 0; i < this->n_block_cols(); ++i)
1795 length += temporary_data.counter_within_block[i];
1804 const std::pair<unsigned int, size_type> row_index =
1805 this->row_block_indices.global_to_local(row);
1806 for (
unsigned int block_col = 0; block_col < n_block_cols(); ++block_col)
1808 if (temporary_data.counter_within_block[block_col] == 0)
1811 block(row_index.first, block_col)
1812 .set(row_index.second,
1813 temporary_data.counter_within_block[block_col],
1814 temporary_data.column_indices[block_col].data(),
1815 temporary_data.column_values[block_col].data(),
1822template <
typename MatrixType>
1826 const value_type value)
1830 prepare_add_operation();
1835 using MatrixTraits =
typename MatrixType::Traits;
1836 if ((MatrixTraits::zero_addition_can_be_elided ==
true) &&
1837 (value == value_type()))
1840 const std::pair<unsigned int, size_type>
1841 row_index = row_block_indices.global_to_local(i),
1842 col_index = column_block_indices.global_to_local(j);
1843 block(row_index.first, col_index.first)
1844 .add(row_index.second, col_index.second, value);
1849template <
typename MatrixType>
1850template <
typename number>
1853 const std::vector<size_type> &col_indices,
1855 const bool elide_zero_values)
1862 for (size_type i = 0; i < row_indices.size(); ++i)
1872template <
typename MatrixType>
1873template <
typename number>
1877 const bool elide_zero_values)
1883 for (size_type i = 0; i < indices.size(); ++i)
1893template <
typename MatrixType>
1894template <
typename number>
1897 const std::vector<size_type> &col_indices,
1898 const std::vector<number> &values,
1899 const bool elide_zero_values)
1916template <
typename MatrixType>
1917template <
typename number>
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)
1926 prepare_add_operation();
1931 if (col_indices_are_sorted ==
true)
1938 for (size_type i = 1; i < n_cols; ++i)
1939 if (col_indices[i] <= before)
1942 ExcMessage(
"Flag col_indices_are_sorted is set, but "
1943 "indices appear to not be sorted."));
1946 before = col_indices[i];
1948 const std::pair<unsigned int, size_type> row_index =
1949 this->row_block_indices.global_to_local(row);
1951 if (this->n_block_cols() > 1)
1955 col_indices + n_cols,
1956 this->column_block_indices.block_start(1));
1958 const size_type n_zero_block_indices = first_block - col_indices;
1959 block(row_index.first, 0)
1960 .add(row_index.second,
1961 n_zero_block_indices,
1965 col_indices_are_sorted);
1967 if (n_zero_block_indices < n_cols)
1969 n_cols - n_zero_block_indices,
1971 values + n_zero_block_indices,
1977 block(row_index.first, 0)
1978 .add(row_index.second,
1983 col_indices_are_sorted);
1990 std::lock_guard<std::mutex> lock(temporary_data.mutex);
1992 if (temporary_data.column_indices.size() < this->n_block_cols())
1994 temporary_data.column_indices.resize(this->n_block_cols());
1995 temporary_data.column_values.resize(this->n_block_cols());
1996 temporary_data.counter_within_block.resize(this->n_block_cols());
2009 if (temporary_data.column_indices[0].size() < n_cols)
2011 for (
unsigned int i = 0; i < this->n_block_cols(); ++i)
2013 temporary_data.column_indices[i].resize(n_cols);
2014 temporary_data.column_values[i].resize(n_cols);
2020 for (
unsigned int i = 0; i < this->n_block_cols(); ++i)
2021 temporary_data.counter_within_block[i] = 0;
2032 for (size_type j = 0; j < n_cols; ++j)
2036 if (value == number() && elide_zero_values ==
true)
2039 const std::pair<unsigned int, size_type> col_index =
2040 this->column_block_indices.global_to_local(col_indices[j]);
2043 temporary_data.counter_within_block[col_index.first]++;
2045 temporary_data.column_indices[col_index.first][local_index] =
2047 temporary_data.column_values[col_index.first][local_index] =
value;
2055 for (
unsigned int i = 0; i < this->n_block_cols(); ++i)
2056 length += temporary_data.counter_within_block[i];
2065 const std::pair<unsigned int, size_type> row_index =
2066 this->row_block_indices.global_to_local(row);
2067 for (
unsigned int block_col = 0; block_col < n_block_cols(); ++block_col)
2069 if (temporary_data.counter_within_block[block_col] == 0)
2072 block(row_index.first, block_col)
2073 .add(row_index.second,
2074 temporary_data.counter_within_block[block_col],
2075 temporary_data.column_indices[block_col].data(),
2076 temporary_data.column_values[block_col].data(),
2078 col_indices_are_sorted);
2084template <
typename MatrixType>
2091 prepare_add_operation();
2096 using MatrixTraits =
typename MatrixType::Traits;
2097 if ((MatrixTraits::zero_addition_can_be_elided ==
true) && (factor == 0))
2100 for (
unsigned int row = 0; row < n_block_rows(); ++row)
2101 for (
unsigned int col = 0; col < n_block_cols(); ++col)
2104 block(row, col).add(factor,
matrix.block(row, col));
2109template <
typename MatrixType>
2112 const size_type j)
const
2114 const std::pair<unsigned int, size_type>
2115 row_index = row_block_indices.global_to_local(i),
2116 col_index = column_block_indices.global_to_local(j);
2117 return block(row_index.first, col_index.first)(row_index.second,
2123template <
typename MatrixType>
2127 const std::pair<unsigned int, size_type>
2128 row_index = row_block_indices.global_to_local(i),
2129 col_index = column_block_indices.global_to_local(j);
2130 return block(row_index.first, col_index.first)
2131 .el(row_index.second, col_index.second);
2136template <
typename MatrixType>
2142 const std::pair<unsigned int, size_type>
index =
2143 row_block_indices.global_to_local(i);
2149template <
typename MatrixType>
2153 for (
unsigned int r = 0; r < n_block_rows(); ++r)
2154 for (
unsigned int c = 0; c < n_block_cols(); ++c)
2155 block(r, c).compress(operation);
2160template <
typename MatrixType>
2167 for (
unsigned int r = 0; r < n_block_rows(); ++r)
2168 for (
unsigned int c = 0; c < n_block_cols(); ++c)
2169 block(r, c) *= factor;
2176template <
typename MatrixType>
2184 const value_type factor_inv = 1. / factor;
2186 for (
unsigned int r = 0; r < n_block_rows(); ++r)
2187 for (
unsigned int c = 0; c < n_block_cols(); ++c)
2188 block(r, c) *= factor_inv;
2195template <
typename MatrixType>
2199 return this->row_block_indices;
2204template <
typename MatrixType>
2208 return this->column_block_indices;
2213template <
typename MatrixType>
2214template <
typename BlockVectorType>
2217 const BlockVectorType &src)
const
2219 Assert(dst.n_blocks() == n_block_rows(),
2221 Assert(src.n_blocks() == n_block_cols(),
2224 for (size_type row = 0; row < n_block_rows(); ++row)
2226 block(row, 0).vmult(dst.block(row), src.block(0));
2227 for (size_type col = 1; col < n_block_cols(); ++col)
2228 block(row, col).vmult_add(dst.block(row), src.block(col));
2234template <
typename MatrixType>
2235template <
typename BlockVectorType,
typename VectorType>
2239 const BlockVectorType &src)
const
2242 Assert(src.n_blocks() == n_block_cols(),
2245 block(0, 0).vmult(dst, src.block(0));
2246 for (size_type col = 1; col < n_block_cols(); ++col)
2247 block(0, col).vmult_add(dst, src.block(col));
2252template <
typename MatrixType>
2253template <
typename BlockVectorType,
typename VectorType>
2256 const VectorType &src)
const
2258 Assert(dst.n_blocks() == n_block_rows(),
2262 for (size_type row = 0; row < n_block_rows(); ++row)
2263 block(row, 0).vmult(dst.block(row), src);
2268template <
typename MatrixType>
2269template <
typename VectorType>
2273 const VectorType &src)
const
2278 block(0, 0).vmult(dst, src);
2283template <
typename MatrixType>
2284template <
typename BlockVectorType>
2287 const BlockVectorType &src)
const
2289 Assert(dst.n_blocks() == n_block_rows(),
2291 Assert(src.n_blocks() == n_block_cols(),
2294 for (
unsigned int row = 0; row < n_block_rows(); ++row)
2295 for (
unsigned int col = 0; col < n_block_cols(); ++col)
2296 block(row, col).vmult_add(dst.block(row), src.block(col));
2301template <
typename MatrixType>
2302template <
typename BlockVectorType>
2305 BlockVectorType &dst,
2306 const BlockVectorType &src)
const
2308 Assert(dst.n_blocks() == n_block_cols(),
2310 Assert(src.n_blocks() == n_block_rows(),
2315 for (
unsigned int row = 0; row < n_block_rows(); ++row)
2317 for (
unsigned int col = 0; col < n_block_cols(); ++col)
2318 block(row, col).Tvmult_add(dst.block(col), src.block(row));
2324template <
typename MatrixType>
2325template <
typename BlockVectorType,
typename VectorType>
2328 const VectorType &src)
const
2330 Assert(dst.n_blocks() == n_block_cols(),
2336 for (
unsigned int col = 0; col < n_block_cols(); ++col)
2337 block(0, col).Tvmult_add(dst.block(col), src);
2342template <
typename MatrixType>
2343template <
typename BlockVectorType,
typename VectorType>
2347 const BlockVectorType &src)
const
2350 Assert(src.n_blocks() == n_block_rows(),
2353 block(0, 0).Tvmult(dst, src.block(0));
2355 for (size_type row = 1; row < n_block_rows(); ++row)
2356 block(row, 0).Tvmult_add(dst, src.block(row));
2361template <
typename MatrixType>
2362template <
typename VectorType>
2366 const VectorType &src)
const
2371 block(0, 0).Tvmult(dst, src);
2376template <
typename MatrixType>
2377template <
typename BlockVectorType>
2380 const BlockVectorType &src)
const
2382 Assert(dst.n_blocks() == n_block_cols(),
2384 Assert(src.n_blocks() == n_block_rows(),
2387 for (
unsigned int row = 0; row < n_block_rows(); ++row)
2388 for (
unsigned int col = 0; col < n_block_cols(); ++col)
2389 block(row, col).Tvmult_add(dst.block(col), src.block(row));
2394template <
typename MatrixType>
2395template <
typename BlockVectorType>
2400 Assert(v.n_blocks() == n_block_rows(),
2403 value_type norm_sqr = 0;
2404 for (
unsigned int row = 0; row < n_block_rows(); ++row)
2405 for (
unsigned int col = 0; col < n_block_cols(); ++col)
2407 norm_sqr += block(row, col).matrix_norm_square(v.block(row));
2410 block(row, col).matrix_scalar_product(v.block(row), v.block(col));
2416template <
typename MatrixType>
2420 value_type norm_sqr = 0;
2424 for (
unsigned int row = 0; row < n_block_rows(); ++row)
2426 for (
unsigned int col = 0; col < n_block_cols(); ++col)
2428 const value_type block_norm = block(row, col).frobenius_norm();
2429 norm_sqr += block_norm * block_norm;
2438template <
typename MatrixType>
2439template <
typename BlockVectorType>
2442 const BlockVectorType &u,
2443 const BlockVectorType &v)
const
2445 Assert(u.n_blocks() == n_block_rows(),
2447 Assert(v.n_blocks() == n_block_cols(),
2450 value_type result = 0;
2451 for (
unsigned int row = 0; row < n_block_rows(); ++row)
2452 for (
unsigned int col = 0; col < n_block_cols(); ++col)
2454 block(row, col).matrix_scalar_product(u.block(row), v.block(col));
2460template <
typename MatrixType>
2461template <
typename BlockVectorType>
2464 const BlockVectorType &x,
2465 const BlockVectorType &b)
const
2467 Assert(dst.n_blocks() == n_block_rows(),
2469 Assert(
b.n_blocks() == n_block_rows(),
2471 Assert(x.n_blocks() == n_block_cols(),
2487 for (
unsigned int row = 0; row < n_block_rows(); ++row)
2489 block(row, 0).residual(dst.block(row), x.block(0),
b.block(row));
2491 for (size_type i = 0; i < dst.block(row).
size(); ++i)
2492 dst.block(row)(i) = -dst.block(row)(i);
2494 for (
unsigned int col = 1; col < n_block_cols(); ++col)
2495 block(row, col).vmult_add(dst.block(row), x.block(col));
2497 for (size_type i = 0; i < dst.block(row).size(); ++i)
2498 dst.block(row)(i) = -dst.block(row)(i);
2502 for (size_type row = 0; row < n_block_rows(); ++row)
2503 res += dst.block(row).norm_sqr();
2509template <
typename MatrixType>
2512 const bool alternative_output)
const
2514 for (
unsigned int row = 0; row < n_block_rows(); ++row)
2515 for (
unsigned int col = 0; col < n_block_cols(); ++col)
2517 if (!alternative_output)
2518 out <<
"Block (" << row <<
", " << col <<
')' << std::endl;
2520 block(row, col).print(out, alternative_output);
2526template <
typename MatrixType>
2530 return const_iterator(
this, 0);
2535template <
typename MatrixType>
2539 return const_iterator(
this, m());
2544template <
typename MatrixType>
2549 return const_iterator(
this, r);
2554template <
typename MatrixType>
2559 return const_iterator(
this, r + 1);
2564template <
typename MatrixType>
2568 return iterator(
this, 0);
2573template <
typename MatrixType>
2577 return iterator(
this, m());
2582template <
typename MatrixType>
2587 return iterator(
this, r);
2592template <
typename MatrixType>
2597 return iterator(
this, r + 1);
2602template <
typename MatrixType>
2606 std::vector<size_type> row_sizes(this->n_block_rows());
2607 std::vector<size_type> col_sizes(this->n_block_cols());
2611 for (
unsigned int r = 0; r < this->n_block_rows(); ++r)
2612 row_sizes[r] = sub_objects[r][0]->m();
2616 for (
unsigned int c = 1; c < this->n_block_cols(); ++c)
2617 for (
unsigned int r = 0; r < this->n_block_rows(); ++r)
2618 Assert(row_sizes[r] == sub_objects[r][c]->m(),
2619 ExcIncompatibleRowNumbers(r, 0, r, c));
2623 this->row_block_indices.reinit(row_sizes);
2627 for (
unsigned int c = 0; c < this->n_block_cols(); ++c)
2628 col_sizes[c] = sub_objects[0][c]->n();
2629 for (
unsigned int r = 1; r < this->n_block_rows(); ++r)
2630 for (
unsigned int c = 0; c < this->n_block_cols(); ++c)
2631 Assert(col_sizes[c] == sub_objects[r][c]->n(),
2632 ExcIncompatibleRowNumbers(0, c, r, c));
2636 this->column_block_indices.reinit(col_sizes);
2641template <
typename MatrixType>
2645 for (
unsigned int row = 0; row < n_block_rows(); ++row)
2646 for (
unsigned int col = 0; col < n_block_cols(); ++col)
2647 block(row, col).prepare_add();
2652template <
typename MatrixType>
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_set();
void Tvmult_nonblock_nonblock(VectorType &dst, const VectorType &src) const
void vmult_nonblock_nonblock(VectorType &dst, const VectorType &src) const
void add(const size_type row, const std::vector< size_type > &col_indices, const std::vector< number > &values, const bool elide_zero_values=true)
void add(const std::vector< size_type > &row_indices, const std::vector< size_type > &col_indices, const FullMatrix< number > &full_matrix, const bool elide_zero_values=true)
void add(const size_type row, const size_type n_cols, const size_type *col_indices, const number *values, const bool elide_zero_values=true, const bool col_indices_are_sorted=false)
void print(std::ostream &out, const bool alternative_output=false) const
friend class BlockMatrixIterators::Accessor
const_iterator end() const
BlockIndices column_block_indices
BlockMatrixBase & operator*=(const value_type factor)
value_type matrix_norm_square(const BlockVectorType &v) const
void Tvmult_block_block(BlockVectorType &dst, const BlockVectorType &src) const
const value_type & const_reference
void Tvmult_nonblock_block(VectorType &dst, const BlockVectorType &src) const
unsigned int n_block_rows() const
value_type operator()(const size_type i, const size_type j) const
value_type el(const size_type i, const size_type j) const
const_iterator begin() const
real_type frobenius_norm() const
void vmult_block_block(BlockVectorType &dst, const BlockVectorType &src) const
void compress(VectorOperation::values operation)
void vmult_add(BlockVectorType &dst, const BlockVectorType &src) const
types::global_dof_index size_type
void Tvmult_add(BlockVectorType &dst, const BlockVectorType &src) const
typename BlockType::value_type value_type
void prepare_add_operation()
void vmult_nonblock_block(VectorType &dst, const BlockVectorType &src) const
value_type residual(BlockVectorType &dst, const BlockVectorType &x, const BlockVectorType &b) const
BlockMatrixBase()=default
void set(const size_type row, const size_type n_cols, const size_type *col_indices, const number *values, const bool elide_zero_values=false)
void prepare_set_operation()
unsigned int n_block_cols() const
const value_type * const_pointer
TemporaryData temporary_data
const BlockIndices & get_row_indices() const
void set(const size_type i, const size_type j, const value_type value)
BlockMatrixBase & copy_from(const BlockMatrixType &source)
std::size_t memory_consumption() const
void add(const value_type factor, const BlockMatrixBase< MatrixType > &matrix)
const BlockType & block(const unsigned int row, const unsigned int column) const
typename numbers::NumberTraits< value_type >::real_type real_type
BlockIndices row_block_indices
~BlockMatrixBase() override
BlockType & block(const unsigned int row, const unsigned int column)
value_type matrix_scalar_product(const BlockVectorType &u, const BlockVectorType &v) const
const BlockIndices & get_column_indices() const
void set(const std::vector< size_type > &indices, const FullMatrix< number > &full_matrix, const bool elide_zero_values=false)
void add(const std::vector< size_type > &indices, const FullMatrix< number > &full_matrix, const bool elide_zero_values=true)
void add(const size_type i, const size_type j, const value_type value)
void set(const size_type row, const std::vector< size_type > &col_indices, const std::vector< number > &values, const bool elide_zero_values=false)
value_type diag_element(const size_type i) const
void vmult_block_nonblock(BlockVectorType &dst, const VectorType &src) const
Table< 2, ObserverPointer< BlockType, BlockMatrixBase< MatrixType > > > sub_objects
iterator begin(const size_type r)
const_iterator begin(const size_type r) const
BlockMatrixBase & operator/=(const value_type factor)
void Tvmult_block_nonblock(BlockVectorType &dst, const VectorType &src) const
iterator end(const size_type r)
const_iterator end(const size_type r) const
void set(const std::vector< size_type > &row_indices, const std::vector< size_type > &col_indices, const FullMatrix< number > &full_matrix, const bool elide_zero_values=false)
unsigned int block_row() const
unsigned int block_column() const
typename BlockMatrixType::value_type value_type
Accessor(BlockMatrixType *m, const size_type row, const size_type col)
BlockMatrixType MatrixType
BlockMatrixType::BlockType::iterator base_iterator
bool operator==(const Accessor &a) const
void set_value(value_type newval) const
typename BlockMatrixType::value_type value_type
const BlockMatrixType MatrixType
typename BlockMatrixType::value_type value_type
bool operator==(const Accessor &a) const
Accessor(const Accessor< BlockMatrixType, false > &)
const BlockMatrixType * matrix
BlockMatrixType::BlockType::const_iterator base_iterator
Accessor(const BlockMatrixType *m, const size_type row, const size_type col)
#define DEAL_II_NAMESPACE_OPEN
constexpr bool running_in_debug_mode()
#define DEAL_II_NAMESPACE_CLOSE
#define DeclException4(Exception4, type1, type2, type3, type4, outsequence)
static ::ExceptionBase & ExcIncompatibleColNumbers(int arg1, int arg2, int arg3, int arg4)
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
static ::ExceptionBase & ExcIncompatibleRowNumbers(int arg1, int arg2, int arg3, int arg4)
static ::ExceptionBase & ExcIteratorPastEnd()
#define AssertIsFinite(number)
#define AssertIndexRange(index, range)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcDivideByZero()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
static ::ExceptionBase & ExcNotInitialized()
static ::ExceptionBase & ExcNotQuadratic()
static ::ExceptionBase & ExcMessage(std::string arg1)
types::global_dof_index size_type
@ matrix
Contents is actually a matrix.
Tpetra::CrsMatrix< Number, LO, GO, NodeType< MemorySpace > > MatrixType
std::enable_if_t< std::is_fundamental_v< T >, std::size_t > memory_consumption(const T &t)
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
VectorType::value_type * begin(VectorType &V)
Iterator lower_bound(Iterator first, Iterator last, const T &val)
constexpr types::global_dof_index invalid_size_type
constexpr unsigned int invalid_unsigned_int
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)
unsigned int global_dof_index
std::vector< std::vector< size_type > > column_indices
std::vector< std::vector< value_type > > column_values
TemporaryData & operator=(const TemporaryData &)
std::vector< size_type > counter_within_block