IMSL, MKL, matrix operators question

[Brian_Murphy]

A code I'm working on uses matrix operations like  A = B .x. C and A = B .ix. C and A = B .tx. C

My understanding is that this is enabled by IMSL.  If that's right, does MKL also enable this, or is there another way?  My sponsor is telling me that distributing our code with IMSL requires buying an IMSL license for every user.  So he's asking me to eliminate IMSL.

[mecej4]

Yes, those matrix/vector operators are available in IMSL but not in MKL.

If you have only a few instances of those, you can remove those expressions and write in the corresponding calls to MKL/Lapack routines. If you follow that route, you may find it easier to use the Lapack95 interfaces (fewer arguments, no need to allocate and pass work arrays).

[Steve_Lionel]

The IMSL run-time licensing issue is a bit complicated. It is explained well in https://software.intel.com/en-us/articles/imsl-fortran-numerical-library-60-end-user-license-agreement/ . This applies only for IMSL purchased from Intel for use with Intel Visual Fortran for Windows.

[Brian_Murphy]

Since it's not my job to worry about licenses, I passed the buck to my sponsor.  He had someone look into it, and all I know is they reported back to me that an IMSL license needs to be bought for each end user, and this was not practical.  So I have to scope out what it will take to ditch IMSL.

A quick survey of the code indicates roughly 350 instances where these operators are used.  How hard is it to create these operators?  In nearly all cases the matrices are in full storage form, with maybe as few as three in type(d_hbc_sparse) form.

[Steve_Lionel]

Is this for use within your organization or for others? If for internal use, you don't need a run-time license if IMSL is bought from Intel. The "someone" may have looked instead at the license for IMSL bought from RogueWave - that DOES require a separate deployment license for each user.

[Brian_Murphy]

The application will be distributed within the university, and to members of their industrial research consortium.  There about three dozen member companies, and roughly 100 non-university users around the world.

[andrew_4619]

look at the example at http://www.mathcs.emory.edu/~cheung/Courses/561/Syllabus/6-Fortran/operators.html

[Brian_Murphy]

Thank you very much Andrew.  With Google I found similar web pages, but that one is much better than any I found.

Defining matrix operators seems like something many fortran users would want.  I'll bet others have already done this with MKL routines.  Does Intel have an include file for MKL that does this?

[mecej4]

Brian, you may find this older thread pertinent:

     https://software.intel.com/en-us/forums/intel-visual-fortran-compiler-for-windows/topic/500652

I have never used the IMSL linear operators modules myself; one of the reasons is perhaps the lack of an operator to read, write and copy matrices (from the standard 2-D array representation to the IMSL matrix type).

In the IMSL manuals, you see that they have one function call to load a single element, so a 1000 X 1000 matrix would entail a million function calls just to read the matrix in. How was this issue handled in your CVF code?

[Brian_Murphy]

Thanks, mecej4.  That thread is interesting.  As I read through it, I expected the discussion to eventually touch on 2D matrix operations, but it didn't.  Unless I'm missing something, 2D matrix operations would be a great place to use that feature.  e.g.  A = PHi .tx. K .x. PHI

The IMSL operators evidently work with regular 2 dimensional arrays.  The code I'm working on does just that in the vast majority of instances where it's used.  The operators also work with IMSL's built-in Harwell sparse matrices defined type.  IMSL, in addition to creating operators like .tx. also overloads at least the = operator so one can use H=A to load a regular matrix into a Harwell.

There was mention of an include file in gcc (gcc/testsuite/gfortran.dg/userdef_operator_1.f90).  I might try to get that and see what's inside.  I am an old dog trying to learn new tricks.  From the little bit that I've learned, I'm surprised that creating .x. for matrix and vector operations is not a widely used feature of new fortran.  Hasn't this capability been around at least since F90?  If it's not widely used, is there some obvious reason for that?

[Brian_Murphy]

What I would like to do is replace IMSL's USE LINEAR_OPERATORS with USE MYOWNOPERATORS so every place where the code presently has K_B8 = BSTAR .tx. K_B .x. BSTAR will work without needing to change anything.  The matrices are all regular 2D real(8) arrays.

Here is my feeble-minded start to MYOWNOPERATORS.  I of course have not tried this, and I expect it needs a lot more logic to work.  I merely edited what was posted earlier for .cross.

module MYOWNOPERATORS
implicit none
 
interface operator(.x.)
module procedure myMult
end interface
 
contains
 
pure function myMult(a, b)
real, dimension(:,:), intent(in) :: a,b
real, dimension(size(a,1),size(b,2)) :: myMult
 
myMult = MATMUL(a,b)
end function myMult
 
end module MYOWNOPERATORS

[Steve_Lionel]

Yes, that's the right approach.

[Brian_Murphy]

Haven't the Intel MKL folks already done this with calls to ?gemm and so forth?  They would be far better at doing this than me.  I posted a question about this in MKL forum.  If I find out anything good, I'll report that here.

By the way, how would gemm compare to MATMUL for matrix multiplication?

[Steve_Lionel]

With MATMUL, the compiler might choose to do it inline or call a specialized MKL routine. If your arrays are of any decent size, I'd suggest calling ?gemm directly - MKL has "optimized the snot" out of it.

[Brian_Murphy]

I finally got a reply in the MKL forum.  Intel does not have a module of interface operators for matrix math.  I must be missing something here.  Because it seems to me that half the fortran programmers on the planet would want like this.

[Brian_Murphy]

The following code gives a compile error.

module mylinear_operators
 implicit none
 interface operator(.x.)
  module procedure myMult
 end interface
 
 contains
  pure function myMult(a, b)
   real, dimension(:,:), intent(in) :: a,b
   real, dimension(size(a,1),size(b,2)) :: myMult
 
   call gemm(a,b, myMult)
  end function myMult
 
end module mylinear_operators

The error is: Error 1  error #7137: Any procedure referenced in a PURE procedure, including one referenced via a defined operation or assignment, must have an explicit interface and be declared PURE.   [GEMM]

Can someone tell what the explicit pure interface for GEMM would look like?

[mecej4]

Note that GEMM is a generic interface name that stands for several subroutines (SGEMM_F95, DGEMM_F95, etc.).

You can find the complete details and the interfaces in the file <installation_directory>\windows\mkl\include\blas.f90.

[andrew_4619]

GEMM is in a module you need the appropriate USE statement that will supply the interface. 

[Steve_Lionel]

True, but the interface would need to declare the specific routines as PURE to allow a call from a pure function. If the MKL module doesn't do this, which is likely, you could make your own interface that calls the procedures pure, assuming of course you are certain they qualify (I think they do.)

[Brian_Murphy]

Below is what I've got.  Did I add the most appropriate include file?  With INCLUDE 'MKL.FI' it compiles fine.  Or should I have USE BLAS95, and does that require adding blas.f90 to the project?

Anyhow, as it is now, linking fails with GEMM as an unresolved external.  In visual studio do I need to tell it to include MLK_BLAS95.LIB?  I tried putting this in Linker/Input/Additional Dependencies, but that didn't work.  I have /Qmkl:sequential turned on in the compiling section.

module mylinear_operators
 implicit none
 INCLUDE 'mkl.fi'
 interface operator(.x.)
  module procedure myMultAB, myMultAv
 end interface
 
 interface operator(.tx.)
  module procedure myTMult, myTMultAv
 end interface
 
 contains
  function myMultAB(a, b)
   real*8, dimension(:,:), intent(in) :: a,b
   real*8, dimension(size(a,1),size(b,2)) :: myMultAB
   call gemm(a,b, myMultAB)
  end function myMultAB
 
  function myMultAv(a, b)
   real*8, dimension(:,:), intent(in) :: a
   real*8, dimension(:),   intent(in) :: b
   real*8, dimension(size(a,1)) :: myMultAv
   call gemm(a,b, myMultAv)
  end function myMultAv
 
  function myTMult(a, b)
   real*8, dimension(:,:), intent(in) :: a,b
   real*8, dimension(size(a,2),size(b,2)) :: myTMult
   call gemm('T','N',a,b, myTMult)
  end function myTMult
 
  function myTMultAv(a, b)
   real*8, dimension(:,:), intent(in) :: a
   real*8, dimension(:),   intent(in) :: b
   real*8, dimension(size(a,2)) :: myTMultAv
   call gemm('T','N',a,b, myTMultAv)
  end function myTMultAv
 
end module mylinear_operators

[mecej4]

Your grasp of the F9X vs. F77 interfaces to MKL is not firm enough. Here are some steps to help get started, and a working example.

1. Decide whether you will use BLAS95 and Lapack95, or the older F77 interfaces to the BLAS and Lapack routines in MKL. I recommend that you use the F95 interfaces, since they are shorter and easier to use.

2. Depending on which version of Intel Fortran and MKL you have, you may have to build the modules f95_precision.mod, blas95.mod, lapack95.mod and libraries mkl_blas95.lib, etc. (there will be one set for targeting IA32, another for LP64 and a third set for ILP64).

3. Make sure that all the F95 modules, libraries are installed in the MKL include and lib directories.

4. Build and run some of the BLAS95 and Lapack95 examples that come with MKL.

After you have finished Steps 1-4 above, you are set to build and run a linear operator example.

Here is a corrected  version of your mylinear_operators module source:

module mylinear_operators
 USE BLAS95
 implicit none
 interface operator(.x.)
  module procedure myMultAB, myMultAv
 end interface

 interface operator(.tx.)
  module procedure myTMult, myTMultAv
 end interface

 contains
  function myMultAB(a, b)
   real*8, dimension(:,:), intent(in) :: a,b
   real*8, dimension(size(a,1),size(b,2)) :: myMultAB
   call gemm(a,b, myMultAB)
  end function myMultAB

  function myMultAv(a, b)
   real*8, dimension(:,:), intent(in) :: a
   real*8, dimension(:),   intent(in) :: b
   real*8, dimension(size(a,1)) :: myMultAv
   call gemv(a,b, myMultAv)
  end function myMultAv

  function myTMult(a, b)
   real*8, dimension(:,:), intent(in) :: a,b
   real*8, dimension(size(a,2),size(b,2)) :: myTMult
   call gemm(a,b, myTMult, transa='T')
  end function myTMult

  function myTMultAv(a, b)
   real*8, dimension(:,:), intent(in) :: a
   real*8, dimension(:),   intent(in) :: b
   real*8, dimension(size(a,2)) :: myTMultAv
   call gemv(a,b, myTMultAv,trans='T')
  end function myTMultAv

end module mylinear_operators

Compile this source file so that the module file gets generated. We are now ready for a test program:

program xgemm
use mylinear_operators
implicit none
double precision :: a(2,3) = [11d0, 21d0, 12d0, 22d0, 13d0, 23d0]
double precision :: b(3,2) = [11d0, 21d0, 12d0, 22d0, 13d0, 23d0]
double precision :: c(2,2)
!
c = a.x.b
print '(1x,2ES14.5)',c(1,:),c(2,:)
end program

Compile and run the example:

s:\mkl> ifort /Qmkl /Od xop.f90 linop.obj mkl_blas95.lib
s:\mkl> xop
    5.29000E+02   9.69000E+02
    6.97000E+02   1.27700E+03