topic Problems with Lapack LS solver in Intel® oneAPI Math Kernel Library

Problems with Lapack LS solver

fercook — Thu, 04 Nov 2010 14:13:30 GMT

Hello,

I am running into a SIGSEGV segmentation fault when I try to use the *GELS* routines for Least Squares problems. I have a small reproducer code in F90 (I am using it on Mac OS X):

[fortran]program testlapack

    use mkl95_lapack
    use mkl95_precision

    implicit none
    integer,parameter :: left=5,right=2
    complex(8),parameter :: II=(0.0d0,1.0d0)
    real(8),parameter :: Tolerance=1.0d-9
    real(8) :: tempR(left,left),tempI(left,left),vecR(left,right),vecI(left,right)
    complex(8) :: matrixdata(left,left),vectordata(left,right)
    integer :: error

    print *,'Prepare data'
    call random_number(tempR)
    call random_number(tempI)
    matrixdata=tempR+II*tempI
    call random_number(vecR)
    call random_number(vecI)
    vectordata=vecR+II*vecI

    print *,'Call MKL Lapack solver'
    call ZGELSD( matrixdata, vectordata, Tolerance )

    print *,'Solution obtained:'
    print *,vectordata

end program testlapack
[/fortran]

which I compile using

[bash]ifort -L$MKLPATH -I$MKLINCLUDE -I$MKLINCLUDE/em64t/lp64 -lmkl_blas95_lp64 $MKLPATH/libmkl_intel_lp64.a $MKLPATH/libmkl_intel_thread.a $MKLPATH/libmkl_core.a $MKLPATH/libmkl_intel_lp64.a $MKLPATH/libmkl_intel_thread.a $MKLPATH/libmkl_core.a -liomp5 -lpthread testlapack.f90 -g -traceback[/bash]

It compiles ok but gives a seg fault when calling ZGELSD.

Am I using the F95 interface correctly? Am I using the right combination for MKL (threaded or not, static or dynamic, etc)? I know there are some routines that need to be linked static but *GELS are not mentioned, so I assumed they could be linked dynamically.

Problems with Lapack LS solver

TimP — Thu, 04 Nov 2010 14:27:09 GMT

It does look risky to link as you have done. All MKL libraries dynamic is a reasonable starting point. For a next step, all libraries static should give you a link time error if any references are missing. You would need to specify the f95 wrapper library just once, ahead of the other MKL static libraries, which should preferably be done as the link advisor suggests.

Problems with Lapack LS solver

fercook — Thu, 04 Nov 2010 15:01:21 GMT

OK, I am going to lift the suggested linker sequences from the MKL User's Guide:

1) Static linking and parallel MKL

[bash]ifort -L$MKLPATH -I$MKLINCLUDE $MKLPATH/libmkl_intel_lp64.a $MKLPATH/libmkl_intel_thread.a $MKLPATH/libmkl_core.a $MKLPATH/libmkl_intel_lp64.a $MKLPATH/libmkl_intel_thread.a $MKLPATH/libmkl_core.a -liomp5 -lpthread testlapack.f90[/bash]

2) All dynamic linking and threaded:

[bash]ifort -L$MKLPATH -I$MKLINCLUDE -lmkl_intel_ilp64 -lmkl_intel_thread -lmkl_core -liomp5 -lpthread testlapack.f90[/bash]

3)Static linking and sequential MKL

[bash]ifort -L$MKLPATH -I$MKLINCLUDE $MKLPATH/libmkl_intel_lp64.a $MKLPATH/libmkl_sequential.a $MKLPATH/libmkl_core.a $MKLPATH/libmkl_intel_lp64.a $MKLPATH/libmkl_sequential.a $MKLPATH/libmkl_core.a -lpthread testlapack.f90[/bash]

4) Dynamic linking and sequential MKL

[bash]ifort -L$MKLPATH -I$MKLINCLUDE -lmkl_intel_lp64 -lmkl_sequential -lmkl_core -lpthread testlapack.f90[/bash]

5) Static linking and parallel MKL

[bash]ifort -L$MKLPATH -I$MKLINCLUDE $MKLPATH/libmkl_intel_ilp64.a $MKLPATH/libmkl_intel_thread.a $MKLPATH/libmkl_core.a $MKLPATH/libmkl_intel_ilp64.a $MKLPATH/libmkl_intel_thread.a $MKLPATH/libmkl_core.a -liomp5 -lpthread testlapack.f90 [/bash]

6) Dynamic linking and parallel MKL with ILP interface

[bash]ifort -L$MKLPATH -I$MKLINCLUDE -lmkl_intel_ilp64 -lmkl_intel_thread -lmkl_core -liomp5 -lpthread testlapack.f90[/bash]

7) Static linking, F95 Lapack interface, and threaded MKL

[bash]ifort -L$MKLPATH -I$MKLINCLUDE -I$MKLINCLUDE/em64t/lp64 -lmkl_lapack95_lp64 $MKLPATH/libmkl_intel_lp64.a $MKLPATH/libmkl_intel_thread.a $MKLPATH/libmkl_core.a $MKLPATH/libmkl_intel_thread.a $MKLPATH/libmkl_core.a $MKLPATH/libmkl_intel_thread.a $MKLPATH/libmkl_core.a -liomp5 -lpthread testlapack.f90[/bash]

They ALL compile ok and give seg fault.

I don't know another way to compile it, and at this point I don't think it is a problem of how I am compiling. Perhaps the code is not as it should be?

Problems with Lapack LS solver

fercook — Thu, 04 Nov 2010 15:23:44 GMT

I like nothing more than answering my own posts.

I switched to the F77 routines and they work fine. Here is the new version:

[fortran]program testlapack

    implicit none
    integer,parameter :: left=5,right=2
    complex(8),parameter :: II=(0.0d0,1.0d0)
    real(8),parameter :: Tolerance=1.0d-8
    real(8) :: tempR(left,left),tempI(left,left),vecR(left,right),vecI(left,right)
    complex(8) :: matrixdata(left,left),vectordata(left,right)
    real(8) :: GELOutputS(left)
    integer :: error, rank,info,tempDim
    complex(8),allocatable :: work(:)
    real(8),allocatable :: rwork(:)
    integer,allocatable :: iwork(:)
    integer :: lwork,lrwork,liwork

    print *,'Prepare data'
    call random_number(tempR)
    call random_number(tempI)
    matrixdata=tempR+II*tempI
    call random_number(vecR)
    call random_number(vecI)
    vectordata=vecR+II*vecI

    print *,'Call MKL Lapack solver'
    print *,'First find out optimum lwork'
    tempDim=1
    allocate (work(tempDim),rwork(tempDim),iwork(tempDim))
    lwork=-1
    call ZGELSD ( left, left, right, matrixdata, left, vectordata, left, GELOutputS, Tolerance, rank,work,lwork,rwork,iwork,info)

    !Reallocate vectors with optimum workspace
    lwork=int(work(1))
    lrwork=int(rwork(1))
    liwork=iwork(1)
    deallocate(work,rwork,iwork)
    allocate(work(lwork),rwork(lrwork),iwork(liwork))

    call ZGELSD ( left, left, right, matrixdata, left, vectordata, left, GELOutputS,Tolerance,rank,work,lwork,rwork,iwork,info)
    print *,'Solution obtained:',info
    print *,vectordata

end program testlapack
[/fortran]

Interestingly, it DID NOT work with the threaded MKL and ILP, only with LP64 interface.

Best,

Fer

Problems with Lapack LS solver

mecej4 — Fri, 05 Nov 2010 01:47:35 GMT

The reason that your F9x version crashed is that you called the F77 subroutine name with improper arguments. The modules that you listed under USE do not contain the F77 interfaces, as a consequence of which no argument checks were made on the ZGELSD call. Here is a corrected version of your program. In this version, I call the F95 interface name GELS, for which the interface is included in mkl95_lapack.mod.

[fortran]program testlapack

    use mkl95_lapack

    implicit none
    integer,parameter :: left=5,right=2
    complex(8),parameter :: II=(0.0d0,1.0d0)
    real(8),parameter :: Tolerance=1.0d-9
    real(8) :: tempR(left,left),tempI(left,left),vecR(left,right),vecI(left,right)
    complex(8) :: matrixdata(left,left),vectordata(left,right)
    integer :: error

    print *,'Prepare data'
    call random_number(tempR)
    call random_number(tempI)
    matrixdata=tempR+II*tempI
    call random_number(vecR)
    call random_number(vecI)
    vectordata=vecR+II*vecI

    print *,'Call MKL Lapack solver'
    call GELS( matrixdata, vectordata )

    print *,'Solution obtained:'
    print *,vectordata

end program testlapack
[/fortran]