ecsort_mix.F90

Go to the documentation of this file.

#ifdef RS6K
@process nocheck
#endif
MODULE ecsort_mix
USE parkind1   , ONLY : jpim     ,jpib     ,jprb     ,jprm    ,jprd
USE yomhook    , ONLY : lhook, dr_hook
USE oml_mod     , ONLY : oml_max_threads, oml_my_thread, oml_in_parallel
USE strhandler_mod , ONLY : toupper
#ifdef SFX_MPI
USE mpl_module , ONLY : mpl_myrank, mpl_nproc
#endif

!$
!..   Author: Sami Saarinen, ECMWF, 10/02/98
!     Fixes : Sami Saarinen, ECMWF, 08/11/99 : Sub-arrays go now correctly (look for addrdiff)
!                                              Genuine real(4) sort "re-habilitated"
!                                              sizeof_int, _real4 & _real8 HARDCODED !
!             Sami Saarinen, ECMWF, 11/10/00 : REAL*4 version included (REAL_M)
!             Sami Saarinen, ECMWF, 28/11/03 : Calls to DR_HOOK added manually (on top of CY28)
!             Sami Saarinen, ECMWF, 18/02/05 : 64-bit integer sorting introduced (for CY30)
!             Sami Saarinen, ECMWF, 22/02/05 : Using genuine 64-bit rsort64() => one-pass through data
!             Sami Saarinen, ECMWF, 06/07/05 : "current_method" made OpenMP-thread aware (for max. # of threads = NTHRDS)
!             Sami Saarinen, ECMWF, 07/07/05 : Quick-sort method finally arrived (and applicable to multikeys, too)
!                                              Quick-sort the default for scalar machines ("non-VPP"), VPPs is radix-sort
!             Sami Saarinen, ECMWF, 03/07/07 : Quick-sort method uses stable approach (was not guaranteed so before)
!             Sami Saarinen, ECMWF, 15/10/07 : Subroutines put into a common file ../include/ecsort_shared.h and
!                                              preprocessed from there
!             Sami Saarinen, ECMWF, 15/10/07 : NTHRDS increased from 32 to 64
!             Sami Saarinen, ECMWF, 16/10/07 : The default sorting method can be overriden via export EC_SORTING_METHOD=[<number>|<string>]
!             Sami Saarinen, ECMWF, 30/10/07 : Support for CountingSort added as part of QuickSort speedup
!             Sami Saarinen, ECMWF, 31/10/07 : CALL SORTING_METHOD() now prints the prevailing method from EC_SORTING_METHOD
!             Sami Saarinen, ECMWF, 01/11/07 : CountingSort implemented as independent method
!             Sami Saarinen, ECMWF, 06/11/07 : index_adj added as an optional argument to init_index ; new routine adjust_index()
!             Sami Saarinen, ECMWF, 07/11/07 : threshold length "nomp" (see below). Override with EC_SORTING_NOMP.
!             Sami Saarinen, ECMWF, 12/11/07 : Gnome-sort -- the easiest sort on Earth (and very slow for large arrays)
!             Sami Saarinen, ECMWF, 15/11/07 : OpenMP-sorting still under development. Do NOT override the EC_SORTING_NOMP yet.
!             Sami Saarinen, ECMWF, 05/12/07 : When export EC_SORTING_INFO=0, then no info messages are printed from CALL sorting_method()
!             Sami Saarinen, ECMWF, 20/12/07 : export EC_SORTING_INFO=0, rather than =1 is now the default --> less hassling output


IMPLICIT NONE
SAVE
PRIVATE

INTEGER(KIND=JPIM), PARAMETER :: nthrds = 64 ! ***Note: A hardcoded max number of threads !!!

INTEGER(KIND=JPIM), PARAMETER :: sizeof_int4  = 4
INTEGER(KIND=JPIM), PARAMETER :: sizeof_int8  = 8
INTEGER(KIND=JPIM), PARAMETER :: sizeof_real4 = 4
INTEGER(KIND=JPIM), PARAMETER :: sizeof_real8 = 8

INTEGER(KIND=JPIM), PARAMETER :: min_method = 1
INTEGER(KIND=JPIM), PARAMETER :: max_method = 5

INTEGER(KIND=JPIM), PARAMETER :: radixsort_method = 1
INTEGER(KIND=JPIM), PARAMETER :: heapsort_method  = 2
INTEGER(KIND=JPIM), PARAMETER :: quicksort_method = 3
INTEGER(KIND=JPIM), PARAMETER :: countingsort_method = 4
INTEGER(KIND=JPIM), PARAMETER :: gnomesort_method = 5

CHARACTER(LEN=12), PARAMETER :: method_name(min_method:max_method) = &
     & (/&
     &   'RADIXSORT   ' &
     &  ,'HEAPSORT    ' &
     &  ,'QUICKSORT   ' &
     &  ,'COUNTINGSORT' &
     &  ,'GNOMESORT   ' &
     &  /)

!-- Select such method for default_method, which also works for multikey sorts
!   Vector machines should choose radixsort_method, others quicksort_method (oh, sorry, countingsort_method!)
!
!   Note: Occasionally radixsort_method may be faster on non-vector machines, too
#if defined(VPP) || defined(NECSX)
INTEGER(KIND=JPIM) :: default_method = radixsort_method
INTEGER(KIND=JPIM) :: current_method(nthrds) = radixsort_method
#else
INTEGER(KIND=JPIM) :: default_method = countingsort_method
INTEGER(KIND=JPIM) :: current_method(nthrds) = countingsort_method
#endif

!-- A threshold length after which OpenMP in sorting, merging, copying may kick in.
!   Override with EC_SORTING_NOMP. Detected while initializing/calling SORTING_METHOD
!   Non-positive values (<= 0) indicate that OpenMP will NOT be attempted at all
INTEGER(KIND=JPIM) :: nomp = -1

!-- EC_SORTING_INFO = MPL-proc id [1..$NPES] to print the info message when CALL sorting_method()
!                   : 0 (no print; the default)
!                   : 1 (print on MPL-task one)
!                   : > 1 and <= $NPES (some other MPL-task than 1 prints)
INTEGER(KIND=JPIM) :: nsinfo = 0

INTERFACE keysort
MODULE PROCEDURE &
     &int4_keysort_1d, int4_keysort_2d, &
     &int8_keysort_1d, int8_keysort_2d, &
     &real8_keysort_1d, real8_keysort_2d, &
     &real4_keysort_1d, real4_keysort_2d
END INTERFACE

INTERFACE sorting_method
MODULE PROCEDURE int_sorting_method, str_sorting_method
END INTERFACE

PUBLIC :: keysort
PUBLIC :: init_index, get_rank, adjust_index
PUBLIC :: sorting_method

CONTAINS

!----------------------------
!--   Public subroutines   --
!----------------------------

SUBROUTINE int_sorting_method(INEW, IOLD)
INTEGER(KIND=JPIM), OPTIONAL, INTENT(IN)  :: INEW
INTEGER(KIND=JPIM), OPTIONAL, INTENT(OUT) :: IOLD
INTEGER(KIND=JPIM) :: ITMP, IMYPROC, INPES
INTEGER(KIND=JPIM) :: ITID
CHARACTER(LEN=20) CLENV
LOGICAL, SAVE :: LLFIRST = .true.
LOGICAL LLOMP, LLHOOK_OK
REAL(KIND=JPRB) :: ZHOOK_HANDLE
!-- This maybe called from the very first call of DR_HOOK_UTIL ...
llhook_ok = lhook .AND. (PRESENT(inew) .OR. PRESENT(iold))
IF (llhook_ok) CALL dr_hook('ECSORT_MIX:INT_SORTING_METHOD',0,zhook_handle)
itid = oml_my_thread()
llomp = oml_in_parallel()
IF (PRESENT(iold)) iold = current_method(itid)
itmp = -1
IF (PRESENT(inew)) THEN
  itmp = inew
ELSE IF (.NOT.llomp) THEN ! Override the default method (only if outside the OpenMP)
  itmp = -1 ! no change
  IF (llfirst) THEN ! Do once per execution only
#ifdef SFX_MPI
    inpes = mpl_nproc()
#else
    inpes = 1
#endif
    CALL get_environment_variable('EC_SORTING_INFO',clenv) ! ../support/env.c
    IF (clenv /= ' ') THEN
      itmp = nsinfo
      READ(clenv,'(i20)',end=89,err=89) itmp
      GOTO 88
89    CONTINUE
      itmp = nsinfo ! no change
88    CONTINUE
      IF (itmp <= 0) THEN
        nsinfo = 0
      ELSE IF (itmp >= 1 .AND. itmp <= inpes) THEN
        nsinfo = itmp
      ENDIF
    ENDIF
#ifdef SFX_MPI
    imyproc = mpl_myrank()
#else
    imyproc = 0
#endif
    CALL get_environment_variable('EC_SORTING_METHOD',clenv) ! ../support/env.c
    IF (clenv /= ' ') THEN
       IF (imyproc == nsinfo) WRITE(0,'(a)')'<EC_SORTING_METHOD='//trim(clenv)
       CALL toupper(clenv)
       SELECT CASE (clenv)
       CASE ('RADIX', 'RADIXSORT')
          itmp = radixsort_method
       CASE ('HEAP', 'HEAPSORT')
          itmp = heapsort_method
       CASE ('QUICK', 'QUICKSORT', 'QSORT')
          itmp = quicksort_method
       CASE ('COUNT', 'COUNTINGSORT', 'COUNTSORT')
          itmp = countingsort_method
       CASE ('GNOME', 'GNOMESORT')
          itmp = gnomesort_method
       CASE ('DEFAULT', 'DEF')
          itmp = -1 ! no change
       CASE DEFAULT
          READ(clenv,'(i20)',end=99,err=99) itmp
          GOTO 98
99        CONTINUE
          itmp = -1 ! no change
98        CONTINUE
       END SELECT
    ENDIF

    IF (itmp < min_method .OR. itmp > max_method) itmp = default_method
    default_method = itmp
    current_method(:) = itmp
    IF (imyproc == nsinfo) &
         & WRITE(0,'(a,i1,a)')'>EC_SORTING_METHOD=',default_method,&
         & ' # '//method_name(default_method)

    CALL get_environment_variable('EC_SORTING_NOMP',clenv)
    IF (clenv /= ' ') THEN
       IF (imyproc == nsinfo) WRITE(0,'(a)')'<EC_SORTING_NOMP='//trim(clenv)
       READ(clenv,'(i20)',end=199,err=199) itmp
       GOTO 198
199    CONTINUE
       itmp = nomp ! no change
198    CONTINUE
       nomp = itmp
    ENDIF

    IF (imyproc == nsinfo) THEN
      WRITE(clenv,'(i20)') nomp
      WRITE(0,'(a)')'>EC_SORTING_NOMP='//trim(adjustl(clenv))
    ENDIF

    itmp = default_method
    llfirst = .false.
  ENDIF
ENDIF

IF (itmp < min_method .OR. itmp > max_method) itmp = default_method
IF (llomp) THEN ! Only this thread sees the change
  current_method(itid) = itmp
ELSE ! All threads see the change
  current_method(:) = itmp
ENDIF
IF (llhook_ok) CALL dr_hook('ECSORT_MIX:INT_SORTING_METHOD',1,zhook_handle)
END SUBROUTINE int_sorting_method


SUBROUTINE str_sorting_method(CDNEW, IOLD)
CHARACTER(LEN=*), INTENT(IN) :: CDNEW
INTEGER(KIND=JPIM), OPTIONAL, INTENT(OUT) :: IOLD
CHARACTER(LEN=LEN(CDNEW)) CLNEW
REAL(KIND=JPRB) :: ZHOOK_HANDLE
IF (lhook) CALL dr_hook('ECSORT_MIX:STR_SORTING_METHOD',0,zhook_handle)
clnew = cdnew
CALL toupper(clnew)
SELECT CASE (clnew)
CASE ('RADIX', 'RADIXSORT')
  CALL sorting_method(radixsort_method, iold)
CASE ('HEAP', 'HEAPSORT')
  CALL sorting_method(heapsort_method, iold)
CASE ('QUICK', 'QUICKSORT', 'QSORT')
  CALL sorting_method(quicksort_method, iold)
CASE ('COUNT', 'COUNTINGSORT', 'COUNTSORT')
  CALL sorting_method(countingsort_method, iold)
CASE ('GNOME', 'GNOMESORT')
  CALL sorting_method(gnomesort_method, iold)
CASE ('DEFAULT', 'DEF')
  CALL sorting_method(default_method, iold)
CASE DEFAULT
  CALL sorting_method(default_method, iold)
END SELECT
IF (lhook) CALL dr_hook('ECSORT_MIX:STR_SORTING_METHOD',1,zhook_handle)
END SUBROUTINE str_sorting_method


SUBROUTINE init_index(INDEX, INDEX_ADJ)
INTEGER(KIND=JPIM), INTENT(OUT):: INDEX(:)
INTEGER(KIND=JPIM), INTENT(IN), OPTIONAL :: INDEX_ADJ
INTEGER(KIND=JPIM) :: I, N
REAL(KIND=JPRB) :: ZHOOK_HANDLE
IF (lhook) CALL dr_hook('ECSORT_MIX:INIT_INDEX',0,zhook_handle)
n = SIZE(index)
IF (PRESENT(index_adj)) THEN
  DO i=1,n
    index(i) = i + index_adj
  ENDDO
ELSE
  DO i=1,n
    index(i) = i
  ENDDO
ENDIF
IF (lhook) CALL dr_hook('ECSORT_MIX:INIT_INDEX',1,zhook_handle)
END SUBROUTINE init_index


SUBROUTINE adjust_index(INDEX, INDEX_ADJ)
INTEGER(KIND=JPIM), INTENT(INOUT):: INDEX(:)
INTEGER(KIND=JPIM), INTENT(IN) :: INDEX_ADJ
INTEGER(KIND=JPIM) :: I, N
REAL(KIND=JPRB) :: ZHOOK_HANDLE
IF (lhook) CALL dr_hook('ECSORT_MIX:ADJUST_INDEX',0,zhook_handle)
IF (index_adj /= 0) THEN
  n = SIZE(index)
  DO i=1,n
    index(i) = index(i) + index_adj
  ENDDO
ENDIF
IF (lhook) CALL dr_hook('ECSORT_MIX:ADJUST_INDEX',1,zhook_handle)
END SUBROUTINE adjust_index


SUBROUTINE get_rank(INDEX, RANK, INDEX_ADJ)
INTEGER(KIND=JPIM), INTENT(IN) :: INDEX(:)
INTEGER(KIND=JPIM), INTENT(OUT):: RANK(:)
INTEGER(KIND=JPIM), INTENT(IN), OPTIONAL :: INDEX_ADJ
INTEGER(KIND=JPIM) :: I, N
REAL(KIND=JPRB) :: ZHOOK_HANDLE
IF (lhook) CALL dr_hook('ECSORT_MIX:GET_RANK',0,zhook_handle)
n = min(SIZE(index),SIZE(rank))
IF (PRESENT(index_adj)) THEN
  DO i=1,n
    rank(index(i)+index_adj) = i
  ENDDO
ELSE
  DO i=1,n
    rank(index(i)) = i
  ENDDO
ENDIF
IF (lhook) CALL dr_hook('ECSORT_MIX:GET_RANK',1,zhook_handle)
END SUBROUTINE get_rank


#undef INT_VERSION
#undef REAL_VERSION

!-- Create version for INTEGER(KIND=JPIM)
#define INT_VERSION   4
#define REAL_VERSION  0
#include "ecsort_shared.h"
#undef INT_VERSION
#undef REAL_VERSION

!-- Create version for INTEGER(KIND=JPIB)
#define INT_VERSION   8
#define REAL_VERSION  0
#include "ecsort_shared.h"
#undef INT_VERSION
#undef REAL_VERSION

!-- Create version for REAL(KIND=JPRM)
#define INT_VERSION   0
#define REAL_VERSION  4
#include "ecsort_shared.h"
#undef INT_VERSION
#undef REAL_VERSION

!-- Create version for REAL(KIND=JPRB)
#define INT_VERSION   0
#define REAL_VERSION  8
#include "ecsort_shared.h"
#undef INT_VERSION
#undef REAL_VERSION

END MODULE ecsort_mix