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?

Awaken

Most people I know would use a subroutine, not a function. Functions that modify their arguments are just evil - they cause confusion, and it is my understanding that it is not always well defined exactly what they do when used in more complex expressions. So if you must use a function use the first way, but I think most people would say a subroutine is the way to go. – Archoplasm

Why should one be "easier" for the caller? – Bain

I agree with @IanBush's comment, up to a point. The point is Why define a function to multiply an array by 2 ? when that capability is built into the language. – Athenian

to multiply an array by two was just to make the example code shorter. The actual question I have is "say I have a subprogram that takes an array of unknown dimension, do a bunch of things and generates another array of the same dimension, how should I do it?". (I'm using function and not subroutine because the purpose of my subprogram is to generate another array based on the first one, not to modify the first one) – Awaken

Missed editing time - in above comment of course the bit about modifying arguments is irrelevant here, I misread the code. – Archoplasm

The actual question I have is .... Usually better to ask the question you want answered. – Athenian

Under the hood the function results are normally implemented by creating a subroutine with a hidden argument. Nevertheless, you may end-up with a temporary on the stack, if you use the function inside a complex expression. But that can happen any time you use complex expressions with arrays. – Hypoacidity

I don't thing transforming it to a subroutine is the best practice @roygvib. I have some reasons to do it as a function. I function can be called in place. sqrt is a function, that means you can evaluate expresions like sqrt(1/(1 + sqrt(x)). Imagine that sqrtwas a subroutine... – Awaken

@VladimirF I understand that you are saying that in terms of memory there is not much difference, right? – Awaken

@Bain I think the first one is easier for the caller because it doesn't have to make sure that x and n are consistent the caller passes x and everything is consistent inside the function. – Awaken

I think I would use an automatic return array for functions like cross(u, v) (for a cross-product of two vectors returning another vector), while use an allocatable return array for functions like matmul(A, B) if A or B can be pretty large (and so they may not fit on the stack). – Hexaemeron

so you would write cross(u,v,n) ? (sorry if I missunderstood you @roygvib) – Awaken

But you can rewrite the second to have y(size(x)) with x(:) still (which I admit is what I thought you'd done: I probably skipped over the actual code based on the most minimal difference between automatic and deferred.) – Bain

thank you @francescalus, I didn't know I could use a function when declaring the shape of an automatic array – Awaken

@ManuelPena I would probably stick to assumed-shape arrays (i.e. cross(u,v) rather than cross(u,v,n)) because one can pass non-contiguous arrays with smaller overheads (vs, copy-in may occur for explicit-shape dummy). But there might be cases where explicit-shape is better (I cannot recall such an example at the moment, sorry...) – Hexaemeron

I think one exceptional (and important) case where explicit-shape dummy arrays are useful is when one wants to use that routine from C (or other languages). – Hexaemeron

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.

Korfonta

Thank you, I didn't know you could specify the shape of the output as size(x), but knowing that I would use that also. I was talking about a non-elemental one, thank you anyway. – Awaken