Matrix inversion and MATMUL

[ssiddarth]

In trying to debug a program and figure out what was happening I wrote a simple program that inverts a matrix and multiplies it by its inverse. The problem is that I don't get back the identity matrix. Thus in the program below T1 is not the identity matrix.

Another issue I discovered is that I get a compilation error when I try to use the general matrix inversion routine "call getri( truelambda, ipiv )". The error is C:\\Allwork\\PROJECTS\\SML\\FORTRAN\\PROGRAMS\\Test3.f90(19): error #6285: There is no matching specific subroutine for this generic subroutine call. [GETRI]

Any help would be appreciated.

Siddarth

PROGRAM TRUE_PARAMS

USE MKL95_LAPACK

USE MKL95_PRECISION

IMPLICIT NONE

REAL*8 :: LBD12(4,4),T(4,4),TRUELAMBDA(4,4),T1(4,4),ipiv(1,4)

INTEGER :: I, J

DATA ((LBD12(I, J), J = 1,4), I = 1,4) / 5, 5, -1, 1, 0, 4, -2, 0, 0, 0, 3, 2, 0, 0,0, 1 /

LBD12=LBD12/2.0D0

TRUELAMBDA=MATMUL(TRANSPOSE(LBD12),LBD12)

T=TRUELAMBDA

CALL potrf(TRUELAMBDA)

CALL potri(TRUELAMBDA)

T1=MATMUL(T,TRUELAMBDA)

!CALL getrf(TRUELAMBDA)

!call getri( truelambda, ipiv )

!T1=MATMUL(T,TRUELAMBDA)

END PROGRAM TRUE_PARAMS

[ArturGuzik]

Hi,

I will try to address the second problem you described (as it's easier to spot), namely, that there is no matching specific subroutine for generic GETRI. Well, your ipiv definition is incorrect. The routine expects array (vector) of length = max(1,n), so in your case you shoudl declare it as ipiv(4). As the expected number of dimensions is 1, and you gave 2, so.....

A.

[ssiddarth]

Thanks for pointing that out. Now GETRI works fine.

However, even when I use the getrf and getri sequence to invert the original matrix the product of the matrix and its inverse is not the identity matrix.

S.

[mecej4]

The call to ?POTRF computes the upper or lower triangular Cholesky factor of the input matrix. The subsequent call to ?POTRI returns only the corresponding triangular part of the inverse.

The MKL documentation is not quite clear about this point.

Furthermore, the documentation of the F95 interface POTRI should state that the argument A contains at subroutine entry the triangular factor resulting from a prior call to POTRF, rather than the original matrix A.

[ssiddarth]

The documentation states that the matrix a is overwritten by the n-by-n matrix inv(A), which I took to mean the whole matrix.

I understood your comment to mean that potri returns the upper triangular matrix of the inverse matrix.
I therefore set the lower diagonal terms to zero and then calculated Utranspose*U (see commands below). However, this does not seem to give me the right inverse. When I multiply the original matrix, by this, I still don't get the identity matrix. In the code below T2 does not contain the identity matrix

How does one get the inverse matrix from what is outputby POTRI?

Thanks!

CALL potrf(TRUELAMBDA)

write(*,*) 'factored'

write(*,10) ((TRUELAMBDA(I,J),J=1,4),I=1,4)

CALL potri(TRUELAMBDA)

write(*,*) 'inverse'

write(*,10) ((TRUELAMBDA(I,J),J=1,4),I=1,4)

forall (j=1:4)

forall (i=4:j+1:-1)

TRUELAMBDA(i,j)=0.0d0

end forall

end forall

write(*,10)((TRUELAMBDA(I,J),J=1,4),I=1,4)

T1=MATMUL(TRANSPOSE(TRUELAMBDA),TRUELAMBDA)

t2=matmul(t,t1)

[mecej4]

I understood your comment to mean that potri returns the upper triangular matrix of the inverse matrix

It returns the upper triangular sub-matrix.

How does one get the inverse matrix from what is outputby POTRI?

By symmetry;

(i) if you called POTRI(A,UPLO) with UPLO='U' or with no UPLO argument, after return from POTRI set a(i,j) = a(j,i) for i > j; in other words, obtain the other triangle by reflection in the diagonal.

(ii) if you used UPLO='L', after return from POTRI set a(i,j) = a(j,i) for i < j

This should become obvious upon a little reflection (pun intended): Cholesky factorization is possible only for positive definite matrices. On the other hand a general matrix multiplication routine knows nothing about symmetric matrices, triangular storage, etc. What it expects are full matrices; therefore, you have to flesh out your symmetric matrices if they are presently in any other form than full-matrix form, before calling MATMULT. Alternatively, if you had a special matrix-multiply routine for symmetric matrices, you could have used that routine.

[ssiddarth]

Thanks!

Sorry, you're right, it wasn't a great question.