Turn on suggestions

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

Showing results for

- Intel Community
- Software Development Tools (Compilers, Debuggers, Profilers & Analyzers)
- Intel® Fortran Compiler
- Inverse matrix with Intel MKL

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

Highlighted
##

Irene_A_

Beginner

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

12-26-2016
05:47 PM

102 Views

Inverse matrix with Intel MKL

Dear All,

I am trying to migrate my code from Compaq Visual Fortran 6.5 to Visual Studio 2015 with Intel Fortran Compiler. I understand that IMSL used to be built in CVF. Now however, IMSL is an additional add-on for Visual Studio. I am having problems trying to change the program for inversing matrix using MKL.

The original code was:

subroutine MatA(ap)

use msimsl

implicit none

include 'eldata.h'

include 'cdata.h'

include 'sdata.h'

include 'iofile.h'

real*8 ap(4,4,numel),iap(4,4),as(4,4),ias(4,4)

integer ii,jj,ik

mct=0

do 15 ik=1,numel

do 15 ii=1,4

do 15 jj=1,4

as(ii,jj)=as(ii,jj)+ap(ii,jj,ik)

15 continue

call dlinrg(4,as,4,ias,4)

do ii=1,4

if(mct.eq.0) then

write(iow,2007)

if(ior.lt.0) write(*,2007)

mct = 50

endif

mct=mct-1

write(iow,2008)(ap(ii,jj,1),jj=1,4)

if(ior.lt.0) write(*,2008)(ap(ii,jj,1),jj=1,4)

enddo

2007 format(' C o e f f i s i e n t o f M a t r i x

1 A p e r E l e m e n t',//

1 ,' comment koef. Matriks A Elemen 1 ')

2008 format(4x,'Yoo',10x,e12.4,',',e12.4,',',e12.4,',',e12.4)

end

instead of use imsl, and call dlinrg(4,as,4,ias,4) should i just change it into:

external DGETRF

external DGETRI

ias=as

n = size(as,1)

call DGETRF(n, n, ias, n, ipiv, info)

call DGETRI(n, ias, n, ipiv, work, n, info)

Would that work? I am new to this and don't understand much about Fortran. Thank you

Best regards

Irene

4 Replies

Highlighted
##

mecej4

Black Belt

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

12-27-2016
08:49 AM

102 Views

When obtaining a numerical solution to a set of simultaneous linear equations, it is almost always a mistake to compute the inverse of the matrix. Likewise, if all you want to do is to solve the equations once, the simplest way is to do:

USE lapack95

...

call gesv( A, b [,ipiv] [,info] )

where A is the matrix, b contains the R.H.S. vector on input and the solution x on output, for the linear system A x = b .

The code fragment that you posted does not show us what you do with the computed inverse. If you simply pre-multiply the vector b by that inverse, the solution that I suggested will suffice.

[Added 29 Dec 2016]

After looking at the code more closely, I am puzzled by the purpose of the call to DLINRG. The returned inverse, in array IAS, is not used in the subroutine at all, unless the array is placed in a common block and used elsewhere. If it turns out that IAS is not used, this call to DLINRG could be just commented out.

Highlighted
##

Irene_A_

Beginner

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

12-27-2016
06:11 PM

102 Views

I've attached another part of the code where I use the inverse of the matrix.

So for example,

MATAX(jj,kj)=MATAX(jj,kj)+iap(jj,ij)*BMATX(ij)

where iap is the inverse of the matrix

I would just have to change iap(jj,ij)*BMATX(ij) as the b in the function you mentioned?

Would that still work? What about looping it like above (jj, ij, ij)? How would I do it?

Thank you

Highlighted
##

mecej4

Black Belt

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

12-27-2016
08:03 PM

102 Views

Thanks for posting the code, which makes it clear that what the subroutine does is to obtain the least squares solution to A.X = B, where A is a 16 X 8 matrix, X is 8 X 5 and B is 16 X 5. The subroutine code is very wasteful and inefficient, and violates many of the standard prescriptions for solving least squares problems:

- It explicitly forms the normal equations, A
^{T}A.X = A^{T}B, which makes the conditioning of the problem much worse - It explicitly calculates the inverse of A
^{T}A, which makes the solution take twice as long. - Summations that occur more than once are executed repeatedly, instead of being evaluated once and the result reused as needed.

It is quite likely that the subroutine was written decades ago, and may have been a reasonable implementation in that period.

The entire subroutine, apart from the heading and declarations, can be replaced by the following lines of code, with some variable names shortened (not tested!). The main part of the solution is performed by calling the Lapack least-squares solver, GELS. The array arithmetic available in Fortran 90 and later allows us to do away with all the DO loops.

USE LAPACK95 ! A is M X N; B is M X NR; M >= N ! M=16; N=8; NR=5 A(:,1) = 1 A(:,2) = xo A(:,3) = yo A(:,4) = xo*xo A(:,5) = xo*yo A(:,6) = yo*yo A(:,7) = xo*xo*yo A(:,8) = xo*yo*yo B(:,1) = mmx B(:,2) = mmy B(:,3) = mmxy B(:,4) = stxz B(:,5) = styz call gels (A, B) AX (:,kj) = B(:N,1) AY (:,kj) = B(:N,2) AXY(:,kj) = B(:N,3) AXZ(:,kj) = B(:N,4) AYZ(:,kj) = B(:N,5)

Highlighted
##

mecej4

Black Belt

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

12-29-2016
03:54 AM

102 Views

As was stated in a private message, the subroutine under discussion is intended to be used with the FEAP finite element package from UC-Berkeley. It appears to be a modification to allow performing least-squares solution of linear equations instead of the usual Gaussian elimination solution.

Here is a replacement for the subroutine, in which vector operations are used and MKL/Lapack-95 is used instead of IMSL. The subroutine compiled without errors when, instead of the include files for FEAP, which are not freely distributed, I used the include files of the downsized version, FEAPPV from http://www.ce.berkeley.edu/projects/feap/feappv/ . I had to add a dummy zero-length file, **upointer.h**. To compile the subroutine, use the /Qmkl flag or, in Visual Studio, select the MKL library in the project properties.

! subroutine least_sq_fit(xo,yo,kj,mmx,mmy,mmxy,stxz,styz,jac,matax &,matay,mataxy,mataxz,matayz) use lapack95, only : gels implicit none include 'eldata.h' include 'cdata.h' include 'sdata.h' include 'iofile.h' include 'pointer.h' include 'upointer.h' include 'comblk.h' real*8 mmx(16),mmy(16),mmxy(16),stxz(16),styz(16) ! all input real*8 xo(16),yo(16),jac(16) ! xo, yo input real*8 matax(8,numnp),matay(8,numnp),mataxy(8,numnp), & mataxz(8,numnp),matayz(8,numnp) ! all output real*8 amat(16,8), bmat(16,5) integer, intent(in) :: kj amat(:,1) = 1d0 amat(:,2) = xo amat(:,3) = yo amat(:,4) = xo*xo amat(:,5) = xo*yo amat(:,6) = yo*yo amat(:,7) = xo*xo*yo amat(:,8) = xo*yo*yo bmat(:,1) = mmx bmat(:,2) = mmy bmat(:,3) = mmxy bmat(:,4) = stxz bmat(:,5) = styz call gels(amat,bmat) matax (:,kj) = bmat(:8,1) matay (:,kj) = bmat(:8,2) mataxy(:,kj) = bmat(:8,3) mataxz(:,kj) = bmat(:8,4) matayz(:,kj) = bmat(:8,5) return end

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