Partially elemental functions

Elemental functions are a neat feature of Fortran but are limited by the requirement that all arguments be scalars. Could this be relaxed to allow a mix of scalar and array arguments? A partially elemental function would be elemental for the scalar arguments, which in the caller would be scalars or conformable arrays as they are now. A function such as

elemental function moment(power,offset,x) result(y)
integer, intent(in) :: power
real   , intent(in) :: offset
real   , intent(in) :: x(:)
real                :: y
y = sum((x-offset)**power)/max(size(x),1)
end function moment

would be allowed. It could be called with moment(1,2.0,x) or moment([1,2],2.0,x) or moment(1,[0.0,1.0],x) or moment([1,2],[0.0,1.0],x) . The result would have the same shape as the expressions passed as arguments power and offset. Currently one can get similar behavior by defining a derived type with an array component and using that as a function argument, since that is considered a scalar. I don’t see why this should be necessary.

3 Likes

Actually, with the assumed rank feature you can simulate this. Okay, it is a trifle contrived, as the essence of the function is repeated, but it works. I only looked at the case where x is an array of one or two dimensions and all other arguments are scalar. Here is my code:

module moments_mod
    implicit none
contains

function moment(power,offset,x) result(y)
    integer, intent(in) :: power
    real   , intent(in) :: offset
    real   , intent(in) :: x(..)
    real                :: y

    select rank (x)
        rank(1)
            y = sum((x-offset)**power)/max(size(x),1)
        rank(2)
            y = sum((x-offset)**power)/max(size(x),1)
        rank(3)
            y = sum((x-offset)**power)/max(size(x),1)
        rank default
            y = -999.0
    end select
end function moment

end module moments_mod

program test_moments
    use moments_mod

    real, dimension(10)  :: x
    real, dimension(2,5) :: z
    real                :: offset

    x = [1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0]
    z = reshape( x, shape(z) )

    write(*,*) moment(1,5.0,x), moment(1,5.0,z)
    write(*,*) moment(2,5.0,x), moment(2,5.0,z)
    write(*,*) moment(3,5.0,x), moment(3,5.0,z)
    write(*,*) moment(4,5.0,x), moment(4,5.0,z)
end program test_moments
1 Like

That is different functionality. I am talking about functions where one or more arguments are arrays of fixed rank, in addition to scalar arguments, and I am proposing that the function be allowed ELEMENTAL in the scalar arguments.

While I understand the scalar argument restriction, I also wish there was someway to also pass arrays. One way to do it is to embed the array in a derived type but I’m not sure what the performance penalty is for that. On the surface an obvious fix might be to allow a finer grain elemental capability where individual dummy arguments could be declared to be elemental. Ie in Arjen’s example, Integer, elemental, intent(in) :: offset etc. but the compiler developers probably have good reasons why that wouldn’t work.

Also, on the subject of elemental function performance, has anyone done a study or know of one where the performance of using elemental functions vs embedding the work inside do loops was measured. I did my own quick test several years ago with a function that did a multiply of two values plus an add for arrays of around 1 million elements and didn’t see much of a difference in response time. (processors these days are so much faster than when some of the features like elemental functions were introduced in the language that doing a meaningful comparison is sometimes difficult).

I’m particularly interested in the case where you are passing arrays of derived types that contain other arrays. An example would be in a Finite Volume code where you have a Cell of Face type/class and you do some operation on a per cell or face basis. Writing the operations for just one cell or face would make things a little less complicated but I’m not sure what the performance hit would be.

What I do in such situations is: define the elemental function using only the scalar arguments, while implicitly passing the fixed rank argument using the scope of the containing unit.

Yes, that sounds more practical but how will you implement it?

e.g., If I define an elemental function in a module, yes, we can use module variables directly in elemental function, but they are static.
if I define an elemental function as type-bound procedure, then it treats class(myType)::this as an elemental argument, not a fixed one. That don’t work.

I just wonder if there is any more flexible way…

Thanks