MPIBAIJ and sor

[Barry Smith]

    I completely forgot about the MatPBRelax_ routines.

    MatRelax() which is used by PCSOR automatically tries  
MatRelax_xxx() for the matrix and then MatPBRelax_xxx()  if the  
MatRelax_xxx() does not exist.

    MatRelax_MPIAIJ() (and MatRelax_MPISBAIJ()) handles the block  
Jacobi and relax on each block automatically and allows for any number  
of local smoothing before updating the ghost points.

    If you need MatPBRelax_SeqBAIJ_N() for N > 7 that is very easily  
added by taking the 7 case and extending it.

    Barry

   This discussion brings up a few issues

1) having the MatRelax() and MatPBRelax() both is a bit confusing. As  
Jed points out it is possible to write a MatRelax() for BAIJ, but that  
is almost always silly since essentially the same number
      of floating point operations gives you the better point block  
relaxation. I originally thought we might support all point-block and  
point relaxation for all matrix types and hence have both the
      MatRelax() and MatPBRelax(). Perhaps that was unneededly general  
and we should merge them. Then I wouldn't have searched for  
MatRelax_SeqBAIJ and given you the wrong answer initially.

  2) Using SOR as the name of the preconditioner but relax as the name  
of the matrix operation is questionable. Better to use the same name  
everywhere?

  3) One could very easily modify the MatRelax_MPIAIJ() to get a  
MatRelax_MPIBAIJ() and we should do that (and test it). But as Jed  
points out it is redundant (as is MatRelax_MPIAIJ()).

On Oct 27, 2009, at 3:39 PM, Jed Brown wrote:

> Stephan Kramer wrote:
>> I'm slightly confused by the fact that if I use my BAIJ matrix in
>> parallel and select the sor preconditioner. it complains about
>> mat->ops->relax and mat->ops->pbrelax not being defined.
>
> There is no parallel SOR, it would make the algorithm nearly serial.
> Use -sub_pc_type sor with -ksp_type asm or bjacobi.
>
>> Also I can't find a MatRelax_MPIBAIJ in
>> src/mat/impls/baij/mpi/mpibaij.c as I would have expected. Is there  
>> no
>> sor type relaxation for parallel block matrices at all?
>
> Since the small dense blocks are available at once, it would be  
> silly to
> relax each component separately (your CPU spends most of its time
> waiting on memory).  Instead, BAIJ relaxes all components  
> simultaneously
> by solving with these small blocks.  This is "point-block" relaxation,
> look for MatPBRelax_SeqBAIJ_N for N in 1..7.
>
> Jed
>