Fortran Allocatable Array Member of a User-Defined Type

Question

I am stuck with segmentation fault at an allocatable array memberof a derived type in the following simple program. This segmentation fault occurs only on one machine (with Intel Fortran 14.0.3 on openSUSE) but not on the other machine (with Intel Fortran 14.0.2 on Ubuntu) that I tried. Also, if I change one of the integer parameters in the program, the program ends normally.

Could anybody reproduce the prolem? Could anybody tell me what is wrong with the code?

Below are the three source code files.

main_dbg.f90 .. whether the segmentation fault occurs or not depends on the values of n1 and n2 in this file.

PROGRAM dbg
  USE tktype
  USE mymodule, ONLY : MyClass, MyClass_constructor
  IMPLICIT NONE

  INTEGER(I4B)                :: n1,n2,n3
  TYPE(MyClass)               :: o_MyClass

  n1=23
  n2=32
  ! .. this does not work.
  ! n2=31 
  ! .. this works.
  n3 = n1*n2
  write(*,'(1X,A,I10)') 'n1=', n1
  write(*,'(1X,A,I10)') 'n2=', n2
  write(*,'(1X,A,I10)') 'n3=', n3

  o_MyClass = MyClass_constructor(n1, n2, n3) 

  call o_MyClass%destructor()
  write(*,*) '***************************'
  write(*,*) '   Normal End :)           '
  write(*,*) '***************************'

END PROGRAM dbg

strange.f90 .. segmentation fault occurs at the forall construct in this file.

!*******************************************************************
MODULE mymodule
!*******************************************************************
  USE tktype
  IMPLICIT NONE
  PRIVATE

  PUBLIC MyClass
  PUBLIC MyClass_constructor

  TYPE :: MyClass
     PRIVATE
     REAL(DP),     DIMENSION(:),     ALLOCATABLE :: arrA
     COMPLEX(DPC), DIMENSION(:,:,:), ALLOCATABLE :: arrB
   CONTAINS
     PROCEDURE :: destructor
  END TYPE MyClass

! ================================================================
CONTAINS
! ================================================================

  ! ****************************************************************
  FUNCTION MyClass_constructor(n1, n2, n3) RESULT(this)
  ! ****************************************************************
    TYPE(MyClass)                :: this
    INTEGER(I4B),    INTENT(IN)  :: n1, n2, n3
    ! local variables
    INTEGER(I4B) :: j1, j2, j3

    write(*,'(1X,A)') 'entered constructor..'

    allocate(this%arrA(n2))
    allocate(this%arrB(n1, n2, n3))

    this%arrA = 1.0_dp

    write(*,*) 'size(this%arrB,1) =', size(this%arrB,1)
    write(*,*) 'n1                = ', n1
    write(*,*) 'size(this%arrB,2) =', size(this%arrB,2)
    write(*,*) 'n2                = ', n2
    write(*,*) 'size(this%arrB,3) =', size(this%arrB,3)
    write(*,*) 'n3                = ', n3

    forall(j1=1:n1, j2=1:n2, j3=1:n3)
       this%arrB(j1,j2,j3)  = this%arrA(j2) 
    end forall

    write(*,'(1X,A)') '..leaving constructor'

  END FUNCTION MyClass_constructor


  ! ****************************************************************
  SUBROUTINE destructor(this)
  ! ****************************************************************
    CLASS(MyClass),             INTENT(INOUT) :: this

    deallocate(this%arrA)
    deallocate(this%arrB)

  END SUBROUTINE destructor

END MODULE mymodule

tktype.f90

! ********************************************************************
MODULE tktype
! ********************************************************************
!   module tktype is an extraction of module nrtype in Numerical Recipes in 
!   Fortran 90.
! ********************************************************************
  !   Symbolic names for kind types of 4-, 2-, and 1-byte integers:
  INTEGER, PARAMETER :: I4B = SELECTED_INT_KIND(9)
  INTEGER, PARAMETER :: I2B = SELECTED_INT_KIND(4)
  INTEGER, PARAMETER :: I1B = SELECTED_INT_KIND(2)
  !   Symbolic names for kind types of single- and double-precision reals:
  INTEGER, PARAMETER :: SP = KIND(1.0)
  INTEGER, PARAMETER :: DP = KIND(1.0D0)
  !   Symbolic names for kind types of single- and double-precision complex:
  INTEGER, PARAMETER :: SPC = KIND((1.0,1.0))
  INTEGER, PARAMETER :: DPC = KIND((1.0D0,1.0D0))
  !   Symbolic name for kind type of default logical:
  INTEGER, PARAMETER :: LGT = KIND(.true.)
END MODULE tktype

Below is a shell script to compile the source codes above and run the generated executable.

compile_run.sh

#!/bin/bash

ifort -v 
echo "compiling.."
ifort -o tktype.o -c -check -g -stand f03 tktype.f90
ifort -o strange.o -c -check -g -stand f03 strange.f90
ifort -o main_dbg.o -c -check -g -stand f03 main_dbg.f90
ifort -o baabaa strange.o tktype.o main_dbg.o
echo "..done"
echo "running.."
./baabaa
echo "..done"

The standard output looked as following.

ifort version 14.0.3
compiling..
..done
running..
 n1=        23
 n2=        32
 n3=       736
 entered constructor..
 size(this%arrB,1) =          23
 n1                =           23
 size(this%arrB,2) =          32
 n2                =           32
 size(this%arrB,3) =         736
 n3                =          736
./compile_run.sh: line 11: 17096 Segmentation fault      ./baabaa
..done

Edit 2016-01-30

I found that adding ulimit -s unlimited at the beginning (after #/bin/bash) of compile_run.sh prevents the segmentation fault. Are the allocatable arrays in fortran stored in stack, not in heap?

Answer 1

This may be a possible duplicate of the similar question (Segmentation fault on 2D array), where some multi-dimensional forall loop causes a problem. The OP of the linked question asked this in the Intel forum (ifort v 14.0 / 15.0 "-g" option causes segFault), and the latest reply is like the following:

Workaround #1 is to increase stack size limit. I was successful with your test case using: ulimit -s unlimited

Workaround #2 is to use DO loops instead of FORALL, as follows:

Also, according to casey's comment in the linked question, this problem does not occur for ifort16, so I guess it may be a compiler issue specific for ifort14/15.

---

More info (just some experiment):

The same problem was reproduced on my computer by limiting the stack size to ulimit -s 4000 and using ifort14.0.1, and it disappeared with the -heap-arrays option. So I initially thought that there may be some automatic arrays or array temporaries of size n1 * n2 * n3, but there seems no such thing in the original code... Attaching -assume realloc_lhs or -check -warn did not help either.

So I have made a test program that performs the same calculation using do or forall:

program main
    implicit none
    integer, parameter :: dp  = KIND(1.0D0)
    integer, parameter :: dpc = KIND((1.0D0,1.0D0))
    type Mytype
        real(dp),     allocatable :: A(:)
        complex(dpc), allocatable :: B(:,:,:)
    endtype
    type(Mytype) :: t
    integer :: n1, n2, n3, j1, j2, j3

    n1 = 23
    n2 = 32
    n3 = n1 * n2   !! = 736

    allocate( t% A( n2 ), t% B( n1, n2, n3 ) )

    t% A(:) = 1.0_dp

    print *, "[1] do (3-dim)"
    do j3 = 1, n3
    do j2 = 1, n2
    do j1 = 1, n1
        t% B( j1, j2, j3 ) = t% A( j2 ) 
    enddo
    enddo
    enddo

    print *, "[2] do (1-dim)"
    do j2 = 1, n2
        t% B( :, j2, : ) = t% A( j2 ) 
    enddo

    print *, "[3] forall (1-dim)"
    forall( j2 = 1:n2 )
        t% B( :, j2, : ) = t% A( j2 ) 
    end forall

    print *, "[4] forall (3-dim)"   ! <-- taken from the original code
    forall( j1 = 1:n1, j2 = 1:n2, j3 = 1:n3 )
        t% B( j1, j2, j3 ) = t% A( j2 )
    end forall

    print *, "all passed."
end program

where pattern [4] corresponds to that used by OP. Limiting the stack size and compiling with no option (ulimit -s 4000 ; ifort test.f90) gives the output

 [1] do (3-dim)
 [2] do (1-dim)
 [3] forall (1-dim)
 [4] forall (3-dim)
Segmentation fault

which means that only the pattern [4] fails when -heap-arrays is not attached. Strangely, the problem disappears when arrays A and B are declared outside the derived type, i.e., the following program works with no options.

program main
    implicit none
    integer, parameter :: dp  = KIND(1.0D0)
    integer, parameter :: dpc = KIND((1.0D0,1.0D0))
    real(dp),     allocatable :: A(:)
    complex(dpc), allocatable :: B(:,:,:)
    integer :: n1, n2, n3, j1, j2, j3

    n1 = 23
    n2 = 32
    n3 = n1 * n2   !! = 736

    allocate( A( n2 ), B( n1, n2, n3 ) )

    A(:) = 1.0_dp

    print *, "[1] do (3-dim)"
    do j3 = 1, n3
    do j2 = 1, n2
    do j1 = 1, n1
        B( j1, j2, j3 ) = A( j2 ) 
    enddo
    enddo
    enddo

    print *, "[2] do (1-dim)"
    do j2 = 1, n2
        B( :, j2, : ) = A( j2 ) 
    enddo

    print *, "[3] forall (1-dim)"
    forall( j2 = 1:n2 )
        B( :, j2, : ) = A( j2 ) 
    end forall

    print *, "[4] forall (3-dim)"
    forall( j1 = 1:n1, j2 = 1:n2, j3 = 1:n3 )
        B( j1, j2, j3 ) = A( j2 )
    end forall

    print *, "all passed."
end program

So it seems that the problem occurs only in some particular case of multi-dimensional forall loop (even without -g option), which might be using an internal temporary array on the stack (though -check -warn option gives no message). FYI, all the above patterns work with gfortran 4.8/5.2 and Oracle fortran 12.4.