Matrix free preconditioning

[Tim Kroeger]

Dear Jed and Matt and all,

On Tue, 14 Oct 2008, Jed Brown wrote:

> I assume your matrix comes from the rank-1 update you brought up
> recently on the libmesh list?  That is
>
>  B = A + u v'
>
> for where A is an explicit matrix and u,v are vectors.

Yes, that's it.  Actually, I didn't know that anyone was reading both 
libmesh-users and petsc-users.

> Did you try the
> rank-1 update I suggested?  In detail, let A^ be an approximate inverse
> of A (i.e. an application of any standard preconditioner) and form a
> preconditioner for B as
>
>  B^ = A^ - alpha w v' A^
>
> where
>
>  w = A^ u
>  alpha = 1/(1 + v' w).

Thank you very much for that suggestion.  I guess this would be the 
best thing to do, but actually I have to admit that I am rather shy of 
implemeting that, in particular with a sensible interface to libMesh.

Hence, I am thinking whether there might be some easier (altough 
certainly less efficient) possibility.  I am thinking abouth the 
following two ideas:

1. I could do MatShellSetOperation(A,MATOP_GET_DIAGONAL,...), which is 
easy in my case and should enable PETSc to do at least JACOBI 
preconditioning for that matrix, shouldn't it?

2. Following Matt's suggestion, I could handle an auxilary matrix C 
for preconditioning.  In my case, it would choose C(i,j)=A(i,j) 
whenever B(i,j)!=0, but C(i,j)=0 else.  I.e., C has the same sparsity 
pattern as the (sparse) matrix B.  This would also be easy to 
construct in my case.

Any opinion on this is welcome.

Best Regards,

Tim

-- 
Dr. Tim Kroeger                                        Phone +49-421-218-7710
tim.kroeger at mevis.de, tim.kroeger at cevis.uni-bremen.de  Fax   +49-421-218-4236

MeVis Research GmbH, Universitaetsallee 29, 28359 Bremen, Germany

Amtsgericht Bremen HRB 16222
Geschaeftsfuehrer: Prof. Dr. H.-O. Peitgen