[petsc-users] Setting/creating Mats for SNES use

[Justin Chang]

Hi all,

In my DMPlex program, I have these lines:

Mat A,J;

...

ierr = DMSetMatType(dm, MATAIJ); CHKERRQ(ierr);
ierr = DMCreateMatrix(dm, &J); CHKERRQ(ierr);
A = J;

ierr = DMSNESSetFunctionLocal(dm, ...); CHKERRQ(ierr);
ierr = DMSNESSetJacobianLocal(dm, ...); CHKERRQ(ierr);
ierr = SNESSetJacobian(snes, A, J, NULL, NULL); CHKERRQ(ierr);
ierr = SNESSetFromOptions(snes); CHKERRQ(ierr);

...

ierr = SNESSolve(snes, NULL, x); CHKERRQ(ierr);

...
ierr = MatDestroy(&J); CHKERRQ(ierr);


For the line "A = J;", what exactly is the difference, if any, between that
and "ierr = MatDuplicate(...)" or "ierr = MatCopy(...)"? Do these different
options somehow affect memory usage/performance? Say I am solving a
standard poisson equation using either GAMG and/or HYPRE.

Thanks,
Justin

[Jed Brown]

Post by Justin Chang
Hi all,
Mat A,J;
...
ierr = DMSetMatType(dm, MATAIJ); CHKERRQ(ierr);
ierr = DMCreateMatrix(dm, &J); CHKERRQ(ierr);
A = J;
ierr = DMSNESSetFunctionLocal(dm, ...); CHKERRQ(ierr);
ierr = DMSNESSetJacobianLocal(dm, ...); CHKERRQ(ierr);
ierr = SNESSetJacobian(snes, A, J, NULL, NULL); CHKERRQ(ierr);
ierr = SNESSetFromOptions(snes); CHKERRQ(ierr);
...
ierr = SNESSolve(snes, NULL, x); CHKERRQ(ierr);
...
ierr = MatDestroy(&J); CHKERRQ(ierr);
For the line "A = J;",

This means you have two handles referring to the same object.

Post by Justin Chang
what exactly is the difference, if any, between that and "ierr =
MatDuplicate(...)"

This creates a new object.

Post by Justin Chang
or "ierr = MatCopy(...)"?

The second argument needs to be a valid Mat to call this function.

Post by Justin Chang
Do these different options somehow affect memory usage/performance?

Yes.

Post by Justin Chang
Say I am solving a standard poisson equation using either GAMG and/or
HYPRE.
Thanks,
Justin

[Justin Chang]

So when would I use one over the other?

- If I wanted to solve a problem using a direct solver or an iterative
solver without a preconditioner, I would use A = J?

- The documentation for SNESSetJacobian() says that AMat and PMat are
usually the same, but if I used something like GAMG would I need to create
two different objects/Mats?

Thanks,
Justin

Post by Jed Brown

Post by Justin Chang
Hi all,
Mat A,J;
...
ierr = DMSetMatType(dm, MATAIJ); CHKERRQ(ierr);
ierr = DMCreateMatrix(dm, &J); CHKERRQ(ierr);
A = J;
ierr = DMSNESSetFunctionLocal(dm, ...); CHKERRQ(ierr);
ierr = DMSNESSetJacobianLocal(dm, ...); CHKERRQ(ierr);
ierr = SNESSetJacobian(snes, A, J, NULL, NULL); CHKERRQ(ierr);
ierr = SNESSetFromOptions(snes); CHKERRQ(ierr);
...
ierr = SNESSolve(snes, NULL, x); CHKERRQ(ierr);
...
ierr = MatDestroy(&J); CHKERRQ(ierr);
For the line "A = J;",

This means you have two handles referring to the same object.

Post by Justin Chang
what exactly is the difference, if any, between that and "ierr =
MatDuplicate(...)"

This creates a new object.

Post by Justin Chang
or "ierr = MatCopy(...)"?

The second argument needs to be a valid Mat to call this function.

Post by Justin Chang
Do these different options somehow affect memory usage/performance?

Yes.

Post by Justin Chang
Say I am solving a standard poisson equation using either GAMG and/or
HYPRE.
Thanks,
Justin

[Jed Brown]

Post by Justin Chang
So when would I use one over the other?
- If I wanted to solve a problem using a direct solver or an iterative
solver without a preconditioner, I would use A = J?
- The documentation for SNESSetJacobian() says that AMat and PMat are
usually the same, but if I used something like GAMG would I need to create
two different objects/Mats?

It is a semantic distinction that has nothing to do with the
preconditioning algorithm.

If you define the operator that you want to solve with using a
matrix-free implementation (could be MFFD, could be some fast evaluation
that doesn't store matrix entries), but want to use an assembled
operator for preconditioning, then you pass a different matrix. If you
use a different discretization to define the operator versus the
preconditioner, then you would pass different matrices. This happens,
for example, when preconditioning using a low-order discretization
and/or dropping some coupling terms that are deemed unimportant for
preconditioning.

[Dave May]

Post by Justin Chang
So when would I use one over the other?
- If I wanted to solve a problem using a direct solver or an iterative
solver without a preconditioner, I would use A = J?

Yes.

Post by Justin Chang
- The documentation for SNESSetJacobian() says that AMat and PMat are
usually the same, but if I used something like GAMG would I need to create
two different objects/Mats?

I would say "maybe".

If the Jacobian (here A) was defined via a matrix-free finite difference
approximation, but you wish to use a non-trivial preconditioner, you might
wish to assemble J. J might be the Picard linearized operator (for example).

Another use case where A != J might arise is if you define A with a high
order spatial discretization (probably matrix free) and you use a low order
discretization to define the preconditioner which will ultimately be passed
to GAMG.

Post by Justin Chang
Thanks,
Justin

Post by Jed Brown

Post by Justin Chang
Hi all,
Mat A,J;
...
ierr = DMSetMatType(dm, MATAIJ); CHKERRQ(ierr);
ierr = DMCreateMatrix(dm, &J); CHKERRQ(ierr);
A = J;
ierr = DMSNESSetFunctionLocal(dm, ...); CHKERRQ(ierr);
ierr = DMSNESSetJacobianLocal(dm, ...); CHKERRQ(ierr);
ierr = SNESSetJacobian(snes, A, J, NULL, NULL); CHKERRQ(ierr);
ierr = SNESSetFromOptions(snes); CHKERRQ(ierr);
...
ierr = SNESSolve(snes, NULL, x); CHKERRQ(ierr);
...
ierr = MatDestroy(&J); CHKERRQ(ierr);
For the line "A = J;",

This means you have two handles referring to the same object.

Post by Justin Chang
what exactly is the difference, if any, between that and "ierr =
MatDuplicate(...)"

This creates a new object.

Post by Justin Chang
or "ierr = MatCopy(...)"?

The second argument needs to be a valid Mat to call this function.

Post by Justin Chang
Do these different options somehow affect memory usage/performance?

Yes.

Post by Justin Chang
Say I am solving a standard poisson equation using either GAMG and/or
HYPRE.
Thanks,
Justin

[Dave May]

Damn - Jed scooped me :D

Post by Dave May

Post by Justin Chang
So when would I use one over the other?
- If I wanted to solve a problem using a direct solver or an iterative
solver without a preconditioner, I would use A = J?

Yes.

Post by Justin Chang
- The documentation for SNESSetJacobian() says that AMat and PMat are
usually the same, but if I used something like GAMG would I need to create
two different objects/Mats?

I would say "maybe".
If the Jacobian (here A) was defined via a matrix-free finite difference
approximation, but you wish to use a non-trivial preconditioner, you might
wish to assemble J. J might be the Picard linearized operator (for example).
Another use case where A != J might arise is if you define A with a high
order spatial discretization (probably matrix free) and you use a low order
discretization to define the preconditioner which will ultimately be passed
to GAMG.

Post by Justin Chang
Thanks,
Justin

Post by Jed Brown

Post by Justin Chang
Hi all,
Mat A,J;
...
ierr = DMSetMatType(dm, MATAIJ); CHKERRQ(ierr);
ierr = DMCreateMatrix(dm, &J); CHKERRQ(ierr);
A = J;
ierr = DMSNESSetFunctionLocal(dm, ...); CHKERRQ(ierr);
ierr = DMSNESSetJacobianLocal(dm, ...); CHKERRQ(ierr);
ierr = SNESSetJacobian(snes, A, J, NULL, NULL); CHKERRQ(ierr);
ierr = SNESSetFromOptions(snes); CHKERRQ(ierr);
...
ierr = SNESSolve(snes, NULL, x); CHKERRQ(ierr);
...
ierr = MatDestroy(&J); CHKERRQ(ierr);
For the line "A = J;",

This means you have two handles referring to the same object.

Post by Justin Chang
what exactly is the difference, if any, between that and "ierr =
MatDuplicate(...)"

This creates a new object.

Post by Justin Chang
or "ierr = MatCopy(...)"?

The second argument needs to be a valid Mat to call this function.

Post by Justin Chang
Do these different options somehow affect memory usage/performance?

Yes.

Post by Justin Chang
Say I am solving a standard poisson equation using either GAMG and/or
HYPRE.
Thanks,
Justin