Hi Mat,
I was investigating a workaround for this ICE when I found some further misbehavior of the numerical results. Consider this test case
program test
implicit none
integer, dimension(3) :: memdims
integer, dimension(3) :: ksdims
real(8), dimension(:), allocatable :: data
real(8), dimension(:,:,:), allocatable :: data_reshaped
integer :: Ngrid, Npoints
integer :: i
memdims = (/1, 2, 4/)
ksdims = (/4, 2, 3/)
Npoints = product (memdims)
Ngrid = product (ksdims)
allocate (data(Ngrid))
allocate (data_reshaped(memdims(3), memdims(2), memdims(1)))
do i = 1, Ngrid
data(i) = real (i)
end do
i = 3
data_reshaped = reshape(data(1:Npoints), memdims(3:1:-1), order=(/3,2,1/))
!!! In the original code, this is part of a loop where i appears as index
data (1:Npoints * ksdims(i)) = pack(spread(data_reshaped, i, ksdims(i)), .true.) !!! BROKEN
!!! If i is replaced by a constant 3, the gfortran results are reproduced
!!!data (1:Npoints * ksdims(3)) = pack(spread(data_reshaped, 3, ksdims(3)), .true.) !!! OK
!!! Put i on different positions. Apparently, the "dim" argument of spread causes trouble
!!!data (1:Npoints * ksdims(3)) = pack(spread(data_reshaped, 3, ksdims(i)), .true.) !!! OK
!!!data (1:Npoints * ksdims(3)) = pack(spread(data_reshaped, i, ksdims(3)), .true.) !!! BROKEN
!!!data (1:Npoints * ksdims(i)) = pack(spread(data_reshaped, 3, ksdims(3)), .true.) !!! OK
do i = 1, Ngrid
print *, data(i)
end do
deallocate (data)
deallocate (data_reshaped)
end program test
I have compiled this both with nvfortran and gfortran. The output of the nvfortran binary contains zeroes at the positions where the array is supposed to be extended. These zeroes are not generated by the gfortran binary. Interestingly, this does not happen when i
is replaced by a numerical constant.
Maybe this example is helpful in fixing the problems with the array operations.
Regards,
Christian