Turn on suggestions

Auto-suggest helps you quickly narrow down your search results by suggesting possible matches as you type.

Showing results for

- Intel Community
- Software Development Tools (Compilers, Debuggers, Profilers & Analyzers)
- Intel® Fortran Compiler
- IMSL, MKL, matrix operators question

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

Highlighted
##

Brian_Murphy

Beginner

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

02-19-2018
10:01 AM

4 Views

IMSL, MKL, matrix operators question

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.

27 Replies

Highlighted
##

mecej4

Black Belt

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

02-19-2018
02:06 PM

4 Views

Yes, those matrix/vector

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

Highlighted
##

Steve_Lionel

Black Belt

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

02-19-2018
04:37 PM

4 Views

The IMSL run-time licensing

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.

Steve (aka "Doctor Fortran") - https://stevelionel.com/drfortran

Highlighted
##

Brian_Murphy

Beginner

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

02-19-2018
05:01 PM

4 Views

Since it's not my job to

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.

Highlighted
##

Steve_Lionel

Black Belt

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

02-19-2018
06:21 PM

4 Views

Is this for use within your

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.

Steve (aka "Doctor Fortran") - https://stevelionel.com/drfortran

Highlighted
##

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.

Brian_Murphy

Beginner

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

02-19-2018
06:50 PM

4 Views

The application will be

Highlighted
##

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

andrew_4619

Valued Contributor II

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

02-20-2018
01:43 AM

4 Views

look at the example at http:/

Highlighted
##

Brian_Murphy

Beginner

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

02-21-2018
05:22 AM

4 Views

Thank you very much Andrew.

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?

Highlighted
##

mecej4

Black Belt

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

02-22-2018
06:17 AM

4 Views

Brian, you may find this

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?

Highlighted
##

Brian_Murphy

Beginner

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

02-22-2018
07:28 AM

4 Views

Thanks, mecej4. That thread

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?

Highlighted
##

Brian_Murphy

Beginner

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

02-22-2018
07:54 AM

4 Views

What I would like to do is

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

Highlighted
##

Steve_Lionel

Black Belt

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

02-22-2018
10:25 AM

4 Views

Yes, that's the right

Yes, that's the right approach.

Steve (aka "Doctor Fortran") - https://stevelionel.com/drfortran

Highlighted
##

Brian_Murphy

Beginner

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

02-22-2018
10:43 AM

4 Views

Haven't the Intel MKL folks

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?

Highlighted
##

Steve_Lionel

Black Belt

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

02-22-2018
01:35 PM

4 Views

With MATMUL, the compiler

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.

Steve (aka "Doctor Fortran") - https://stevelionel.com/drfortran

Highlighted
##

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.

Brian_Murphy

Beginner

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

02-24-2018
05:24 AM

4 Views

I finally got a reply in the

Highlighted
##

Brian_Murphy

Beginner

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

02-24-2018
06:15 AM

4 Views

The following code gives

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?

Highlighted
##

mecej4

Black Belt

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

02-24-2018
07:37 AM

4 Views

Note that GEMM is a generic

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.

Highlighted
##

andrew_4619

Valued Contributor II

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

02-24-2018
07:40 AM

4 Views

GEMM is in a module you need

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

Highlighted
##

Steve_Lionel

Black Belt

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

02-24-2018
08:23 AM

4 Views

True, but the interface would

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

Steve (aka "Doctor Fortran") - https://stevelionel.com/drfortran

Highlighted
##

Brian_Murphy

Beginner

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

02-24-2018
08:25 AM

4 Views

Below is what I've got. Did

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

For more complete information about compiler optimizations, see our Optimization Notice.