Silverfrost Forums

Welcome to our forums

Code from source with SPREAD causes memory corruption

18 Dec 2018 6:23 (Edited: 18 Dec 2018 6:59) #23028

The test code below arose from an investigation of bugs in the venerable Eispack package at Netlib (test CHTEST). The old code has some undefined variables, and some subroutine calls with aliased arguments. FTN95 helped with the undefined variables, but not with the aliased arguments.

Refactored code from the investigation showed up an FTN95 bug. The bug is present in the 32 and 64 bit Version 8.30.290 compiler, but not in the 7.20 compiler. To keep the test code short, I have placed the data in an unformatted file (696 bytes), which can be downloaded from Dropbox: https://www.dropbox.com/s/i0loyu5ip2ya02q/sbug.sbn?dl=0 .

The code:

      program sbug
      implicit none
      integer, parameter :: double = kind(0.0d0)
      real(double) :: ar(5,5),ai(5,5),tau(2,5),zr(5,5),zi(5,5)
      integer :: n
      open(20,file='sbug.sbn',form='unformatted',status='old')
      read(20)n
      if(n /= 5)stop 'Incorrect n for data file'
      read(20)ar,ai,tau,zr
      close(20)
      call htribk (n, ar, ai, tau, zr, zi)
      end

      subroutine htribk(n, ar, ai, tau, zr, zi)
      implicit none
      integer, parameter :: double = kind(0.0d0)
      integer , intent(in) :: n
      real(double) , intent(in) :: ar(n,n), ai(n,n), tau(2,n)
      real(double) , intent(in out) :: zr(n,n)
      real(double) , intent(out)    :: zi(n,n)
      integer :: i, j, k
      real(double) :: h
      real(kind(0.0d0)) :: s1u(n), si1u(n)
!
      zi(:n,:) = -zr(:n,:)*spread(tau(2,:),dim = 2,ncopies = n)
      zr(:n,:) =  zr(:n,:)*spread(tau(1,:),dim = 2,ncopies = n)
      i = 2
      k = i - 1
      h = ai(i,i)
      s1u = sum(spread(ar(i,:k),dim = 2,ncopies = n)*zr(:k,:)- &
                spread(ai(i,:k),dim = 2,ncopies = n)*zi(:k,:), dim=1)
      si1u = sum(spread(ar(i,:k),dim = 2,ncopies = n)*zi(:k,:)+&
                spread(ai(i,:k),dim = 2,ncopies = n)*zr(:k,:), dim=1)
      s1u  = (s1u/h)/h
      si1u = (si1u/h)/h
      PRINT *,' i      S1u       Si1u:'
      print '(1x,i2,2ES13.4)',(j,s1u(j),si1u(j),j=1,n)
      stop
      end subroutine htribk

The correct output results, which FTN95 7.20 gives (and also Gfortran, etc., provided a suitable unformatted data file is used):

  i      S1u       Si1u:
  1   2.3405E-01   0.0000E+00
  2   8.5779E-02   0.0000E+00
  3   1.3610E-01   0.0000E+00
  4   1.0824E-01   0.0000E+00
  5   8.4260E-04   0.0000E+00

The results from FTN95 8.30.279:

32-bit

Default, /lgo
  i      S1u       Si1u:
  1   0.0000E+00   0.0000E+00
  2   9.8118-294   9.8089-294
  3   0.0000E+00   0.0000E+00
  4   9.8118-294   9.8089-294
  5   0.0000E+00   0.0000E+00

/checkmate /lgo
     Invalid floating point operation

/checkmate
  1   9.8040-294   9.8033-294
  2  -4.9760+234  -4.9760+234
  3   0.0000E+00   0.0000E+00
  4   9.8107-294   9.8078-294
  5   9.8040-294   9.8033-294

/opt
  i      S1u       Si1u:
  1   2.3405E-01   0.0000E+00
  2   8.5779E-02   0.0000E+00
  3   1.3610E-01   0.0000E+00
  4   1.0824E-01   0.0000E+00
  5   8.4260E-04   0.0000E+00      (all correct!)

[Message limit reached, to be continued in next post in thread]

18 Dec 2018 6:25 #23029

...CONTINUATION...

64-bit results

Default, /lgo
  i      S1u       Si1u:
  1   2.3405E-01   0.0000E+00      (correct)
  2   8.5779E-02   0.0000E+00         '
  3   1.3610E-01   0.0000E+00         '
  4   1.0824E-01   0.0000E+00         '
  5   2.3405E-01   0.0000E+00      (incorrect)

/checkmate
  same as for Default

/opt
  i      S1u       Si1u:
  1   2.3405E-01   0.0000E+00
  2   8.5779E-02   0.0000E+00
  3   1.3610E-01   0.0000E+00
  4   1.0824E-01   0.0000E+00
  5   8.4260E-04   0.0000E+00      (all correct)

NOTE: I rarely use the SPREAD intrinsic when writing code myself, and I noted that replacing the four statements containing SPREAD with the corresponding DO loops or other array expressions gives correctly functioning code. The SPREAD calls were introduced by an automatic Fortran clean-up utility; the original code was in Fortran 66.

18 Dec 2018 7:13 #23030

It's a glitch in The Matrix, surely, when one gets that creepy sense of Déjà vu about a post in The Forum.

18 Dec 2018 7:21 #23031

If you are interested in knowing how this got started, see the following thread by Jack Wolfman on Comp. Lang. Fortran:

 https://groups.google.com/forum/#!topic/comp.lang.fortran/trzQK_bMipY

In the past, I had filed a couple of bug reports involving SPREAD. One of them involved the compiler rejecting valid code, and another a performance hit as a result of using SPREAD. I do not remember having seen wrong results caused by its use.

What is quite a surprise with the current bug report is that one gets correct results (mostly) with /OPT and really bad results with /CHECKMATE.

19 Dec 2018 8:22 #23035

mecej4

Thank you for the feedback.

There are similarities to earlier outstanding bug reports and hopefully when they are fixed this one will also come out in the wash.

4 Apr 2019 6:50 #23429

This has now been fixed for the next release of FTN95.

Please login to reply.