16 #ifndef dealii_block_matrix_base_h
17 #define dealii_block_matrix_base_h
62 template <
typename BlockMatrixType>
110 template <
typename BlockMatrixType,
bool Constness>
117 template <
typename BlockMatrixType>
194 friend class ::MatrixIterator;
196 friend class Accessor<BlockMatrixType, true>;
205 template <
typename BlockMatrixType>
284 friend class ::MatrixIterator;
349 template <
typename MatrixType>
403 template <
typename BlockMatrixType>
411 block(
const unsigned int row,
const unsigned int column);
419 block(
const unsigned int row,
const unsigned int column)
const;
473 template <
typename number>
475 set(
const std::vector<size_type> &indices,
477 const bool elide_zero_values =
false);
484 template <
typename number>
486 set(
const std::vector<size_type> &row_indices,
487 const std::vector<size_type> &col_indices,
489 const bool elide_zero_values =
false);
501 template <
typename number>
504 const std::vector<size_type> &col_indices,
505 const std::vector<number> &
values,
506 const bool elide_zero_values =
false);
517 template <
typename number>
523 const bool elide_zero_values =
false);
547 template <
typename number>
549 add(
const std::vector<size_type> &indices,
551 const bool elide_zero_values =
true);
558 template <
typename number>
560 add(
const std::vector<size_type> &row_indices,
561 const std::vector<size_type> &col_indices,
563 const bool elide_zero_values =
true);
574 template <
typename number>
577 const std::vector<size_type> &col_indices,
578 const std::vector<number> &
values,
579 const bool elide_zero_values =
true);
590 template <
typename number>
596 const bool elide_zero_values =
true,
597 const bool col_indices_are_sorted =
false);
673 template <
typename BlockVectorType>
675 vmult_add(BlockVectorType &dst,
const BlockVectorType &src)
const;
682 template <
typename BlockVectorType>
684 Tvmult_add(BlockVectorType &dst,
const BlockVectorType &src)
const;
698 template <
typename BlockVectorType>
712 template <
typename BlockVectorType>
715 const BlockVectorType &v)
const;
720 template <
typename BlockVectorType>
723 const BlockVectorType &x,
724 const BlockVectorType &
b)
const;
733 print(std::ostream &out,
const bool alternative_output =
false)
const;
815 <<
"The blocks [" << arg1 <<
',' << arg2 <<
"] and [" << arg3
816 <<
',' << arg4 <<
"] have differing row numbers.");
825 <<
"The blocks [" << arg1 <<
',' << arg2 <<
"] and [" << arg3
826 <<
',' << arg4 <<
"] have differing column numbers.");
886 template <
typename BlockVectorType>
900 template <
typename BlockVectorType,
typename VectorType>
914 template <
typename BlockVectorType,
typename VectorType>
928 template <
typename VectorType>
943 template <
typename BlockVectorType>
957 template <
typename BlockVectorType,
typename VectorType>
971 template <
typename BlockVectorType,
typename VectorType>
985 template <
typename VectorType>
1069 template <
typename,
bool>
1085 template <
typename BlockMatrixType>
1092 template <
typename BlockMatrixType>
1094 AccessorBase<BlockMatrixType>::block_row()
const
1102 template <
typename BlockMatrixType>
1104 AccessorBase<BlockMatrixType>::block_column()
const
1112 template <
typename BlockMatrixType>
1113 inline Accessor<BlockMatrixType, true>::Accessor(
1114 const BlockMatrixType *
matrix,
1125 if (row < matrix->
m())
1127 const std::pair<unsigned int, size_type> indices =
1128 matrix->row_block_indices.global_to_local(row);
1132 for (
unsigned int bc = 0; bc <
matrix->n_block_cols(); ++bc)
1135 matrix->block(indices.first, bc).begin(indices.second);
1136 if (base_iterator !=
1137 matrix->block(indices.first, bc).end(indices.second))
1139 this->row_block = indices.first;
1140 this->col_block = bc;
1151 *
this = Accessor(
matrix, row + 1, 0);
1176 template <
typename BlockMatrixType>
1177 inline Accessor<BlockMatrixType, true>::Accessor(
1178 const Accessor<BlockMatrixType, false> &other)
1180 , base_iterator(other.base_iterator)
1182 this->row_block = other.row_block;
1183 this->col_block = other.col_block;
1187 template <
typename BlockMatrixType>
1189 Accessor<BlockMatrixType, true>::row()
const
1194 return (
matrix->row_block_indices.local_to_global(this->row_block, 0) +
1195 base_iterator->row());
1199 template <
typename BlockMatrixType>
1201 Accessor<BlockMatrixType, true>::column()
const
1206 return (
matrix->column_block_indices.local_to_global(this->col_block, 0) +
1207 base_iterator->column());
1211 template <
typename BlockMatrixType>
1212 inline typename Accessor<BlockMatrixType, true>::value_type
1213 Accessor<BlockMatrixType, true>::value()
const
1220 return base_iterator->value();
1225 template <
typename BlockMatrixType>
1235 size_type local_row = base_iterator->row();
1246 while (base_iterator ==
1247 matrix->block(this->row_block, this->col_block).end(local_row))
1252 if (this->col_block < matrix->n_block_cols() - 1)
1256 matrix->block(this->row_block, this->col_block).begin(local_row);
1262 this->col_block = 0;
1272 matrix->block(this->row_block, this->col_block).m())
1276 if (this->row_block ==
matrix->n_block_rows())
1285 matrix->block(this->row_block, this->col_block).begin(local_row);
1291 template <
typename BlockMatrixType>
1298 if (this->row_block == a.row_block && this->col_block == a.col_block)
1305 (base_iterator == a.base_iterator));
1313 template <
typename BlockMatrixType>
1314 inline Accessor<BlockMatrixType, false>::Accessor(BlockMatrixType *
matrix,
1324 if (row < matrix->
m())
1326 const std::pair<unsigned int, size_type> indices =
1327 matrix->row_block_indices.global_to_local(row);
1334 matrix->block(indices.first, bc).begin(indices.second);
1335 if (base_iterator !=
1336 matrix->block(indices.first, bc).end(indices.second))
1338 this->row_block = indices.first;
1339 this->col_block = bc;
1350 *
this = Accessor(
matrix, row + 1, 0);
1362 template <
typename BlockMatrixType>
1364 Accessor<BlockMatrixType, false>::row()
const
1368 return (
matrix->row_block_indices.local_to_global(this->row_block, 0) +
1369 base_iterator->row());
1373 template <
typename BlockMatrixType>
1375 Accessor<BlockMatrixType, false>::column()
const
1379 return (
matrix->column_block_indices.local_to_global(this->col_block, 0) +
1380 base_iterator->column());
1384 template <
typename BlockMatrixType>
1385 inline typename Accessor<BlockMatrixType, false>::value_type
1386 Accessor<BlockMatrixType, false>::value()
const
1391 return base_iterator->value();
1396 template <
typename BlockMatrixType>
1398 Accessor<BlockMatrixType, false>::set_value(
1399 typename Accessor<BlockMatrixType, false>::value_type newval)
const
1404 base_iterator->value() = newval;
1409 template <
typename BlockMatrixType>
1417 size_type local_row = base_iterator->row();
1428 while (base_iterator ==
1429 matrix->block(this->row_block, this->col_block).end(local_row))
1434 if (this->col_block < matrix->n_block_cols() - 1)
1438 matrix->block(this->row_block, this->col_block).begin(local_row);
1444 this->col_block = 0;
1454 matrix->block(this->row_block, this->col_block).m())
1458 if (this->row_block ==
matrix->n_block_rows())
1467 matrix->block(this->row_block, this->col_block).begin(local_row);
1474 template <
typename BlockMatrixType>
1481 if (this->row_block == a.row_block && this->col_block == a.col_block)
1488 (base_iterator == a.base_iterator));
1497 template <
typename MatrixType>
1509 template <
typename MatrixType>
1510 template <
typename BlockMatrixType>
1514 for (
unsigned int r = 0; r < n_block_rows(); ++r)
1515 for (
unsigned int c = 0; c < n_block_cols(); ++c)
1516 block(r, c).
copy_from(source.block(r, c));
1522 template <
typename MatrixType>
1533 sizeof(temporary_data.mutex);
1535 for (
unsigned int r = 0; r < n_block_rows(); ++r)
1536 for (
unsigned int c = 0; c < n_block_cols(); ++c)
1538 MatrixType *p = this->sub_objects[r][c];
1547 template <
typename MatrixType>
1551 for (
unsigned int r = 0; r < n_block_rows(); ++r)
1552 for (
unsigned int c = 0; c < n_block_cols(); ++c)
1554 MatrixType *p = this->sub_objects[r][c];
1555 this->sub_objects[r][c] =
nullptr;
1558 sub_objects.reinit(0, 0);
1561 row_block_indices = column_block_indices =
BlockIndices();
1566 template <
typename MatrixType>
1569 const unsigned int column)
1574 return *sub_objects[row][column];
1579 template <
typename MatrixType>
1582 const unsigned int column)
const
1587 return *sub_objects[row][column];
1591 template <
typename MatrixType>
1595 return row_block_indices.total_size();
1600 template <
typename MatrixType>
1604 return column_block_indices.total_size();
1609 template <
typename MatrixType>
1613 return column_block_indices.size();
1618 template <
typename MatrixType>
1622 return row_block_indices.size();
1630 template <
typename MatrixType>
1634 const value_type value)
1636 prepare_set_operation();
1640 const std::pair<unsigned int, size_type>
1641 row_index = row_block_indices.global_to_local(i),
1642 col_index = column_block_indices.global_to_local(j);
1643 block(row_index.first, col_index.first)
1644 .set(row_index.second, col_index.second, value);
1649 template <
typename MatrixType>
1650 template <
typename number>
1653 const std::vector<size_type> &col_indices,
1655 const bool elide_zero_values)
1662 for (
size_type i = 0; i < row_indices.size(); ++i)
1672 template <
typename MatrixType>
1673 template <
typename number>
1677 const bool elide_zero_values)
1683 for (
size_type i = 0; i < indices.size(); ++i)
1693 template <
typename MatrixType>
1694 template <
typename number>
1697 const std::vector<size_type> &col_indices,
1698 const std::vector<number> &
values,
1699 const bool elide_zero_values)
1716 template <
typename MatrixType>
1717 template <
typename number>
1723 const bool elide_zero_values)
1725 prepare_set_operation();
1729 std::lock_guard<std::mutex> lock(temporary_data.mutex);
1732 if (temporary_data.column_indices.size() < this->n_block_cols())
1734 temporary_data.column_indices.resize(this->n_block_cols());
1735 temporary_data.column_values.resize(this->n_block_cols());
1736 temporary_data.counter_within_block.resize(this->n_block_cols());
1749 if (temporary_data.column_indices[0].size() < n_cols)
1751 for (
unsigned int i = 0; i < this->n_block_cols(); ++i)
1753 temporary_data.column_indices[i].resize(n_cols);
1754 temporary_data.column_values[i].resize(n_cols);
1760 for (
unsigned int i = 0; i < this->n_block_cols(); ++i)
1761 temporary_data.counter_within_block[i] = 0;
1776 if (value == number() && elide_zero_values ==
true)
1779 const std::pair<unsigned int, size_type> col_index =
1780 this->column_block_indices.global_to_local(col_indices[j]);
1783 temporary_data.counter_within_block[col_index.first]++;
1785 temporary_data.column_indices[col_index.first][local_index] =
1787 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(),
1822 template <
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);
1849 template <
typename MatrixType>
1850 template <
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)
1872 template <
typename MatrixType>
1873 template <
typename number>
1877 const bool elide_zero_values)
1883 for (
size_type i = 0; i < indices.size(); ++i)
1893 template <
typename MatrixType>
1894 template <
typename number>
1897 const std::vector<size_type> &col_indices,
1898 const std::vector<number> &
values,
1899 const bool elide_zero_values)
1916 template <
typename MatrixType>
1917 template <
typename number>
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 if (col_indices[i] <= before)
1941 ExcMessage(
"Flag col_indices_are_sorted is set, but "
1942 "indices appear to not be sorted."));
1945 before = col_indices[i];
1947 const std::pair<unsigned int, size_type> row_index =
1948 this->row_block_indices.global_to_local(row);
1950 if (this->n_block_cols() > 1)
1954 col_indices + n_cols,
1955 this->column_block_indices.block_start(1));
1957 const size_type n_zero_block_indices = first_block - col_indices;
1958 block(row_index.first, 0)
1959 .add(row_index.second,
1960 n_zero_block_indices,
1964 col_indices_are_sorted);
1966 if (n_zero_block_indices < n_cols)
1968 n_cols - n_zero_block_indices,
1970 values + n_zero_block_indices,
1976 block(row_index.first, 0)
1977 .add(row_index.second,
1982 col_indices_are_sorted);
1989 std::lock_guard<std::mutex> lock(temporary_data.mutex);
1991 if (temporary_data.column_indices.size() < this->n_block_cols())
1993 temporary_data.column_indices.resize(this->n_block_cols());
1994 temporary_data.column_values.resize(this->n_block_cols());
1995 temporary_data.counter_within_block.resize(this->n_block_cols());
2008 if (temporary_data.column_indices[0].size() < n_cols)
2010 for (
unsigned int i = 0; i < this->n_block_cols(); ++i)
2012 temporary_data.column_indices[i].resize(n_cols);
2013 temporary_data.column_values[i].resize(n_cols);
2019 for (
unsigned int i = 0; i < this->n_block_cols(); ++i)
2020 temporary_data.counter_within_block[i] = 0;
2035 if (value == number() && elide_zero_values ==
true)
2038 const std::pair<unsigned int, size_type> col_index =
2039 this->column_block_indices.global_to_local(col_indices[j]);
2042 temporary_data.counter_within_block[col_index.first]++;
2044 temporary_data.column_indices[col_index.first][local_index] =
2046 temporary_data.column_values[col_index.first][local_index] =
value;
2053 for (
unsigned int i = 0; i < this->n_block_cols(); ++i)
2054 length += temporary_data.counter_within_block[i];
2063 const std::pair<unsigned int, size_type> row_index =
2064 this->row_block_indices.global_to_local(row);
2065 for (
unsigned int block_col = 0; block_col < n_block_cols(); ++block_col)
2067 if (temporary_data.counter_within_block[block_col] == 0)
2070 block(row_index.first, block_col)
2071 .add(row_index.second,
2072 temporary_data.counter_within_block[block_col],
2073 temporary_data.column_indices[block_col].data(),
2074 temporary_data.column_values[block_col].data(),
2076 col_indices_are_sorted);
2082 template <
typename MatrixType>
2089 prepare_add_operation();
2094 using MatrixTraits =
typename MatrixType::Traits;
2095 if ((MatrixTraits::zero_addition_can_be_elided ==
true) && (factor == 0))
2098 for (
unsigned int row = 0; row < n_block_rows(); ++row)
2099 for (
unsigned int col = 0; col < n_block_cols(); ++col)
2102 block(row, col).add(factor,
matrix.block(row, col));
2107 template <
typename MatrixType>
2112 const std::pair<unsigned int, size_type>
2113 row_index = row_block_indices.global_to_local(i),
2114 col_index = column_block_indices.global_to_local(j);
2115 return block(row_index.first, col_index.first)(row_index.second,
2121 template <
typename MatrixType>
2125 const std::pair<unsigned int, size_type>
2126 row_index = row_block_indices.global_to_local(i),
2127 col_index = column_block_indices.global_to_local(j);
2128 return block(row_index.first, col_index.first)
2129 .el(row_index.second, col_index.second);
2134 template <
typename MatrixType>
2140 const std::pair<unsigned int, size_type>
index =
2141 row_block_indices.global_to_local(i);
2147 template <
typename MatrixType>
2151 for (
unsigned int r = 0; r < n_block_rows(); ++r)
2152 for (
unsigned int c = 0; c < n_block_cols(); ++c)
2153 block(r, c).compress(operation);
2158 template <
typename MatrixType>
2165 for (
unsigned int r = 0; r < n_block_rows(); ++r)
2166 for (
unsigned int c = 0; c < n_block_cols(); ++c)
2167 block(r, c) *= factor;
2174 template <
typename MatrixType>
2182 const value_type factor_inv = 1. / factor;
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;
2193 template <
typename MatrixType>
2197 return this->row_block_indices;
2202 template <
typename MatrixType>
2206 return this->column_block_indices;
2211 template <
typename MatrixType>
2212 template <
typename BlockVectorType>
2215 const BlockVectorType &src)
const
2217 Assert(dst.n_blocks() == n_block_rows(),
2219 Assert(src.n_blocks() == n_block_cols(),
2222 for (
size_type row = 0; row < n_block_rows(); ++row)
2224 block(row, 0).vmult(dst.block(row), src.block(0));
2225 for (
size_type col = 1; col < n_block_cols(); ++col)
2226 block(row, col).vmult_add(dst.block(row), src.block(col));
2232 template <
typename MatrixType>
2233 template <
typename BlockVectorType,
typename VectorType>
2237 const BlockVectorType &src)
const
2240 Assert(src.n_blocks() == n_block_cols(),
2243 block(0, 0).vmult(dst, src.block(0));
2244 for (
size_type col = 1; col < n_block_cols(); ++col)
2245 block(0, col).vmult_add(dst, src.block(col));
2250 template <
typename MatrixType>
2251 template <
typename BlockVectorType,
typename VectorType>
2254 const VectorType &src)
const
2256 Assert(dst.n_blocks() == n_block_rows(),
2260 for (
size_type row = 0; row < n_block_rows(); ++row)
2261 block(row, 0).vmult(dst.block(row), src);
2266 template <
typename MatrixType>
2267 template <
typename VectorType>
2271 const VectorType &src)
const
2276 block(0, 0).vmult(dst, src);
2281 template <
typename MatrixType>
2282 template <
typename BlockVectorType>
2285 const BlockVectorType &src)
const
2287 Assert(dst.n_blocks() == n_block_rows(),
2289 Assert(src.n_blocks() == n_block_cols(),
2292 for (
unsigned int row = 0; row < n_block_rows(); ++row)
2293 for (
unsigned int col = 0; col < n_block_cols(); ++col)
2294 block(row, col).vmult_add(dst.block(row), src.block(col));
2299 template <
typename MatrixType>
2300 template <
typename BlockVectorType>
2303 BlockVectorType &dst,
2304 const BlockVectorType &src)
const
2306 Assert(dst.n_blocks() == n_block_cols(),
2308 Assert(src.n_blocks() == n_block_rows(),
2313 for (
unsigned int row = 0; row < n_block_rows(); ++row)
2315 for (
unsigned int col = 0; col < n_block_cols(); ++col)
2316 block(row, col).Tvmult_add(dst.block(col), src.block(row));
2322 template <
typename MatrixType>
2323 template <
typename BlockVectorType,
typename VectorType>
2326 const VectorType &src)
const
2328 Assert(dst.n_blocks() == n_block_cols(),
2334 for (
unsigned int col = 0; col < n_block_cols(); ++col)
2335 block(0, col).Tvmult_add(dst.block(col), src);
2340 template <
typename MatrixType>
2341 template <
typename BlockVectorType,
typename VectorType>
2345 const BlockVectorType &src)
const
2348 Assert(src.n_blocks() == n_block_rows(),
2351 block(0, 0).Tvmult(dst, src.block(0));
2353 for (
size_type row = 1; row < n_block_rows(); ++row)
2354 block(row, 0).Tvmult_add(dst, src.block(row));
2359 template <
typename MatrixType>
2360 template <
typename VectorType>
2364 const VectorType &src)
const
2369 block(0, 0).Tvmult(dst, src);
2374 template <
typename MatrixType>
2375 template <
typename BlockVectorType>
2378 const BlockVectorType &src)
const
2380 Assert(dst.n_blocks() == n_block_cols(),
2382 Assert(src.n_blocks() == n_block_rows(),
2385 for (
unsigned int row = 0; row < n_block_rows(); ++row)
2386 for (
unsigned int col = 0; col < n_block_cols(); ++col)
2387 block(row, col).Tvmult_add(dst.block(col), src.block(row));
2392 template <
typename MatrixType>
2393 template <
typename BlockVectorType>
2398 Assert(v.n_blocks() == n_block_rows(),
2401 value_type norm_sqr = 0;
2402 for (
unsigned int row = 0; row < n_block_rows(); ++row)
2403 for (
unsigned int col = 0; col < n_block_cols(); ++col)
2405 norm_sqr += block(row, col).matrix_norm_square(v.block(row));
2408 block(row, col).matrix_scalar_product(v.block(row), v.block(col));
2414 template <
typename MatrixType>
2418 value_type norm_sqr = 0;
2422 for (
unsigned int row = 0; row < n_block_rows(); ++row)
2424 for (
unsigned int col = 0; col < n_block_cols(); ++col)
2426 const value_type block_norm = block(row, col).frobenius_norm();
2427 norm_sqr += block_norm * block_norm;
2436 template <
typename MatrixType>
2437 template <
typename BlockVectorType>
2440 const BlockVectorType &u,
2441 const BlockVectorType &v)
const
2443 Assert(u.n_blocks() == n_block_rows(),
2445 Assert(v.n_blocks() == n_block_cols(),
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)
2452 block(row, col).matrix_scalar_product(u.block(row), v.block(col));
2458 template <
typename MatrixType>
2459 template <
typename BlockVectorType>
2462 const BlockVectorType &x,
2463 const BlockVectorType &
b)
const
2465 Assert(dst.n_blocks() == n_block_rows(),
2467 Assert(
b.n_blocks() == n_block_rows(),
2469 Assert(x.n_blocks() == n_block_cols(),
2485 for (
unsigned int row = 0; row < n_block_rows(); ++row)
2487 block(row, 0).residual(dst.block(row), x.block(0),
b.block(row));
2489 for (
size_type i = 0; i < dst.block(row).size(); ++i)
2490 dst.block(row)(i) = -dst.block(row)(i);
2492 for (
unsigned int col = 1; col < n_block_cols(); ++col)
2493 block(row, col).vmult_add(dst.block(row), x.block(col));
2495 for (
size_type i = 0; i < dst.block(row).size(); ++i)
2496 dst.block(row)(i) = -dst.block(row)(i);
2500 for (
size_type row = 0; row < n_block_rows(); ++row)
2501 res += dst.block(row).norm_sqr();
2507 template <
typename MatrixType>
2510 const bool alternative_output)
const
2512 for (
unsigned int row = 0; row < n_block_rows(); ++row)
2513 for (
unsigned int col = 0; col < n_block_cols(); ++col)
2515 if (!alternative_output)
2516 out <<
"Block (" << row <<
", " << col <<
')' << std::endl;
2518 block(row, col).print(out, alternative_output);
2524 template <
typename MatrixType>
2528 return const_iterator(
this, 0);
2533 template <
typename MatrixType>
2537 return const_iterator(
this, m());
2542 template <
typename MatrixType>
2547 return const_iterator(
this, r);
2552 template <
typename MatrixType>
2557 return const_iterator(
this, r + 1);
2562 template <
typename MatrixType>
2566 return iterator(
this, 0);
2571 template <
typename MatrixType>
2575 return iterator(
this, m());
2580 template <
typename MatrixType>
2585 return iterator(
this, r);
2590 template <
typename MatrixType>
2595 return iterator(
this, r + 1);
2600 template <
typename MatrixType>
2604 std::vector<size_type> row_sizes(this->n_block_rows());
2605 std::vector<size_type> col_sizes(this->n_block_cols());
2609 for (
unsigned int r = 0; r < this->n_block_rows(); ++r)
2610 row_sizes[r] = sub_objects[r][0]->m();
2614 for (
unsigned int c = 1; c < this->n_block_cols(); ++c)
2615 for (
unsigned int r = 0; r < this->n_block_rows(); ++r)
2616 Assert(row_sizes[r] == sub_objects[r][c]->m(),
2617 ExcIncompatibleRowNumbers(r, 0, r, c));
2621 this->row_block_indices.reinit(row_sizes);
2625 for (
unsigned int c = 0; c < this->n_block_cols(); ++c)
2626 col_sizes[c] = sub_objects[0][c]->n();
2627 for (
unsigned int r = 1; r < this->n_block_rows(); ++r)
2628 for (
unsigned int c = 0; c < this->n_block_cols(); ++c)
2629 Assert(col_sizes[c] == sub_objects[r][c]->n(),
2630 ExcIncompatibleRowNumbers(0, c, r, c));
2634 this->column_block_indices.reinit(col_sizes);
2639 template <
typename MatrixType>
2643 for (
unsigned int row = 0; row < n_block_rows(); ++row)
2644 for (
unsigned int col = 0; col < n_block_cols(); ++col)
2645 block(row, col).prepare_add();
2650 template <
typename MatrixType>
2654 for (
unsigned int row = 0; row < n_block_rows(); ++row)
2655 for (
unsigned int col = 0; col < n_block_cols(); ++col)
2656 block(row, col).prepare_set();
bool operator==(const AlignedVector< T > &lhs, const AlignedVector< T > &rhs)
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
const_iterator end() const
BlockIndices column_block_indices
value_type matrix_norm_square(const BlockVectorType &v) const
void Tvmult_block_block(BlockVectorType &dst, const BlockVectorType &src) const
BlockType & block(const unsigned int row, const unsigned int column)
const value_type & const_reference
void Tvmult_nonblock_block(VectorType &dst, const BlockVectorType &src) const
BlockMatrixBase & operator*=(const value_type factor)
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
BlockMatrixBase & operator/=(const value_type factor)
types::global_dof_index size_type
Table< 2, SmartPointer< BlockType, BlockMatrixBase< MatrixType > > > sub_objects
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
void set(const size_type i, const size_type j, const value_type value)
std::size_t memory_consumption() const
void add(const value_type factor, const BlockMatrixBase< MatrixType > &matrix)
typename numbers::NumberTraits< value_type >::real_type real_type
BlockIndices row_block_indices
const BlockIndices & get_row_indices() const
~BlockMatrixBase() override
value_type matrix_scalar_product(const BlockVectorType &u, const BlockVectorType &v) const
BlockMatrixBase & copy_from(const BlockMatrixType &source)
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
iterator begin(const size_type r)
const_iterator begin(const size_type r) const
const BlockType & block(const unsigned int row, const unsigned int column) const
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)
const BlockIndices & get_column_indices() const
unsigned int block_row() const
types::global_dof_index size_type
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
const BlockMatrixType MatrixType
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
#define DEAL_II_NAMESPACE_CLOSE
__global__ void set(Number *val, const Number s, const size_type N)
static ::ExceptionBase & ExcDivideByZero()
static ::ExceptionBase & ExcInternalError()
#define DeclException4(Exception4, type1, type2, type3, type4, outsequence)
static ::ExceptionBase & ExcNotInitialized()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
static ::ExceptionBase & ExcIteratorPastEnd()
#define Assert(cond, exc)
static ::ExceptionBase & ExcNotImplemented()
#define AssertIsFinite(number)
static ::ExceptionBase & ExcIncompatibleColNumbers(int arg1, int arg2, int arg3, int arg4)
#define AssertIndexRange(index, range)
static ::ExceptionBase & ExcIncompatibleRowNumbers(int arg1, int arg2, int arg3, int arg4)
static ::ExceptionBase & ExcMessage(std::string arg1)
static ::ExceptionBase & ExcNotQuadratic()
@ matrix
Contents is actually a matrix.
types::global_dof_index size_type
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)
static const unsigned int invalid_unsigned_int
const types::global_dof_index invalid_size_type
unsigned int global_dof_index
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)
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
void advance(std::tuple< I1, I2 > &t, const unsigned int n)