[petsc-users] Using both DMDA and matrix-free methods with KSP

[Smith, Barry F.]

   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 <jte1 at rice.edu> 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 at 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 petsc-maint at mcs.anl.gov----------
> 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>