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
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)
For more complete information about compiler optimizations, see our Optimization Notice.