some sor questions

[Stephan Kramer]

Thanks for your answers

Barry Smith wrote:
> On Sep 22, 2009, at 8:47 AM, Stephan Kramer wrote:
> 
>> Hi all,
>>
>> I have some questions basically about the MatRelax_SeqAIJ routine:
>>
>> If I understand correctly there are 2 versions of the sor routine  
>> depending on whether or not there is a zero guess, so that with a  
>> zero guess in the forward sweep you don't need to multiply the upper  
>> diagonal part U with the part of the x vector that is still zero.  
>> Why then does it look like that both versions log the same number of  
>> flops? I would have expected that the full forward sweep (i.e. no  
>> zero guess) takes 2*a->nz flops (i.e. the same as a matvec) and not  
>> a->nz.
> 
>     You are right. This is an error in our code. It will be in the  
> next patch.
> 
>> Why does the Richardson iteration with sor not take the zero guess  
>> into account, i.e. why does PCApplyRichardson_SOR not set  
>> SOR_ZERO_INIT_GUESS in the call to MatRelax if the Richardson ksp  
>> has a zero initial guess set?
> 
>     This appears to be a design limitation. There is no mechanism to  
> pass the information that the initial guess is zero into  
> PCApplyRichardson(). We could add support for this by adding one more  
> argument to PCApplyRichardson() for this information. I don't see a  
> simpler way.  If one is running, say 2 iterations of Richardson then  
> this would be a measurable improvement in time. If one is running many  
> iterations then the savings is tiny. Perhaps this support should be  
> added.

I'm thinking of the application of ksp richardson with sor as a smoother 
in pcmg. In which case the down smoother will have zero initial guess 
(as it only acts on the residual), and there will be typicaly only 1 or 
2 iterations, so the saving would be significant. Is there another way I 
should set this up instead?

> 
> 
> 
>> In parallel if you specify SOR_LOCAL_FORWARD_SWEEP or  
>> SOR_LOCAL_BACKWARD_SWEEP it
> 
> 
>> calls MatRelax on the local part of the matrix, mat->A, with  
>> its=lits and lits=PETSC_NULL (?).
> 
>> However the first line of MatRelax_SeqAIJ then says: its = its*lits.  
>> Is that right?
> 
>     This is all wrong. It should be passing 1 in, not PETSC_NULL. This  
> was fixed in petsc-dev but not in petsc-3.0.0 I will fix it in  
> petsc-3.0.0 and it will be in the next patch.
> 
>     Thanks for pointing out the problems.
> 
>     If you plan to use SOR a lot, you might consider switching to  
> petsc-dev since I have made some improvements there. Also consider the  
> Eisenstat trick preconditioner.
> 
>     Barry
> 
> 
>> Please tell me if I'm totally misunderstanding how the routine  
>> works, thanks for any help.
>>
>> Cheers
>> Stephan
>>
>> -- 
>> Stephan Kramer  <s.kramer at imperial.ac.uk>
>> Applied Modelling and Computation Group,
>> Department of Earth Science and Engineering,
>> Imperial College London
> 
>