Just pass mJac, mJac instead of mJac, nullptr and it will be happy. In
your case, the second mJac won't be used in your preconditioner it is just a
place holder so other parts of SNES won't try to create a matrix.
Barry
> On Nov 6, 2024, at 11:36 AM, Weide, Edwin van der (UT-ET) via petsc-users
> <petsc-users@mcs.anl.gov> wrote:
>
> 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
