Aliasing issue with zero-length array

Indeed. Many codes on Netlib share this problem and inspection of the documentation (if any) or the procedure statements to infer the expected subscript range is needed to judge what is (L)ength supposed to be.

I think the typical example is in the ODE integrators, where you can provide either a scalar absolute tolerance, or a vector. Since these are typically intent(in) you can usually fix this today with an array expression call odeslv(...,[atol],iopt). An example is shown here: Has anyone used DVODE and is it good? - #27 by CRquantum


This code was published in 1991, but probably it originated a decade earlier, around 1980 or so. In the accompanying articles, the authors mention running it on an Amdahl 5850. I found a digital copy of an Amdahl FORTRAN 77 Reference Manual, and in Section 3.4 it describes assumed-size semantics. In Section 10.7 it also says,

If an actual argument is an array element name, that element is associated with the initial dummy array element. In this case, the size declared for the dummy array argument must satisfy the formula size <= (ds+1-as)*(e1) where size is the storage required by the dummy array argument, ds is the maximum subscript value that an actual array element can attain, as is the specified subscript value for the actual array element, and e1 is the storage size occupied by one actual array element.

The array in the subroutine must not attempt to access a subscript value that beyond the limit of the array in the invoking program (see chapter 3, “Dummy Array and Actual Array” under paragraph 3.4). […]

According to these rules, the actual argument V(IV) would amount to a 1-element array. The authors of the code knew this, because in a manual, they say

Note that in the case when there are no ODEs coupled to the PDEs, V and VDOT will be dummy vectors of length one.

but they stop short of saying that the element V(1) should not be referenced when there are no ODEs (NV = 0).

Does that mean with Fortran 90, I can correct the code like shown below?

      IV = NPDE*NPTS
      CALL CSET(NPDE,NPTS,U,WK(I6),WK,WK(I2),WK(I5),NEL,NPTL,WK(I4),
     *          WK(I12),XBK,IBK,WK(I3),U(IV+1:IV+NV),NV)

This way when NV=0, the slice U(IV+1:IV+NV) would be a zero-length array?

EDIT: running a quick check on the program below,

      subroutine zerolen(nv,v)
      integer nv
      double precision v(nv)
      print *, nv
      v(1) = 3.0
      end

C main program
      external zerolen
      double precision v(6)

C Fortran 90 Solution
      nv = 0
      call zerolen(nv,v(5+1:5+nv))
      end

compiled with -fcheck=bounds shows the output I wanted to achieve:

At line 5 of file /app/example.f90
Fortran runtime error: Index '1' of dimension 1 of array 'v' above upper bound of 0