[Bug fortran/91778] New: gfortran GCC9 optimizer bug

[mark.wieczorek at oca dot eu]

https://gcc.gnu.org/bugzilla/show_bug.cgi?id=91778

            Bug ID: 91778
           Summary: gfortran GCC9 optimizer bug
           Product: gcc
           Version: 9.2.0
            Status: UNCONFIRMED
          Severity: normal
          Priority: P3
         Component: fortran
          Assignee: unassigned at gcc dot gnu.org
          Reporter: mark.wieczorek at oca dot eu
  Target Milestone: ---

I am writing about a possible bug in the gfortran GCC9 optimizer on macOS
(installed via brew).

Before going into the details, I note that my code (SHTOOLS/pyshtools) is
widely used on many platforms and compilers. My code works with GCC8 compiled
with optimizations "-O" or "-O3", and it works fine with GCC9 when compiled
_without_ optimizations. I was able to "fix" my code to work with GCC9, but I
feel that what I am doing is avoiding a bug in the GCC9 optimizer, and that I
am not in fact "fixing" my code (perhaps I am wrong...).

The problem is related to using the FFTW3 library, which is the most widely
used FFT library for scientific computing. If this is a bug, then others will
probably encounter similar problems. As my code is somewhat long (and given the
lack of time I have now), I will just give you a summary of two problems. If
necessary, I could try to write a "small" example that reproduces these
problems when I have more free time later.

I start by describing how FFTW routines are use. First, you initialize the FFT
operation and get pointers to all the input and output arrays, which are stored
in the variable "plan":

    call dfftw_plan_dft_c2r_1d(plan, nlong, coef, grid)

Then you perform the FFT simply by calling

    call dfftw_execute(plan)

The first problem boils down to this:

    call dfftw_plan_dft_c2r_1d(plan, nlong, coef, grid)

    coef(1) = dcmplx(coef0,0.0d0) ! A
    coef(2:lmax_comp+1) = coef(2:lmax_comp+1) / 2.0d0

    call dfftw_execute(plan) ! AA
    gridglq(i,1:nlong) = grid(1:nlong)

    coef(1) = dcmplx(coef0s,0.0d0) ! B
    coef(2:lmax_comp+1) = coefs(2:lmax_comp+1)/2.0d0

    call dfftw_execute(plan) ! BB
    gridglq(i_s,1:nlong) = grid(1:nlong)


The problem is that the optimizer thinks the line A is redundant with line B
(the same variable is being defined twice). Thus, the optimizer sets line A to
that of line B and deletes line B. I have verified this by doing so in my code.
However, line A is necessary to execute line AA, and line B is necessary to
execute line BB. The optimizer probably doesn't realize this because the
variable "coef" is not explicitly included when calling the function
dfftw_execute(plan).

The second problem I encountered is a little more mysterious. These are the
_last_ 4 lines of the subroutine:

    coef(lmax_comp+1) = coef(lmax_comp+1) + cilm(1,lmax_comp+1,lmax_comp+1)
    coef(nlong-(lmax_comp-1)) = coef(nlong-(lmax_comp-1)) &
                                + cilm(2,lmax_comp+1,lmax_comp+1)

    call dfftw_execute(plan)

    griddh(i_eq,1:nlong) = grid(1:nlong)

The problem is that the optimizer ignores the first two lines. The reason for
this is probably because (1) the variable coef is not explicitly noted in the
fftw call, and (2) the variable coef is not output in the subroutine. Thus, the
optimizer probably thinks that it doesn't need to compute the first two lines 

So, in summary, I believe that the GCC9 optimizer is not working correctly
because it doesn't realize that the function call dfftw_execute(plan) actually
depends on the variables coef and grid. Given that my code has worked well with
all other versions of GCC, I suspect that there has been a change in how the
optimizer works.