tem_mergesort_module.f90 Source File

! Copyright (c) 2012 Metin Cakircali <m.cakircali@grs-sim.de>
! Copyright (c) 2013, 2018 Harald Klimach <harald.klimach@uni-siegen.de>
! Copyright (c) 2013 Simon Zimny <s.zimny@grs-sim.de>
! Copyright (c) 2016 Tobias Schneider <tobias1.schneider@student.uni-siegen.de>
!
! Redistribution and use in source and binary forms, with or without
! modification, are permitted provided that the following conditions are met:
!
! 1. Redistributions of source code must retain the above copyright notice, this
! list of conditions and the following disclaimer.
!
! 2. Redistributions in binary form must reproduce the above copyright notice,
! this list of conditions and the following disclaimer in the documentation
! and/or other materials provided with the distribution.
!
! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
! AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
! IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
! DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
! FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
! DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
! SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
! CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
! OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
! OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
! ****************************************************************************** !
!> description missing
!!
Module tem_mergesort_module

  ! include treelm modules
  use env_module, only : double_k, single_k, int_k, long_k

  ! Module taken from OrderPack 2.0
  ! available at http://www.fortran-2000.com/rank/index.html
  ! from Michel Olagnon

  public :: mrgrnk
  private :: R_mrgrnk, I_mrgrnk, D_mrgrnk, L_mrgrnk

  interface mrgrnk
    module procedure D_mrgrnk, R_mrgrnk, I_mrgrnk, L_mrgrnk
  end interface mrgrnk

  contains

! ****************************************************************************** !
  !> description missing
  !!
  Subroutine D_mrgrnk (XDONT, IRNGT)
    ! ---------------------------------------------------------------------------
    !   MRGRNK = Merge-sort ranking of an array
    !   For performance reasons, the first 2 passes are taken
    !   out of the standard loop, and use dedicated coding.
    ! ---------------------------------------------------------------------------
    Real (kind=double_k), Dimension (:), Intent (In) :: XDONT
    Integer, Dimension (:), Intent (Out) :: IRNGT
    ! ---------------------------------------------------------------------------
    Real (kind=double_k) :: XVALA, XVALB
    Integer, Dimension (SIZE(IRNGT)) :: JWRKT
    Integer :: LMTNA, LMTNC, IRNG1, IRNG2
    Integer :: NVAL, IIND, IWRKD, IWRK, IWRKF, JINDA, IINDA, IINDB
    ! ---------------------------------------------------------------------------

    NVAL = Min (SIZE(XDONT), SIZE(IRNGT))
    Select Case (NVAL)
      Case (:0)
        Return
      Case (1)
        IRNGT (1) = 1
        Return
      Case Default
        Continue
    End Select
    !
    !  Fill-in the index array, creating ordered couples
    !
    Do IIND = 2, NVAL, 2
      If (XDONT(IIND-1) <= XDONT(IIND)) Then
        IRNGT (IIND-1) = IIND - 1
        IRNGT (IIND) = IIND
      Else
        IRNGT (IIND-1) = IIND
        IRNGT (IIND) = IIND - 1
      End If
    End Do
    If (Modulo(NVAL, 2) /= 0) Then
      IRNGT (NVAL) = NVAL
    End If
    !
    !  We will now have ordered subsets A - B - A - B - ...
    !  and merge A and B couples into     C   -   C   - ...
    !
    LMTNA = 2
    LMTNC = 4
    !
    !  First iteration. The length of the ordered subsets goes from 2 to 4
    !
    Do
      If (NVAL <= 2) Exit
      !
      !   Loop on merges of A and B into C
      !
      Do IWRKD = 0, NVAL - 1, 4
        If ((IWRKD+4) > NVAL) Then
          If ((IWRKD+2) >= NVAL) Exit
          !
          !   1 2 3
          !
          If (XDONT(IRNGT(IWRKD+2)) <= XDONT(IRNGT(IWRKD+3))) Exit
          !
          !   1 3 2
          !
          If (XDONT(IRNGT(IWRKD+1)) <= XDONT(IRNGT(IWRKD+3))) Then
            IRNG2 = IRNGT (IWRKD+2)
            IRNGT (IWRKD+2) = IRNGT (IWRKD+3)
            IRNGT (IWRKD+3) = IRNG2
          !
          !   3 1 2
          !
          Else
            IRNG1 = IRNGT (IWRKD+1)
            IRNGT (IWRKD+1) = IRNGT (IWRKD+3)
            IRNGT (IWRKD+3) = IRNGT (IWRKD+2)
            IRNGT (IWRKD+2) = IRNG1
          End If
          Exit
        End If
        !
        !   1 2 3 4
        !
        If (XDONT(IRNGT(IWRKD+2)) <= XDONT(IRNGT(IWRKD+3))) Cycle
        !
        !   1 3 x x
        !
        If (XDONT(IRNGT(IWRKD+1)) <= XDONT(IRNGT(IWRKD+3))) Then
          IRNG2 = IRNGT (IWRKD+2)
          IRNGT (IWRKD+2) = IRNGT (IWRKD+3)
          If (XDONT(IRNG2) <= XDONT(IRNGT(IWRKD+4))) Then
            !   1 3 2 4
            IRNGT (IWRKD+3) = IRNG2
          Else
            !   1 3 4 2
            IRNGT (IWRKD+3) = IRNGT (IWRKD+4)
            IRNGT (IWRKD+4) = IRNG2
          End If
        !
        !   3 x x x
        !
        Else
          IRNG1 = IRNGT (IWRKD+1)
          IRNG2 = IRNGT (IWRKD+2)
          IRNGT (IWRKD+1) = IRNGT (IWRKD+3)
          If (XDONT(IRNG1) <= XDONT(IRNGT(IWRKD+4))) Then
            IRNGT (IWRKD+2) = IRNG1
            If (XDONT(IRNG2) <= XDONT(IRNGT(IWRKD+4))) Then
              !   3 1 2 4
              IRNGT (IWRKD+3) = IRNG2
            Else
              !   3 1 4 2
              IRNGT (IWRKD+3) = IRNGT (IWRKD+4)
              IRNGT (IWRKD+4) = IRNG2
            End If
          Else
            !   3 4 1 2
            IRNGT (IWRKD+2) = IRNGT (IWRKD+4)
            IRNGT (IWRKD+3) = IRNG1
            IRNGT (IWRKD+4) = IRNG2
          End If
        End If
      End Do
      !
      !  The Cs become As and Bs
      !
      LMTNA = 4
      Exit
    End Do
    !
    !  Iteration loop. Each time, the length of the ordered subsets
    !  is doubled.
    !
    Do
      If (LMTNA >= NVAL) Exit
      IWRKF = 0
      LMTNC = 2 * LMTNC
      !
      !   Loop on merges of A and B into C
      !
      Do
        IWRK = IWRKF
        IWRKD = IWRKF + 1
        JINDA = IWRKF + LMTNA
        IWRKF = IWRKF + LMTNC
        If (IWRKF >= NVAL) Then
          If (JINDA >= NVAL) Exit
          IWRKF = NVAL
        End If
        IINDA = 1
        IINDB = JINDA + 1
        !
        !   Shortcut for the case when the max of A is smaller
        !   than the min of B. This line may be activated when the
        !   initial set is already close to sorted.
        !
        !          IF (XDONT(IRNGT(JINDA)) <= XDONT(IRNGT(IINDB))) CYCLE
        !
        !  One steps in the C subset, that we build in the final rank array
        !
        !  Make a copy of the rank array for the merge iteration
        !
        JWRKT (1:LMTNA) = IRNGT (IWRKD:JINDA)
        !
        XVALA = XDONT (JWRKT(IINDA))
        XVALB = XDONT (IRNGT(IINDB))
        !
        Do
          IWRK = IWRK + 1
          !
          !  We still have unprocessed values in both A and B
          !
          If (XVALA > XVALB) Then
            IRNGT (IWRK) = IRNGT (IINDB)
            IINDB = IINDB + 1
            If (IINDB > IWRKF) Then
              !  Only A still with unprocessed values
              IRNGT (IWRK+1:IWRKF) = JWRKT (IINDA:LMTNA)
              Exit
            End If
            XVALB = XDONT (IRNGT(IINDB))
          Else
            IRNGT (IWRK) = JWRKT (IINDA)
            IINDA = IINDA + 1
            If (IINDA > LMTNA) Exit! Only B still with unprocessed values
            XVALA = XDONT (JWRKT(IINDA))
          End If
        !
        End Do
      End Do
      !
      !  The Cs become As and Bs
      !
      LMTNA = 2 * LMTNA
    End Do
    !
    Return
  !
  End Subroutine D_mrgrnk
! ****************************************************************************** !


! ****************************************************************************** !
  !>
  !!
  Subroutine R_mrgrnk (XDONT, IRNGT)
    ! ---------------------------------------------------------------------------
    !   MRGRNK = Merge-sort ranking of an array
    !   For performance reasons, the first 2 passes are taken
    !   out of the standard loop, and use dedicated coding.
    ! ---------------------------------------------------------------------------
    Real(kind=single_k), Dimension (:), Intent (In) :: XDONT
    Integer, Dimension (:), Intent (Out) :: IRNGT
    Real(kind=single_k) :: XVALA, XVALB
    Integer, Dimension (SIZE(IRNGT)) :: JWRKT
    Integer :: LMTNA, LMTNC, IRNG1, IRNG2
    Integer :: NVAL, IIND, IWRKD, IWRK, IWRKF, JINDA, IINDA, IINDB
    ! ---------------------------------------------------------------------------

    NVAL = Min (SIZE(XDONT), SIZE(IRNGT))
    Select Case (NVAL)
      Case (:0)
        Return
      Case (1)
        IRNGT (1) = 1
        Return
      Case Default
        Continue
    End Select
    !
    !  Fill-in the index array, creating ordered couples
    !
    Do IIND = 2, NVAL, 2
      If (XDONT(IIND-1) <= XDONT(IIND)) Then
        IRNGT (IIND-1) = IIND - 1
        IRNGT (IIND) = IIND
      Else
        IRNGT (IIND-1) = IIND
        IRNGT (IIND) = IIND - 1
      End If
    End Do
    If (Modulo(NVAL, 2) /= 0) Then
      IRNGT (NVAL) = NVAL
    End If
    !
    !  We will now have ordered subsets A - B - A - B - ...
    !  and merge A and B couples into     C   -   C   - ...
    !
    LMTNA = 2
    LMTNC = 4
    !
    !  First iteration. The length of the ordered subsets goes from 2 to 4
    !
    Do
      If (NVAL <= 2) Exit
      !
      !   Loop on merges of A and B into C
      !
      Do IWRKD = 0, NVAL - 1, 4
        If ((IWRKD+4) > NVAL) Then
          If ((IWRKD+2) >= NVAL) Exit
          !
          !   1 2 3
          !
          If (XDONT(IRNGT(IWRKD+2)) <= XDONT(IRNGT(IWRKD+3))) Exit
          !
          !   1 3 2
          !
          If (XDONT(IRNGT(IWRKD+1)) <= XDONT(IRNGT(IWRKD+3))) Then
            IRNG2 = IRNGT (IWRKD+2)
            IRNGT (IWRKD+2) = IRNGT (IWRKD+3)
            IRNGT (IWRKD+3) = IRNG2
          !
          !   3 1 2
          !
          Else
            IRNG1 = IRNGT (IWRKD+1)
            IRNGT (IWRKD+1) = IRNGT (IWRKD+3)
            IRNGT (IWRKD+3) = IRNGT (IWRKD+2)
            IRNGT (IWRKD+2) = IRNG1
          End If
          Exit
        End If
        !
        !   1 2 3 4
        !
        If (XDONT(IRNGT(IWRKD+2)) <= XDONT(IRNGT(IWRKD+3))) Cycle
        !
        !   1 3 x x
        !
        If (XDONT(IRNGT(IWRKD+1)) <= XDONT(IRNGT(IWRKD+3))) Then
          IRNG2 = IRNGT (IWRKD+2)
          IRNGT (IWRKD+2) = IRNGT (IWRKD+3)
          If (XDONT(IRNG2) <= XDONT(IRNGT(IWRKD+4))) Then
            !   1 3 2 4
            IRNGT (IWRKD+3) = IRNG2
          Else
            !   1 3 4 2
            IRNGT (IWRKD+3) = IRNGT (IWRKD+4)
            IRNGT (IWRKD+4) = IRNG2
          End If
        !
        !   3 x x x
        !
        Else
          IRNG1 = IRNGT (IWRKD+1)
          IRNG2 = IRNGT (IWRKD+2)
          IRNGT (IWRKD+1) = IRNGT (IWRKD+3)
          If (XDONT(IRNG1) <= XDONT(IRNGT(IWRKD+4))) Then
            IRNGT (IWRKD+2) = IRNG1
            If (XDONT(IRNG2) <= XDONT(IRNGT(IWRKD+4))) Then
              !   3 1 2 4
              IRNGT (IWRKD+3) = IRNG2
            Else
              !   3 1 4 2
              IRNGT (IWRKD+3) = IRNGT (IWRKD+4)
              IRNGT (IWRKD+4) = IRNG2
            End If
          Else
            !   3 4 1 2
            IRNGT (IWRKD+2) = IRNGT (IWRKD+4)
            IRNGT (IWRKD+3) = IRNG1
            IRNGT (IWRKD+4) = IRNG2
          End If
        End If
      End Do
      !
      !  The Cs become As and Bs
      !
      LMTNA = 4
      Exit
    End Do
    !
    !  Iteration loop. Each time, the length of the ordered subsets
    !  is doubled.
    !
    Do
      If (LMTNA >= NVAL) Exit
      IWRKF = 0
      LMTNC = 2 * LMTNC
      !
      !   Loop on merges of A and B into C
      !
      Do
        IWRK = IWRKF
        IWRKD = IWRKF + 1
        JINDA = IWRKF + LMTNA
        IWRKF = IWRKF + LMTNC
        If (IWRKF >= NVAL) Then
          If (JINDA >= NVAL) Exit
          IWRKF = NVAL
        End If
        IINDA = 1
        IINDB = JINDA + 1
        !
        !   Shortcut for the case when the max of A is smaller
        !   than the min of B. This line may be activated when the
        !   initial set is already close to sorted.
        !
        !          IF (XDONT(IRNGT(JINDA)) <= XDONT(IRNGT(IINDB))) CYCLE
        !
        !  One steps in the C subset, that we build in the final rank array
        !
        !  Make a copy of the rank array for the merge iteration
        !
        JWRKT (1:LMTNA) = IRNGT (IWRKD:JINDA)
        !
        XVALA = XDONT (JWRKT(IINDA))
        XVALB = XDONT (IRNGT(IINDB))
        !
        Do
          IWRK = IWRK + 1
          !
          !  We still have unprocessed values in both A and B
          !
          If (XVALA > XVALB) Then
            IRNGT (IWRK) = IRNGT (IINDB)
            IINDB = IINDB + 1
            If (IINDB > IWRKF) Then
              !  Only A still with unprocessed values
              IRNGT (IWRK+1:IWRKF) = JWRKT (IINDA:LMTNA)
              Exit
            End If
            XVALB = XDONT (IRNGT(IINDB))
          Else
            IRNGT (IWRK) = JWRKT (IINDA)
            IINDA = IINDA + 1
            If (IINDA > LMTNA) Exit! Only B still with unprocessed values
            XVALA = XDONT (JWRKT(IINDA))
          End If
        !
        End Do
      End Do
      !
      !  The Cs become As and Bs
      !
      LMTNA = 2 * LMTNA
    End Do
    !
    Return
  !
  End Subroutine R_mrgrnk
! ****************************************************************************** !


! ****************************************************************************** !
  !>
  !!
  Subroutine L_mrgrnk (XDONT, IRNGT)
    ! ---------------------------------------------------------------------------
    !   MRGRNK = Merge-sort ranking of an array
    !   For performance reasons, the first 2 passes are taken
    !   out of the standard loop, and use dedicated coding.
    ! ---------------------------------------------------------------------------
    Integer(kind=long_k), Dimension (:), Intent (In) :: XDONT
    Integer, Dimension (:), Intent (Out) :: IRNGT
    ! ---------------------------------------------------------------------------
    Integer(kind=long_k) :: XVALA, XVALB
    Integer, Dimension (SIZE(IRNGT)) :: JWRKT
    Integer :: LMTNA, LMTNC, IRNG1, IRNG2
    Integer :: NVAL, IIND, IWRKD, IWRK, IWRKF, JINDA, IINDA, IINDB
    ! ---------------------------------------------------------------------------

    NVAL = Min (SIZE(XDONT), SIZE(IRNGT))
    Select Case (NVAL)
      Case (:0)
        Return
      Case (1)
        IRNGT (1) = 1
        Return
      Case Default
        Continue
    End Select
    !
    !  Fill-in the index array, creating ordered couples
    !
    Do IIND = 2, NVAL, 2
      If (XDONT(IIND-1) <= XDONT(IIND)) Then
        IRNGT (IIND-1) = IIND - 1
        IRNGT (IIND) = IIND
      Else
        IRNGT (IIND-1) = IIND
        IRNGT (IIND) = IIND - 1
      End If
    End Do
    If (Modulo(NVAL, 2) /= 0) Then
      IRNGT (NVAL) = NVAL
    End If
    !
    !  We will now have ordered subsets A - B - A - B - ...
    !  and merge A and B couples into     C   -   C   - ...
    !
    LMTNA = 2
    LMTNC = 4
    !
    !  First iteration. The length of the ordered subsets goes from 2 to 4
    !
    Do
      If (NVAL <= 2) Exit
      !
      !   Loop on merges of A and B into C
      !
      Do IWRKD = 0, NVAL - 1, 4
        If ((IWRKD+4) > NVAL) Then
          If ((IWRKD+2) >= NVAL) Exit
          !
          !   1 2 3
          !
          If (XDONT(IRNGT(IWRKD+2)) <= XDONT(IRNGT(IWRKD+3))) Exit
          !
          !   1 3 2
          !
          If (XDONT(IRNGT(IWRKD+1)) <= XDONT(IRNGT(IWRKD+3))) Then
            IRNG2 = IRNGT (IWRKD+2)
            IRNGT (IWRKD+2) = IRNGT (IWRKD+3)
            IRNGT (IWRKD+3) = IRNG2
          !
          !   3 1 2
          !
          Else
            IRNG1 = IRNGT (IWRKD+1)
            IRNGT (IWRKD+1) = IRNGT (IWRKD+3)
            IRNGT (IWRKD+3) = IRNGT (IWRKD+2)
            IRNGT (IWRKD+2) = IRNG1
          End If
          Exit
        End If
        !
        !   1 2 3 4
        !
        If (XDONT(IRNGT(IWRKD+2)) <= XDONT(IRNGT(IWRKD+3))) Cycle
        !
        !   1 3 x x
        !
        If (XDONT(IRNGT(IWRKD+1)) <= XDONT(IRNGT(IWRKD+3))) Then
          IRNG2 = IRNGT (IWRKD+2)
          IRNGT (IWRKD+2) = IRNGT (IWRKD+3)
          If (XDONT(IRNG2) <= XDONT(IRNGT(IWRKD+4))) Then
            !   1 3 2 4
            IRNGT (IWRKD+3) = IRNG2
          Else
            !   1 3 4 2
            IRNGT (IWRKD+3) = IRNGT (IWRKD+4)
            IRNGT (IWRKD+4) = IRNG2
          End If
        !
        !   3 x x x
        !
        Else
          IRNG1 = IRNGT (IWRKD+1)
          IRNG2 = IRNGT (IWRKD+2)
          IRNGT (IWRKD+1) = IRNGT (IWRKD+3)
          If (XDONT(IRNG1) <= XDONT(IRNGT(IWRKD+4))) Then
            IRNGT (IWRKD+2) = IRNG1
            If (XDONT(IRNG2) <= XDONT(IRNGT(IWRKD+4))) Then
              !   3 1 2 4
              IRNGT (IWRKD+3) = IRNG2
            Else
              !   3 1 4 2
              IRNGT (IWRKD+3) = IRNGT (IWRKD+4)
              IRNGT (IWRKD+4) = IRNG2
            End If
          Else
            !   3 4 1 2
            IRNGT (IWRKD+2) = IRNGT (IWRKD+4)
            IRNGT (IWRKD+3) = IRNG1
            IRNGT (IWRKD+4) = IRNG2
          End If
        End If
      End Do
      !
      !  The Cs become As and Bs
      !
      LMTNA = 4
      Exit
    End Do
    !
    !  Iteration loop. Each time, the length of the ordered subsets
    !  is doubled.
    !
    Do
      If (LMTNA >= NVAL) Exit
      IWRKF = 0
      LMTNC = 2 * LMTNC
      !
      !   Loop on merges of A and B into C
      !
      Do
        IWRK = IWRKF
        IWRKD = IWRKF + 1
        JINDA = IWRKF + LMTNA
        IWRKF = IWRKF + LMTNC
        If (IWRKF >= NVAL) Then
          If (JINDA >= NVAL) Exit
          IWRKF = NVAL
        End If
        IINDA = 1
        IINDB = JINDA + 1
        !
        !   Shortcut for the case when the max of A is smaller
        !   than the min of B. This line may be activated when the
        !   initial set is already close to sorted.
        !
        !          IF (XDONT(IRNGT(JINDA)) <= XDONT(IRNGT(IINDB))) CYCLE
        !
        !  One steps in the C subset, that we build in the final rank array
        !
        !  Make a copy of the rank array for the merge iteration
        !
        JWRKT (1:LMTNA) = IRNGT (IWRKD:JINDA)
        !
        XVALA = XDONT (JWRKT(IINDA))
        XVALB = XDONT (IRNGT(IINDB))
        !
        Do
          IWRK = IWRK + 1
          !
          !  We still have unprocessed values in both A and B
          !
          If (XVALA > XVALB) Then
            IRNGT (IWRK) = IRNGT (IINDB)
            IINDB = IINDB + 1
            If (IINDB > IWRKF) Then
              !  Only A still with unprocessed values
              IRNGT (IWRK+1:IWRKF) = JWRKT (IINDA:LMTNA)
              Exit
            End If
            XVALB = XDONT (IRNGT(IINDB))
          Else
            IRNGT (IWRK) = JWRKT (IINDA)
            IINDA = IINDA + 1
            If (IINDA > LMTNA) Exit! Only B still with unprocessed values
            XVALA = XDONT (JWRKT(IINDA))
          End If
        !
        End Do
      End Do
      !
      !  The Cs become As and Bs
      !
      LMTNA = 2 * LMTNA
    End Do
    !
    Return
  !
  End Subroutine L_mrgrnk
! ****************************************************************************** !


! ****************************************************************************** !
  Subroutine I_mrgrnk (XDONT, IRNGT)
    ! ---------------------------------------------------------------------------
    !   MRGRNK = Merge-sort ranking of an array
    !   For performance reasons, the first 2 passes are taken
    !   out of the standard loop, and use dedicated coding.
    ! ---------------------------------------------------------------------------
    Integer(kind=int_k), Dimension (:), Intent (In)  :: XDONT
    Integer, Dimension (:), Intent (Out) :: IRNGT
    Integer(kind=int_k) :: XVALA, XVALB
    Integer, Dimension (SIZE(IRNGT)) :: JWRKT
    Integer :: LMTNA, LMTNC, IRNG1, IRNG2
    Integer :: NVAL, IIND, IWRKD, IWRK, IWRKF, JINDA, IINDA, IINDB
    ! ---------------------------------------------------------------------------

    NVAL = Min (SIZE(XDONT), SIZE(IRNGT))
    Select Case (NVAL)
      Case (:0)
        Return
      Case (1)
        IRNGT (1) = 1
        Return
      Case Default
        Continue
    End Select
    !
    !  Fill-in the index array, creating ordered couples
    !
    Do IIND = 2, NVAL, 2
      If (XDONT(IIND-1) <= XDONT(IIND)) Then
        IRNGT (IIND-1) = IIND - 1
        IRNGT (IIND) = IIND
      Else
        IRNGT (IIND-1) = IIND
        IRNGT (IIND) = IIND - 1
      End If
    End Do
    If (Modulo(NVAL, 2) /= 0) Then
      IRNGT (NVAL) = NVAL
    End If
    !
    !  We will now have ordered subsets A - B - A - B - ...
    !  and merge A and B couples into     C   -   C   - ...
    !
    LMTNA = 2
    LMTNC = 4
    !
    !  First iteration. The length of the ordered subsets goes from 2 to 4
    !
    Do
      If (NVAL <= 2) Exit
      !
      !   Loop on merges of A and B into C
      !
      Do IWRKD = 0, NVAL - 1, 4
        If ((IWRKD+4) > NVAL) Then
          If ((IWRKD+2) >= NVAL) Exit
          !
          !   1 2 3
          !
          If (XDONT(IRNGT(IWRKD+2)) <= XDONT(IRNGT(IWRKD+3))) Exit
          !
          !   1 3 2
          !
          If (XDONT(IRNGT(IWRKD+1)) <= XDONT(IRNGT(IWRKD+3))) Then
            IRNG2 = IRNGT (IWRKD+2)
            IRNGT (IWRKD+2) = IRNGT (IWRKD+3)
            IRNGT (IWRKD+3) = IRNG2
          !
          !   3 1 2
          !
          Else
            IRNG1 = IRNGT (IWRKD+1)
            IRNGT (IWRKD+1) = IRNGT (IWRKD+3)
            IRNGT (IWRKD+3) = IRNGT (IWRKD+2)
            IRNGT (IWRKD+2) = IRNG1
          End If
          Exit
        End If
        !
        !   1 2 3 4
        !
        If (XDONT(IRNGT(IWRKD+2)) <= XDONT(IRNGT(IWRKD+3))) Cycle
        !
        !   1 3 x x
        !
        If (XDONT(IRNGT(IWRKD+1)) <= XDONT(IRNGT(IWRKD+3))) Then
          IRNG2 = IRNGT (IWRKD+2)
          IRNGT (IWRKD+2) = IRNGT (IWRKD+3)
          If (XDONT(IRNG2) <= XDONT(IRNGT(IWRKD+4))) Then
            !   1 3 2 4
            IRNGT (IWRKD+3) = IRNG2
          Else
            !   1 3 4 2
            IRNGT (IWRKD+3) = IRNGT (IWRKD+4)
            IRNGT (IWRKD+4) = IRNG2
          End If
        !
        !   3 x x x
        !
        Else
          IRNG1 = IRNGT (IWRKD+1)
          IRNG2 = IRNGT (IWRKD+2)
          IRNGT (IWRKD+1) = IRNGT (IWRKD+3)
          If (XDONT(IRNG1) <= XDONT(IRNGT(IWRKD+4))) Then
            IRNGT (IWRKD+2) = IRNG1
            If (XDONT(IRNG2) <= XDONT(IRNGT(IWRKD+4))) Then
              !   3 1 2 4
              IRNGT (IWRKD+3) = IRNG2
            Else
              !   3 1 4 2
              IRNGT (IWRKD+3) = IRNGT (IWRKD+4)
              IRNGT (IWRKD+4) = IRNG2
            End If
          Else
            !   3 4 1 2
            IRNGT (IWRKD+2) = IRNGT (IWRKD+4)
            IRNGT (IWRKD+3) = IRNG1
            IRNGT (IWRKD+4) = IRNG2
          End If
        End If
      End Do
      !
      !  The Cs become As and Bs
      !
      LMTNA = 4
      Exit
    End Do
    !
    !  Iteration loop. Each time, the length of the ordered subsets
    !  is doubled.
    !
    Do
      If (LMTNA >= NVAL) Exit
      IWRKF = 0
      LMTNC = 2 * LMTNC
      !
      !   Loop on merges of A and B into C
      !
      Do
        IWRK = IWRKF
        IWRKD = IWRKF + 1
        JINDA = IWRKF + LMTNA
        IWRKF = IWRKF + LMTNC
        If (IWRKF >= NVAL) Then
          If (JINDA >= NVAL) Exit
          IWRKF = NVAL
        End If
        IINDA = 1
        IINDB = JINDA + 1
        !
        !   Shortcut for the case when the max of A is smaller
        !   than the min of B. This line may be activated when the
        !   initial set is already close to sorted.
        !
        !          IF (XDONT(IRNGT(JINDA)) <= XDONT(IRNGT(IINDB))) CYCLE
        !
        !  One steps in the C subset, that we build in the final rank array
        !
        !  Make a copy of the rank array for the merge iteration
        !
        JWRKT (1:LMTNA) = IRNGT (IWRKD:JINDA)
        !
        XVALA = XDONT (JWRKT(IINDA))
        XVALB = XDONT (IRNGT(IINDB))
        !
        Do
          IWRK = IWRK + 1
          !
          !  We still have unprocessed values in both A and B
          !
          If (XVALA > XVALB) Then
            IRNGT (IWRK) = IRNGT (IINDB)
            IINDB = IINDB + 1
            If (IINDB > IWRKF) Then
              !  Only A still with unprocessed values
              IRNGT (IWRK+1:IWRKF) = JWRKT (IINDA:LMTNA)
              Exit
            End If
            XVALB = XDONT (IRNGT(IINDB))
          Else
            IRNGT (IWRK) = JWRKT (IINDA)
            IINDA = IINDA + 1
            If (IINDA > LMTNA) Exit! Only B still with unprocessed values
            XVALA = XDONT (JWRKT(IINDA))
          End If
        !
        End Do
      End Do
      !
      !  The Cs become As and Bs
      !
      LMTNA = 2 * LMTNA
    End Do
    !
    Return
  !
  End Subroutine I_mrgrnk
! ****************************************************************************** !


end module tem_mergesort_module
! ****************************************************************************** !