inverting independent, small dense matrices

[Barry Smith]

On Feb 13, 2009, at 3:01 PM, Shao-Ching Huang wrote:

> Hi
>
> I need to invert a bunch of independent, dense 4x4 matrices on each
> MPI process. (This is separate from the main computation, a parallel
> Ax=b solve.) I am thinking of using LAPACK's DGESV to invert these 4x4
> matrices.
>

      You would not want to use LAPACK's solvers for such tiny  
matrices. LAPACK is optimized
for large dense matrices; the overhead is high for small matrices so  
that actual performance will
be very slow calling LAPACK.

> In my code, is it necessary to call LAPACK directly? Or does PETSc
> already have a convenient interface to LAPACK?
>
      For PETSc's own internal use we have written custom inverters  
for small matrices (these are used by
the factorization routines for BAIJ matrices).  They are called  
Kernel_A_gets_inverse_A_* where * is the block size.
The prototypes for the functions are in src/inline/ilu.h and the  
actual source code in files in src/mat/impls/baij/seq/dgef*.c

    Barry


> The size 4x4 may be increased to 8x8 or so later, so I am not trying
> to use the closed-form inversion formula.
>
> Thanks,
>
> Shao-Ching