Matrix free preconditioning

Tim Kroeger tim.kroeger at
Wed Oct 15 01:49:28 CDT 2008

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,


Dr. Tim Kroeger                                        Phone +49-421-218-7710
tim.kroeger at, tim.kroeger at  Fax   +49-421-218-4236

MeVis Research GmbH, Universitaetsallee 29, 28359 Bremen, Germany

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

More information about the petsc-users mailing list