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:
do 15 ik=1,numel
do 15 ii=1,4
do 15 jj=1,4
mct = 50
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 ')
instead of use imsl, and call dlinrg(4,as,4,ias,4) should i just change it into:
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
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:
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.
I've attached another part of the code where I use the inverse of the matrix.
So for example,
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?
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)
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