Implied-do array constructor, type specs, and differences between GFortran, Intel, and LFortran

[wspector]

@Beliavsky

One can define a binary operator .to. so that i.to.j gives what you intend with i:j

That is a neat solution.

I have a little routine I call ‘iota’ - after the APL operator of the same name. Besides a simple number of items, it has optional base and stride arguments. This is a “data parallel” version that uses my homegrown version of the HPF sum_prefix (“scan” in APL) function:

pure function iota (n, base, stride) result (res)
integer, intent(in) :: n
integer, intent(in), optional :: base
integer, intent(in), optional :: stride
integer :: res(n)

if (present (stride)) then
  res = stride
else
  res = 1
end if

if (present (base)) then
  res(1) = base
else
  res(1) = 1
end if

res = sum_prefix (res)

end function iota

[GregB]

Yes, I agree that APL’s iota is a good inspiration for this. Can we write a rank independent version in Fortran? Seems like reshape does what we need, but I can’t quite figure out the type declaration…

function iota32_rank_independent(array_shape) result(iota)
  ! reminiscent of APL ι operator
  ! iota([4]) == [1,2,3,4]
  ! iota([3,3]) == ⎡1 4 7⎤
  !                ⎢2 5 8⎥
  !                ⎣3 6 9⎦
  ! copy/paste for different integer kinds
  integer(int32), intent(in) :: array_shape(:)
  integer(int32), dimension(array_shape) :: iota  ! <-- What to put here?
  integer(int32) :: temp(product(array_shape))
  integer(int32) :: elements
  integer(int32) :: i
  elements = product(array_shape)

  do concurrent (i = 1:elements)
    temp(i) = i
  end do

  iota = reshape(temp,array_shape)
end function iota32_rank_independent

[PierU]

@everythingfunctional

@PierU

Just

[n:m:p] (:p being optional)

I get the feeling there’s some syntactic ambiguity hiding in here somewhere, but if not it seems reasonable.

I have been thinking about possible ambiguities for some time, but couldn’t really find one. An implied do is beween parenthesis (...), and a list of scalars and/or arrays does have the columns :

Breaking news: This syntax is accepted in the Intel compiler:

print*, [1:10:3]
print*, [5:3:-1, 100]
end

Gives

Program returned: 0
Program stdout

           1           4           7          10
           5           4           3         100

I’m going to submit a proposal

[snano]

by Flang also:

print*, [1:10:3]
print*, [5:3:-1]
end

output:

Program stdout

 1 4 7 10
 5 4 3

[PierU]

@snano

by Flang also:

print*, [1:10:3]
print*, [5:3:-1]

flang is apparently promoting this constructor to integer(kind=8), that’s why print*, [5:3:-1,100] won’t work. print*, [5:3:-1,100_8] is OK. Looks like a bug in the extension (if it’s intentional, I don’t think it’s a good idea).

[wspector]

@GregB

Yes, I agree that APL’s iota is a good inspiration for this. Can we write a rank independent version in Fortran? Seems like reshape does what we need, but I can’t quite figure out the type declaration…

In classic APL, one would use iota, then use reshape (rho) to turn the vector into a matrix. And yes, reshape is more or less the Fortran equivalent to APL’s reshape.

In Fortran, you’d need to write a function for each rank you want to support, then tie them together with a generic name. The dummy argument list is used to distinguish which version is selected. There are several ways it could be done. An easy way to write them would be to have a separate integer argument for each array dimension.

[RonShepard]

@wspector

In Fortran, you’d need to write a function for each rank you want to support, then tie them together with a generic name.

The RESHAPE() intrinsic already works with arrays of arbitrary type, arbitrary kind, and arbitrary rank. As you note, it is difficult to write user functions that are that flexible, so at least with the current language features, it seems like using directly the intrinsic RESHAPE() is the best approach.

There is another set of related gotchas that lurks in this swamp. The only way to allocate an array with nondefault bounds is with an explicit ALLOCATE() statement, so in these cases a simple assignment to the allocatable array is not sufficient, it takes at least two fortran statements. But further, even if the programmer allocates the target array with the right shape and bounds, the compiler is allowed to reallocate the lhs array within a simple array assignment statement with the default bounds. In order to avoid that, the assignment must be like the second assignment below.

allocate( array(lb1:ub1,lb2:ub2,lb3:ub3) )
array =  reshape( vector, shape=shape(array) )  ! possible error, reallocation of lhs is allowed.
array(:,:,:) = reshape( vector, shape=shape(array) )  ! correct.

There are, of course, many ways to write that last statement to achieve the same result.

There are three possible changes to the language standard that would affect this predicament. First, the language standard could specify that reallocation cannot occur when the shapes of the lhs and rhs match; this change would require the above two assignments to always be exactly the same. The downside for this is that the ability to reallocate is sometimes a useful low-level optimization, one that would no long be allowed with that change. Second, if there were a way to change the bounds of an existing allocatable array (without a copy, just with metadata modifications), then this would not be an issue. The move_alloc() intrinsic, which is the current way to transfer metadata between allocatable arrays, might be generalized to achieve this, but there are other ways it might be achieved too. Third, the reshape() intinsic could be generalized to allow nondefault bounds. This solves this particular problem involving nondefault bounds, but it would not affect the other cases where the programmer needs to change the bounds of existing allocatable arrays.

[urbanjost]

https://www.intel.com/content/www/us/en/docs/fortran-compiler/developer-guide-reference/2024-1/array-constructors.html

documents the ‘low:high[:stride]’ array constructor extension for ifx/ifort. I was curious if there were any discussion of the extension that might be useful in inspiring a proposal but I see no other features or extensions related to it, simply a few examples. No discussion of whether only default integer ranges are allowed, support of negative strides, result of M<N, etc. that I see.

[GregB]

Some empirical testing with ifx 2024.1.0:

program triplet
    use iso_fortran_env

    print *, [1:5]
    print *, [0:9:2]
    print *, [9:1:-2]
    print *, "size([7:1]): ", size([7:1])

    print *, [-3_Int64 : -1_Int64]
    print *, [1_Int64 : 3_Int16]
    print *, [1_Int16 : 3_Int32 : 1_Int64]
end program triplet

…gives the outputs…

           1           2           3           4           5
           0           2           4           6           8
           9           7           5           3           1
 size([7:1]):            0
                    -3                    -2                    -1
                     1                     2                     3
      1      2      3

[urbanjost]

integer :: a(11)
character(len=*),parameter :: all='(*(g0,1x))'
!print all, [1:10:3]*2
!print all, [5:3:-1, 100]
a=-999;FORALL (I = 1:11) A (I) = I;      print all,a
a=-999;FORALL (I = 5:-5:-1) A (I+6) = I; print all,a
a=-999;FORALL (I = 11:10) A (I) = I;     print all,a
end

gfortran xx.f90
./a.out
1 2 3 4 5 6 7 8 9 10 11
-5 -4 -3 -2 -1 0 1 2 3 4 5
-999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999

So that is similar behavior to the triplex syntax of FORALL (which is the only place in Fortran I can think of that uses something like the proposed syntax) and at least to me seems like the most intuitive behavior.

[everythingfunctional]

@RonShepard

allocate( array(lb1:ub1,lb2:ub2,lb3:ub3) )
array =  reshape( vector, shape=shape(array) )  ! possible error, reallocation of lhs is allowed.
array(:,:,:) = reshape( vector, shape=shape(array) )  ! correct.

This can be written in a single statement as

allocate( &
    array(lb1:ub1,lb2:ub2,lb3:ub3), &
    source = reshape( vector, shape=shape(array) ) &
    )

@urbanjost

which is the only place in Fortran I can think of that uses something like the proposed syntax

DO CONCURRENT also uses that syntax, i.e. do concurrent (i = 1:10:2).

[urbanjost]

Hah. I actually use that more than FORALL! Corrected.

[sblionel]

@RonShepard

First, the language standard could specify that reallocation cannot occur when the shapes of the lhs and rhs match;

Confused - the standard already says that.

[PierU]

@sblionel The paragraph you are citing says that “…[the variable] is deallocated if expr is an array of different shape…”: this is an “if”, not an “if and only if”.

[sblionel]

The standard specifies what happens for intrinsic assignment. Doing something that the standard doesn’t specify is non-conforming. I have serious doubts anyone would implement a Fortran compiler that chose to deallocate and reallocate in the absence of a shape difference (or the other things the standard calls out), as that would have additional side-effects.

[RonShepard]

@sblionel

I have serious doubts anyone would implement a Fortran compiler that chose to deallocate and reallocate in the absence of a shape difference (or the other things the standard calls out), as that would have additional side-effects.

@PierU has pointed out the problem with the language standard. I do not know of any compiler that does reallocation in this case, but it would be allowed according to the text of the standard. To give an example of when this might be done, suppose the lhs is an allocatable array with intent(out) or intent(inout) semantics. In this case, the compiler would “know” that reallocation would be allowed. Then suppose that the rhs is itself stored in some anonymous intermediate result allocated on the heap. The assignment could then occur either by copying the results from the rhs to the lhs, or it could do the equivalent of a move_alloc() and assign the metadata of the lhs to the memory on the rhs. That latter would be more efficient, so a compiler could do this as a low-level operation.

From our previous discussions on this topic in c.l.f., I had the impression that this possible optimization (reallocation of the lhs) was a feature, not a bug, and that the ambiguity in the language standard text was intentional. As I noted in my previous comment, and as @everythingfunctional also showed, the programmer has ways to avoid this reallocation, but the programmer must be aware of the situation in order to know to do the right thing. It does seem like an obscure gotcha to me.

As to when this could cause problems for the programmer, if the lhs has the target attribute, then any previously associated pointers would no longer work as expected, they would need to be reassociated with new pointer assignments.

[sblionel]

It would not be allowed - deallocation could cause finalization, as well as invalidation of any pointers to the LHS as a target. Hard no.

[RonShepard]

Wouldn’t finalization be allowed anyway if the lhs has intent(out) or intent(inout) attributes? In any case, if it is not allowed, then the language standard should say so explicitly and plainly. Similarly, if it is allowed, then I think the language standard should also say so explicitly and plainly.

But I also have a question about the text “the shape of expr with each lower bound equal to the corresponding element of LBOUND(expr) if expr is an array.” In the case of reallocation of the lhs, does this REQUIRE the bounds on the rhs side to be copied to the lhs?

[sblionel]

If no deallocation is done as part of the assignment, then no finalization. intent(out) would deallocate on procedure entry, intent(inout) would have no effect on this.

On reallocation, the bounds depend on what the rhs is. If it’s a whole array, then the bounds are copied, otherwise the shape is copied with the lower bound being 1 for each dimension. Example:

D:\Projects>type t.f90
integer, allocatable :: a(:,:), b(:,:)

allocate (a(2:4,3:5))
allocate (b(10:20,5:10))
print *, lbound(a)
print *, ubound(a)
a = b
print *, lbound(a)
print *, ubound(a)
end

D:\Projects>ifx /nologo t.f90

D:\Projects>t.exe
           2           3
           4           5
          10           5
          20          10

[PierU]

@sblionel

It would not be allowed - deallocation could cause finalization, as well as invalidation of any pointers to the LHS as a target. Hard no.

I think we all agree that it would be extremely weird and misleading if a compiler was doing that, but the question is about if some very liberal interpretation of the standard would allow it or not…

Admittedly, the standard neither explicitely says that a processor is disallowed to rellocate a variable (even not an allocatable one) that is in use… But obviously no compiler will do that.