Intel® Fortran Compiler
Build applications that can scale for the future with optimized code designed for Intel® Xeon® and compatible processors.

Inverse matrix with Intel MKL

Irene_A_
Beginner
1,736 Views

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 

0 Kudos
4 Replies
mecej4
Honored Contributor III
1,736 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.

0 Kudos
Irene_A_
Beginner
1,736 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

0 Kudos
mecej4
Honored Contributor III
1,736 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, ATA.X = ATB, which makes the conditioning of the problem much worse
  • It explicitly calculates the inverse of ATA, 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)

 

0 Kudos
mecej4
Honored Contributor III
1,736 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

 

0 Kudos
Reply