CODEWITHSUNDEEP
fortranintel-fortran
JuS
JuS

Reputation: 31

Derived types with allocatable arrays and overload operators in Fortran

In my code, I am using an allocatable derived data type (say, type data) in which I store multidimensional allocatable arrays (x and y). In the same module, I also define routines to allocate/deallocate the whole thing, the assignment operator (=), and additional overload operators (*) and (+). Now, I allocate data1 (of type data), as well as data1%x and data1%y in my main program, initialize them, and perform a simple operation using the overload operators (let's say a simple multiplication of all the elements of data1%x and data1%y by a constant). Here is a minimal code that compiles and reproduces what I have just described:

program minimal

  USE dimensions
  USE typedef

  IMPLICIT NONE

  integer :: i, k
  type(data), dimension(:), allocatable :: data1, data2

  call alloc ( data1 )
  call alloc ( data2 )

  do k = 1 , ndat
    data1(k)%x = real(k)
    data1(k)%y = -real(k)
    data2(k)%x = 0.
    data2(k)%y = 0.
  enddo

  do i = 1, 10
    data2 = data2 + 2.*data1
  enddo

  do k = 1, ndat
    print*, k, maxval(data2(k)%x), maxval(data2(k)%y)
  enddo

  call dealloc ( data1 )
  call dealloc ( data2 )

end program

and the modules:

module dimensions
  integer :: ndat=2
  integer :: m1=10, m2=50
  integer :: n1=10, n2=50
end module dimensions


module typedef

  USE dimensions

  type :: data
    real, dimension(:,:), allocatable :: x
    real, dimension(:,:), allocatable :: y
  end type data

  interface alloc
    module procedure alloc_data
  end interface alloc

  interface dealloc
    module procedure dealloc_data
  end interface dealloc

  interface assignment (=)
    module procedure data_to_data
  end interface

  interface operator (*)
    module procedure const_times_data
  end interface

  interface operator (+)
    module procedure data_plus_data
  end interface

  CONTAINS

  subroutine alloc_data (data1)
    type(data), dimension(:), allocatable, intent(inout) :: data1
    integer :: i

    allocate ( data1(1:ndat) )
    do i = 1, ndat
      allocate ( data1(i)%x(m1:m2,n1:n2) )
      allocate ( data1(i)%y(m1:m2,n1:n2) )
    enddo

  end subroutine alloc_data

  subroutine dealloc_data (data1)
    type(data), dimension(:), allocatable, intent(inout) :: data1
    integer :: i

    do i = 1, ndat
      deallocate ( data1(i)%x )
      deallocate ( data1(i)%y )
    enddo
    deallocate ( data1 )

  end subroutine dealloc_data

  subroutine data_to_data (data2,data1)
    type(data), dimension(:), intent(in) :: data1
    type(data), dimension(1:ndat), intent(out) :: data2
    integer :: i

    do i = 1, ndat
      data2(i)%x = data1(i)%x
      data2(i)%y = data1(i)%y
    enddo

  end subroutine data_to_data

  function const_times_data (c,data1) result(data2)
    type(data), dimension(:), intent(in) :: data1
    real, intent(in) :: c
    type(data), dimension(1:ndat) :: data2
    integer :: i

    do i = 1, ndat
      data2(i)%x = c*data1(i)%x
      data2(i)%y = c*data1(i)%y
    enddo

  end function const_times_data

  function data_plus_data (data1,data2) result(data3)
    type(data), dimension(:), intent(in) :: data1, data2
    type(data), dimension(1:ndat) :: data3
    integer :: i

    do i = 1, ndat
      data3(i)%x = data1(i)%x + data2(i)%x
      data3(i)%y = data1(i)%y + data2(i)%y
    enddo

  end function data_plus_data

end module typedef

Compiling the code with ifort 17.0 (the recommended version on our machine) and the -O0 option for debug does not return any problem. Using optimization levels -O2 or -O3 however produces a segmentation fault. I have tried the same procedure with ifort 18.0 with the same result, whereas ifort 19.0 seems to work.

I have also played a little bit with this minimal code and found out that, for example, it works with optimized ifort 17 if the data structure "data" contains a single element x, or if it is not an allocatable array itself.

The question is very simple: was there a problem in earlier versions of the ifort compiler, or am I doing something wrong? For now, I have found a very simple workaround (which consists in redefining the operator (*) to work on single elements of data, i.e. without any loop in function data_times_data), but I would be interested in knowing a clean way to rewrite the code above to avoid the current issue while taking full advantage of the overload operator functionality.

Many thanks.

Upvotes: 3

Views: 566

Answers (1)

Raul Laasner
Raul Laasner

Reputation: 1585

I can confirm the segfault with ifort 18.0. For some reason the compiler doesn't like the dummy arguments to be arrays when overloading the + or * operators. I suggest keeping the arguments scalar and making the functions elemental instead:

module dimensions
  integer :: ndat=2
  integer :: m1=10, m2=50
  integer :: n1=10, n2=50
end module dimensions


module typedef

  USE dimensions

  type :: data
    real, dimension(:,:), allocatable :: x
    real, dimension(:,:), allocatable :: y
  end type data

  interface alloc
    module procedure alloc_data
  end interface alloc

  interface dealloc
    module procedure dealloc_data
  end interface dealloc

  interface assignment (=)
    module procedure data_to_data
  end interface

  interface operator (*)
    module procedure const_times_data
  end interface

  interface operator (+)
    module procedure data_plus_data
  end interface

  CONTAINS

  subroutine alloc_data (data1)
    type(data), dimension(:), allocatable, intent(inout) :: data1
    integer :: i

    allocate ( data1(1:ndat) )
    do i = 1, ndat
      allocate ( data1(i)%x(m1:m2,n1:n2) )
      allocate ( data1(i)%y(m1:m2,n1:n2) )
    enddo

  end subroutine alloc_data

  subroutine dealloc_data (data1)
    type(data), dimension(:), allocatable, intent(inout) :: data1
    integer :: i

    do i = 1, ndat
      deallocate ( data1(i)%x )
      deallocate ( data1(i)%y )
    enddo
    deallocate ( data1 )

  end subroutine dealloc_data

  elemental subroutine data_to_data (data2,data1)
    type(data), intent(in) :: data1
    type(data), intent(out) :: data2
    integer :: i

    data2%x = data1%x
    data2%y = data1%y

  end subroutine data_to_data

  elemental function const_times_data (c,data1) result(data2)
    type(data), intent(in) :: data1
    real, intent(in) :: c
    type(data) :: data2
    integer :: i

    data2%x = c*data1%x
    data2%y = c*data1%y

  end function const_times_data

  elemental function data_plus_data (data1,data2) result(data3)
    type(data), intent(in) :: data1, data2
    type(data) :: data3
    integer :: i

    data3%x = data1%x + data2%x
    data3%y = data1%y + data2%y

  end function data_plus_data

end module typedef

I think using elemental is better style anyway as opposed to hard-coding the dimensions into the functions, although looking into the Fortran standard I can't immediately find anything that directly prohibits what you were trying to do.

Upvotes: 2

Related Questions