- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

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.

Link Copied

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

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).

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

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.

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

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.

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

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.

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

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.

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

look at the example at http://www.mathcs.emory.edu/~cheung/Courses/561/Syllabus/6-Fortran/operators.html

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

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?

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

Brian, you may find this older thread pertinent:

https://software.intel.com/en-us/forums/intel-visual-fortran-compiler-for-windows/topic/500652

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?

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

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?

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

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

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

Yes, that's the right approach.

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

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?

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

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.

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

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.

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

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?

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

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.

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

GEMM is in a module you need the appropriate USE statement that will supply the interface.

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

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.)

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

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

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

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

- Subscribe to RSS Feed
- Mark Topic as New
- Mark Topic as Read
- Float this Topic for Current User
- Bookmark
- Subscribe
- Printer Friendly Page