Fortran - Function vs subroutine performance

Question

Since a couple of years ago or so I was completely new to Fortran, I overused SUBROUTINEs with no arguments, together with shared data, so that these procedures have made calculations on the actual arguments, available through USE statements. Now that I need to reuse some of these procedures (think of computing the divergence in a volume, a big DIMENSION(:,:,:) array, from a vector field in that volume, 3 big DIMENSION(:,:,:) arrays glued together in a derived type), I'd like either to

• keep them SUBROUTINEs but remove USE statements and use IN/OUT/INOUT dummy arguments (easy), or
• transform them in FUNCTIONs (little harder since I've to study a bit)

I imagine there could be a difference in performance in the two approaches, that I'd like to understand. In the following MWE I wrote 3 procedures to do the same computation, but I have no idea of how I should choose one or the others; neither I know whether some other approach would be preferable.

As a note, all rank-3 actual arrays in my program are ALLOCATABLE and must be so.

PROGRAM mymod

    IMPLICIT NONE

    TYPE blk3d
        REAL,    DIMENSION(:,:,:), ALLOCATABLE :: values
    END TYPE blk3d
    TYPE(blk3d) :: A, B

    INTEGER, PARAMETER :: n = 2
    INTEGER :: i

    ALLOCATE(A%values(n,n,n))
    A%values = RESHAPE([(i/2.0, i = 1, PRODUCT(SHAPE(A%values)))], SHAPE(A%values))
    print *, A%values

    ! 1st way
    B = myfun(A)
    print *, B%values

    DEALLOCATE(B%values)

    ! 2nd way
    ALLOCATE(B%values(n,n,n))
    CALL mysub(A, B)
    print *, B%values

    DEALLOCATE(B%values)

    ! 3rd way
    ALLOCATE(B%values(n,n,n))
    CALL mysub2(A, B%values)
    print *, B%values

CONTAINS

  FUNCTION myfun(Adummy) RESULT(Bdummy)                                                                                                              
    IMPLICIT NONE

    TYPE(blk3d), INTENT(IN) :: Adummy
    TYPE(blk3d)             :: Bdummy

    ALLOCATE(Bdummy%values, mold = Adummy%values)
    Bdummy%values(:,:,:) = 2*Adummy%values

  END FUNCTION myfun

  SUBROUTINE mysub(Adummy, Bdummy)

    IMPLICIT NONE

    TYPE(blk3d), INTENT(IN)    :: Adummy
    TYPE(blk3d), INTENT(INOUT) :: Bdummy

    Bdummy%values(:,:,:) = 2*Adummy%values

  END SUBROUTINE mysub

  SUBROUTINE mysub2(Adummy, Bdummy)

    IMPLICIT NONE

    TYPE(blk3d),            INTENT(IN)  :: Adummy
    REAL, DIMENSION(:,:,:), INTENT(OUT) :: Bdummy

    Bdummy(:,:,:) = 2*Adummy%values

  END SUBROUTINE mysub2

END PROGRAM mymod

EDIT In my program to do CFD, I use several rank-3 big-sized arrays. These arrays interact with each other in that computations (not simply pointwise +/-/*,...) is performed on some of them to obtain others of them. Think of B which is computed based on A through one of the 4 procedures in the example and then used to upgrade A itself by A = A + B. Am I wrong to think that the 4 options above could accomplish the same task? In this sense I could simply call A = A + myfun(A) if I chose the function approach.

EDIT2 The actual derived type, along with that big rank-3 array, has half a dozen of other fields (scalars, and small static arrays); additionally, most of the variables of this type are arrays, e.g. TYPE(blk3d), DIMENSION(n) :: A.

Floorer

"Performance" is perhaps not the thing to worry about at the moment. Those subroutines mean different things regarding the effect/requirements on arguments. So, before we answer on performance (don't forget you are better placed to analyse according to your specific setting) are you sure you understand the differences? – Finnegan

@francescalus, I've edited the question. Maybe I should in some cases replace OUT with INOUT to deal with arrays allocated outside the procedures? To answer to your question: no, the differences are not clear to me. What I want to do is obviously clear(er), but your objections could help me to have more insight and lead me to rethink the problem. – Floorer

@francescalus, I've updated the question with a more realistic example (it's no more minimal I fear, but maybe is the only way for it to be not too simplistic). – Floorer

Why have a type with a single field? Just use the array directly. Unless you are planning to write some type-specific methods in the future. – Hectogram

Perhaps your functions can be elemental instead of copying the entire array to process. – Hectogram

No, the actual relation between A and B is much more than B = 2*A, which is a toy line of code. – Floorer

Answer 1

Choosing between subroutine or function should generally be based on how you will be using the result and how clear it is to the reader.

What you do want to be concerned with is how many times data gets unnecessarily copied. With big arrays, you would want to reduce that.

Ignoring the actual work in the procedures, myfun copies the data a second time and (potentially) does two allocations. First the function result variable gets allocated and the data copied to it. Then back in the caller, B%values gets reallocated if it isn't already the same shape as the result and the data copied again, then the function result is deallocated.

mysub and mysub2 don't have this extra allocate/copy and are pretty much equivalent, though the call to mysub2 may have some extra work setting up a descriptor on the stack. I expect that to be noise if the subroutine does any real work.

Choosing between mysub and mysub2 really depends on how clear it is in your real application. Having a derived type with only one component seems unrealistic unless you are looking to have an array of these.

Perspiratory

Thank you, @SteveLionel! I'll comment the paragraphs. (1) It'd be much clearer, in my opinion, if I used FUNCTIONs, but clarity is not the sole concern, thus the question. (2) Yes, unnecessary copies is what I'm worried about, since those arrays can be very large in size. (3) So I could avoid the second allocation for the FUNCTION approach by pre-allocationg B%values, right? With this done, would this approach be still "inferior" to the other two (as I think you hint at)? – Floorer

(4-5) What do you mean by it in how clear is it? The actual derived type, along with that big rank-3 array, has half a dozen of other fields (scalars, and small static arrays); additionally, most of the variables of this type are indeed, as you have envisaged, arrays (when I started the project, that seemed the only way to emulate MATLAB's cells, which I was used/addicted to use, and it still seems the only way to me). – Floorer

The problem here with the function is that the semantics are that the function is evaluated and then the results copied to the variable on the left of the assignment, There isn't anything you can do in advance to avoid that copy. You'll have two allocates no matter what you do there - one for the function result and then one for the array in the variable. Given what else you have said, I suggest a subroutine where you pass the derived type that has the various components you want to reference or update. – Perspiratory

By "clear" I meant that if someone else read the code, would they understand it? Using a function just for function's sake is no help if the rest of the code is more complex. – Perspiratory