> On Jan 7, 2016, at 7:58 PM, Timothée Nicolas <[email protected]>
> 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 <[email protected]>:
>
> 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 <[email protected]>
> > 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 <[email protected]>:
> > On Thu, Jan 7, 2016 at 7:49 AM, Timothée Nicolas
> > <[email protected]> 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
> >
>
>
