How to implement a block version of MatDiagonalScale?

[Barry Smith]

   Stef,

     You do not want to do it either way you have outlined below.  
Since the blocks are of size 5 or 7 you DO NOT EVER IN A BILLION YEARS  
WANT
to use LAPACK/BLAS to do the little factorizations and solves. The  
overhead of the LAPACK/BLAS will kill performance for that size  
array.  Similarly
using the PETSc Mat objects  and PC objects are not suitable for those  
tiny matrices.

    Here is what I would do. Take a look at the function  
MatLUFactorNumeric_SeqBAIJ_5() in the file src/mat/impls/baij/seq/ 
baijfact9.c Note it uses
Kernel_A_gets_inverse_A_5() to invert the little 5 by 5 blocks. Also  
it uses inlined code to do the little 5 by 5 matrix matrix  
multiplies.  (Yes for this
size problem you do want to invert the little 5 by 5 matrices and do  
matrix matrix products instead of only doing 5 by 5 LU factorizations  
and lots
of triangular solves; it is much faster this way.) I think by looking  
at this subroutine you can figure out how to loop over the nonzero  
blocks of A
and multiply by the appropriate diagonal blocks you have inverted to  
obtain the new block diagonal scaled A matrix.

    Barry

We would like to include the resulting code in PETSc if you would like  
to donate it. Thanks




On May 6, 2008, at 12:20 PM, Stephane Aubert wrote:

> Hello,
> 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
>     PCCreate/PCSetType/PCSetOperators/PCSetUp.
>  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,
> Stef.

>
>
> -- 
> ___________________________________________________________
> Dr. Stephane AUBERT, CEO & CTO
> FLUOREM s.a.s
> Centre Scientifique Auguste MOIROUX
> 64 chemin des MOUILLES
> F-69130 ECULLY, FRANCE
> International:      fax: +33 4.78.33.99.39      tel: +33 4.78.33.99.35
> France:             fax: 04.78.33.99.39         tel: 04.78.33.99.35
> email: stephane.aubert at fluorem.com
> web: www.fluorem.com
>
>