- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

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

Link Copied

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

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.

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

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.

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

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.

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

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)

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

*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.

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

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

- Subscribe to RSS Feed
- Mark Topic as New
- Mark Topic as Read
- Float this Topic for Current User
- Bookmark
- Subscribe
- Printer Friendly Page