CODEWITHSUNDEEP
iteratorfortran
Stein
Stein

Reputation: 3277

Replace nested loops to be n-dimensional

I have a nested loop which is supposed to iterate over a an N-dim cube (line, square, cube,...)

So far I do this using nested loops, eg. for 2D:

do i = 1,Ni
   do j = 1,Nj
      ! do stuff with f(i,j)
   enddo
enddo

or 3D:

do i = 1,Ni
   do j = 1,Nj
      do k = 1,Nk
      ! do stuff with f(i,j,k)
      enddo
   enddo
enddo

I would like to replace these nested loops with a construct that works for every possible N-dimensional case. How can I do this in Fortran? In C++ one can use iterators. Is there something like this in Fortran?

I am considering writing a class, which kind of works like a mechanical n-dimensional counter (i.e. if we count 3 per dimension):

0,0,0,0,0
1,0,0,0,0
2,0,0,0,0
0,1,0,0,0
1,1,0,0,0

Would this be somewhat fast, compared to nested for loops? Do you guys have better ideas?

Upvotes: 1

Views: 251

Answers (2)

Stein
Stein

Reputation: 3277

I decided to use a super index. The size of my dimensions is limited to powers of two. This allows me to use bitshift operations and regular addition. Hopefully this is faster than a object-oriented counter.

  pure function supI(j, k, n) result(vec)
    implicit none

    integer*8, intent(in)         :: j
    integer*8                     :: copy
    integer*4, intent(in)         :: k, n
    integer*4                     :: vec(n), x, mask
    ! transforms super index j to a n-dimensional iteration,
    ! where each dimension has a range of 2^k

    ! mask to extract left most indicies
    mask = ISHFT(1, k) -1
    copy = j

    do x = 1,n
      ! extract current value using mask
      vec(x) = and(copy, mask)

      ! shift to next value
      copy = ISHFT(copy, -k)
    enddo
  end function supI

A 3d counter with 4 in each dimension would look like this:

  do j = 0,64
    write(*,*) supI(j, 2, 3)
  enddo

Edit: Compared to loop speed without optimization it is not to bad: ~20% more time, but If I use the optimizer it is much much slower than nested loops.

Upvotes: 1

Holmz
Holmz

Reputation: 724

This may be the easiest...

Also you will note that the index order for the array is reversed. Fortran in fastest in the left, so it is opposite of C.

When you call it it uses which ever one matches the array dimensions that you sent to it automagically, by checking the interface of the module.

MODULE Stuff
IMPLICIT NONE
PUBLIC MY$Work
INTERFACE MY$Work
  MODULE PROCEDURE MY$Work1D, MY$Work2D, MY$Work3D
END INTERFACE
PRIVATE 

CONTAINS   

!~~~~~~~~~~
SUBROUTINE MY$Work1d(Array, F_of_Array)
REAL, DIMENSION(:), INTENT(IN   ) :: Array
REAL, DIMENSION(:), INTENT(INOUT) :: F_of_Array
INTEGER                           :: I

!DIR$ DO SIMD
do i = LBOUND(Array,1) ,UBOUND(Array,1)
  ! do stuff with f(i)
enddo

RETURN
END SUBROUTINE MY$Work1d

!~~~~~~~~~~
SUBROUTINE MY$Work2d(Array, F_of_Array)
REAL, DIMENSION(:,:), INTENT(IN   ) :: Array
REAL, DIMENSION(:,:), INTENT(INOUT) :: F_of_Array
INTEGER                             :: I,J

!DEC$ UNROLL_AND_JAM
 do j = LBOUND(Array,2) ,UBOUND(Array,2)
   do i = LBOUND(Array,1) ,UBOUND(Array,1)
      ! do stuff with f(i,j)
   enddo
enddo

RETURN
END SUBROUTINE MY$Work2d

!~~~~~~~~~~
SUBROUTINE MY$Work3d(Array, F_of_Array)
REAL, DIMENSION(:,:,:), INTENT(IN   ) :: Array
REAL, DIMENSION(:,:,:), INTENT(INOUT) :: F_of_Array
INTEGER                               :: I,J,K

!DEC$ UNROLL_AND_JAM
do k = LBOUND(Array,3) , UBOUND(Array,3)
   do j = LBOUND(Array,2) ,UBOUND(Array,2)
      do i = LBOUND(Array,1) ,UBOUND(Array,1)
      ! do stuff with f(i,j,k)
      enddo
   enddo
enddo

RETURN
END SUBROUTINE MY$Work3d(Array, F_of_Array)

END MODULE Stuff


PROGRAM XX
USE STUFF
...

So you could have a 3D array and do this

DO K = 1, nk
  CALL MY$WORK(Array(:,:,k), ResultArray(:,:,K) )
ENDDO

That would use the 2D version and loop over it for all the K-cases...

Upvotes: 1

Related Questions