Question on a matrix with a particular structure

[sihv]

Dear LAPACK Forum

We have just implemented LAPACK into our software for solving our matrices and it has speed up our processing from 3 hours to 50minutes! I would really like to thank you for that! This will save us a lot of CPU time and also wait time. I only wish I had found out about this many many years ago!!!

Besides these thanks I would also like to take this opportunity for a small question. I need to invert (or at least solve) a symmetrical positive definite matrix A of the following form:

Code: Select all

A = ( a11 a12 )   
    ( a21 a22 )


where:
- a11 is a square (n by n) full (more or less) matrix, and
- a22 is square (m by m, m<n) (band) diagonal matrix, and
- a12 (and its transposed a21) is a (more or less) full matrix.


Is there any particular routine/driver in LAPACK to invert this matrix
more smartly than with the DPOTRx routines which I am using right now!? My Linear Algebra classes are a bit long in the past to remember the details which makes it difficult to understand (parts of) the documentation. I have been and am reading the user guide but did not find something applicable for this. Can you help out a bit!?

Many thanks in advance,
Tim

[Julien Langou]

Hello,

The system you describe is generally called 'Saddle Point Systems'. In general a22 is 0, but it is happens that it is nonzero, as in your case.

So first 'no' lapack as nothing to deal directly and efficiently with those kind of matrices.

From the lapack point of view, a solution I'll see might be to try to form
the schur complement then to do the factorization of the Shur complement.

I'll assume that
- a11 is full rank
- you want to solve a linear system
- n1 (size of a11) is bigger than n2 (size of a22) (not necessary, this is better)

The idea of the Schur complement is basically to solve a 2x2 linear systems. If you want to solve

Code: Select all

 a11 * x1 + a12 * x2 = b1
 a21 * x2 + a22 * x2 = b2


then you'll first find x2 via the equation

Code: Select all

(- a21 * a11^{-1} * a12 + a22) * x2 = - a21 * a11^{-1} * b1 + b2


with (- a21 * a11^{-1} * a12 + a22) that is symmetric (and it shall be positive definite as well) of size n2-by-n2. (The matrix (- a21 * a11^{-1} * a12 + a22) is called Schur complement.)

If one of the two dimensions is really biggest than the other one, this is really faster. However it will requires a bit of coding.

Julien Langou

[Julien Langou]

Jim Demmel wrote:

Just a small addition to Julien's good reply:

(1) If n2 is bigger than n1, you can just reorder the equations to swap the positions of the a11 and a22 blocks before applying the formula. This will be faster.

(2) In the formula, this is probably clear to you, but the multiplication by a11^(-1) is just multiplication by a diagonal matrix and very cheap to do by a pair of nested loops.

(3) To form a22 - a21 * a11^{-1} * a12 it may be cheaper to form X = a11^(-1/2)*a12 (again just multiplying each entry of a12 by a number) call the BLAS routine xSYRK to form a22 - X^T*X

Depending on the size of your matrix and the BLAS available, this may be faster.


Jim Demmel