Error in Lapack routine SGELS using an interface and allocatable WORK array

Question

Problem:

Using an explicit interface for Lapack (http://www.netlib.org/lapack/) routines simplifies coding. There is a problem with routines like SGELS (https://www.netlib.org/lapack/explore-3.1.1-html/sgels.f.html) which require a WORK (https://www.netlib.org/lapack/lug/node117.html) array of variable size. Fixing the size of the WORK array to a large number is successful, but defeats the economy of Lapack.

But I am unable to use an allocatable array for WORK.

The usual error message is something like

Program received signal SIGSEGV: Segmentation fault - invalid memory reference.

Fortran query routines on the WORK array such as shape, size, is_contiguous and bounds checking are successful. But any attempt to print a value from the array produces a segmentation fault.

After several permutations, here is a minimal working example of the code.

Code

Main program

program mwe

    use mLapackInterfaceSGELS, only : sgels
    implicit none

    integer, parameter :: n = 2

    real, dimension ( : ),     allocatable :: work
    real, dimension ( : , : ), allocatable :: AstarA, Astarb

    integer :: nrhs = 0, lda = 0, ldb = 0, lwork = 0, info = 0

        lda = n
        ldb = n
        nrhs = 1
        lwork = -1

        allocate ( Astarb ( 1 : n, 1 : 1 ) )
        allocate ( AstarA ( 1 : n, 1 : n ) )
        allocate ( work ( 1 : 2 ) )

        Astarb ( 1 : n, 1 )     = [ 466.7, 2898.]
        AstarA ( 1 : n, 1 : n ) = reshape ( [ [ 9., 45. ], [ 45., 285.] ], [ n, n ] )
        work ( 1 : 2 ) = 0.0

        write ( * , * ) "work array before query = ", work

        call sgels( trans = 'No transpose', m = n, n = n, nrhs = nrhs, A = AstarA, lda = lda, b = Astarb, ldb = ldb, &
                    work = work, lwork = -1, info = info )

        write ( * , * ) "out: shape ( work ) = ", shape ( work )
        write ( * , * ) "lbound ( work ) = ", lbound ( work )
        write ( * , * ) "ubound ( work ) = ", ubound ( work )
        write ( * , * ) "size ( work ) = ", size ( work )
        write ( * , * ) "is_contiguous ( work ) = ", is_contiguous ( work )
        write ( * , * ) "work array after query = ", work

    stop
end program mwe

Interface module

(As seen in Numerical Methods with Modern Fortran, Richard Hanson, Tim Hopkins, Listing 2.5. (https://my.siam.org/Store/Product/viewproduct/?ProductId=24372445))

module mLapackInterfaceSGELS

    implicit none

    interface lapack_sgels

        subroutine sgels ( trans, m, n, nrhs, A, lda, b, ldb, work, lwork, info )
            integer,           intent ( in )    :: m, n, nrhs, lda, ldb, lwork
            integer,           intent ( out )   :: info
            real, allocatable, intent ( out )   :: work ( : )
            real,              intent ( inout ) :: A ( 1 : lda , 1 : n ), b ( 1 : ldb , 1 : nrhs )
            character,         intent ( in )    :: trans
        end subroutine sgels

    end interface lapack_sgels

end module mLapackInterfaceSGELS

Output from code

The line which causes the error is:

write ( * , * ) "work array after query = ", work

The output from the code is:

 work array before query =    0.00000000       0.00000000    
 out: shape ( work ) =            2
 lbound ( work ) =            1
 ubound ( work ) =            2
 size ( work ) =            2
 is_contiguous ( work ) =  T

Program received signal SIGSEGV: Segmentation fault - invalid memory reference.

Compiler version

% gcc --version                                                                        (1c01947)fortran-alpha
gcc (Spack GCC) 10.2.0

Compilation

gfortran -g -c -Og -pedantic -Wall -Warray-temporaries -Wextra -Waliasing -Wsurprising -Wimplicit-procedure -Wintrinsics-std -Wfunction-elimination -Wc-binding-type -Wrealloc-lhs-all -Wuse-without-only -Wconversion-extra -fno-realloc-lhs -ffpe-trap=denormal,invalid,zero -fbacktrace -fmax-errors=5 -fcheck=all -fcheck=do -fcheck=pointer -fno-protect-parens -faggressive-function-elimination -fdiagnostics-color=auto -finit-derived -o m-lapack-interface-sgels.o m-lapack-interface-sgels.f08
gfortran -g -c -Og -pedantic -Wall -Warray-temporaries -Wextra -Waliasing -Wsurprising -Wimplicit-procedure -Wintrinsics-std -Wfunction-elimination -Wc-binding-type -Wrealloc-lhs-all -Wuse-without-only -Wconversion-extra -fno-realloc-lhs -ffpe-trap=denormal,invalid,zero -fbacktrace -fmax-errors=5 -fcheck=all -fcheck=do -fcheck=pointer -fno-protect-parens -faggressive-function-elimination -fdiagnostics-color=auto -finit-derived -o mwe.o mwe.f08

Execution

Two different Lapack libraries were used. Both produced the same failure.

/usr/local/lib

gfortran -g -L/usr/local/lib -llapack -lblas -o mwe m-lapack-interface-sgels.o mwe.o

Apple Accelerate framework

gfortran -g -framework Accelerate -o mwe m-lapack-interface-sgels.o mwe.o

Tip: Consider redirecting the output to avoid overrunning the screen buffer:

./mwe 2>&1 | tee output.txt

Question

How do I fix the interface and the call to allow the use of an allocatable WORK array?

Reference:

Parameter corruption when interfacing FORTRAN lapack routine

Segmentation error in fortran LAPACK routine

user_1280439 · Answer

You seem to be somewhat confused about how LAPACK routines work. Showing their vintage they never do any memory management for you, you have to allocate all memory yourself before you call the routines. Thus a dummy argument can never be allocatable for a LAPACK routine - and I have to say I would find it hard to trust any book that claims the above is the correct interface.

So your interface is wrong. But also are you confused about how workspace queries work in LAPACK? The routine doesn't allocate any memory for you - you call it once with the workspace argument sent to -1 to indicate that this is a workspace query, as documented at sgles documentation page, and this call will return the required size of the workspace array in work( 1 ). You then allocate the work array to the correct size, and finally call the routine again with the correct arguments for the workspace. It's all a bit clunky, again showing LAPACK's age, but I hope the above is clear

Below is what I think is a correct version of what you want to do. Note I have made a couple of other small changes to your LAPACK interface to make it what I think is the fully correct version, and I have also removed the unnecessary initialisation of a lot of your variables - this is a bad habit to get into, it can mask bugs, or at least make them harder to find, and worse is often more or less incompatible with threaded programming, an important form of parallelism that I personally think all Fortran programmers should have some familiarity with.

ijb@ijb-Latitude-5410:~/work/stack$ cat sgles.f90
module mLapackInterfaceSGELS

  implicit none

  interface lapack_sgels

     subroutine sgels ( trans, m, n, nrhs, A, lda, b, ldb, work, lwork, info )
       integer,   intent ( in )    :: m, n, nrhs, lda, ldb, lwork
       integer,   intent ( out )   :: info
       real,      intent ( out )   :: work ( * )
       real,      intent ( inout ) :: A ( 1 : lda , * ), b ( 1 : ldb , * )
       character, intent ( in )    :: trans
     end subroutine sgels

  end interface lapack_sgels

end module mLapackInterfaceSGELS

program mwe

  use mLapackInterfaceSGELS, only : sgels
  implicit none

  integer, parameter :: n = 2

  real, dimension ( : ),     allocatable :: work
  real, dimension ( : , : ), allocatable :: AstarA, Astarb

  Real, Dimension( 1:1 ) :: r_lwork 
  
  integer :: nrhs, lda, ldb, lwork, info

  lda = n
  ldb = n
  nrhs = 1
  lwork = -1

  allocate ( Astarb ( 1 : n, 1 : 1 ) )
  allocate ( AstarA ( 1 : n, 1 : n ) )

  Astarb ( 1 : n, 1 )     = [ 466.7, 2898.]
  AstarA ( 1 : n, 1 : n ) = reshape ( [ [ 9., 45. ], [ 45., 285.] ], [ n, n ] )

  call sgels( trans = 'No transpose', m = n, n = n, nrhs = nrhs, A = AstarA, lda = lda, b = Astarb, ldb = ldb, &
       work = r_lwork, lwork = -1, info = info )
  lwork = Nint( r_lwork( 1 ) )
  Allocate( work( 1:lwork ) )
  call sgels( trans = 'No transpose', m = n, n = n, nrhs = nrhs, A = AstarA, lda = lda, b = Astarb, ldb = ldb, &
       work = work, lwork = lwork, info = info )

  write ( * , * ) "out: shape ( work ) = ", shape ( work )
  write ( * , * ) "lbound ( work ) = ", lbound ( work )
  write ( * , * ) "ubound ( work ) = ", ubound ( work )
  write ( * , * ) "size ( work ) = ", size ( work )
  write ( * , * ) "is_contiguous ( work ) = ", is_contiguous ( work )
  write ( * , * ) "work array after query = ", work

  stop
end program mwe
ijb@ijb-Latitude-5410:~/work/stack$ gfortran --version
GNU Fortran (Ubuntu 9.3.0-17ubuntu1~20.04) 9.3.0
Copyright (C) 2019 Free Software Foundation, Inc.
This is free software; see the source for copying conditions.  There is NO
warranty; not even for MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.

ijb@ijb-Latitude-5410:~/work/stack$ gfortran -Wall -Wextra -fcheck=all -std=f2018 -O -g sgles.f90 -lopenblas
ijb@ijb-Latitude-5410:~/work/stack$ ./a.out
 out: shape ( work ) =           66
 lbound ( work ) =            1
 ubound ( work ) =           66
 size ( work ) =           66
 is_contiguous ( work ) =  T
 work array after query =    66.0000000       0.00000000       4192.00000       0.00000000       0.00000000       0.00000000       0.00000000       0.00000000       0.00000000       0.00000000       0.00000000       0.00000000       0.00000000       0.00000000       0.00000000       0.00000000       0.00000000       0.00000000       0.00000000       0.00000000       0.00000000       0.00000000       0.00000000       0.00000000       0.00000000       0.00000000       0.00000000       0.00000000       0.00000000       0.00000000       0.00000000       0.00000000       0.00000000       0.00000000       0.00000000       0.00000000       0.00000000       0.00000000       0.00000000       0.00000000       0.00000000       0.00000000       0.00000000       0.00000000       0.00000000       0.00000000       0.00000000       0.00000000       0.00000000       0.00000000       0.00000000       0.00000000       0.00000000       0.00000000       0.00000000       0.00000000       0.00000000       0.00000000       0.00000000       0.00000000       0.00000000       0.00000000       0.00000000       0.00000000       0.00000000       0.00000000