15#ifndef dealii_petsc_matrix_base_h
16#define dealii_petsc_matrix_base_h
21#ifdef DEAL_II_WITH_PETSC
31# include <boost/container/small_vector.hpp>
44template <
typename Matrix>
54 namespace MatrixIterators
134 <<
"You tried to access row " << arg1
135 <<
" of a distributed matrix, but only rows " << arg2
136 <<
" through " << arg3
137 <<
" are stored locally and can be accessed.");
250 <<
"Attempt to access element " << arg2 <<
" of row "
251 << arg1 <<
" which doesn't have that many elements.");
417 set(
const std::vector<size_type> &indices,
419 const bool elide_zero_values =
false);
427 set(
const std::vector<size_type> &row_indices,
428 const std::vector<size_type> &col_indices,
430 const bool elide_zero_values =
false);
448 const std::vector<size_type> &col_indices,
449 const std::vector<PetscScalar> &values,
450 const bool elide_zero_values =
false);
470 const PetscScalar *values,
471 const bool elide_zero_values =
false);
505 add(
const std::vector<size_type> &indices,
507 const bool elide_zero_values =
true);
515 add(
const std::vector<size_type> &row_indices,
516 const std::vector<size_type> &col_indices,
518 const bool elide_zero_values =
true);
536 const std::vector<size_type> &col_indices,
537 const std::vector<PetscScalar> &values,
538 const bool elide_zero_values =
true);
558 const PetscScalar *values,
559 const bool elide_zero_values =
true,
560 const bool col_indices_are_sorted =
false);
591 const PetscScalar new_diag_value = 0);
598 const PetscScalar new_diag_value = 0);
681 std::pair<size_type, size_type>
706 std::pair<size_type, size_type>
941 operator Mat()
const;
979 write_ascii(
const PetscViewerFormat format = PETSC_VIEWER_DEFAULT);
989 print(std::ostream &out,
const bool alternative_output =
false)
const;
1001 "You are attempting an operation on two vectors that "
1002 "are the same object, but the operation requires that the "
1003 "two objects are in fact different.");
1011 <<
"You tried to do a "
1012 << (arg1 == 1 ?
"'set'" : (arg1 == 2 ?
"'add'" :
"???"))
1013 <<
" operation but the matrix is currently in "
1014 << (arg2 == 1 ?
"'set'" : (arg2 == 2 ?
"'add'" :
"???"))
1015 <<
" mode. You first have to call 'compress()'.");
1105 friend class ::BlockMatrixBase;
1116 const PetscScalar *values,
1117 const bool elide_zero_values);
1126 namespace MatrixIterators
1129 const size_type row,
1130 const size_type index)
1135 visit_present_row();
1140 inline const_iterator::Accessor::size_type
1141 const_iterator::Accessor::row()
const
1143 Assert(a_row < matrix->m(), ExcBeyondEndOfMatrix());
1148 inline const_iterator::Accessor::size_type
1149 const_iterator::Accessor::column()
const
1151 Assert(a_row < matrix->m(), ExcBeyondEndOfMatrix());
1152 return (*colnum_cache)[a_index];
1156 inline const_iterator::Accessor::size_type
1157 const_iterator::Accessor::index()
const
1159 Assert(a_row < matrix->m(), ExcBeyondEndOfMatrix());
1165 const_iterator::Accessor::value()
const
1167 Assert(a_row < matrix->m(), ExcBeyondEndOfMatrix());
1168 return (*value_cache)[a_index];
1172 inline const_iterator::const_iterator(
const MatrixBase *matrix,
1173 const size_type row,
1174 const size_type index)
1180 inline const_iterator &
1181 const_iterator::operator++()
1189 if (accessor.a_index >= accessor.colnum_cache->size())
1191 accessor.a_index = 0;
1194 while ((accessor.a_row < accessor.matrix->m()) &&
1195 (accessor.a_row < accessor.matrix->local_range().second) &&
1196 (accessor.matrix->row_length(accessor.a_row) == 0))
1199 accessor.visit_present_row();
1205 inline const_iterator
1206 const_iterator::operator++(
int)
1208 const const_iterator old_state = *
this;
1214 inline const const_iterator::Accessor &
1215 const_iterator::operator*()
const
1221 inline const const_iterator::Accessor *
1222 const_iterator::operator->()
const
1229 const_iterator::operator==(
const const_iterator &other)
const
1231 return (accessor.a_row == other.accessor.a_row &&
1232 accessor.a_index == other.accessor.a_index);
1237 const_iterator::operator!=(
const const_iterator &other)
const
1239 return !(*
this == other);
1244 const_iterator::operator<(
const const_iterator &other)
const
1246 return (accessor.row() < other.accessor.row() ||
1247 (accessor.row() == other.accessor.row() &&
1248 accessor.index() < other.accessor.index()));
1263 MatrixBase::set(
const size_type i,
const size_type j,
const PetscScalar value)
1267 set(i, 1, &j, &value,
false);
1273 MatrixBase::set(
const std::vector<size_type> &indices,
1275 const bool elide_zero_values)
1281 for (size_type i = 0; i < indices.size(); ++i)
1292 MatrixBase::set(
const std::vector<size_type> &row_indices,
1293 const std::vector<size_type> &col_indices,
1295 const bool elide_zero_values)
1302 for (size_type i = 0; i < row_indices.size(); ++i)
1313 MatrixBase::set(
const size_type row,
1314 const std::vector<size_type> &col_indices,
1315 const std::vector<PetscScalar> &values,
1316 const bool elide_zero_values)
1331 MatrixBase::set(
const size_type row,
1332 const size_type n_cols,
1333 const size_type *col_indices,
1334 const PetscScalar *values,
1335 const bool elide_zero_values)
1348 MatrixBase::add(
const size_type i,
const size_type j,
const PetscScalar value)
1352 if (value == PetscScalar())
1363 add(i, 1, &j, &value,
false);
1369 MatrixBase::add(
const std::vector<size_type> &indices,
1371 const bool elide_zero_values)
1377 for (size_type i = 0; i < indices.size(); ++i)
1388 MatrixBase::add(
const std::vector<size_type> &row_indices,
1389 const std::vector<size_type> &col_indices,
1391 const bool elide_zero_values)
1398 for (size_type i = 0; i < row_indices.size(); ++i)
1409 MatrixBase::add(
const size_type row,
1410 const std::vector<size_type> &col_indices,
1411 const std::vector<PetscScalar> &values,
1412 const bool elide_zero_values)
1427 MatrixBase::add(
const size_type row,
1428 const size_type n_cols,
1429 const size_type *col_indices,
1430 const PetscScalar *values,
1431 const bool elide_zero_values,
1446 const size_type row,
1447 const size_type n_cols,
1448 const size_type *col_indices,
1449 const PetscScalar *values,
1450 const bool elide_zero_values)
1452 prepare_action(operation);
1454 const auto petsc_row =
static_cast<PetscInt
>(row);
1464 boost::container::small_vector<PetscInt, 100> column_indices;
1465 std::optional<boost::container::small_vector<PetscScalar, 100>>
1468 const PetscScalar *values_ptr =
nullptr;
1469 if (elide_zero_values ==
false)
1471 column_indices.resize(n_cols);
1473 for (size_type j = 0; j < n_cols; ++j)
1476 column_indices[j] =
static_cast<PetscInt
>(col_indices[j]);
1486 column_values.emplace();
1487 for (size_type j = 0; j < n_cols; ++j)
1491 if (value != PetscScalar())
1493 column_indices.push_back(
static_cast<PetscInt
>(col_indices[j]));
1495 column_values->push_back(value);
1498 values_ptr = column_values->data();
1501 const auto petsc_n_columns =
static_cast<PetscInt
>(column_indices.size());
1507 const PetscErrorCode ierr =
1508 MatSetValues(matrix,
1512 column_indices.data(),
1522 MatrixBase::operator()(
const size_type i,
const size_type j)
const
1529 inline MatrixBase::const_iterator
1530 MatrixBase::begin()
const
1533 (in_local_range(0) && in_local_range(m() - 1)),
1535 "begin() and end() can only be called on a processor owning the entire matrix. If this is a distributed matrix, use begin(row) and end(row) instead."));
1540 while ((first_nonempty_row < m()) && (row_length(first_nonempty_row) == 0))
1541 ++first_nonempty_row;
1543 return const_iterator(
this, first_nonempty_row, 0);
1547 inline MatrixBase::const_iterator
1548 MatrixBase::end()
const
1551 (in_local_range(0) && in_local_range(m() - 1)),
1553 "begin() and end() can only be called on a processor owning the entire matrix. If this is a distributed matrix, use begin(row) and end(row) instead."));
1555 return const_iterator(
this, m(), 0);
1559 inline MatrixBase::const_iterator
1560 MatrixBase::begin(
const size_type r)
const
1562 Assert(in_local_range(r),
1565 if (row_length(r) > 0)
1566 return const_iterator(
this, r, 0);
1572 inline MatrixBase::const_iterator
1573 MatrixBase::end(
const size_type r)
const
1575 Assert(in_local_range(r),
1588 for (size_type i = r + 1; i < m(); ++i)
1589 if (i == local_range().second || (row_length(i) > 0))
1590 return const_iterator(
this, i, 0);
1596 return {
this, m(), 0};
1602 MatrixBase::in_local_range(
const size_type index)
const
1604 PetscInt petsc_begin, petsc_end;
1606 const PetscErrorCode ierr =
1607 MatGetOwnershipRange(
static_cast<const Mat &
>(matrix),
1617 return ((index >= begin) && (index < end));
1626 last_action = new_action;
1628 Assert(last_action == new_action, ExcWrongMode(last_action, new_action));
1634 MatrixBase::assert_is_compressed()
1639 ExcMessage(
"Error: missing compress() call."));
1645 MatrixBase::prepare_add()
1653 MatrixBase::prepare_set()
1659 MatrixBase::get_mpi_communicator()
const
1661 return PetscObjectComm(
reinterpret_cast<PetscObject
>(matrix));
void add(const size_type i, const size_type j, const PetscScalar value)
void add(const std::vector< size_type > &row_indices, const std::vector< size_type > &col_indices, const FullMatrix< PetscScalar > &full_matrix, const bool elide_zero_values=true)
const_iterator end(const size_type r) const
size_type row_length(const size_type row) const
std::size_t memory_consumption() const
PetscReal l1_norm() const
VectorOperation::values last_action
void vmult(VectorBase &dst, const VectorBase &src) const
MPI_Comm get_mpi_communicator() const
PetscScalar diag_element(const size_type i) const
MatrixBase(const MatrixBase &)=delete
size_type local_domain_size() const
const_iterator begin() const
MatrixBase & operator/=(const PetscScalar factor)
void mmult(MatrixBase &C, const MatrixBase &B, const VectorBase &V) const
PetscBool is_symmetric(const double tolerance=1.e-12)
PetscReal frobenius_norm() const
virtual ~MatrixBase() override
std::pair< size_type, size_type > local_domain() const
const_iterator end() const
void Tvmult_add(VectorBase &dst, const VectorBase &src) const
void print(std::ostream &out, const bool alternative_output=false) const
const_iterator begin(const size_type r) const
void add_or_set(const VectorOperation::values &operation, const size_type row, const size_type n_cols, const size_type *col_indices, const PetscScalar *values, const bool elide_zero_values)
void add(const size_type row, const std::vector< size_type > &col_indices, const std::vector< PetscScalar > &values, const bool elide_zero_values=true)
PetscScalar el(const size_type i, const size_type j) const
size_type local_size() const
PetscScalar operator()(const size_type i, const size_type j) const
MatrixBase & operator=(const MatrixBase &)=delete
void set(const std::vector< size_type > &row_indices, const std::vector< size_type > &col_indices, const FullMatrix< PetscScalar > &full_matrix, const bool elide_zero_values=false)
PetscScalar trace() const
void prepare_action(const VectorOperation::values new_action)
void add(const std::vector< size_type > &indices, const FullMatrix< PetscScalar > &full_matrix, const bool elide_zero_values=true)
bool in_local_range(const size_type index) const
PetscBool is_hermitian(const double tolerance=1.e-12)
PetscScalar matrix_scalar_product(const VectorBase &u, const VectorBase &v) const
void assert_is_compressed()
void Tvmult(VectorBase &dst, const VectorBase &src) const
void clear_rows_columns(const std::vector< size_type > &row_and_column_indices, const PetscScalar new_diag_value=0)
MatrixBase & operator*=(const PetscScalar factor)
PetscScalar residual(VectorBase &dst, const VectorBase &x, const VectorBase &b) const
void Tmmult(MatrixBase &C, const MatrixBase &B, const VectorBase &V) const
void set(const size_type row, const std::vector< size_type > &col_indices, const std::vector< PetscScalar > &values, const bool elide_zero_values=false)
void add(const size_type row, const size_type n_cols, const size_type *col_indices, const PetscScalar *values, const bool elide_zero_values=true, const bool col_indices_are_sorted=false)
void write_ascii(const PetscViewerFormat format=PETSC_VIEWER_DEFAULT)
void vmult_add(VectorBase &dst, const VectorBase &src) const
std::pair< size_type, size_type > local_range() const
void compress(const VectorOperation::values operation)
PetscScalar matrix_norm_square(const VectorBase &v) const
void set(const size_type row, const size_type n_cols, const size_type *col_indices, const PetscScalar *values, const bool elide_zero_values=false)
void clear_rows(const ArrayView< const size_type > &rows, const PetscScalar new_diag_value=0)
void set(const std::vector< size_type > &indices, const FullMatrix< PetscScalar > &full_matrix, const bool elide_zero_values=false)
void set(const size_type i, const size_type j, const PetscScalar value)
std::uint64_t n_nonzero_elements() const
void clear_row(const size_type row, const PetscScalar new_diag_value=0)
PetscReal linfty_norm() const
std::shared_ptr< const std::vector< PetscScalar > > value_cache
PetscScalar value() const
std::shared_ptr< const std::vector< size_type > > colnum_cache
Accessor(const MatrixBase *matrix, const size_type row, const size_type index)
const_iterator operator++(int)
bool operator!=(const const_iterator &) const
const_iterator(const MatrixBase *matrix, const size_type row, const size_type index)
bool operator==(const const_iterator &) const
bool operator<(const const_iterator &) const
const_iterator & operator++()
const Accessor * operator->() const
const Accessor & operator*() const
#define DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_CLOSE
#define DeclException0(Exception0)
static ::ExceptionBase & ExcAccessToNonlocalRow(int arg1, int arg2, int arg3)
#define Assert(cond, exc)
static ::ExceptionBase & ExcIteratorPastEnd()
static ::ExceptionBase & ExcInvalidIndexWithinRow(int arg1, int arg2)
#define AssertIsFinite(number)
#define DeclException2(Exception2, type1, type2, outsequence)
static ::ExceptionBase & ExcBeyondEndOfMatrix()
#define DeclExceptionMsg(Exception, defaulttext)
static ::ExceptionBase & ExcSourceEqualsDestination()
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcWrongMode(int arg1, int arg2)
#define DeclException3(Exception3, type1, type2, type3, outsequence)
static ::ExceptionBase & ExcIndexRange(std::size_t arg1, std::size_t arg2, std::size_t arg3)
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
static ::ExceptionBase & ExcNotQuadratic()
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
#define AssertIntegerConversion(index1, index2)
types::global_dof_index size_type
@ matrix
Contents is actually a matrix.
VectorType::value_type * end(VectorType &V)
VectorType::value_type * begin(VectorType &V)
unsigned int global_dof_index