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.

[Brian_Murphy]

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.

[mecej4]

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.

[Brian_Murphy]

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.

[FortranFan]

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

[Brian_Murphy]

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.

[Brian_Murphy]

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

[FortranFan]

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.