MPIBAIJ and sor

[Barry Smith]

    I have added the support for SOR with MPIBAIJ matrices  
(MatSOR_MPIBAIJ()) to petsc-dev. (Note as in the discussion this is  
actually block Jacobi with SOR on each process).

    I have changed the name of MatRelax() to MatSOR() in petsc-dev so  
that a consistent name is used, instead of sometimes using relax.

    I have also merged the MatPBRelax() into MatRelax() in petsc-dev  
(now MatSOR()) to simplify and clarify the code.

    Barry


On Oct 27, 2009, at 4:04 PM, Barry Smith wrote:

>
>   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
>>
>