Standard Algorithms: iota

[ashermancinelli]

Motivation

General algorithms have become ubiquitous in GPU and parallel programming, leading many developers to write code in terms of a set of standard algorithms that are broadly present in programming languages today. For example, the C++ standard library provides the function std::iota to make a sequence out of an array, which might be called from Fortran like so:

program main
  use std_algorithm
  implicit none
  integer, parameter    :: N = 5
  integer               :: i
  integer, dimension(N) :: a
  a = 0
  call iota(a)
  print *, (a(i), i=1, N)
end program main

$ ./example
           0           1           2           3           4

A draft of this algorithm may be found below.

module std_algorithm
  implicit none

  interface iota
    module procedure iota_i, iota_r
  end interface

contains

  subroutine iota_i(array, init)
    integer, dimension(:), intent(out) :: array
    integer, intent(in), optional :: init

    integer :: i, init_value = 0

    if (present(init)) init_value = init

    do i=1, size(array)
      array(i) = i + init_value - 1
    end do
  end subroutine iota_i

  subroutine iota_r(array, init)
    real, dimension(:), intent(out) :: array
    real, intent(in), optional :: init

    integer :: i
    real :: init_value = 0

    if (present(init)) init_value = init

    do i=1, size(array)
      array(i) = i + init_value - 1
    end do
  end subroutine iota_r

end module

Providing a set of standard algorithms should lower Fortran's barrier to entry for programmers already familiar with these concepts. In my opinion, the highest-priority algorithms after iota are scan and reduce. iota was chosen because of its simplicity and ubiquity.

Prior Art

• C++ STL's std::iota
• Thrust's make_sequence
• Python's range function

Additional Information

No response

[ashermancinelli]

I am happy to work on this, scan, and reduce if this idea is found acceptable by project maintainers.

[awvwgk]

I find the name iota somewhat unintuitive, especially if one has to look up the origin from the APL function ι, the letter iota, to make sense of the name. Therefore, make_sequence or range seem preferable, since they are also descriptive about the provided functionality.

Note that we already have the function arange.

[14NGiestas]

You can do the the following in fortran with array constructors, how different a function like iota would be?

program main
  implicit none
  integer, parameter    :: N = 5
  integer               :: i
  integer, allocatable  :: a(:)
  a = [(i,i=1,N,2)] ! i = start, end, steps
  print*, a
end program

[awvwgk]

As a side note, the implicit saved variable init_value might result in surprising behavior for the proposed implementation:

use std_algorithm, only : iota
implicit none
integer :: a(5), b(5), c(5)
call iota(a)
call iota(b, 1)
call iota(c)

print *, all(c == a)  ! F
print *, all(c == b)  ! T
end

[ashermancinelli]

@awvwgk Ah, I handn't found arange, thanks for pointing that out. Iota may not be the most accessible name, and the saved variable is certainly an issue in my naïve module.

@14NGiestas My example iota can work on array slices which may be more helpful than the array constructors.

Are there existing versions of scan or reduce? If not, are these worth creating another issue for? I appreciate both of your responses either way.

[awvwgk]

Are there existing versions of scan or reduce? If not, are these worth creating another issue for? I appreciate both of your responses either way.

No, I don't think we have those yet. Feel free to open an issue and we can discuss the API.

[ivan-pi]

I also find the name iota kind of cryptic. My favorite range command is probably the MATLAB colon where you just use 1:n. I use a similar function in some of my own codes which is a just a wrapper of the implied do loop (i,i=1,n).

In general I'm very supportive of "mining" the C++ standard template library for algorithms. Granted, the lack of generics makes this feasible only for the intrinsic types. Although with some preprocessing, you could perhaps extend this also to derived types. Arjen Markus posted an interesting template for reversing an array at Discourse. It gives a hint of what a generic Fortran template library might look like.

If a compiler vendor would ever want to re-use their C++ STL implementation, they could do so easily by making an STL-compatible wrapper of the Fortran array descriptor (I've made a few steps in this direction in #325).

[ashermancinelli]

Great. I'll open another issue for reduce and we can continue the conversation there.