Module ARRAYS
integer, parameter :: N=5e7, Nbin =8
Real*8, DIMENSION(9,N) ::   X, Y
Integer, DIMENSION(N)  :: LL, ListLev, Level, Level2
Integer L(0:Nbin)

end module ARRAYS
!-----------------------------------------------------------------
PROGRAM MEMTEST
use ARRAYS

!$OMP PARALLEL DO DEFAULT(SHARED) 
do i=1,N
   Y(:,i) =0.
   do j=1,9
      X(j,i) = sin(float(i)/n)+j
   enddo
enddo
Do i=0,Nbin
   !L(i) =N*(1.-1./2**i)
   L(i) =N*i/Nbin
EndDo
write (*,'(a,9i10)')' L=',L
CALL SimpleADD
CALL SHUFFLE
CALL SortedADD

      stop
end PROGRAM MEMTEST

!-----------------------------------------------------------------
SUBROUTINE SimpleADD
use ARRAYS

call cpu_time(t0)
call system_clock(i1time,iratetime)
do k=1,Nbin
   do j=1,4 !2**k
!$OMP PARALLEL DO DEFAULT(SHARED) 
     do i=L(k-1)+1,L(k)
       CALL ADDss(k,i)
     enddo
   enddo
enddo
call cpu_time(t1)
call system_clock(i2time,iratetime)

write(*,'(a,12g14.6)')' time (cpu, wall)=',t1-t0,float(i2time-i1time)/ iratetime


end SUBROUTINE SimpleADD
!----------------------------------
SUBROUTINE SortedADD
use ARRAYS

call cpu_time(t0)
call system_clock(i1time,iratetime)
do k=1,Nbin
   CALL SelectCells(nLev,k)
   do j=1,4 !2**k
!$OMP PARALLEL DO DEFAULT(SHARED) 
     do i=1,nLev
        ic = ListLev(i) 
       CALL ADDss(k,ic)
     enddo
   enddo
enddo
call cpu_time(t1)
call system_clock(i2time,iratetime)

write(*,'(a,12g14.6)')' time (cpu, wall)=',t1-t0,float(i2time-i1time)/ iratetime

end SUBROUTINE SortedADD
!----------------------------------
SUBROUTINE ADDss(k,i)
use ARRAYS
real*8 ww(9)
integer*8 ii

ww =0.
do m =1,9
do j =1,9
    ii = i
    ii = 12701*(ii+m*9+j)+5407
    i0 =mod(ii,N)+1
   ww(m)  =ww(m) + X(j,i)*j*m
enddo 
enddo 
Y(:,i) = ww

end SUBROUTINE ADDss

!----------------------------------
Integer Function Im(i,m,ii)
integer*8 jj,kk
kk = ii
jj= kk*i+5039
Im = MOD(jj,m)
i  = Im

end Function Im
!----------------------------------
SUBROUTINE SHUFFLE
use ARRAYS

LL =0
ii = N/10*4+1
write(*,*) ' ii=',ii
iseed =121071
nrep =0
klast =1
do i=1,N
   j =Im(iseed,N,ii)+1
   If(j==0.or. j>N)Then
          write(*,*) 'error j->',i,j
          stop
   Else
      If(LL(j) /= 0)Then
          Do k =klast,N
            if(LL(k)==0)Then
              j = k
              exit
            endIf
          EndDo
          klast =k
          nrep =nrep +1
      EndIf
      LL(j) =i
   EndIf
endDo
if(nrep/=0)write(*,*) ' Nrepeat=',nrep,klast
do j=1,N
  if(LL(j) == 0)write(*,*) ' gap=',j
Enddo
!$OMP PARALLEL DO DEFAULT(SHARED) 
Do k=1,Nbin
     Do i=L(k-1)+1,L(k)
          Level(i) = k
     EndDo       
EndDo
!$OMP PARALLEL DO DEFAULT(SHARED) 
Do i=1,N
        Y(:,LL(i)) = X(:,i)
         Level2(LL(i)) =Level(i)
EndDo
!$OMP PARALLEL DO DEFAULT(SHARED) 
Do i=1,N
        X(:,i) = Y(:,i)
        Level(i) =Level2(i)
        Y(:,i) =0.
EndDo

end SUBROUTINE SHUFFLE
!-------------------------------------------------------
SUBROUTINE SelectCells( nLev,k)
use ARRAYS

nLev  =0
Do i=1,N
    If(Level(i)==k)Then
             nLev =nLev +1
             ListLev(nLev) =i
    EndIf
EndDo


end SUBROUTINE SelectCells