Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

Automatic array allocation upon assignment in Fortran

We recently discovered we were making an assignment to an unallocated array in Fortran. The GNU gfortran compiler didn't catch the error, and the code runs under both OSX and Linux. However, the same code segmentation faults on a IBM Power PC.

My question is, is the following code correct? It seems that the array assigned to array is automatically allocating memory, on some architectures, but not on others. Are there implementation specific details at work here?

The code is a mixed C/Fortran code :

#include <stdlib.h>

void assign_array_(double x[], int* n);
void print_array_();

int main()
{
    int n,i;
    double *x;

    n = 5;
    x = (double*) malloc(sizeof(double)*n);

    for (i = 0; i < n; i++)
        x[i] = (double) i;

    assign_array_(x,&n);
    print_array_();

    return 0;
}

And the Fortran code :

MODULE test_mod
  DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:) :: array
  integer :: nsize
END MODULE test_mod

SUBROUTINE assign_array(x,n)
  USE test_mod
  IMPLICIT NONE

  INTEGER :: n
  DOUBLE PRECISION :: x(n)

  CALL test_allocated()
  array = x
  CALL test_allocated()
  nsize = n

END SUBROUTINE assign_array


SUBROUTINE print_array()
  USE test_mod, ONLY: nsize, array
  IMPLICIT NONE

  INTEGER :: i

  DO i = 1,nsize
     WRITE(6,'(F24.16)') array(i)
  END DO

END SUBROUTINE print_array

SUBROUTINE test_allocated()
  USE test_mod
  IMPLICIT NONE

  IF (ALLOCATED(array)) THEN
     WRITE(6,*) 'Array is allocated'
     WRITE(6,*) 'size is ', SIZE(array)
  ELSE
     WRITE(6,*) 'Array is NOT allocated'
  END IF
END SUBROUTINE test_allocated

The output (when it runs) is :

Array is NOT allocated
Array is allocated
size is            5
  0.0000000000000000
  1.0000000000000000
  2.0000000000000000
  3.0000000000000000
  4.0000000000000000

And here is the output on the Power PC :

Array is NOT allocated
Segmentation fault (core dumped)

In summary: It runs when compiled under GNU (GNU Fortran (MacPorts gcc5 5.4.0_0) 5.4.0) gfortran on OSX (arch : x86_64h) and Linux (in a virtual machine hosted on OSX, GNU Fortran (Ubuntu 4.9.4-2ubuntu1~14.04.1) 4.9.4), but fails to run when compiled on the Power PC (arch: ppc64) compiled using GNU Fortran (GCC) 4.4.7 20120313 (Red Hat 4.4.7-17). In our original code, the Power PC implementation only segfaulted much later in the code, where entries assigned array was referenced, making our 'bug' (if it is in fact a bug) really hard to track down.

What is the correct behavior for the above code?

like image 480
Donna Avatar asked Feb 09 '17 15:02

Donna


1 Answers

The validity of code like

integer, allocatable :: array(:)
array = (/1,2,3/)
end

depends on the standard of Fortran used to interpret it.

Fortran 2003 introduced the concept of automatic allocation on intrinsic assignment. Before Fortran 2003 the array on the left-hand side of such an assignment statement must be allocated, and of the same shape as the array on the right-hand side.

From Fortran 2003 only the rank needs match. If there is a shape mismatch, the array would be first deallocated and then reallocated to the correct shape. If not allocated initially, it would be allocated.

So, the program above isn't valid Fortran 90, but is valid Fortran 2003.

The difference in the real-world code, then, is from what language syntax the compilers support.

For gfortran, Fortran 2003's assignment to an allocatable array was introduced in 4.6, 2011-01-28.

As also commented, the command line option -fno-realloc-lhs1 disables this automatic (re-)allocation, making the compiler not Fortran 2003+ compliant.


1 Other compilers have similar behaviour: adding a required check for whether reallocation is necessary is a performance hit which is redundant in Fortran 90 compliant code and may be a feature not used by many even in modern code. For Intel's compiler, for example, in some versions which do support F2003 the default is to ignore this.

One can always suppress the (re-)allocation checks/actions of an array in modern code by using an array section

array(:) = (/1,2,3/)

In this case, array (if allocatable) must be allocated, of rank 1 and of size 3 for the assignment statement to be valid. This is as would be under a Fortran 90 interpretation of the assignment with the whole array array=(/1,2,3/).

The reason for this is that with the array section of this footnote, the left-hand side is not allocatable, even though the array itself is.

like image 170
francescalus Avatar answered Sep 30 '22 20:09

francescalus