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>
36 DEAL_II_NAMESPACE_OPEN
41 namespace PETScWrappers
47 namespace MatrixIterators
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;
259 <<
"Attempt to access element " << arg2
260 <<
" of row " << arg1
261 <<
" which doesn't have that many elements.");
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);
459 const std::vector<size_type > &col_indices,
460 const std::vector<PetscScalar> &values,
461 const bool elide_zero_values =
false);
492 const PetscScalar *values,
493 const bool elide_zero_values =
false);
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);
595 const std::vector<size_type> &col_indices,
596 const std::vector<PetscScalar> &values,
597 const bool elide_zero_values =
true);
629 const PetscScalar *values,
630 const bool elide_zero_values =
true,
631 const bool col_indices_are_sorted =
false);
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,
831 size_type
row_length (const size_type row) const;
936 #if DEAL_II_PETSC_VERSION_GTE(3,1,0)
943 PetscScalar trace ()
const;
963 const PetscScalar factor);
1137 operator Mat ()
const;
1152 #if DEAL_II_PETSC_VERSION_LT(3,2,0)
1168 #if DEAL_II_PETSC_VERSION_LT(3,2,0)
1183 void write_ascii (
const PetscViewerFormat format = PETSC_VIEWER_DEFAULT);
1189 void print (std::ostream &out,
1190 const bool alternative_output =
false)
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 };
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)),
1368 visit_present_row ();
1379 colnum_cache (a.colnum_cache),
1380 value_cache (a.value_cache)
1385 const_iterator::Accessor::size_type
1386 const_iterator::Accessor::row()
const
1388 Assert (a_row < matrix->m(), ExcBeyondEndOfMatrix());
1394 const_iterator::Accessor::size_type
1395 const_iterator::Accessor::column()
const
1397 Assert (a_row < matrix->m(), ExcBeyondEndOfMatrix());
1398 return (*colnum_cache)[a_index];
1403 const_iterator::Accessor::size_type
1404 const_iterator::Accessor::index()
const
1406 Assert (a_row < matrix->m(), ExcBeyondEndOfMatrix());
1413 const_iterator::Accessor::value()
const
1415 Assert (a_row < matrix->m(), ExcBeyondEndOfMatrix());
1416 return (*value_cache)[a_index];
1422 const_iterator(
const MatrixBase *matrix,
1423 const size_type row,
1424 const size_type index)
1426 accessor(matrix, row, index)
1433 const_iterator::operator++ ()
1442 if (accessor.a_index >= accessor.colnum_cache->size())
1444 accessor.a_index = 0;
1447 while ((accessor.a_row < accessor.matrix->m())
1449 (accessor.matrix->row_length(accessor.a_row) == 0))
1452 accessor.visit_present_row();
1460 const_iterator::operator++ (
int)
1462 const const_iterator old_state = *
this;
1469 const const_iterator::Accessor &
1470 const_iterator::operator* ()
const
1477 const const_iterator::Accessor *
1478 const_iterator::operator-> ()
const
1487 operator == (
const const_iterator &other)
const
1489 return (accessor.a_row == other.accessor.a_row &&
1490 accessor.a_index == other.accessor.a_index);
1497 operator != (
const const_iterator &other)
const
1499 return ! (*
this == other);
1506 operator < (
const const_iterator &other)
const
1508 return (accessor.row() < other.accessor.row() ||
1509 (accessor.row() == other.accessor.row() &&
1510 accessor.index() < other.accessor.index()));
1526 MatrixBase::set (
const size_type i,
1528 const PetscScalar value)
1532 set (i, 1, &j, &value,
false);
1539 MatrixBase::set (
const std::vector<size_type> &indices,
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),
1556 MatrixBase::set (
const std::vector<size_type> &row_indices,
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),
1575 MatrixBase::set (
const size_type row,
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],
1591 MatrixBase::set (
const size_type row,
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);
1653 MatrixBase::add (
const size_type i,
1655 const PetscScalar value)
1660 if (value == PetscScalar())
1672 prepare_action(LastAction::add);
1677 add (i, 1, &j, &value,
false);
1684 MatrixBase::add (
const std::vector<size_type> &indices,
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),
1701 MatrixBase::add (
const std::vector<size_type> &row_indices,
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),
1720 MatrixBase::add (
const size_type row,
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],
1736 MatrixBase::add (
const size_type row,
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);
1802 MatrixBase::operator() (
const size_type i,
1803 const size_type j)
const
1811 MatrixBase::const_iterator
1812 MatrixBase::begin()
const
1814 return const_iterator(
this, 0, 0);
1819 MatrixBase::const_iterator
1820 MatrixBase::end()
const
1822 return const_iterator(
this, m(), 0);
1827 MatrixBase::const_iterator
1828 MatrixBase::begin(
const size_type r)
const
1831 if (row_length(r) > 0)
1832 return const_iterator(
this, r, 0);
1839 MatrixBase::const_iterator
1840 MatrixBase::end(
const size_type r)
const
1847 for (size_type i=r+1; i<m(); ++i)
1848 if (row_length(i) > 0)
1849 return const_iterator(
this, i, 0);
1860 MatrixBase::in_local_range (
const size_type index)
const
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)));
1876 MatrixBase::prepare_action(
const LastAction::Values new_action)
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));
1890 MatrixBase::prepare_add()
1892 prepare_action(LastAction::add);
1899 MatrixBase::prepare_set()
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
PetscScalar matrix_scalar_product(const VectorBase &u, const VectorBase &v) const
void add(const size_type i, const size_type j, const PetscScalar value)
PetscScalar residual(VectorBase &dst, const VectorBase &x, const VectorBase &b) const
void prepare_action(const LastAction::Values new_action)
void compress() DEAL_II_DEPRECATED
void Tvmult(VectorBase &dst, const VectorBase &src) const
PetscScalar value() const
#define AssertThrow(cond, exc)
bool operator<(const const_iterator &) const
PetscReal linfty_norm() const
size_type local_size() const
DeclException2(ExcWrongMode, int, int,<< "You tried to do a "<< (arg1==1?"'set'":(arg1==2?"'add'":"???"))<< " operation but the matrix is currently in "<< (arg2==1?"'set'":(arg2==2?"'add'":"???"))<< " mode. You first have to call 'compress()'.")
PetscBool is_symmetric(const double tolerance=1.e-12)
PetscBool is_hermitian(const double tolerance=1.e-12)
const_iterator begin() const
bool is_finite(const double x)
void clear_rows(const std::vector< size_type > &rows, const PetscScalar new_diag_value=0)
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
std::size_t memory_consumption() const
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
::ExceptionBase & ExcIndexRange(int arg1, int arg2, int arg3)
types::global_dof_index size_type
BlockCompressedSparsityPattern CompressedBlockSparsityPattern DEAL_II_DEPRECATED
PetscReal l1_norm() const
PetscScalar diag_element(const size_type i) const
PetscScalar el(const size_type i, const size_type j) const
::ExceptionBase & ExcIteratorPastEnd()
void write_ascii(const PetscViewerFormat format=PETSC_VIEWER_DEFAULT)
void print(std::ostream &out, const bool alternative_output=false) const
PetscScalar matrix_norm_square(const VectorBase &v) const
DeclException0(ExcSourceEqualsDestination)
const_iterator & operator++()
::ExceptionBase & ExcNumberNotFinite()
MatrixBase & operator/=(const PetscScalar factor)
std::pair< size_type, size_type > local_range() const
void Tvmult_add(VectorBase &dst, const VectorBase &src) const
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.")
PetscReal frobenius_norm() const
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)
MatrixBase & operator=(const value_type d)
DeclException1(ExcPETScError, int,<< "An error with error number "<< arg1<< " occurred while calling a PETSc function")
void vmult(VectorBase &dst, const VectorBase &src) const
::ExceptionBase & ExcInternalError()
virtual const MPI_Comm & get_mpi_communicator() const =0
std_cxx1x::shared_ptr< const std::vector< PetscScalar > > value_cache
size_type n_nonzero_elements() const
void clear_row(const size_type row, const PetscScalar new_diag_value=0)
MatrixBase & operator*=(const PetscScalar factor)
void vmult_add(VectorBase &dst, const VectorBase &src) const
const Accessor * operator->() const