To: J3 J3/19-110r1 From: Van Snyder Subject: Supporting rank genericity -- NOT VECTOR SUBSCRIPTS! Date: 2019-January-16 References: 04-195, 18-247 19-110r1 adds a use case related to generic programming. Introduction ============ 1. In support of facilities in Fortran 2018 ------------------------------------------- MAXLOC(A) returns a rank-1 array whose extent is equal to the rank of A. It cannot be used directly as a subscript for A. One must assign it to a variable of the appropriate rank, and use elements of that array as subscripts. For example integer :: S(3) real :: A(10,10,10) ... s = maxloc(a) print *, 'Maximum value in A is ', a(s(1),s(2),s(3)) LBOUND(A) and UBOUND(A) return rank-1 arrays whose extents are equal to the rank of A. They cannot be used directly as bounds for an automatic or allocatable variable. For example real, intent(in) :: A(:,:,:) real :: X(ubound(a,1),ubound(a,2),ubound(a,3)) or real, allocatable :: X(:,:,:) ... allocate ( x(ubound(a,1),ubound(a,2),ubound(a,3)) ) or, if you need a halo around all dimensions of X allocate ( x(lbound(a,1)-1:ubound(a,1)+1, & & lbound(a,2)-1:ubound(a,2)+1, & & lbound(a,3)-1:ubound(a,3)+1) ) These are tedious if you know the rank, and impossible if we develop facilities for rank-generic programming, except perhaps by using a MACRO DO, which results in unreadable code. History ======= At meeting 217 on Monday 15 October /JOR announced without plenary discussion that there would be no further discussion of 18-247. Meeting 217 minutes state "Some features interesting; essentially vector subscripts." The proposal is NOT ESSENTIALLY VECTOR SUBSCRIPTS. This is EXPLICITLY STATED in 18-247, notwithstanding that the specifications concerning the relationship of vector subscripts to component initialization, association, finalization, contiguity, and I/O, would also apply. Read this paper all the way to the end, and read the referenced paper 04-195 for more details. Proposal ======== 1. Array bounds use case ------------------------ Allow a single array B of rank one and extent R to specify the upper or lower bounds of a rank R array A. The extent of B shall be a constant, except maybe for the A(..) case. For example real, intent(in) :: A(:,:,:) real :: X(ubound(a)) or real, rank(A), allocatable :: X ... allocate ( x(ubound(a)) ) or (for the halo case) allocate ( x(lbound(a)-1:ubound(a)+1) ) The case of one bound being a scalar and the other not is obvious: allocate ( x(0:ubound(a)) ) is the same as allocate ( x([(0,i=1,rank(a))]:ubound(a)) ) Indeed, the default case (above) of lower bounds being all 1 allocate ( x(ubound(a)) ) is the same as allocate ( x(1:ubound(a)) ) The halo case above becomes allocate ( x(lbound(a)-1:ubound(a)+1) ) If rank genericity is provided by the generic programming facility, there is no other way to specify bounds for a generic-rank array, except by something equivalent to a MACRO DO. This results in unreadable code. This would also be useful within a BLOCK within a SELECT RANK construct. 2. In support of generic programming ------------------------------------ This illustration presupposes generic programming is provided by parameterized modules. module M ( R ) parameters integer :: R ! a scalar integer used to specify array ranks end parameters contains subroutine P ( A ) real, intent(inout), rank(r) :: A ! an assumed-shape array, for ! which a declaration using existing syntax could in principle be ! created by a MACRO DO statement or construct. A subsequent ! reader might need to do some reverse engineering to understand ! the effect. integer :: S(r) ! S can be used as a subscript for A real :: X([shape(A),2]) ! rank(X) = 1 + rank(A) integer :: S2(r+1,2) ! S2 can be used as a subscript for X .... If you don't understand the implications of this, PLEASE ASK FOR MORE EXPLANATION! There is nothing in DIMENSION(..) and SELECT RANK that provides any support in this direction. 3. Subscripts ------------- Allow a single array of rank one and extent R as the subscript of a rank R array. For example real :: A(10,10,10) ... print *, 'Maximum value in A is ', a(maxloc(a)) If S = [ 3, 4, 5 ], A(S) is a scalar, the same as A(3,4,5), not A(3:3,4:4,5:5). The obvious generalization here is to allow a single subscript S of rank K+1 and extents (R,N1, ..., NK) as a subscript for a rank-R array A. The elements of the rank-one sections in the first dimension of S are used consecutively as subscripts for A, as in the rank one case, resulting in a rank K object with shape [ N1, ..., NK ]. If K == 0, the result is a scalar. This provides a MORE GENERAL SCATTER/GATHER FACILITY THAN THE PRESENT VECTOR SUBSCRIPT FACILITY. This is NOT THE SAME AS USING THE ELEMENTS OF THE RANK-ONE SECTIONS IN S AS VECTOR SUBSCRIPTS FOR A, which would result in a rectangular section of shape [n,n] (in the case A is of rank 2 and S is of shape [2,n]). Example: Suppose we have arrays A3 with dimensions (10,10,10) and S3 with dimensions (3,2). If we assume [ 3 6 ] S3 = reshape( [3, 4, 5, 6, 7, 8], [3,2] ) = [ 4 7 ] [ 5 8 ] then A3(S3) is a rank-1 extent-2 array that can appear in a variable- definition context (except as an actual argument associated with a dummy argument having INTENT(OUT) or INTENT(INOUT)); it specifies the same array as [ A3(3,4,5), A3(6,7,8) ], which cannot appear in a variable-definition context because it is not a variable. This is different from A3(S3(1,:),S3(2,:),S3(3,:)), which can appear in a variable-definition context, but is an object with shape [2,2,2], not [2]. The former is an arbitrary collection of elements of A3, while the latter is a rectangular section of A3. The degenerate case (same as above) is real :: A(10,10,10) integer :: S(3) = [ 3, 4, 5 ] Then A(S) is the same as A(3,4,5), not A(3:3,4:4,5:5). This illustrates that a rank K+1 subscript with first extent R for a rank R array produces a rank K result. In this example, K = 0, so A(S) is a scalar, not a rank-3 array with shape [1,1,1]. I.e., S is NOT A VECTOR SUBSCRIPT. But many of the same restrictions for vector subscripts remain, e.g. not being an actual argument corresponding to an INTENT(OUT) dummy argument. This can probably be made to work with DIMENSION(..) arrays. A run-time check is probably required to make sure that in A(S), the first extent of S is the same as the (dynamic) rank of A. This might be especially useful to create automatic variables with the same shape as DIMENSION(..) arrays. If rank genericity is provided by the generic programming facility, there is no other way to specify subscripts for a generic-rank array, except by something equivalent to a MACRO DO, which would appear ubiquitously, in every reference, not only in declarations. This results in unreadable code. Let's not do kludges and half measures, which is how we got nonpointer optional arguments corresponding to disassociated pointers being absent. The original solution proposed for the combinatorial-explosion problem was a special case of the conditional expression proposal, with an empty ELSE part meaning "if the logical expression is true, the argument is the consequent (not its value), else the argument is absent." Use Case ======== Compute a list of sets of subscripts where an array has nonzero values. Here is an elegant rank-generic parallel solution that also depends upon generalizing the CRITICAL construct to apply to the relationship between iterations of a DO CONCURRENT construct (or provide a MONITOR procedure for the same purpose). integer :: Subs(rank(a),size(a)) integer :: Num_NZ num_nz = 0 ! The I below is an array having the same extent as ! the rank of A. It takes on all values of sets of subscripts of A. do concurrent ( integer :: I(rank(a)) = lbound(a) : ubound(a) ) & & shared ( num_nz ) if ( a(i) /= 0 ) then critical num_nz = num_nz + 1 ! Presently prohibited subs(:,num_nz) = i end critical end if end do print '(a/1p5g15.8)', 'Nonzero elements of A:', a(subs(:,:num_nz)) or write ( *, '(a)' ) 'Nonzero elements of A:' do i = 1, num_nz write ( *, '(a,i0,*(",",i0:))', advance='no' ) "(", subs(:,i) write ( *, '(a, 1pg15.8)' ) ") ", a(subs(:,i)) end do Here is an equivalent rank-generic but not parallel solution using an extension of non-concurrent DO. The extension allows the to be an array with the same extent as the loop bounds. ! The I below is an array having the same extent as the ! rank of A. It takes on all values of sets of subscripts of A, in ! array-element order. do ( integer :: I(rank(a)) = lbound(a), ubound(a) ) if ( a(i) /= 0 ) then num_nz = num_nz + 1 subs(:,num_nz) = i end if end do If these could be generated by a MACRO DO, the result would be truly unreadable, and fragile. The DO nesting depth would depend upon the rank of the array.