16 #ifndef dealii_petsc_matrix_base_h
17 #define dealii_petsc_matrix_base_h
22 #ifdef DEAL_II_WITH_PETSC
32 # include <petscmat.h>
42 template <
typename Matrix>
52 namespace MatrixIterators
124 <<
"You tried to access row " << arg1
125 <<
" of a distributed matrix, but only rows " << arg2
126 <<
" through " << arg3
127 <<
" are stored locally and can be accessed.");
240 <<
"Attempt to access element " << arg2 <<
" of row "
241 << arg1 <<
" which doesn't have that many elements.");
407 set(
const std::vector<size_type> &indices,
409 const bool elide_zero_values =
false);
417 set(
const std::vector<size_type> &row_indices,
418 const std::vector<size_type> &col_indices,
420 const bool elide_zero_values =
false);
438 const std::vector<size_type> &col_indices,
439 const std::vector<PetscScalar> &
values,
440 const bool elide_zero_values =
false);
460 const PetscScalar *
values,
461 const bool elide_zero_values =
false);
495 add(
const std::vector<size_type> &indices,
497 const bool elide_zero_values =
true);
505 add(
const std::vector<size_type> &row_indices,
506 const std::vector<size_type> &col_indices,
508 const bool elide_zero_values =
true);
526 const std::vector<size_type> &col_indices,
527 const std::vector<PetscScalar> &
values,
528 const bool elide_zero_values =
true);
548 const PetscScalar *
values,
549 const bool elide_zero_values =
true,
550 const bool col_indices_are_sorted =
false);
580 clear_rows(
const std::vector<size_type> &rows,
581 const PetscScalar new_diag_value = 0);
588 const PetscScalar new_diag_value = 0);
671 std::pair<size_type, size_type>
696 std::pair<size_type, size_type>
931 operator Mat()
const;
969 write_ascii(
const PetscViewerFormat format = PETSC_VIEWER_DEFAULT);
979 print(std::ostream &out,
const bool alternative_output =
false)
const;
991 "You are attempting an operation on two vectors that "
992 "are the same object, but the operation requires that the "
993 "two objects are in fact different.");
1001 <<
"You tried to do a "
1002 << (arg1 == 1 ?
"'set'" : (arg1 == 2 ?
"'add'" :
"???"))
1003 <<
" operation but the matrix is currently in "
1004 << (arg2 == 1 ?
"'set'" : (arg2 == 2 ?
"'add'" :
"???"))
1005 <<
" mode. You first have to call 'compress()'.");
1122 friend class ::BlockMatrixBase;
1131 namespace MatrixIterators
1140 visit_present_row();
1157 return (*colnum_cache)[a_index];
1173 return (*value_cache)[a_index];
1185 inline const_iterator &
1194 if (accessor.a_index >= accessor.colnum_cache->size())
1196 accessor.a_index = 0;
1199 while ((accessor.a_row < accessor.matrix->m()) &&
1200 (accessor.a_row < accessor.matrix->local_range().second) &&
1201 (accessor.matrix->row_length(accessor.a_row) == 0))
1204 accessor.visit_present_row();
1210 inline const_iterator
1213 const const_iterator old_state = *
this;
1219 inline const const_iterator::Accessor &
1226 inline const const_iterator::Accessor *
1227 const_iterator::operator->()
const
1236 return (accessor.a_row == other.accessor.a_row &&
1237 accessor.a_index == other.accessor.a_index);
1244 return !(*
this == other);
1251 return (accessor.row() < other.accessor.row() ||
1252 (accessor.row() == other.accessor.row() &&
1253 accessor.index() < other.accessor.index()));
1272 set(i, 1, &j, &value,
false);
1280 const bool elide_zero_values)
1286 for (
size_type i = 0; i < indices.size(); ++i)
1298 const std::vector<size_type> &col_indices,
1300 const bool elide_zero_values)
1307 for (
size_type i = 0; i < row_indices.size(); ++i)
1319 const std::vector<size_type> &col_indices,
1320 const std::vector<PetscScalar> &
values,
1321 const bool elide_zero_values)
1339 const PetscScalar *
values,
1340 const bool elide_zero_values)
1344 const PetscInt petsc_i = row;
1345 const PetscInt *col_index_ptr;
1347 const PetscScalar *col_value_ptr;
1351 if (elide_zero_values ==
false)
1353 col_index_ptr =
reinterpret_cast<const PetscInt *
>(col_indices);
1361 if (column_indices.size() < n_cols)
1363 column_indices.resize(n_cols);
1364 column_values.resize(n_cols);
1372 if (value != PetscScalar())
1374 column_indices[n_columns] = col_indices[j];
1375 column_values[n_columns] =
value;
1381 col_index_ptr = column_indices.data();
1382 col_value_ptr = column_values.data();
1385 const PetscErrorCode ierr = MatSetValues(
matrix,
1402 if (value == PetscScalar())
1413 add(i, 1, &j, &value,
false);
1419 MatrixBase::add(
const std::vector<size_type> &indices,
1421 const bool elide_zero_values)
1427 for (
size_type i = 0; i < indices.size(); ++i)
1438 MatrixBase::add(
const std::vector<size_type> &row_indices,
1439 const std::vector<size_type> &col_indices,
1441 const bool elide_zero_values)
1448 for (
size_type i = 0; i < row_indices.size(); ++i)
1460 const std::vector<size_type> &col_indices,
1461 const std::vector<PetscScalar> &
values,
1462 const bool elide_zero_values)
1480 const PetscScalar *
values,
1481 const bool elide_zero_values,
1484 (void)elide_zero_values;
1488 const PetscInt petsc_i = row;
1489 const PetscInt *col_index_ptr;
1491 const PetscScalar *col_value_ptr;
1495 if (elide_zero_values ==
false)
1497 col_index_ptr =
reinterpret_cast<const PetscInt *
>(col_indices);
1505 if (column_indices.size() < n_cols)
1507 column_indices.resize(n_cols);
1508 column_values.resize(n_cols);
1516 if (value != PetscScalar())
1518 column_indices[n_columns] = col_indices[j];
1519 column_values[n_columns] =
value;
1525 col_index_ptr = column_indices.data();
1526 col_value_ptr = column_values.data();
1529 const PetscErrorCode ierr = MatSetValues(
1530 matrix, 1, &petsc_i, n_columns, col_index_ptr, col_value_ptr, ADD_VALUES);
1544 inline MatrixBase::const_iterator
1548 (in_local_range(0) && in_local_range(m() - 1)),
1550 "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 while ((first_nonempty_row < m()) && (row_length(first_nonempty_row) == 0))
1556 ++first_nonempty_row;
1558 return const_iterator(
this, first_nonempty_row, 0);
1562 inline MatrixBase::const_iterator
1566 (in_local_range(0) && in_local_range(m() - 1)),
1568 "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."));
1570 return const_iterator(
this, m(), 0);
1574 inline MatrixBase::const_iterator
1577 Assert(in_local_range(r),
1580 if (row_length(r) > 0)
1581 return const_iterator(
this, r, 0);
1587 inline MatrixBase::const_iterator
1590 Assert(in_local_range(r),
1604 if (i == local_range().second || (row_length(i) > 0))
1605 return const_iterator(
this, i, 0);
1611 return {
this, m(), 0};
1617 MatrixBase::in_local_range(
const size_type index)
const
1621 const PetscErrorCode ierr =
1622 MatGetOwnershipRange(
static_cast<const Mat &
>(
matrix), &
begin, &
end);
1635 last_action = new_action;
1637 Assert(last_action == new_action, ExcWrongMode(last_action, new_action));
1643 MatrixBase::assert_is_compressed()
1648 ExcMessage(
"Error: missing compress() call."));
1654 MatrixBase::prepare_add()
1662 MatrixBase::prepare_set()
1668 MatrixBase::get_mpi_communicator()
const
1670 return PetscObjectComm(
reinterpret_cast<PetscObject
>(
matrix));
bool operator!=(const AlignedVector< T > &lhs, const AlignedVector< T > &rhs)
bool operator==(const AlignedVector< T > &lhs, const AlignedVector< T > &rhs)
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
void clear_rows(const std::vector< size_type > &rows, const PetscScalar new_diag_value=0)
std::vector< PetscScalar > column_values
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)
types::global_dof_index size_type
PetscReal frobenius_norm() const
MatrixBase & operator=(const MatrixBase &)=delete
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(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
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()
std::vector< PetscInt > column_indices
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 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
types::global_dof_index size_type
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)
types::global_dof_index size_type
bool operator==(const const_iterator &) const
const_iterator & operator++()
bool operator<(const const_iterator &) const
const Accessor * operator->() const
const Accessor & operator*() const
std::enable_if_t< std::is_floating_point_v< T > &&std::is_floating_point_v< U >, typename ProductType< std::complex< T >, std::complex< U > >::type > operator*(const std::complex< T > &left, const std::complex< U > &right)
#define DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_CLOSE
__global__ void set(Number *val, const Number s, const size_type N)
#define DeclException0(Exception0)
static ::ExceptionBase & ExcInvalidIndexWithinRow(int arg1, int arg2)
static ::ExceptionBase & ExcBeyondEndOfMatrix()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
static ::ExceptionBase & ExcIteratorPastEnd()
static ::ExceptionBase & ExcBeyondEndOfMatrix()
#define Assert(cond, exc)
static ::ExceptionBase & ExcAccessToNonlocalRow(int arg1, int arg2, int arg3)
#define AssertIsFinite(number)
#define DeclException2(Exception2, type1, type2, outsequence)
#define AssertIndexRange(index, range)
#define DeclExceptionMsg(Exception, defaulttext)
static ::ExceptionBase & ExcWrongMode(int arg1, int arg2)
static ::ExceptionBase & ExcIndexRange(std::size_t arg1, std::size_t arg2, std::size_t arg3)
#define DeclException3(Exception3, type1, type2, type3, outsequence)
static ::ExceptionBase & ExcSourceEqualsDestination()
static ::ExceptionBase & ExcMessage(std::string arg1)
static ::ExceptionBase & ExcNotQuadratic()
#define AssertThrow(cond, exc)
@ matrix
Contents is actually a matrix.
types::global_dof_index size_type
SymmetricTensor< 2, dim, Number > C(const Tensor< 2, dim, Number > &F)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
VectorType::value_type * begin(VectorType &V)
VectorType::value_type * end(VectorType &V)
unsigned int global_dof_index
SynchronousIterators< Iterators > operator++(SynchronousIterators< Iterators > &a)
bool operator<(const SynchronousIterators< Iterators > &a, const SynchronousIterators< Iterators > &b)