18 #ifdef DEAL_II_WITH_PETSC
30 namespace MatrixIterators
39 if (
matrix->in_local_range(this->a_row) ==
false)
49 const PetscInt *colnums;
53 MatGetRow(*
matrix, this->a_row, &ncols, &colnums, &
values);
64 std::make_shared<std::vector<size_type>>(colnums, colnums + ncols);
66 std::make_shared<std::vector<PetscScalar>>(
values,
values + ncols);
69 ierr = MatRestoreRow(*
matrix, this->a_row, &ncols, &colnums, &
values);
87 const PetscErrorCode ierr =
88 PetscObjectReference(
reinterpret_cast<PetscObject
>(
matrix));
98 PetscObjectReference(
reinterpret_cast<PetscObject
>(
A));
100 ierr = MatDestroy(&
matrix);
107 PetscErrorCode ierr = MatDestroy(&
matrix);
117 const PetscErrorCode ierr = MatDestroy(&
matrix);
123 const int m = 0,
n = 0, n_nonzero_per_row = 0;
124 const PetscErrorCode ierr = MatCreateSeqAIJ(
125 PETSC_COMM_SELF,
m,
n, n_nonzero_per_row,
nullptr, &
matrix);
139 const PetscErrorCode ierr = MatZeroEntries(
matrix);
150 std::vector<size_type> rows(1, row);
158 const PetscScalar new_diag_value)
164 const std::vector<PetscInt> petsc_rows(rows.begin(), rows.end());
179 ierr = MatZeroRowsIS(
matrix, index_set, new_diag_value,
nullptr,
nullptr);
181 ierr = ISDestroy(&index_set);
187 const PetscScalar new_diag_value)
193 const std::vector<PetscInt> petsc_rows(rows.begin(), rows.end());
209 MatZeroRowsColumnsIS(
matrix, index_set, new_diag_value,
nullptr,
nullptr);
211 ierr = ISDestroy(&index_set);
220 PetscInt petsc_i = i, petsc_j = j;
224 const PetscErrorCode ierr =
225 MatGetValues(
matrix, 1, &petsc_i, 1, &petsc_j, &value);
253 int all_int_last_action;
255 const int ierr = MPI_Allreduce(&my_int_last_action,
256 &all_int_last_action,
265 ExcMessage(
"Error: not all processors agree on the last "
266 "VectorOperation before this compress() call."));
273 "Missing compress() or calling with wrong VectorOperation argument."));
276 PetscErrorCode ierr = MatAssemblyBegin(
matrix, MAT_FINAL_ASSEMBLY);
279 ierr = MatAssemblyEnd(
matrix, MAT_FINAL_ASSEMBLY);
290 PetscInt n_rows, n_cols;
292 const PetscErrorCode ierr = MatGetSize(
matrix, &n_rows, &n_cols);
303 PetscInt n_rows, n_cols;
305 const PetscErrorCode ierr = MatGetSize(
matrix, &n_rows, &n_cols);
318 const PetscErrorCode ierr = MatGetLocalSize(
matrix, &n_rows,
nullptr);
326 std::pair<MatrixBase::size_type, MatrixBase::size_type>
331 const PetscErrorCode ierr =
332 MatGetOwnershipRange(
static_cast<const Mat &
>(
matrix), &
begin, &
end);
345 const PetscErrorCode ierr = MatGetLocalSize(
matrix,
nullptr, &n_cols);
353 std::pair<MatrixBase::size_type, MatrixBase::size_type>
358 const PetscErrorCode ierr =
359 MatGetOwnershipRangeColumn(
static_cast<const Mat &
>(
matrix),
373 const PetscErrorCode ierr = MatGetInfo(
matrix, MAT_GLOBAL_SUM, &mat_info);
378 return static_cast<std::uint64_t
>(mat_info.nz_used);
398 const PetscInt *colnums;
399 const PetscScalar *
values;
403 PetscErrorCode ierr = MatGetRow(*
this, row, &ncols, &colnums, &
values);
411 const PetscInt ncols_saved = ncols;
412 ierr = MatRestoreRow(*
this, row, &ncols, &colnums, &
values);
424 const PetscErrorCode ierr = MatNorm(
matrix, NORM_1, &result);
437 const PetscErrorCode ierr = MatNorm(
matrix, NORM_INFINITY, &result);
450 const PetscErrorCode ierr = MatNorm(
matrix, NORM_FROBENIUS, &result);
481 const PetscErrorCode ierr = MatGetTrace(
matrix, &result);
492 const PetscErrorCode ierr = MatScale(
matrix, a);
503 const PetscScalar factor = 1. / a;
504 const PetscErrorCode ierr = MatScale(
matrix, factor);
515 const PetscErrorCode ierr =
516 MatAXPY(
matrix, factor, other, DIFFERENT_NONZERO_PATTERN);
528 const PetscErrorCode ierr = MatMult(
matrix, src, dst);
539 const PetscErrorCode ierr = MatMultTranspose(
matrix, src, dst);
550 const PetscErrorCode ierr = MatMultAdd(
matrix, src, dst, dst);
561 const PetscErrorCode ierr = MatMultTransposeAdd(
matrix, src, dst, dst);
573 const bool transpose_left)
575 const bool use_vector = (
V.size() == inputright.
m() ? true :
false);
576 if (transpose_left ==
false)
578 Assert(inputleft.
n() == inputright.
m(),
583 Assert(inputleft.
m() == inputright.
m(),
591 if (use_vector ==
false)
595 ierr = MatTransposeMatMult(inputleft,
604 ierr = MatMatMult(inputleft,
615 ierr = MatDuplicate(inputleft, MAT_COPY_VALUES, &tmp);
619 # if DEAL_II_PETSC_VERSION_LT(3, 8, 0)
620 ierr = MatTranspose(tmp, MAT_REUSE_MATRIX, &tmp);
622 ierr = MatTranspose(tmp, MAT_INPLACE_MATRIX, &tmp);
626 ierr = MatDiagonalScale(tmp,
nullptr,
V);
628 ierr = MatMatMult(tmp,
634 ierr = MatDestroy(&tmp);
673 MatrixBase::operator Mat()
const
687 # if DEAL_II_PETSC_VERSION_LT(3, 8, 0)
688 const PetscErrorCode ierr = MatTranspose(
matrix, MAT_REUSE_MATRIX, &
matrix);
690 const PetscErrorCode ierr =
701 const PetscErrorCode ierr = MatIsSymmetric(
matrix, tolerance, &truth);
712 const PetscErrorCode ierr = MatIsHermitian(
matrix, tolerance, &truth);
722 MPI_Comm
comm = PetscObjectComm(
reinterpret_cast<PetscObject
>(
matrix));
725 PetscErrorCode ierr =
726 PetscViewerSetFormat(PETSC_VIEWER_STDOUT_(
comm), format);
730 ierr = MatView(
matrix, PETSC_VIEWER_STDOUT_(
comm));
739 PetscErrorCode ierr = MatHasOperation(
matrix, MATOP_GET_ROW, &has);
745 ierr = MatConvert(
matrix, MATAIJ, MAT_INITIAL_MATRIX, &vmatrix);
749 std::pair<MatrixBase::size_type, MatrixBase::size_type> loc_range =
753 const PetscInt *colnums;
754 const PetscScalar *
values;
757 for (row = loc_range.first; row < loc_range.second; ++row)
759 ierr = MatGetRow(vmatrix, row, &ncols, &colnums, &
values);
762 for (PetscInt col = 0; col < ncols; ++col)
764 out <<
"(" << row <<
"," << colnums[col] <<
") " <<
values[col]
768 ierr = MatRestoreRow(vmatrix, row, &ncols, &colnums, &
values);
773 ierr = MatDestroy(&vmatrix);
785 const PetscErrorCode ierr = MatGetInfo(
matrix, MAT_LOCAL, &info);
788 return sizeof(*this) +
static_cast<size_type>(info.memory);
void add(const size_type i, const size_type j, const PetscScalar value)
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)
MPI_Comm get_mpi_communicator() const
PetscScalar diag_element(const size_type i) const
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
PetscScalar el(const size_type i, const size_type j) const
size_type local_size() const
PetscScalar trace() 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 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
std::uint64_t n_nonzero_elements() const
void clear_row(const size_type row, const PetscScalar new_diag_value=0)
PetscReal linfty_norm() const
real_type l2_norm() const
#define DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_CLOSE
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcScalarAssignmentOnlyForZeroValue()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
#define Assert(cond, exc)
#define AssertThrowMPI(error_code)
#define AssertNothrow(cond, exc)
static ::ExceptionBase & ExcIO()
static ::ExceptionBase & ExcSourceEqualsDestination()
static ::ExceptionBase & ExcMessage(std::string arg1)
static ::ExceptionBase & ExcNotQuadratic()
#define AssertThrow(cond, exc)
@ matrix
Contents is actually a matrix.
void perform_mmult(const MatrixBase &inputleft, const MatrixBase &inputright, MatrixBase &result, const VectorBase &V, const bool transpose_left)
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)