Intel® Fortran Compiler
Build applications that can scale for the future with optimized code designed for Intel® Xeon® and compatible processors.
28749 Discussions

More trouble with MKL routines getrf, getri

dboggs
New Contributor I
1,512 Views

I have a working (or so I thought) program that uses these mkl routines to invert a matrix. The first few lines are

PROGRAM MKLinvert
USE LAPACK95, ONLY: GETRF, GETRI
USE F95_PRECISION, ONLY: WP => SP ! or DP
:

The Fortran libraries are
Debug: Debug QuickWin (/libs:qwin /dbglibs) /Qmkl:parallel
Release: QuickWin (/libs:qwin) /Qmkl:parallel

The linker dependencies are:
mkl_Lapack95.lib

The program builds and runs OK from within the VS environment. However, when I try to run the stand-alone exe files, I get a rte "libiomp5md.dll not found." In the past I have fixed a similar error (report of a missing dll file) by observing that the Fortran libray was "../libs:dll" (by mistake) and changing it to "../libs:static", but that's not the case here.

I DO have the file libiomp5md.dll on my system, in several different locations. Adding this to my link line does not help. Would it help to include the path to one of these? (I suspect that the reported missing dll is a white elephant and the problem is really somewhere else.)

Any ideas?
 

 

0 Kudos
9 Replies
mecej4
Honored Contributor III
1,512 Views

If your build configuration is set to use OpenMP, or you use a library such as MKL, the EXE is going to depend on libiomp5md.dll at run time. Years ago, you could link an EXE that would be completely self-contained, but that is no longer true now. When you installed Ifort, the path to that DLL would have been added to the system path, but that path could have been altered later for one reason or another. You can check your path variable to settle this question.

The MS linker does not link to DLLs directly, as some other linkers do. It uses an import library, instead. The DLL has to be available at run time.

When you use a library such as MKL/Lapack95, three or four paths are significant: (i) the path to included files and module files (usually both the same), (ii) the path to libraries, and (iii) the path to DLLs. These are used during compiling, linking and running, respectively. 

A side comment: except in some very specific applications, one should never need to compute the inverse of a matrix explicitly. There is a major performance penalty associated with forming the inverse. so that should be avoided if possible. If you are solving simultaneous linear equations, for example, you can call GETRF and then GETRS or, even more simply, call the single routine GESV.

0 Kudos
dboggs
New Contributor I
1,512 Views

So, how can I determine which path to use for this application? How do I easily provide the path--explicitly on the link line? It is a very long path. Can I copy the dll to another location that is more convenient to this particular project? e.g., can I put it in the project's home directory and then not have to worry about it?

Incidentally about your side comment: I actually started this subject a year or so ago and explained in this forum why I want the inverse. I am well aware of the overhead in general and the advantage of using an equation solver instead of the coefficient matrix inverse.

I have very little experience with dll's in the IVF environment. It might help if I knew the hierarchy of search locations used. Is it the same as the hierarchy used for include files?

 

0 Kudos
mecej4
Honored Contributor III
1,512 Views

Please see these MSDN articles: https://msdn.microsoft.com/en-us/library/7d83bc18.aspx ; https://docs.microsoft.com/en-us/windows/desktop/Dlls/dynamic-link-library-search-order ;

The MS and Intel compilers use the MS Link program to build an EXE or a DLL. LINK does not read any DLL files. It is only when you run an EXE or invoke a DLL from an EXE that the search path for executables becomes relevant. In fact, if the development and deployment are performed on different machines,  the development machine need not have the DLLs on it at all, and the deployment machine need not have any of the libraries, include and module files.

Yes, you can place copies of the DLL in the same directory as the EXE in your project or somewhere nearby, but in general it is not good to have multiple copies and multiple versions of DLLs lying around in user directories. Bad things can happen if the DLL that is found is not the same as the one that matches the import library that was used at link time.

0 Kudos
gib
New Contributor II
1,512 Views

Once you know which version of the DLL is needed by the program, you can ensure that the path to it appears in the PATH environment variable ahead of the paths to any other versions.

0 Kudos
dboggs
New Contributor I
1,512 Views

Thanks for the valuable info mecej4. I have fixed the problem by placing a copy of libiomp5md.dll in the same directory as the application exe. I am well aware of the risks in doing this, as you describe, but in this case it is the most expedient way. For the record, I found the dll here:

c:\Program Files (x86)\Common Files\Intel\Shared Libraries\redist\ia32\compiler

Now, my problem is that I would really like to use the routines potrf, potri in addition to getrf, getri (just to broaden my toolbag). I thought I had this working in a small demo program, but it doesn't work now. The program:

   PROGRAM MKLinvert
      !include 'lapack.f90'
      !USE IFQWIN
      USE LAPACK95 !, only: getrf, getri, potrf, potri
      !use mkl95_lapack ! appears in Intel mkl samples; does same thing
      USE F95_PRECISION, ONLY: WP => SP ! or DP
      IMPLICIT NONE
      INTEGER ipiv(2) ! NOTE must have length >= A,B,C
      !REAL(WP)  :: A(2, 2), B(2, 2), C(2, 2)
      real(4) a(2,2), b(2,2), c(2,2) ! same result
     ! potrf,potri?
      character(1) uplo
      integer info
      
!---- Define an input matrix.      
      
      A(1,1) = 10.
      A(1,2) = 3.
      A(2,1) = A(1,2)
      A(2,2) = 5.
      print *, 'Input matrix:'
      print *, 'A0 =          ', A ! OK
      print *
      
!---- Call MKL functions to factor, then inverse.      

      ! Can invert in place with following calls,
      !CALL getrf (A, ipiv, info)
      !CALL getri (A, ipiv)
      ! But we will work on a copy to verify the result.
      ! Also demonstrate absence of optional parameters.
      !B = A_inv
      
      B = A
      print *, 'A1 =          ', A
      CALL getrf (B, ipiv)
      print *, 'A2 =          ', A
      CALL getri (B, ipiv)
      print *, 'A3 =          ', A 
      ! Note that ipiv is optional in getrf according to the documentation, but
      ! the result is wrong unless it is included. In getri it is (correctly)
      ! shown to be necessary. IPIV IS THE PIVOT INFORMATION, WHICH COMES FROM
      ! GETRF AND IS NEEDED BY GETRI. So it is optional when calling getrf only
      ! to factor the matrix, but is required when followed by getri to do the
      ! inversion.
      
!---- Show the result and check it.      

      PRINT *
      PRINT *, 'Matrix inversion using MKL functions getrf, getri'
      PRINT *, 'A =        ', A
      PRINT *, 'B = A_inv =', B
      PRINT *, 'A * B =    ', MATMUL (A, B)
      PRINT *, 'Identiy matrix result indicates inversion was correct.'
      
!---- Repeat the problem using routines potrf and potri. These are for a
!     positive-definite symmetric matrix. According to the documentation the no.
!     of floating-point operations for potri is about half of getri.

      C = A
      CALL potrf (C) ! Upper tri part of C is returned as Cholesky factor U
      CALL potri (C)
      !CALL potri (C, 'U', info) ! No help. info returns 0 (ok).
      
      PRINT *
      PRINT *, 'Result using potrf, potri'
      PRINT *, 'A =        ', A
      PRINT *, 'C = A_inv =', C
      PRINT *, 'A * C =    ', MATMUL (A, C)
      print *, 'Result using potrf, potri may be wrong. Don''t know why.'
      

   END PROGRAM MKLinvert

The result from getri are correct, as seen by the identity matrix resulting from multiplying the original matrix by its inverse:

A * Ainv = 1.    0.
                 0.    1.
The result from potri is not correct--not even close:

A* Ainv = 10.2     15.4
                  ~0         1.

 

0 Kudos
mecej4
Honored Contributor III
1,512 Views

You need to pass only the upper triangle to POTRF/POTRI. The strictly lower half may be used as a scratch pad by Lapack, and any values in that part that you pass are not used.

Only the upper triangle of the inverse is filled in and returned by POTRI. You should not use any of the elements in the strictly lower half, which may have their original values, or may have been overwritten with temporary results (scratchpad values). If you wish to use the whole inverse matrix, do so after copying the strictly upper half to the strictly lower half.

In your program, after calling POTRI, add " C(2,1)=C(1,2) " before passing matrix C to MATMUL or printing matrix C.

Questions regarding MKL/Lapack/BLAS tend to get better responses in the MKL forum, https://software.intel.com/en-us/forums/intel-math-kernel-library .

0 Kudos
dboggs
New Contributor I
1,512 Views

Thank you mecej4, your explanation and instruction are perfect.

But wow, how is a new user supposed to know this? The documentation clearly says that potri(A) returns the inverse of A, which is simply not true. It could easily have said that it returns the upper triangular part of Ainv, with junk in the lower triangular part. There isn't even anything that hints it is the upper triangle that is good, not the lower triangle! This is disastrous documentation. 

0 Kudos
mecej4
Honored Contributor III
1,512 Views

The MKL documentation for POTRF says, "If uplo = 'U', the array a stores the upper triangular part of the matrix A, and the strictly lower triangular part of the matrix is not referenced.", and the page for POTRI does say, "a     Overwritten by the upper or lower triangle of the inverse of A."

It is usually necessary to read the documentation in detail, especially at the first occasion of using a routine. Furthermore, given the number of different matrix types that are covered in Lapack, it is not surprising that different representations of a matrix are used in different sections of Lapack. The representation of a matrix as a rectangular array is just one such representation, and this representation is rarely used for symmetric matrices in the routines of the ?POxxx family.

In early mainframes, memory was very limited, and triangular storage consumes about 1/2 of the memory that would be consumed if the full matrix were stored.

0 Kudos
dboggs
New Contributor I
1,513 Views

When I complained about the inaccuracy of "the documentation," I refer to what you get when you click on "potri" in your source and press F1, bringing up a documentation brief which, allegedly (as in all of the usual Fortran subroutines) gives sufficient information to actually use the thing. In this case, it clearly describes the return value "a" as "overwritten by the n-by-n matrix inv(A). This tells me that a will contain--well--the complete inverse matrix. It does not. This documentation description is not only incomplete, it is dangerously wrong. Would have been a very simple correction to inform the user that a contains only a triangular part of the matrix and to ignore the other triangular part. 

But I have learned my lesson. If I want to use MKL for additional routines in the future, I must either (1) spend hours studying the rest of the "complete documentation" or (2) take my best understanding of the subroutine brief, try it, test it to see if it gives the right answer in a known case, and then go to this forum if it does not. Thankfully there are people like mecej4 who are willing to help and in quick order too.

0 Kudos
Reply