时间戳

首页 > 解决方案 > 使用接口和可分配 WORK 数组的 Lapack 例程 SGELS 中的错误

fortran - 使用接口和可分配 WORK 数组的 Lapack 例程 SGELS 中的错误

问题描述

问题：

为 Lapack ( http://www.netlib.org/lapack/ ) 例程使用显式接口可以简化编码。像 SGELS ( https://www.netlib.org/lapack/explore-3.1.1-html/sgels.f.html ) 这样需要WORK( https://www.netlib.org/lapack ) 的例程存在问题/lug/node117.html ) 可变大小的数组。将数组的大小固定WORK为一个大数是成功的，但会破坏 Lapack 的经济性。

但我无法将allocatable数组用于WORK.

通常的错误信息类似于

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

数组上的 Fortran 查询例程（WORK例如shape、size和is_contiguous边界检查）是成功的。但是任何从数组中打印值的尝试都会产生分段错误。

经过几次排列后，这是代码的最小工作示例。

代码

主程序

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

接口模块

（如现代 Fortran 的数值方法中所见，Richard Hanson，Tim Hopkins，清单 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

代码输出

导致错误的行是：

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

代码的输出是：

 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.

编译器版本

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

汇编

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

执行

使用了两个不同的 Lapack 库。两者都产生了同样的失败。

/usr/local/lib 
gfortran -g -L/usr/local/lib -llapack -lblas -o mwe m-lapack-interface-sgels.o mwe.o
苹果加速框架 
gfortran -g -framework Accelerate -o mwe m-lapack-interface-sgels.o mwe.o

提示：考虑重定向输出以避免超出屏幕缓冲区：

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

问题

如何修复接口和调用以允许使用allocatable WORK数组？

参考：

连接 FORTRAN lapack 例程时的参数损坏

fortran LAPACK 例程中的分段错误

标签： fortrandynamic-memory-allocationallocationlapack

解决方案

您似乎对 LAPACK 例程的工作方式有些困惑。显示他们的年份，他们从不为您做任何内存管理，您必须在调用例程之前自己分配所有内存。因此，虚拟参数永远不能分配给 LAPACK 例程——我不得不说，我很难相信任何声称上述是正确接口的书。

所以你的界面是错误的。但是您是否对 LAPACK 中的工作区查询如何工作感到困惑？该例程不会为您分配任何内存 - 您调用它一次并将工作空间参数发送到 -1 以指示这是一个工作空间查询，如sgles 文档页面中所述，并且此调用将返回所需的工作空间大小工作中的数组（1）。然后将工作数组分配到正确的大小，最后使用工作区的正确参数再次调用例程。都是有点笨拙，再次显示LAPACK的年龄，但希望上面说清楚

以下是我认为您想要做的正确版本。请注意，我对您的 LAPACK 界面进行了一些其他小的更改，以使其成为我认为完全正确的版本，并且我还删除了许多变量的不必要初始化 - 这是一个不好的习惯，它可以掩盖错误，或者至少使它们更难找到，更糟糕的是通常或多或少与线程编程不兼容，线程编程是我个人认为所有 Fortran 程序员都应该熟悉的一种重要的并行形式。

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

推荐阅读