Array comprehensions in Fortran

Help
1

Hi, I have another fairly noob question hope y’all don’t mind!

tldr is — what’s the best way to do array comprehensions in Fortran?

Background:
In Python and Haskell (for example) you have list comprehensions and Julia has the same idea but array comprehensions. For example you can do something like

[a**3 for a in 1:10 if a**2 % 2 != 0]

to compute all cubes of a such that a squared is not even (completely arbitrary example). I’ve just used a generic syntax but basic idea works in all three.

I view these as essentially alternatives to either more explicit for/do loops or higher order functions like map or filter.

It seems to me that array comprehensions, a la Julia’s version of list comprehension with arrays, are reasonably compatible with the ‘spirit’ of modern Fortran in which eg whole array operations or implied do loops are often used in place of perfectly fine but somewhat verbose explicit do loops.

So, my question is — can this be done similarly in Fortran?

I know there are implied do loops, eg

[(i**3, i=1,10)]

but unfortunately I don’t think these can include a masking condition?

I also know there is pack, which is essentially the functional programming alternative to a list comprehension, but I think this requires constructing the whole array to filter first eg

a = [(i, i=1,10)]
b = pack(a**3, mod(a**2,2) /= 2)

I think this is required so that the same variable can be referenced in the two arguments of pack. This is mostly fine but possibly memory inefficient?

There is also the apparently obsolete forall array assignment which almost works:

forall (i=1:10, mod(i**2,2)/=0) a(i) = i**3

but the indexing of a is broken (only need to have a of size given by subset of i elements). Fixing this with a dummy incrementing by one for each unmasked value just becomes a loop I think…

The forall ‘replacement’, do concurrent, similarly allows masking but this makes indexing of the retained elements awkward and essentially reduces back to an explicit do loop.

So! I guess my question is — have I missed anything, or is array comprehension with masking not quite possible in Fortran currently and instead requires an explicit loop or a functional pack statement with a pre-constructed indexing vector?

Thanks!

4 Likes
2

I think that you have it right. The most common approach would be to expand out your example and use a explicit loop and a counter. PACK can be very useful but as you note, will require additional storage and depending on the compiler, can be relatively slow.

1 Like
3

Thanks! A little bit of a shame — an implied do that allows a mask (or perhaps better to say allows a filter/cycle condition as the goal is to avoid allocating an array cf pack) would be very cool imo! :wink:

5

The problem with building the array on the fly to avoid creating a full array first and then filtering out the unwanted values is that - inevitably - you will have to repeatedly allocate memory. So, on the one hand you do not spill memory (temporarily) and on the other hand you have a series of fairly costly operations.
I would say: use the approach that is simplest in terms of coding, which in my view is the use of pack as you show it. Only when that turns out to give problems - huge temporary arrays for instance, then go a more complicated route.
Consider this maxim: “Premature optimisation is the root of all evil” - 1691: Optimization - explain xkcd

2 Likes
6

True. I guess I was also thinking that the implied do loop construct (cf array based) doesn’t explicitly build an array on the fly but it is rather just that eg

(i, i=1,3)

is syntactic sugar for

1, 2, 3

which can then be passed to eg an array constructor. So then thinking about the similar eg

(i, i=1,3, i**2 /=4)

as syntactic sugar for

1, 3

(Again just making up a toy example).

But I’m also not 100 percent sure that what I wrote above is correct — if it’s syntactic sugar for the loop itself then I guess you’d need to build the the list up piece by piece in a temporary array or something anyway and reallocate? No idea really!

Regardless though I wouldn’t strictly say what I’m going for is premature optimisation, rather premature elegance — I happen to like the aesthetics of list/array comprehensions in Python/Haskell/Julia but I suppose there might be good (performance etc) reasons like those discussed above for not having the exact same construct in Fortran :slight_smile:

7

Gfortran can compile implied do-loops very slowly when they are used to set large arrays. On my PC it takes 9.1s on Windows and 6.4s on WSL2 to compile

program list
implicit none
integer, parameter :: n = 10**7
integer, allocatable :: vec(:)
integer :: i
vec = [(i,i=1,n)]
print*,vec(n)
end program list

The Windows version is GNU Fortran (GCC) 12.0.0 20211024 from equation.com, and the WSL2 version is 9.3.0 . Intel Fortran compiles the code very quickly. Gfortran compiles an analogous code with forall or a do loop quickly.

2 Likes
8

Well, it may take more programming underneath, but I have written sample programs that mimic lambdas (anonymous functions) in Fortran or that enumerate an integer data set without explicitly constructing the entire set first. You can find several at http://flibs.sourceforge.net/.

3 Likes
9

One could use do concurrent construct to achieve the goal and also automatically benefit from auto-parallelization features of compilers,

do concurrent(i = 1:10, mod(i,2) == 0)
    Array(i) = i**3
end do

But here is the follow-up question: if performance really matters here, is list-comprehension really the best solution (not just in Fortran, any language)? If performance does not matter, then pack is a good solution as you mentioned.

3 Likes
10

Any language that supports this feature must do one of the following:

  • Construct the entire mask array first
  • Use a linked-list structure to dynamically add elements via iteration
  • re-allocate the array for every included element

I suspect the most efficient method is the first one, as it would allow the minimal amount of allocations. I would be in favor of adding the capability to the language. It would essentially be “syntactic sugar” like the following.

lhs = [(some_expression(i), i = ub, lb, stp where condition(i))]

would be equivalent to

block
  logical :: mask(ub:lb:stp)
  type(typeof(lhs)), allocatable :: tmp(:)
  integer :: i, j, num_elements

  do i = ub, lb, stp
    mask(i) = condition(i)
  end do
  num_elements = count(mask)
  allocate(tmp(num_elements))
  j = 1
  do i = ub, lb, stp
    if (mask(i)) then
      tmp(j) = some_expression(i)
      j = j + 1
      if (j > num_elements) exit
    end if
  end do
  lhs = tmp
end block

I’m sure there’s nuances in the standard that need to be considered.

3 Likes
11

Does this work? What are the indices of the array? Is it preallocated of some size? I don’t think it’s quite right but I may be wrong…

12

Yeah I think this would some really nice syntactic sugar too :slight_smile:

Essentially just extending implicit do loops to allow a where clause

13

This would not be equivalent, because the semantics of the language are such that the [(..., i = lb, ub <, stp>)] is a temporary array whose lower bound is one, and whose upper bound is (ub - lb) / (present(stp) ? stp : 1). Adding the where clause should mean that it has upper bound count([(condition(i), i = lb, ub <, stp>)]), but in that example Arr must have bounds lb:ub<:stp>, and will contain elements that are undefined (or at least not redefined) in the places that condition(i) .eqv. .false.

With some rearranging I think I could change my example to use do concurrent, which is implied (I think) by the semantics of the language, but I’d have to double check, and it would no longer be quite as clear.

1 Like
14

I think you’re replying to my reply re: do concurrent :slight_smile: Just to be clear, I don’t think this works.

I’m not 100% sure but I I’ve come to the tentative view that do concurrent isn’t actually the right approach here due to the need to keep track of how many cases have been accepted so far (so it’s not a set of fully independent operations if you want proper indexing, though it can preserve the ordering with a bunch of extra unneeded indices. But speaking as a noob here)

So now mainly thinking of an ordinary implied do but with where condition

15

Looking at your use case more carefully now, I think there is a slight difference between what you desired and what do concurrent does in my comment, in the sense that it does not construct an array, but selectively fill elements that meet the condition. What I intended to say is well detailed in everythingfunctional's comment (that there is likely no escape from the performance penalties of masks in any language). But regarding do concurrent, here is an implementation,

integer :: Array(10) = 0
#if __GFORTRAN__
integer :: i
do concurrent(i = 1:10, mod(i,2) == 0)
    Array(i) = i**3
end do
#else
do concurrent(integer :: i = 1:10, mod(i,2) == 0)
    Array(i) = i**3
end do
#endif
print *, Array
end

Intel ifort supports the full do concurrent syntax, which is why I separated it from the gfortran code. Compile it with,

ifort -fpp concurrent.f90
gfortran -cpp concurrent.f90

Based on my experience, Intel ifort does indeed automatically parallelize do concurrent constructs over threads if it finds it beneficial. NVIDIA compiler can offload the work to GPUs as well. I have no experience with auto parallelization features of gfortran.

1 Like
16

Re: “this requires constructing the whole array,” what do you think Python, Haskell, Julia, etc. do “under the hood”?

By the way, if having to declare and define variables (and that too in a section preceding the executable one in a program) are your concern, you can consider ASSOCIATE construct with your PACK intrinsic solution:

   associate( a => [( i, integer :: i = 1, 10 )] )
      print *, pack(a**3, mod(a**2,2) /= 0)
   end associate
end

C:\temp>ifort p.f90
Intel(R) Fortran Intel(R) 64 Compiler Classic for applications running on Intel(R) 64, Version 2021.5.0 Build 20211109_000000
Copyright (C) 1985-2021 Intel Corporation. All rights reserved.

Microsoft (R) Incremental Linker Version 14.30.30706.0
Copyright (C) Microsoft Corporation. All rights reserved.

-out:p.exe
-subsystem:console
p.obj

C:\temp>p.exe
1 27 125 343 729

C:\temp>

2 Likes
17

I think Haskell explicitly does some sort of lazy evaluation thing?

18

Yes, Haskell uses lazy evaluation, so you can actually construct an infinitely long array (list in Haskell lingo) using list comprehension syntax like [x**2 -> x in 1.. if (x % 2) == 0] (I may have messed up the syntax a bit there) and it won’t execute forever.

Edit: the correct syntax in Haskell would be [ x^3 | x <- [1..], even x ]

1 Like
19

Yeah it’s fun to be able to work with infinite lists in Haskell :slight_smile:

And I guess this further highlights some differences between ‘iterate through as needed and filter as you go’ vs ‘pre-construct then filter’?

20

Well, Haskell is known for that, isn’t it? But it’s not some rocket science, besides the quality of compiler implementation of that functional programming thought-process varies quite a bit, particularly with real-world needs, so much so that there are more cons to consider with “under the hood” attempts at lazy evaluation.

Under the circumstances, with Fortran the practitioners can do better with DIY approaches, say for the case at hand:

   print *, "a(2) = ", a(2)
   print *, "a(5) = ", a(5)
   print *, "a(10) = ", a(10)
   print *, "a(20) = ", a(20)
contains
   function a(n) result(r)
      integer, intent(in) :: n
      integer, allocatable :: r(:)
      integer, save :: m = 0
      integer, allocatable, save :: dat(:)
      associate ( x => [( i, integer :: i = 1, n )] )
         if ( n > m ) then
            m = max( 2*n, 25 )
            associate ( y => [ x, ( i, integer :: i = n+1, m )] )
               dat = pack( y**3, mod(y**2, 2) /= 0 )
            end associate
         end if
         r = dat(1:count( mod(x**2,2) /= 0 ))
      end associate
   end function
end

C:\temp>ifort /standard-semantics p.f90
Intel(R) Fortran Intel(R) 64 Compiler Classic for applications running on Intel(R) 64, Version 2021.5.0 Build 20211109_000000
Copyright (C) 1985-2021 Intel Corporation. All rights reserved.

Microsoft (R) Incremental Linker Version 14.30.30706.0
Copyright (C) Microsoft Corporation. All rights reserved.

-out:p.exe
-subsystem:console
p.obj

C:\temp>p.exe
a(2) = 1
a(5) = 1 27 125
a(10) = 1 27 125 343 729
a(20) = 1 27 125 343 729
1331 2197 3375 4913 6859

C:\temp>

21

Sure, fair enough :slight_smile:

Again, my original motivation was — I like the simplicity and elegance of one line list comprehensions in other languages (they track math language closely), though I understand they’re not necessarily the most efficient etc. So was wondering if this was possible in terms of a simple one liner in Fortran, especially when I learned about implied do loops.

Implied do loops seemed promising for this use case but don’t quite work as they don’t (currently) allow a filter condition. A two line pack approach with a fully allocated iteration array is probably the cleanest looking, though slightly less aesthetically pleasing to me and not quite the same as a list comprehension. More sophisticated approaches like yours are certainly also nice, but IMO relatively verbose compared to the math concept motivating comprehensions. The associate idea seems cool though

1 Like