Hi Dave Thanks for the response. I'm actually using fortran and I wasn't sure that PetscObjectComposeFunction would be accessible, and if so, what sort of fortran magic I might need to call this function (possibly an interface, with or without c_opt).
Do you know if it is possible to call that routine directly from fortran? Thanks Steven On Fri, Aug 26, 2016 at 1:35 PM, Dave May <[email protected]> wrote: > > > On 26 August 2016 at 14:34, Dave May <[email protected]> wrote: > >> >> >> On 26 August 2016 at 14:14, Steven Dargaville < >> [email protected]> wrote: >> >>> Hi all >>> >>> I'm just wondering if there is any plans in the future for >>> MatGetDiagonalBlock to support shell matrices by registering a >>> user-implemented MATOP? MatGetDiagonal supports MATOP, but the block >>> version of this does not. >>> >>> I found a previous query on the user list which touched on this and >>> mentioned that it would be easy to add: >>> >>> http://lists.mcs.anl.gov/pipermail/petsc-users/2011-May/008700.html >>> >>> I have implemented a matrix-free multigrid algorithm using shell >>> operations in PETSc, and it would be very convenient to be able to provide >>> a local shell Mat so I could apply local GMRES (or other matvec-based >>> solvers) as a local block smoother. >>> >> >> It looks like the thing you need to do is use >> PetscObjectComposeFunction() and not MatShellSetOperation() >> >> http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/ >> Sys/PetscObjectComposeFunction.html#PetscObjectComposeFunction >> >> with your MatShell object. >> >> That is, instead of calling MatShellSetOperation(), call >> ierr = PetscObjectComposeFunction(myshell,"MatGetDiagonalBlock_C", >> MatGetDiagonalBlock_MyShell);CHKERRQ(ierr); >> > > Oops - I forgot the cast, the above should be > > ierr = PetscObjectComposeFunction((PetscObject)myshell,"MatGetDia > gonalBlock_C", MatGetDiagonalBlock_MyShell);CHKERRQ(ierr); > > >> >> where MatGetDiagonalBlock_MyShell is a function pointer to your method to >> get the diagonal block. >> >> Important detail is that you don't change the string >> "MatGetDiagonalBlock_C". The method MatGetDiagonalBlock() does a function >> pointer query using this string. See >> http://www.mcs.anl.gov/petsc/petsc-current/src/mat/interface >> /matrix.c.html#MatGetDiagonalBlock >> >> >> Thanks >> Dave >> >> >> >> >> >> >>> >>> Thanks! >>> Steven >>> >>> >>> >> >
