Auto-suggest helps you quickly narrow down your search results by suggesting possible matches as you type.

Showing results for

- Intel Community
- Software Development SDKs and Libraries
- Intel® oneAPI Math Kernel Library & Intel® Math Kernel Library
- Matrix inversion and MATMUL

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

ssiddarth

Beginner

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

08-20-2010
02:04 PM

86 Views

Matrix inversion and MATMUL

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

6 Replies

ArturGuzik

Valued Contributor I

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

08-20-2010
08:42 PM

86 Views

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

Beginner

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

08-22-2010
02:16 PM

86 Views

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

Black Belt

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

08-22-2010
07:55 PM

86 Views

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

Beginner

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

08-23-2010
10:03 AM

86 Views

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

Black Belt

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

08-23-2010
11:01 AM

86 Views

It returns the upper triangular sub-matrix.

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

Beginner

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

08-23-2010
11:14 AM

86 Views

Thanks!

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

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

For more complete information about compiler optimizations, see our Optimization Notice.