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.
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
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
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.
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.
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.
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
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.
@FortranFan - what is the motivation for replacing DP with WP? Is it to keep DP free for something else?
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.
For more complete information about compiler optimizations, see our Optimization Notice.