topic Thanks for posting the code, in IntelĀ® Fortran Compiler
https://community.intel.com/t5/Intel-Fortran-Compiler/Inverse-matrix-with-Intel-MKL/m-p/1122859#M132112
<P>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:</P>
<UL>
<LI>It explicitly forms the normal equations, A<SUP>T</SUP>A.X = A<SUP>T</SUP>B, which makes the conditioning of the problem much worse</LI>
<LI>It explicitly calculates the inverse of A<SUP>T</SUP>A, which makes the solution take twice as long. </LI>
<LI>Summations that occur more than once are executed repeatedly, instead of being evaluated once and the result reused as needed.</LI>
</UL>
<P>It is quite likely that the subroutine was written decades ago, and may have been a reasonable implementation in that period.</P>
<P>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.</P>
<PRE class="brush:fortran;">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)</PRE>
<P> </P>Wed, 28 Dec 2016 04:03:00 GMTmecej42016-12-28T04:03:00ZInverse matrix with Intel MKL
https://community.intel.com/t5/Intel-Fortran-Compiler/Inverse-matrix-with-Intel-MKL/m-p/1122856#M132109
<P>Dear All,</P>
<P>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.</P>
<P>The original code was:</P>
<P> subroutine MatA(ap)<BR />
use msimsl<BR />
implicit none<BR />
include 'eldata.h'<BR />
include 'cdata.h'<BR />
include 'sdata.h'<BR />
include 'iofile.h'<BR />
<BR />
real*8 ap(4,4,numel),iap(4,4),as(4,4),ias(4,4)<BR />
integer ii,jj,ik<BR />
mct=0<BR />
do 15 ik=1,numel<BR />
do 15 ii=1,4<BR />
do 15 jj=1,4<BR />
as(ii,jj)=as(ii,jj)+ap(ii,jj,ik)<BR />
15 continue<BR />
call dlinrg(4,as,4,ias,4)<BR />
do ii=1,4<BR />
if(mct.eq.0) then<BR />
write(iow,2007)<BR />
if(ior.lt.0) write(*,2007)<BR />
mct = 50<BR />
endif<BR />
mct=mct-1<BR />
write(iow,2008)(ap(ii,jj,1),jj=1,4)<BR />
if(ior.lt.0) write(*,2008)(ap(ii,jj,1),jj=1,4)<BR />
enddo<BR />
2007 format(' C o e f f i s i e n t o f M a t r i x <BR />
1 A p e r E l e m e n t',//<BR />
1 ,' comment koef. Matriks A Elemen 1 ')<BR />
2008 format(4x,'Yoo',10x,e12.4,',',e12.4,',',e12.4,',',e12.4)</P>
<P> end</P>
<P>instead of use imsl, and<SPAN style="font-size: 13.008px;"> call dlinrg(4,as,4,ias,4) should i just change it into:</SPAN></P>
<P style="margin-left:36.0pt;"> external DGETRF<BR />
<SPAN style="font-size: 1em;">external DGETRI</SPAN></P>
<P style="margin-left:36.0pt;"> ias=as</P>
<P style="margin-left:36.0pt;"> n = size(as,1)</P>
<P style="margin-left:36.0pt;">call DGETRF(n, n, ias, n, ipiv, info)<BR />
<SPAN style="font-size: 1em;"> call DGETRI(n, ias, n, ipiv, work, n, info)</SPAN></P>
<P style="margin-left:36.0pt;">Would that work? I am new to this and don't understand much about Fortran. Thank you</P>
<P style="margin-left:36.0pt;">Best regards</P>
<P style="margin-left:36.0pt;">Irene </P>Tue, 27 Dec 2016 01:47:31 GMThttps://community.intel.com/t5/Intel-Fortran-Compiler/Inverse-matrix-with-Intel-MKL/m-p/1122856#M132109Irene_A_2016-12-27T01:47:31ZWhen obtaining a numerical
https://community.intel.com/t5/Intel-Fortran-Compiler/Inverse-matrix-with-Intel-MKL/m-p/1122857#M132110
<P>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:</P>
<P> USE lapack95</P>
<P> ...</P>
<P> call gesv( A, b [,ipiv] [,info] )</P>
<P>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 .</P>
<P>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.</P>
<P>[Added 29 Dec 2016] </P>
<P>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.</P>Tue, 27 Dec 2016 16:49:00 GMThttps://community.intel.com/t5/Intel-Fortran-Compiler/Inverse-matrix-with-Intel-MKL/m-p/1122857#M132110mecej42016-12-27T16:49:00ZI've attached another part of
https://community.intel.com/t5/Intel-Fortran-Compiler/Inverse-matrix-with-Intel-MKL/m-p/1122858#M132111
<P>I've attached another part of the code where I use the inverse of the matrix.</P>
<P>So for example,</P>
<P>MATAX(jj,kj)=MATAX(jj,kj)+iap(jj,ij)*BMATX(ij)</P>
<P>where iap is the inverse of the matrix</P>
<P>I would just have to change <SPAN style="font-size: 13.008px;">iap(jj,ij)*BMATX(ij) as the b in the function you mentioned?</SPAN></P>
<P><SPAN style="font-size: 13.008px;">Would that still work? What about looping it like above (jj, ij, ij)? How would I do it?</SPAN></P>
<P><SPAN style="font-size: 13.008px;">Thank you</SPAN></P>Wed, 28 Dec 2016 02:11:49 GMThttps://community.intel.com/t5/Intel-Fortran-Compiler/Inverse-matrix-with-Intel-MKL/m-p/1122858#M132111Irene_A_2016-12-28T02:11:49ZThanks for posting the code,
https://community.intel.com/t5/Intel-Fortran-Compiler/Inverse-matrix-with-Intel-MKL/m-p/1122859#M132112
<P>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:</P>
<UL>
<LI>It explicitly forms the normal equations, A<SUP>T</SUP>A.X = A<SUP>T</SUP>B, which makes the conditioning of the problem much worse</LI>
<LI>It explicitly calculates the inverse of A<SUP>T</SUP>A, which makes the solution take twice as long. </LI>
<LI>Summations that occur more than once are executed repeatedly, instead of being evaluated once and the result reused as needed.</LI>
</UL>
<P>It is quite likely that the subroutine was written decades ago, and may have been a reasonable implementation in that period.</P>
<P>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.</P>
<PRE class="brush:fortran;">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)</PRE>
<P> </P>Wed, 28 Dec 2016 04:03:00 GMThttps://community.intel.com/t5/Intel-Fortran-Compiler/Inverse-matrix-with-Intel-MKL/m-p/1122859#M132112mecej42016-12-28T04:03:00ZAs was stated in a private
https://community.intel.com/t5/Intel-Fortran-Compiler/Inverse-matrix-with-Intel-MKL/m-p/1122860#M132113
<P>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.</P>
<P>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, <STRONG>upointer.h</STRONG>. To compile the subroutine, use the /Qmkl flag or, in Visual Studio, select the MKL library in the project properties.</P>
<PRE class="brush:fortran;">!
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</PRE>
<P> </P>Thu, 29 Dec 2016 11:54:00 GMThttps://community.intel.com/t5/Intel-Fortran-Compiler/Inverse-matrix-with-Intel-MKL/m-p/1122860#M132113mecej42016-12-29T11:54:00Z