Intel® Fortran Compiler
Build applications that can scale for the future with optimized code designed for Intel® Xeon® and compatible processors.

IMSL, MKL, matrix operators question

Brian_Murphy
New Contributor II
872 Views

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
Brian_Murphy
New Contributor II
103 Views

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
mecej4
Honored Contributor III
103 Views

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
Brian_Murphy
New Contributor II
103 Views

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
FortranFan
Honored Contributor II
103 Views

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
Brian_Murphy
New Contributor II
103 Views

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
Brian_Murphy
New Contributor II
103 Views

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

0 Kudos
FortranFan
Honored Contributor II
103 Views

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
Reply