[petsc-users] PCMGSetResidual and fortran

[Barry Smith]

  Steven,

    Could you please send your test code. It is possible there is a bug in our Fortran interface since we do not test it for this functionality. 

  Barry


> On Nov 4, 2022, at 10:50 AM, Steven Dargaville <dargaville.steven at gmail.com> wrote:
> 
> Hi all
> 
> I have a quick question regarding the use of PCMGSetResidual within fortran code. I'm calling PCMGSetResidual from within fortran: 
> 
> call PCMGSetResidual(pc_mg, petsc_level, mg_residual, coarse_matrix ierr) 
> 
> and just for testing purposes I've written a trivial residual evaluation routine:
> 
>     subroutine mg_residual(mat, b, x, r, ierr)
>       !!< Compute the residual
> 
>       ! ~~~~~~
>       type(tMat) :: mat
>       type(tVec) :: b, x, r
>       PetscErrorCode :: ierr
>       ! ~~~~~~
> 
>       print *, "inside residual evaluation"
>       call MatResidual(mat, b, x, r, ierr)
>       
>    end subroutine
> 
> The problem I am having is that this segfaults when the residual routine is called. Valgrind shows that it is failing in the fortran interface in ftn-custom/zmgfuncf.c, with the message:
> 
> ==24742== Invalid read of size 8
> ==24742==    at 0x5B7CBC0: ourresidualfunction (in /home/projects/dependencies/petsc_main/arch-linux-c-opt/lib/libpetsc.so.3.015.0)
> ==24742==    by 0x5B6D804: PCMGMCycle_Private (in /home/projects/dependencies/petsc_main/arch-linux-c-opt/lib/libpetsc.so.3.015.0)
> 
> ==24742== Process terminating with default action of signal 11 (SIGSEGV)
> ==24742==  Access not within mapped region at address 0x0
> ==24742==    at 0x5B7CBC0: ourresidualfunction (in /home/sdargavi/projects/dependencies/petsc_main/arch-linux-c-opt/lib/libpetsc.so.3.015.0)
> ==24742==    by 0x5B6D804: PCMGMCycle_Private (in /home/sdargavi/projects/dependencies/petsc_main/arch-linux-c-opt/lib/libpetsc.so.3.015.0)
> 
> I'm guessing this is because the fortran_func_pointers isn't pointing to the mg_residual routine, but I am not sure why. I noticed that in the C code of PCMGSetResidual that it calls MatDestroy on the A matrix in mg_levels and replaces it with the mat passed in:
> 
> if (mat) PetscObjectReference((PetscObject)mat);
> MatDestroy(&mglevels[l]->A);
> mglevels[l]->A = mat
> 
> so I modified my code to call PCMGSetResidual either before the operators are set, or after but passing in an extra copy, but this doesn't seem to help.
> 
> I'm guessing I'm doing something silly, but just wondering if anyone had any ideas? Thanks for your help
> Steven