How to implement a block version of MatDiagonalScale?

Stephane Aubert stephane.aubert at
Tue May 6 12:20:07 CDT 2008

I'm "playing" with badly conditioned block matrices (MPIBAIJ type) 
arising from turbulent compressible fluid dynamics (RANS eqs, block size 
= 5 (laminar) or 7 (turbulent)).
I tested that using bi-normalization approach to build l and r vectors, 
then calling MatDiagonalScale(A,l,r), improve the accuracy and the 
convergence of GMRES+ILU(0), basically by reducing the dependence to the 
mesh cells size and to the variables magnitudes (in particular for k-w 
turbulence model).
Now, I would like to go further and use L and R block diagonal matrices, 
instead of vectors, for example to get block identity along the diagonal 
of L.A.R. The sparcity of A is preserved, from a block point of view.
My first idea is to use a LU factorization of Aii, with L=Li^(-1), 
R=Ui^(-1), Aii=Li.Ui.

My question is: How to compute L and R using PETSC available functions 
in a clever way, instead of calling some LAPACK functions after getting 
individual diagonal blocks in my own piece of code?

My actual guess is:

   1. Create a new matrix D with MPIBDIAG type, from the block-diagonal
      of A using MatGetSubMatrices(), to do something like MatGetDiagonal().
   2. Create a PC of type PCLU using D as operators, with the sequence
   3. Get Li and Ui. This is where I'm glued: I can't figure out how to
      use PCGetFactoredMatrix() to get Li and Ui... Is MatLUFactor() a
      better candidate?
   4. Do I need to compute Li^(-1) and Ui^(-1), or is there some
      "backsubstitution" functions available?

With many thanks for your suggestions,

Dr. Stephane AUBERT, CEO & CTO
Centre Scientifique Auguste MOIROUX
64 chemin des MOUILLES
International:      fax: +33      tel: +33
France:             fax:         tel:
email: stephane.aubert at

More information about the petsc-users mailing list