[petsc-users] Block Jacobi for Matrix-free
Barry Smith
bsmith at mcs.anl.gov
Thu Jan 7 11:38:58 CST 2016
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
>
More information about the petsc-users
mailing list