I have an N-dimensional dataset (say real numbers) which is stored as a 1D array with an additional dimension array which specifies the original dimensions.
In addition, the functions to deduce the 1-D index from an N-D index and vice-versa are given.
I'm trying to figure out how do I make a do-loop (or equivalent) for the generic N-dimensional index (which will be transformed to 1D index of course) from some set of limiting lower indices to a set of upper indices. So I need an "N-dimensional" loop which does not go over all of the values - only a portion of the array, therefore doing the linear index of an equivalent 1D array is not relevant (at least without modifications).
This is a schematic of my problem:
subroutine Test(Array,Dims,MinIndex,MaxIndex)
implicit none
real , dimension(1:), intent(inout) :: Array
integer, dimension(1:), intent(in) :: Dims,MinIndex,MaxIndex
integer, dimension(size(Dims)) :: CurrInd
integer :: 1Dindex
! size(Dims) can be 1, 2, 3 ,..., N
! size(MinIndex)==size(MaxIndex)==size(Dims)
! size(Array)==Product(Dims)
! 1Dindex=Get1dInd(NDindex,Dims)
! NDindex=GetNdInd(1Dindex,Dims)
! How do I actually preform this?
do CurrInd=MinIndex,MaxIndex
1Dindex=Get1dInd(CurrInd,Dims)
<Some operation>
enddo
end subroutine
I figured it is possible to loop over the Dims array and use an inner loop but I don't manage to write the procedure down properly.
Another option that didn't work for me (maybe because I use it incorrectly?) is FORALL, as that requires each index to be specified separately.
If you knew the dimensions of your arrays at compilation time you could do a series of nested DO loops, each of them running between pairs of components of MinIndex and MaxIndex. Since you don't know the dimensions, that's not possible.
The easiest strategy I can think of is to loop with a single DO loop over all the 1D indices. For each of them, compute the N-dimensional index, and check if it is within the bounds provided by MinIndex and MaxIndex: if it is, continue doing what you need; if it is not, discard that 1D index and go to the next one. If the indices are sequential you may be able to do something smarter that skips blocks of indices that you know you will not be interested in.
do OneDindex = 1, size(Array)
CurrInd = GetNDInd(OneDindex, Dims)
if ((any(CurrInd<MinIndex)) .or. (any(CurrInd>MaxIndex))) cycle
! <Some operation>
end do
Note that, as far as the index manipulation is concerned, this strategy is compatible with parallelising the loop.
Also note that Fortran variables must start with a letter: 1Dindex is not a valid Fortran variable name.