!--------------------------------------------------------------------
      Module numerics
!
      Integer, Parameter  ::  sp = kind(1.0  )
      Integer, Parameter  ::  dp = kind(1.0D0)
!
      Integer, Parameter  ::  rp = dp
!
      End Module numerics
!--------------------------------------------------------------------
      PROGRAM fake_adi
!
! Guy.Moebs@crihan.fr
! exemple tire de la documentation SGI
!
! Probleme de cache et de TLB
!
!==========================
! Declaration des variables
!==========================
!
      use numerics
!
      implicit none
!
      integer, parameter  ::  ldx = 128
      integer, parameter  ::  ldy = 128
      integer, parameter  ::  ldz = 128
      integer, parameter  ::  nx  = 128
      integer, parameter  ::  ny  = 128
      integer, parameter  ::  nz  = 128
      integer, parameter  ::  maxsteps = 2
!
      real(rp), dimension(ldx, ldy, ldz)  ::  data
!
      integer   ::  i, j, k, istep
      real(rp)  ::  checksum
      real(sp), external  ::  dtime
      real(rp), external  ::  rand
      real(sp)  ::  t, t2(2)
!
!===================
! Debut du programme
!===================
!
      do k = 1, nz
         do j = 1, ny
            do i = 1, nx
               data(i,j,k) = rand()
            end do
         end do
      end do
!
      t = dtime(t2)
!
      do istep = 1, maxsteps
!
!*$* assert concurrent call
         do k = 1, nz
            do j = 1, ny
               call xsweep(data(1,j,k),1,nx)
            end do
         end do
!
!*$* assert concurrent call
         do k = 1, nz
            do i = 1, nx
               call ysweep(data(i,1,k),ldx,ny)
            end do
         end do
!
!*$* assert concurrent call
         do j = 1, ny
            do i = 1, nx
               call zsweep(data(i,j,1),ldx*ldy,nz)
            end do
         end do
!
      end do
!
      t = dtime(t2)
      write(6,100) t
  100 format(1x,'Time:     ',f6.3,' seconds')
!
      checksum = 0.0_rp
      do k = 1, nz
         do j = 1, ny
            do i = 1, nx
               checksum = checksum + data(i,j,k)
            end do
         end do
      end do
!
      write(6,200) checksum
  200 format(1x,'Checksum: ',1pe17.10)
!
      end program fake_adi
!---------------------------------------------------------------------
      subroutine xsweep(v,is,n)
!
      use numerics
      implicit none
      real(rp), dimension(1+is*(n-1)) ::  v
      integer  ::  is, n
      integer  ::  i
      real(rp), parameter  ::  half = 0.5_rp
!
      do i = 2, n
         v(1+is*(i-1)) = v(1+is*(i-1)) + half*v(1+is*(i-2))
      enddo
!
      do i = n-1, 1, -1
         v(1+is*(i-1)) = v(1+is*(i-1)) - half*v(1+is*i)
      enddo
!
      return
      end subroutine xsweep
!---------------------------------------------------------------------
      subroutine ysweep(v,is,n)
!
      use numerics
      implicit none
      real(rp), dimension(1+is*(n-1)) ::  v
      integer  ::  is, n
      integer  ::  i
      real(rp), parameter  ::  half = 0.5_rp
!
      do i = 2, n
         v(1+is*(i-1)) = v(1+is*(i-1)) + half*v(1+is*(i-2))
      end do
!
      do i = n-1, 1, -1
         v(1+is*(i-1)) = v(1+is*(i-1)) - half*v(1+is*i)
      end do
!
      return
      end subroutine ysweep
!---------------------------------------------------------------------
      subroutine zsweep(v,is,n)
!
      use numerics
      implicit none
      real(rp), dimension(1+is*(n-1)) ::  v
      integer  ::  is, n
      integer  ::  i
      real(rp), parameter  ::  half = 0.5_rp
!
      do i = 2, n
         v(1+is*(i-1)) = v(1+is*(i-1)) + half*v(1+is*(i-2))
      end do
!
      do i = n-1, 1, -1
         v(1+is*(i-1)) = v(1+is*(i-1)) - half*v(1+is*i)
      end do
!
      return
      end subroutine zsweep
