Community
cancel
Showing results for 
Search instead for 
Did you mean: 
Highlighted
Beginner
17 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
1 View

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

 

 

0 Kudos
Highlighted
Beginner
1 View

Thanks, mecej4.  I was trying to go by what I saw in mklman_f_11.3_0.pdf.  Doing the calls they way you show got it compiling with USE BLAS95.  However, my code still won't link.  I will try one of the MKL\examples\BLAS95 and follow the instructions at:

Intel® Math Kernel Library 11.0.5 User's Guide

Creating, Configuring, and Running the Intel Visual Fortran Project

This section demonstrates how to create an Intel Visual Fortran project running an Intel MKL example in Microsoft Visual Studio 2008

I've already tried applying those instructions to my project, and linking still fails.  I will give it a go with the blas95 example and hopefully find what did wrong.

0 Kudos
Highlighted
Black Belt
1 View

You have to add mkl_blas95.lib (and mkl_lapack95.lib, if relevant) explicitly; the /Qmkl option does not pull these in. Please consult the MKL Link Line Advisor at https://software.intel.com/en-us/articles/intel-mkl-link-line-advisor .

Please report the linking command used and the linker's error messages.Or, if you are using Visual Studio to build, attach the build log of a failed build attempt. Simply saying "won't link" is too vague to be of use in diagnosing the problem.

0 Kudos
Highlighted
Beginner
1 View

Many thanks again, mecej4.  I added C:\Program Files (x86)\Intel\Composer XE 2013\mkl\lib\ia32\mkl_blas95.lib to the list of source file for the visual studio project, and it now links successfully.  I have not actually run the code yet.  I will do that as soon as I can in the next few days.

0 Kudos
Highlighted
Valued Contributor III
1 View

mecej4 wrote:

Your grasp of the F9X vs. F77 interfaces to MKL is not firm enough. ..

@Brian Murphy:

Depending on one's computing circumstances, the following might appear a real (pun intended!) nit; or it can turn out helpful from a portability aspect.  See this blog:

https://software.intel.com/en-us/blogs/2017/03/27/doctor-fortran-in-it-takes-all-kinds

Considering your code along with @mecej4's solution for you and also how BLAS95 module is setup, you may want to consider using consistent KINDs of your floating-point (REAL) variables and move away from non-standard REAL*8 declarations for sure but also standard but somewhat ambiguous DOUBLE PRECISION:

module mylinear_operators

   use F95_PRECISION, only: WP => DP
   use BLAS95, only : GEMM, GEMV

   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(WP), dimension(:,:), intent(in) :: a,b
      real(WP), dimension(size(a,1),size(b,2)) :: myMultAB
      call gemm(a,b, myMultAB)
   end function myMultAB

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

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

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

end module mylinear_operators

program xgemm

   use F95_PRECISION, only: WP => DP
   use mylinear_operators, only : operator(.x.)

   implicit none

   real(WP) :: a(2,3) = [11.0_wp, 21.0_wp, 12.0_wp, 22.0_wp, 13.0_wp, 23.0_wp]
   real(WP) :: b(3,2) = [11.0_wp, 21.0_wp, 12.0_wp, 22.0_wp, 13.0_wp, 23.0_wp]
   real(WP) :: c(2,2)

   c = a.x.b

   print '(1x,2ES14.5)',c(1,:),c(2,:)

end program

If you plan to work more with Fortran, consider the references in this Dr Fortran blog too.

https://software.intel.com/en-us/blogs/2013/12/30/doctor-fortran-in-its-a-modern-fortran-world

0 Kudos
Highlighted
Beginner
1 View

I appreciate the tips on modernizing the code.  Portability is not much of an issue at this point, but anything that makes it easier to maintain going forward is good.  I have run some tests and compared output to the old code, and so far it is looking very good.  This interface operators thing is really nice.

Sorry for not being more clear about the linking errors.  I was still getting "unresolved external [GEMM]" which I mentioned in post #20.

0 Kudos
Highlighted
Beginner
1 View

@FortranFan - what is the motivation for replacing DP with WP?  Is it to keep DP free for something else?

0 Kudos
Highlighted
Valued Contributor III
1 View

Brian Murphy wrote:

@FortranFan - what is the motivation for replacing DP with WP?  Is it to keep DP free for something else?

A couple of things: first look in blas.f90 on your system as part of your Intel Fortran installation and see how the generic interface to GEMM is setup.  Secondly you will realize WP, as mnemonic for working precision, is also reminding you an arbitrary precision is in effect for your computations rather than the one with double (or single) precision whose definition on some systems (or certain circumstances) may not match what you presume it to be and thus might prove misleading.

0 Kudos