sort.F90

Go to the documentation of this file.

!! Copyright (C) 2002-2006 M. Marques, A. Castro, A. Rubio, G. Bertsch
!! Copyright (C) 2021 S. Ohlmann
!!
!! This program is free software; you can redistribute it and/or modify
!! it under the terms of the GNU General Public License as published by
!! the Free Software Foundation; either version 2, or (at your option)
!! any later version.
!!
!! This program is distributed in the hope that it will be useful,
!! but WITHOUT ANY WARRANTY; without even the implied warranty of
!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!! GNU General Public License for more details.
!!
!! You should have received a copy of the GNU General Public License
!! along with this program; if not, write to the Free Software
!! Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
!! 02110-1301, USA.
!!
 
#include "global.h"
 
module sort_oct_m
  use debug_oct_m
  use global_oct_m
 
  implicit none
 
  private
  public ::                     &
    sort,                       &
    robust_sort_by_abs
 
  interface robust_sort_by_abs
    procedure irobust_sort_by_abs, drobust_sort_by_abs
  end interface robust_sort_by_abs
 
  ! ------------------------------------------------------------------------------
  interface sort
    module procedure dsort, isort, lsort
    module procedure dshellsort1, zshellsort1, ishellsort1
    module procedure dshellsort2, zshellsort2, ishellsort2
  end interface sort
 
  interface
    subroutine dsort1(size, array)
      use, intrinsic :: iso_fortran_env
      implicit none
      integer, intent(in)    :: size
      real(real64),   intent(inout) :: array(*)
    end subroutine dsort1
 
    subroutine dsort2(size, array, indices)
      use, intrinsic :: iso_fortran_env
      implicit none
      integer, intent(in)    :: size
      real(real64),   intent(inout) :: array(*)
      integer, intent(out)   :: indices(*)
    end subroutine dsort2
 
    subroutine isort1(size, array)
      implicit none
      integer, intent(in)    :: size
      integer, intent(inout) :: array(*)
    end subroutine isort1
 
    subroutine isort2(size, array, indices)
      implicit none
      integer, intent(in)    :: size
      integer, intent(inout) :: array(*)
      integer, intent(out)   :: indices(*)
    end subroutine isort2
 
    subroutine lsort1(size, array)
      use, intrinsic :: iso_fortran_env
      implicit none
      integer,    intent(in)    :: size
      integer(int64), intent(inout) :: array(*)
    end subroutine lsort1
 
    subroutine lsort2(size, array, indices)
      use, intrinsic :: iso_fortran_env
      implicit none
      integer,    intent(in)    :: size
      integer(int64), intent(inout) :: array(*)
      integer,    intent(out)   :: indices(*)
    end subroutine lsort2
  end interface
 
contains
 
  ! ---------------------------------------------------------
  subroutine dsort(a, ind)
    real(real64),      intent(inout) :: a(:)
    integer, optional, intent(out)   :: ind(:)
 
    push_sub(dsort)
 
    if (size(a) > 0) then
 
      if (.not. present(ind)) then
        call dsort1(size(a), a)
      else
        call dsort2(size(a), a, ind)
      end if
 
    end if
 
    pop_sub(dsort)
  end subroutine dsort
 
 
  ! ---------------------------------------------------------
  subroutine isort(a, ind)
    integer,           intent(inout) :: a(:)
    integer, optional, intent(out)   :: ind(:)
 
    push_sub(isort)
 
    if (size(a) > 0) then
 
      if (.not. present(ind)) then
        call isort1(size(a), a)
      else
        call isort2(size(a), a, ind)
      end if
 
    end if
 
    pop_sub(isort)
  end subroutine isort
 
  ! ---------------------------------------------------------
  subroutine lsort(a, ind)
    integer(int64),       intent(inout) :: a(:)
    integer, optional, intent(out)   :: ind(:)
 
    push_sub(lsort)
 
    if (size(a) > 0) then
 
      if (.not. present(ind)) then
        call lsort1(size(a), a)
      else
        call lsort2(size(a), a, ind)
      end if
 
    end if
 
    pop_sub(lsort)
  end subroutine lsort
 
 
  pure logical function iless_idx(i, j, off, kabs, ksgn) result(less)
    integer, intent(in)     :: i, j
    integer, intent(in)     :: off(:, :)          ! (ndim, n)
    real(int64), intent(in) :: kabs(:)
    integer, intent(in)     :: ksgn(:)
 
    integer :: d, ndim
 
    if (are_different(kabs(i), kabs(j))) then
      less = kabs(i) < kabs(j)
      return
    end if
 
    if (ksgn(i) /= ksgn(j)) then
      less = ksgn(i) < ksgn(j)    ! 0 before 1
      return
    end if
 
    ndim = size(off, 1)
    do d = 1, ndim
      if (off(d,i) /= off(d,j)) then
        less = off(d,i) < off(d,j)
        return
      end if
    end do
 
    less = i < j                  ! final stability (should not matter if offsets are unique)
  end function iless_idx
 
  pure logical function dless_idx(i, j, off, kabs, ksgn) result(less)
    integer, intent(in)     :: i, j
    real(int64), intent(in) :: off(:, :)          ! (ndim, n)
    real(int64), intent(in) :: kabs(:)
    integer, intent(in)     :: ksgn(:)
 
    integer :: d, ndim
 
    if (are_different(kabs(i), kabs(j))) then
      less = kabs(i) < kabs(j)
      return
    end if
 
    if (ksgn(i) /= ksgn(j)) then
      less = ksgn(i) < ksgn(j)    ! 0 before 1
      return
    end if
 
    ndim = size(off, 1)
    do d = 1, ndim
      if (are_different(off(d,i), off(d,j))) then
        less = off(d,i) < off(d,j)
        return
      end if
    end do
 
    less = i < j                  ! final stability (should not matter if offsets are unique)
  end function dless_idx
 
  logical pure function are_different(x, y)
    real(real64), intent(in) :: x,y
    are_different = abs(x - y) > abs(y) * 1.0e-14_real64 + 1.0e-14_real64
  end function are_different
 
  recursive subroutine imergesort_perm(perm, tmp, l, r, off, kabs, ksgn)
    integer, intent(inout)   :: perm(:), tmp(:)
    integer, intent(in)      :: l, r
    integer, intent(in)      :: off(:, :)
    real(real64), intent(in) :: kabs(:)
    integer, intent(in)      :: ksgn(:)
 
    integer :: m, i, j, k
 
    if (r <= l) return
    m = (l + r) / 2
 
    call imergesort_perm(perm, tmp, l,   m, off, kabs, ksgn)
    call imergesort_perm(perm, tmp, m+1, r, off, kabs, ksgn)
 
    i = l; j = m+1; k = l
    do while (i <= m .and. j <= r)
      if (iless_idx(perm(i), perm(j), off, kabs, ksgn)) then
        tmp(k) = perm(i); i = i+1
      else
        tmp(k) = perm(j); j = j+1
      end if
      k = k+1
    end do
    do while (i <= m); tmp(k) = perm(i); i=i+1; k=k+1; end do
    do while (j <= r); tmp(k) = perm(j); j=j+1; k=k+1; end do
 
    perm(l:r) = tmp(l:r)
  end subroutine imergesort_perm
 
  recursive subroutine dmergesort_perm(perm, tmp, l, r, off, kabs, ksgn)
    integer, intent(inout)   :: perm(:), tmp(:)
    integer, intent(in)      :: l, r
    real(real64), intent(in) :: off(:, :)
    real(real64), intent(in) :: kabs(:)
    integer, intent(in)      :: ksgn(:)
 
    integer :: m, i, j, k
 
    if (r <= l) return
    m = (l + r) / 2
 
    call dmergesort_perm(perm, tmp, l,   m, off, kabs, ksgn)
    call dmergesort_perm(perm, tmp, m+1, r, off, kabs, ksgn)
 
    i = l; j = m+1; k = l
    do while (i <= m .and. j <= r)
      if (dless_idx(perm(i), perm(j), off, kabs, ksgn)) then
        tmp(k) = perm(i); i = i+1
      else
        tmp(k) = perm(j); j = j+1
      end if
      k = k+1
    end do
    do while (i <= m); tmp(k) = perm(i); i=i+1; k=k+1; end do
    do while (j <= r); tmp(k) = perm(j); j=j+1; k=k+1; end do
 
    perm(l:r) = tmp(l:r)
  end subroutine dmergesort_perm
 
 
  subroutine irobust_sort_by_abs(v, off, perm, negative_first)
    real(real64), intent(in) :: v(:)
    integer, intent(in)      :: off(:, :)              ! (ndim, n)
    integer, intent(out)     :: perm(size(v))
    logical, intent(in), optional :: negative_first
 
    integer :: n, i
    logical :: neg_first
    integer, allocatable :: tmp(:)
    real(real64), allocatable :: kabs(:)
    integer, allocatable :: ksgn(:)
 
    push_sub(irobust_sort_by_abs)
 
    n = size(v)
    assert(size(off, dim=2) == n)
 
    neg_first = optional_default(negative_first, .true.)
 
    allocate(tmp(n), kabs(n), ksgn(n))
 
    do i = 1, n
      perm(i) = i
      kabs(i) = abs(v(i))
 
      if (neg_first) then
        ksgn(i) = merge(0, 1, v(i) < 0.0_real64)   ! neg (0) then pos (1)
      else
        ksgn(i) = merge(0, 1, v(i) >= 0.0_real64)  ! pos (0) then neg (1)
      end if
    end do
 
    call imergesort_perm(perm, tmp, 1, n, off, kabs, ksgn)
 
    deallocate(tmp, kabs, ksgn)
    pop_sub(irobust_sort_by_abs)
  end subroutine irobust_sort_by_abs
 
  subroutine drobust_sort_by_abs(v, off, perm, negative_first)
    real(real64), intent(in) :: v(:)
    real(real64), intent(in) :: off(:, :)              ! (ndim, n)
    integer, intent(out)     :: perm(size(v))
    logical, intent(in), optional :: negative_first
 
    integer :: n, i
    logical :: neg_first
    integer, allocatable :: tmp(:)
    real(real64), allocatable :: kabs(:)
    integer, allocatable :: ksgn(:)
 
    push_sub(drobust_sort_by_abs)
 
    n = size(v)
    assert(size(off, dim=2) == n)
 
    neg_first = optional_default(negative_first, .true.)
 
    allocate(tmp(n), kabs(n), ksgn(n))
 
    do i = 1, n
      perm(i) = i
      kabs(i) = abs(v(i))
 
      if (neg_first) then
        ksgn(i) = merge(0, 1, v(i) < 0.0_real64)   ! neg (0) then pos (1)
      else
        ksgn(i) = merge(0, 1, v(i) >= 0.0_real64)  ! pos (0) then neg (1)
      end if
    end do
 
    call dmergesort_perm(perm, tmp, 1, n, off, kabs, ksgn)
 
    deallocate(tmp, kabs, ksgn)
    pop_sub(drobust_sort_by_abs)
  end subroutine drobust_sort_by_abs
 
 
#include "undef.F90"
#include "complex.F90"
#include "sort_inc.F90"
 
#include "undef.F90"
#include "real.F90"
#include "sort_inc.F90"
 
#include "undef.F90"
#include "integer.F90"
#include "sort_inc.F90"
 
end module sort_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: