Problems with Lapack LS solver

[fercook]

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.

[TimP]

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.

[fercook]

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?

[fercook]

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

[mecej4]

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]