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>
914 operator Mat()
const;
952 write_ascii(
const PetscViewerFormat format = PETSC_VIEWER_DEFAULT);
962 print(std::ostream &out,
const bool alternative_output =
false)
const;
974 "You are attempting an operation on two vectors that "
975 "are the same object, but the operation requires that the "
976 "two objects are in fact different.");
984 <<
"You tried to do a "
985 << (arg1 == 1 ?
"'set'" : (arg1 == 2 ?
"'add'" :
"???"))
986 <<
" operation but the matrix is currently in "
987 << (arg2 == 1 ?
"'set'" : (arg2 == 2 ?
"'add'" :
"???"))
988 <<
" mode. You first have to call 'compress()'.");
1105 friend class ::BlockMatrixBase;
1114 namespace MatrixIterators
1123 visit_present_row();
1140 return (*colnum_cache)[a_index];
1156 return (*value_cache)[a_index];
1168 inline const_iterator &
1177 if (accessor.a_index >= accessor.colnum_cache->size())
1179 accessor.a_index = 0;
1182 while ((accessor.a_row < accessor.matrix->m()) &&
1183 (accessor.a_row < accessor.matrix->local_range().second) &&
1184 (accessor.matrix->row_length(accessor.a_row) == 0))
1187 accessor.visit_present_row();
1193 inline const_iterator
1196 const const_iterator old_state = *
this;
1202 inline const const_iterator::Accessor &
1209 inline const const_iterator::Accessor *
1210 const_iterator::operator->()
const
1219 return (accessor.a_row == other.accessor.a_row &&
1220 accessor.a_index == other.accessor.a_index);
1227 return !(*
this == other);
1234 return (accessor.row() < other.accessor.row() ||
1235 (accessor.row() == other.accessor.row() &&
1236 accessor.index() < other.accessor.index()));
1255 set(i, 1, &j, &value,
false);
1263 const bool elide_zero_values)
1269 for (
size_type i = 0; i < indices.size(); ++i)
1281 const std::vector<size_type> & col_indices,
1283 const bool elide_zero_values)
1290 for (
size_type i = 0; i < row_indices.size(); ++i)
1302 const std::vector<size_type> & col_indices,
1303 const std::vector<PetscScalar> &
values,
1304 const bool elide_zero_values)
1322 const PetscScalar *
values,
1323 const bool elide_zero_values)
1327 const PetscInt petsc_i = row;
1328 PetscInt
const *col_index_ptr;
1330 PetscScalar
const *col_value_ptr;
1334 if (elide_zero_values ==
false)
1336 col_index_ptr =
reinterpret_cast<const PetscInt *
>(col_indices);
1344 if (column_indices.size() < n_cols)
1346 column_indices.resize(n_cols);
1347 column_values.resize(n_cols);
1355 if (value != PetscScalar())
1357 column_indices[n_columns] = col_indices[j];
1358 column_values[n_columns] =
value;
1364 col_index_ptr = column_indices.data();
1365 col_value_ptr = column_values.data();
1368 const PetscErrorCode ierr = MatSetValues(
matrix,
1385 if (value == PetscScalar())
1396 add(i, 1, &j, &value,
false);
1402 MatrixBase::add(
const std::vector<size_type> & indices,
1404 const bool elide_zero_values)
1410 for (
size_type i = 0; i < indices.size(); ++i)
1421 MatrixBase::add(
const std::vector<size_type> & row_indices,
1422 const std::vector<size_type> & col_indices,
1424 const bool elide_zero_values)
1431 for (
size_type i = 0; i < row_indices.size(); ++i)
1443 const std::vector<size_type> & col_indices,
1444 const std::vector<PetscScalar> &
values,
1445 const bool elide_zero_values)
1463 const PetscScalar *
values,
1464 const bool elide_zero_values,
1467 (void)elide_zero_values;
1471 const PetscInt petsc_i = row;
1472 PetscInt
const *col_index_ptr;
1474 PetscScalar
const *col_value_ptr;
1478 if (elide_zero_values ==
false)
1480 col_index_ptr =
reinterpret_cast<const PetscInt *
>(col_indices);
1488 if (column_indices.size() < n_cols)
1490 column_indices.resize(n_cols);
1491 column_values.resize(n_cols);
1499 if (value != PetscScalar())
1501 column_indices[n_columns] = col_indices[j];
1502 column_values[n_columns] =
value;
1508 col_index_ptr = column_indices.data();
1509 col_value_ptr = column_values.data();
1512 const PetscErrorCode ierr = MatSetValues(
1513 matrix, 1, &petsc_i, n_columns, col_index_ptr, col_value_ptr, ADD_VALUES);
1527 inline MatrixBase::const_iterator
1531 (in_local_range(0) && in_local_range(m() - 1)),
1533 "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."));
1538 while ((first_nonempty_row < m()) && (row_length(first_nonempty_row) == 0))
1539 ++first_nonempty_row;
1541 return const_iterator(
this, first_nonempty_row, 0);
1545 inline MatrixBase::const_iterator
1549 (in_local_range(0) && in_local_range(m() - 1)),
1551 "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."));
1553 return const_iterator(
this, m(), 0);
1557 inline MatrixBase::const_iterator
1560 Assert(in_local_range(r),
1563 if (row_length(r) > 0)
1564 return const_iterator(
this, r, 0);
1570 inline MatrixBase::const_iterator
1573 Assert(in_local_range(r),
1587 if (i == local_range().second || (row_length(i) > 0))
1588 return const_iterator(
this, i, 0);
1594 return {
this, m(), 0};
1600 MatrixBase::in_local_range(
const size_type index)
const
1604 const PetscErrorCode ierr =
1605 MatGetOwnershipRange(
static_cast<const Mat &
>(
matrix), &
begin, &
end);
1618 last_action = new_action;
1620 Assert(last_action == new_action, ExcWrongMode(last_action, new_action));
1626 MatrixBase::assert_is_compressed()
1631 ExcMessage(
"Error: missing compress() call."));
1637 MatrixBase::prepare_add()
1645 MatrixBase::prepare_set()
1650 inline const MPI_Comm &
1651 MatrixBase::get_mpi_communicator()
const
1653 static MPI_Comm
comm = PETSC_COMM_SELF;
1654 MPI_Comm pcomm = PetscObjectComm(
reinterpret_cast<PetscObject
>(
matrix));
1655 if (pcomm != MPI_COMM_NULL)
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
PetscScalar diag_element(const size_type i) const
MatrixBase(const MatrixBase &)=delete
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
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)
const MPI_Comm & get_mpi_communicator() const
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< T >::value &&std::is_floating_point< U >::value, 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)