Clarification on elemental subroutines

1

Hello,

I need some help on code that I wrote a while ago :wink: The function below samples from a normal distribution. It is elemental (and impure due to the call to random_number). At least for 1 dimensional arrays x it seems to does what I expect: x is sampled from a normal distribution with given mean and standard deviation. What I donā€™t understand is why real, dimension(2) :: rnd is ā€œbroadcastedā€ to something like dimension(2,N).

impure elemental subroutine math_normal(x,mu,sigma)

  real, intent(out) :: x
  real, intent(in), optional :: mu, sigma

  real, dimension(2) :: rnd


  call random_number(rnd)
  x = misc_optional(mu,0.0) &
    + misc_optional(sigma,1.0) * sqrt(-2.0*log(1.0-rnd(1)))*cos(TAU*(1.0 - rnd(2)))

end subroutine math_normal

is this ok? Or should I use the modification below?

impure elemental subroutine math_normal(x,mu,sigma)

  real, intent(out) :: x
  real, intent(in), optional :: mu, sigma

  real :: rnd1, rnd2


  call random_number(rnd1)
  call random_number(rnd2)
  x = misc_optional(mu,0.0) &
    + misc_optional(sigma,1.0) * sqrt(-2.0*log(1.0-rnd1))*cos(TAU*(1.0 - rnd2))

end subroutine math_normal

The general question is then: How are array variables treated in elemental functions/subroutines?

2

How did you detect that rnd becomes a two-dimensional array? I can imagine something like that happening, although only from a formal point of view, not from an implementation pont.

3

Maybe Iā€™m wrong, but my understanding of such an elemental routine is that the result obtained with an array has to be equivalent to the result obtained by looping on the elements:

real :: x(10)
call math_normal(x)
! gives the same result as
do i = 1, 10
   call math_normal(x(i))
end do

But the compiler is free to do whatever he wants (including replacing the unique call by the above loop) as long as this constraint is honored, and I donā€™t think that the standard specifies implementations details.

1 Like
4

The 2D array was just my interpretation. But @PierUā€™s interpretation makes more sense.

5

that makes sense and explains the observed behavior.

I thought more in terms of array evaluation, but then arrayed variables would not work at all because the dummy arguments can have any shape.