16 #ifndef dealii__petsc_matrix_base_h 17 #define dealii__petsc_matrix_base_h 20 #include <deal.II/base/config.h> 22 #ifdef DEAL_II_WITH_PETSC 24 # include <deal.II/base/subscriptor.h> 25 # include <deal.II/lac/exceptions.h> 26 # include <deal.II/lac/full_matrix.h> 27 # include <deal.II/lac/petsc_compatibility.h> 28 # include <deal.II/lac/vector.h> 30 # include <petscmat.h> 31 # include <deal.II/base/std_cxx11/shared_ptr.h> 36 DEAL_II_NAMESPACE_OPEN
47 namespace MatrixIterators
83 const size_type
index);
88 size_type
row()
const;
93 size_type
index()
const;
103 PetscScalar
value()
const;
114 <<
"You tried to access row " << arg1
115 <<
" of a distributed matrix, but only rows " 116 << arg2 <<
" through " << arg3
117 <<
" are stored locally and can be accessed.");
152 std_cxx11::shared_ptr<const std::vector<PetscScalar> >
value_cache;
179 const size_type
index);
223 <<
"Attempt to access element " << arg2
224 <<
" of row " << arg1
225 <<
" which doesn't have that many elements.");
307 operator = (
const value_type d);
323 void set (
const size_type i,
325 const PetscScalar
value);
346 void set (
const std::vector<size_type> &indices,
348 const bool elide_zero_values =
false);
355 void set (
const std::vector<size_type> &row_indices,
356 const std::vector<size_type> &col_indices,
358 const bool elide_zero_values =
false);
374 void set (
const size_type
row,
375 const std::vector<size_type > &col_indices,
376 const std::vector<PetscScalar> &values,
377 const bool elide_zero_values =
false);
393 void set (
const size_type
row,
394 const size_type n_cols,
395 const size_type *col_indices,
396 const PetscScalar *values,
397 const bool elide_zero_values =
false);
408 void add (
const size_type i,
410 const PetscScalar
value);
431 void add (
const std::vector<size_type> &indices,
433 const bool elide_zero_values =
true);
440 void add (
const std::vector<size_type> &row_indices,
441 const std::vector<size_type> &col_indices,
443 const bool elide_zero_values =
true);
459 void add (
const size_type row,
460 const std::vector<size_type> &col_indices,
461 const std::vector<PetscScalar> &values,
462 const bool elide_zero_values =
true);
478 void add (
const size_type row,
479 const size_type n_cols,
480 const size_type *col_indices,
481 const PetscScalar *values,
482 const bool elide_zero_values =
true,
483 const bool col_indices_are_sorted =
false);
501 void clear_row (
const size_type row,
502 const PetscScalar new_diag_value = 0);
512 void clear_rows (
const std::vector<size_type> &rows,
513 const PetscScalar new_diag_value = 0);
539 PetscScalar operator () (
const size_type i,
540 const size_type j)
const;
549 PetscScalar el (
const size_type i,
550 const size_type j)
const;
561 PetscScalar diag_element (
const size_type i)
const;
566 size_type m ()
const;
571 size_type n ()
const;
581 size_type local_size ()
const;
591 std::pair<size_type, size_type>
592 local_range ()
const;
598 bool in_local_range (
const size_type
index)
const;
604 virtual const MPI_Comm &get_mpi_communicator ()
const = 0;
611 size_type n_nonzero_elements ()
const;
616 size_type row_length (
const size_type row)
const;
625 PetscReal l1_norm ()
const;
634 PetscReal linfty_norm ()
const;
640 PetscReal frobenius_norm ()
const;
662 PetscScalar matrix_norm_square (
const VectorBase &v)
const;
678 PetscScalar matrix_scalar_product (
const VectorBase &u,
682 #if DEAL_II_PETSC_VERSION_GTE(3,1,0) 687 PetscScalar trace ()
const;
693 MatrixBase &operator *= (
const PetscScalar factor);
698 MatrixBase &operator /= (
const PetscScalar factor);
715 const PetscScalar factor) DEAL_II_DEPRECATED;
795 const_iterator begin ()
const;
800 const_iterator end ()
const;
810 const_iterator begin (
const size_type r)
const;
820 const_iterator end (
const size_type r)
const;
829 operator Mat ()
const;
841 is_symmetric (
const double tolerance = 1.e-12);
849 is_hermitian (
const double tolerance = 1.e-12);
857 void write_ascii (
const PetscViewerFormat format = PETSC_VIEWER_DEFAULT);
866 void print (std::ostream &out,
867 const bool alternative_output =
false)
const;
872 std::size_t memory_consumption()
const;
884 <<
"You tried to do a " 889 <<
" operation but the matrix is currently in " 894 <<
" mode. You first have to call 'compress()'.");
920 void assert_is_compressed();
971 template <
class>
friend class ::BlockMatrixBase;
980 namespace MatrixIterators
987 const size_type
index)
989 matrix(const_cast<MatrixBase *>(matrix)),
1037 const size_type row,
1038 const size_type index)
1076 const const_iterator old_state = *
this;
1113 return ! (*
this == other);
1120 operator < (
const const_iterator &other)
const 1142 const PetscScalar
value)
1146 set (i, 1, &j, &
value,
false);
1155 const bool elide_zero_values)
1157 Assert (indices.size() == values.
m(),
1161 for (size_type i=0; i<indices.size(); ++i)
1162 set (indices[i], indices.size(), &indices[0], &values(i,0),
1171 const std::vector<size_type> &col_indices,
1173 const bool elide_zero_values)
1175 Assert (row_indices.size() == values.
m(),
1177 Assert (col_indices.size() == values.
n(),
1180 for (size_type i=0; i<row_indices.size(); ++i)
1181 set (row_indices[i], col_indices.size(), &col_indices[0], &values(i,0),
1190 const std::vector<size_type> &col_indices,
1191 const std::vector<PetscScalar> &values,
1192 const bool elide_zero_values)
1194 Assert (col_indices.size() == values.size(),
1197 set (
row, col_indices.size(), &col_indices[0], &values[0],
1206 const size_type n_cols,
1207 const size_type *col_indices,
1208 const PetscScalar *values,
1209 const bool elide_zero_values)
1211 (void)elide_zero_values;
1215 const PetscInt petsc_i =
row;
1216 PetscInt *col_index_ptr;
1218 PetscScalar
const *col_value_ptr;
1222 #ifndef PETSC_USE_64BIT_INDICES 1223 if (elide_zero_values ==
false)
1225 col_index_ptr = (
int *)col_indices;
1226 col_value_ptr = values;
1234 if (column_indices.size() < n_cols)
1236 column_indices.resize(n_cols);
1237 column_values.resize(n_cols);
1241 for (size_type j=0; j<n_cols; ++j)
1243 const PetscScalar value = values[j];
1245 if (value != PetscScalar())
1247 column_indices[n_columns] = col_indices[j];
1248 column_values[n_columns] =
value;
1254 col_index_ptr = &column_indices[0];
1255 col_value_ptr = &column_values[0];
1258 const PetscErrorCode ierr = MatSetValues (
matrix, 1, &petsc_i, n_columns,
1260 col_value_ptr, INSERT_VALUES);
1270 const PetscScalar value)
1275 if (value == PetscScalar())
1287 add (i, 1, &j, &value,
false);
1296 const bool elide_zero_values)
1298 Assert (indices.size() == values.
m(),
1302 for (size_type i=0; i<indices.size(); ++i)
1303 add (indices[i], indices.size(), &indices[0], &values(i,0),
1312 const std::vector<size_type> &col_indices,
1314 const bool elide_zero_values)
1316 Assert (row_indices.size() == values.
m(),
1318 Assert (col_indices.size() == values.
n(),
1321 for (size_type i=0; i<row_indices.size(); ++i)
1322 add (row_indices[i], col_indices.size(), &col_indices[0], &values(i,0),
1331 const std::vector<size_type> &col_indices,
1332 const std::vector<PetscScalar> &values,
1333 const bool elide_zero_values)
1335 Assert (col_indices.size() == values.size(),
1338 add (row, col_indices.size(), &col_indices[0], &values[0],
1347 const size_type n_cols,
1348 const size_type *col_indices,
1349 const PetscScalar *values,
1350 const bool elide_zero_values,
1353 (void)elide_zero_values;
1357 const PetscInt petsc_i =
row;
1358 PetscInt *col_index_ptr;
1360 PetscScalar
const *col_value_ptr;
1364 #ifndef PETSC_USE_64BIT_INDICES 1365 if (elide_zero_values ==
false)
1367 col_index_ptr = (
int *)col_indices;
1368 col_value_ptr = values;
1376 if (column_indices.size() < n_cols)
1378 column_indices.resize(n_cols);
1379 column_values.resize(n_cols);
1383 for (size_type j=0; j<n_cols; ++j)
1385 const PetscScalar value = values[j];
1387 if (value != PetscScalar())
1389 column_indices[n_columns] = col_indices[j];
1390 column_values[n_columns] =
value;
1396 col_index_ptr = &column_indices[0];
1397 col_value_ptr = &column_values[0];
1400 const PetscErrorCode ierr = MatSetValues (
matrix, 1, &petsc_i, n_columns,
1401 col_index_ptr, col_value_ptr,
1414 const size_type j)
const 1442 if (row_length(r) > 0)
1458 for (size_type i=r+1; i<m(); ++i)
1459 if (row_length(i) > 0)
1473 PetscInt begin, end;
1475 const PetscErrorCode ierr = MatGetOwnershipRange (static_cast<const Mat &>(
matrix),
1479 return ((index >= static_cast<size_type>(begin)) &&
1480 (index < static_cast<size_type>(end)));
1490 last_action = new_action;
1492 Assert (last_action == new_action, ExcWrongMode (last_action, new_action));
1504 ExcMessage(
"Error: missing compress() call."));
1529 DEAL_II_NAMESPACE_CLOSE
1532 #endif // DEAL_II_WITH_PETSC types::global_dof_index size_type
#define DeclException2(Exception2, type1, type2, outsequence)
void add(const size_type i, const size_type j, const PetscScalar value)
static::ExceptionBase & ExcInvalidIndexWithinRow(int arg1, int arg2)
std_cxx11::shared_ptr< const std::vector< PetscScalar > > value_cache
PetscScalar value() const
#define AssertThrow(cond, exc)
bool operator<(const const_iterator &) const
static::ExceptionBase & ExcIndexRange(int arg1, int arg2, int arg3)
const_iterator begin() const
std::vector< PetscInt > column_indices
const_iterator(const MatrixBase *matrix, const size_type row, const size_type index)
types::global_dof_index size_type
std::vector< PetscScalar > column_values
MatrixIterators::const_iterator const_iterator
static::ExceptionBase & ExcMessage(std::string arg1)
const_iterator end() const
unsigned int global_dof_index
Accessor(const MatrixBase *matrix, const size_type row, const size_type index)
bool operator==(const const_iterator &) const
#define Assert(cond, exc)
static::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
bool in_local_range(const size_type index) const
bool operator!=(const const_iterator &) const
PetscScalar operator()(const size_type i, const size_type j) const
#define DeclException0(Exception0)
size_type row_length(const size_type row) const
types::global_dof_index size_type
std_cxx11::shared_ptr< const std::vector< size_type > > colnum_cache
static::ExceptionBase & ExcIteratorPastEnd()
friend class const_iterator
static::ExceptionBase & ExcNotQuadratic()
VectorOperation::values last_action
void prepare_action(const VectorOperation::values new_action)
const_iterator & operator++()
static::ExceptionBase & ExcAccessToNonlocalRow(int arg1, int arg2, int arg3)
void set(const size_type i, const size_type j, const PetscScalar value)
const Accessor & operator*() const
#define DeclException3(Exception3, type1, type2, type3, outsequence)
static::ExceptionBase & ExcBeyondEndOfMatrix()
#define AssertIsFinite(number)
void assert_is_compressed()
static::ExceptionBase & ExcInternalError()
const Accessor * operator->() const