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.
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).
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.
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.
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.
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.
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?
Brian, you may find this older thread pertinent:
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?
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?
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
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?
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.
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.
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?
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.
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.)
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