This is wrong
ierr = MatCreateShell(PETSC_COMM_WORLD,ctx->high-ctx->low,ctx->high-ctx->low,\
PETSC_DETERMINE, PETSC_DETERMINE, (void*)ctx, &A);CHKERRQ(ierr);
It changes the copy of A inside this routine to a shell matrix; it doesn't
affect the A that is passed in so it is not a setup shell matrix like you need.
Instead you need to use
ierr = MatSetSizes(A,ctx->high-ctx->low,ctx->high-ctx->low,\
PETSC_DETERMINE, PETSC_DETERMINE);CHKERRQ(ierr);
ierr = MatSetType(A,MATSHELL);CHKERRQ(ierr);
ierr = MatShellSetContext(A,ctx);CHKERRQ(ierr);
ierr = MatSetUp(A);CHKERRQ(ierr);
Also in your main program immediately after creating the DMDA you must call
DMSetMatType(da,MATSHELL); otherwise it will create an AIJ matrix which you
don't want.
Barry
> On Jul 19, 2018, at 5:58 PM, Joseph Essman <[email protected]> wrote:
>
> Hello all,
>
> I'm a relatively new PETSc user and am currently attempting to apply
> matrix-free methods to a program I wrote that uses the KSP and DMDA
> libraries.
>
> However, I've come across an issue. When I use KSPSetDM, it requires that I
> set operators with KSPSetComputeOperators rather than simply setting an
> operator. However, the primary issue here is that when I use
> KSPSetComputeOperators, I attempt to compute my operators with the following
> function:
>
> PetscErrorCode ComputeMatrix(KSP ksp, Mat A, Mat jac, void *user){
> PetscErrorCode ierr;
> AppCtx *ctx = (AppCtx*)user;
>
> PetscFunctionBegin;
> ierr = MatSetType(A, MATSHELL);CHKERRQ(ierr);
> ierr =
> MatCreateShell(PETSC_COMM_WORLD,ctx->high-ctx->low,ctx->high-ctx->low,\
> PETSC_DETERMINE, PETSC_DETERMINE, (void*)ctx, &A);CHKERRQ(ierr);
> ierr = MatShellSetOperation(A,MATOP_MULT,(void(*)(void))MyMatMult);
> CHKERRQ(ierr);
> ierr = MatSetUp(A);CHKERRQ(ierr);
>
> ierr = MatSetType(jac, MATSHELL);CHKERRQ(ierr);
> ierr =
> MatCreateShell(PETSC_COMM_WORLD,ctx->high-ctx->low,ctx->high-ctx->low,\
> PETSC_DETERMINE, PETSC_DETERMINE, (void*)ctx, &jac);CHKERRQ(ierr);
> ierr = MatShellSetOperation(jac,MATOP_MULT,(void(*)(void))MyMatMult);
> CHKERRQ(ierr);
> ierr = MatSetUp(jac);CHKERRQ(ierr);
>
> PetscFunctionReturn(ierr);
> }
>
> When I try to run KSPSolve, I get the following error message:
>
> joseph@Otto ~/petsc-3.9.3/Diffusion $ mpirun -n 1 ./diffu2matfree -x 3 -y 3
> shell
> [0]PETSC ERROR: --------------------- Error Message
> --------------------------------------------------------------
> [0]PETSC ERROR: Object is in wrong state
> [0]PETSC ERROR: Must call MatXXXSetPreallocation() or MatSetUp() on argument
> 1 "mat" before MatMult()
> [0]PETSC ERROR: See http://www.mcs.anl.gov/petsc/documentation/faq.html for
> trouble shooting.
> [0]PETSC ERROR: Petsc Release Version 3.9.3, Jul, 02, 2018
> [0]PETSC ERROR: ./diffu2matfree on a arch-linux2-c-debug named Otto by joseph
> Thu Jul 19 17:53:00 2018
> [0]PETSC ERROR: Configure options --with-cc=gcc --with-cxx=g++
> --with-fc=gfortran --download-mpich --download-fblaslapack
> [0]PETSC ERROR: #1 MatMult() line 2306 in
> /home/joseph/petsc-3.9.3/src/mat/interface/matrix.c
> [0]PETSC ERROR: #2 PCApplyBAorAB() line 682 in
> /home/joseph/petsc-3.9.3/src/ksp/pc/interface/precon.c
> [0]PETSC ERROR: #3 KSP_PCApplyBAorAB() line 304 in
> /home/joseph/petsc-3.9.3/include/petsc/private/kspimpl.h
> [0]PETSC ERROR: #4 KSPGMRESCycle() line 152 in
> /home/joseph/petsc-3.9.3/src/ksp/ksp/impls/gmres/gmres.c
> [0]PETSC ERROR: #5 KSPSolve_GMRES() line 234 in
> /home/joseph/petsc-3.9.3/src/ksp/ksp/impls/gmres/gmres.c
> [0]PETSC ERROR: #6 KSPSolve() line 669 in
> /home/joseph/petsc-3.9.3/src/ksp/ksp/interface/itfunc.c
> [0]PETSC ERROR: #7 main() line 110 in
> /home/joseph/petsc-3.9.3/Diffusion/diffu2matfree.cpp
> [0]PETSC ERROR: PETSc Option Table entries:
> [0]PETSC ERROR: -x 3
> [0]PETSC ERROR: -y 3
> [0]PETSC ERROR: ----------------End of Error Message -------send entire error
> message to [email protected]
> application called MPI_Abort(MPI_COMM_WORLD, 73) - process 0
>
>
> Does anybody know what's going on here? I've also tried to make the algorithm
> work without KSPSetDM and KSPSetComputeOperators, but my MatMult step is
> dependent on DMDA, which is no longer used for the vectors in the KSP if I
> don't use KSPSetDM. I've attached my code in full if there are any questions
> as to the structure of my code.
>
> Thanks,
>
> Joseph Essman
> <diffu2matfree.cpp>
