# How to implement a block version of MatDiagonalScale?

Barry Smith bsmith at mcs.anl.gov
Tue May 6 15:57:35 CDT 2008

```   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
>
>

```