Bug in zheevd #43

There is an issue with the routine zheevd when compiled with Intel fortran using the optimization -O2 or -O3.

The problem occurs at line 285 of dlaed4:

DELTA( J ) = ( D( J )-D( I ) ) - TAU

which evaluates to zero because of the compiler not respecting the parenthesis. This results in a division-by-zero later on. Including -assume protect_parens fixes the problem.

Here is a simple program which produces the error:

program test
implicit none

integer, parameter :: n=27
integer i,j
integer liwork,lrwork,lwork,info
integer, allocatable :: iwork(:)
real(8), allocatable :: w(:),rwork(:)
complex(8), allocatable :: a(:,:),work(:)
liwork=5*n+3
lrwork=2*n**2+5*n+1
lwork=n**2+2*n
allocate(w(n),a(n,n))
allocate(iwork(liwork),rwork(lrwork),work(lwork))

a(:,:)=0.d0
do i=1,n
  a(i,i)=1.d0
  do j=i+1,n
    a(i,j)=cos(dble(i))
  end do
end do

call zheevd('V','U',n,a,n,w,work,lwork,rwork,lrwork,iwork,liwork,info)

print *,'info',info

end program

I'm not sure whether this is a compiler, LAPACK or user issue. I've added -assume protect_parens to our code (elk.sourceforge.net) and made zheev as the default eigenvalue solver for the moment.