[petsc-users] Block Jacobi for Matrix-free

[Barry Smith]

> On Jan 7, 2016, at 8:31 PM, Timothée Nicolas <timothee.nicolas at gmail.com> wrote:
> 
> Ah, I understand, so by allocating this BAIJ in an intelligent way (allocating only the diagonal 3x3 blocks), I can still be basically memory efficient, and use matrix-free formulation for the first matrix in KSPSetOperator, right ?

  Exactly

> 
> Timothee
> 
> 2016-01-08 11:06 GMT+09:00 Barry Smith <bsmith at mcs.anl.gov>:
> 
> > On Jan 7, 2016, at 7:58 PM, Timothée Nicolas <timothee.nicolas at gmail.com> wrote:
> >
> > I see, I just tested PCPBJACOBI, which is better than PCJACOBI, but still I may need PCBJACOBI.
> 
>    Note that using PCBJACOBI means you are providing big blocks of the Jacobian. If you do provide big blocks of the Jacobian you might as well just provide the entire Jacobin IMHO.
> 
>    Anyways the easiest way to do either PCPBJACOBI or PCBJACOBI is to explicitly construct the portion of the Jacobian you need, in a AIJ or BAIJ matrix and pass that as the SECOND matrix argument to KSPSetOperator() or SNESSetJacobian() then PETSc will use the piece you provide to build the preconditioner. So for example if you want PBJACOBI you would create a BAIJ matrix with block size 3 and only fill up the 3 by 3 block diagonal with Jacobian entries.
> 
> 
>    Barry
> 
> 
> > The problem is that I don't seem to be allowed to define the matrix operation for MatGetDiagonalBlock... Indeed, I don't find
> >
> > MATOP_GET_DIAGONAL_BLOCK in ${PETSC_DIR}/include/petscmat.h
> >
> > Therefore, when I try to define it, I get the following error at compilation (quite logically)
> >
> > matrices.F90(174): error #6404: This name does not have a type, and must have an explicit type.   [MATOP_GET_DIAGONAL_BLOCK]
> >      call MatShellSetOperation(lctx(1)%PSmat,MATOP_GET_DIAGONAL_BLOCK,PSmatGetDiagonalBlock,ierr)
> > ---------------------------------------------^
> >
> > Also, if I change my mind and instead decide to go for PCPBJACOBI, I still have a problem because the manual says that the routine you talk about, MatInvertBlockDiagonal, is not available from FORTRAN. Indeed I cannot call it. I still cannot call it after I provide a routine corresponding to MATOP_INVERT_BLOCK_DIAGONAL.
> >
> > So, it seems to mean that if I want to use this kind of algorithms, I will have to hard code them, which would be too bad. Is that right, or is there an other way around these two issues ?
> >
> > Best
> >
> > Timothee
> >
> >
> >
> >
> > 2016-01-08 2:38 GMT+09:00 Barry Smith <bsmith at mcs.anl.gov>:
> >
> >   Timothee,
> >
> >       You are mixing up block Jacobi PCBJACOBI (which in PETSc generally uses "big" blocks) and point block Jacobi PCPBJACOBI (which generally means all the degrees of freedom associated with a single grid point -- in your case 3).
> >
> >    If you are doing matrix free with a shell matrix then you need to provide your own MatInvertBlockDiagonal() which in your case would invert each of your little 3 by 3 blocks and store the result in a 1d array; each little block in column major order followed by the next one. See for example MatInvertBlockDiagonal_SeqAIJ(). You also need you matrix to return a block size of 3.
> >
> >
> >    Barry
> >
> >
> >
> > > On Jan 7, 2016, at 8:08 AM, Timothée Nicolas <timothee.nicolas at gmail.com> wrote:
> > >
> > > Ok, so it should be sufficient. Great, I think I can do it.
> > >
> > > Best
> > >
> > > Timothée
> > >
> > > 2016-01-07 23:06 GMT+09:00 Matthew Knepley <knepley at gmail.com>:
> > > On Thu, Jan 7, 2016 at 7:49 AM, Timothée Nicolas <timothee.nicolas at gmail.com> wrote:
> > > Hello everyone,
> > >
> > > I have discovered that I need to use Block Jacobi, rather than Jacobi, as a preconditioner/smoother. The linear problem I am solving at this stage lives in a subspace with 3 degrees of freedom, which represent the 3 components of a 3D vector. In particular for multigrid, using BJACOBI instead of JACOBI as a smoother changes everything in terms of efficiency. I know it because I have tested with the actual matrix in matrix format for my problem. However, eventually, I want to be matrix free.
> > >
> > > My question is, what are the operations I need to provide for the matrix-free approach to accept BJACOBI ? I am confused because when I try to apply BJACOBI to my matrix-free operator; the code asks for MatGetDiagonalBlock (see error below). But MatGetDiagonalBlock, in my understanding, returns a uniprocessor matrix representing the diagonal part of the matrix on this processor (as defined in the manual). Instead, I would expect that what is needed is a routine which returns a 3x3 matrix at the grid point (that is, the block associated with this grid point, coupling the 3 components of the vector together). How does this work ? Do I simply need to code MatGetDiagonalBlock ?
> > >
> > > Just like Jacobi does not request one diagonal element at a time, Block-Jacobi does not request one diagonal block at a time. You
> > > would need to implement that function, or write a custom block Jacobi for this matrix.
> > >
> > >   Thanks,
> > >
> > >     Matt
> > >
> > >
> > > Thx
> > > Best
> > >
> > > Timothée
> > >
> > > [0]PETSC ERROR: --------------------- Error Message --------------------------------------------------------------
> > > [0]PETSC ERROR: No support for this operation for this object type
> > > [0]PETSC ERROR: Matrix type shell does not support getting diagonal block
> > > [0]PETSC ERROR: See http://www.mcs.anl.gov/petsc/documentation/faq.html for trouble shooting.
> > > [0]PETSC ERROR: Petsc Release Version 3.6.1, Jul, 22, 2015
> > > [0]PETSC ERROR: ./miips on a arch-linux2-c-debug named Carl-9000 by timothee Thu Jan  7 22:41:13 2016
> > > [0]PETSC ERROR: Configure options --with-cc=gcc --with-cxx=g++ --with-fc=gfortran --download-fblaslapack --download-mpich
> > > [0]PETSC ERROR: #1 MatGetDiagonalBlock() line 166 in /home/timothee/Documents/petsc-3.6.1/src/mat/interface/matrix.c
> > > [0]PETSC ERROR: #2 PCSetUp_BJacobi() line 126 in /home/timothee/Documents/petsc-3.6.1/src/ksp/pc/impls/bjacobi/bjacobi.c
> > > [0]PETSC ERROR: #3 PCSetUp() line 982 in /home/timothee/Documents/petsc-3.6.1/src/ksp/pc/interface/precon.c
> > > [0]PETSC ERROR: #4 KSPSetUp() line 332 in /home/timothee/Documents/petsc-3.6.1/src/ksp/ksp/interface/itfunc.c
> > > [0]PETSC ERROR: #5 KSPSolve() line 546 in /home/timothee/Documents/petsc-3.6.1/src/ksp/ksp/interface/itfunc.c
> > >
> > >
> > >
> > > --
> > > What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead.
> > > -- Norbert Wiener
> > >
> >
> >
> 
>