invert SPD matrix

[v-nijs]

I wrote a test program to learn the basics of MKL like inverting a 2x2 matrix

The routines for a general matrix work fine:

M(1,1) = 3.0
M(1,2) = 2.0
M(2,1) = 2.0
M(2,2) = 3.0

call getrf(M,ipiv)
call getri(M,ipiv)

gives the correct result. However, if i pass the marix M to routines for a symmetric positive-definite matrix I seem to get an incorrect result

call potrf(M)
call potri(M)

gives

0.6 2.0
-0.4 0.6

Is this normal or an error?

Thanks!

Vincent

[bill_woessner]

Assuming you called potrf and potri with uplo == "L", this is the correct output. Because M is assumed the be SPD, the routines only operate on one traingle of M.

If you're using the F95 interface, the uplo and info parameters are optional. Unfortunately, the documentation doesn't say what the default is for uplo. Given your results, if you're calling these routines without passing uplo, it appears the default is "L".

If you like, I can provide you with sample C++ code that demonstrates all this more explicitly.

--Bill

[v-nijs]

Thanks for the reply Bill. I am using the F95 interface. The default is actually 'U'. The output is actually the transpose of the matrix I originally posted.

You comment suggests that i am getting the correct result. The question I have then is how to get the output matrix to be the same as from getrf/getri (i.e., the complete, symmetric, inverse matrix).

I am also curious what the speed benefit of potrf/potri will be once you include this extra operation.

I'd appreciate seeing some the example code you mentioned.

Best,

Vincent

[Anonymous70]

Part of programming with LAPACK is learning to take advantage of any special properties of your matrix. In your case, you're dealing with a symmetric matrix. For the sake of efficiency, you should only dealwith one half of the matrix (though you must allocate storage for the other half, a fact which continues to irk me).

So to answer your question, you shouldn't try to force the output from potrf/potri to look like the output from getrf/getri. Instead, you should try to use LAPACK and BLAS functions that will continue to take advantage of the symmetric nature of your matrix. Here's an example:

#include 
#include "mkl_cblas.h"
#include "mkl_lapack.h"
int main()
{  double A[4], Ainv[4], b[2], x[2];  int info, N(2);
  // A is stored in a flat array in column-major order.
  // Elements 0, 1 and 3 in the array correspond to elements (1, 1), (2, 1) and (2, 2) in the matrix.  // Setting A[2] is unnecessary because we're assuming A is symmetric (i.e. A[2] = A[1]).
  A[0] = 3;
  A[1] = 2;
  A[3] = 3;
  b[0] = 3.1415926535;
  b[1] = 1.618;
  // Ainv <- A
  std::copy(A, A + 4, Ainv);
  // Ainv <- A^-1
  dpotrf("L", &N, Ainv, &N, &info);
  dpotri("L", &N, Ainv, &N, &info);
  // x <- A^-1 * b
  cblas_dsymv(CblasColMajor, CblasLower, 2, 1, Ainv, 2, b, 1, 0, x, 1);  // b <- b - A * x
  cblas_dsymv(CblasColMajor, CblasLower, 2, -1, A, 2, x, 1, 1, b, 1);
  // print ||b - A * x||
  std::cout << cblas_dnrm2(2, b, 1) << '
';
  return EXIT_SUCCESS;
}

As you can see, I didn't even set the upper triangle of A. After the factorization, I used BLAS functions that operate on symmetric matrices (that's the "sy" in "dsymv" stands for symmetric). Of course, this code solves Ax = b by explicitly forming the A^-1 which is a TERRIBLE idea. But it gets the point across.

So that's my advice to you: exploit the structure of your matrix. LAPACK and BLAS have routines for dealing with every sort of matrix you can think of (banded, diagonal, triangular, symmetric, hermetian, etc.). Take adv antage of that.

--Bill