What is the best practice when a Fortran function returns an array?

Question

Say I want to write a function that takes as input an array x of one dimension and returns another array y of the same dimension based on it (to ilustrate it I'm using a function that multiplies it by 2). I have two options for that code:

function times2(x) result(y)
    real, intent(in) :: x(:)
    real, allocatable :: y(:)

    allocate(y(size(x))
    y = 2*x
end function

or

function times2(x,n) result(y)
    real, intent(in) :: x(n)
    integer, intent(in) :: n
    real  :: y(n)

    y = 2*x
end function

Personally I prefer the first one because is easier to use for the caller, but I'm not sure which is better in terms of memory, suppose that the array x can be huge, I don't know if it is better to be a deferred array or an automatic array. In any case, which is the good way to do it in modern Fortran?

Answer 1

Probably neither, though as usual with these things, the answer depends on specific circumstances.

Assuming a non-elemental operation, I would tend to write such a function (in a module) as:

function times2(x) result(y)
  real, intent(in) :: x(:)
  real :: y(size(x))

  y = 2*x
end function

The above has an assumed shape dummy argument, with automatic function result. It:

• is available when writing to the Fortran 95 standard;

• specifies explicitly in source the dependence of the size of the function result on the function arguments, which may (or may not) help readers of your code understand what is going on (one such reader is the compiler itself, which might help it with optimisation);

• may (or may not) avoid intermediate copies of the value of the array;

• may (or may not) require space for the function result or equivalent temporary on the stack.

If the operation was elemental (i.e. the same operation on each element, as per the actual example given), I would write an elemental function. The source for such a function takes a scalar argument and provides a non-allocatable, non-pointer scalar result.

elemental function times2(x) result(y)
  real, intent(in) :: x
  real :: y

  y = 2*x
end function

I typically use deferred shape allocatable function results when the shape (or some other attribute) of the function result is not able to be described by a simple specification expression. Allocatable function results:

• requires writing to at least the Fortran 2003 standard;

• may require an extra heap memory allocation/deallocation pair above what is strictly necessary, which may (or may not) have performance implications;

• may not require the same use of the stack as the automatic result case, which may (or may not) avoid stack overflow problems at execution time.

Compiler implementation details (including compiler options) affect the comparison. In particular, differences in how a compiler manages the memory for temporaries may make the two approaches converge in terms of their requirements for stack and heap allocations.

Avoid explicit shape array dummy arguments, unless you have particular requirements.