      program DISTRINDIRECT3
!    Testing DISTRIBUTE and REDISTRIBUTE directives
!            INDIRECT  distribution
	print *,'=== START OF distrindirect3 ========================'
      call distrindirect31
      print *,'=== END OF distrindirect3 ========================= '    
      end

        subroutine distrindirect31
        parameter    (L=10,  ER=100000)
        integer::     A(L,L,L), B(L,L,L),AS(L,L,L), BS(L,L,L)
        integer,dimension(:,:,:),allocatable:: ib1, ib2, ib3, ib4, ib5, ib6
        integer,dimension(L,L,L):: indir_x, indir_y, indir_z
        integer MAP1(L), MAP2(L), MAP3(L)
        integer:: erri=ER
        character*15:: tname="distrindirect31"
!DVM$   TEMPLATE E(L,L,L)
!DVM$   DISTRIBUTE :: E
!DVM$   ALIGN :: A,B
!DVM$   ALIGN :: indir_x, indir_y,indir_z
!DVM$   ALIGN :: ib1,ib2,ib3,ib4,ib5,ib6

      call distrindirect31_s (AS, BS)
      call fillMap(MAP1,L,1)
      call fillMap(MAP2,L,2)
      call fillMap(MAP3,L,3)
      allocate( ib1(L,L,L), ib2(L,L,L), ib3(L,L,L),  &
     &          ib4(L,L,L), ib5(L,L,L), ib6(L,L,L) )

!DVM$ REDISTRIBUTE E(INDIRECT(MAP1),INDIRECT(MAP2),INDIRECT(MAP3))
!DVM$ REALIGN (I,J,K)  WITH  E(I,J,K) :: A,B
!DVM$ REALIGN (I,J,K)  WITH  E(I,J,K) :: indir_x, indir_y,indir_z
!DVM$ REALIGN (I,J,K)  WITH  E(I,J,K) :: ib1,ib2,ib3,ib4,ib5,ib6
      do i = 1,L
        do j = 1,L
          do k = 1,L

            indir_x(i,j,k) = i
            indir_y(i,j,k) = j
            indir_z(i,j,k) = k

            if (i.gt.1) then
              ib1(i,j,k) = i - 1
            else
              ib1(i,j,k) = 0
            endif
            if (i.lt.L) then
              ib2(i,j,k) = i + 1
            else
              ib2(i,j,k) = 0
            endif
            if (j.gt.1) then
              ib3(i,j,k) = j - 1
            else
              ib3(i,j,k) = 0
            endif
            if (j.lt.L) then
              ib4(i,j,k) = j + 1
            else
              ib4(i,j,k) = 0
            endif
            if (k.gt.1) then
              ib5(i,j,k) = k - 1
            else
              ib5(i,j,k) = 0
            endif
            if (k.lt.L) then
              ib6(i,j,k) = k + 1
            else
              ib6(i,j,k) = 0
            endif
          enddo
        enddo
      enddo

!DVM$ SHADOW_ADD (E((ib1(i,j,k)) with E(@i,@j,@k),*,*) = "nei1") include_to A
!DVM$ SHADOW_ADD (E((ib2(i,j,k)) with E(@i,@j,@k),*,*) = "nei2") include_to A
!DVM$ SHADOW_ADD (E(*,(ib3(i,j,k)) with E(@i,@j,@k),*) = "nei3") include_to A
!DVM$ SHADOW_ADD (E(*,(ib4(i,j,k)) with E(@i,@j,@k),*) = "nei4") include_to A
!DVM$ SHADOW_ADD (E(*,*,(ib5(i,j,k)) with E(@i,@j,@k)) = "nei5") include_to A
!DVM$ SHADOW_ADD (E(*,*,(ib6(i,j,k)) with E(@i,@j,@k)) = "nei6") include_to A

!DVM$ LOCALIZE(ib1 => A(:,*,*))
!DVM$ LOCALIZE(ib2 => A(:,*,*))
!DVM$ LOCALIZE(ib3 => A(*,:,*))
!DVM$ LOCALIZE(ib4 => A(*,:,*))
!DVM$ LOCALIZE(ib5 => A(*,*,:))
!DVM$ LOCALIZE(ib6 => A(*,*,:))

!DVM$ REGION

!DVM$ PARALLEL (k,j,i) ON B(i,j,k)
      do k = 1,L
        do j = 1,L
          do i = 1,L

            A(i,j,k) = 0

            if (indir_x(i,j,k) == 1 .or. indir_x(i,j,k) == L .or.   &
     &          indir_y(i,j,k) == 1 .or. indir_y(i,j,k) == L .or.   &
     &          indir_z(i,j,k) == 1 .or. indir_z(i,j,k) == L) then
              B(i,j,k) = 0
            else
              B(i,j,k) = 1 + indir_x(i,j,k) + indir_y(i,j,k) + indir_z(i,j,k)
            endif

          enddo
        enddo
      enddo

!DVM$ PARALLEL (k,j,i) ON B(i,j,k), SHADOW_RENEW (A)
      do k = 2,L-1
      do j = 2,L-1
      do i = 2,L-1
            if (indir_x(i,j,k) /= 1 .and. indir_x(i,j,k) /= L .and.    &
     &          indir_y(i,j,k) /= 1 .and. indir_y(i,j,k) /= L .and.    &
     &          indir_z(i,j,k) /= 1 .and. indir_z(i,j,k) /= L) then
 
        B(i,j,k) = (A(ib1(i,j,k),j,k) + A(ib2(i,j,k),j,k) +            &
     & A(i,ib3(i,j,k),k) + A(i,ib4(i,j,k),k) + A(i,j,ib5(i,j,k)) +     &
     & A(i,j,ib6(i,j,k))) / 6.0
        endif
      enddo
      enddo
      enddo
!DVM$ PARALLEL (k,j,i) ON B(i,j,k), REDUCTION(min(erri))
      do k = 2,L-1
      do j = 2,L-1
      do i = 2,L-1
             if (indir_x(i,j,k) /= 1 .and. indir_x(i,j,k) /= L .and.    &
     &           indir_y(i,j,k) /= 1 .and. indir_y(i,j,k) /= L .and.    &
     &           indir_z(i,j,k) /= 1 .and. indir_z(i,j,k) /= L) then

            if(B(i,j,k) .ne. BS(i,j,k)) erri =  min(erri, ABS(B(i,j,k)-BS(i,j,k)))
        endif
      enddo
      enddo
      enddo

!DVM$ END REGION
!DVM$ GET_ACTUAL(erri)
      if (erri .eq. ER) then     
          call ansyes(tname)
      else
          call ansno(tname)
      endif
      
      deallocate (ib1,ib2,ib3,ib4,ib5,ib6)
      end subroutine

!---------------------------------------------------------------
      subroutine fillMap(MAP,L,dim)
      integer  numproc 
      integer i,L,dim
      real x
      integer MAP(L)      
      PROCESSORS_SIZE(i) = 1
      numproc = PROCESSORS_SIZE(dim) ! dvmh_get_num_procs(1) 
        do i=1,L
           call RANDOM_NUMBER(x) 
           MAP(i) = MOD(INT(x*10), numproc)     !rand()
        enddo
      end subroutine    
!---------------------------------------------------------------
        subroutine distrindirect31_s (A,B)
        parameter    (L=10)
        integer::     A(L,L,L), B(L,L,L)
        integer,dimension(:,:,:),allocatable:: ib1, ib2, ib3, ib4, ib5, ib6
        integer,dimension(L,L,L):: indir_x, indir_y, indir_z
      allocate( ib1(L,L,L), ib2(L,L,L), ib3(L,L,L),  &
     &          ib4(L,L,L), ib5(L,L,L), ib6(L,L,L) )

      do i = 1,L
        do j = 1,L
          do k = 1,L

            indir_x(i,j,k) = i
            indir_y(i,j,k) = j
            indir_z(i,j,k) = k

            if (i.gt.1) then
              ib1(i,j,k) = i - 1
            else
              ib1(i,j,k) = 0
            endif
            if (i.lt.L) then
              ib2(i,j,k) = i + 1
            else
              ib2(i,j,k) = 0
            endif
            if (j.gt.1) then
              ib3(i,j,k) = j - 1
            else
              ib3(i,j,k) = 0
            endif
            if (j.lt.L) then
              ib4(i,j,k) = j + 1
            else
              ib4(i,j,k) = 0
            endif
            if (k.gt.1) then
              ib5(i,j,k) = k - 1
            else
              ib5(i,j,k) = 0
            endif
            if (k.lt.L) then
              ib6(i,j,k) = k + 1
            else
              ib6(i,j,k) = 0
            endif
          enddo
        enddo
      enddo

      do k = 1,L
        do j = 1,L
          do i = 1,L

            A(i,j,k) = 0

            if (indir_x(i,j,k) == 1 .or. indir_x(i,j,k) == L .or.   &
     &          indir_y(i,j,k) == 1 .or. indir_y(i,j,k) == L .or.   &
     &          indir_z(i,j,k) == 1 .or. indir_z(i,j,k) == L) then
              B(i,j,k) = 0
            else
              B(i,j,k) = 1 + indir_x(i,j,k) + indir_y(i,j,k) + indir_z(i,j,k)
            endif

          enddo
        enddo
      enddo

      do k = 2,L-1
      do j = 2,L-1
      do i = 2,L-1
            if (indir_x(i,j,k) /= 1 .and. indir_x(i,j,k) /= L .and.    &
     &          indir_y(i,j,k) /= 1 .and. indir_y(i,j,k) /= L .and.    &
     &          indir_z(i,j,k) /= 1 .and. indir_z(i,j,k) /= L) then
 
        B(i,j,k) = (A(ib1(i,j,k),j,k) + A(ib2(i,j,k),j,k) +            &
     & A(i,ib3(i,j,k),k) + A(i,ib4(i,j,k),k) + A(i,j,ib5(i,j,k)) +     &
     & A(i,j,ib6(i,j,k))) / 6.0
        endif
      enddo
      enddo
      enddo

      deallocate (ib1,ib2,ib3,ib4,ib5,ib6)
      end subroutine



!---------------------------------------------------------------
      subroutine ansyes(name)
      character*14 name
      print *,name,'  -  complete'
      end
      subroutine ansno(name)
      character*9 name
      print *,name,'  -  ***error'
      end

