!---------------------------------------------------
!
!  Cluster analysis of a 3d array
!  FoF of data in matrix A(N,N,N)
!  Returns iCl(N,N,N) - index of cluster
!             = 0 if cell is not in a cluster
!             = m if cell belongs to cluster number m
!      Nclusters - total number of clusters
!      iMap(1:Nclusters) - number of points in each cluster
!---------------------------------------------------
Module Clusters

integer*4, allocatable, dimension(:,:,:) :: A,iCl
integer*4, allocatable, dimension(:)     :: iMap,jMap
integer*4 :: N
Contains
!---------------------------------------------------
!    Main routine of FoF
!
SUBROUTINE MakeClusters(nClusters)

   Call InitCl(nClusters)
   Do
      Call IterateCl
      Call ClusterCount(nClustersnew)
      if(nClusters==nClustersnew)exit ! no change in clustering
      nClusters =nClustersnew
   EndDo

end SUBROUTINE MakeClusters
!---------------------------------------------------
!             Count number of clusters
!             jMap(i) = # of iCl cells with given iCl value i
!             iMap
SUBROUTINE ClusterStat(nClusters)

    maxM = MAXVAL(iCl)
    allocate (jMap(maxM))
    jMap(:) = 0
    Do k=1,N
       Do j=1,N
          Do i=1,N
                If(iCl(i,j,k)/=0)jMap(iCl(i,j,k)) = jMap(iCl(i,j,k)) +1
          end Do
       end Do
    end Do
    nClusters2 = 0
    jmax      = MAXVAL(iCl)
    Do i=1,MaxM
       If(jMap(i)/=0)nclusters2 =nclusters2 +1
    endDo
    If(nClusters2 /= nClusters)Stop 'Error in cluster numbers'
    imax = MAXVAL(jMap)
    imin = MINVAL(jMap)
    write(*,'(3(a,i8/))') ' Number of clusters= ',nClusters, &
                          ' Cells in smallest = ',imin, &
                          ' Cells in largest  = ',imax
    write(*,*) ' Cluster       Ncells'
    Do i=1,min(MaxM,200),5
       write(*,'(i6,3x,10i7)') i,(jMap(j),j=i,min(i+4,MaxM))
    end Do
    deallocate (jMap)
end SUBROUTINE ClusterStat
!---------------------------------------------------
!             Count number of clusters
!             jMap(i) = # of iCl cells with given iCl value i
!             iMap
SUBROUTINE ClusterCount(nClusters)

    maxM = MAXVAL(iCl)
    allocate (jMap(maxM))
    jMap(:) = 0
    Do k=1,N
       Do j=1,N
          Do i=1,N
                If(iCl(i,j,k)/=0)jMap(iCl(i,j,k)) = jMap(iCl(i,j,k)) +1
          end Do
       end Do
    end Do
    nClusters = 0
    jmax      = MAXVAL(iCl)
    Do i=1,MaxM
       If(jMap(i)/=0)nclusters =nclusters +1
    endDo

    allocate (iMap(jmax))
    iMap(:) = 0
 
    write(*,*) ' Max iCl =',maxM,jmax
    write(*,*) ' N clusters= ',nClusters
    ii =0
    Do i=1,MaxM             ! find cluster index iMap
       If(jMap(i)/=0)Then
           ii = ii +1
           iMap(i) = ii
          ! write(*,'(a,3i10)') '  Npoints= ',jMap(i),iMap(i),i
        end If
    endDo
    !                    assign in new index iMap 
    !                    for old cluster iCl
    !                    New indexes do not have gaps     
!$OMP PARALLEL DO DEFAULT(SHARED)   &                                    
!$OMP PRIVATE ( i,j,k )       
    Do k=1,N
       Do j=1,N
          Do i=1,N
                If(iCl(i,j,k)/=0)Then
                  iCl(i,j,k) =iMap(iCl(i,j,k))
                EndIf
          end Do
       end Do
    end Do
                      ! output for testing
    Do k=1,6
       write(*,*) ' k=',k
     Do j=1,20
        write(*,'(20i6)') (iCl(i,j,k),i=1,20)
     EndDo
    Enddo
    deallocate(jMap,iMap) 
end SUBROUTINE ClusterCount
!---------------------------------------------------
!             initilaize iCl 
!             every non-zero A(i,j,k) gets iCl= number
SUBROUTINE InitCl(nClusters)

myindex =0
Do k=1,N
   Do j=1,N
      Do i=1,N
         If(A(i,j,k)/=0)Then
            myindex = myindex +1
            If(iCl(i,j,k)/=0)Stop ' Error in Init: cell is already taken'
            iCl(i,j,k) = myindex 
         End If
      end Do
   end Do
end Do
nClusters = myindex
write(*,*) ' Init: total number of occupied cells= ',myindex
end SUBROUTINE InitCl
!---------------------------------------------------
!      select smallest of indexes of neighbors
!              
SUBROUTINE IterateCl

!$OMP PARALLEL DO DEFAULT(SHARED)   &                                    
!$OMP PRIVATE ( i,j,k,kp,km,jp,jm,ip,im,k2,j2,i2,mymin )       
Do k=1,N
   kp = min(k+1,N)
   km = max(k-1,1)
   Do j=1,N
      jp = min(j+1,N)
      jm = max(j-1,1)
      Do i=1,N
         ip = min(i+1,N)
         im = max(i-1,1)
         If(iCl(i,j,k)/=0)Then
            Do k2 =km,kp
            Do j2 =jm,jp
            Do i2 =im,ip
               If(iCl(i2,j2,k2)/=0)Then
                  mymin = min(iCl(i,j,k),iCl(i2,j2,k2))
                  iCl(i,j,k)    = mymin
                  iCl(i2,j2,k2) = mymin
               End If
            End Do
            End Do
            End Do
         End If
      end Do
   end Do
end Do

end SUBROUTINE IterateCl

end Module Clusters
!---------------------------------------------------
!
Program ClusterAnalysis
use Clusters

N = 352   ! size of matrix
ALLOCATE (A(N,N,N),iCl(N,N,N))
   A(:,:,:)   = 0
   iCL(:,:,:) = 0

!Do k=9,N,5  ! test
!   Do j=1,N,5
!      Do i=1,N,5
!         A(i,j,k) = 1
!         A(i+1,j,k) = 1
!         A(i,j+1,k) = 1
!         A(i+1,j+1,k) = 1
!         A(i,j,k+1) = 1
!         A(i+1,j,k+1) = 1
!         A(i,j+1,k+1) = 1
!         A(i+1,j+1,k+1) = 1
!      end Do
!   end Do
!end Do

!          test
!A(1,1,1) = 1; A(2,1,1) = 1; A(1,2,1) = 1; A(3,1,1) = 1;A(4,1,1) = 1;A(5,1,1) = 1;
!A(1,1,5) = 1; A(2,1,5) = 1; A(1,2,5) = 1; A(3,1,5) = 1; A(2,2,5) = 1; 
!A(2,3,5) = 1;A(2,4,5) = 1;A(3,4,5) = 1;A(4,4,5) = 1;A(4,3,5) = 1;
!A(3,3,3) = 1

!          test
A(1,1,1) = 1; A(2,1,1) = 1; A(1,2,1) = 1; A(3,1,1) = 1;A(4,1,1) = 1;A(5,1,1) = 1;
A(1,1,5) = 1; A(2,1,5) = 1; A(1,2,5) = 1; A(3,1,5) = 1; A(2,2,5) = 1; 
A(2,3,5) = 1;A(2,4,5) = 1;A(3,4,5) = 1;A(4,4,5) = 1;A(5,3,5) = 1;A(5,2,5) = 1;
A(3,3,3) = 1
A(:,:,107) =100.
A(:,:,207) =100.
A(:,:,307) =100.
A(:,:,308) =100.


Call MakeClusters(nClusters)

If(nClusters==0)Stop ' No clusters are found'

Call ClusterStat(nClusters)


end Program ClusterAnalysis