Problem compiling Fortran code with MKL

raposo__henrique · ‎02-04-2018

Hi all,

I'm trying to replace my own LU decomposition code with the fortran subroutine gbtrf from MKL. Please bear with me as I am somewhat new to Fortran.

I successfully compiled, run and validated the example (sgbtrfx.f) provided by Intel (2018 release). The compilation was done with:

gfortran -w source/sgbtrfx.f -Wl,--start-group /opt/intel/compilers_and_libraries_2018/linux/mkl/lib/intel64/libmkl_intel_lp64.a /opt/intel/compilers_and_libraries_2018/linux/mkl/lib/intel64/libmkl_intel_thread.a /opt/intel/compilers_and_libraries_2018/linux/mkl/lib/intel64/libmkl_core.a -Wl,--end-group lib/libaux_intel64_gfortran.a -L/opt/intel/compilers_and_libraries_2018/linux/mkl/../compiler/lib/intel64 -liomp5 -lpthread -ldl -o _results/intel_lp64_omp_intel_iomp5_intel64_lib/sgbtrfx.out

And it also required me to set:

export LD_LIBRARY_PATH=/opt/intel/compilers_and_libraries_2018.1.163/linux/compiler/lib/intel64_lin

I am now finding some difficulties getting the same subroutine to work with some code I previously had. My compilation and linking were done with:

1. Compilation

gfortran -O1 -fopenmp -ffast-math -funroll-loops --param max-unroll-times=2 -Ofast -march=native -mcmodel=medium -fno-range-check -ffree-line-length-0 -fdefault-double-8 -fdefault-real-8 -fdefault-integer-8 -w -c <FILENAME>.f90

2. Linking

gfortran -O1 -fopenmp -ffast-math -funroll-loops --param max-unroll-times=2 -Ofast -march=native -mcmodel=medium -fno-range-check -ffree-line-length-0 -fdefault-double-8 -fdefault-real-8 -fdefault-integer-8 -o ......

And I got it to link to the MKL library using:

1. Compilation

gfortran -O1 -fopenmp -m64 -I/opt/intel/parallel_studio_xe_2018/compilers_and_libraries_2018/linux/mkl/include -ffast-math -funroll-loops --param max-unroll-times=2 -Ofast -march=native -mcmodel=medium -fno-range-check -ffree-line-length-0 -fdefault-double-8 -fdefault-real-8 -fdefault-integer-8 -w -c <FILENAME>.f90

2. Linking

gfortran -O1 -fopenmp -L/opt/intel/parallel_studio_xe_2018/compilers_and_libraries_2018/linux/mkl/lib/intel64 -Wl,--no-as-needed -lmkl_gf_ilp64 -lmkl_gnu_thread -lmkl_core -lgomp -lpthread -lm -ldl -ffast-math -funroll-loops --param max-unroll-times=2 -Ofast -march=native -mcmodel=medium -fno-range-check -ffree-line-length-0 -fdefault-double-8 -fdefault-real-8 -fdefault-integer-8 -o ......

as well as

export LD_LIBRARY_PATH=/opt/intel/parallel_studio_xe_2018/compilers_and_libraries_2018/linux/mkl/lib/intel64

But when I call the subroutine sgbtrf with the exact same input as in the example file (sgbtrfx.f) I don't obtain the correct output.

Can anyone help me figure out what might be wrong here?

Many thanks,

Henrique

mecej4 · ‎02-04-2018

Almost all Lapack and BLAS routines, such as ?gbtrf, are deterministic and use no global variables. In other words, all information is passed through the argument list, and the codepath chosen through the MKL routines is typically determined by the CPU family ID.

My guess is that you are either not correctly passing the arguments or have used incompatible versions of the MKL libraries. To say more, I would need to see example code that does what you say it does. For starters, how did you verify that the correct arguments were passed?

raposo__henrique · ‎02-04-2018

Ok so here's an altered version of the code provided in MKL's example:

!SUBROUTINE INVERT()
PROGRAM MATRIXINV

      INTEGER          NIN, NOUT
      PARAMETER        (NIN=5,NOUT=6)
      INTEGER          MMAX, NMAX, KLMAX, KUMAX, LDAB
      PARAMETER        (MMAX=5,NMAX=5,KLMAX=3,KUMAX=2,LDAB=2*KLMAX+KUMAX+1)
!     .. Local Scalars ..
      INTEGER          I, INFO, J, K, KL, KU, M, N,num_el,j_min,j_max,LDB
!     .. Local Arrays ..
      REAL             AB(LDAB,NMAX), B(NMAX)
      REAL	       AUX(KLMAX+KUMAX+1)
      INTEGER          IPIV(NMAX)
!     .. External Subroutines ..
      EXTERNAL         SGBTRF
      EXTERNAL         SGBTRS
!     .. Intrinsic Functions ..
      INTRINSIC        MAX, MIN
!     .. Executable Statements ..
      WRITE (NOUT,*) 'SGBTRF Example Program Results'
!     Skip heading in data file
      OPEN(UNIT=10, FILE="sgbtrfx.d", STATUS='OLD')


      READ(10,*,END=12,IOSTAT=JFK,ERR=1099)
      READ(10,*,END=12,IOSTAT=JFK,ERR=1099) M, N, KL, KU

      IF (M.LE.MMAX .AND. N.LE.NMAX .AND. KL.LE.KLMAX .AND. KU.LE.KUMAX) THEN
!        Read A from data file

         K = KL + KU + 1
	  DO I=1,M
           num_el=K
	       j_min=i-kl
	       j_max=i+ku
	       IF(i+KU>N) then
             num_el=num_el-(i+ku-n)
	         j_max=N
           end if
	       IF(i-KL<=0) then
              num_el=num_el+i-kl-1
	          j_min=1
           end if
           
           READ(10,*,END=12,IOSTAT=JFK,ERR=1099) AUX(1:num_el)

           DO j=j_min,j_max
             AB(kl+ku+1+i-j,j)=aux(j-j_min+1)
           END DO
      END DO
	 

!        Factorize A
       CALL SGBTRF(M,N,KL,KU,AB,LDAB,IPIV,INFO)

         ! Fill in the RHS vector
       DO i=1,N
         b(i)=float(i)
	   end do
	   LDB=NMAX
!		 ! Solve linear system
       CALL SGBTRS('N',M,kl,ku,1,ab,ldab,ipiv,b,ldb,info)
       print*,b
!        Print details of factorization
       WRITE (NOUT,*)

       CALL PRINT_BAND_MATRIX( 'Details of factorization', M, N, KL,KL+KU, AB, LDAB )

!        Print pivot indices
       WRITE (NOUT,*)
       WRITE (NOUT,*) 'IPIV'
       WRITE (NOUT,99999) (IPIV(I),I=1,MIN(M,N))

       IF (INFO.NE.0) WRITE (NOUT,*) 'The factor U is singular'

      END IF
      STOP

99999 FORMAT ((3X,7I11))
1099  CONTINUE
12    CONTINUE
      END

!     End of SGBTRF Example

!  =============================================================================

!     Auxiliary routine: printing a banded matrix stored in packed form.

      SUBROUTINE PRINT_BAND_MATRIX( DESC, M, N, KL, KU, A, LDA )
      CHARACTER*(*)    DESC
      INTEGER          M, N, KL, KU, LDA
      REAL             A( LDA, * )
      REAL             A_B( M, N )
      REAL             ZERO
      PARAMETER        (ZERO=0.0)

      INTEGER          I, J, LB, UB

      DO J = 1, N
         LB = MAX(1, J - KU)
         UB = MIN(M, J + KL)
         DO I = 1, M
            IF ((LB.LE.I).AND.(I.LE.UB)) THEN
               A_B( I, J ) = A( KU + 1 + I - J, J )
            ELSE
               A_B( I, J) = ZERO
            END IF
         END DO
      END DO

      WRITE(*,*) DESC
      WRITE(*, 9999) ( J, J = 1, N)
      DO I = 1, M
         WRITE(*, 9998) I, ( A_B( I, J ), J = 1, N )
      END DO

 9998 FORMAT( I2, ' ', 11(:,1X,F10.4) )
 9999 FORMAT( '   ', 11(:,1X,I10) )

      RETURN
      END

This works fine when compiling with:

gfortran -w MATRIXINV.f90 -Wl,--start-group /opt/intel/compilers_and_libraries_2018/linux/mkl/lib/intel64/libmkl_intel_lp64.a /opt/intel/compilers_and_libraries_2018/linux/mkl/lib/intel64/libmkl_intel_thread.a /opt/intel/compilers_and_libraries_2018/linux/mkl/lib/intel64/libmkl_core.a -Wl,--end-group lib/libaux_intel64_gfortran.a -L/opt/intel/compilers_and_libraries_2018/linux/mkl/../compiler/lib/intel64 -liomp5 -lpthread -ldl -o matrixinv.out

The file sgbtrfx.d is the one provided in MKL's examples. It contains a 4x4 matrix (I can check wether the input is correct with a simple print since the matrix is so small):

SGBTRF (F07BDE) Example Program Data
  4  4  1  2                  :Values of M, N, KL and KU
 -0.23   2.54  -3.66
 -6.98   2.46  -2.73  -2.13
         2.56   2.46   4.07
               -4.78  -3.82   :End of matrix A

My problem really boils down to not knowing how to compile and link this bit of code to a program. Suppose the first line is uncommented (and the second commented) and that I have another file with

PROGRAM MASTER

CALL INVERT()

END

I basically just want to know what the correct way of compiling and linking these two files is! Should be a simple matter and yet I can't seem to make it work. When using the lines for compilation and linking that I provided in my original post I get the wrong result.

Thanks again

mecej4 · ‎02-04-2018

You show different compiler options in #1 and #3, and I am not sure which one you actually used. For example, you specified 8-byte integers in #1 and default (4-byte) integers in #3. You have to be quite correct with such options.

Furthermore, I suggest that for the first few runs you do not use threading, optimization, etc. For linking, please consult the MKL link line advisor at https://software.intel.com/en-us/articles/intel-mkl-link-line-advisor .

On my openSuse Linux system, I set MKLROOT and LD_LIBRARY_PATH, placed the source codes into files rmain.f90 and rpso.f90, and built the program with the command

gfortran -I${MKLROOT}/include rmain.f90 rpso.f90   -L${MKLROOT}/lib/intel64 -Wl,--no-as-needed -lmkl_gf_lp64 -lmkl_sequential -lmkl_core -lpthread -lm -ldl

The program ran fine and gave correct results.

TimP · ‎02-05-2018

I agree with Sham's suggestion you start with minimal command line options and later check the effect of adding complications incrementally. It's hard to see through your forest, but you must choose either the libiomp5 or the libgomp and use exclusively the corresponding mkl libraries.

raposo__henrique · ‎02-05-2018

Thanks a lot, I now understand what exactly is going wrong.

When adding the flag -fdefault-real-8 , the MKL subroutine no longer gives me the right output. This actually means that whenever I declare 8 byte reals and provide them as inputs, I don't get the correct output.

Do you know why this is? Shouldn't MKL work with different precisions?

Thanks for your patience!

raposo__henrique · ‎02-05-2018

Thanks a lot, I now understand what exactly is going wrong.

When adding the flag -fdefault-real-8 , the MKL subroutine no longer gives me the right output. This actually means that whenever I declare 8 byte reals and provide them as inputs, I don't get the correct output.

Do you know why this is? Shouldn't MKL work with different precisions?

Thanks for your patience!

mecej4 · ‎02-05-2018

No, I think that you are jumping to incorrect conclusions.

Most MKL routines (especially those that are part of BLAS and Lapack) have versions for single- and double-precision real, and single- and double precision complex arguments. These have slightly different names, e.g., SDOT, DDOT, CDOT and ZDOT, commonly designated by ?DOT. With Pardiso, you have to specify whether you are passing single- or double -precision REAL arguments by setting appropriate elements of IPARM.

When you call a routine using its implicit interface, be careful about compiler options (such as -fdefault-integer-8) that change the types of default reals and integer variables. The MKL library routines will work only if they are given arguments of the expected types. If you are compiling for 64-bit targets, integer arguments may be 4- or 8-byte integers, but you have to use the *lp64 libraries in the former case and *ilp64 libraries in the latter case. If you do not do this matching correctly, the resulting program will often crash with no warning messages.

In short, I think that your errors are more probably related to integer arguments than real arguments.

prabha_d_ · ‎02-05-2018

Hi ,

I bought license for Intel parallel studio XE composer edition 2018 for windows . I have installed the same with math kernel library 2018 update 1. I have my code in Fortran and using it in Visual Studio 2017 with the intel fortran compiler. But I am getting error while build a program :

Unresolved exernal Symbol_DIVPAG, UMACH , UMPOL etc . which I understand MKL is not liking with that.

Kindly help me to resolve this problem .

Thanks in advance

Prabha

mecej4 · ‎02-06-2018

The routines named in #9 are from the IMSL library, and they do not exist in MKL. You will need to consider purchasing an IMSL license or replacing the calls to the IMSL routines with alternatives from other sources.

Prabha, it is not good etiquette to post unrelated responses in an existing thread. Please start a new thread for your question.

prabha_d_ · ‎02-06-2018

Thank you so much for you advice

raposo__henrique · ‎02-07-2018

Yes I should have been using the *ilp64 libraries. But most of all I had to use DGBTRF and not SGBTRF for the double precision case to work.

Anyway, all is working smoothly now. Thank you for enlightening me!