Search code examples
fortranlapackintel-mkl

lapack stemr segmentation fault with a particular matrix

I am trying to find the first (smallest) k eigenvalues of a real symmetric tridiagonal matrix using the appropriate lapack routine. I am new to both Fortran and lapack libraries, but (d)stemr seemd to me a good choice so I tried calling it but keep getting segmentation faults.

After some trials I noticed the problem was my input matrix, which has:

  • diagonal = 2 * (1 + order 1e-5 to 1e-3 small variable correction)
  • subdiagonal all equal -1 (if I use e.g. 0.95 instead everything works)

I reduced the code to a single M(not)WE program shown below.

So the question are:

  1. why is stemr failing with such a matrix, while e.g. stev works?
  2. why a segmentation fault?
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

Compilers are invoked via:

  • 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
Solution

  • According to the manual, one issue seems that the argument TRYRAC needs to be a variable (rather than a constant) because it can be overwritten by 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.

    So, for example, a modified code may look like:

    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)