Community
cancel
Showing results for 
Search instead for 
Did you mean: 
Highlighted
Beginner
4 Views

IMSL, MKL, matrix operators question

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.

0 Kudos
27 Replies
Highlighted
Black Belt
4 Views

Yes, those matrix/vector

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).

0 Kudos
Highlighted
Black Belt
4 Views

The IMSL run-time licensing

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.

Steve (aka "Doctor Fortran") - https://stevelionel.com/drfortran
0 Kudos
Highlighted
Beginner
4 Views

Since it's not my job to

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.

0 Kudos
Highlighted
Black Belt
4 Views

Is this for use within your

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.

Steve (aka "Doctor Fortran") - https://stevelionel.com/drfortran
0 Kudos
Highlighted
Beginner
4 Views

The application will be

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.

0 Kudos
Highlighted
Valued Contributor II
4 Views

look at the example at http:/

0 Kudos
Highlighted
Beginner
4 Views

Thank you very much Andrew. 

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?

0 Kudos
Highlighted
Black Belt
4 Views

Brian, you may find this

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?

0 Kudos
Highlighted
Beginner
4 Views

Thanks, mecej4.  That thread

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?

0 Kudos
Highlighted
Beginner
4 Views

What I would like to do is

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

 

0 Kudos
Highlighted
Black Belt
4 Views

Yes, that's the right

Yes, that's the right approach.

Steve (aka "Doctor Fortran") - https://stevelionel.com/drfortran
0 Kudos
Highlighted
Beginner
4 Views

Haven't the Intel MKL folks

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?

0 Kudos
Highlighted
Black Belt
4 Views

With MATMUL, the compiler

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.

Steve (aka "Doctor Fortran") - https://stevelionel.com/drfortran
0 Kudos
Highlighted
Beginner
4 Views

I finally got a reply in the

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.

0 Kudos
Highlighted
Beginner
4 Views

The following code gives

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?

 

0 Kudos
Highlighted
Black Belt
4 Views

Note that GEMM is a generic

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.

0 Kudos
Highlighted
Valued Contributor II
4 Views

GEMM is in a module you need

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

0 Kudos
Highlighted
Black Belt
4 Views

True, but the interface would

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.)

Steve (aka "Doctor Fortran") - https://stevelionel.com/drfortran
0 Kudos
Highlighted
Beginner
4 Views

Below is what I've got.  Did

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

 

0 Kudos