17 #ifndef __deal2__petsc_matrix_base_h 18 #define __deal2__petsc_matrix_base_h 21 #include <deal.II/base/config.h> 23 #ifdef DEAL_II_WITH_PETSC 25 # include <deal.II/base/subscriptor.h> 26 # include <deal.II/lac/full_matrix.h> 27 # include <deal.II/lac/exceptions.h> 28 # include <deal.II/lac/vector.h> 30 # include <petscmat.h> 31 # include <deal.II/base/std_cxx1x/shared_ptr.h> 47 namespace MatrixIterators
85 const size_type
index);
97 size_type
row()
const;
104 size_type
index()
const;
116 PetscScalar
value()
const;
127 <<
"You tried to access row " << arg1
128 <<
" of a distributed matrix, but only rows " 129 << arg2 <<
" through " << arg3
130 <<
" are stored locally and can be accessed.");
178 std_cxx1x::shared_ptr<const std::vector<PetscScalar> >
value_cache;
209 const size_type
index);
259 <<
"Attempt to access element " << arg2
260 <<
" of row " << arg1
261 <<
" which doesn't have that many elements.");
353 operator = (
const value_type d);
377 void set (
const size_type i,
379 const PetscScalar
value);
416 void set (
const std::vector<size_type> &indices,
418 const bool elide_zero_values =
false);
427 void set (
const std::vector<size_type> &row_indices,
428 const std::vector<size_type> &col_indices,
430 const bool elide_zero_values =
false);
458 void set (
const size_type
row,
459 const std::vector<size_type > &col_indices,
460 const std::vector<PetscScalar> &values,
461 const bool elide_zero_values =
false);
489 void set (
const size_type
row,
490 const size_type n_cols,
491 const size_type *col_indices,
492 const PetscScalar *values,
493 const bool elide_zero_values =
false);
510 void add (
const size_type i,
512 const PetscScalar
value);
551 void add (
const std::vector<size_type> &indices,
553 const bool elide_zero_values =
true);
562 void add (
const std::vector<size_type> &row_indices,
563 const std::vector<size_type> &col_indices,
565 const bool elide_zero_values =
true);
594 void add (
const size_type row,
595 const std::vector<size_type> &col_indices,
596 const std::vector<PetscScalar> &values,
597 const bool elide_zero_values =
true);
626 void add (
const size_type row,
627 const size_type n_cols,
628 const size_type *col_indices,
629 const PetscScalar *values,
630 const bool elide_zero_values =
true,
631 const bool col_indices_are_sorted =
false);
663 void clear_row (
const size_type row,
664 const PetscScalar new_diag_value = 0);
679 void clear_rows (
const std::vector<size_type> &rows,
680 const PetscScalar new_diag_value = 0);
699 void compress (::VectorOperation::values operation);
722 PetscScalar operator () (const size_type i,
723 const size_type j) const;
736 PetscScalar el (const size_type i,
737 const size_type j) const;
754 PetscScalar diag_element (const size_type i) const;
760 size_type m () const;
766 size_type n () const;
781 size_type local_size () const;
799 std::pair<size_type, size_type>
800 local_range () const;
807 bool in_local_range (const size_type
index) const;
815 virtual const MPI_Comm &get_mpi_communicator () const = 0;
826 size_type n_nonzero_elements () const;
831 size_type row_length (const size_type row) const;
845 PetscReal l1_norm () const;
860 PetscReal linfty_norm () const;
868 PetscReal frobenius_norm () const;
907 PetscScalar matrix_norm_square (const
VectorBase &v) const;
932 PetscScalar matrix_scalar_product (const
VectorBase &u,
936 #if DEAL_II_PETSC_VERSION_GTE(3,1,0) 943 PetscScalar trace ()
const;
950 MatrixBase &operator *= (
const PetscScalar factor);
956 MatrixBase &operator /= (
const PetscScalar factor);
963 const PetscScalar factor);
1092 const_iterator begin ()
const;
1097 const_iterator end ()
const;
1111 const_iterator begin (
const size_type r)
const;
1124 const_iterator end (
const size_type r)
const;
1137 operator Mat ()
const;
1152 #if DEAL_II_PETSC_VERSION_LT(3,2,0) 1157 is_symmetric (
const double tolerance = 1.e-12);
1168 #if DEAL_II_PETSC_VERSION_LT(3,2,0) 1173 is_hermitian (
const double tolerance = 1.e-12);
1183 void write_ascii (
const PetscViewerFormat format = PETSC_VIEWER_DEFAULT);
1189 void print (std::ostream &out,
1190 const bool alternative_output =
false)
const;
1196 std::size_t memory_consumption()
const;
1203 <<
"An error with error number " << arg1
1204 <<
" occurred while calling a PETSc function");
1215 <<
"You tried to do a " 1220 <<
" operation but the matrix is currently in " 1225 <<
" mode. You first have to call 'compress()'.");
1258 enum Values { none, insert, add };
1275 void prepare_action(
const LastAction::Values new_action);
1346 template <
class>
friend class ::BlockMatrixBase;
1355 namespace MatrixIterators
1361 const size_type row,
1362 const size_type
index)
1364 matrix(const_cast<MatrixBase *>(matrix)),
1388 Assert (a_row < matrix->m(), ExcBeyondEndOfMatrix());
1397 Assert (a_row < matrix->m(), ExcBeyondEndOfMatrix());
1406 Assert (a_row < matrix->m(), ExcBeyondEndOfMatrix());
1415 Assert (a_row < matrix->m(), ExcBeyondEndOfMatrix());
1423 const size_type row,
1424 const size_type index)
1462 const const_iterator old_state = *
this;
1499 return ! (*
this == other);
1506 operator < (
const const_iterator &other)
const 1528 const PetscScalar
value)
1532 set (i, 1, &j, &
value,
false);
1541 const bool elide_zero_values)
1543 Assert (indices.size() == values.
m(),
1545 Assert (values.
m() == values.
n(), ExcNotQuadratic());
1547 for (size_type i=0; i<indices.size(); ++i)
1548 set (indices[i], indices.size(), &indices[0], &values(i,0),
1557 const std::vector<size_type> &col_indices,
1559 const bool elide_zero_values)
1561 Assert (row_indices.size() == values.
m(),
1563 Assert (col_indices.size() == values.
n(),
1566 for (size_type i=0; i<row_indices.size(); ++i)
1567 set (row_indices[i], col_indices.size(), &col_indices[0], &values(i,0),
1576 const std::vector<size_type> &col_indices,
1577 const std::vector<PetscScalar> &values,
1578 const bool elide_zero_values)
1580 Assert (col_indices.size() == values.size(),
1583 set (
row, col_indices.size(), &col_indices[0], &values[0],
1592 const size_type n_cols,
1593 const size_type *col_indices,
1594 const PetscScalar *values,
1595 const bool elide_zero_values)
1597 prepare_action(LastAction::insert);
1599 const PetscInt petsc_i =
row;
1600 PetscInt *col_index_ptr;
1602 PetscScalar
const *col_value_ptr;
1607 #ifndef PETSC_USE_64BIT_INDICES 1608 if (elide_zero_values ==
false)
1610 col_index_ptr = (
int *)col_indices;
1611 col_value_ptr = values;
1619 if (column_indices.size() < n_cols)
1621 column_indices.resize(n_cols);
1622 column_values.resize(n_cols);
1626 for (size_type j=0; j<n_cols; ++j)
1628 const PetscScalar value = values[j];
1630 if (value != PetscScalar())
1632 column_indices[n_columns] = col_indices[j];
1633 column_values[n_columns] =
value;
1639 col_index_ptr = &column_indices[0];
1640 col_value_ptr = &column_values[0];
1644 = MatSetValues (
matrix, 1, &petsc_i, n_columns, col_index_ptr,
1645 col_value_ptr, INSERT_VALUES);
1655 const PetscScalar value)
1660 if (value == PetscScalar())
1672 prepare_action(LastAction::add);
1677 add (i, 1, &j, &value,
false);
1686 const bool elide_zero_values)
1688 Assert (indices.size() == values.
m(),
1690 Assert (values.
m() == values.
n(), ExcNotQuadratic());
1692 for (size_type i=0; i<indices.size(); ++i)
1693 add (indices[i], indices.size(), &indices[0], &values(i,0),
1702 const std::vector<size_type> &col_indices,
1704 const bool elide_zero_values)
1706 Assert (row_indices.size() == values.
m(),
1708 Assert (col_indices.size() == values.
n(),
1711 for (size_type i=0; i<row_indices.size(); ++i)
1712 add (row_indices[i], col_indices.size(), &col_indices[0], &values(i,0),
1721 const std::vector<size_type> &col_indices,
1722 const std::vector<PetscScalar> &values,
1723 const bool elide_zero_values)
1725 Assert (col_indices.size() == values.size(),
1728 add (row, col_indices.size(), &col_indices[0], &values[0],
1737 const size_type n_cols,
1738 const size_type *col_indices,
1739 const PetscScalar *values,
1740 const bool elide_zero_values,
1743 prepare_action(LastAction::add);
1745 const PetscInt petsc_i =
row;
1746 PetscInt *col_index_ptr;
1748 PetscScalar
const *col_value_ptr;
1753 #ifndef PETSC_USE_64BIT_INDICES 1754 if (elide_zero_values ==
false)
1756 col_index_ptr = (
int *)col_indices;
1757 col_value_ptr = values;
1765 if (column_indices.size() < n_cols)
1767 column_indices.resize(n_cols);
1768 column_values.resize(n_cols);
1772 for (size_type j=0; j<n_cols; ++j)
1774 const PetscScalar value = values[j];
1776 if (value != PetscScalar())
1778 column_indices[n_columns] = col_indices[j];
1779 column_values[n_columns] =
value;
1785 col_index_ptr = &column_indices[0];
1786 col_value_ptr = &column_values[0];
1790 = MatSetValues (
matrix, 1, &petsc_i, n_columns, col_index_ptr,
1791 col_value_ptr, ADD_VALUES);
1803 const size_type j)
const 1831 if (row_length(r) > 0)
1847 for (size_type i=r+1; i<m(); ++i)
1848 if (row_length(i) > 0)
1862 PetscInt begin, end;
1864 const int ierr = MatGetOwnershipRange (static_cast<const Mat &>(
matrix),
1868 return ((index >= static_cast<size_type>(begin)) &&
1869 (index < static_cast<size_type>(end)));
1878 if (last_action == new_action)
1880 else if (last_action == LastAction::none)
1881 last_action = new_action;
1883 Assert (
false, ExcWrongMode (last_action, new_action));
1892 prepare_action(LastAction::add);
1901 prepare_action(LastAction::insert);
1908 DEAL_II_NAMESPACE_CLOSE
1911 #endif // DEAL_II_WITH_PETSC types::global_dof_index size_type
std_cxx1x::shared_ptr< const std::vector< size_type > > colnum_cache
void add(const size_type i, const size_type j, const PetscScalar value)
void prepare_action(const LastAction::Values new_action)
PetscScalar value() const
#define AssertThrow(cond, exc)
bool operator<(const const_iterator &) const
const_iterator begin() const
bool is_finite(const double x)
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
const_iterator end() const
#define DeclException1(Exception1, type1, outsequence)
DeclException0(ExcBeyondEndOfMatrix)
unsigned int global_dof_index
size_type row_length(const size_type row) const
Accessor(const MatrixBase *matrix, const size_type row, const size_type index)
bool operator==(const const_iterator &) const
#define Assert(cond, exc)
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
::ExceptionBase & ExcIndexRange(int arg1, int arg2, int arg3)
types::global_dof_index size_type
BlockCompressedSparsityPattern CompressedBlockSparsityPattern DEAL_II_DEPRECATED
friend class const_iterator
::ExceptionBase & ExcIteratorPastEnd()
const_iterator & operator++()
::ExceptionBase & ExcNumberNotFinite()
LastAction::Values last_action
void set(const size_type i, const size_type j, const PetscScalar value)
DeclException3(ExcAccessToNonlocalRow, int, int, int,<< "You tried to access row "<< arg1<< " of a distributed matrix, but only rows "<< arg2<< " through "<< arg3<< " are stored locally and can be accessed.")
const Accessor & operator*() const
DeclException2(ExcInvalidIndexWithinRow, int, int,<< "Attempt to access element "<< arg2<< " of row "<< arg1<< " which doesn't have that many elements.")
::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
::ExceptionBase & ExcInternalError()
std_cxx1x::shared_ptr< const std::vector< PetscScalar > > value_cache
const Accessor * operator->() const