MODULE Space_Filling_Curves
      USE Precision
      USE Error_Handling
      USE System_Monitors
      USE Sorting_Ranking
      PUBLIC::SFCOrder
      PRIVATE
      
      INTERFACE
      
      FUNCTION fhilbert_c2i(n_dim,n_bits,coords)RESULT(SFC_index)
      USE Precision
      INTEGER,INTENT(IN)::n_dim,n_bits
      INTEGER,DIMENSION(n_dim),INTENT(IN)::coords
      INTEGER(KIND=i_dp)::SFC_index
      END FUNCTION fhilbert_c2i
      
      FUNCTION fmorton_c2i(n_dim,n_bits,coords)RESULT(SFC_index)
      USE Precision
      INTEGER,INTENT(IN)::n_dim,n_bits
      INTEGER,DIMENSION(n_dim),INTENT(IN)::coords
      INTEGER(KIND=i_dp)::SFC_index
      END FUNCTION fmorton_c2i
      END INTERFACE
      
      
      INTERFACE SFCOrder
      MODULE PROCEDURE SFCOrder_i_sp
      MODULE PROCEDURE SFCOrder_i_dp
      MODULE PROCEDURE SFCOrder_r_sp
      MODULE PROCEDURE SFCOrder_r_dp
      END INTERFACE
      
      
      
      CONTAINS
      
      
      
      SUBROUTINE SFCOrder_i_sp(points_coords,SFC_order,&
      SFC,coords_min,coords_max)
      IMPLICIT NONE
      INTEGER(KIND=i_sp),DIMENSION(:,:),INTENT(IN)::points_coords
      INTEGER(KIND=i_wp),DIMENSION(:),INTENT(OUT)::SFC_order
      CHARACTER(LEN=*),INTENT(IN),OPTIONAL::SFC
      INTEGER(KIND=i_sp),DIMENSION(:),INTENT(IN),OPTIONAL::coords_min,coords_max
      
      INTEGER(KIND=i_dp),DIMENSION(:),ALLOCATABLE::SFC_keys
      INTEGER(KIND=i_dp)::max_key
      INTEGER(KIND=i_wp)::n_points,point
      INTEGER::dim,n_dim,n_bits,alloc_status
      INTEGER(KIND=i_32)::max_int_coord
      INTEGER,DIMENSION(SIZE(points_coords,DIM=1))::i_point_coords
      INTEGER(KIND=i_sp),DIMENSION(SIZE(points_coords,DIM=1))::coords_lb,coords_ub,coords_range
      REAL::scale_factor
      LOGICAL::Hilbert_SFC
      
      n_dim=SIZE(points_coords,DIM=1)
      n_points=INT(SIZE(points_coords,DIM=2),KIND=i_wp)
      
      Hilbert_SFC=.TRUE.
      IF(PRESENT(SFC))THEN
      IF((SFC(1:1)=='M').OR.(SFC(1:1)=='m'))Hilbert_SFC=.FALSE.
      END IF
      
      IF(PRESENT(coords_min))THEN
      coords_lb=coords_min
      ELSE
      coords_lb=MINVAL(points_coords,DIM=2)
      END IF
      IF(PRESENT(coords_max))THEN
      coords_ub=coords_max
      ELSE
      coords_ub=MAXVAL(points_coords,DIM=2)
      END IF
      coords_range=coords_ub-coords_lb
      
      ALLOCATE(SFC_keys(n_points),STAT=alloc_status)
      CALL RecordAllocation(n_elements=n_points,mold=1_i_dp,&
      caller="SFCOrder",alloc_status=alloc_status)
      
      n_bits=MIN(INT(BIT_SIZE(1_i_dp)-1)/n_dim,31)
      max_int_coord=IBITS(HUGE(1_i_32),0,n_bits)
      scale_factor=REAL(max_int_coord)/REAL(MAXVAL(coords_range))
      
      DO point=1,n_points
      i_point_coords=&
      INT(scale_factor*REAL(points_coords(1:n_dim,point)-coords_lb),i_wp)
      IF(Hilbert_SFC)THEN
      SFC_keys(point)=fhilbert_c2i(n_dim,n_bits,i_point_coords)
      ELSE
      SFC_keys(point)=fmorton_c2i(n_dim,n_bits,i_point_coords)
      END IF
      END DO
      max_key=2_i_dp**INT(n_bits*n_dim-1)
      CALL QuickRank(array=SFC_keys,permutation=SFC_order,pivot_selection='U',&
      mean_value=max_key,standard_deviation=INT(max_key/SQRT(3.0),KIND=i_dp))
      CALL RecordAllocation(n_elements=-INT(SIZE(SFC_keys),KIND=i_wp),mold=1_i_dp)
      DEALLOCATE(SFC_keys)
      
      END SUBROUTINE
      
      SUBROUTINE SFCOrder_i_dp(points_coords,SFC_order,&
      SFC,coords_min,coords_max)
      IMPLICIT NONE
      INTEGER(KIND=i_dp),DIMENSION(:,:),INTENT(IN)::points_coords
      INTEGER(KIND=i_wp),DIMENSION(:),INTENT(OUT)::SFC_order
      CHARACTER(LEN=*),INTENT(IN),OPTIONAL::SFC
      INTEGER(KIND=i_dp),DIMENSION(:),INTENT(IN),OPTIONAL::coords_min,coords_max
      
      INTEGER(KIND=i_dp),DIMENSION(:),ALLOCATABLE::SFC_keys
      INTEGER(KIND=i_dp)::max_key
      INTEGER(KIND=i_wp)::n_points,point
      INTEGER::dim,n_dim,n_bits,alloc_status
      INTEGER(KIND=i_32)::max_int_coord
      INTEGER,DIMENSION(SIZE(points_coords,DIM=1))::i_point_coords
      INTEGER(KIND=i_dp),DIMENSION(SIZE(points_coords,DIM=1))::coords_lb,coords_ub,coords_range
      REAL::scale_factor
      LOGICAL::Hilbert_SFC
      
      n_dim=SIZE(points_coords,DIM=1)
      n_points=INT(SIZE(points_coords,DIM=2),KIND=i_wp)
      
      Hilbert_SFC=.TRUE.
      IF(PRESENT(SFC))THEN
      IF((SFC(1:1)=='M').OR.(SFC(1:1)=='m'))Hilbert_SFC=.FALSE.
      END IF
      
      IF(PRESENT(coords_min))THEN
      coords_lb=coords_min
      ELSE
      coords_lb=MINVAL(points_coords,DIM=2)
      END IF
      IF(PRESENT(coords_max))THEN
      coords_ub=coords_max
      ELSE
      coords_ub=MAXVAL(points_coords,DIM=2)
      END IF
      coords_range=coords_ub-coords_lb
      
      ALLOCATE(SFC_keys(n_points),STAT=alloc_status)
      CALL RecordAllocation(n_elements=n_points,mold=1_i_dp,&
      caller="SFCOrder",alloc_status=alloc_status)
      
      n_bits=MIN(INT(BIT_SIZE(1_i_dp)-1)/n_dim,31)
      max_int_coord=IBITS(HUGE(1_i_32),0,n_bits)
      scale_factor=REAL(max_int_coord)/REAL(MAXVAL(coords_range))
      
      DO point=1,n_points
      i_point_coords=&
      INT(scale_factor*REAL(points_coords(1:n_dim,point)-coords_lb),i_wp)
      IF(Hilbert_SFC)THEN
      SFC_keys(point)=fhilbert_c2i(n_dim,n_bits,i_point_coords)
      ELSE
      SFC_keys(point)=fmorton_c2i(n_dim,n_bits,i_point_coords)
      END IF
      END DO
      max_key=2_i_dp**INT(n_bits*n_dim-1)
      CALL QuickRank(array=SFC_keys,permutation=SFC_order,pivot_selection='U',&
      mean_value=max_key,standard_deviation=INT(max_key/SQRT(3.0),KIND=i_dp))
      CALL RecordAllocation(n_elements=-INT(SIZE(SFC_keys),KIND=i_wp),mold=1_i_dp)
      DEALLOCATE(SFC_keys)
      
      END SUBROUTINE
      
      SUBROUTINE SFCOrder_r_sp(points_coords,SFC_order,&
      SFC,coords_min,coords_max)
      IMPLICIT NONE
      REAL(KIND=r_sp),DIMENSION(:,:),INTENT(IN)::points_coords
      INTEGER(KIND=i_wp),DIMENSION(:),INTENT(OUT)::SFC_order
      CHARACTER(LEN=*),INTENT(IN),OPTIONAL::SFC
      REAL(KIND=r_sp),DIMENSION(:),INTENT(IN),OPTIONAL::coords_min,coords_max
      
      INTEGER(KIND=i_dp),DIMENSION(:),ALLOCATABLE::SFC_keys
      INTEGER(KIND=i_dp)::max_key
      INTEGER(KIND=i_wp)::n_points,point
      INTEGER::dim,n_dim,n_bits,alloc_status
      INTEGER(KIND=i_32)::max_int_coord
      INTEGER,DIMENSION(SIZE(points_coords,DIM=1))::i_point_coords
      REAL(KIND=r_sp),DIMENSION(SIZE(points_coords,DIM=1))::coords_lb,coords_ub,coords_range
      REAL::scale_factor
      LOGICAL::Hilbert_SFC
      
      n_dim=SIZE(points_coords,DIM=1)
      n_points=INT(SIZE(points_coords,DIM=2),KIND=i_wp)
      
      Hilbert_SFC=.TRUE.
      IF(PRESENT(SFC))THEN
      IF((SFC(1:1)=='M').OR.(SFC(1:1)=='m'))Hilbert_SFC=.FALSE.
      END IF
      
      IF(PRESENT(coords_min))THEN
      coords_lb=coords_min
      ELSE
      coords_lb=MINVAL(points_coords,DIM=2)
      END IF
      IF(PRESENT(coords_max))THEN
      coords_ub=coords_max
      ELSE
      coords_ub=MAXVAL(points_coords,DIM=2)
      END IF
      coords_range=coords_ub-coords_lb
      
      ALLOCATE(SFC_keys(n_points),STAT=alloc_status)
      CALL RecordAllocation(n_elements=n_points,mold=1_i_dp,&
      caller="SFCOrder",alloc_status=alloc_status)
      
      n_bits=MIN(INT(BIT_SIZE(1_i_dp)-1)/n_dim,31)
      max_int_coord=IBITS(HUGE(1_i_32),0,n_bits)
      scale_factor=REAL(max_int_coord)/REAL(MAXVAL(coords_range))
      
      DO point=1,n_points
      i_point_coords=&
      INT(scale_factor*REAL(points_coords(1:n_dim,point)-coords_lb),i_wp)
      IF(Hilbert_SFC)THEN
      SFC_keys(point)=fhilbert_c2i(n_dim,n_bits,i_point_coords)
      ELSE
      SFC_keys(point)=fmorton_c2i(n_dim,n_bits,i_point_coords)
      END IF
      END DO
      max_key=2_i_dp**INT(n_bits*n_dim-1)
      CALL QuickRank(array=SFC_keys,permutation=SFC_order,pivot_selection='U',&
      mean_value=max_key,standard_deviation=INT(max_key/SQRT(3.0),KIND=i_dp))
      CALL RecordAllocation(n_elements=-INT(SIZE(SFC_keys),KIND=i_wp),mold=1_i_dp)
      DEALLOCATE(SFC_keys)
      
      END SUBROUTINE
      
      SUBROUTINE SFCOrder_r_dp(points_coords,SFC_order,&
      SFC,coords_min,coords_max)
      IMPLICIT NONE
      REAL(KIND=r_dp),DIMENSION(:,:),INTENT(IN)::points_coords
      INTEGER(KIND=i_wp),DIMENSION(:),INTENT(OUT)::SFC_order
      CHARACTER(LEN=*),INTENT(IN),OPTIONAL::SFC
      REAL(KIND=r_dp),DIMENSION(:),INTENT(IN),OPTIONAL::coords_min,coords_max
      
      INTEGER(KIND=i_dp),DIMENSION(:),ALLOCATABLE::SFC_keys
      INTEGER(KIND=i_dp)::max_key
      INTEGER(KIND=i_wp)::n_points,point
      INTEGER::dim,n_dim,n_bits,alloc_status
      INTEGER(KIND=i_32)::max_int_coord
      INTEGER,DIMENSION(SIZE(points_coords,DIM=1))::i_point_coords
      REAL(KIND=r_dp),DIMENSION(SIZE(points_coords,DIM=1))::coords_lb,coords_ub,coords_range
      REAL::scale_factor
      LOGICAL::Hilbert_SFC
      
      n_dim=SIZE(points_coords,DIM=1)
      n_points=INT(SIZE(points_coords,DIM=2),KIND=i_wp)
      
      Hilbert_SFC=.TRUE.
      IF(PRESENT(SFC))THEN
      IF((SFC(1:1)=='M').OR.(SFC(1:1)=='m'))Hilbert_SFC=.FALSE.
      END IF
      
      IF(PRESENT(coords_min))THEN
      coords_lb=coords_min
      ELSE
      coords_lb=MINVAL(points_coords,DIM=2)
      END IF
      IF(PRESENT(coords_max))THEN
      coords_ub=coords_max
      ELSE
      coords_ub=MAXVAL(points_coords,DIM=2)
      END IF
      coords_range=coords_ub-coords_lb
      
      ALLOCATE(SFC_keys(n_points),STAT=alloc_status)
      CALL RecordAllocation(n_elements=n_points,mold=1_i_dp,&
      caller="SFCOrder",alloc_status=alloc_status)
      
      n_bits=MIN(INT(BIT_SIZE(1_i_dp)-1)/n_dim,31)
      max_int_coord=IBITS(HUGE(1_i_32),0,n_bits)
      scale_factor=REAL(max_int_coord)/REAL(MAXVAL(coords_range))
      
      DO point=1,n_points
      i_point_coords=&
      INT(scale_factor*REAL(points_coords(1:n_dim,point)-coords_lb),i_wp)
      IF(Hilbert_SFC)THEN
      SFC_keys(point)=fhilbert_c2i(n_dim,n_bits,i_point_coords)
      ELSE
      SFC_keys(point)=fmorton_c2i(n_dim,n_bits,i_point_coords)
      END IF
      END DO
      max_key=2_i_dp**INT(n_bits*n_dim-1)
      CALL QuickRank(array=SFC_keys,permutation=SFC_order,pivot_selection='U',&
      mean_value=max_key,standard_deviation=INT(max_key/SQRT(3.0),KIND=i_dp))
      CALL RecordAllocation(n_elements=-INT(SIZE(SFC_keys),KIND=i_wp),mold=1_i_dp)
      DEALLOCATE(SFC_keys)
      
      END SUBROUTINE
      
      
      
      END MODULE Space_Filling_Curves