Hi,
I am trying to solve a nonlinear problem with matrix-free SNES where I would
like to provide both the matrix vector product and the preconditioner myself.
For that purpose, I use the following construction.
// Set up the matrix free evaluation of the Jacobian times a vector
// by setting the appropriate function in snes.
PetscCall(MatCreateShell(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE,
nEqns, nEqns, this, &mJac));
PetscCall(MatShellSetOperation(mJac, MATOP_MULT,
(void (*)(void))JacobianTimesVector));
PetscCall(SNESSetJacobian(mSnes, mJac, nullptr, nullptr, nullptr));
// Set the function to be used as preconditioner for the krylov solver.
KSP ksp;
PC pc;
PetscCall(SNESGetKSP(mSnes, &ksp));
PetscCall(KSPGetPC(ksp, &pc));
PetscCall(PCSetType(pc, PCSHELL));
PetscCall(PCSetApplicationContext(pc, this));
PetscCall(PCShellSetApply(pc, Preconditioner));
For small problems this construction works, and it does exactly what I expect
it to do. However, when I increase the problem size, I get a memory allocation
failure in SNESSolve, because it looks like SNES attempts to allocate memory
for a full dense matrix for the preconditioner, which is not used. This is the
call stack when the error occurs.
[0]PETSC ERROR: #1 PetscMallocAlign() at
/home/vdweide/petsc/src/sys/memory/mal.c:53
[0]PETSC ERROR: #2 PetscTrMallocDefault() at
/home/vdweide/petsc/src/sys/memory/mtr.c:175
[0]PETSC ERROR: #3 PetscMallocA() at
/home/vdweide/petsc/src/sys/memory/mal.c:421
[0]PETSC ERROR: #4 MatSeqDenseSetPreallocation_SeqDense() at
/home/vdweide/petsc/src/mat/impls/dense/seq/dense.c:3357
[0]PETSC ERROR: #5 MatSeqDenseSetPreallocation() at
/home/vdweide/petsc/src/mat/impls/dense/seq/dense.c:3338
[0]PETSC ERROR: #6 MatDuplicateNoCreate_SeqDense() at
/home/vdweide/petsc/src/mat/impls/dense/seq/dense.c:372
[0]PETSC ERROR: #7 MatDuplicate_SeqDense() at
/home/vdweide/petsc/src/mat/impls/dense/seq/dense.c:399
[0]PETSC ERROR: #8 MatDuplicate() at
/home/vdweide/petsc/src/mat/interface/matrix.c:4964
[0]PETSC ERROR: #9 DMCreateMatrix_Shell() at
/home/vdweide/petsc/src/dm/impls/shell/dmshell.c:195
[0]PETSC ERROR: #10 DMCreateMatrix() at
/home/vdweide/petsc/src/dm/interface/dm.c:1501
[0]PETSC ERROR: #11 SNESSetUpMatrices() at
/home/vdweide/petsc/src/snes/interface/snes.c:794
[0]PETSC ERROR: #12 SNESSetUp_NEWTONLS() at
/home/vdweide/petsc/src/snes/impls/ls/ls.c:290
[0]PETSC ERROR: #13 SNESSetUp() at
/home/vdweide/petsc/src/snes/interface/snes.c:3395
[0]PETSC ERROR: #14 SNESSolve() at
/home/vdweide/petsc/src/snes/interface/snes.c:4831
[0]PETSC ERROR: #15 SolveCurrentStage() at SolverClass.cpp:502
In the function SNESSetUpMatrices the source looks as follows
784 } else if (!snes->jacobian_pre) {
785 PetscDS prob;
786 Mat J, B;
787 PetscBool hasPrec = PETSC_FALSE;
788
789 J = snes->jacobian;
790 PetscCall(DMGetDS(dm, &prob));
791 if (prob) PetscCall(PetscDSHasJacobianPreconditioner(prob, &hasPrec));
792 if (J) PetscCall(PetscObjectReference((PetscObject)J));
793 else if (hasPrec) PetscCall(DMCreateMatrix(snes->dm, &J));
794 PetscCall(DMCreateMatrix(snes->dm, &B));
795 PetscCall(SNESSetJacobian(snes, J ? J : B, B, NULL, NULL));
796 PetscCall(MatDestroy(&J));
797 PetscCall(MatDestroy(&B));
798 }
It looks like in line 794 it is attempted to create the preconditioner, because
it was (intentionally) not provided.
Hence my question. Is it possible to use matrix-free SNES with a user provided
matrix vector product (via MatShell) and a user provided preconditioner
operation without SNES allocating the memory for a dense matrix? If so, what do
I need to change in the construction above to make it work?
If needed, I can provide the source code for which this problem occurs.
Thanks,
Edwin
---------------------------------------------------
Edwin van der Weide
Department of Mechanical Engineering
University of Twente
Enschede, the Netherlands
