!  Demonstration programs for scans (parallel prefix operations).
!  These programs are intended only to illustrate how various operations
!  may be performed in O(log(N)) parallel steps; they should not be used
!  in actual practice.
!
!  Source file:  scans.f90       Parameters:  scans.hp
!  Common data:  scans.hc        Interfaces:  scans.hI
!
!  For convenience, all operations are implemented along in
!  an M-by-N array.

!  Parameters and common items in scans.h, scans.hc :
!     M, N - array dimensions, each is a power of 2
!     Him = M/2,      Hin = N/2     (high-order coordinate bit)
!     Km  = log2(M),  Kn  = log2(N)
!     Mkm - M - 1,    Mkn = N - 1   (mask for bit processing)
!     MaxId - the most negative integer (identity)
!     Prn - controls printing (true for progress reports)
!
!     Ix(0:M-1), Jx(0:N-1)               1-D coordinate arrays for subscripts
!     I(0:M-1, 0:N-1), J(0:M-1, 0:N-1)   2-D coordinate arrays

!  Initialize scan system.  Call this routine (once) before performing
!  any scans.

      subroutine scaninit()
      INCLUDE 'scans.hp'
      INCLUDE 'scans.hc'

      integer ii

      Ix = (/ (ii, ii=0,M-1) /);  Jx = (/ (ii, ii=0,N-1) /)
      I  = SPREAD(Ix, 2, N);      J  = SPREAD(Jx, 1, M)
      MaxId = -2147473648           ! The most negative integer
      print *, "Initialized for ", M, " by ", N, &
               " arrays with identity ", MaxId
      end subroutine scaninit


!  Print integer array - intended for small arrays, small values.
!  Chops extreme negative values

      subroutine pparray(X)
      INCLUDE 'scans.hp'             ! For parameters
      integer :: X(0:M-1, 0:N-1)

      integer k

      do, k = 0, M-1
        print 21, max(X(k,:), -99)
      end do
  21  FORMAT (16(1X, I3))       ! 16 is about as big as we can handle
      end subroutine pparray


!  Scan functions use max, but any associative operation will do.
!  The algorithms are also suitable for operations which are only
!  approximately associative, such as floating-point addition.

! =======================================================================
!
!                    Scans without segmentation

!  Exclusive max scan in axis 1 using 2*Kn steps, without segmentation,
!  Scans either upward or downward.

      subroutine XmaxScan1(Rslt, Src, upwd)
      INCLUDE 'scans.hp'
      integer, dimension(0:M-1, 0:N-1) :: Rslt, Src
      logical upwd        ! .true. for upward scan

      INCLUDE 'scans.hc'

      integer, dimension(0:M-1, 0:N-1) :: V, W, Z, Id
      integer U(0:Km-1, 0:M-1, 0:N-1)
      logical act(0:M-1, 0:N-1)
      integer k, h, mask

      V = Src;  act = .true.
      Id = I;  if ( .not. upwd ) Id = not(Id)     ! Logical address
    ! Reduce to last position.
      h = 1       ! h = 2**k
      do, k = 0, Km-1
        W = MaxId
        WHERE ( act .and. iand(Id, h) == 0 )
          W(ieor(Ix, h), :) = V          ! Send to "partner"
          act = .false.                  ! Drop out
        END WHERE
        U(k, :, :) = W;  V = max(W, V);  h = 2*h
        if ( Prn ) then
          print *, "Ph1, k = ", k, "   Uk:";  call pparray(W);
        end if
      end do
    ! Phase 1 is finished
    ! Propagate down the tree from last processor, the only one active
      W = MaxId;  mask = Mkm;  h = Him   ! h = 2^k
      do k = Km-1, 0, -1
        WHERE ( act )
          W(ieor(Ix, h), :) = W          ! Send to partner
          W = max(W, U(k, :, :))
        END WHERE
        if ( Prn ) then
          print *, "Ph2, k, h = ", k, h, "   W:";  call pparray(W)
        end if
        mask = ishft(mask, -1);  h = h/2
        act = (iand(Id, mask) == mask)
      end do
      Rslt = W                           ! Report result
      end subroutine XmaxScan1


!  Inclusive maximum scans in axis 1.  Arguments are like XmaxScan1

      subroutine ImaxScan1(Rslt, Src, upwd)
      INCLUDE 'scans.hp'
      integer, dimension(0:M-1, 0:N-1) :: Rslt, Src
      logical upwd        ! .true. for upward scan

      INTERFACE
      subroutine XmaxScan1(Rslt, Src, upwd)
      INCLUDE 'scans.hp'
      integer, dimension(0:M-1, 0:N-1) :: Rslt, Src
      logical upwd        ! .true. for upward scan
      END INTERFACE

      integer W(0:M-1, 0:N-1)

      call XmaxScan1(W, Src, upwd)
      Rslt = max(W, Src)
      end subroutine ImaxScan1


!  Exclusive max scan in axis 2 using 2*Kn steps, without segmentation,
!  Scans either upward or downward.

      subroutine XmaxScan2(Rslt, Src, upwd)
      INCLUDE 'scans.hp'
      integer, dimension(0:M-1, 0:N-1) :: Rslt, Src
      logical upwd        ! .true. for upward scan

      INCLUDE 'scans.hc'

      integer, dimension(0:M-1, 0:N-1) :: V, W, Z, Jd
      integer U(0:Kn-1, 0:M-1, 0:N-1)
      logical act(0:M-1, 0:N-1)
      integer k, h, mask

      V = Src;  act = .true.
      Jd = J;  if ( .not. upwd ) Jd = not(Jd)     ! Logical address
    ! Reduce to last position.
      h = 1       ! h = 2**k
      do, k = 0, Kn-1
        W = MaxId
        WHERE ( act .and. iand(Jd, h) == 0 )
          W(:, ieor(Jx, h)) = V          ! Send to "partner"
          act = .false.                  ! Drop out
        END WHERE
        U(k, :, :) = W;  V = max(W, V);  h = 2*h
        if ( Prn ) then
          print *, "Ph1, k = ", k, "   Uk:";  call pparray(W);
        end if
      end do
    ! Phase 1 is finished
    ! Propagate down the tree from last processor, the only one active
      W = MaxId;  mask = N-1;  h = N/2   ! h = 2^k
      do k = Kn-1, 0, -1
        WHERE ( act )
          W(:, ieor(Jx, h)) = W          ! Send to partner
          W = max(W, U(k, :, :))
        END WHERE
        if ( Prn ) then
          print *, "Ph2, k, h = ", k, h, "   W:";  call pparray(W)
        end if
        mask = ishft(mask, -1);  h = h/2
        act = (iand(Jd, mask) == mask)
      end do
      Rslt = W                           ! Report result
      end subroutine XmaxScan2


!  Inclusive maximum scans in axis 2.  Arguments are like XmaxScan2

      subroutine ImaxScan2(Rslt, Src, upwd)
      INCLUDE 'scans.hp'
      integer, dimension(0:M-1, 0:N-1) :: Rslt, Src
      logical upwd        ! .true. for upward scan

      INTERFACE
      subroutine XmaxScan2(Rslt, Src, upwd)
      INCLUDE 'scans.hp'
      integer, dimension(0:M-1, 0:N-1) :: Rslt, Src
      logical upwd        ! .true. for upward scan
      END INTERFACE

      integer W(0:M-1, 0:N-1)

      call XmaxScan2(W, Src, upwd)
      Rslt = max(W, Src)
      end subroutine ImaxScan2

! ==========================================================================
!
!                       Scans with segmentation


!  Exclusive max scan in axis 1 using 2*Kn steps, with segmentation.
!  Scans either upward or downward.

      subroutine XmaxSeg1(Rslt, Src, Sgmt, upwd)
      INCLUDE 'scans.hp'
      integer, dimension(0:M-1, 0:N-1) :: Rslt, Src
      logical, dimension(0:M-1, 0:N-1) :: Sgmt       ! Segmentation flags
      logical upwd        ! .true. for upward scan

      INCLUDE 'scans.hc'

      integer, dimension(0:M-1, 0:N-1) :: V,  W,  Z, Id
      logical, dimension(0:M-1, 0:N-1) :: Vf, Wf, Zf
      integer U(0:Km-1, 0:M-1, 0:N-1)
      logical Uf(0:Km-1, 0:M-1, 0:N-1)
      logical act(0:M-1, 0:N-1)
      integer k, h, mask

      V = Src;  Vf = Sgmt;  act = .true.
      Id = I;  if ( .not. upwd ) Id = not(Id)     ! Logical address
    ! Reduce to last position.
      h = 1       ! h = 2**k
      do, k = 0, Km-1
        W = MaxId;  Wf = .false.
        WHERE ( act .and. iand(Id, h) == 0 )
           W(ieor(Ix, h), :) = V         ! Send to "partner"
          Wf(ieor(Ix, h), :) = Vf
          act = .false.                  ! Drop out
        END WHERE
        U(k, :, :) = W;  Uf(k, :, :) = Wf;  h = 2*h
        WHERE ( .not. Vf )          ! <Vf,V> := <Wf,W> max' <Vf,V>
          V = max(W, V)
        END WHERE
        Vf = Wf .or. Vf
        if ( Prn ) then
          print *, "Ph1, k = ", k, "   Uk:";  call pparray(W);
        end if
      end do
    ! Phase 1 is finished
    ! Propagate down the tree from last processor, the only one active
      W = MaxId;  Wf = .false.;  mask = Mkm;  h = Him   ! h = 2^k
      do k = Km-1, 0, -1
        WHERE ( act )
           W(ieor(Ix, h), :) = W          ! Send to partner
          Wf(ieor(Ix, h), :) = Wf
          ! <Wf,W> := <Wf,W> max' <Ufk,Uk>
          W = MERGE( U(k, :, :), max(W, U(k, :, :)), Uf(k, :, :) )
          Wf = Wf .or. Uf(k, :, :)
        END WHERE
        if ( Prn ) then
          print *, "Ph2, k, h = ", k, h, "   W:";  call pparray(W)
        end if
        mask = ishft(mask, -1);  h = h/2
        act = (iand(Id, mask) == mask)
      end do
      Rslt = W                           ! Report result (value)
      end subroutine XmaxSeg1


!  Inclusive segmented scans.  Arguments are like XmaxSeg1 above.

      subroutine ImaxSeg1(Rslt, Src, Sgmt, upwd)
      INCLUDE 'scans.hp'
      integer, dimension(0:M-1, 0:N-1) :: Rslt, Src
      logical, dimension(0:M-1, 0:N-1) :: Sgmt       ! Segmentation flags
      logical upwd        ! .true. for upward scan

      INTERFACE
      subroutine XmaxSeg1(Rslt, Src, Sgmt, upwd)
      INCLUDE 'scans.hp'
      integer, dimension(0:M-1, 0:N-1) :: Rslt, Src
      logical, dimension(0:M-1, 0:N-1) :: Sgmt       ! Segmentation flags
      logical upwd        ! .true. for upward scan
      END INTERFACE

      INCLUDE 'scans.hc'

      integer, dimension(0:M-1, 0:N-1) :: W

      call XmaxSeg1(W, Src, Sgmt, upwd)
      WHERE ( Sgmt ) W = MaxId
      Rslt = max(W, Src)
      end subroutine ImaxSeg1


!  Exclusive max scan in axis 2 using 2*Kn steps, with segmentation.
!  Scans either upward or downward.

      subroutine XmaxSeg2(Rslt, Src, Sgmt, upwd)
      INCLUDE 'scans.hp'
      integer, dimension(0:M-1, 0:N-1) :: Rslt, Src
      logical, dimension(0:M-1, 0:N-1) :: Sgmt       ! Segmentation flags
      logical upwd        ! .true. for upward scan

      INCLUDE 'scans.hc'

      integer, dimension(0:M-1, 0:N-1) :: V,  W,  Z, Jd
      logical, dimension(0:M-1, 0:N-1) :: Vf, Wf, Zf
      integer U(0:Kn-1, 0:M-1, 0:N-1)
      logical Uf(0:Kn-1, 0:M-1, 0:N-1)
      logical act(0:M-1, 0:N-1)
      integer k, h, mask

      V = Src;  Vf = Sgmt;  act = .true.
      if ( Prn ) then
        print *, "XmaxSeg2,  MaxId = ", MaxId
      end if
      Jd = J;  if ( .not. upwd ) Jd = not(Jd)     ! Logical address
    ! Reduce to last position.
      h = 1       ! h = 2**k
      do, k = 0, Kn-1
        W = MaxId;  Wf = .false.
        WHERE ( act .and. iand(Jd, h) == 0 )
           W(:, ieor(Jx, h)) = V         ! Send to "partner"
          Wf(:, ieor(Jx, h)) = Vf
          act = .false.                  ! Drop out
        END WHERE
        U(k, :, :) = W;  Uf(k, :, :) = Wf;  h = 2*h
        WHERE ( .not. Vf )          ! <Vf,V> := <Wf,W> max' <Vf,V>
          V = max(W, V)
        END WHERE
        Vf = Wf .or. Vf
        if ( Prn ) then
          print *, "Ph1, k = ", k, "   Uk:";  call pparray(W);
        end if
      end do
    ! Phase 1 is finished
    ! Propagate down the tree from last processor, the only one active
      W = MaxId;  Wf = .false.;  mask = Mkn;  h = Hin   ! h = 2^k
      do k = Kn-1, 0, -1
        WHERE ( act )
           W(:, ieor(Jx, h)) = W          ! Send to partner
          Wf(:, ieor(Jx, h)) = Wf
          ! <Wf,W> := <Wf,W> max' <Ufk,Uk>
          W = MERGE( U(k, :, :), max(W, U(k, :, :)), Uf(k, :, :) )
          Wf = Wf .or. Uf(k, :, :)
        END WHERE
        if ( Prn ) then
          print *, "Ph2, k, h = ", k, h, "   W:";  call pparray(W)
        end if
        mask = ishft(mask, -1);  h = h/2
        act = (iand(Jd, mask) == mask)
      end do
      Rslt = W                           ! Report result (value)
      end subroutine XmaxSeg2


!  Inclusive segmented scans.  Arguments are like XmaxSeg2 above.

      subroutine ImaxSeg2(Rslt, Src, Sgmt, upwd)
      INCLUDE 'scans.hp'
      integer, dimension(0:M-1, 0:N-1) :: Rslt, Src
      logical, dimension(0:M-1, 0:N-1) :: Sgmt       ! Segmentation flags
      logical upwd        ! .true. for upward scan

      INTERFACE
      subroutine XmaxSeg2(Rslt, Src, Sgmt, upwd)
      INCLUDE 'scans.hp'
      integer, dimension(0:M-1, 0:N-1) :: Rslt, Src
      logical, dimension(0:M-1, 0:N-1) :: Sgmt       ! Segmentation flags
      logical upwd        ! .true. for upward scan
      END INTERFACE

      INCLUDE 'scans.hc'

      integer, dimension(0:M-1, 0:N-1) :: W

      call XmaxSeg2(W, Src, Sgmt, upwd)
      WHERE ( Sgmt ) W = MaxId
      Rslt = max(W, Src)
      end subroutine ImaxSeg2

! ===========================================================================
!
!  General scan routines

!  Unsegmented scans of Src along axis dim,
!                       Inclusive/Exclusive, Upward/Downward

      function MaxScan(Src, dim, ixcl, dirn)
      INCLUDE 'scans.hp'
      integer, dimension(0:M-1, 0:N-1) :: Src, MaxScan
      integer dim, ixcl
      logical dirn

      INCLUDE 'scans.hI'

      integer, dimension(0:M-1, 0:N-1) :: W

      if ( dim == 1 ) then
        if ( ixcl == Exclusive ) then
          call XmaxScan1(W, Src, dirn)
        else
          call ImaxScan1(W, Src, dirn)
        end if
      else
        if ( ixcl == Exclusive ) then
          call XmaxScan2(W, Src, dirn)
        else
          call ImaxScan2(W, Src, dirn)
        end if
      end if
      MaxScan = W
      end function MaxScan

      function MinScan(Src, dim, ixcl, dirn)
      INCLUDE 'scans.hp'
      integer, dimension(0:M-1, 0:N-1) :: Src, MinScan
      integer dim, ixcl
      logical dirn

      INTERFACE
      function MaxScan(Src, dim, ixcl, dirn)
      INCLUDE 'scans.hp'
      integer, dimension(0:M-1, 0:N-1) :: Src, MaxScan
      integer dim, ixcl
      logical dirn
      END INTERFACE

      Minscan = - MaxScan(-Src, dim, ixcl, dirn)
      end function MinScan


!  Segmented scans of Src along axis dim, Segmentation
!                       Inclusive/Exclusive, Upward/Downward

      function SegMaxScan(Src, dim, Sgmt, ixcl, dirn)
      INCLUDE 'scans.hp'
      integer, dimension(0:M-1, 0:N-1) :: Src, SegMaxScan
      logical, dimension(0:M-1, 0:N-1) :: Sgmt
      integer dim, ixcl
      logical dirn
      INCLUDE 'scans.hI'

      integer, dimension(0:M-1, 0:N-1) :: W

      if ( dim == 1 ) then
        if ( ixcl == Exclusive ) then
          call XmaxSeg1(W, Src, Sgmt, dirn)
        else
          call ImaxSeg1(W, Src, Sgmt, dirn)
        end if
      else
        if ( ixcl == Exclusive ) then
          call XmaxSeg2(W, Src, Sgmt, dirn)
        else
          call ImaxSeg2(W, Src, Sgmt, dirn)
        end if
      end if
      SegMaxScan = W
      end function SegMaxScan


      function SegMinScan(Src, dim, Sgmt, ixcl, dirn)
      INCLUDE 'scans.hp'
      integer, dimension(0:M-1, 0:N-1) :: Src, SegMinScan
      logical, dimension(0:M-1, 0:N-1) :: Sgmt
      integer dim, ixcl
      logical dirn
      INCLUDE 'scans.hI'

      INTERFACE
      function SegMaxScan(Src, dim, Sgmt, ixcl, dirn)
      INCLUDE 'scans.hp'
      integer, dimension(0:M-1, 0:N-1) :: Src, SegMaxScan
      logical, dimension(0:M-1, 0:N-1) :: Sgmt
      integer dim, ixcl
      logical dirn
      INCLUDE 'scans.hI'
      END INTERFACE

      SegMinScan = - SegMaxScan(-Src, dim, Sgmt, ixcl, dirn)
      end function SegMinScan
