Skip to content

Expand PetscMatrix interface #2169

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 3 commits into from
Jul 2, 2019
Merged

Expand PetscMatrix interface #2169

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
b45e739
Expand PetscMatrix interface
lindsayad Jul 2, 2019
93a2181
PETSC_NULL -> NULL
lindsayad Jul 2, 2019
fd629fe
Guard PETSc version
lindsayad Jul 2, 2019
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
22 changes: 22 additions & 0 deletions include/numerics/petsc_matrix.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -151,6 +151,11 @@ class PetscMatrix final : public SparseMatrix<T>
*/
void update_preallocation_and_zero();

/**
* Reset matrix to use the original nonzero pattern provided by users
*/
void reset_preallocation();

virtual void clear () override;

virtual void zero () override;
Expand All @@ -163,8 +168,25 @@ class PetscMatrix final : public SparseMatrix<T>

virtual numeric_index_type m () const override;

/**
* Get the number of rows owned by this process
*/
numeric_index_type local_m () const;
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Can we bump this up to the parent class and make the default implementation row_stop()-row_start()?


virtual numeric_index_type n () const override;

/**
* Get the number of columns owned by this process
*/
numeric_index_type local_n () const;
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'm confused by this one. My first interpretation here is "the total number of columns with non-empty entries in the sparsity pattern in rows owned by this process", but I'm not actually sure how PETSc would determine that - we give them n_nz and n_oz rather than complete sparsity info, right? And without the latter there's no way to immediately determine which non-empty columns from which rows overlap? So is this method only well-defined for matrices that have already been assembled??

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It is not related to nnz, noz or matrix assembly. It is a layout thing. It is something we tell PETSc to use in libMesh.

In PetscMatrix::init(), we have something like

ierr = MatSetSizes(_mat, m_local, n_local, m_global, n_global);
  LIBMESH_CHKERR(ierr);

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

"n_local" is "local_n()"

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks @fdkong. This is what I had determined through experimentation this morning 😄

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

That certainly clears it up! Thanks!


/**
* Get the number of rows and columns owned by this process
* @param i Row size
* @param j Column size
*/
void get_local_size (numeric_index_type & m, numeric_index_type & n) const;

virtual numeric_index_type row_start () const override;

virtual numeric_index_type row_stop () const override;
Expand Down
73 changes: 62 additions & 11 deletions src/numerics/petsc_matrix.C
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -172,11 +172,11 @@ void PetscMatrix<T>::init (const numeric_index_type m_in,

ierr = MatSetType(_mat, MATBAIJ); // Automatically chooses seqbaij or mpibaij
LIBMESH_CHKERR(ierr);
ierr = MatSeqBAIJSetPreallocation(_mat, blocksize, n_nz/blocksize, PETSC_NULL);
ierr = MatSeqBAIJSetPreallocation(_mat, blocksize, n_nz/blocksize, NULL);
LIBMESH_CHKERR(ierr);
ierr = MatMPIBAIJSetPreallocation(_mat, blocksize,
n_nz/blocksize, PETSC_NULL,
n_oz/blocksize, PETSC_NULL);
n_nz/blocksize, NULL,
n_oz/blocksize, NULL);
LIBMESH_CHKERR(ierr);
}
else
Expand All @@ -186,17 +186,17 @@ void PetscMatrix<T>::init (const numeric_index_type m_in,
case AIJ:
ierr = MatSetType(_mat, MATAIJ); // Automatically chooses seqaij or mpiaij
LIBMESH_CHKERR(ierr);
ierr = MatSeqAIJSetPreallocation(_mat, n_nz, PETSC_NULL);
ierr = MatSeqAIJSetPreallocation(_mat, n_nz, NULL);
LIBMESH_CHKERR(ierr);
ierr = MatMPIAIJSetPreallocation(_mat, n_nz, PETSC_NULL, n_oz, PETSC_NULL);
ierr = MatMPIAIJSetPreallocation(_mat, n_nz, NULL, n_oz, NULL);
LIBMESH_CHKERR(ierr);
break;

case HYPRE:
#if !PETSC_VERSION_LESS_THAN(3,9,4) && LIBMESH_HAVE_PETSC_HYPRE
ierr = MatSetType(_mat, MATHYPRE);
LIBMESH_CHKERR(ierr);
ierr = MatHYPRESetPreallocation(_mat, n_nz, PETSC_NULL, n_oz, PETSC_NULL);
ierr = MatHYPRESetPreallocation(_mat, n_nz, NULL, n_oz, NULL);
LIBMESH_CHKERR(ierr);
#else
libmesh_error_msg("PETSc 3.9.4 or higher with hypre is required for MatHypre");
Expand Down Expand Up @@ -470,6 +470,20 @@ void PetscMatrix<T>::update_preallocation_and_zero ()
libmesh_not_implemented();
}

template <typename T>
void PetscMatrix<T>::reset_preallocation()
{
#if !PETSC_VERSION_LESS_THAN(3,8,0)
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Oh I see.

#if !PETSC_VERSION_LESS_THAN(3,9,0)

8-->9

I was thinking about 9. Why I gave you 8

libmesh_assert (this->initialized());

auto ierr = MatResetPreallocation(_mat);
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Version control

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Ugh @fdkong you told me the wrong version number here. Hence the test failures on idaholab/moose#13682. Looks like it should be 3.9.0 based on git tag --contains. I'm guessing you just made a typo sigh.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

https://www.mcs.anl.gov/petsc/petsc-3.9/docs/manualpages/Mat/MatResetPreallocation.html

It is in 3.9. I am confused here.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

!PETSC_VERSION_LESS_THAN(3,8,0) means >= 3.8.0, i.e. including 3.8.0

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Right.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

We need >=3.9.0

LIBMESH_CHKERR(ierr);
#else
libmesh_warning("Your version of PETSc doesn't support resetting of "
"preallocation, so we will use your most recent sparsity "
"pattern. This may result in a degradation of performance\n");
#endif
}

template <typename T>
void PetscMatrix<T>::zero ()
Expand Down Expand Up @@ -506,7 +520,7 @@ void PetscMatrix<T>::zero_rows (std::vector<numeric_index_type> & rows, T diag_v
ierr = MatZeroRows(_mat, rows.size(),
numeric_petsc_cast(rows.data()), diag_value);
else
ierr = MatZeroRows(_mat, 0, PETSC_NULL, diag_value);
ierr = MatZeroRows(_mat, 0, NULL, diag_value);
#else
// As of petsc-dev at the time of 3.1.0, MatZeroRows now takes two additional
// optional arguments. The optional arguments (x,b) can be used to specify the
Expand All @@ -515,10 +529,10 @@ void PetscMatrix<T>::zero_rows (std::vector<numeric_index_type> & rows, T diag_v
if (!rows.empty())
ierr = MatZeroRows(_mat, cast_int<PetscInt>(rows.size()),
numeric_petsc_cast(rows.data()), diag_value,
PETSC_NULL, PETSC_NULL);
NULL, NULL);
else
ierr = MatZeroRows(_mat, 0, PETSC_NULL, diag_value, PETSC_NULL,
PETSC_NULL);
ierr = MatZeroRows(_mat, 0, NULL, diag_value, NULL,
NULL);
#endif

LIBMESH_CHKERR(ierr);
Expand Down Expand Up @@ -931,7 +945,7 @@ void PetscMatrix<T>::get_transpose (SparseMatrix<T> & dest) const
PetscErrorCode ierr;
#if PETSC_VERSION_LESS_THAN(3,0,0)
if (&petsc_dest == this)
ierr = MatTranspose(_mat,PETSC_NULL);
ierr = MatTranspose(_mat,NULL);
else
ierr = MatTranspose(_mat,&petsc_dest._mat);
LIBMESH_CHKERR(ierr);
Expand Down Expand Up @@ -1002,7 +1016,18 @@ numeric_index_type PetscMatrix<T>::m () const
return static_cast<numeric_index_type>(petsc_m);
}

template <typename T>
numeric_index_type PetscMatrix<T>::local_m () const
{
libmesh_assert (this->initialized());

PetscInt m = 0;

auto ierr = MatGetLocalSize (_mat, &m, NULL);
LIBMESH_CHKERR(ierr);

return static_cast<numeric_index_type>(m);
}

template <typename T>
numeric_index_type PetscMatrix<T>::n () const
Expand All @@ -1018,7 +1043,33 @@ numeric_index_type PetscMatrix<T>::n () const
return static_cast<numeric_index_type>(petsc_n);
}

template <typename T>
numeric_index_type PetscMatrix<T>::local_n () const
{
libmesh_assert (this->initialized());

PetscInt n = 0;

auto ierr = MatGetLocalSize (_mat, NULL, &n);
LIBMESH_CHKERR(ierr);

return static_cast<numeric_index_type>(n);
}

template <typename T>
void PetscMatrix<T>::get_local_size (numeric_index_type & m,
numeric_index_type & n) const
{
libmesh_assert (this->initialized());

PetscInt petsc_m = 0, petsc_n = 0;

auto ierr = MatGetLocalSize (_mat, &petsc_m, &petsc_n);
LIBMESH_CHKERR(ierr);

m = static_cast<numeric_index_type>(petsc_m);
n = static_cast<numeric_index_type>(petsc_n);
}

template <typename T>
numeric_index_type PetscMatrix<T>::row_start () const
Expand Down