!  Copyright (C) 2015-2019 Richard Weed.
!  All rights reserved.
  
!  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, in whole or in part, must retain the  
!  above copyright notice, this list of conditions and the following 
!  disclaimer.
  
!  2. Redistributions in binary form, in whole or in part, 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.
  
!  3. The names of the contributors may not be used to endorse or promote from 
!  products derived from this software without specific prior written 
!  permission.

!  4. Redistributions of this software, in whole or in part, in any form, 
!  must be freely available and licensed under this original License. The 
!  U.S. Government may add additional restrictions to their modified and 
!  redistributed software as required by Law. However, these restrictions 
!  do not apply to the original software distribution.
 
!  5. Redistribution of this source code, including any modifications, may 
!  not be intentionally obfuscated.
  
!  6. Other code may make use of this software, in whole or in part, without 
!  restriction, provided that it does not apply any restriction to this 
!  software other than outlined above.
  
!  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 HOLDERS AND
!  CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, 
!  EXEMPLARARY 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.


! Original Licence:

! This code is patterned after code taken from Hansons and Hopkins, 
! "Numerical Computation with Modern Fortran", SIAM, 2013. and the
! associated Web site, http://www.siam.org/books/ot134. SIAM's copyright
! license information, and restrictions are given as follows:

! Copyright (C) 2013 by the Society of Industrial and Applied Mathematics

! SIAM grants a royalty-free license to copy and distribute the sofware code
! posted on the book's supplemental Web page, provided the source is
! acknowledged:

! No warranties, expressed or implied are made by the publisher, authors, and
! their employers, that the programs contained in this volume are free of 
! error. They should not be relied on as the sole basis to solve a problem 
! whose incorrect solution could result in injury to person or property.
! if the programs are employed in such a manner, it is at the user's own risk
! and the publisher, authors, and their employers disclaim all liability for
! such misuse. 

Module quickSort

! Implements type explicit quicksort routines for real, integer, and character
! arrays along with a user defined type derived from the User_t base class.
! REAL32, REAL64, INT32, INT64, Character(LEN=1) arrays, and strings are 
! supported along with the user class.

! Arrays default to sorting into ascending order. This can be
! changed by setting the optional argument order to order=DESCENDING
! for any of the overloaded quicksort interfaces. Two forms of the
! quicksort routines are provided, one sorts the arrays in place
! and changes the values in the array, the second doesn't change
! the values but returns a permutation (index) array that can
! be used to reorder the input array outside of qsort

! Note access to the different versions are through the overloaded
! qsort interfaces

! Examples:

!  Real(REAL32) :: a(50)
!  Integer      :: index(50)

!  Call qsort(a)                  ! sorts in place in ascending order 
!  Call qsort(a, PERMUTE=index)   ! returns permutation index, a is not changed
!  Call qsort(a  ORDER=DESCENDING) ! sorts in decending order

  Use ISO_FORTRAN_ENV, ONLY: REAL32, REAL64, INT8, INT16, INT32, INT64
  Use userClass,       ONLY: User_t

  Implicit NONE

  Private

! Define crossover point for switching from qsort to insertion sort

  Integer, PARAMETER, PRIVATE :: SWITCHSORTS = 10 

! Define some global logical variables to define/select sort order

  Logical, Parameter  :: ASCENDING  = .TRUE. 
  Logical, Parameter  :: DESCENDING = .FALSE. 

  Interface qsort
    Module Procedure qsortR32
    Module Procedure qsortR64
    Module Procedure qsortI8
    Module Procedure qsortI16
    Module Procedure qsortI32
    Module Procedure qsortI64
    Module Procedure qsortCharArray
    Module Procedure qsortCharString
    Module Procedure qsortUser

  End Interface

  Public :: qsort, ASCENDING, DESCENDING

Contains

! Sort 32 bit real array. Wraps sort in place or by permutation routines

  Subroutine qsortR32(a, left, right, permute, order)

    Implicit NONE

    Real(REAL32), Contiguous, Intent(INOUT) :: a(:)
 
    Integer,      Optional,   Intent(IN)    :: left 
    Integer,      Optional,   Intent(IN)    :: right
    Integer,      Optional,   Intent(INOUT) :: permute(:)
    Logical,      Optional,   Intent(IN)    :: order

    Integer :: is, ie
    Logical :: sortorder

    is        = 1
    ie        = SIZE(a)
    sortorder = ASCENDING
    If (PRESENT(left))  is        = left
    If (PRESENT(right)) ie        = right
    If (PRESENT(order)) sortorder = order  

    If (PRESENT(permute)) Then
      Call qsortIndexR32(a, is, ie, permute, sortorder)
    Else
      Call qsortInPlaceR32(a, is, ie, sortorder)
    EndIf

  End Subroutine qsortR32

! Sort 64 bit real array. Wraps sort in place or by permutation routines
  
  Subroutine qsortR64(a, left, right, permute, order)

    Implicit NONE

    Real(REAL64), Contiguous, Intent(INOUT) :: a(:)
 
    Integer,      Optional,   Intent(IN)    :: left 
    Integer,      Optional,   Intent(IN)    :: right
    Integer,      Optional,   Intent(INOUT) :: permute(:)
    Logical,      Optional,   Intent(IN)    :: order

    Integer :: is, ie
    Logical :: sortorder

    is        = 1
    ie        = SIZE(a)
    sortorder = ASCENDING
    If (PRESENT(left))  is        = left
    If (PRESENT(right)) ie        = right
    If (PRESENT(order)) sortorder = order  

    If (PRESENT(permute)) Then
      Call qsortIndexR64(a, is, ie, permute, sortorder)
    Else
      Call qsortInPlaceR64(a, is, ie, sortorder)
    EndIf

  End Subroutine qsortR64

! Sort 8 bit integer array. Wraps sort in place or by permutation routines

  Subroutine qsortI8(a, left, right, permute, order)

    Implicit NONE

    Integer(INT8), Contiguous, Intent(INOUT) :: a(:)
    Integer,       Optional,   Intent(IN)    :: left 
    Integer,       Optional,   Intent(IN)    :: right
    Integer,       Optional,   Intent(INOUT) :: permute(:)
    Logical,       Optional,   Intent(IN)    :: order

    Integer :: is, ie
    Logical :: sortorder

    is        = 1
    ie        = SIZE(a)
    sortorder = ASCENDING
    If (PRESENT(left))  is        = left
    If (PRESENT(right)) ie        = right
    If (PRESENT(order)) sortorder = order  

    If (PRESENT(permute)) Then
      Call qsortIndexI8(a, is, ie, permute, sortorder)
    Else
      Call qsortInPlaceI8(a, is, ie, sortorder)
    EndIf

  End Subroutine qsortI8

! Sort 16 bit integer array. Wraps sort in place or by permutation routines

  Subroutine qsortI16(a, left, right, permute, order)

    Implicit NONE

    Integer(INT16), Contiguous, Intent(INOUT) :: a(:)
    Integer,        Optional,   Intent(IN)    :: left 
    Integer,        Optional,   Intent(IN)    :: right
    Integer,        Optional,   Intent(INOUT) :: permute(:)
    Logical,        Optional,   Intent(IN)    :: order

    Integer :: is, ie
    Logical :: sortorder

    is        = 1
    ie        = SIZE(a)
    sortorder = ASCENDING
    If (PRESENT(left))  is        = left
    If (PRESENT(right)) ie        = right
    If (PRESENT(order)) sortorder = order  

    If (PRESENT(permute)) Then
      Call qsortIndexI16(a, is, ie, permute, sortorder)
    Else
      Call qsortInPlaceI16(a, is, ie, sortorder)
    EndIf

  End Subroutine qsortI16

! Sort 32 bit integer array. Wraps sort in place or by permutation routines

  Subroutine qsortI32(a, left, right, permute, order)

    Implicit NONE

    Integer(INT32), Contiguous, Intent(INOUT) :: a(:)
    Integer,        Optional,   Intent(IN)    :: left 
    Integer,        Optional,   Intent(IN)    :: right
    Integer,        Optional,   Intent(INOUT) :: permute(:)
    Logical,        Optional,   Intent(IN)    :: order

    Integer :: is, ie
    Logical :: sortorder

    is        = 1
    ie        = SIZE(a)
    sortorder = ASCENDING
    If (PRESENT(left))  is        = left
    If (PRESENT(right)) ie        = right
    If (PRESENT(order)) sortorder = order  

    If (PRESENT(permute)) Then
      Call qsortIndexI32(a, is, ie, permute, sortorder)
    Else
      Call qsortInPlaceI32(a, is, ie, sortorder)
    EndIf

  End Subroutine qsortI32

! Sort 64 bit integer array. Wraps sort in place or by permutation routines

  Subroutine qsortI64(a, left, right, permute, order)

    Implicit NONE

    Integer(INT64), Contiguous, Intent(INOUT) :: a(:)
    Integer,        Optional,   Intent(IN)    :: left 
    Integer,        Optional,   Intent(IN)    :: right
    Integer,        Optional,   Intent(INOUT) :: permute(:)
    Logical,        Optional,   Intent(IN)    :: order

    Integer :: is, ie
    Logical :: sortorder

    is        = 1
    ie        = SIZE(a)
    sortorder = ASCENDING
    If (PRESENT(left))  is        = left
    If (PRESENT(right)) ie        = right
    If (PRESENT(order)) sortorder = order  

    If (PRESENT(permute)) Then
      Call qsortIndexI64(a, is, ie, permute, sortorder)
    Else
      Call qsortInPlaceI64(a, is, ie, sortorder)
    EndIf

  End Subroutine qsortI64

! Sort Character strings or array. Wraps sort in place or by permutation routines

  Subroutine qsortCharArray(a, left, right, permute, order)

    Implicit NONE

    Character(LEN=1), Contiguous, Intent(INOUT) :: a(:)
    Integer,          Optional,   Intent(IN)    :: left 
    Integer,          Optional,   Intent(IN)    :: right
    Integer,          Optional,   Intent(INOUT) :: permute(:)
    Logical,          Optional,   Intent(IN)    :: order

    Integer :: is, ie
    Logical :: sortorder

    is        = 1
    ie        = SIZE(a)
    sortorder = ASCENDING
    If (PRESENT(left))  is        = left
    If (PRESENT(right)) ie        = right
    If (PRESENT(order)) sortorder = order  

    If (PRESENT(permute)) Then
      Call qsortIndexChar(a, is, ie, permute, sortorder)
    Else
      Call qsortInPlaceChar(a, is, ie, sortorder)
    EndIf

  End Subroutine qsortCharArray

  Subroutine qsortCharString(astring, left, right, permute, order)

    Implicit NONE

    Character(LEN=*),           Intent(INOUT) :: astring
    Integer,          Optional, Intent(IN)    :: left 
    Integer,          Optional, Intent(IN)    :: right
    Integer,          Optional, Intent(INOUT) :: permute(:)
    Logical,          Optional, Intent(IN)    :: order

    Character(LEN=1), ALLOCATABLE :: a(:)

    Integer :: is, ie, i
    Logical :: sortorder

    ALLOCATE(a(LEN(astring)))
    Do i=1,LEN(astring)
      a(i) = astring(i:i)
    EndDo

    is        = 1
    ie        = SIZE(a)
    sortorder = ASCENDING
    If (PRESENT(left))  is        = left
    If (PRESENT(right)) ie        = right
    If (PRESENT(order)) sortorder = order  

! convert string to array


    If (PRESENT(permute)) Then
      Call qsortIndexChar(a, is, ie, permute, sortorder)
    Else
      Call qsortInPlaceChar(a, is, ie, sortorder)
      astring = REPEAT(" ", LEN(astring))
      Do i=1,LEN(astring)
        astring(i:i) = a(i)
      EndDo
    EndIf

    If (ALLOCATED(a)) DEALLOCATE(a)

  End Subroutine qsortCharString

! Sort User_t array. Wraps sort in place or by permutation routines

  Subroutine qsortUser(a, left, right, permute, order)

    Implicit NONE

    Class(User_t), Contiguous, Target, Intent(INOUT) :: a(:)
    Integer,       Optional,           Intent(IN)    :: left 
    Integer,       Optional,           Intent(IN)    :: right
    Integer,       Optional,           Intent(INOUT) :: permute(:)
    Logical,       Optional,           Intent(IN)    :: order

    Integer :: is, ie
    Logical :: sortorder

    is        = 1
    ie        = SIZE(a)
    sortorder = ASCENDING
    If (PRESENT(left))  is        = left
    If (PRESENT(right)) ie        = right
    If (PRESENT(order)) sortorder = order  

    If (PRESENT(permute)) Then
      Call qsortIndexUser(a, is, ie, permute, sortorder)
    Else
      Call qsortInPlaceUser(a, is, ie, sortorder)
    EndIf

!    If (ALLOCATED(saveUser)) Deallocate(saveUser)

  End Subroutine qsortUser

! Sort 32 bit real arrays

  Subroutine qsortInPlaceR32(a, left, right, order) 

    Implicit NONE

    Integer,      Intent(IN)             :: left, right
    Real(REAL32), Intent(INOUT)          :: a(left:right)
    Logical,      Intent(IN),   OPTIONAL :: order

    Logical      :: sortorder
    Integer      :: savedIndex
    Real(REAL32) :: savedR32
    
    savedIndex = 0 
    savedR32   = 0.0_REAL32 
    sortorder = ASCENDING
    If (PRESENT(order)) sortorder=order
    
    Call quicksortIP(left, right)

    Call insertionSort

  Contains
      
    Subroutine exchange(i, j)

      Implicit NONE

      Integer, Intent(IN) :: i, j

      Real(REAL32)   :: tr32

      tr32 = a(i)
      a(i) = a(j)
      a(j) = tr32

    End Subroutine exchange

    Logical Function compare(i, j)

      Implicit NONE

      Integer, Intent(IN) :: i, j

      If (sortorder .EQV. ASCENDING) Then
        compare = a(i) < a(j)
      Else
        compare = a(i) > a(j)
      EndIf

    End Function compare

    Subroutine compex(i,j)

      Implicit NONE

      Integer, Intent(IN) :: i, j

      If (compare(j,i)) Call exchange(i,j)

    End Subroutine compex

    Subroutine moveValue(i, j)

      Implicit NONE

      Integer, Intent(IN) :: i, j

      a(i) = a(j)

    End Subroutine moveValue
    Logical Function compareValue(j)

      Implicit NONE

      Integer, Intent(IN) :: j

      If (sortorder .EQV. ASCENDING) Then
        compareValue = savedR32 < a(j)
      Else
        compareValue = savedR32 > a(j)
      EndIf

    End Function compareValue

    Subroutine saveValue(i)

      Implicit NONE

      Integer, Intent(IN) :: i

      savedR32 = a(i)

    End Subroutine saveValue

    Subroutine restoreValue(i)

      Implicit NONE

      Integer, Intent(IN) :: i

      a(i) = savedR32

    End Subroutine restoreValue

    Subroutine insertionSort()

      Implicit NONE

      Integer :: i, j

      Do i=left+1, right
        Call compex(left,i)
      EndDo

      Do i=left+2, right

        j=i
        Call saveValue(i)

        Do While(compareValue(j-1))
          Call moveValue(j,j-1)
          j=j-1
        EndDo
     
        Call restoreValue(j)

      EndDo

    End Subroutine insertionSort

    Function partition(left, right) RESULT(i)

      Implicit NONE

      Integer, Intent(IN) :: left, right
      Integer :: i, j

      i = left - 1
      j = right

      Call saveValue(right)
   
      Loop1 : Do

        Loop2 : Do

          i = i+1
      
          If (i>right) EXIT Loop2
          If (compareValue(i)) EXIT Loop2
        
        EndDo Loop2

        Loop3: Do

          j = j-1
          If (.NOT. compareValue(j) .OR. j==left) EXIT LOOP3

        EndDo Loop3

        If (i>=j) EXIT Loop1
        Call exchange(i, j)

      EndDo Loop1

      Call exchange(i, right)

    End Function partition

    Recursive Subroutine quickSortIP(left, right)

      Implicit NONE

      Integer, Intent(IN) :: left, right

      Integer :: i

      If ((right-left) > SWITCHSORTS) Then

        Call exchange((right+left)/2, (right-1))
        Call compex(left, right-1) 
        Call compex(right, left) 
        Call compex(right-1, right)
        i    = partition(left+1, right-1)
        Call quicksortIP(left, i-1)
        Call quicksortIP(i+1, right)
      EndIf

    End Subroutine quicksortIP

  End Subroutine qsortInPlaceR32

  Subroutine qsortIndexR32(a, left, right, permute, order)

    Implicit NONE

    Integer,      Intent(IN)           :: left, right
    Real(REAL32), Intent(INOUT)        :: a(left:right)
    Integer,      Intent(INOUT)        :: permute(left:right)
    Logical,      Intent(IN), OPTIONAL :: order

    Integer      :: j
    Integer      :: savedIndex
    Logical      :: sortorder
    Real(REAL32) :: savedR32
   
    savedIndex = 0 
    savedR32   = 0.0_REAL32 

    sortorder = ASCENDING
    If (PRESENT(order)) sortorder = order
  
    Do j=left,right

      permute(j) = j

    EndDo 

    Call quicksortIndex(left, right)

    Call insertionSort

  Contains
 
    Subroutine exchange(i, j)

      Implicit NONE

      Integer, Intent(IN) :: i, j

      Integer :: tndx 
      
      tndx = permute(i)
      permute(i) = permute(j)
      permute(j) = tndx

    End Subroutine exchange

    Logical Function compare(i, j)

      Implicit NONE

      Integer, Intent(IN) :: i, j

      Integer :: ii, ij

      ii = permute(i)
      ij = permute(j)

      If (sortorder .EQV. ASCENDING) Then
        compare = a(ii) < a(ij)
      Else
        compare = a(ii) > a(ij)
      EndIf

    End Function compare

    Subroutine compex(i,j)

      Implicit NONE

      Integer, Intent(IN) :: i, j

      If (compare(j,i)) Call exchange(i,j)

    End Subroutine compex

    Subroutine moveValue(i, j)

      Implicit NONE

      Integer, Intent(IN) :: i, j

      permute(i) = permute(j)

    End Subroutine moveValue

    Logical Function compareValue(j)

      Implicit NONE

      Integer, Intent(IN) :: j

      Integer :: ij

      ij = permute(j)
      If (sortorder .EQV. ASCENDING) Then
        compareValue = savedR32 < a(ij)
      Else
        compareValue = savedR32 > a(ij)
      EndIf

    End Function compareValue

    Subroutine saveValue(i)

      Implicit NONE

      Integer, Intent(IN) :: i
      
      savedIndex = permute(i)
      savedR32   = a(savedIndex)
      
    End Subroutine saveValue
    Subroutine restoreValue(i)

      Implicit NONE

      Integer, Intent(IN) :: i

      permute(i) = savedIndex

    End Subroutine restoreValue

    Subroutine insertionSort()

      Implicit NONE

      Integer :: i, j

      Do i=left+1, right
        Call compex(left,i)
      EndDo

      Do i=left+2, right

        j=i
        Call saveValue(i)

        Do While(compareValue(j-1))
          Call moveValue(j,j-1)
          j=j-1
        EndDo
     
        Call restoreValue(j)

      EndDo

    End Subroutine insertionSort

    Function partition(left, right) RESULT(i)

      Implicit NONE

      Integer, Intent(IN) :: left, right
      Integer :: i, j

      i = left - 1
      j = right

      Call saveValue(right)
   
      Loop1 : Do

        Loop2 : Do

          i = i+1
      
          If (i>right) EXIT Loop2
          If (compareValue(i)) EXIT Loop2
        
        EndDo Loop2

        Loop3: Do

          j = j-1
          If (.NOT. compareValue(j) .OR. j==left) EXIT LOOP3

        EndDo Loop3

        If (i>=j) EXIT Loop1
        Call exchange(i, j)

      EndDo Loop1

      Call exchange(i, right)

    End Function partition

    Recursive Subroutine quickSortIndex(left, right)

      Implicit NONE

      Integer, Intent(IN) :: left, right

      Integer :: i

      If ((right-left) > SWITCHSORTS) Then

        Call exchange((right+left)/2, (right-1))
        Call compex(left, right-1) 
        Call compex(right, left) 
        Call compex(right-1, right)
        i    = partition(left+1, right-1)
        Call quicksortIndex(left, i-1)
        Call quicksortIndex(i+1, right)
   
      EndIf

    End Subroutine quickSortIndex

  End Subroutine qsortIndexR32

! Sort 64 bit real arrays

  Subroutine qsortInPlaceR64(a, left, right, order)

    Implicit NONE

    Integer,      Intent(IN)             :: left, right
    Real(REAL64), Intent(INOUT)          :: a(left:right)
    Logical,      Intent(IN),   OPTIONAL :: order

    Logical      :: sortorder
    Integer      :: savedIndex
    Real(REAL64) :: savedR64
   
    savedIndex = 0 
    savedR64   = 0.0_REAL32 

    sortorder = ASCENDING
    If (PRESENT(order)) sortorder = order

    Call quicksortIP(left, right)

    Call insertionSort

  Contains
      
    Subroutine exchange(i, j)

      Implicit NONE

      Integer, Intent(IN) :: i, j

      Real(REAL64)   :: tr64

      tr64 = a(i)
      a(i) = a(j)
      a(j) = tr64

    End Subroutine exchange

    Logical Function compare(i, j)

      Implicit NONE

      Integer, Intent(IN) :: i, j

      If (sortorder .EQV. ASCENDING) Then
        compare = a(i) < a(j)
      Else
        compare = a(i) > a(j)
      EndIf

    End Function compare

    Subroutine compex(i,j)

      Implicit NONE

      Integer, Intent(IN) :: i, j

      If (compare(j,i)) Call exchange(i,j)

    End Subroutine compex

    Subroutine moveValue(i, j)

      Implicit NONE

      Integer, Intent(IN) :: i, j

      a(i) = a(j)

    End Subroutine moveValue
    Logical Function compareValue(j)

      Implicit NONE

      Integer, Intent(IN) :: j

      If (sortorder .EQV. ASCENDING) Then
        compareValue = savedR64 < a(j)
      Else
        compareValue = savedR64 > a(j)
      EndIf

    End Function compareValue

    Subroutine saveValue(i)

      Implicit NONE

      Integer, Intent(IN) :: i

      savedR64 = a(i)

    End Subroutine saveValue

    Subroutine restoreValue(i)

      Implicit NONE

      Integer, Intent(IN) :: i

      a(i) = savedR64

    End Subroutine restoreValue

    Subroutine insertionSort()

      Implicit NONE

      Integer :: i, j

      Do i=left+1, right
        Call compex(left,i)
      EndDo

      Do i=left+2, right

        j=i
        Call saveValue(i)

        Do While(compareValue(j-1))
          Call moveValue(j,j-1)
          j=j-1
        EndDo
     
        Call restoreValue(j)

      EndDo

    End Subroutine insertionSort

    Function partition(left, right) RESULT(i)

      Implicit NONE

      Integer, Intent(IN) :: left, right
      Integer :: i, j

      i = left - 1
      j = right

      Call saveValue(right)
   
      Loop1 : Do

        Loop2 : Do

          i = i+1
      
          If (i>right) EXIT Loop2
          If (compareValue(i)) EXIT Loop2
        
        EndDo Loop2

        Loop3: Do

          j = j-1
          If (.NOT. compareValue(j) .OR. j==left) EXIT LOOP3

        EndDo Loop3

        If (i>=j) EXIT Loop1
        Call exchange(i, j)

      EndDo Loop1

      Call exchange(i, right)

    End Function partition

    Recursive Subroutine quickSortIP(left, right)

      Implicit NONE

      Integer, Intent(IN) :: left, right

      Integer :: i

      If ((right-left) > SWITCHSORTS) Then

        Call exchange((right+left)/2, (right-1))
        Call compex(left, right-1) 
        Call compex(right, left) 
        Call compex(right-1, right)
        i    = partition(left+1, right-1)
        Call quicksortIP(left, i-1)
        Call quicksortIP(i+1, right)
      EndIf

    End Subroutine quicksortIP

  End Subroutine qsortInPlaceR64

  Subroutine qsortIndexR64(a, left, right, permute, order)

    Implicit NONE

    Integer,      Intent(IN)            :: left, right
    Real(REAL64), Intent(INOUT)         :: a(left:right)
    Integer,      Intent(INOUT)         :: permute(left:right)
    Logical,      Intent(IN),  OPTIONAL :: order

    Integer      :: savedIndex
    Integer      :: j
    Logical      :: sortorder
    Real(REAL64) :: savedR64
    
    savedIndex = 0
    savedR64   = 0.0_REAL64 

    sortorder = ASCENDING
    If(PRESENT(order)) sortorder = order

    Do j=left,right

      permute(j) = j

    EndDo 

    Call quicksortIndex(left, right)

    Call insertionSort

  Contains
 
    Subroutine exchange(i, j)

      Implicit NONE

      Integer, Intent(IN) :: i, j

      Integer :: tndx 
      
      tndx       = permute(i)
      permute(i) = permute(j)
      permute(j) = tndx

    End Subroutine exchange

    Logical Function compare(i, j)

      Implicit NONE

      Integer, Intent(IN) :: i, j

      Integer :: ii, ij

      ii = permute(i)
      ij = permute(j)

      If (sortorder .EQV. ASCENDING) Then
        compare = a(ii) < a(ij)
      Else
        compare = a(ii) > a(ij)
      EndIf

    End Function compare

    Subroutine compex(i,j)

      Implicit NONE

      Integer, Intent(IN) :: i, j

      If (compare(j,i)) Call exchange(i,j)

    End Subroutine compex

    Subroutine moveValue(i, j)

      Implicit NONE

      Integer, Intent(IN) :: i, j

      permute(i) = permute(j)

    End Subroutine moveValue

    Logical Function compareValue(j)

      Implicit NONE

      Integer, Intent(IN) :: j

      Integer :: ij

      ij           = permute(j)
      If (sortorder .EQV. ASCENDING) Then
        compareValue = savedR64 < a(ij)
      Else
        compareValue = savedR64 > a(ij)
      EndIf

    End Function compareValue

    Subroutine saveValue(i)

      Implicit NONE

      Integer, Intent(IN) :: i
      
      savedIndex = permute(i)
      savedR64 = a(savedIndex)
      
    End Subroutine saveValue

    Subroutine restoreValue(i)

      Implicit NONE

      Integer, Intent(IN) :: i

      permute(i) = savedIndex

    End Subroutine restoreValue

    Subroutine insertionSort()

      Implicit NONE

      Integer :: i, j

      Do i=left+1, right
        Call compex(left,i)
      EndDo

      Do i=left+2, right

        j=i
        Call saveValue(i)

        Do While(compareValue(j-1))
          Call moveValue(j,j-1)
          j=j-1
        EndDo
     
        Call restoreValue(j)

      EndDo

    End Subroutine insertionSort

    Function partition(left, right) RESULT(i)

      Implicit NONE

      Integer, Intent(IN) :: left, right
      Integer :: i, j

      i = left - 1
      j = right

      Call saveValue(right)
   
      Loop1 : Do

        Loop2 : Do

          i = i+1
      
          If (i>right) EXIT Loop2
          If (compareValue(i)) EXIT Loop2
        
        EndDo Loop2

        Loop3: Do

          j = j-1
          If (.NOT. compareValue(j) .OR. j==left) EXIT LOOP3

        EndDo Loop3

        If (i>=j) EXIT Loop1
        Call exchange(i, j)

      EndDo Loop1

      Call exchange(i, right)

    End Function partition

    Recursive Subroutine quickSortIndex(left, right)

      Implicit NONE

      Integer, Intent(IN) :: left, right

      Integer :: i

      If ((right-left) > SWITCHSORTS) Then

        Call exchange((right+left)/2, (right-1))
        Call compex(left, right-1) 
        Call compex(right, left) 
        Call compex(right-1, right)
        i    = partition(left+1, right-1)
        Call quicksortIndex(left, i-1)
        Call quicksortIndex(i+1, right)
   
      EndIf

    End Subroutine quickSortIndex

  End Subroutine qsortIndexR64

! Sort 8 bit integer arrays

  Subroutine qsortInPlaceI8(a, left, right, order)

    Implicit NONE

    Integer,       Intent(IN)             :: left, right
    Integer(INT8), Intent(INOUT)          :: a(left:right)
    Logical,       Intent(IN),   OPTIONAL :: order

    Integer       :: savedIndex
    Logical       :: sortorder
    Integer(INT8) :: savedI8
   
    savedIndex = 0 
    savedI8    = 0_INT8

    sortorder = ASCENDING
    If (PRESENT(order)) sortorder = order

    Call quicksortIP(left, right)

    Call insertionSort

  Contains
      
    Subroutine exchange(i, j)

      Implicit NONE

      Integer, Intent(IN) :: i, j

      Integer(INT8) :: ti8

      ti8 = a(i)
      a(i) = a(j)
      a(j) = ti8

    End Subroutine exchange

    Logical Function compare(i, j)

      Implicit NONE

      Integer, Intent(IN) :: i, j

      If (sortorder .EQV. ASCENDING) Then
        compare = a(i) < a(j)
      Else
        compare = a(i) > a(j)
      EndIf

    End Function compare

    Subroutine compex(i,j)

      Implicit NONE

      Integer, Intent(IN) :: i, j

      If (compare(j,i)) Call exchange(i,j)

    End Subroutine compex

    Subroutine moveValue(i, j)

      Implicit NONE

      Integer, Intent(IN) :: i, j

      a(i) = a(j)

    End Subroutine moveValue
    Logical Function compareValue(j)

      Implicit NONE

      Integer, Intent(IN) :: j

      If (sortorder .EQV. ASCENDING) Then
        compareValue = savedI8 < a(j)
      Else
        compareValue = savedI8 > a(j)
      EndIf

    End Function compareValue

    Subroutine saveValue(i)

      Implicit NONE

      Integer, Intent(IN) :: i

      savedI8 = a(i)

    End Subroutine saveValue

    Subroutine restoreValue(i)

      Implicit NONE

      Integer, Intent(IN) :: i

      a(i) = savedI8

    End Subroutine restoreValue

    Subroutine insertionSort()

      Implicit NONE

      Integer :: i, j

      Do i=left+1, right
        Call compex(left,i)
      EndDo

      Do i=left+2, right

        j=i
        Call saveValue(i)

        Do While(compareValue(j-1))
          Call moveValue(j,j-1)
          j=j-1
        EndDo
     
        Call restoreValue(j)

      EndDo

    End Subroutine insertionSort

    Function partition(left, right) RESULT(i)

      Implicit NONE

      Integer, Intent(IN) :: left, right
      Integer :: i, j

      i = left - 1
      j = right

      Call saveValue(right)
   
      Loop1 : Do

        Loop2 : Do

          i = i+1
      
          If (i>right) EXIT Loop2
          If (compareValue(i)) EXIT Loop2
        
        EndDo Loop2

        Loop3: Do

          j = j-1
          If (.NOT. compareValue(j) .OR. j==left) EXIT LOOP3

        EndDo Loop3

        If (i>=j) EXIT Loop1
        Call exchange(i, j)

      EndDo Loop1

      Call exchange(i, right)

    End Function partition

    Recursive Subroutine quickSortIP(left, right)

      Implicit NONE

      Integer, Intent(IN) :: left, right

      Integer :: i

      If ((right-left) > SWITCHSORTS) Then

        Call exchange((right+left)/2, (right-1))
        Call compex(left, right-1) 
        Call compex(right, left) 
        Call compex(right-1, right)
        i    = partition(left+1, right-1)
        Call quicksortIP(left, i-1)
        Call quicksortIP(i+1, right)
      EndIf

    End Subroutine quicksortIP

  End Subroutine qsortInPlaceI8

  Subroutine qsortIndexI8(a, left, right, permute, order)

    Implicit NONE

    Integer,       Intent(IN)             :: left, right
    Integer(INT8), Intent(INOUT)          :: a(left:right)
    Integer,       Intent(INOUT)          :: permute(left:right)
    Logical,       Intent(IN),   OPTIONAL :: order

    Integer       :: j
    Integer       :: savedIndex
    Logical       :: sortorder
    Integer(INT8) :: savedI8
   
    savedIndex = 0 
    savedI8    = 0_INT8 

    sortorder = ASCENDING
    If (PRESENT(order)) sortorder = order

    Call quicksortIndex(left, right)

    Call insertionSort

  Contains
 
    Subroutine exchange(i, j)

      Implicit NONE

      Integer, Intent(IN) :: i, j

      Integer :: tndx 
      
      tndx = permute(i)
      permute(i) = permute(j)
      permute(j) = tndx

    End Subroutine exchange

    Logical Function compare(i, j)

      Implicit NONE

      Integer, Intent(IN) :: i, j

      Integer :: ii, ij

      ii = permute(i)
      ij = permute(j)

      If (sortorder .EQV. ASCENDING) Then
        compare = a(ii) < a(ij)
      Else
        compare = a(ii) > a(ij)
      EndIf

    End Function compare

    Subroutine compex(i,j)

      Implicit NONE

      Integer, Intent(IN) :: i, j

      If (compare(j,i)) Call exchange(i,j)

    End Subroutine compex

    Subroutine moveValue(i, j)

      Implicit NONE

      Integer, Intent(IN) :: i, j

      permute(i) = permute(j)

    End Subroutine moveValue

    Logical Function compareValue(j)

      Implicit NONE

      Integer, Intent(IN) :: j

      Integer :: ij

      ij = permute(j)
      If (sortorder .EQV. ASCENDING) Then
        compareValue = savedI8 < a(ij)
      Else
        compareValue = savedI8 > a(ij)
      EndIf

    End Function compareValue

    Subroutine saveValue(i)

      Implicit NONE

      Integer, Intent(IN) :: i
      
      savedIndex = permute(i)
      savedI8 = a(savedIndex)
      
    End Subroutine saveValue
    Subroutine restoreValue(i)

      Implicit NONE

      Integer, Intent(IN) :: i

      permute(i) = savedIndex

    End Subroutine restoreValue

    Subroutine insertionSort()

      Implicit NONE

      Integer :: i, j

      Do i=left+1, right
        Call compex(left,i)
      EndDo

      Do i=left+2, right

        j=i
        Call saveValue(i)

        Do While(compareValue(j-1))
          Call moveValue(j,j-1)
          j=j-1
        EndDo
     
        Call restoreValue(j)

      EndDo

    End Subroutine insertionSort

    Function partition(left, right) RESULT(i)

      Implicit NONE

      Integer, Intent(IN) :: left, right
      Integer :: i, j

      i = left - 1
      j = right

      Call saveValue(right)
   
      Loop1 : Do

        Loop2 : Do

          i = i+1
      
          If (i>right) EXIT Loop2
          If (compareValue(i)) EXIT Loop2
        
        EndDo Loop2

        Loop3: Do

          j = j-1
          If (.NOT. compareValue(j) .OR. j==left) EXIT LOOP3

        EndDo Loop3

        If (i>=j) EXIT Loop1
        Call exchange(i, j)

      EndDo Loop1

      Call exchange(i, right)

    End Function partition

    Recursive Subroutine quickSortIndex(left, right)

      Implicit NONE

      Integer, Intent(IN) :: left, right

      Integer :: i

      If ((right-left) > SWITCHSORTS) Then

        Call exchange((right+left)/2, (right-1))
        Call compex(left, right-1) 
        Call compex(right, left) 
        Call compex(right-1, right)
        i    = partition(left+1, right-1)
        Call quicksortIndex(left, i-1)
        Call quicksortIndex(i+1, right)
   
      EndIf

    End Subroutine quickSortIndex

  End Subroutine qsortIndexI8

! 16 bit integer in place or by permutation

  Subroutine qsortInPlaceI16(a, left, right, order)

    Implicit NONE

    Integer,        Intent(IN)             :: left, right
    Integer(INT16), Intent(INOUT)          :: a(left:right)
    Logical,        Intent(IN),   OPTIONAL :: order

    Integer        :: savedIndex
    Logical        :: sortorder
    Integer(INT16) :: savedI16
   
    savedIndex = 0 
    savedI16   = 0_INT16

    sortorder = ASCENDING
    If (PRESENT(order)) sortorder = order

    Call quicksortIP(left, right)

    Call insertionSort

  Contains
      
    Subroutine exchange(i, j)

      Implicit NONE

      Integer, Intent(IN) :: i, j

      Integer(INT16) :: ti16

      ti16 = a(i)
      a(i) = a(j)
      a(j) = ti16

    End Subroutine exchange

    Logical Function compare(i, j)

      Implicit NONE

      Integer, Intent(IN) :: i, j

      If (sortorder .EQV. ASCENDING) Then
        compare = a(i) < a(j)
      Else
        compare = a(i) > a(j)
      EndIf

    End Function compare

    Subroutine compex(i,j)

      Implicit NONE

      Integer, Intent(IN) :: i, j

      If (compare(j,i)) Call exchange(i,j)

    End Subroutine compex

    Subroutine moveValue(i, j)

      Implicit NONE

      Integer, Intent(IN) :: i, j

      a(i) = a(j)

    End Subroutine moveValue

    Logical Function compareValue(j)

      Implicit NONE

      Integer, Intent(IN) :: j

      If (sortorder .EQV. ASCENDING) Then
        compareValue = savedI16 < a(j)
      Else
        compareValue = savedI16 > a(j)
      EndIf

    End Function compareValue

    Subroutine saveValue(i)

      Implicit NONE

      Integer, Intent(IN) :: i

      savedI16 = a(i)

    End Subroutine saveValue

    Subroutine restoreValue(i)

      Implicit NONE

      Integer, Intent(IN) :: i

      a(i) = savedI16

    End Subroutine restoreValue

    Subroutine insertionSort()

      Implicit NONE

      Integer :: i, j

      Do i=left+1, right
        Call compex(left,i)
      EndDo

      Do i=left+2, right

        j=i
        Call saveValue(i)

        Do While(compareValue(j-1))
          Call moveValue(j,j-1)
          j=j-1
        EndDo
     
        Call restoreValue(j)

      EndDo

    End Subroutine insertionSort

    Function partition(left, right) RESULT(i)

      Implicit NONE

      Integer, Intent(IN) :: left, right
      Integer :: i, j

      i = left - 1
      j = right

      Call saveValue(right)
   
      Loop1 : Do

        Loop2 : Do

          i = i+1
      
          If (i>right) EXIT Loop2
          If (compareValue(i)) EXIT Loop2
        
        EndDo Loop2

        Loop3: Do

          j = j-1
          If (.NOT. compareValue(j) .OR. j==left) EXIT LOOP3

        EndDo Loop3

        If (i>=j) EXIT Loop1
        Call exchange(i, j)

      EndDo Loop1

      Call exchange(i, right)

    End Function partition

    Recursive Subroutine quickSortIP(left, right)

      Implicit NONE

      Integer, Intent(IN) :: left, right

      Integer :: i

      If ((right-left) > SWITCHSORTS) Then

        Call exchange((right+left)/2, (right-1))
        Call compex(left, right-1) 
        Call compex(right, left) 
        Call compex(right-1, right)
        i    = partition(left+1, right-1)
        Call quicksortIP(left, i-1)
        Call quicksortIP(i+1, right)
      EndIf

    End Subroutine quicksortIP

  End Subroutine qsortInPlaceI16

  Subroutine qsortIndexI16(a, left, right, permute, order)

    Implicit NONE

    Integer,        Intent(IN)             :: left, right
    Integer(INT16), Intent(INOUT)          :: a(left:right)
    Integer,        Intent(INOUT)          :: permute(left:right)
    Logical,        Intent(IN),   OPTIONAL :: order

    Integer        :: j
    Integer        :: savedIndex
    Logical        :: sortorder
    Integer(INT16) :: savedI16
   
    savedIndex = 0 
    savedI16   = 0_INT16 

    sortorder = ASCENDING
    If (PRESENT(order)) sortorder = order

    Call quicksortIndex(left, right)

    Call insertionSort

  Contains
 
    Subroutine exchange(i, j)

      Implicit NONE

      Integer, Intent(IN) :: i, j

      Integer :: tndx 
      
      tndx = permute(i)
      permute(i) = permute(j)
      permute(j) = tndx

    End Subroutine exchange

    Logical Function compare(i, j)

      Implicit NONE

      Integer, Intent(IN) :: i, j

      Integer :: ii, ij

      ii = permute(i)
      ij = permute(j)

      If (sortorder .EQV. ASCENDING) Then
        compare = a(ii) < a(ij)
      Else
        compare = a(ii) > a(ij)
      EndIf

    End Function compare

    Subroutine compex(i,j)

      Implicit NONE

      Integer, Intent(IN) :: i, j

      If (compare(j,i)) Call exchange(i,j)

    End Subroutine compex

    Subroutine moveValue(i, j)

      Implicit NONE

      Integer, Intent(IN) :: i, j

      permute(i) = permute(j)

    End Subroutine moveValue

    Logical Function compareValue(j)

      Implicit NONE

      Integer, Intent(IN) :: j

      Integer :: ij

      ij = permute(j)
      If (sortorder .EQV. ASCENDING) Then
        compareValue = savedI16 < a(ij)
      Else
        compareValue = savedI16 > a(ij)
      EndIf

    End Function compareValue

    Subroutine saveValue(i)

      Implicit NONE

      Integer, Intent(IN) :: i
      
      savedIndex = permute(i)
      savedI16 = a(savedIndex)
      
    End Subroutine saveValue
    Subroutine restoreValue(i)

      Implicit NONE

      Integer, Intent(IN) :: i

      permute(i) = savedIndex

    End Subroutine restoreValue

    Subroutine insertionSort()

      Implicit NONE

      Integer :: i, j

      Do i=left+1, right
        Call compex(left,i)
      EndDo

      Do i=left+2, right

        j=i
        Call saveValue(i)

        Do While(compareValue(j-1))
          Call moveValue(j,j-1)
          j=j-1
        EndDo
     
        Call restoreValue(j)

      EndDo

    End Subroutine insertionSort

    Function partition(left, right) RESULT(i)

      Implicit NONE

      Integer, Intent(IN) :: left, right
      Integer :: i, j

      i = left - 1
      j = right

      Call saveValue(right)
   
      Loop1 : Do

        Loop2 : Do

          i = i+1
      
          If (i>right) EXIT Loop2
          If (compareValue(i)) EXIT Loop2
        
        EndDo Loop2

        Loop3: Do

          j = j-1
          If (.NOT. compareValue(j) .OR. j==left) EXIT LOOP3

        EndDo Loop3

        If (i>=j) EXIT Loop1
        Call exchange(i, j)

      EndDo Loop1

      Call exchange(i, right)

    End Function partition

    Recursive Subroutine quickSortIndex(left, right)

      Implicit NONE

      Integer, Intent(IN) :: left, right

      Integer :: i

      If ((right-left) > SWITCHSORTS) Then

        Call exchange((right+left)/2, (right-1))
        Call compex(left, right-1) 
        Call compex(right, left) 
        Call compex(right-1, right)
        i    = partition(left+1, right-1)
        Call quicksortIndex(left, i-1)
        Call quicksortIndex(i+1, right)
   
      EndIf

    End Subroutine quickSortIndex

  End Subroutine qsortIndexI16

! Sort 32 bit integer arrays

  Subroutine qsortInPlaceI32(a, left, right, order)

    Implicit NONE

    Integer,        Intent(IN)             :: left, right
    Integer(INT32), Intent(INOUT)          :: a(left:right)
    Logical,        Intent(IN),   OPTIONAL :: order

    Integer        :: savedIndex
    Logical        :: sortorder
    Integer(INT32) :: savedI32
   
    savedIndex = 0 
    savedI32   = 0_INT32 

    sortorder = ASCENDING
    If (PRESENT(order)) sortorder = order

    Call quicksortIP(left, right)

    Call insertionSort

  Contains
      
    Subroutine exchange(i, j)

      Implicit NONE

      Integer, Intent(IN) :: i, j

      Integer(INT32) :: ti32

      ti32 = a(i)
      a(i) = a(j)
      a(j) = ti32

    End Subroutine exchange

    Logical Function compare(i, j)

      Implicit NONE

      Integer, Intent(IN) :: i, j

      If (sortorder .EQV. ASCENDING) Then
        compare = a(i) < a(j)
      Else
        compare = a(i) > a(j)
      EndIf

    End Function compare

    Subroutine compex(i,j)

      Implicit NONE

      Integer, Intent(IN) :: i, j

      If (compare(j,i)) Call exchange(i,j)

    End Subroutine compex

    Subroutine moveValue(i, j)

      Implicit NONE

      Integer, Intent(IN) :: i, j

      a(i) = a(j)

    End Subroutine moveValue
    Logical Function compareValue(j)

      Implicit NONE

      Integer, Intent(IN) :: j

      If (sortorder .EQV. ASCENDING) Then
        compareValue = savedI32 < a(j)
      Else
        compareValue = savedI32 > a(j)
      EndIf

    End Function compareValue

    Subroutine saveValue(i)

      Implicit NONE

      Integer, Intent(IN) :: i

      savedI32 = a(i)

    End Subroutine saveValue

    Subroutine restoreValue(i)

      Implicit NONE

      Integer, Intent(IN) :: i

      a(i) = savedI32

    End Subroutine restoreValue

    Subroutine insertionSort()

      Implicit NONE

      Integer :: i, j

      Do i=left+1, right
        Call compex(left,i)
      EndDo

      Do i=left+2, right

        j=i
        Call saveValue(i)

        Do While(compareValue(j-1))
          Call moveValue(j,j-1)
          j=j-1
        EndDo
     
        Call restoreValue(j)

      EndDo

    End Subroutine insertionSort

    Function partition(left, right) RESULT(i)

      Implicit NONE

      Integer, Intent(IN) :: left, right
      Integer :: i, j

      i = left - 1
      j = right

      Call saveValue(right)
   
      Loop1 : Do

        Loop2 : Do

          i = i+1
      
          If (i>right) EXIT Loop2
          If (compareValue(i)) EXIT Loop2
        
        EndDo Loop2

        Loop3: Do

          j = j-1
          If (.NOT. compareValue(j) .OR. j==left) EXIT LOOP3

        EndDo Loop3

        If (i>=j) EXIT Loop1
        Call exchange(i, j)

      EndDo Loop1

      Call exchange(i, right)

    End Function partition

    Recursive Subroutine quickSortIP(left, right)

      Implicit NONE

      Integer, Intent(IN) :: left, right

      Integer :: i

      If ((right-left) > SWITCHSORTS) Then

        Call exchange((right+left)/2, (right-1))
        Call compex(left, right-1) 
        Call compex(right, left) 
        Call compex(right-1, right)
        i    = partition(left+1, right-1)
        Call quicksortIP(left, i-1)
        Call quicksortIP(i+1, right)
      EndIf

    End Subroutine quicksortIP

  End Subroutine qsortInPlaceI32

  Subroutine qsortIndexI32(a, left, right, permute, order)

    Implicit NONE

    Integer,        Intent(IN)             :: left, right
    Integer(INT32), Intent(INOUT)          :: a(left:right)
    Integer,        Intent(INOUT)          :: permute(left:right)
    Logical,        Intent(IN),   OPTIONAL :: order

    Integer        :: j
    Integer        :: savedIndex
    Logical        :: sortorder
    Integer(INT32) :: savedI32
   
    savedIndex = 0 
    savedI32   = 0_INT32 

    sortorder = ASCENDING
    If (PRESENT(order)) sortorder = order

    Call quicksortIndex(left, right)

    Call insertionSort

  Contains
 
    Subroutine exchange(i, j)

      Implicit NONE

      Integer, Intent(IN) :: i, j

      Integer :: tndx 
      
      tndx = permute(i)
      permute(i) = permute(j)
      permute(j) = tndx

    End Subroutine exchange

    Logical Function compare(i, j)

      Implicit NONE

      Integer, Intent(IN) :: i, j

      Integer :: ii, ij

      ii = permute(i)
      ij = permute(j)

      If (sortorder .EQV. ASCENDING) Then
        compare = a(ii) < a(ij)
      Else
        compare = a(ii) > a(ij)
      EndIf

    End Function compare

    Subroutine compex(i,j)

      Implicit NONE

      Integer, Intent(IN) :: i, j

      If (compare(j,i)) Call exchange(i,j)

    End Subroutine compex

    Subroutine moveValue(i, j)

      Implicit NONE

      Integer, Intent(IN) :: i, j

      permute(i) = permute(j)

    End Subroutine moveValue

    Logical Function compareValue(j)

      Implicit NONE

      Integer, Intent(IN) :: j

      Integer :: ij

      ij = permute(j)
      If (sortorder .EQV. ASCENDING) Then
        compareValue = savedI32 < a(ij)
      Else
        compareValue = savedI32 > a(ij)
      EndIf

    End Function compareValue

    Subroutine saveValue(i)

      Implicit NONE

      Integer, Intent(IN) :: i
      
      savedIndex = permute(i)
      savedI32 = a(savedIndex)
      
    End Subroutine saveValue
    Subroutine restoreValue(i)

      Implicit NONE

      Integer, Intent(IN) :: i

      permute(i) = savedIndex

    End Subroutine restoreValue

    Subroutine insertionSort()

      Implicit NONE

      Integer :: i, j

      Do i=left+1, right
        Call compex(left,i)
      EndDo

      Do i=left+2, right

        j=i
        Call saveValue(i)

        Do While(compareValue(j-1))
          Call moveValue(j,j-1)
          j=j-1
        EndDo
     
        Call restoreValue(j)

      EndDo

    End Subroutine insertionSort

    Function partition(left, right) RESULT(i)

      Implicit NONE

      Integer, Intent(IN) :: left, right
      Integer :: i, j

      i = left - 1
      j = right

      Call saveValue(right)
   
      Loop1 : Do

        Loop2 : Do

          i = i+1
      
          If (i>right) EXIT Loop2
          If (compareValue(i)) EXIT Loop2
        
        EndDo Loop2

        Loop3: Do

          j = j-1
          If (.NOT. compareValue(j) .OR. j==left) EXIT LOOP3

        EndDo Loop3

        If (i>=j) EXIT Loop1
        Call exchange(i, j)

      EndDo Loop1

      Call exchange(i, right)

    End Function partition

    Recursive Subroutine quickSortIndex(left, right)

      Implicit NONE

      Integer, Intent(IN) :: left, right

      Integer :: i

      If ((right-left) > SWITCHSORTS) Then

        Call exchange((right+left)/2, (right-1))
        Call compex(left, right-1) 
        Call compex(right, left) 
        Call compex(right-1, right)
        i    = partition(left+1, right-1)
        Call quicksortIndex(left, i-1)
        Call quicksortIndex(i+1, right)
   
      EndIf

    End Subroutine quickSortIndex

  End Subroutine qsortIndexI32

! Sort 64 bit integer arrays

  Subroutine qsortInPlaceI64(a, left, right, order)

    Implicit NONE

    Integer,        Intent(IN)             :: left, right
    Integer(INT64), Intent(INOUT)          :: a(left:right)
    Logical,        Intent(IN),   OPTIONAL :: order

    Logical        :: sortorder
    Integer        :: savedIndex
    Integer(INT64) :: savedI64
   
    savedIndex = 0 
    savedI64    = 0_INT64 

    sortorder = ASCENDING
    If (PRESENT(order)) sortorder = order

    Call quicksortIP(left, right)

    Call insertionSort

  Contains
      
    Subroutine exchange(i, j)

      Implicit NONE

      Integer, Intent(IN) :: i, j

      Integer(INT64) :: ti64

      ti64 = a(i)
      a(i) = a(j)
      a(j) = ti64

    End Subroutine exchange

    Logical Function compare(i, j)

      Implicit NONE

      Integer, Intent(IN) :: i, j

      If (sortorder .EQV. ASCENDING) Then
        compare = a(i) < a(j)
      Else
        compare = a(i) > a(j)
      EndIf

    End Function compare

    Subroutine compex(i,j)

      Implicit NONE

      Integer, Intent(IN) :: i, j

      If (compare(j,i)) Call exchange(i,j)

    End Subroutine compex

    Subroutine moveValue(i, j)

      Implicit NONE

      Integer, Intent(IN) :: i, j

      a(i) = a(j)

    End Subroutine moveValue
    Logical Function compareValue(j)

      Implicit NONE

      Integer, Intent(IN) :: j

      If (sortorder .EQV. ASCENDING) Then
        compareValue = savedI64 < a(j)
      Else
        compareValue = savedI64 > a(j)
      EndIf

    End Function compareValue

    Subroutine saveValue(i)

      Implicit NONE

      Integer, Intent(IN) :: i

      savedI64 = a(i)

    End Subroutine saveValue

    Subroutine restoreValue(i)

      Implicit NONE

      Integer, Intent(IN) :: i

      a(i) = savedI64

    End Subroutine restoreValue

    Subroutine insertionSort()

      Implicit NONE

      Integer :: i, j

      Do i=left+1, right
        Call compex(left,i)
      EndDo

      Do i=left+2, right

        j=i
        Call saveValue(i)

        Do While(compareValue(j-1))
          Call moveValue(j,j-1)
          j=j-1
        EndDo
     
        Call restoreValue(j)

      EndDo

    End Subroutine insertionSort

    Function partition(left, right) RESULT(i)

      Implicit NONE

      Integer, Intent(IN) :: left, right
      Integer :: i, j

      i = left - 1
      j = right

      Call saveValue(right)
   
      Loop1 : Do

        Loop2 : Do

          i = i+1
      
          If (i>right) EXIT Loop2
          If (compareValue(i)) EXIT Loop2
        
        EndDo Loop2

        Loop3: Do

          j = j-1
          If (.NOT. compareValue(j) .OR. j==left) EXIT LOOP3

        EndDo Loop3

        If (i>=j) EXIT Loop1
        Call exchange(i, j)

      EndDo Loop1

      Call exchange(i, right)

    End Function partition

    Recursive Subroutine quickSortIP(left, right)

      Implicit NONE

      Integer, Intent(IN) :: left, right

      Integer :: i

      If ((right-left) > SWITCHSORTS) Then

        Call exchange((right+left)/2, (right-1))
        Call compex(left, right-1) 
        Call compex(right, left) 
        Call compex(right-1, right)
        i    = partition(left+1, right-1)
        Call quicksortIP(left, i-1)
        Call quicksortIP(i+1, right)
      EndIf

    End Subroutine quicksortIP

  End Subroutine qsortInPlaceI64

  Subroutine qsortIndexI64(a, left, right, permute, order)

    Implicit NONE

    Integer,        Intent(IN)             :: left, right
    Integer(INT64), Intent(INOUT)          :: a(left:right)
    Integer,        Intent(INOUT)          :: permute(left:right)
    Logical,        Intent(IN),   OPTIONAL :: order

    Integer        :: j
    Integer        :: savedIndex
    Logical        :: sortorder
    Integer(INT64) :: savedI64
   
    savedIndex = 0 
    savedI64   = 0_INT64 

    sortorder = ASCENDING
    If (PRESENT(order)) sortorder = order

    Do j=left,right

      permute(j) = j

    EndDo 

    Call quicksortIndex(left, right)

    Call insertionSort

  Contains
 
    Subroutine exchange(i, j)

      Implicit NONE

      Integer, Intent(IN) :: i, j

      Integer :: tndx 
      
      tndx = permute(i)
      permute(i) = permute(j)
      permute(j) = tndx

    End Subroutine exchange

    Logical Function compare(i, j)

      Implicit NONE

      Integer, Intent(IN) :: i, j

      Integer :: ii, ij

      ii = permute(i)
      ij = permute(j)

      If (sortorder .EQV. ASCENDING) Then
        compare = a(ii) < a(ij)
      Else
        compare = a(ii) > a(ij)
      EndIf

    End Function compare

    Subroutine compex(i,j)

      Implicit NONE

      Integer, Intent(IN) :: i, j

      If (compare(j,i)) Call exchange(i,j)

    End Subroutine compex

    Subroutine moveValue(i, j)

      Implicit NONE

      Integer, Intent(IN) :: i, j

      permute(i) = permute(j)

    End Subroutine moveValue

    Logical Function compareValue(j)

      Implicit NONE

      Integer, Intent(IN) :: j

      Integer :: ij

      ij = permute(j)
      If (sortorder .EQV. ASCENDING) Then
        compareValue = savedI64 < a(ij)
      Else
        compareValue = savedI64 > a(ij)
      EndIf

    End Function compareValue

    Subroutine saveValue(i)

      Implicit NONE

      Integer, Intent(IN) :: i
      
      savedIndex = permute(i)
      savedI64 = a(savedIndex)
      
    End Subroutine saveValue
    Subroutine restoreValue(i)

      Implicit NONE

      Integer, Intent(IN) :: i

      permute(i) = savedIndex

    End Subroutine restoreValue

    Subroutine insertionSort()

      Implicit NONE

      Integer :: i, j

      Do i=left+1, right
        Call compex(left,i)
      EndDo

      Do i=left+2, right

        j=i
        Call saveValue(i)

        Do While(compareValue(j-1))
          Call moveValue(j,j-1)
          j=j-1
        EndDo
     
        Call restoreValue(j)

      EndDo

    End Subroutine insertionSort

    Function partition(left, right) RESULT(i)

      Implicit NONE

      Integer, Intent(IN) :: left, right
      Integer :: i, j

      i = left - 1
      j = right

      Call saveValue(right)
   
      Loop1 : Do

        Loop2 : Do

          i = i+1
      
          If (i>right) EXIT Loop2
          If (compareValue(i)) EXIT Loop2
        
        EndDo Loop2

        Loop3: Do

          j = j-1
          If (.NOT. compareValue(j) .OR. j==left) EXIT LOOP3

        EndDo Loop3

        If (i>=j) EXIT Loop1
        Call exchange(i, j)

      EndDo Loop1

      Call exchange(i, right)

    End Function partition

    Recursive Subroutine quickSortIndex(left, right)

      Implicit NONE

      Integer, Intent(IN) :: left, right

      Integer :: i

      If ((right-left) > SWITCHSORTS) Then

        Call exchange((right+left)/2, (right-1))
        Call compex(left, right-1) 
        Call compex(right, left) 
        Call compex(right-1, right)
        i    = partition(left+1, right-1)
        Call quicksortIndex(left, i-1)
        Call quicksortIndex(i+1, right)
   
      EndIf

    End Subroutine quickSortIndex

  End Subroutine qsortIndexI64

! Sort character arrays

  Subroutine qsortInPlaceChar(a, left, right, order)

    Implicit NONE

    Integer,          Intent(IN)             :: left, right
    Character(LEN=1), Intent(INOUT)          :: a(left:right)
    Logical,          Intent(IN),   OPTIONAL :: order

    Integer          :: savedIndex
    Logical          :: sortorder
    Character(LEN=1) :: savedChar
   
    savedIndex = 0 
    savedChar  = ' '
    sortorder  = ASCENDING
    If (PRESENT(order)) sortorder = order

    Call quicksortIP(left, right)

    Call insertionSort

  Contains
      
    Subroutine exchange(i, j)

      Implicit NONE

      Integer, Intent(IN) :: i, j

      Character(LEN=1) :: tchar

      tchar = a(i)
      a(i)  = a(j)
      a(j)  = tchar 

    End Subroutine exchange

    Logical Function compare(i, j)

      Implicit NONE

      Integer, Intent(IN) :: i, j

      If (sortorder .EQV. ASCENDING) Then
        compare = LLT(a(i), a(j))
      Else
        compare = LGT(a(i),a(j))
      EndIf

    End Function compare

    Subroutine compex(i,j)

      Implicit NONE

      Integer, Intent(IN) :: i, j

      If (compare(j,i)) Call exchange(i,j)

    End Subroutine compex

    Subroutine moveValue(i, j)

      Implicit NONE

      Integer, Intent(IN) :: i, j

      a(i) = a(j)

    End Subroutine moveValue
    Logical Function compareValue(j)

      Implicit NONE

      Integer, Intent(IN) :: j

      If (sortorder .EQV. ASCENDING) Then
        compareValue = LLT(savedChar, a(j))
      Else
        compareValue = LGT(savedChar, a(j))
      EndIf

    End Function compareValue

    Subroutine saveValue(i)

      Implicit NONE

      Integer, Intent(IN) :: i

      savedChar = a(i)

    End Subroutine saveValue

    Subroutine restoreValue(i)

      Implicit NONE

      Integer, Intent(IN) :: i

      a(i) = savedChar

    End Subroutine restoreValue

    Subroutine insertionSort()

      Implicit NONE

      Integer :: i, j

      Do i=left+1, right
        Call compex(left,i)
      EndDo

      Do i=left+2, right

        j=i
        Call saveValue(i)

        Do While(compareValue(j-1))
          Call moveValue(j,j-1)
          j=j-1
        EndDo
     
        Call restoreValue(j)

      EndDo

    End Subroutine insertionSort

    Function partition(left, right) RESULT(i)

      Implicit NONE

      Integer, Intent(IN) :: left, right
      Integer :: i, j

      i = left - 1
      j = right

      Call saveValue(right)
   
      Loop1 : Do

        Loop2 : Do

          i = i+1
      
          If (i>right) EXIT Loop2
          If (compareValue(i)) EXIT Loop2
        
        EndDo Loop2

        Loop3: Do

          j = j-1
          If (.NOT. compareValue(j) .OR. j==left) EXIT LOOP3

        EndDo Loop3

        If (i>=j) EXIT Loop1
        Call exchange(i, j)

      EndDo Loop1

      Call exchange(i, right)

    End Function partition

    Recursive Subroutine quickSortIP(left, right)

      Implicit NONE

      Integer, Intent(IN) :: left, right

      Integer :: i

      If ((right-left) > SWITCHSORTS) Then

        Call exchange((right+left)/2, (right-1))
        Call compex(left, right-1) 
        Call compex(right, left) 
        Call compex(right-1, right)
        i    = partition(left+1, right-1)
        Call quicksortIP(left, i-1)
        Call quicksortIP(i+1, right)
      EndIf

    End Subroutine quicksortIP

  End Subroutine qsortInPlaceChar

  Subroutine qsortIndexChar(a, left, right, permute, order)

    Implicit NONE

    Integer,          Intent(IN)             :: left, right
    Character(LEN=1), Intent(INOUT)          :: a(left:right)
    Integer,          Intent(INOUT)          :: permute(left:right)
    Logical,          Intent(IN),   OPTIONAL :: order

    Integer          :: j
    Integer          :: savedIndex
    Logical          :: sortorder
    Character(LEN=1) :: savedChar
   
    savedIndex = 0 
    savedChar  = ' '
    sortorder = ASCENDING
    If (PRESENT(order)) sortorder = order

    Do j=left,right

      permute(j) = j

    EndDo 

    Call quicksortIndex(left, right)

    Call insertionSort

  Contains
 
    Subroutine exchange(i, j)

      Implicit NONE

      Integer, Intent(IN) :: i, j

      Integer :: tndx
      
      tndx = permute(i)
      permute(i) = permute(j)
      permute(j) = tndx

    End Subroutine exchange

    Logical Function compare(i, j)

      Implicit NONE

      Integer, Intent(IN) :: i, j

      Integer :: ii, ij

      ii = permute(i)
      ij = permute(j)

      If (sortorder .EQV. ASCENDING) Then
        compare = LLT(a(ii), a(ij))
      Else
        compare = LGT(a(ii), a(ij))
      EndIf

    End Function compare

    Subroutine compex(i,j)

      Implicit NONE

      Integer, Intent(IN) :: i, j

      If (compare(j,i)) Call exchange(i,j)

    End Subroutine compex

    Subroutine moveValue(i, j)

      Implicit NONE

      Integer, Intent(IN) :: i, j

      permute(i) = permute(j)

    End Subroutine moveValue

    Logical Function compareValue(j)

      Implicit NONE

      Integer, Intent(IN) :: j

      Integer :: ij

      ij = permute(j)
      If (sortorder .EQV. ASCENDING) Then
        compareValue = LLT(savedChar, a(ij))
      Else
        compareValue = LGT(savedChar, a(ij))
      EndIf

    End Function compareValue

    Subroutine saveValue(i)

      Implicit NONE

      Integer, Intent(IN) :: i
      
      savedIndex = permute(i)
      savedChar = a(savedIndex)
      
    End Subroutine saveValue
    Subroutine restoreValue(i)

      Implicit NONE

      Integer, Intent(IN) :: i

      permute(i) = savedIndex

    End Subroutine restoreValue

    Subroutine insertionSort()

      Implicit NONE

      Integer :: i, j

      Do i=left+1, right
        Call compex(left,i)
      EndDo

      Do i=left+2, right

        j=i
        Call saveValue(i)

        Do While(compareValue(j-1))
          Call moveValue(j,j-1)
          j=j-1
        EndDo
     
        Call restoreValue(j)

      EndDo

    End Subroutine insertionSort

    Function partition(left, right) RESULT(i)

      Implicit NONE

      Integer, Intent(IN) :: left, right
      Integer :: i, j

      i = left - 1
      j = right

      Call saveValue(right)
   
      Loop1 : Do

        Loop2 : Do

          i = i+1
      
          If (i>right) EXIT Loop2
          If (compareValue(i)) EXIT Loop2
        
        EndDo Loop2

        Loop3: Do

          j = j-1
          If (.NOT. compareValue(j) .OR. j==left) EXIT LOOP3

        EndDo Loop3

        If (i>=j) EXIT Loop1
        Call exchange(i, j)

      EndDo Loop1

      Call exchange(i, right)

    End Function partition

    Recursive Subroutine quickSortIndex(left, right)

      Implicit NONE

      Integer, Intent(IN) :: left, right

      Integer :: i

      If ((right-left) > SWITCHSORTS) Then

        Call exchange((right+left)/2, (right-1))
        Call compex(left, right-1) 
        Call compex(right, left) 
        Call compex(right-1, right)
        i    = partition(left+1, right-1)
        Call quicksortIndex(left, i-1)
        Call quicksortIndex(i+1, right)
   
      EndIf

    End Subroutine quickSortIndex

  End Subroutine qsortIndexChar

! Sort Extensions of User_t user type

  Subroutine qsortInPlaceUser(a, left, right, order)

    Implicit NONE

    Integer,       Intent(IN)              :: left, right
    Class(User_t), Intent(INOUT), TARGET   :: a(left:right)
    Logical,       Intent(IN),    OPTIONAL :: order

    Integer                    :: savedIndex
    Logical                    :: sortorder
    Class(User_t), ALLOCATABLE :: savedUser
    Class(User_t), ALLOCATABLE :: tUser

    savedIndex = 0
    sortorder  = ASCENDING
    If (PRESENT(order)) sortorder = order

    Call quicksortIP(left, right)

    Call insertionSort

    If (ALLOCATED(tUser))     DEALLOCATE(tUser)
    If (ALLOCATED(savedUser)) DEALLOCATE(savedUser)

  Contains
      
    Subroutine exchange(i, j)

      Implicit NONE

      Integer, Intent(IN) :: i, j

      If (ALLOCATED(tUser)) Then
        Call tUser%assignValue(a(i)) 
      Else
        ALLOCATE(tUser, SOURCE=a(i))
      EndIf
      Call a(i)%assignValue(a(j))
      Call a(j)%assignValue(tUser)

    End Subroutine exchange

    Logical Function compare(i, j)

      Implicit NONE

      Integer, Intent(IN) :: i, j

      If (sortorder .EQV. ASCENDING) Then
        compare = (a(i) < a(j))
      Else
        compare = (a(i) > a(j))
      EndIf

    End Function compare

    Subroutine compex(i,j)

      Implicit NONE

      Integer, Intent(IN) :: i, j

      If (compare(j,i)) Call exchange(i,j)

    End Subroutine compex

    Subroutine moveValue(i, j)

      Implicit NONE

      Integer, Intent(IN) :: i, j

      Call a(i)%assignValue(a(j))

    End Subroutine moveValue

    Logical Function compareValue(j)

      Implicit NONE

      Integer, Intent(IN) :: j

      If (sortorder .EQV. ASCENDING) Then
        compareValue = (savedUser < a(j))
      Else
        compareValue = (savedUser > a(j))
      EndIf

    End Function compareValue

    Subroutine saveValue(i)

      Implicit NONE

      Integer, Intent(IN) :: i

      If (ALLOCATED(savedUser)) Then
        Call savedUser%assignValue(a(i))
      Else
        ALLOCATE(savedUser,SOURCE=a(i))
      EndIf

    End Subroutine saveValue

    Subroutine restoreValue(i)

      Implicit NONE

      Integer, Intent(IN) :: i

      Call a(i)%assignValue(savedUser)

    End Subroutine restoreValue

    Subroutine insertionSort()

      Implicit NONE

      Integer :: i, j

      Do i=left+1, right
        Call compex(left,i)
      EndDo

      Do i=left+2, right

        j=i
        Call saveValue(i)

        Do While(compareValue(j-1))
          Call moveValue(j,j-1)
          j=j-1
        EndDo
     
        Call restoreValue(j)

      EndDo

    End Subroutine insertionSort

    Function partition(left, right) RESULT(i)

      Implicit NONE

      Integer, Intent(IN) :: left, right
      Integer :: i, j

      i = left - 1
      j = right

      Call saveValue(right)
   
      Loop1 : Do

        Loop2 : Do

          i = i+1
      
          If (i>right) EXIT Loop2
          If (compareValue(i)) EXIT Loop2
        
        EndDo Loop2

        Loop3: Do

          j = j-1
          If (.NOT. compareValue(j) .OR. j==left) EXIT LOOP3

        EndDo Loop3

        If (i>=j) EXIT Loop1
        Call exchange(i, j)

      EndDo Loop1

      Call exchange(i, right)

    End Function partition

    Recursive Subroutine quickSortIP(left, right)

      Implicit NONE

      Integer, Intent(IN) :: left, right

      Integer :: i

      If ((right-left) > SWITCHSORTS) Then

        Call exchange((right+left)/2, (right-1))
        Call compex(left, right-1) 
        Call compex(right, left) 
        Call compex(right-1, right)
        i    = partition(left+1, right-1)
        Call quicksortIP(left, i-1)
        Call quicksortIP(i+1, right)
      EndIf

    End Subroutine quicksortIP

  End Subroutine qsortInPlaceUser

  Subroutine qsortIndexUser(a, left, right, permute, order)

    Implicit NONE

    Integer,        Intent(IN)             :: left, right
    Class(User_t),  Intent(INOUT), Target  :: a(left:right)
    Integer,        Intent(INOUT)          :: permute(left:right)
    Logical,        Intent(IN),   OPTIONAL :: order

    Integer                    :: j
    Integer                    :: savedIndex
    Logical                    :: sortorder
    Class(User_t), ALLOCATABLE :: savedUser

    savedIndex = 0
    sortorder = ASCENDING
    If (PRESENT(order)) sortorder = order

    Do j=left,right

      permute(j) = j

    EndDo 

    Call quicksortIndex(left, right)

    Call insertionSort

    If (ALLOCATED(savedUser)) DEALLOCATE(savedUser)

  Contains
 
    Subroutine exchange(i, j)

      Implicit NONE

      Integer, Intent(IN) :: i, j

      Integer :: tndx
      
      tndx = permute(i)
      permute(i) = permute(j)
      permute(j) = tndx

    End Subroutine exchange

    Logical Function compare(i, j)

      Implicit NONE

      Integer, Intent(IN) :: i, j

      Integer :: ii, ij

      ii = permute(i)
      ij = permute(j)

      If (sortorder .EQV. ASCENDING) Then
        compare = (a(ii) < a(ij))
      Else
        compare = (a(ii) > a(ij))
      EndIf

    End Function compare

    Subroutine compex(i,j)

      Implicit NONE

      Integer, Intent(IN) :: i, j

      If (compare(j,i)) Call exchange(i,j)

    End Subroutine compex

    Subroutine moveValue(i, j)

      Implicit NONE

      Integer, Intent(IN) :: i, j

      permute(i) = permute(j)

    End Subroutine moveValue

    Logical Function compareValue(j)

      Implicit NONE

      Integer, Intent(IN) :: j

      Integer :: ij

      ij = permute(j)
      If (sortorder .EQV. ASCENDING) Then
        compareValue = (savedUser < a(ij))
      Else
        compareValue = (savedUser > a(ij))
      EndIf

    End Function compareValue

    Subroutine saveValue(i)

      Implicit NONE

      Integer, Intent(IN) :: i
      
      savedIndex = permute(i)
      If(ALLOCATED(savedUser)) Then
        Call savedUser%assignValue(a(savedIndex))
      Else
        ALLOCATE(savedUser, SOURCE=a(savedIndex))
      EndIf
      
    End Subroutine saveValue

    Subroutine restoreValue(i)

      Implicit NONE

      Integer, Intent(IN) :: i

      permute(i) = savedIndex

    End Subroutine restoreValue

    Subroutine insertionSort()

      Implicit NONE

      Integer :: i, j

      Do i=left+1, right
        Call compex(left,i)
      EndDo

      Do i=left+2, right

        j=i
        Call saveValue(i)

        Do While(compareValue(j-1))
          Call moveValue(j,j-1)
          j=j-1
        EndDo
     
        Call restoreValue(j)

      EndDo

    End Subroutine insertionSort

    Function partition(left, right) RESULT(i)

      Implicit NONE

      Integer, Intent(IN) :: left, right
      Integer :: i, j

      i = left - 1
      j = right

      Call saveValue(right)
   
      Loop1 : Do

        Loop2 : Do

          i = i+1
      
          If (i>right) EXIT Loop2
          If (compareValue(i)) EXIT Loop2
        
        EndDo Loop2

        Loop3: Do

          j = j-1
          If (.NOT. compareValue(j) .OR. j==left) EXIT LOOP3

        EndDo Loop3

        If (i>=j) EXIT Loop1
        Call exchange(i, j)

      EndDo Loop1

      Call exchange(i, right)

    End Function partition

    Recursive Subroutine quickSortIndex(left, right)

      Implicit NONE

      Integer, Intent(IN) :: left, right

      Integer :: i

      If ((right-left) > SWITCHSORTS) Then

        Call exchange((right+left)/2, (right-1))
        Call compex(left, right-1) 
        Call compex(right, left) 
        Call compex(right-1, right)
        i    = partition(left+1, right-1)
        Call quicksortIndex(left, i-1)
        Call quicksortIndex(i+1, right)
   
      EndIf

    End Subroutine quickSortIndex

  End Subroutine qsortIndexUser

End Module quicksort