0

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
rgbmrc
  • 3
  • 2
  • a) This isn't a complete program, could you please provide one b) providing the exact compile and link lines as well would be very useful as well c) real(8) is not portable – Ian Bush Nov 15 '19 at 17:01
  • I updated the code with the missing program/end program and added the compiling commands. I also changed (8) to real*8, didn't know about the portability issues, thank you for letting me know. – rgbmrc Nov 15 '19 at 17:20
  • real*8 is much, MUCH WORSE - it is not portable and not standard! See https://stackoverflow.com/questions/838310/fortran-90-kind-parameter. – Ian Bush Nov 15 '19 at 17:25

1 Answers1

3

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)
roygvib
  • 7,218
  • 2
  • 19
  • 36
  • I hope there are some explicit interface header or something for this kind of routines... (lapack95 etc might include it but not very sure) – roygvib Nov 15 '19 at 18:56
  • I'm not sure too. Intel mkl documentation lists lapack95 for many routines but not for this one. Anyway this is not a problem for me as this is part of an uni assignment that I have to submit, I prefer to use to as few dependencies as possible. Your answer solved my problem, I read the [in,out] in the documentation but then forgot it, thank you very much for pointing it out! – rgbmrc Nov 15 '19 at 19:19