Returning array from function

[leonfoks]

I was going to ask this in #113 but didn't want to derail that thread.
This is more of a personal question... since I've actually done very little with returning an allocatable array from a function. I may have flawed logic, so please correct me if i'm wrong.
I've had compile issues in the past returning allocatable from functions (cant remember all the details but perhaps others have had these too) and I "intuitively" thought there might be a slow down if those functions need calling multiple times.

Ill use the snippet from #113

function mean(mat, dim) result(res)
    real(sp), intent(in) :: mat(:,:)
    integer, intent(in), optional :: dim !dim = 1 or dim = 2 
    real(sp), allocatable :: res(:)
end function

If I have a large array, and I call this multiple times, does it reallocate memory for that variable every time? Or is there some cleverness that happens?

e.g. 1D result is already allocated, function allocates its own result, on return those results are copied to pre-allocated array.

real(sp), allocatable :: x(:), A(:, :)

! A = full of numbers
allocate(x(shape(A)(2)))

x(:) = mean(A, dim=1)

Is it known what happens specifically, memory wise? Is it compiler dependent?

I'm asking because it will be slower if allocations keep happening, we might want to be careful and always make subroutines available where preallocated memory can be accessed without any implicit allocations and then copies. Am I worrying about this too much haha.

[certik]

This is really important. I don't think we should use allocatable at all in the lowest level API, unless the size of the output is truly not known ahead of time, like in loadtxt. Also, in the past I have found that if a function returns an array (non allocatable), it can be slower than an equivalent subroutine with intent(out) array. Apparently there is an extra copy happening. This is something where it feels Fortran compilers must do better. I feel ideally there shouldn't be an overhead of using functions over subroutines, but unfortunately it seems there is.

…

On Sun, Jan 19, 2020, at 12:07 PM, Leon Foks wrote: I was going to ask this in #113 <#113> but didn't want to derail that thread. This is more of a personal question... since I've actually done very little with returning an allocatable array from a function. I may have flawed logic, so please correct me if i'm wrong. I've had compile issues in the past returning allocatable from functions (cant remember all the details but perhaps others have had these too) and I "intuitively" thought there might be a slow down if those functions need calling multiple times. Ill use the snippet from #113 <#113> function mean(mat, dim) result(res) real(sp), intent(in) :: mat(:,:) integer, intent(in), optional :: dim !dim = 1 or dim = 2 real(sp), allocatable :: res(:) end function If I have a large array, and I call this multiple times, does it reallocate memory for that variable every time? Or is there some cleverness that happens? e.g. 1D result is already allocated, function allocates its own result, on return those results are copied to pre-allocated array. real(sp), allocatable :: x(:), A(:, :) ! A = full of numbers allocate(x(shape(A)(2))) x(:) = mean(A, dim=1) Is it known what happens specifically, memory wise? Is it compiler dependent? I'm asking because it will be slower if allocations keep happening, we might want to be careful and always make subroutines available where preallocated memory can be accessed without any implicit allocations and then copies. Am I worrying about this too much haha. — You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub <#114?email_source=notifications&email_token=AAAFAWC6TNXY3SQSYYYBK3DQ6SQI7A5CNFSM4KI2655KYY3PNVWWK3TUL52HS4DFUVEXG43VMWVGG33NNVSW45C7NFSM4IHGL63Q>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AAAFAWBI63LLCUGKPWVAQNTQ6SQI7ANCNFSM4KI2655A>.

[leonfoks]

Phew, I was hoping I wasn't crazy haha

[nncarlson]

I don't believe there is any performance difference between an automatically-sized array valued function, for example

function foo(array) result(res)
real, intent(in) :: array(:,:)
real :: res(size(array,dim=2))

versus one whose result is allocatable. In either case the size of the result is determined at run time. So there is absolutely no reason to avoid allocatable results if one is determined to use a function.

Also, in the past I have found that if a function returns an array (non allocatable), it can be slower than an equivalent subroutine with intent(out) array. Apparently there is an extra copy happening.

This is the real problem with array-valued functions, allocatable or not, and has nothing to do with allocatable.

[milancurcic]

So there is absolutely no reason to avoid allocatable results if one is determined to use a function.

How about the allocation on stack vs. heap? Yes, this isn't defined in the language and is implementation specific, but I found that in practice there is a difference:

• Non-allocatable function result may be allocated on the stack, which is faster;
• Allocatable function result is always allocated on the heap, which is slower;

I never measured this, but I have anecdotal personal experience of WRF being noticeably slower when compiled with -heap-arrays (ifort flag to make all arrays be allocated on heap).

This is also the caveat with non-allocatable function (or subroutine) results (output arguments): For large arrays, stack can overflow. This problem doesn't occur with allocatable arrays.

[nncarlson]

I can certainly believe there is a performance difference between stack and heap allocated. But in the case of an automatic array result that must be allocated at run time, do you think the compiler inserts code to decide where to allocate (stack or heap) based on the size it determines for a particular invocation of the function? (Note I'm not talking about a size known at compile time.) I suppose that could be the case; do you know? But even if so, I think the far bigger cost is the one @certik pointed out with using a function at all.

[certik]

We should setup a project to carefully benchmark all these options. Note that gfortran and Intel can put an allocatable on a stack as well.

…

On Mon, Jan 20, 2020, at 9:59 AM, Milan Curcic wrote: > > So there is absolutely no reason to avoid allocatable results if one is determined to use a function. How about the allocation on stack vs. heap? Yes, this isn't defined in the language and is implementation specific, but I found that in practice there is a difference: * Non-allocatable function result may be allocated on the stack, which is faster; * Allocatable function result is *always* allocated on the heap, which is slower; I never measured this, but I have anecdotal personal experience of WRF being noticeably slower when compiled with `-heap-arrays` (ifort flag to make all arrays be allocated on heap). This is also the caveat with non-allocatable function (or subroutine) results (output arguments): For large arrays, stack can overflow. This problem doesn't occur with allocatable arrays. — You are receiving this because you commented. Reply to this email directly, view it on GitHub <#114?email_source=notifications&email_token=AAAFAWCIENGWONBPWR32MH3Q6XJ5JA5CNFSM4KI2655KYY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGOEJNIWIY#issuecomment-576359203>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AAAFAWAEIKYND6XWHNJ4S3TQ6XJ5JANCNFSM4KI2655A>.

[nncarlson]

We should setup a project to carefully benchmark all these options.

That would be very interesting to see. And there are many other specific performance questions like this. Conclusions from such a project could feed into a "Fortran best practices" guide.

Note that gfortran and Intel can put an allocatable on a stack as well.

Yep. My experience with intel is that I routinely have to unlimit the stack size to avoid hitting a segfault. I have no such issues with NAG, for example.

[certik]

Also I want to setup benchmarks for operations on arrays (y(:) = a(:) + b(:)) compared to an explicit loop. The loop seems to be faster some times and I have seen several codes moving away from arrays to explicit loops because of that.

Those are things where I want Fortran compilers to do better. And a first step is to have reliable benchmarks, showing where the "performance bug" is.

[leonfoks]

In that benchmark, a manually unrolled loop version would be cool to see. The compilers should already do a good job of this. It would be nice to see in a benchmark though if it’s worth it?

[milancurcic]

This looks like a good candidate for a fortran-lang/benchmarks repo. This could include both benchmarks of stdlib, and more fundamental basic operations discussed in this thread. It would be nice to have a pipeline there that automatically generates a report for all procedures/programs in the suite. Later on when we change implementations of some stdlib procedures that are used by other procedures, it would be useful to verify that there's no significant performance regression when changing the implementation of a building block.

Besides these, I'm also curious about what's the overhead of calling a generic procedure that overloads some 50 or 100 specific procedures.

[zerothi]

Additionally one should be careful with allocatable statements, the below is a perfectly reasonable code:

real, allocatable :: a(:)

a = [1, 2, 3]
a = [1, 2, 3, 4]

However, the compiler have to add checks for sizes and do deallocation/allocation in case they aren't commensurate.
To force the compiler to not check one has to add this:

a(:) = [1, 2, 3]

in which case it will segfault if the dimensions doesn't match.
If this library is low-level enough it should not rely on automatic allocation AND it should disallow array dimension checks (i.e. use a(:) = )