Bug in zheevd

[jkd2016]

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.

[poulson]

This appears to be an issue with Intel's Fortran compiler aggressively employing FMA's: https://software.intel.com/en-us/forums/intel-c-compiler/topic/579190

Notice that gfortran does not allow parentheses to be ignored by default: see the -fno-protect-parens discussion at https://gcc.gnu.org/onlinedocs/gfortran.pdf

EDIT: I am frightened to realize that this aggressive optimization also effects C and C++.

[jeffhammond]

I cannot reproduce this.

$ ifort -mkl -O2 fma-parens.f90 && ./a.out
 info           0
$ ifort -mkl -O3 fma-parens.f90 && ./a.out
 info           0
$ ifort -mkl -Ofast -xCORE-AVX2 -fma -fp-model fast=2 fma-parens.f90 && ./a.out
 info           0
$ ifort --version
ifort (IFORT) 16.0.2 20160204
Copyright (C) 1985-2016 Intel Corporation.  All rights reserved.

What compiler version are you using?

[jeffhammond]

Oh, I see, the reproducer requires me to build LAPACK from source. Sorry.

At least it is good to know that MKL doesn't have any issue here.

[langou]

Thanks Jay and Jeff. As Kay suggested, I have included the -assume protect_parens in the make.inc.ifort file in the INSTALL directory. Hopefully, that should do the trick. Cheers, Julien.

[jeffhammond]

For what it's worth, this page discusses this in quite a bit of detail, including multiple pages on reassociation.

I attached the latest version of the PDF (FP_Consistency_070816.pdf), which was updated for version 17 of the Intel compilers.

[langou]

Here are some comments from Sven Hammarling.

"Perhaps not surprisingly, I would phrase it the other way round that compiler optimizations should not get in the way of code writers.

It is challenging enough to write quality numerical software without having to worry about what the compiler might do to compromise numerical stability. I would hope that users would prefer correct answers over speed.

That is most certainly not to say that code writers should not strive for efficiency, but that compilers should respect the integrity of their code.

Under "Integrity of parentheses" the Fortran standard explicitly says that any expression in parentheses shall be treated as a data entity. The following example is given:

NOTE 7.16 For example, in evaluating the expression A + (B – C) where A, B, and C are of numeric types, the difference of B and C shall be evaluated before the addition operation is performed; the processor shall not evaluate the mathematically equivalent expression (A + B) – C.

So, I still feel strongly that I should not need to use a switch to ensure that parentheses are not discarded. Parentheses being discarded should only be allowed by a positive action, that is a switch that allows such an action."

[jeffhammond]

I would hope that users would prefer correct answers over speed.

It depends on how you define correct. We are talking about floating point, so accuracy is not a boolean but a threshold.

Scientific computing users want the code to run fast and not give incorrect answers. The vast majority of scientific applications are not rendered incorrect by changes in the last few bits of the significand, particularly since no scientific application I have ever seen is rigorous to all the bits anyways.

It so happens that in this particular case, reassociation causes noticeable incorrectness, but this is exceptional. Thus, it makes sense to cure the exception rather than disable productize optimizations (performance-increasing but not correctness-damaging code transformations).

And you can't ask the compiler to know when reassociation is dangerous, because it is strictly data-dependent, and compilation of Fortran is done offline. In a dynamic language like Python or Julia, one might be able to do something better.

[jeffhammond]

In any case, even if the Intel Fortran compiler turns off productive optimizations everywhere because of a small number of reports where it causes a problem [1], there will be many versions of the compiler that contain this feature, and LAPACK will compile incorrectly with all of them. You have a way to solve this problem, if you choose to use it.

1. If compilers disabled every optimization that broke at least one code, I suspect no compilers would optimize any code ever.

[jeffhammond]

I raised these issues with the Intel compiler team and there are open support issues. Will update as appropriate.

[weslleyspereira]

I couldn't reproduce this issue using ifort version 2021.9.0. Maybe this was solved? I am inclined to close this issue since, as I understand, it is related to compiler optimizations and not to heevd.