Fortran's interface operator behavior on allocatable arrays

Question

I have a derived type that contains a set allocatable arrays and I'm trying to overload some operators. As my arrays can become very large I don't want Fortran to do implicit copy of my arrays, but I can't figure out what is it actually done. Here is a minimal program that illustrate my point:

module type_vector_field
  implicit none


  type vector_field
     double precision, dimension(:), allocatable :: f
  end type vector_field

  interface operator(+)
     module procedure vector_field_add
  end interface


contains


  function vector_field_add(vf1, vf2) result(vf3)
    type(vector_field), intent(in) :: vf1, vf2
    type(vector_field) :: vf3

    print*, "vf* addresses in add:  ", loc(vf1), loc(vf2), loc(vf3)
    vf3%f = vf1%f + vf2%f

  end function vector_field_add


  subroutine vector_field_add2(vf1, vf2, vf3)
    type(vector_field), intent(in) :: vf1, vf2
    type(vector_field), intent(out) :: vf3

    print*, "vf* addresses in add2: ", loc(vf1), loc(vf2), loc(vf3)
    vf3%f = vf1%f + vf2%f

  end subroutine vector_field_add2


end module type_vector_field

The addition of two vector_field can be done either by the + operator or by the vector_field_add2 subroutine. I tested this with the following sample program:

program main
  use type_vector_field
  implicit none

  integer :: n = 6283
  type(vector_field) :: vf1, vf2, vf3

  allocate(vf1%f(n), vf2%f(n))
  vf1%f = 3.1415
  vf2%f = 3.1415

  !! first version

  allocate(vf3%f(n))
  print*, "vf* addresses in main: ", loc(vf1), loc(vf2), loc(vf3)

  vf3 = vf1 + vf2
  print*, "vf* addresses in main: ", loc(vf1), loc(vf2), loc(vf3)
  print*, vf3%f(n)

  deallocate(vf3%f)

  !! second version

  allocate(vf3%f(n))
  print*, "vf* addresses in main: ", loc(vf1), loc(vf2), loc(vf3)

  call vector_field_add2(vf1, vf2, vf3)
  print*, "vf* addresses in main: ", loc(vf1), loc(vf2), loc(vf3)
  print*, vf3%f(n)

end program main

And this is what it outputs when compiled with gfortran-4.8.2:

$ gfortran -g -Wall test.f90 -o test && test
 vf* addresses in main:           6303968          6304032          6304096
 vf* addresses in add:            6303968          6304032  140733663003232
 vf* addresses in main:           6303968          6304032          6304096
   6.2829999923706055     
 vf* addresses in main:           6303968          6304032          6304096
 vf* addresses in add2:           6303968          6304032          6304096
 vf* addresses in main:           6303968          6304032          6304096
   6.2829999923706055

It seems that the subroutine works directly with the vf3 instance defined in the main scope, while the operator create an other vf3 instance in its scope and copy back the result. When I print the arrays addresses instead of the derived type addresses, I get:

 vf*%f addresses in main:          39440240         39490512         39540784
 vf*%f addresses in add:           39440240         39490512                0
 vf*%f addresses in main:          39440240         39490512         39591056
   6.2829999923706055     
 vf*%f addresses in main:          39440240         39490512         39540784
 vf*%f addresses in add2:          39440240         39490512                0
 vf*%f addresses in main:          39440240         39490512         39540784
   6.2829999923706055

Hence the operator seems to copy the address of its vf3%f in the vf3%f instance of the main scope, while the subroutine works directly with the vf3%f of the main scope. In both cases, I don't understand where points the vf3%f of both the subroutine and the operator.

My questions:

Is there a way to know for sure in which array the operator and the subroutine write their results?
Is there a way to avoid copies while using the + operator, or should I drop this idea?

I did not thought of this case. I would expect the value of `vf1` to be added into `vf3` which is a mess from a memory addresses point of vue, and my subroutine `vector_field_add2` crashes. _invalid memory reference_ Are you suggesting me that I cannot implement `+` without creating a temporary array because of your case? — jorispilot, Feb 26 '15 at 14:14
That is precisely it. In particular, for `c=a+b` the result `a+b` is fully evaluated and the result assigned to `c`. I won't write a full answer, as it's perhaps included in my answer to http://stackoverflow.com/q/28241473/3157076. [Though I wouldn't call this a duplicate.] — francescalus, Feb 26 '15 at 14:17

score 1 · Answer 1 · 2015-02-26T15:28:50.750

1

The problem you encounter has, as far as I can see, nothing to do with operator overloading.

More or less by accident (I presume) you had INTENT(OUT) in your subroutine vector_field_add2, and you don't have it in vector_field_add.

One often-overlooked effect of INTENT(OUT) is that all allocatable arrays get deallocated if the dummy arguments have the allocatable attribute upon entry into vector_field_add2. The same goes for types having components with the allocatable attribute.

If you have

   real, dimension(:), allocatable: my_array

   allocate(my_array(some:thing))

and then do

 call sub1(my_array)

with

  sub1(a)
    real, dimension(:), intent(out) :: a

nothing surprising happens.

If you do instead

   call sub2(a)

where sub2 is

   subroutine sub2(a)
   real, dimension(:), allocatable, intent(out) :: a
   end subroutine sub2

the call to sub2 is equivalent to a call to deallocate. If you leave out the intent(out) or replace it, e.g. with intent(inout), nothing happens.

The same thing happened to you, only you had the allocatable component inside your type, so it was even less visible.

Otherwise, there is no difference if you use an operator or if you call a subroutine directly. It only matters for clarity of the code. The compiler translates the assignment statement with the overloaded operator to the subroutine call.

If you want to see how it is done, translate your program using the -fdump-tree-original option. This will generate a file with the main name like your source file, with an added .003.original, where you can see a C-like representation of your Fortran code.

edited Feb 26 '15 at 15:28

answered Feb 26 '15 at 14:29

Having `intent(out)` in `vector_field_add` makes little sense (of course), but it may be worth adding that for the function result variable `vf3` we have that `vf3%f` is initially unallocated, much as in the subroutine? [And, to be clear, this is answering just the part about "vf*%f addresses..."?] – francescalus Feb 26 '15 at 14:51
Effectively, when I leave out the `intent(out)` or replace it by `intent(inout)` I got a non-zero address for `vf3%f` in `vector_field_add2`. But this dosen't help me on whether I should better use an operator or a subroutine. – jorispilot Feb 26 '15 at 15:12
Perhaps I'm confused, but you mention subroutine and "assignment statement". The question is about the defined operation (overloading `+`) rather than defined assignment. The focus must, then, be on the function not the subroutine? That is, "The compiler translates the assignment statement with the overloaded operator to the subroutine call." is unclear to me. Are you saying that the compiler will treat `vf3 = vf1+vf2` as `call vector_field_add2(vf1, vf2, vf3)`? – francescalus Feb 26 '15 at 15:42
Almost, but not not quite (or near enough so it doesn't make a difference) What it does is `vf3 = vector_field_add (&vf1, &vf2);` where vf3 is a struct. There is a slight difference in calling conventions, but not enough to matter. – Feb 26 '15 at 16:41
I was indeed confused. I assumed you were talking about Fortran, rather than C. I now see what you mean. – francescalus Feb 26 '15 at 16:51
As a matter of fact, I am talkling about the intermediate, C-like code that gfortran uses "under the hood". – Feb 26 '15 at 17:24

Fortran's interface operator behavior on allocatable arrays

1 Answers1