Functions returning arrays, derived types and the creation of temporaries.

[Todor_K_]

Dear all,

Consider the following situation followed up by a working test case.

Let there be a derived type wrapping a multidimensional array. Let there further be a function "sum" which adds array components of the objects of a derived type and returns the resulting array. Furthermore, there is a "from_array" assignment subroutine which takes an object and an array, then fills the array component of the object with its array argument.

I am interested how does the compiler interpret an expression such as c = sum(a,b), where = corresponds to "from_array". Is this going to be  optimised into c%arr=a%arr + b%arr, or is an (unnecessary) temporary going to be produced first, e.g: temp=a%arr+b%arr;c%arr=temp?

In the following example, the derived type is "board" (line 3), its objects a,b,c and the statement in question: c=sum(a,b) (line 53). Function "from_array" is also defined as an assignment (lines 7 and 27), and "sum" is defined on line 14.

module board_module
  implicit none
  type board
     real,dimension(:,:), allocatable :: arr
  end type board

  interface assignment(=)
     ! Call from_arr "="
     module procedure from_arr
  end interface assignment(=)


contains
  function sum(x,y)
    ! Add x%arr and y%arr, then return the resulting _array_
    type(board), intent(in) :: x,y
    real,dimension(size(x%arr,1),size(x%arr,2)) :: sum !Is there a
                                                       !more concise
                                                       !way of
                                                       !specifying an
                                                       !automatic
                                                       !multidimensional
                                                       !array object?
    sum = x%arr+y%arr
  end function sum

  subroutine from_arr(dest,src)
    ! Take the src array, and fill dest%arr with it.
    type(board),intent(inout) :: dest
    real,dimension(:,:),intent(in) :: src
    dest%arr = src
  end subroutine from_arr

  subroutine make_board(dest,d1,d2)
    ! Allocate arr
    type(board),intent(out) :: dest
    integer,intent(in) :: d1,d2
    allocate(dest%arr(d1,d2))
  end subroutine make_board
end module board_module

program add
  use board_module
  implicit none
  type(board) :: a,b,c

  call make_board(a,10,20)      !Allocations
  call make_board(b,10,20)
  call make_board(c,10,20)
  a%arr=2                       !Assign some stuff
  b%arr=3

  c = sum(a,b)                  !Is this "temp=a%arr+b%arr;c%arr=temp", or "c%arr=a%arr+b%arr" for the compiler?

  print *, c%arr
end program add

[FortranFan]

In Intel Fortran, you should be able to turn on the run-time option to check for array temporaries being created and investigate this.  But I think for a defined operation that you are trying to use, a temporary array is unavoidable.  Experts on this forum can guide you better.

[FortranFan · Accepted solution]

I forgot to add.  If you rewrite your line 31 as:

[fortran]

   dest%arr(:,:) = src

[/fortran]

then you might be able to avoid the temporary.  You can try it out with the run-time option I suggested.  But then you may lose some resizing benefits with the Fortran array notation.  See Dr Fortran's blog.

[TimP]

Recent versions of ifort have overcome the tendency in some cases to allocate a temporary result, when it is possible to store directly to a declared array.  If you don't want to examine the generated code, the vecanalysis tool

http://software.intel.com/en-us/articles/vecanalysis-python-script-for-annotating-intelr-compiler-vectorization-report

should allow you to infer whether this happened, by searching for memcpy substitutions.

Overloading the name of a Fortran intrinsic seems undesirable, at least from the point of view of making it easier for us to comment.

[Todor_K_]

FortranFan wrote:

I forgot to add.  If you rewrite your line 31 as:

   dest%arr(:,:) = src

then you might be able to avoid the temporary.  You can try it out with the run-time option I suggested.  But then you may lose some resizing benefits with the Fortran array notation.  See Dr Fortran's blog.

That's precisely what I want to avoid. The grand idea behind all of this is to abstract away the types, sizes and dimensions of the arrays. Ideally, I'd like to work with the derived types and let the mechanics bound to them take care of adding, subtracting at all. E.g c=a+b, c=a/b, .... The trouble is allocatable components within the derived types would, in my view, certainly require explicit temporaries. I thought there is some chance that, if I use a simple, contiguous array, this might encourage inlining.

It's not that I want to hide the implementation just for the sake of "writing beautiful code" , the stance I know attracts condescending glares of hardcore number crunching crowd :) . I have a couple of different libraries all working with different intrinsic types , but the algorithm which uses them has the same general structure. As templating does not exist for Fortran, I guess I'd need to resort to some kind of preprocessing If I want to avoid redundancies.

[Todor_K_]

Just to add, I guess something like : c%assignadd(b,c) where "assignadd" is

subroutine assignadd(c,a,b)
type(board), intent(out) :: c
type(board), intent(in) :: a,b

c%arr=a%arr+b%arr
end subroutine assignadd

would do the trick.

It's just that overloading operators is more convenient, because one would not need  to write a separate procedure for any of the more complicated combinations of the array operations. This would also result in a less understandable code, at least, for humans.

What I'd love the most is to do away with any derived types and focus on the array syntax with some means to parametrise dimensions and types of arrays. Fortran Committee, does anybody hear? :-)

[FortranFan]

Todor K. wrote:

....

What I'd love the most is to do away with any derived types and focus on the array syntax with some means to parametrise dimensions and types of arrays. Fortran Committee, does anybody hear? :-)

Have you looked at the parametrized derived type feature in Fortran 2003?  That's the main outstanding aspect from the 2003 standard that Intel Fortran hasn't covered yet.  But I understand Intel wants to gets all Fortran 2003 features implemented in the near future, perhaps even this year.  So one may be able to try it out in an Intel Fortran version in the not too distant future.

Also, have you considered object-oriented methodology with methods i.e., type-bound procedures (TBP)?  With TBP, you can bring in more control of what gets allocated and when and this may help you with your desire to avoid temporary arrays being created.

[Todor_K_]

FortranFan wrote:

Quote:

Todor K. wrote:

....

What I'd love the most is to do away with any derived types and focus on the array syntax with some means to parametrise dimensions and types of arrays. Fortran Committee, does anybody hear? :-)

Have you looked at the parametrized derived type feature in Fortran 2003?  That's the main outstanding aspect from the 2003 standard that Intel Fortran hasn't covered yet.  But I understand Intel wants to gets all Fortran 2003 features implemented in the near future, perhaps even this year.  So one may be able to try it out in an Intel Fortran version in the not too distant future.

Also, have you considered object-oriented methodology with methods i.e., type-bound procedures (TBP)?  With TBP, you can bring in more control of what gets allocated and when and this may help you with your desire to avoid temporary arrays being created.

Yes, I  know about the parameterised derived types. But they offer only "kind" and "len" parametrisation. I was mostly moaning about the parametrisation of  the intrinsic type. Something along the lines (in bastard C++/Fortran pseudocode) of, template<intrinsic> type A; intrinsic comp; end type A . Then at the point of object definition: type(A)<complex> obj_a; type(A)<integer> obj_b, type(A)<real> obj_c, ...

But, I'll admit, this might not be the Fortran way. After all, one can patch up something similar using a preprocessor.

The Fortran way, as I see it, would be to wrap arrays in a derived type like I did with "board", then pass around the board object instead of the naked array, thus eliminating the need to write a separate subroutine with the same code for each complex, real, integer etc. array. Now we come to another point -- we could perform numerical operations on the array components directly, just as you suggested: c%arr=a%arr+b%arr. However, Fortran standard allows a nice way to avoid such ugly code through overriding operators like =, +, - , etc. Moreover, you can define functions that return derived objects (and this is where all the temporaries might rise their ugly heads). Even so, It is much cleaner to write c = a+b, or c = somefun(a,b) + someotherfun(d). Abstracting away the native types is one of the most important usage patterns of types/classes systems in strongly typed languages. For that reason, I think that a compiler should see through this kind of standard code and transparently substitute an expression of the form c=lhs(purefun1(a),purefun2(b),purefun3(c,d),...) by c%arr = lhs(a%arr,b%arr,c%arr,d%arr, ...) . Of course, implementing such optimisations may not  be really high up on the priority list of the compiler designers. Just to be clear, I haven't tested this with the intel compiler yet. Maybe it is smart enough to figure all that out.

All right, this rant is now way off topic for a support forum. My apologies to the Intel crew.

Also, thank you FortranFan for your contribution to my thread!

[jimdempseyatthecove]

It may be likely that with full optimizations on, in particular IPO that the abstract expression/assignment may fully optimize to remove temporaries. The basis I use to make this assumption is many years ago (24) I wrote an optimizer whereby a C++ compiler (either Borland or MS) would be use to generate an assembler file. The assembler file was then parsed by a series of TECO macros to "MUNG" optimize and rework the code. It converted from large and huge model to flat model, and optimized the code in the process. It was noted early on that the optimized code could be optimized with a different set of rules. (TECO was used on PDP-10, PDP-8, PDP11 and MUNG is a recursive acronym for "MUNG Until No Good")

My educated guess is history will repeat itself and this is the same experience as today.

I haven't tried this with the latest compiler, but the earlier versions of IVF with IPO would not generate a .ASM or .s assembler source file with IPO enabled (since this comes in at Link time). What you can to though, is set a break point at a convenient spot such as at a PRINT statement you insert into the code, then open a disassembly window. Reading disassembly is not all that hard to get the jist of things.

Jim Dempsey