Inverting the following positive semidefinite matrix

[Olumide]

I'm trying to invert the following positive semidefinite matrix, that
has zeros along its diagonal. Matrices of this sort that arise from
scattered data interpolation with radial basis functions.

[ 0.000000 529.831737 1198.292909 529.831737 1.000000
10.000000 0.000000 ]
[ 529.831737 0.000000 529.831737 1198.292909 1.000000
0.000000 10.000000 ]
[ 1198.292909 529.831737 0.000000 529.831737 1.000000
-10.000000 0.000000 ]
[ 529.831737 1198.292909 529.831737 0.000000 1.000000
0.000000 -10.000000 ]
[ 1.000000 1.000000 1.000000 1.000000 0.000000
0.000000 0.000000 ]
[ 10.000000 0.000000 -10.000000 0.000000 0.000000
0.000000 0.000000 ]
[ 0.000000 10.000000 0.000000 -10.000000 0.000000
0.000000 0.000000 ]

The LAPACK routines DPOTRI and DGETRI fail because the matrix has zero
elements along its diagonal. However, GNU Octave successfully computes
the inverts the matrix to be

[ 1.8034e-003 -1.8034e-003 1.8034e-003 -1.8034e-003
2.5000e-001 5.0000e-002 -2.2854e-018 ]
[ -1.8034e-003 1.8034e-003 -1.8034e-003 1.8034e-003
2.5000e-001 -1.4014e-017 5.0000e-002 ]
[ 1.8034e-003 -1.8034e-003 1.8034e-003 -1.8034e-003
2.5000e-001 -5.0000e-002 -1.1356e-017 ]
[ -1.8034e-003 1.8034e-003 -1.8034e-003 1.8034e-003
2.5000e-001 -1.2975e-017 -5.0000e-002 ]
[ 2.5000e-001 2.5000e-001 2.5000e-001 2.5000e-001 -5.6449e
+002 1.6031e-015 -2.2001e-015 ]
[ 5.0000e-002 1.5053e-017 -5.0000e-002 1.1470e-017
-1.6031e-015 5.9915e+000 2.3211e-016 ]
[ 1.3878e-017 5.0000e-002 5.8566e-018 -5.0000e-002
1.4668e-015 -1.6985e-017 5.9915e+000 ]

I've verified that the product of both matrices is the identify
matrix. Unfortunately, there is no way to integrate Octave into my
standalone application. As such I'm still looking for a way to compute
the inverse of the above matrix, preferably using LAPACK .

Thanks,

- Olumide

[Mark Hoemmen]

On Mon, Apr 4, 2011 at 07:38, Olumide <50...@web.de> wrote:
> I'm trying to invert the following  positive semidefinite matrix, that
> has zeros along its diagonal. Matrices of this sort that arise from
> scattered data interpolation with radial basis functions.
>

> ...

> The LAPACK routines DPOTRI and DGETRI fail because the matrix has zero
> elements along its diagonal.

Have you tried DSYTRI?

http://www.netlib.org/lapack/double/dsytri.f

mfh

[Bill Greene]

I called DGETRF with your matrix as input followed by DGETRI and it
returns the inverse you show below. Did you call DGETRF before
your call to DGETRI?

--
You received this message because you are subscribed to the Google Groups "matrixprogramming" group.
To post to this group, send email to matrixpr...@googlegroups.com.
To unsubscribe from this group, send email to matrixprogramm...@googlegroups.com.
For more options, visit this group at http://groups.google.com/group/matrixprogramming?hl=en.

[Bill Greene]

As Mark Hoemmen pointed out, DSYTRI is the better choice in this case because
the input matrix is symmetric. If works fine, too, but in that case you have to
remember to call DSYTRF first.

[Olumide]

On Apr 4, 5:49 pm, Bill Greene <w.h.gre...@gmail.com> wrote:
> As Mark Hoemmen pointed out, DSYTRI is the better choice in this case
> because
> the input matrix is symmetric. If works fine, too, but in that case you have
> to
> remember to call DSYTRF first.
>

> ...

>
> On Mon, Apr 4, 2011 at 12:42 PM, Bill Greene <w.h.gre...@gmail.com> wrote:
> > I called DGETRF with your matrix as input followed by DGETRI and it
> > returns the inverse you show below. Did you call DGETRF before
> > your call to DGETRI?

Thanks for your feedback. I'm having trouble linking to DSYTRI. (I'm
programming in C++ on Windows.) I've examined the symbols in
LAPACKd.lib as follows

nm LAPACKd.lib | grep DSYTRI
00000000 T _DSYTRI

I've also added the declaration, extern "C" void _DSYTRI( ...) to the
source file, however I'm getting the error message: unresolved
external symbol __DSYTRI at the linking stage.

Why is the linker appending an extra underscore to _DSYTRI?

[Evgenii Rudnyi]

On 05.04.2011 02:43 Olumide said the following:

The right declaration would be extern "C" void DSYTRI( ...). The first
underscore is added automatically. This is the way the linker works.

[Olumide]

On Apr 5, 6:12 pm, Evgenii Rudnyi <use...@rudnyi.ru> wrote:
> On 05.04.2011 02:43 Olumide said the following:

> > Why is the linker appending an extra underscore to _DSYTRI?
>
> The right declaration would be extern "C" void DSYTRI( ...). The first
> underscore is added automatically. This is the way the linker works.

Just as I thought. The problem is that when I use the declaration
extern "C" void DSYTRI( ...), the application crashes. The only
information I get when I debug is

MyApplication.exe!_mkl_blas_p4_dswap() + 0x86 Asm

It would seem that my application is linking against MKL BLAS. This is
very unusual because (i) I've never had this problem (ii) LAPACK and
GotoBLAS libraries appear before any MKL libraries in my project. (MLK
BLAS is there in the first pace because some of the classes in my
project are wrappers for MUMPS, which is why I've got MKL. The point
I'm trying to make is that I cannot remove MKL libraries from my
dependencies.)

[Evgenii Rudnyi]

On 05.04.2011 19:31 Olumide said the following:

First check parameters to the function. It well might be that something
is wrong and then there is a crash within the function. Alternatively
make a small test to call it some simpler setup. Then you can try to
link it independently to different BLAS.

[Olumide]

On Apr 5, 6:31 pm, Olumide <50...@web.de> wrote:
> Just as I thought. The problem is that when I use the declaration
> extern "C" void DSYTRI( ...), the application crashes. The only
> information I get when I debug is
>
> MyApplication.exe!_mkl_blas_p4_dswap()  + 0x86  Asm
>

I think this problem was caused by the fact that I made some non-const
arguments const. Anyway, here's what my implementation (DSYTRF +
DSYTRI) looks like:

// _dim is the dimension of the matrix
// _uplo is a pointer to the character 'L'
// _A is a pointer the array of matrix data

int info;
int* ipiv = new int [ _dim ];
double* work = new double [ _dim ];

DSYTRF( _uplo , &_dim , _A , &_dim , ipiv , work , &_dim , &info );

if( info != 0 )
{
cout << "Step 1 | DSYTRF failed ..." << endl;
return info;
}

int lwork = (int) work[0];
//work = (double*) realloc( work , sizeof(double) * lwork );
double* nwork = new double [ lwork ];

DSYTRI( _uplo , &_dim , _A , &_dim , ipiv , nwork , &info );

return info;

Here's the the lower matrix after DSYTRF
[ 0.000000 ]
[ 529.831737 0.000000 ]
[ 1198.292909 529.831737 0.000000 ]
[ 529.831737 1198.292909 529.831737 0.000000 ]
[ 1.000000 1.000000 1.000000 1.000000 0.000000 ]

[ 10.000000 0.000000 -10.000000 0.000000 0.000000
0.000000 ]

[ 0.000000 10.000000 0.000000 -10.000000 0.000000
0.000000 0.000000 ]

And here's the result after DSYTRI

[ 0.000000 ]
[ -0.000000 0.000000 ]
[ -0.000002 0.000007 0.000000 ]
[ 0.001891 0.000833 529.831737 0.000000 ]
[ -0.000000 -0.001891 1.000000 -0.000000 0.000000 ]
[ 0.000000 0.000000 -0.000368 0.000000 0.000000
-0.000000 ]
[ -0.000000 -0.000000 -0.000002 -10.000000 0.001895
0.000000 -0.000000 ]

It looks a very suspect. Is it correct? ...

And how do I get the inverse from this lower triangular matrix? My
matrix theory isn't very strong.

Thanks.

[Olumide]

Hooray, DGETRF + DGETRI works!!! ... The inverse is:

[ 0.001803 -0.001803 0.001803 -0.001803 0.250000
0.050000 -0.000000 ]
[ -0.001803 0.001803 -0.001803 0.001803 0.250000
-0.000000 0.050000 ]
[ 0.001803 -0.001803 0.001803 -0.001803 0.250000
-0.050000 -0.000000 ]
[ -0.001803 0.001803 -0.001803 0.001803 0.250000
-0.000000 -0.050000 ]
[ 0.250000 0.250000 0.250000 0.250000 -564.489096
0.000000 0.000000 ]
[ 0.050000 0.000000 -0.050000 0.000000 -0.000000
5.991465 0.000000 ]
[ 0.000000 0.050000 0.000000 -0.050000 0.000000
0.000000 5.991465 ]

[Olumide]

Here's related question. Would DSPSV be the correct routine to use in
solving linear systems Ax = b, where A is a positive semidefinite
matrix -- the sort I'm trying to invert, that has zeros along its
diagonal e.g. in the original post.

[Mark Hoemmen]

DSPSV is for packed symmetric storage. DSYSV is for unpacked
symmetric storage -- the usual way that you would store a nonsymmetric
matrix (the routine only reads and writes the upper or lower triangle,
depending on the UPLO argument). It looks like you are using unpacked
symmetric storage.

mfh

[Olumide]

I've just taken another look at my code. I used DGESV in another case
where the matrix A is not symmetrical but had zeros along its
diagonal. My concern is the appropriateness of the subroutines in
cases where the matrix A has zeros along its diagonal.

[Mark Hoemmen]

On Tue, Apr 5, 2011 at 16:51, Olumide <50...@web.de> wrote:
> I've just taken another look at my code. I used DGESV in another case
> where the matrix A is not symmetrical but had zeros along its
> diagonal. My concern is the appropriateness of the subroutines in
> cases where the matrix A has zeros along its diagonal.

_SY__(_) routines are for any symmetric matrix. As long as your
matrix is nonsingular, the solver driver routines should work just
fine.

If your matrices are as small as the example problem you enclosed, you
should just use DGESV and stop worrying ;-)

mfh

[Olumide]

Thanks Mark. The matrix in my original post is just a toy example. The
actual matrices I'm working with are huge, up to order 5000 and always
have zeros along their diagonals, and are not always symmetric. What
concerns me is that DPOTRI and DGETRI had returned errors that
indicated that zeros along the diagonal of the matrix. I just wanted
to be sure that I'm using the correct subroutines.

[Mark Hoemmen]

On Tue, Apr 5, 2011 at 17:41, Olumide <50...@web.de> wrote:
> Thanks Mark. The matrix in my original post is just a toy example. The
> actual matrices I'm working with are huge, up to order 5000 and always
> have zeros along their diagonals, and are not always symmetric. What
> concerns me is that DPOTRI and DGETRI had returned errors that
> indicated that zeros along the diagonal of the matrix. I just wanted
> to be sure that I'm using the correct subroutines.

DPOTRI should return an error, since that routine should only work for
symmetric positive definite ("PO") matrices. "SY" routines are for
symmetric, but not necessarily positive definite matrices.

mfh