Partially elemental functions

1

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
2

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
3

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.

4

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.

5

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.

6

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