具有特定矩阵的 lapack stemr 分段错误

Question

我正在尝试使用适当的 lapack 例程找到实数对称三对角矩阵的第一个（最小）k 个特征值。我是 Fortran 和 lapack 库的新手，但 (d)stemr 在我看来是个不错的选择，所以我尝试调用它，但一直出现分段错误。

经过一些试验，我发现问题出在我的输入矩阵上，它有：

对角线=2*（1+阶1e-5到1e-3小变量修正）
次对角线都等于 -1（如果我使用例如 0.95 而不是一切正常）

我将代码缩减为单个 M(not)WE 程序，如下所示。

所以问题是：

为什么 stemmr 使用这样的矩阵会失败，而例如史蒂夫在工作吗？
为什么会出现分段错误？

program mwe

    implicit none

    integer, parameter :: n = 10
    integer, parameter :: iu = 3
    integer :: k
    double precision :: d(n), e(n)
    double precision :: vals(n), vecs(n,iu)

    integer :: m, ldz, nzc, lwk, liwk, info
    integer, allocatable :: isuppz(:), iwk(:)
    double precision, allocatable :: wk(:)

    do k = 1, n
        d(k) = +2d0 + ((k-5.5d0)*1d-2)**2
        e(k) = -1d0 ! e(n) not really needed
    end do

    ldz = n
    nzc = iu
    allocate(wk(1), iwk(1), isuppz(2*iu))

    ! ifort -mkl gives SIGSEGV at this call  <----------------
    call dstemr( &                              
        'V', 'I', n, d, e, 0d0, 0d0, 1, iu, &
        m, vals, vecs, ldz, -1, isuppz, .true., &
        wk, -1, iwk, -1, info)
    lwk = ceiling(wk(1)); deallocate(wk); allocate(wk(lwk))
    liwk = iwk(1); deallocate(iwk); allocate(iwk(liwk))
    print *, info, lwk, liwk ! ok with gfortran

    ! gfortran -llapack gives SIGSEGV at this call  <---------
    call dstemr( &
        'V', 'I', n, d, e, 0d0, 0d0, 1, iu, &
        m, vals, vecs, ldz, nzc, isuppz, .true., &
        wk, lwk, iwk, liwk, info)

end program

编译器通过以下方式调用：

gfortran [(GCC) 9.2.0]: gfortran -llapack -o o.x mwe.f90
ifort [(IFORT) 19.0.5.281 20190815]: ifort -mkl -o o.x mwe.f90

Answer 1

根据 manual，一个问题似乎是参数 TRYRAC 需要是一个变量（而不是常量），因为它可以被 dstemr() 覆盖：

[in,out] TRYRAC : ... On exit, a .TRUE. TRYRAC will be set to .FALSE. if the matrix does not define its eigenvalues to high relative accuracy.

因此，例如，修改后的代码可能如下所示：

logical :: tryrac
...
tryrac = .true.
call dstemr( &                              
    'V', 'I', n, d, e, 0d0, 0d0, 1, iu, &
    m, vals, vecs, ldz, -1, isuppz, tryrac, &  !<--
    wk, -1, iwk, -1, info)
...
tryrac = .true.
call dstemr( &
    'V', 'I', n, d, e, 0d0, 0d0, 1, iu, &
    m, vals, vecs, ldz, nzc, isuppz, tryrac, &  !<--
    wk, lwk, iwk, liwk, info)

具有特定矩阵的 lapack stemr 分段错误

lapack stemr segmentation fault with a particular matrix

fortran

lapack

intel-mkl