From 744731b4ed91e774a4da6b2bc8691ad443569e43 Mon Sep 17 00:00:00 2001 From: Salvatore Filippone Date: Mon, 5 Mar 2018 17:12:54 +0000 Subject: [PATCH] Reworked IP_REORD. --- base/modules/Makefile | 6 +- base/modules/auxil/psb_c_ip_reord_mod.f90 | 7 +- base/modules/auxil/psb_c_sort_mod.f90 | 608 ++++++++++++++++++++++ base/modules/auxil/psb_d_ip_reord_mod.f90 | 7 +- base/modules/auxil/psb_d_sort_mod.f90 | 570 ++++++++++++++++++++ base/modules/auxil/psb_e_ip_reord_mod.f90 | 186 +++++++ base/modules/auxil/psb_i_sort_mod.f90 | 577 ++++++++++++++++++++ base/modules/auxil/psb_ip_reord_mod.f90 | 4 +- base/modules/auxil/psb_l_sort_mod.f90 | 577 ++++++++++++++++++++ base/modules/auxil/psb_m_ip_reord_mod.f90 | 186 +++++++ base/modules/auxil/psb_s_ip_reord_mod.f90 | 7 +- base/modules/auxil/psb_s_sort_mod.f90 | 570 ++++++++++++++++++++ base/modules/auxil/psb_z_ip_reord_mod.f90 | 7 +- base/modules/auxil/psb_z_sort_mod.f90 | 608 ++++++++++++++++++++++ base/serial/sort/psb_e_msort_impl.f90 | 3 +- base/serial/sort/psb_m_msort_impl.f90 | 1 - 16 files changed, 3901 insertions(+), 23 deletions(-) create mode 100644 base/modules/auxil/psb_c_sort_mod.f90 create mode 100644 base/modules/auxil/psb_d_sort_mod.f90 create mode 100644 base/modules/auxil/psb_e_ip_reord_mod.f90 create mode 100644 base/modules/auxil/psb_i_sort_mod.f90 create mode 100644 base/modules/auxil/psb_l_sort_mod.f90 create mode 100644 base/modules/auxil/psb_m_ip_reord_mod.f90 create mode 100644 base/modules/auxil/psb_s_sort_mod.f90 create mode 100644 base/modules/auxil/psb_z_sort_mod.f90 diff --git a/base/modules/Makefile b/base/modules/Makefile index 0aba1cb13..a603fbed0 100644 --- a/base/modules/Makefile +++ b/base/modules/Makefile @@ -55,7 +55,7 @@ UTIL_MODS = auxil/psb_string_mod.o desc/psb_desc_const_mod.o desc/psb_indx_map_m auxil/psi_c_serial_mod.o auxil/psi_z_serial_mod.o \ psi_mod.o psi_i_mod.o psi_s_mod.o psi_d_mod.o psi_c_mod.o psi_z_mod.o\ auxil/psb_ip_reord_mod.o\ - auxil/psb_i_ip_reord_mod.o auxil/psb_l_ip_reord_mod.o \ + auxil/psb_m_ip_reord_mod.o auxil/psb_e_ip_reord_mod.o \ auxil/psb_s_ip_reord_mod.o auxil/psb_d_ip_reord_mod.o \ auxil/psb_c_ip_reord_mod.o auxil/psb_z_ip_reord_mod.o \ auxil/psb_m_hsort_mod.o auxil/psb_m_isort_mod.o \ @@ -84,6 +84,7 @@ UTIL_MODS = auxil/psb_string_mod.o desc/psb_desc_const_mod.o desc/psb_indx_map_m serial/psb_z_base_mat_mod.o serial/psb_z_csr_mat_mod.o serial/psb_z_csc_mat_mod.o serial/psb_z_mat_mod.o + MODULES=$(BASIC_MODS) $(UTIL_MODS) OBJS = error.o psb_base_mod.o $(EXTRA_COBJS) cutil.o LIBDIR=.. @@ -173,10 +174,11 @@ auxil/psi_serial_mod.o: auxil/psi_i_serial_mod.o \ auxil/psi_s_serial_mod.o auxil/psi_d_serial_mod.o auxil/psi_c_serial_mod.o auxil/psi_z_serial_mod.o auxil/psi_i_serial_mod.o auxil/psi_s_serial_mod.o auxil/psi_d_serial_mod.o auxil/psi_c_serial_mod.o auxil/psi_z_serial_mod.o: psb_const_mod.o -auxil/psb_ip_reord_mod.o: auxil/psb_i_ip_reord_mod.o auxil/psb_l_ip_reord_mod.o \ +auxil/psb_ip_reord_mod.o: auxil/psb_m_ip_reord_mod.o auxil/psb_e_ip_reord_mod.o \ auxil/psb_s_ip_reord_mod.o auxil/psb_d_ip_reord_mod.o \ auxil/psb_c_ip_reord_mod.o auxil/psb_z_ip_reord_mod.o + serial/psb_base_mat_mod.o: auxil/psi_serial_mod.o serial/psb_s_base_mat_mod.o serial/psb_d_base_mat_mod.o serial/psb_c_base_mat_mod.o serial/psb_z_base_mat_mod.o: serial/psb_base_mat_mod.o serial/psb_s_base_mat_mod.o: serial/psb_s_base_vect_mod.o diff --git a/base/modules/auxil/psb_c_ip_reord_mod.f90 b/base/modules/auxil/psb_c_ip_reord_mod.f90 index c31643fbf..9c58b9ace 100644 --- a/base/modules/auxil/psb_c_ip_reord_mod.f90 +++ b/base/modules/auxil/psb_c_ip_reord_mod.f90 @@ -53,10 +53,9 @@ contains integer(psb_ipk_), intent(in) :: n integer(psb_ipk_) :: iaux(0:*) complex(psb_spk_) :: x(*) - integer(psb_ipk_) :: lswap, lp, k complex(psb_spk_) :: swap - + lp = iaux(0) k = 1 do @@ -82,10 +81,11 @@ contains integer(psb_ipk_) :: iaux(0:*) complex(psb_spk_) :: x(*) integer(psb_ipk_) :: indx(*) - integer(psb_ipk_) :: lswap, lp, k, ixswap complex(psb_spk_) :: swap + + lp = iaux(0) k = 1 do @@ -150,7 +150,6 @@ contains integer(psb_ipk_) :: iaux(0:*) complex(psb_spk_) :: x(*) integer(psb_ipk_) :: i1(*), i2(*), i3(*) - integer(psb_ipk_) :: lswap, lp, k, isw1, isw2, isw3 complex(psb_spk_) :: swap diff --git a/base/modules/auxil/psb_c_sort_mod.f90 b/base/modules/auxil/psb_c_sort_mod.f90 new file mode 100644 index 000000000..a2f0f3c6b --- /dev/null +++ b/base/modules/auxil/psb_c_sort_mod.f90 @@ -0,0 +1,608 @@ +! +! Parallel Sparse BLAS version 3.5 +! (C) Copyright 2006-2018 +! Salvatore Filippone +! Alfredo Buttari +! +! 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. +! 3. The name of the PSBLAS group or the names of its contributors may +! not be used to endorse or promote products derived from this +! software without specific written permission. +! +! 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 PSBLAS GROUP OR ITS 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. +! +! +! +! Sorting routines +! References: +! D. Knuth +! The Art of Computer Programming, vol. 3 +! Addison-Wesley +! +! Aho, Hopcroft, Ullman +! Data Structures and Algorithms +! Addison-Wesley +! +module psb_c_sort_mod + use psb_const_mod + + + interface psb_msort_unique + subroutine psb_cmsort_u(x,nout,dir) + import + complex(psb_spk_), intent(inout) :: x(:) + integer(psb_ipk_), intent(out) :: nout + integer(psb_ipk_), optional, intent(in) :: dir + end subroutine psb_cmsort_u + end interface psb_msort_unique + + type psb_c_heap + integer(psb_ipk_) :: last, dir + complex(psb_spk_), allocatable :: keys(:) + contains + procedure, pass(heap) :: init => psb_c_init_heap + procedure, pass(heap) :: howmany => psb_c_howmany + procedure, pass(heap) :: insert => psb_c_insert_heap + procedure, pass(heap) :: get_first => psb_c_heap_get_first + procedure, pass(heap) :: dump => psb_c_dump_heap + procedure, pass(heap) :: free => psb_c_free_heap + end type psb_c_heap + + type psb_c_idx_heap + integer(psb_ipk_) :: last, dir + complex(psb_spk_), allocatable :: keys(:) + integer(psb_ipk_), allocatable :: idxs(:) + contains + procedure, pass(heap) :: init => psb_c_idx_init_heap + procedure, pass(heap) :: howmany => psb_c_idx_howmany + procedure, pass(heap) :: insert => psb_c_idx_insert_heap + procedure, pass(heap) :: get_first => psb_c_idx_heap_get_first + procedure, pass(heap) :: dump => psb_c_idx_dump_heap + procedure, pass(heap) :: free => psb_c_idx_free_heap + end type psb_c_idx_heap + + + interface psb_msort + subroutine psb_cmsort(x,ix,dir,flag) + import + complex(psb_spk_), intent(inout) :: x(:) + integer(psb_ipk_), optional, intent(in) :: dir, flag + integer(psb_ipk_), optional, intent(inout) :: ix(:) + end subroutine psb_cmsort + end interface psb_msort + + interface + subroutine psi_c_lmsort_up(n,k,l,iret) + import + implicit none + integer(psb_ipk_) :: n, iret + complex(psb_spk_) :: k(n) + integer(psb_ipk_) :: l(0:n+1) + end subroutine psi_c_lmsort_up + subroutine psi_c_lmsort_dw(n,k,l,iret) + import + implicit none + integer(psb_ipk_) :: n, iret + complex(psb_spk_) :: k(n) + integer(psb_ipk_) :: l(0:n+1) + end subroutine psi_c_lmsort_dw + subroutine psi_c_almsort_up(n,k,l,iret) + import + implicit none + integer(psb_ipk_) :: n, iret + complex(psb_spk_) :: k(n) + integer(psb_ipk_) :: l(0:n+1) + end subroutine psi_c_almsort_up + subroutine psi_c_almsort_dw(n,k,l,iret) + import + implicit none + integer(psb_ipk_) :: n, iret + complex(psb_spk_) :: k(n) + integer(psb_ipk_) :: l(0:n+1) + end subroutine psi_c_almsort_dw + end interface + interface + subroutine psi_c_amsort_up(n,k,l,iret) + import + implicit none + integer(psb_ipk_) :: n, iret + complex(psb_spk_) :: k(n) + integer(psb_ipk_) :: l(0:n+1) + end subroutine psi_c_amsort_up + subroutine psi_c_amsort_dw(n,k,l,iret) + import + implicit none + integer(psb_ipk_) :: n, iret + complex(psb_spk_) :: k(n) + integer(psb_ipk_) :: l(0:n+1) + end subroutine psi_c_amsort_dw + end interface + + + interface psb_qsort + subroutine psb_cqsort(x,ix,dir,flag) + import + complex(psb_spk_), intent(inout) :: x(:) + integer(psb_ipk_), optional, intent(in) :: dir, flag + integer(psb_ipk_), optional, intent(inout) :: ix(:) + end subroutine psb_cqsort + end interface psb_qsort + + interface psb_isort + subroutine psb_cisort(x,ix,dir,flag) + import + complex(psb_spk_), intent(inout) :: x(:) + integer(psb_ipk_), optional, intent(in) :: dir, flag + integer(psb_ipk_), optional, intent(inout) :: ix(:) + end subroutine psb_cisort + end interface psb_isort + + + interface psb_hsort + subroutine psb_chsort(x,ix,dir,flag) + import + complex(psb_spk_), intent(inout) :: x(:) + integer(psb_ipk_), optional, intent(in) :: dir, flag + integer(psb_ipk_), optional, intent(inout) :: ix(:) + end subroutine psb_chsort + end interface psb_hsort + + + interface + subroutine psi_c_insert_heap(key,last,heap,dir,info) + import + implicit none + + ! + ! Input: + ! key: the new value + ! last: pointer to the last occupied element in heap + ! heap: the heap + ! dir: sorting direction + + complex(psb_spk_), intent(in) :: key + complex(psb_spk_), intent(inout) :: heap(:) + integer(psb_ipk_), intent(in) :: dir + integer(psb_ipk_), intent(inout) :: last + integer(psb_ipk_), intent(out) :: info + end subroutine psi_c_insert_heap + end interface + + interface + subroutine psi_c_idx_insert_heap(key,index,last,heap,idxs,dir,info) + import + implicit none + + ! + ! Input: + ! key: the new value + ! last: pointer to the last occupied element in heap + ! heap: the heap + ! dir: sorting direction + + complex(psb_spk_), intent(in) :: key + complex(psb_spk_), intent(inout) :: heap(:) + integer(psb_ipk_), intent(in) :: index + integer(psb_ipk_), intent(in) :: dir + integer(psb_ipk_), intent(inout) :: idxs(:) + integer(psb_ipk_), intent(inout) :: last + integer(psb_ipk_), intent(out) :: info + end subroutine psi_c_idx_insert_heap + end interface + + + interface + subroutine psi_c_heap_get_first(key,last,heap,dir,info) + import + implicit none + complex(psb_spk_), intent(inout) :: key + integer(psb_ipk_), intent(inout) :: last + integer(psb_ipk_), intent(in) :: dir + complex(psb_spk_), intent(inout) :: heap(:) + integer(psb_ipk_), intent(out) :: info + end subroutine psi_c_heap_get_first + end interface + + interface + subroutine psi_c_idx_heap_get_first(key,index,last,heap,idxs,dir,info) + import + complex(psb_spk_), intent(inout) :: key + integer(psb_ipk_), intent(out) :: index + complex(psb_spk_), intent(inout) :: heap(:) + integer(psb_ipk_), intent(in) :: dir + integer(psb_ipk_), intent(inout) :: last + integer(psb_ipk_), intent(inout) :: idxs(:) + integer(psb_ipk_), intent(out) :: info + end subroutine psi_c_idx_heap_get_first + end interface + + interface + subroutine psi_clisrx_up(n,x,ix) + import + complex(psb_spk_), intent(inout) :: x(:) + integer(psb_ipk_), intent(inout) :: ix(:) + integer(psb_ipk_), intent(in) :: n + end subroutine psi_clisrx_up + subroutine psi_clisrx_dw(n,x,ix) + import + complex(psb_spk_), intent(inout) :: x(:) + integer(psb_ipk_), intent(inout) :: ix(:) + integer(psb_ipk_), intent(in) :: n + end subroutine psi_clisrx_dw + subroutine psi_clisr_up(n,x) + import + complex(psb_spk_), intent(inout) :: x(:) + integer(psb_ipk_), intent(in) :: n + end subroutine psi_clisr_up + subroutine psi_clisr_dw(n,x) + import + complex(psb_spk_), intent(inout) :: x(:) + integer(psb_ipk_), intent(in) :: n + end subroutine psi_clisr_dw + subroutine psi_calisrx_up(n,x,ix) + import + complex(psb_spk_), intent(inout) :: x(:) + integer(psb_ipk_), intent(inout) :: ix(:) + integer(psb_ipk_), intent(in) :: n + end subroutine psi_calisrx_up + subroutine psi_calisrx_dw(n,x,ix) + import + complex(psb_spk_), intent(inout) :: x(:) + integer(psb_ipk_), intent(inout) :: ix(:) + integer(psb_ipk_), intent(in) :: n + end subroutine psi_calisrx_dw + subroutine psi_calisr_up(n,x) + import + complex(psb_spk_), intent(inout) :: x(:) + integer(psb_ipk_), intent(in) :: n + end subroutine psi_calisr_up + subroutine psi_calisr_dw(n,x) + import + complex(psb_spk_), intent(inout) :: x(:) + integer(psb_ipk_), intent(in) :: n + end subroutine psi_calisr_dw + subroutine psi_caisrx_up(n,x,ix) + import + complex(psb_spk_), intent(inout) :: x(:) + integer(psb_ipk_), intent(inout) :: ix(:) + integer(psb_ipk_), intent(in) :: n + end subroutine psi_caisrx_up + subroutine psi_caisrx_dw(n,x,ix) + import + complex(psb_spk_), intent(inout) :: x(:) + integer(psb_ipk_), intent(inout) :: ix(:) + integer(psb_ipk_), intent(in) :: n + end subroutine psi_caisrx_dw + subroutine psi_caisr_up(n,x) + import + complex(psb_spk_), intent(inout) :: x(:) + integer(psb_ipk_), intent(in) :: n + end subroutine psi_caisr_up + subroutine psi_caisr_dw(n,x) + import + complex(psb_spk_), intent(inout) :: x(:) + integer(psb_ipk_), intent(in) :: n + end subroutine psi_caisr_dw + end interface + + interface + subroutine psi_clqsrx_up(n,x,ix) + import + complex(psb_spk_), intent(inout) :: x(:) + integer(psb_ipk_), intent(inout) :: ix(:) + integer(psb_ipk_), intent(in) :: n + end subroutine psi_clqsrx_up + subroutine psi_clqsrx_dw(n,x,ix) + import + complex(psb_spk_), intent(inout) :: x(:) + integer(psb_ipk_), intent(inout) :: ix(:) + integer(psb_ipk_), intent(in) :: n + end subroutine psi_clqsrx_dw + subroutine psi_clqsr_up(n,x) + import + complex(psb_spk_), intent(inout) :: x(:) + integer(psb_ipk_), intent(in) :: n + end subroutine psi_clqsr_up + subroutine psi_clqsr_dw(n,x) + import + complex(psb_spk_), intent(inout) :: x(:) + integer(psb_ipk_), intent(in) :: n + end subroutine psi_clqsr_dw + subroutine psi_calqsrx_up(n,x,ix) + import + complex(psb_spk_), intent(inout) :: x(:) + integer(psb_ipk_), intent(inout) :: ix(:) + integer(psb_ipk_), intent(in) :: n + end subroutine psi_calqsrx_up + subroutine psi_calqsrx_dw(n,x,ix) + import + complex(psb_spk_), intent(inout) :: x(:) + integer(psb_ipk_), intent(inout) :: ix(:) + integer(psb_ipk_), intent(in) :: n + end subroutine psi_calqsrx_dw + subroutine psi_calqsr_up(n,x) + import + complex(psb_spk_), intent(inout) :: x(:) + integer(psb_ipk_), intent(in) :: n + end subroutine psi_calqsr_up + subroutine psi_calqsr_dw(n,x) + import + complex(psb_spk_), intent(inout) :: x(:) + integer(psb_ipk_), intent(in) :: n + end subroutine psi_calqsr_dw + subroutine psi_caqsrx_up(n,x,ix) + import + complex(psb_spk_), intent(inout) :: x(:) + integer(psb_ipk_), intent(inout) :: ix(:) + integer(psb_ipk_), intent(in) :: n + end subroutine psi_caqsrx_up + subroutine psi_caqsrx_dw(n,x,ix) + import + complex(psb_spk_), intent(inout) :: x(:) + integer(psb_ipk_), intent(inout) :: ix(:) + integer(psb_ipk_), intent(in) :: n + end subroutine psi_caqsrx_dw + subroutine psi_caqsr_up(n,x) + import + complex(psb_spk_), intent(inout) :: x(:) + integer(psb_ipk_), intent(in) :: n + end subroutine psi_caqsr_up + subroutine psi_caqsr_dw(n,x) + import + complex(psb_spk_), intent(inout) :: x(:) + integer(psb_ipk_), intent(in) :: n + end subroutine psi_caqsr_dw + end interface + +contains + + subroutine psb_c_init_heap(heap,info,dir) + use psb_realloc_mod, only : psb_ensure_size + implicit none + class(psb_c_heap), intent(inout) :: heap + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_), intent(in), optional :: dir + + info = psb_success_ + heap%last=0 + if (present(dir)) then + heap%dir = dir + else + heap%dir = psb_asort_up_ + endif + select case(heap%dir) + case (psb_asort_up_,psb_asort_down_) + ! ok, do nothing + case default + write(psb_err_unit,*) 'Invalid direction, defaulting to psb_asort_up_' + heap%dir = psb_asort_up_ + end select + call psb_ensure_size(psb_heap_resize,heap%keys,info) + + return + end subroutine psb_c_init_heap + + + function psb_c_howmany(heap) result(res) + implicit none + class(psb_c_heap), intent(in) :: heap + integer(psb_ipk_) :: res + res = heap%last + end function psb_c_howmany + + subroutine psb_c_insert_heap(key,heap,info) + use psb_realloc_mod, only : psb_ensure_size + implicit none + + complex(psb_spk_), intent(in) :: key + class(psb_c_heap), intent(inout) :: heap + integer(psb_ipk_), intent(out) :: info + + info = psb_success_ + if (heap%last < 0) then + write(psb_err_unit,*) 'Invalid last in heap ',heap%last + info = heap%last + return + endif + + call psb_ensure_size(heap%last+1,heap%keys,info,addsz=psb_heap_resize) + if (info /= psb_success_) then + write(psb_err_unit,*) 'Memory allocation failure in heap_insert' + info = -5 + return + end if + call psi_c_insert_heap(key,& + & heap%last,heap%keys,heap%dir,info) + + return + end subroutine psb_c_insert_heap + + subroutine psb_c_heap_get_first(key,heap,info) + implicit none + + class(psb_c_heap), intent(inout) :: heap + integer(psb_ipk_), intent(out) :: info + complex(psb_spk_), intent(out) :: key + + + info = psb_success_ + + call psi_c_heap_get_first(key,& + & heap%last,heap%keys,heap%dir,info) + + return + end subroutine psb_c_heap_get_first + + subroutine psb_c_dump_heap(iout,heap,info) + + implicit none + class(psb_c_heap), intent(in) :: heap + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_), intent(in) :: iout + + info = psb_success_ + if (iout < 0) then + write(psb_err_unit,*) 'Invalid file ' + info =-1 + return + end if + + write(iout,*) 'Heap direction ',heap%dir + write(iout,*) 'Heap size ',heap%last + if ((heap%last > 0).and.((.not.allocated(heap%keys)).or.& + & (size(heap%keys) 0).and.((.not.allocated(heap%keys)).or.& + & (size(heap%keys) 0).and.((.not.allocated(heap%idxs)).or.& + & (size(heap%idxs) psb_d_init_heap + procedure, pass(heap) :: howmany => psb_d_howmany + procedure, pass(heap) :: insert => psb_d_insert_heap + procedure, pass(heap) :: get_first => psb_d_heap_get_first + procedure, pass(heap) :: dump => psb_d_dump_heap + procedure, pass(heap) :: free => psb_d_free_heap + end type psb_d_heap + + type psb_d_idx_heap + integer(psb_ipk_) :: last, dir + real(psb_dpk_), allocatable :: keys(:) + integer(psb_ipk_), allocatable :: idxs(:) + contains + procedure, pass(heap) :: init => psb_d_idx_init_heap + procedure, pass(heap) :: howmany => psb_d_idx_howmany + procedure, pass(heap) :: insert => psb_d_idx_insert_heap + procedure, pass(heap) :: get_first => psb_d_idx_heap_get_first + procedure, pass(heap) :: dump => psb_d_idx_dump_heap + procedure, pass(heap) :: free => psb_d_idx_free_heap + end type psb_d_idx_heap + + + interface psb_msort + subroutine psb_dmsort(x,ix,dir,flag) + import + real(psb_dpk_), intent(inout) :: x(:) + integer(psb_ipk_), optional, intent(in) :: dir, flag + integer(psb_ipk_), optional, intent(inout) :: ix(:) + end subroutine psb_dmsort + end interface psb_msort + + + interface psb_bsrch + function psb_dbsrch(key,n,v) result(ipos) + import + integer(psb_ipk_) :: ipos, n + real(psb_dpk_) :: key + real(psb_dpk_) :: v(:) + end function psb_dbsrch + end interface psb_bsrch + + interface psb_ssrch + function psb_dssrch(key,n,v) result(ipos) + import + implicit none + integer(psb_ipk_) :: ipos, n + real(psb_dpk_) :: key + real(psb_dpk_) :: v(:) + end function psb_dssrch + end interface psb_ssrch + + interface + subroutine psi_d_msort_up(n,k,l,iret) + import + implicit none + integer(psb_ipk_) :: n, iret + real(psb_dpk_) :: k(n) + integer(psb_ipk_) :: l(0:n+1) + end subroutine psi_d_msort_up + subroutine psi_d_msort_dw(n,k,l,iret) + import + implicit none + integer(psb_ipk_) :: n, iret + real(psb_dpk_) :: k(n) + integer(psb_ipk_) :: l(0:n+1) + end subroutine psi_d_msort_dw + end interface + interface + subroutine psi_d_amsort_up(n,k,l,iret) + import + implicit none + integer(psb_ipk_) :: n, iret + real(psb_dpk_) :: k(n) + integer(psb_ipk_) :: l(0:n+1) + end subroutine psi_d_amsort_up + subroutine psi_d_amsort_dw(n,k,l,iret) + import + implicit none + integer(psb_ipk_) :: n, iret + real(psb_dpk_) :: k(n) + integer(psb_ipk_) :: l(0:n+1) + end subroutine psi_d_amsort_dw + end interface + + + interface psb_qsort + subroutine psb_dqsort(x,ix,dir,flag) + import + real(psb_dpk_), intent(inout) :: x(:) + integer(psb_ipk_), optional, intent(in) :: dir, flag + integer(psb_ipk_), optional, intent(inout) :: ix(:) + end subroutine psb_dqsort + end interface psb_qsort + + interface psb_isort + subroutine psb_disort(x,ix,dir,flag) + import + real(psb_dpk_), intent(inout) :: x(:) + integer(psb_ipk_), optional, intent(in) :: dir, flag + integer(psb_ipk_), optional, intent(inout) :: ix(:) + end subroutine psb_disort + end interface psb_isort + + + interface psb_hsort + subroutine psb_dhsort(x,ix,dir,flag) + import + real(psb_dpk_), intent(inout) :: x(:) + integer(psb_ipk_), optional, intent(in) :: dir, flag + integer(psb_ipk_), optional, intent(inout) :: ix(:) + end subroutine psb_dhsort + end interface psb_hsort + + + interface + subroutine psi_d_insert_heap(key,last,heap,dir,info) + import + implicit none + + ! + ! Input: + ! key: the new value + ! last: pointer to the last occupied element in heap + ! heap: the heap + ! dir: sorting direction + + real(psb_dpk_), intent(in) :: key + real(psb_dpk_), intent(inout) :: heap(:) + integer(psb_ipk_), intent(in) :: dir + integer(psb_ipk_), intent(inout) :: last + integer(psb_ipk_), intent(out) :: info + end subroutine psi_d_insert_heap + end interface + + interface + subroutine psi_d_idx_insert_heap(key,index,last,heap,idxs,dir,info) + import + implicit none + + ! + ! Input: + ! key: the new value + ! last: pointer to the last occupied element in heap + ! heap: the heap + ! dir: sorting direction + + real(psb_dpk_), intent(in) :: key + real(psb_dpk_), intent(inout) :: heap(:) + integer(psb_ipk_), intent(in) :: index + integer(psb_ipk_), intent(in) :: dir + integer(psb_ipk_), intent(inout) :: idxs(:) + integer(psb_ipk_), intent(inout) :: last + integer(psb_ipk_), intent(out) :: info + end subroutine psi_d_idx_insert_heap + end interface + + + interface + subroutine psi_d_heap_get_first(key,last,heap,dir,info) + import + implicit none + real(psb_dpk_), intent(inout) :: key + integer(psb_ipk_), intent(inout) :: last + integer(psb_ipk_), intent(in) :: dir + real(psb_dpk_), intent(inout) :: heap(:) + integer(psb_ipk_), intent(out) :: info + end subroutine psi_d_heap_get_first + end interface + + interface + subroutine psi_d_idx_heap_get_first(key,index,last,heap,idxs,dir,info) + import + real(psb_dpk_), intent(inout) :: key + integer(psb_ipk_), intent(out) :: index + real(psb_dpk_), intent(inout) :: heap(:) + integer(psb_ipk_), intent(in) :: dir + integer(psb_ipk_), intent(inout) :: last + integer(psb_ipk_), intent(inout) :: idxs(:) + integer(psb_ipk_), intent(out) :: info + end subroutine psi_d_idx_heap_get_first + end interface + + interface + subroutine psi_disrx_up(n,x,ix) + import + real(psb_dpk_), intent(inout) :: x(:) + integer(psb_ipk_), intent(inout) :: ix(:) + integer(psb_ipk_), intent(in) :: n + end subroutine psi_disrx_up + subroutine psi_disrx_dw(n,x,ix) + import + real(psb_dpk_), intent(inout) :: x(:) + integer(psb_ipk_), intent(inout) :: ix(:) + integer(psb_ipk_), intent(in) :: n + end subroutine psi_disrx_dw + subroutine psi_disr_up(n,x) + import + real(psb_dpk_), intent(inout) :: x(:) + integer(psb_ipk_), intent(in) :: n + end subroutine psi_disr_up + subroutine psi_disr_dw(n,x) + import + real(psb_dpk_), intent(inout) :: x(:) + integer(psb_ipk_), intent(in) :: n + end subroutine psi_disr_dw + subroutine psi_daisrx_up(n,x,ix) + import + real(psb_dpk_), intent(inout) :: x(:) + integer(psb_ipk_), intent(inout) :: ix(:) + integer(psb_ipk_), intent(in) :: n + end subroutine psi_daisrx_up + subroutine psi_daisrx_dw(n,x,ix) + import + real(psb_dpk_), intent(inout) :: x(:) + integer(psb_ipk_), intent(inout) :: ix(:) + integer(psb_ipk_), intent(in) :: n + end subroutine psi_daisrx_dw + subroutine psi_daisr_up(n,x) + import + real(psb_dpk_), intent(inout) :: x(:) + integer(psb_ipk_), intent(in) :: n + end subroutine psi_daisr_up + subroutine psi_daisr_dw(n,x) + import + real(psb_dpk_), intent(inout) :: x(:) + integer(psb_ipk_), intent(in) :: n + end subroutine psi_daisr_dw + end interface + + interface + subroutine psi_dqsrx_up(n,x,ix) + import + real(psb_dpk_), intent(inout) :: x(:) + integer(psb_ipk_), intent(inout) :: ix(:) + integer(psb_ipk_), intent(in) :: n + end subroutine psi_dqsrx_up + subroutine psi_dqsrx_dw(n,x,ix) + import + real(psb_dpk_), intent(inout) :: x(:) + integer(psb_ipk_), intent(inout) :: ix(:) + integer(psb_ipk_), intent(in) :: n + end subroutine psi_dqsrx_dw + subroutine psi_dqsr_up(n,x) + import + real(psb_dpk_), intent(inout) :: x(:) + integer(psb_ipk_), intent(in) :: n + end subroutine psi_dqsr_up + subroutine psi_dqsr_dw(n,x) + import + real(psb_dpk_), intent(inout) :: x(:) + integer(psb_ipk_), intent(in) :: n + end subroutine psi_dqsr_dw + subroutine psi_daqsrx_up(n,x,ix) + import + real(psb_dpk_), intent(inout) :: x(:) + integer(psb_ipk_), intent(inout) :: ix(:) + integer(psb_ipk_), intent(in) :: n + end subroutine psi_daqsrx_up + subroutine psi_daqsrx_dw(n,x,ix) + import + real(psb_dpk_), intent(inout) :: x(:) + integer(psb_ipk_), intent(inout) :: ix(:) + integer(psb_ipk_), intent(in) :: n + end subroutine psi_daqsrx_dw + subroutine psi_daqsr_up(n,x) + import + real(psb_dpk_), intent(inout) :: x(:) + integer(psb_ipk_), intent(in) :: n + end subroutine psi_daqsr_up + subroutine psi_daqsr_dw(n,x) + import + real(psb_dpk_), intent(inout) :: x(:) + integer(psb_ipk_), intent(in) :: n + end subroutine psi_daqsr_dw + end interface + +contains + + subroutine psb_d_init_heap(heap,info,dir) + use psb_realloc_mod, only : psb_ensure_size + implicit none + class(psb_d_heap), intent(inout) :: heap + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_), intent(in), optional :: dir + + info = psb_success_ + heap%last=0 + if (present(dir)) then + heap%dir = dir + else + heap%dir = psb_sort_up_ + endif + select case(heap%dir) + case (psb_sort_up_,psb_sort_down_,psb_asort_up_,psb_asort_down_) + ! ok, do nothing + case default + write(psb_err_unit,*) 'Invalid direction, defaulting to psb_sort_up_' + heap%dir = psb_sort_up_ + end select + call psb_ensure_size(psb_heap_resize,heap%keys,info) + + return + end subroutine psb_d_init_heap + + + function psb_d_howmany(heap) result(res) + implicit none + class(psb_d_heap), intent(in) :: heap + integer(psb_ipk_) :: res + res = heap%last + end function psb_d_howmany + + subroutine psb_d_insert_heap(key,heap,info) + use psb_realloc_mod, only : psb_ensure_size + implicit none + + real(psb_dpk_), intent(in) :: key + class(psb_d_heap), intent(inout) :: heap + integer(psb_ipk_), intent(out) :: info + + info = psb_success_ + if (heap%last < 0) then + write(psb_err_unit,*) 'Invalid last in heap ',heap%last + info = heap%last + return + endif + + call psb_ensure_size(heap%last+1,heap%keys,info,addsz=psb_heap_resize) + if (info /= psb_success_) then + write(psb_err_unit,*) 'Memory allocation failure in heap_insert' + info = -5 + return + end if + call psi_d_insert_heap(key,& + & heap%last,heap%keys,heap%dir,info) + + return + end subroutine psb_d_insert_heap + + subroutine psb_d_heap_get_first(key,heap,info) + implicit none + + class(psb_d_heap), intent(inout) :: heap + integer(psb_ipk_), intent(out) :: info + real(psb_dpk_), intent(out) :: key + + + info = psb_success_ + + call psi_d_heap_get_first(key,& + & heap%last,heap%keys,heap%dir,info) + + return + end subroutine psb_d_heap_get_first + + subroutine psb_d_dump_heap(iout,heap,info) + + implicit none + class(psb_d_heap), intent(in) :: heap + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_), intent(in) :: iout + + info = psb_success_ + if (iout < 0) then + write(psb_err_unit,*) 'Invalid file ' + info =-1 + return + end if + + write(iout,*) 'Heap direction ',heap%dir + write(iout,*) 'Heap size ',heap%last + if ((heap%last > 0).and.((.not.allocated(heap%keys)).or.& + & (size(heap%keys) 0).and.((.not.allocated(heap%keys)).or.& + & (size(heap%keys) 0).and.((.not.allocated(heap%idxs)).or.& + & (size(heap%idxs)n)) exit + do + if (lp >= k) exit + lp = iaux(lp) + end do + swap = x(lp) + x(lp) = x(k) + x(k) = swap + lswap = iaux(lp) + iaux(lp) = iaux(k) + iaux(k) = lp + lp = lswap + k = k + 1 + enddo + return + end subroutine psb_ip_reord_e1 + + subroutine psb_ip_reord_e1i1(n,x,indx,iaux) + integer(psb_ipk_), intent(in) :: n + integer(psb_epk_) :: iaux(0:*) + integer(psb_epk_) :: x(*) + integer(psb_epk_) :: indx(*) + integer(psb_epk_) :: lswap, lp, k, ixswap + integer(psb_epk_) :: swap + + + + lp = iaux(0) + k = 1 + do + if ((lp == 0).or.(k>n)) exit + do + if (lp >= k) exit + lp = iaux(lp) + end do + swap = x(lp) + x(lp) = x(k) + x(k) = swap + ixswap = indx(lp) + indx(lp) = indx(k) + indx(k) = ixswap + lswap = iaux(lp) + iaux(lp) = iaux(k) + iaux(k) = lp + lp = lswap + k = k + 1 + enddo + return + end subroutine psb_ip_reord_e1i1 + + subroutine psb_ip_reord_e1i2(n,x,i1,i2,iaux) + integer(psb_ipk_), intent(in) :: n + integer(psb_epk_) :: iaux(0:*) + integer(psb_epk_) :: x(*) + integer(psb_epk_) :: i1(*), i2(*) + + + integer(psb_epk_) :: lswap, lp, k, isw1, isw2 + integer(psb_epk_) :: swap + + lp = iaux(0) + k = 1 + do + if ((lp == 0).or.(k>n)) exit + do + if (lp >= k) exit + lp = iaux(lp) + end do + swap = x(lp) + x(lp) = x(k) + x(k) = swap + isw1 = i1(lp) + i1(lp) = i1(k) + i1(k) = isw1 + isw2 = i2(lp) + i2(lp) = i2(k) + i2(k) = isw2 + lswap = iaux(lp) + iaux(lp) = iaux(k) + iaux(k) = lp + lp = lswap + k = k + 1 + enddo + return + end subroutine psb_ip_reord_e1i2 + + subroutine psb_ip_reord_e1i3(n,x,i1,i2,i3,iaux) + integer(psb_ipk_), intent(in) :: n + integer(psb_epk_) :: iaux(0:*) + integer(psb_epk_) :: x(*) + integer(psb_epk_) :: i1(*), i2(*), i3(*) + + integer(psb_epk_) :: lswap, lp, k, isw1, isw2, isw3 + integer(psb_epk_) :: swap + + lp = iaux(0) + k = 1 + do + if ((lp == 0).or.(k>n)) exit + do + if (lp >= k) exit + lp = iaux(lp) + end do + swap = x(lp) + x(lp) = x(k) + x(k) = swap + isw1 = i1(lp) + i1(lp) = i1(k) + i1(k) = isw1 + isw2 = i2(lp) + i2(lp) = i2(k) + i2(k) = isw2 + isw3 = i3(lp) + i3(lp) = i3(k) + i3(k) = isw3 + lswap = iaux(lp) + iaux(lp) = iaux(k) + iaux(k) = lp + lp = lswap + k = k + 1 + enddo + return + end subroutine psb_ip_reord_e1i3 + +end module psb_e_ip_reord_mod diff --git a/base/modules/auxil/psb_i_sort_mod.f90 b/base/modules/auxil/psb_i_sort_mod.f90 new file mode 100644 index 000000000..a7b13682a --- /dev/null +++ b/base/modules/auxil/psb_i_sort_mod.f90 @@ -0,0 +1,577 @@ +! +! Parallel Sparse BLAS version 3.5 +! (C) Copyright 2006-2018 +! Salvatore Filippone +! Alfredo Buttari +! +! 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. +! 3. The name of the PSBLAS group or the names of its contributors may +! not be used to endorse or promote products derived from this +! software without specific written permission. +! +! 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 PSBLAS GROUP OR ITS 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. +! +! +! +! Sorting routines +! References: +! D. Knuth +! The Art of Computer Programming, vol. 3 +! Addison-Wesley +! +! Aho, Hopcroft, Ullman +! Data Structures and Algorithms +! Addison-Wesley +! +module psb_i_sort_mod + use psb_const_mod + + interface psb_isaperm + logical function psb_iisaperm(n,eip) + import + integer(psb_ipk_), intent(in) :: n + integer(psb_ipk_), intent(in) :: eip(n) + end function psb_iisaperm + end interface psb_isaperm + + interface psb_msort_unique + subroutine psb_imsort_u(x,nout,dir) + import + integer(psb_ipk_), intent(inout) :: x(:) + integer(psb_ipk_), intent(out) :: nout + integer(psb_ipk_), optional, intent(in) :: dir + end subroutine psb_imsort_u + end interface psb_msort_unique + + type psb_i_heap + integer(psb_ipk_) :: last, dir + integer(psb_ipk_), allocatable :: keys(:) + contains + procedure, pass(heap) :: init => psb_i_init_heap + procedure, pass(heap) :: howmany => psb_i_howmany + procedure, pass(heap) :: insert => psb_i_insert_heap + procedure, pass(heap) :: get_first => psb_i_heap_get_first + procedure, pass(heap) :: dump => psb_i_dump_heap + procedure, pass(heap) :: free => psb_i_free_heap + end type psb_i_heap + + type psb_i_idx_heap + integer(psb_ipk_) :: last, dir + integer(psb_ipk_), allocatable :: keys(:) + integer(psb_ipk_), allocatable :: idxs(:) + contains + procedure, pass(heap) :: init => psb_i_idx_init_heap + procedure, pass(heap) :: howmany => psb_i_idx_howmany + procedure, pass(heap) :: insert => psb_i_idx_insert_heap + procedure, pass(heap) :: get_first => psb_i_idx_heap_get_first + procedure, pass(heap) :: dump => psb_i_idx_dump_heap + procedure, pass(heap) :: free => psb_i_idx_free_heap + end type psb_i_idx_heap + + + interface psb_msort + subroutine psb_imsort(x,ix,dir,flag) + import + integer(psb_ipk_), intent(inout) :: x(:) + integer(psb_ipk_), optional, intent(in) :: dir, flag + integer(psb_ipk_), optional, intent(inout) :: ix(:) + end subroutine psb_imsort + end interface psb_msort + + + interface psb_bsrch + function psb_ibsrch(key,n,v) result(ipos) + import + integer(psb_ipk_) :: ipos, n + integer(psb_ipk_) :: key + integer(psb_ipk_) :: v(:) + end function psb_ibsrch + end interface psb_bsrch + + interface psb_ssrch + function psb_issrch(key,n,v) result(ipos) + import + implicit none + integer(psb_ipk_) :: ipos, n + integer(psb_ipk_) :: key + integer(psb_ipk_) :: v(:) + end function psb_issrch + end interface psb_ssrch + + interface + subroutine psi_i_msort_up(n,k,l,iret) + import + implicit none + integer(psb_ipk_) :: n, iret + integer(psb_ipk_) :: k(n) + integer(psb_ipk_) :: l(0:n+1) + end subroutine psi_i_msort_up + subroutine psi_i_msort_dw(n,k,l,iret) + import + implicit none + integer(psb_ipk_) :: n, iret + integer(psb_ipk_) :: k(n) + integer(psb_ipk_) :: l(0:n+1) + end subroutine psi_i_msort_dw + end interface + interface + subroutine psi_i_amsort_up(n,k,l,iret) + import + implicit none + integer(psb_ipk_) :: n, iret + integer(psb_ipk_) :: k(n) + integer(psb_ipk_) :: l(0:n+1) + end subroutine psi_i_amsort_up + subroutine psi_i_amsort_dw(n,k,l,iret) + import + implicit none + integer(psb_ipk_) :: n, iret + integer(psb_ipk_) :: k(n) + integer(psb_ipk_) :: l(0:n+1) + end subroutine psi_i_amsort_dw + end interface + + + interface psb_qsort + subroutine psb_iqsort(x,ix,dir,flag) + import + integer(psb_ipk_), intent(inout) :: x(:) + integer(psb_ipk_), optional, intent(in) :: dir, flag + integer(psb_ipk_), optional, intent(inout) :: ix(:) + end subroutine psb_iqsort + end interface psb_qsort + + interface psb_isort + subroutine psb_iisort(x,ix,dir,flag) + import + integer(psb_ipk_), intent(inout) :: x(:) + integer(psb_ipk_), optional, intent(in) :: dir, flag + integer(psb_ipk_), optional, intent(inout) :: ix(:) + end subroutine psb_iisort + end interface psb_isort + + + interface psb_hsort + subroutine psb_ihsort(x,ix,dir,flag) + import + integer(psb_ipk_), intent(inout) :: x(:) + integer(psb_ipk_), optional, intent(in) :: dir, flag + integer(psb_ipk_), optional, intent(inout) :: ix(:) + end subroutine psb_ihsort + end interface psb_hsort + + + interface + subroutine psi_i_insert_heap(key,last,heap,dir,info) + import + implicit none + + ! + ! Input: + ! key: the new value + ! last: pointer to the last occupied element in heap + ! heap: the heap + ! dir: sorting direction + + integer(psb_ipk_), intent(in) :: key + integer(psb_ipk_), intent(inout) :: heap(:) + integer(psb_ipk_), intent(in) :: dir + integer(psb_ipk_), intent(inout) :: last + integer(psb_ipk_), intent(out) :: info + end subroutine psi_i_insert_heap + end interface + + interface + subroutine psi_i_idx_insert_heap(key,index,last,heap,idxs,dir,info) + import + implicit none + + ! + ! Input: + ! key: the new value + ! last: pointer to the last occupied element in heap + ! heap: the heap + ! dir: sorting direction + + integer(psb_ipk_), intent(in) :: key + integer(psb_ipk_), intent(inout) :: heap(:) + integer(psb_ipk_), intent(in) :: index + integer(psb_ipk_), intent(in) :: dir + integer(psb_ipk_), intent(inout) :: idxs(:) + integer(psb_ipk_), intent(inout) :: last + integer(psb_ipk_), intent(out) :: info + end subroutine psi_i_idx_insert_heap + end interface + + + interface + subroutine psi_i_heap_get_first(key,last,heap,dir,info) + import + implicit none + integer(psb_ipk_), intent(inout) :: key + integer(psb_ipk_), intent(inout) :: last + integer(psb_ipk_), intent(in) :: dir + integer(psb_ipk_), intent(inout) :: heap(:) + integer(psb_ipk_), intent(out) :: info + end subroutine psi_i_heap_get_first + end interface + + interface + subroutine psi_i_idx_heap_get_first(key,index,last,heap,idxs,dir,info) + import + integer(psb_ipk_), intent(inout) :: key + integer(psb_ipk_), intent(out) :: index + integer(psb_ipk_), intent(inout) :: heap(:) + integer(psb_ipk_), intent(in) :: dir + integer(psb_ipk_), intent(inout) :: last + integer(psb_ipk_), intent(inout) :: idxs(:) + integer(psb_ipk_), intent(out) :: info + end subroutine psi_i_idx_heap_get_first + end interface + + interface + subroutine psi_iisrx_up(n,x,ix) + import + integer(psb_ipk_), intent(inout) :: x(:) + integer(psb_ipk_), intent(inout) :: ix(:) + integer(psb_ipk_), intent(in) :: n + end subroutine psi_iisrx_up + subroutine psi_iisrx_dw(n,x,ix) + import + integer(psb_ipk_), intent(inout) :: x(:) + integer(psb_ipk_), intent(inout) :: ix(:) + integer(psb_ipk_), intent(in) :: n + end subroutine psi_iisrx_dw + subroutine psi_iisr_up(n,x) + import + integer(psb_ipk_), intent(inout) :: x(:) + integer(psb_ipk_), intent(in) :: n + end subroutine psi_iisr_up + subroutine psi_iisr_dw(n,x) + import + integer(psb_ipk_), intent(inout) :: x(:) + integer(psb_ipk_), intent(in) :: n + end subroutine psi_iisr_dw + subroutine psi_iaisrx_up(n,x,ix) + import + integer(psb_ipk_), intent(inout) :: x(:) + integer(psb_ipk_), intent(inout) :: ix(:) + integer(psb_ipk_), intent(in) :: n + end subroutine psi_iaisrx_up + subroutine psi_iaisrx_dw(n,x,ix) + import + integer(psb_ipk_), intent(inout) :: x(:) + integer(psb_ipk_), intent(inout) :: ix(:) + integer(psb_ipk_), intent(in) :: n + end subroutine psi_iaisrx_dw + subroutine psi_iaisr_up(n,x) + import + integer(psb_ipk_), intent(inout) :: x(:) + integer(psb_ipk_), intent(in) :: n + end subroutine psi_iaisr_up + subroutine psi_iaisr_dw(n,x) + import + integer(psb_ipk_), intent(inout) :: x(:) + integer(psb_ipk_), intent(in) :: n + end subroutine psi_iaisr_dw + end interface + + interface + subroutine psi_iqsrx_up(n,x,ix) + import + integer(psb_ipk_), intent(inout) :: x(:) + integer(psb_ipk_), intent(inout) :: ix(:) + integer(psb_ipk_), intent(in) :: n + end subroutine psi_iqsrx_up + subroutine psi_iqsrx_dw(n,x,ix) + import + integer(psb_ipk_), intent(inout) :: x(:) + integer(psb_ipk_), intent(inout) :: ix(:) + integer(psb_ipk_), intent(in) :: n + end subroutine psi_iqsrx_dw + subroutine psi_iqsr_up(n,x) + import + integer(psb_ipk_), intent(inout) :: x(:) + integer(psb_ipk_), intent(in) :: n + end subroutine psi_iqsr_up + subroutine psi_iqsr_dw(n,x) + import + integer(psb_ipk_), intent(inout) :: x(:) + integer(psb_ipk_), intent(in) :: n + end subroutine psi_iqsr_dw + subroutine psi_iaqsrx_up(n,x,ix) + import + integer(psb_ipk_), intent(inout) :: x(:) + integer(psb_ipk_), intent(inout) :: ix(:) + integer(psb_ipk_), intent(in) :: n + end subroutine psi_iaqsrx_up + subroutine psi_iaqsrx_dw(n,x,ix) + import + integer(psb_ipk_), intent(inout) :: x(:) + integer(psb_ipk_), intent(inout) :: ix(:) + integer(psb_ipk_), intent(in) :: n + end subroutine psi_iaqsrx_dw + subroutine psi_iaqsr_up(n,x) + import + integer(psb_ipk_), intent(inout) :: x(:) + integer(psb_ipk_), intent(in) :: n + end subroutine psi_iaqsr_up + subroutine psi_iaqsr_dw(n,x) + import + integer(psb_ipk_), intent(inout) :: x(:) + integer(psb_ipk_), intent(in) :: n + end subroutine psi_iaqsr_dw + end interface + +contains + + subroutine psb_i_init_heap(heap,info,dir) + use psb_realloc_mod, only : psb_ensure_size + implicit none + class(psb_i_heap), intent(inout) :: heap + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_), intent(in), optional :: dir + + info = psb_success_ + heap%last=0 + if (present(dir)) then + heap%dir = dir + else + heap%dir = psb_sort_up_ + endif + select case(heap%dir) + case (psb_sort_up_,psb_sort_down_,psb_asort_up_,psb_asort_down_) + ! ok, do nothing + case default + write(psb_err_unit,*) 'Invalid direction, defaulting to psb_sort_up_' + heap%dir = psb_sort_up_ + end select + call psb_ensure_size(psb_heap_resize,heap%keys,info) + + return + end subroutine psb_i_init_heap + + + function psb_i_howmany(heap) result(res) + implicit none + class(psb_i_heap), intent(in) :: heap + integer(psb_ipk_) :: res + res = heap%last + end function psb_i_howmany + + subroutine psb_i_insert_heap(key,heap,info) + use psb_realloc_mod, only : psb_ensure_size + implicit none + + integer(psb_ipk_), intent(in) :: key + class(psb_i_heap), intent(inout) :: heap + integer(psb_ipk_), intent(out) :: info + + info = psb_success_ + if (heap%last < 0) then + write(psb_err_unit,*) 'Invalid last in heap ',heap%last + info = heap%last + return + endif + + call psb_ensure_size(heap%last+1,heap%keys,info,addsz=psb_heap_resize) + if (info /= psb_success_) then + write(psb_err_unit,*) 'Memory allocation failure in heap_insert' + info = -5 + return + end if + call psi_i_insert_heap(key,& + & heap%last,heap%keys,heap%dir,info) + + return + end subroutine psb_i_insert_heap + + subroutine psb_i_heap_get_first(key,heap,info) + implicit none + + class(psb_i_heap), intent(inout) :: heap + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_), intent(out) :: key + + + info = psb_success_ + + call psi_i_heap_get_first(key,& + & heap%last,heap%keys,heap%dir,info) + + return + end subroutine psb_i_heap_get_first + + subroutine psb_i_dump_heap(iout,heap,info) + + implicit none + class(psb_i_heap), intent(in) :: heap + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_), intent(in) :: iout + + info = psb_success_ + if (iout < 0) then + write(psb_err_unit,*) 'Invalid file ' + info =-1 + return + end if + + write(iout,*) 'Heap direction ',heap%dir + write(iout,*) 'Heap size ',heap%last + if ((heap%last > 0).and.((.not.allocated(heap%keys)).or.& + & (size(heap%keys) 0).and.((.not.allocated(heap%keys)).or.& + & (size(heap%keys) 0).and.((.not.allocated(heap%idxs)).or.& + & (size(heap%idxs) psb_l_init_heap + procedure, pass(heap) :: howmany => psb_l_howmany + procedure, pass(heap) :: insert => psb_l_insert_heap + procedure, pass(heap) :: get_first => psb_l_heap_get_first + procedure, pass(heap) :: dump => psb_l_dump_heap + procedure, pass(heap) :: free => psb_l_free_heap + end type psb_l_heap + + type psb_l_idx_heap + integer(psb_ipk_) :: last, dir + integer(psb_lpk_), allocatable :: keys(:) + integer(psb_lpk_), allocatable :: idxs(:) + contains + procedure, pass(heap) :: init => psb_l_idx_init_heap + procedure, pass(heap) :: howmany => psb_l_idx_howmany + procedure, pass(heap) :: insert => psb_l_idx_insert_heap + procedure, pass(heap) :: get_first => psb_l_idx_heap_get_first + procedure, pass(heap) :: dump => psb_l_idx_dump_heap + procedure, pass(heap) :: free => psb_l_idx_free_heap + end type psb_l_idx_heap + + + interface psb_msort + subroutine psb_lmsort(x,ix,dir,flag) + import + integer(psb_lpk_), intent(inout) :: x(:) + integer(psb_ipk_), optional, intent(in) :: dir, flag + integer(psb_lpk_), optional, intent(inout) :: ix(:) + end subroutine psb_lmsort + end interface psb_msort + + + interface psb_bsrch + function psb_lbsrch(key,n,v) result(ipos) + import + integer(psb_ipk_) :: ipos, n + integer(psb_lpk_) :: key + integer(psb_lpk_) :: v(:) + end function psb_lbsrch + end interface psb_bsrch + + interface psb_ssrch + function psb_lssrch(key,n,v) result(ipos) + import + implicit none + integer(psb_ipk_) :: ipos, n + integer(psb_lpk_) :: key + integer(psb_lpk_) :: v(:) + end function psb_lssrch + end interface psb_ssrch + + interface + subroutine psi_l_msort_up(n,k,l,iret) + import + implicit none + integer(psb_ipk_) :: n, iret + integer(psb_lpk_) :: k(n) + integer(psb_lpk_) :: l(0:n+1) + end subroutine psi_l_msort_up + subroutine psi_l_msort_dw(n,k,l,iret) + import + implicit none + integer(psb_ipk_) :: n, iret + integer(psb_lpk_) :: k(n) + integer(psb_lpk_) :: l(0:n+1) + end subroutine psi_l_msort_dw + end interface + interface + subroutine psi_l_amsort_up(n,k,l,iret) + import + implicit none + integer(psb_ipk_) :: n, iret + integer(psb_lpk_) :: k(n) + integer(psb_lpk_) :: l(0:n+1) + end subroutine psi_l_amsort_up + subroutine psi_l_amsort_dw(n,k,l,iret) + import + implicit none + integer(psb_ipk_) :: n, iret + integer(psb_lpk_) :: k(n) + integer(psb_lpk_) :: l(0:n+1) + end subroutine psi_l_amsort_dw + end interface + + + interface psb_qsort + subroutine psb_lqsort(x,ix,dir,flag) + import + integer(psb_lpk_), intent(inout) :: x(:) + integer(psb_ipk_), optional, intent(in) :: dir, flag + integer(psb_lpk_), optional, intent(inout) :: ix(:) + end subroutine psb_lqsort + end interface psb_qsort + + interface psb_isort + subroutine psb_lisort(x,ix,dir,flag) + import + integer(psb_lpk_), intent(inout) :: x(:) + integer(psb_ipk_), optional, intent(in) :: dir, flag + integer(psb_lpk_), optional, intent(inout) :: ix(:) + end subroutine psb_lisort + end interface psb_isort + + + interface psb_hsort + subroutine psb_lhsort(x,ix,dir,flag) + import + integer(psb_lpk_), intent(inout) :: x(:) + integer(psb_ipk_), optional, intent(in) :: dir, flag + integer(psb_lpk_), optional, intent(inout) :: ix(:) + end subroutine psb_lhsort + end interface psb_hsort + + + interface + subroutine psi_l_insert_heap(key,last,heap,dir,info) + import + implicit none + + ! + ! Input: + ! key: the new value + ! last: pointer to the last occupied element in heap + ! heap: the heap + ! dir: sorting direction + + integer(psb_lpk_), intent(in) :: key + integer(psb_lpk_), intent(inout) :: heap(:) + integer(psb_ipk_), intent(in) :: dir + integer(psb_ipk_), intent(inout) :: last + integer(psb_ipk_), intent(out) :: info + end subroutine psi_l_insert_heap + end interface + + interface + subroutine psi_l_idx_insert_heap(key,index,last,heap,idxs,dir,info) + import + implicit none + + ! + ! Input: + ! key: the new value + ! last: pointer to the last occupied element in heap + ! heap: the heap + ! dir: sorting direction + + integer(psb_lpk_), intent(in) :: key + integer(psb_lpk_), intent(inout) :: heap(:) + integer(psb_lpk_), intent(in) :: index + integer(psb_ipk_), intent(in) :: dir + integer(psb_lpk_), intent(inout) :: idxs(:) + integer(psb_ipk_), intent(inout) :: last + integer(psb_ipk_), intent(out) :: info + end subroutine psi_l_idx_insert_heap + end interface + + + interface + subroutine psi_l_heap_get_first(key,last,heap,dir,info) + import + implicit none + integer(psb_lpk_), intent(inout) :: key + integer(psb_ipk_), intent(inout) :: last + integer(psb_ipk_), intent(in) :: dir + integer(psb_lpk_), intent(inout) :: heap(:) + integer(psb_ipk_), intent(out) :: info + end subroutine psi_l_heap_get_first + end interface + + interface + subroutine psi_l_idx_heap_get_first(key,index,last,heap,idxs,dir,info) + import + integer(psb_lpk_), intent(inout) :: key + integer(psb_lpk_), intent(out) :: index + integer(psb_lpk_), intent(inout) :: heap(:) + integer(psb_ipk_), intent(in) :: dir + integer(psb_ipk_), intent(inout) :: last + integer(psb_lpk_), intent(inout) :: idxs(:) + integer(psb_ipk_), intent(out) :: info + end subroutine psi_l_idx_heap_get_first + end interface + + interface + subroutine psi_lisrx_up(n,x,ix) + import + integer(psb_lpk_), intent(inout) :: x(:) + integer(psb_lpk_), intent(inout) :: ix(:) + integer(psb_lpk_), intent(in) :: n + end subroutine psi_lisrx_up + subroutine psi_lisrx_dw(n,x,ix) + import + integer(psb_lpk_), intent(inout) :: x(:) + integer(psb_lpk_), intent(inout) :: ix(:) + integer(psb_lpk_), intent(in) :: n + end subroutine psi_lisrx_dw + subroutine psi_lisr_up(n,x) + import + integer(psb_lpk_), intent(inout) :: x(:) + integer(psb_lpk_), intent(in) :: n + end subroutine psi_lisr_up + subroutine psi_lisr_dw(n,x) + import + integer(psb_lpk_), intent(inout) :: x(:) + integer(psb_lpk_), intent(in) :: n + end subroutine psi_lisr_dw + subroutine psi_laisrx_up(n,x,ix) + import + integer(psb_lpk_), intent(inout) :: x(:) + integer(psb_lpk_), intent(inout) :: ix(:) + integer(psb_lpk_), intent(in) :: n + end subroutine psi_laisrx_up + subroutine psi_laisrx_dw(n,x,ix) + import + integer(psb_lpk_), intent(inout) :: x(:) + integer(psb_lpk_), intent(inout) :: ix(:) + integer(psb_lpk_), intent(in) :: n + end subroutine psi_laisrx_dw + subroutine psi_laisr_up(n,x) + import + integer(psb_lpk_), intent(inout) :: x(:) + integer(psb_lpk_), intent(in) :: n + end subroutine psi_laisr_up + subroutine psi_laisr_dw(n,x) + import + integer(psb_lpk_), intent(inout) :: x(:) + integer(psb_lpk_), intent(in) :: n + end subroutine psi_laisr_dw + end interface + + interface + subroutine psi_lqsrx_up(n,x,ix) + import + integer(psb_lpk_), intent(inout) :: x(:) + integer(psb_lpk_), intent(inout) :: ix(:) + integer(psb_lpk_), intent(in) :: n + end subroutine psi_lqsrx_up + subroutine psi_lqsrx_dw(n,x,ix) + import + integer(psb_lpk_), intent(inout) :: x(:) + integer(psb_lpk_), intent(inout) :: ix(:) + integer(psb_lpk_), intent(in) :: n + end subroutine psi_lqsrx_dw + subroutine psi_lqsr_up(n,x) + import + integer(psb_lpk_), intent(inout) :: x(:) + integer(psb_lpk_), intent(in) :: n + end subroutine psi_lqsr_up + subroutine psi_lqsr_dw(n,x) + import + integer(psb_lpk_), intent(inout) :: x(:) + integer(psb_lpk_), intent(in) :: n + end subroutine psi_lqsr_dw + subroutine psi_laqsrx_up(n,x,ix) + import + integer(psb_lpk_), intent(inout) :: x(:) + integer(psb_lpk_), intent(inout) :: ix(:) + integer(psb_lpk_), intent(in) :: n + end subroutine psi_laqsrx_up + subroutine psi_laqsrx_dw(n,x,ix) + import + integer(psb_lpk_), intent(inout) :: x(:) + integer(psb_lpk_), intent(inout) :: ix(:) + integer(psb_lpk_), intent(in) :: n + end subroutine psi_laqsrx_dw + subroutine psi_laqsr_up(n,x) + import + integer(psb_lpk_), intent(inout) :: x(:) + integer(psb_lpk_), intent(in) :: n + end subroutine psi_laqsr_up + subroutine psi_laqsr_dw(n,x) + import + integer(psb_lpk_), intent(inout) :: x(:) + integer(psb_lpk_), intent(in) :: n + end subroutine psi_laqsr_dw + end interface + +contains + + subroutine psb_l_init_heap(heap,info,dir) + use psb_realloc_mod, only : psb_ensure_size + implicit none + class(psb_l_heap), intent(inout) :: heap + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_), intent(in), optional :: dir + + info = psb_success_ + heap%last=0 + if (present(dir)) then + heap%dir = dir + else + heap%dir = psb_sort_up_ + endif + select case(heap%dir) + case (psb_sort_up_,psb_sort_down_,psb_asort_up_,psb_asort_down_) + ! ok, do nothing + case default + write(psb_err_unit,*) 'Invalid direction, defaulting to psb_sort_up_' + heap%dir = psb_sort_up_ + end select + call psb_ensure_size(psb_heap_resize,heap%keys,info) + + return + end subroutine psb_l_init_heap + + + function psb_l_howmany(heap) result(res) + implicit none + class(psb_l_heap), intent(in) :: heap + integer(psb_ipk_) :: res + res = heap%last + end function psb_l_howmany + + subroutine psb_l_insert_heap(key,heap,info) + use psb_realloc_mod, only : psb_ensure_size + implicit none + + integer(psb_lpk_), intent(in) :: key + class(psb_l_heap), intent(inout) :: heap + integer(psb_ipk_), intent(out) :: info + + info = psb_success_ + if (heap%last < 0) then + write(psb_err_unit,*) 'Invalid last in heap ',heap%last + info = heap%last + return + endif + + call psb_ensure_size(heap%last+1,heap%keys,info,addsz=psb_heap_resize) + if (info /= psb_success_) then + write(psb_err_unit,*) 'Memory allocation failure in heap_insert' + info = -5 + return + end if + call psi_l_insert_heap(key,& + & heap%last,heap%keys,heap%dir,info) + + return + end subroutine psb_l_insert_heap + + subroutine psb_l_heap_get_first(key,heap,info) + implicit none + + class(psb_l_heap), intent(inout) :: heap + integer(psb_ipk_), intent(out) :: info + integer(psb_lpk_), intent(out) :: key + + + info = psb_success_ + + call psi_l_heap_get_first(key,& + & heap%last,heap%keys,heap%dir,info) + + return + end subroutine psb_l_heap_get_first + + subroutine psb_l_dump_heap(iout,heap,info) + + implicit none + class(psb_l_heap), intent(in) :: heap + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_), intent(in) :: iout + + info = psb_success_ + if (iout < 0) then + write(psb_err_unit,*) 'Invalid file ' + info =-1 + return + end if + + write(iout,*) 'Heap direction ',heap%dir + write(iout,*) 'Heap size ',heap%last + if ((heap%last > 0).and.((.not.allocated(heap%keys)).or.& + & (size(heap%keys) 0).and.((.not.allocated(heap%keys)).or.& + & (size(heap%keys) 0).and.((.not.allocated(heap%idxs)).or.& + & (size(heap%idxs)n)) exit + do + if (lp >= k) exit + lp = iaux(lp) + end do + swap = x(lp) + x(lp) = x(k) + x(k) = swap + lswap = iaux(lp) + iaux(lp) = iaux(k) + iaux(k) = lp + lp = lswap + k = k + 1 + enddo + return + end subroutine psb_ip_reord_m1 + + subroutine psb_ip_reord_m1i1(n,x,indx,iaux) + integer(psb_ipk_), intent(in) :: n + integer(psb_mpk_) :: iaux(0:*) + integer(psb_mpk_) :: x(*) + integer(psb_mpk_) :: indx(*) + integer(psb_mpk_) :: lswap, lp, k, ixswap + integer(psb_mpk_) :: swap + + + + lp = iaux(0) + k = 1 + do + if ((lp == 0).or.(k>n)) exit + do + if (lp >= k) exit + lp = iaux(lp) + end do + swap = x(lp) + x(lp) = x(k) + x(k) = swap + ixswap = indx(lp) + indx(lp) = indx(k) + indx(k) = ixswap + lswap = iaux(lp) + iaux(lp) = iaux(k) + iaux(k) = lp + lp = lswap + k = k + 1 + enddo + return + end subroutine psb_ip_reord_m1i1 + + subroutine psb_ip_reord_m1i2(n,x,i1,i2,iaux) + integer(psb_ipk_), intent(in) :: n + integer(psb_ipk_) :: iaux(0:*) + integer(psb_mpk_) :: x(*) + integer(psb_ipk_) :: i1(*), i2(*) + + + integer(psb_ipk_) :: lswap, lp, k, isw1, isw2 + integer(psb_mpk_) :: swap + + lp = iaux(0) + k = 1 + do + if ((lp == 0).or.(k>n)) exit + do + if (lp >= k) exit + lp = iaux(lp) + end do + swap = x(lp) + x(lp) = x(k) + x(k) = swap + isw1 = i1(lp) + i1(lp) = i1(k) + i1(k) = isw1 + isw2 = i2(lp) + i2(lp) = i2(k) + i2(k) = isw2 + lswap = iaux(lp) + iaux(lp) = iaux(k) + iaux(k) = lp + lp = lswap + k = k + 1 + enddo + return + end subroutine psb_ip_reord_m1i2 + + subroutine psb_ip_reord_m1i3(n,x,i1,i2,i3,iaux) + integer(psb_ipk_), intent(in) :: n + integer(psb_ipk_) :: iaux(0:*) + integer(psb_mpk_) :: x(*) + integer(psb_ipk_) :: i1(*), i2(*), i3(*) + + integer(psb_ipk_) :: lswap, lp, k, isw1, isw2, isw3 + integer(psb_mpk_) :: swap + + lp = iaux(0) + k = 1 + do + if ((lp == 0).or.(k>n)) exit + do + if (lp >= k) exit + lp = iaux(lp) + end do + swap = x(lp) + x(lp) = x(k) + x(k) = swap + isw1 = i1(lp) + i1(lp) = i1(k) + i1(k) = isw1 + isw2 = i2(lp) + i2(lp) = i2(k) + i2(k) = isw2 + isw3 = i3(lp) + i3(lp) = i3(k) + i3(k) = isw3 + lswap = iaux(lp) + iaux(lp) = iaux(k) + iaux(k) = lp + lp = lswap + k = k + 1 + enddo + return + end subroutine psb_ip_reord_m1i3 + +end module psb_m_ip_reord_mod diff --git a/base/modules/auxil/psb_s_ip_reord_mod.f90 b/base/modules/auxil/psb_s_ip_reord_mod.f90 index 20c5f3601..9ab4c5c4f 100644 --- a/base/modules/auxil/psb_s_ip_reord_mod.f90 +++ b/base/modules/auxil/psb_s_ip_reord_mod.f90 @@ -53,10 +53,9 @@ contains integer(psb_ipk_), intent(in) :: n integer(psb_ipk_) :: iaux(0:*) real(psb_spk_) :: x(*) - integer(psb_ipk_) :: lswap, lp, k real(psb_spk_) :: swap - + lp = iaux(0) k = 1 do @@ -82,10 +81,11 @@ contains integer(psb_ipk_) :: iaux(0:*) real(psb_spk_) :: x(*) integer(psb_ipk_) :: indx(*) - integer(psb_ipk_) :: lswap, lp, k, ixswap real(psb_spk_) :: swap + + lp = iaux(0) k = 1 do @@ -150,7 +150,6 @@ contains integer(psb_ipk_) :: iaux(0:*) real(psb_spk_) :: x(*) integer(psb_ipk_) :: i1(*), i2(*), i3(*) - integer(psb_ipk_) :: lswap, lp, k, isw1, isw2, isw3 real(psb_spk_) :: swap diff --git a/base/modules/auxil/psb_s_sort_mod.f90 b/base/modules/auxil/psb_s_sort_mod.f90 new file mode 100644 index 000000000..bd67ac068 --- /dev/null +++ b/base/modules/auxil/psb_s_sort_mod.f90 @@ -0,0 +1,570 @@ +! +! Parallel Sparse BLAS version 3.5 +! (C) Copyright 2006-2018 +! Salvatore Filippone +! Alfredo Buttari +! +! 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. +! 3. The name of the PSBLAS group or the names of its contributors may +! not be used to endorse or promote products derived from this +! software without specific written permission. +! +! 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 PSBLAS GROUP OR ITS 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. +! +! +! +! Sorting routines +! References: +! D. Knuth +! The Art of Computer Programming, vol. 3 +! Addison-Wesley +! +! Aho, Hopcroft, Ullman +! Data Structures and Algorithms +! Addison-Wesley +! +module psb_s_sort_mod + use psb_const_mod + + + interface psb_msort_unique + subroutine psb_smsort_u(x,nout,dir) + import + real(psb_spk_), intent(inout) :: x(:) + integer(psb_ipk_), intent(out) :: nout + integer(psb_ipk_), optional, intent(in) :: dir + end subroutine psb_smsort_u + end interface psb_msort_unique + + type psb_s_heap + integer(psb_ipk_) :: last, dir + real(psb_spk_), allocatable :: keys(:) + contains + procedure, pass(heap) :: init => psb_s_init_heap + procedure, pass(heap) :: howmany => psb_s_howmany + procedure, pass(heap) :: insert => psb_s_insert_heap + procedure, pass(heap) :: get_first => psb_s_heap_get_first + procedure, pass(heap) :: dump => psb_s_dump_heap + procedure, pass(heap) :: free => psb_s_free_heap + end type psb_s_heap + + type psb_s_idx_heap + integer(psb_ipk_) :: last, dir + real(psb_spk_), allocatable :: keys(:) + integer(psb_ipk_), allocatable :: idxs(:) + contains + procedure, pass(heap) :: init => psb_s_idx_init_heap + procedure, pass(heap) :: howmany => psb_s_idx_howmany + procedure, pass(heap) :: insert => psb_s_idx_insert_heap + procedure, pass(heap) :: get_first => psb_s_idx_heap_get_first + procedure, pass(heap) :: dump => psb_s_idx_dump_heap + procedure, pass(heap) :: free => psb_s_idx_free_heap + end type psb_s_idx_heap + + + interface psb_msort + subroutine psb_smsort(x,ix,dir,flag) + import + real(psb_spk_), intent(inout) :: x(:) + integer(psb_ipk_), optional, intent(in) :: dir, flag + integer(psb_ipk_), optional, intent(inout) :: ix(:) + end subroutine psb_smsort + end interface psb_msort + + + interface psb_bsrch + function psb_sbsrch(key,n,v) result(ipos) + import + integer(psb_ipk_) :: ipos, n + real(psb_spk_) :: key + real(psb_spk_) :: v(:) + end function psb_sbsrch + end interface psb_bsrch + + interface psb_ssrch + function psb_sssrch(key,n,v) result(ipos) + import + implicit none + integer(psb_ipk_) :: ipos, n + real(psb_spk_) :: key + real(psb_spk_) :: v(:) + end function psb_sssrch + end interface psb_ssrch + + interface + subroutine psi_s_msort_up(n,k,l,iret) + import + implicit none + integer(psb_ipk_) :: n, iret + real(psb_spk_) :: k(n) + integer(psb_ipk_) :: l(0:n+1) + end subroutine psi_s_msort_up + subroutine psi_s_msort_dw(n,k,l,iret) + import + implicit none + integer(psb_ipk_) :: n, iret + real(psb_spk_) :: k(n) + integer(psb_ipk_) :: l(0:n+1) + end subroutine psi_s_msort_dw + end interface + interface + subroutine psi_s_amsort_up(n,k,l,iret) + import + implicit none + integer(psb_ipk_) :: n, iret + real(psb_spk_) :: k(n) + integer(psb_ipk_) :: l(0:n+1) + end subroutine psi_s_amsort_up + subroutine psi_s_amsort_dw(n,k,l,iret) + import + implicit none + integer(psb_ipk_) :: n, iret + real(psb_spk_) :: k(n) + integer(psb_ipk_) :: l(0:n+1) + end subroutine psi_s_amsort_dw + end interface + + + interface psb_qsort + subroutine psb_sqsort(x,ix,dir,flag) + import + real(psb_spk_), intent(inout) :: x(:) + integer(psb_ipk_), optional, intent(in) :: dir, flag + integer(psb_ipk_), optional, intent(inout) :: ix(:) + end subroutine psb_sqsort + end interface psb_qsort + + interface psb_isort + subroutine psb_sisort(x,ix,dir,flag) + import + real(psb_spk_), intent(inout) :: x(:) + integer(psb_ipk_), optional, intent(in) :: dir, flag + integer(psb_ipk_), optional, intent(inout) :: ix(:) + end subroutine psb_sisort + end interface psb_isort + + + interface psb_hsort + subroutine psb_shsort(x,ix,dir,flag) + import + real(psb_spk_), intent(inout) :: x(:) + integer(psb_ipk_), optional, intent(in) :: dir, flag + integer(psb_ipk_), optional, intent(inout) :: ix(:) + end subroutine psb_shsort + end interface psb_hsort + + + interface + subroutine psi_s_insert_heap(key,last,heap,dir,info) + import + implicit none + + ! + ! Input: + ! key: the new value + ! last: pointer to the last occupied element in heap + ! heap: the heap + ! dir: sorting direction + + real(psb_spk_), intent(in) :: key + real(psb_spk_), intent(inout) :: heap(:) + integer(psb_ipk_), intent(in) :: dir + integer(psb_ipk_), intent(inout) :: last + integer(psb_ipk_), intent(out) :: info + end subroutine psi_s_insert_heap + end interface + + interface + subroutine psi_s_idx_insert_heap(key,index,last,heap,idxs,dir,info) + import + implicit none + + ! + ! Input: + ! key: the new value + ! last: pointer to the last occupied element in heap + ! heap: the heap + ! dir: sorting direction + + real(psb_spk_), intent(in) :: key + real(psb_spk_), intent(inout) :: heap(:) + integer(psb_ipk_), intent(in) :: index + integer(psb_ipk_), intent(in) :: dir + integer(psb_ipk_), intent(inout) :: idxs(:) + integer(psb_ipk_), intent(inout) :: last + integer(psb_ipk_), intent(out) :: info + end subroutine psi_s_idx_insert_heap + end interface + + + interface + subroutine psi_s_heap_get_first(key,last,heap,dir,info) + import + implicit none + real(psb_spk_), intent(inout) :: key + integer(psb_ipk_), intent(inout) :: last + integer(psb_ipk_), intent(in) :: dir + real(psb_spk_), intent(inout) :: heap(:) + integer(psb_ipk_), intent(out) :: info + end subroutine psi_s_heap_get_first + end interface + + interface + subroutine psi_s_idx_heap_get_first(key,index,last,heap,idxs,dir,info) + import + real(psb_spk_), intent(inout) :: key + integer(psb_ipk_), intent(out) :: index + real(psb_spk_), intent(inout) :: heap(:) + integer(psb_ipk_), intent(in) :: dir + integer(psb_ipk_), intent(inout) :: last + integer(psb_ipk_), intent(inout) :: idxs(:) + integer(psb_ipk_), intent(out) :: info + end subroutine psi_s_idx_heap_get_first + end interface + + interface + subroutine psi_sisrx_up(n,x,ix) + import + real(psb_spk_), intent(inout) :: x(:) + integer(psb_ipk_), intent(inout) :: ix(:) + integer(psb_ipk_), intent(in) :: n + end subroutine psi_sisrx_up + subroutine psi_sisrx_dw(n,x,ix) + import + real(psb_spk_), intent(inout) :: x(:) + integer(psb_ipk_), intent(inout) :: ix(:) + integer(psb_ipk_), intent(in) :: n + end subroutine psi_sisrx_dw + subroutine psi_sisr_up(n,x) + import + real(psb_spk_), intent(inout) :: x(:) + integer(psb_ipk_), intent(in) :: n + end subroutine psi_sisr_up + subroutine psi_sisr_dw(n,x) + import + real(psb_spk_), intent(inout) :: x(:) + integer(psb_ipk_), intent(in) :: n + end subroutine psi_sisr_dw + subroutine psi_saisrx_up(n,x,ix) + import + real(psb_spk_), intent(inout) :: x(:) + integer(psb_ipk_), intent(inout) :: ix(:) + integer(psb_ipk_), intent(in) :: n + end subroutine psi_saisrx_up + subroutine psi_saisrx_dw(n,x,ix) + import + real(psb_spk_), intent(inout) :: x(:) + integer(psb_ipk_), intent(inout) :: ix(:) + integer(psb_ipk_), intent(in) :: n + end subroutine psi_saisrx_dw + subroutine psi_saisr_up(n,x) + import + real(psb_spk_), intent(inout) :: x(:) + integer(psb_ipk_), intent(in) :: n + end subroutine psi_saisr_up + subroutine psi_saisr_dw(n,x) + import + real(psb_spk_), intent(inout) :: x(:) + integer(psb_ipk_), intent(in) :: n + end subroutine psi_saisr_dw + end interface + + interface + subroutine psi_sqsrx_up(n,x,ix) + import + real(psb_spk_), intent(inout) :: x(:) + integer(psb_ipk_), intent(inout) :: ix(:) + integer(psb_ipk_), intent(in) :: n + end subroutine psi_sqsrx_up + subroutine psi_sqsrx_dw(n,x,ix) + import + real(psb_spk_), intent(inout) :: x(:) + integer(psb_ipk_), intent(inout) :: ix(:) + integer(psb_ipk_), intent(in) :: n + end subroutine psi_sqsrx_dw + subroutine psi_sqsr_up(n,x) + import + real(psb_spk_), intent(inout) :: x(:) + integer(psb_ipk_), intent(in) :: n + end subroutine psi_sqsr_up + subroutine psi_sqsr_dw(n,x) + import + real(psb_spk_), intent(inout) :: x(:) + integer(psb_ipk_), intent(in) :: n + end subroutine psi_sqsr_dw + subroutine psi_saqsrx_up(n,x,ix) + import + real(psb_spk_), intent(inout) :: x(:) + integer(psb_ipk_), intent(inout) :: ix(:) + integer(psb_ipk_), intent(in) :: n + end subroutine psi_saqsrx_up + subroutine psi_saqsrx_dw(n,x,ix) + import + real(psb_spk_), intent(inout) :: x(:) + integer(psb_ipk_), intent(inout) :: ix(:) + integer(psb_ipk_), intent(in) :: n + end subroutine psi_saqsrx_dw + subroutine psi_saqsr_up(n,x) + import + real(psb_spk_), intent(inout) :: x(:) + integer(psb_ipk_), intent(in) :: n + end subroutine psi_saqsr_up + subroutine psi_saqsr_dw(n,x) + import + real(psb_spk_), intent(inout) :: x(:) + integer(psb_ipk_), intent(in) :: n + end subroutine psi_saqsr_dw + end interface + +contains + + subroutine psb_s_init_heap(heap,info,dir) + use psb_realloc_mod, only : psb_ensure_size + implicit none + class(psb_s_heap), intent(inout) :: heap + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_), intent(in), optional :: dir + + info = psb_success_ + heap%last=0 + if (present(dir)) then + heap%dir = dir + else + heap%dir = psb_sort_up_ + endif + select case(heap%dir) + case (psb_sort_up_,psb_sort_down_,psb_asort_up_,psb_asort_down_) + ! ok, do nothing + case default + write(psb_err_unit,*) 'Invalid direction, defaulting to psb_sort_up_' + heap%dir = psb_sort_up_ + end select + call psb_ensure_size(psb_heap_resize,heap%keys,info) + + return + end subroutine psb_s_init_heap + + + function psb_s_howmany(heap) result(res) + implicit none + class(psb_s_heap), intent(in) :: heap + integer(psb_ipk_) :: res + res = heap%last + end function psb_s_howmany + + subroutine psb_s_insert_heap(key,heap,info) + use psb_realloc_mod, only : psb_ensure_size + implicit none + + real(psb_spk_), intent(in) :: key + class(psb_s_heap), intent(inout) :: heap + integer(psb_ipk_), intent(out) :: info + + info = psb_success_ + if (heap%last < 0) then + write(psb_err_unit,*) 'Invalid last in heap ',heap%last + info = heap%last + return + endif + + call psb_ensure_size(heap%last+1,heap%keys,info,addsz=psb_heap_resize) + if (info /= psb_success_) then + write(psb_err_unit,*) 'Memory allocation failure in heap_insert' + info = -5 + return + end if + call psi_s_insert_heap(key,& + & heap%last,heap%keys,heap%dir,info) + + return + end subroutine psb_s_insert_heap + + subroutine psb_s_heap_get_first(key,heap,info) + implicit none + + class(psb_s_heap), intent(inout) :: heap + integer(psb_ipk_), intent(out) :: info + real(psb_spk_), intent(out) :: key + + + info = psb_success_ + + call psi_s_heap_get_first(key,& + & heap%last,heap%keys,heap%dir,info) + + return + end subroutine psb_s_heap_get_first + + subroutine psb_s_dump_heap(iout,heap,info) + + implicit none + class(psb_s_heap), intent(in) :: heap + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_), intent(in) :: iout + + info = psb_success_ + if (iout < 0) then + write(psb_err_unit,*) 'Invalid file ' + info =-1 + return + end if + + write(iout,*) 'Heap direction ',heap%dir + write(iout,*) 'Heap size ',heap%last + if ((heap%last > 0).and.((.not.allocated(heap%keys)).or.& + & (size(heap%keys) 0).and.((.not.allocated(heap%keys)).or.& + & (size(heap%keys) 0).and.((.not.allocated(heap%idxs)).or.& + & (size(heap%idxs) psb_z_init_heap + procedure, pass(heap) :: howmany => psb_z_howmany + procedure, pass(heap) :: insert => psb_z_insert_heap + procedure, pass(heap) :: get_first => psb_z_heap_get_first + procedure, pass(heap) :: dump => psb_z_dump_heap + procedure, pass(heap) :: free => psb_z_free_heap + end type psb_z_heap + + type psb_z_idx_heap + integer(psb_ipk_) :: last, dir + complex(psb_dpk_), allocatable :: keys(:) + integer(psb_ipk_), allocatable :: idxs(:) + contains + procedure, pass(heap) :: init => psb_z_idx_init_heap + procedure, pass(heap) :: howmany => psb_z_idx_howmany + procedure, pass(heap) :: insert => psb_z_idx_insert_heap + procedure, pass(heap) :: get_first => psb_z_idx_heap_get_first + procedure, pass(heap) :: dump => psb_z_idx_dump_heap + procedure, pass(heap) :: free => psb_z_idx_free_heap + end type psb_z_idx_heap + + + interface psb_msort + subroutine psb_zmsort(x,ix,dir,flag) + import + complex(psb_dpk_), intent(inout) :: x(:) + integer(psb_ipk_), optional, intent(in) :: dir, flag + integer(psb_ipk_), optional, intent(inout) :: ix(:) + end subroutine psb_zmsort + end interface psb_msort + + interface + subroutine psi_z_lmsort_up(n,k,l,iret) + import + implicit none + integer(psb_ipk_) :: n, iret + complex(psb_dpk_) :: k(n) + integer(psb_ipk_) :: l(0:n+1) + end subroutine psi_z_lmsort_up + subroutine psi_z_lmsort_dw(n,k,l,iret) + import + implicit none + integer(psb_ipk_) :: n, iret + complex(psb_dpk_) :: k(n) + integer(psb_ipk_) :: l(0:n+1) + end subroutine psi_z_lmsort_dw + subroutine psi_z_almsort_up(n,k,l,iret) + import + implicit none + integer(psb_ipk_) :: n, iret + complex(psb_dpk_) :: k(n) + integer(psb_ipk_) :: l(0:n+1) + end subroutine psi_z_almsort_up + subroutine psi_z_almsort_dw(n,k,l,iret) + import + implicit none + integer(psb_ipk_) :: n, iret + complex(psb_dpk_) :: k(n) + integer(psb_ipk_) :: l(0:n+1) + end subroutine psi_z_almsort_dw + end interface + interface + subroutine psi_z_amsort_up(n,k,l,iret) + import + implicit none + integer(psb_ipk_) :: n, iret + complex(psb_dpk_) :: k(n) + integer(psb_ipk_) :: l(0:n+1) + end subroutine psi_z_amsort_up + subroutine psi_z_amsort_dw(n,k,l,iret) + import + implicit none + integer(psb_ipk_) :: n, iret + complex(psb_dpk_) :: k(n) + integer(psb_ipk_) :: l(0:n+1) + end subroutine psi_z_amsort_dw + end interface + + + interface psb_qsort + subroutine psb_zqsort(x,ix,dir,flag) + import + complex(psb_dpk_), intent(inout) :: x(:) + integer(psb_ipk_), optional, intent(in) :: dir, flag + integer(psb_ipk_), optional, intent(inout) :: ix(:) + end subroutine psb_zqsort + end interface psb_qsort + + interface psb_isort + subroutine psb_zisort(x,ix,dir,flag) + import + complex(psb_dpk_), intent(inout) :: x(:) + integer(psb_ipk_), optional, intent(in) :: dir, flag + integer(psb_ipk_), optional, intent(inout) :: ix(:) + end subroutine psb_zisort + end interface psb_isort + + + interface psb_hsort + subroutine psb_zhsort(x,ix,dir,flag) + import + complex(psb_dpk_), intent(inout) :: x(:) + integer(psb_ipk_), optional, intent(in) :: dir, flag + integer(psb_ipk_), optional, intent(inout) :: ix(:) + end subroutine psb_zhsort + end interface psb_hsort + + + interface + subroutine psi_z_insert_heap(key,last,heap,dir,info) + import + implicit none + + ! + ! Input: + ! key: the new value + ! last: pointer to the last occupied element in heap + ! heap: the heap + ! dir: sorting direction + + complex(psb_dpk_), intent(in) :: key + complex(psb_dpk_), intent(inout) :: heap(:) + integer(psb_ipk_), intent(in) :: dir + integer(psb_ipk_), intent(inout) :: last + integer(psb_ipk_), intent(out) :: info + end subroutine psi_z_insert_heap + end interface + + interface + subroutine psi_z_idx_insert_heap(key,index,last,heap,idxs,dir,info) + import + implicit none + + ! + ! Input: + ! key: the new value + ! last: pointer to the last occupied element in heap + ! heap: the heap + ! dir: sorting direction + + complex(psb_dpk_), intent(in) :: key + complex(psb_dpk_), intent(inout) :: heap(:) + integer(psb_ipk_), intent(in) :: index + integer(psb_ipk_), intent(in) :: dir + integer(psb_ipk_), intent(inout) :: idxs(:) + integer(psb_ipk_), intent(inout) :: last + integer(psb_ipk_), intent(out) :: info + end subroutine psi_z_idx_insert_heap + end interface + + + interface + subroutine psi_z_heap_get_first(key,last,heap,dir,info) + import + implicit none + complex(psb_dpk_), intent(inout) :: key + integer(psb_ipk_), intent(inout) :: last + integer(psb_ipk_), intent(in) :: dir + complex(psb_dpk_), intent(inout) :: heap(:) + integer(psb_ipk_), intent(out) :: info + end subroutine psi_z_heap_get_first + end interface + + interface + subroutine psi_z_idx_heap_get_first(key,index,last,heap,idxs,dir,info) + import + complex(psb_dpk_), intent(inout) :: key + integer(psb_ipk_), intent(out) :: index + complex(psb_dpk_), intent(inout) :: heap(:) + integer(psb_ipk_), intent(in) :: dir + integer(psb_ipk_), intent(inout) :: last + integer(psb_ipk_), intent(inout) :: idxs(:) + integer(psb_ipk_), intent(out) :: info + end subroutine psi_z_idx_heap_get_first + end interface + + interface + subroutine psi_zlisrx_up(n,x,ix) + import + complex(psb_dpk_), intent(inout) :: x(:) + integer(psb_ipk_), intent(inout) :: ix(:) + integer(psb_ipk_), intent(in) :: n + end subroutine psi_zlisrx_up + subroutine psi_zlisrx_dw(n,x,ix) + import + complex(psb_dpk_), intent(inout) :: x(:) + integer(psb_ipk_), intent(inout) :: ix(:) + integer(psb_ipk_), intent(in) :: n + end subroutine psi_zlisrx_dw + subroutine psi_zlisr_up(n,x) + import + complex(psb_dpk_), intent(inout) :: x(:) + integer(psb_ipk_), intent(in) :: n + end subroutine psi_zlisr_up + subroutine psi_zlisr_dw(n,x) + import + complex(psb_dpk_), intent(inout) :: x(:) + integer(psb_ipk_), intent(in) :: n + end subroutine psi_zlisr_dw + subroutine psi_zalisrx_up(n,x,ix) + import + complex(psb_dpk_), intent(inout) :: x(:) + integer(psb_ipk_), intent(inout) :: ix(:) + integer(psb_ipk_), intent(in) :: n + end subroutine psi_zalisrx_up + subroutine psi_zalisrx_dw(n,x,ix) + import + complex(psb_dpk_), intent(inout) :: x(:) + integer(psb_ipk_), intent(inout) :: ix(:) + integer(psb_ipk_), intent(in) :: n + end subroutine psi_zalisrx_dw + subroutine psi_zalisr_up(n,x) + import + complex(psb_dpk_), intent(inout) :: x(:) + integer(psb_ipk_), intent(in) :: n + end subroutine psi_zalisr_up + subroutine psi_zalisr_dw(n,x) + import + complex(psb_dpk_), intent(inout) :: x(:) + integer(psb_ipk_), intent(in) :: n + end subroutine psi_zalisr_dw + subroutine psi_zaisrx_up(n,x,ix) + import + complex(psb_dpk_), intent(inout) :: x(:) + integer(psb_ipk_), intent(inout) :: ix(:) + integer(psb_ipk_), intent(in) :: n + end subroutine psi_zaisrx_up + subroutine psi_zaisrx_dw(n,x,ix) + import + complex(psb_dpk_), intent(inout) :: x(:) + integer(psb_ipk_), intent(inout) :: ix(:) + integer(psb_ipk_), intent(in) :: n + end subroutine psi_zaisrx_dw + subroutine psi_zaisr_up(n,x) + import + complex(psb_dpk_), intent(inout) :: x(:) + integer(psb_ipk_), intent(in) :: n + end subroutine psi_zaisr_up + subroutine psi_zaisr_dw(n,x) + import + complex(psb_dpk_), intent(inout) :: x(:) + integer(psb_ipk_), intent(in) :: n + end subroutine psi_zaisr_dw + end interface + + interface + subroutine psi_zlqsrx_up(n,x,ix) + import + complex(psb_dpk_), intent(inout) :: x(:) + integer(psb_ipk_), intent(inout) :: ix(:) + integer(psb_ipk_), intent(in) :: n + end subroutine psi_zlqsrx_up + subroutine psi_zlqsrx_dw(n,x,ix) + import + complex(psb_dpk_), intent(inout) :: x(:) + integer(psb_ipk_), intent(inout) :: ix(:) + integer(psb_ipk_), intent(in) :: n + end subroutine psi_zlqsrx_dw + subroutine psi_zlqsr_up(n,x) + import + complex(psb_dpk_), intent(inout) :: x(:) + integer(psb_ipk_), intent(in) :: n + end subroutine psi_zlqsr_up + subroutine psi_zlqsr_dw(n,x) + import + complex(psb_dpk_), intent(inout) :: x(:) + integer(psb_ipk_), intent(in) :: n + end subroutine psi_zlqsr_dw + subroutine psi_zalqsrx_up(n,x,ix) + import + complex(psb_dpk_), intent(inout) :: x(:) + integer(psb_ipk_), intent(inout) :: ix(:) + integer(psb_ipk_), intent(in) :: n + end subroutine psi_zalqsrx_up + subroutine psi_zalqsrx_dw(n,x,ix) + import + complex(psb_dpk_), intent(inout) :: x(:) + integer(psb_ipk_), intent(inout) :: ix(:) + integer(psb_ipk_), intent(in) :: n + end subroutine psi_zalqsrx_dw + subroutine psi_zalqsr_up(n,x) + import + complex(psb_dpk_), intent(inout) :: x(:) + integer(psb_ipk_), intent(in) :: n + end subroutine psi_zalqsr_up + subroutine psi_zalqsr_dw(n,x) + import + complex(psb_dpk_), intent(inout) :: x(:) + integer(psb_ipk_), intent(in) :: n + end subroutine psi_zalqsr_dw + subroutine psi_zaqsrx_up(n,x,ix) + import + complex(psb_dpk_), intent(inout) :: x(:) + integer(psb_ipk_), intent(inout) :: ix(:) + integer(psb_ipk_), intent(in) :: n + end subroutine psi_zaqsrx_up + subroutine psi_zaqsrx_dw(n,x,ix) + import + complex(psb_dpk_), intent(inout) :: x(:) + integer(psb_ipk_), intent(inout) :: ix(:) + integer(psb_ipk_), intent(in) :: n + end subroutine psi_zaqsrx_dw + subroutine psi_zaqsr_up(n,x) + import + complex(psb_dpk_), intent(inout) :: x(:) + integer(psb_ipk_), intent(in) :: n + end subroutine psi_zaqsr_up + subroutine psi_zaqsr_dw(n,x) + import + complex(psb_dpk_), intent(inout) :: x(:) + integer(psb_ipk_), intent(in) :: n + end subroutine psi_zaqsr_dw + end interface + +contains + + subroutine psb_z_init_heap(heap,info,dir) + use psb_realloc_mod, only : psb_ensure_size + implicit none + class(psb_z_heap), intent(inout) :: heap + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_), intent(in), optional :: dir + + info = psb_success_ + heap%last=0 + if (present(dir)) then + heap%dir = dir + else + heap%dir = psb_asort_up_ + endif + select case(heap%dir) + case (psb_asort_up_,psb_asort_down_) + ! ok, do nothing + case default + write(psb_err_unit,*) 'Invalid direction, defaulting to psb_asort_up_' + heap%dir = psb_asort_up_ + end select + call psb_ensure_size(psb_heap_resize,heap%keys,info) + + return + end subroutine psb_z_init_heap + + + function psb_z_howmany(heap) result(res) + implicit none + class(psb_z_heap), intent(in) :: heap + integer(psb_ipk_) :: res + res = heap%last + end function psb_z_howmany + + subroutine psb_z_insert_heap(key,heap,info) + use psb_realloc_mod, only : psb_ensure_size + implicit none + + complex(psb_dpk_), intent(in) :: key + class(psb_z_heap), intent(inout) :: heap + integer(psb_ipk_), intent(out) :: info + + info = psb_success_ + if (heap%last < 0) then + write(psb_err_unit,*) 'Invalid last in heap ',heap%last + info = heap%last + return + endif + + call psb_ensure_size(heap%last+1,heap%keys,info,addsz=psb_heap_resize) + if (info /= psb_success_) then + write(psb_err_unit,*) 'Memory allocation failure in heap_insert' + info = -5 + return + end if + call psi_z_insert_heap(key,& + & heap%last,heap%keys,heap%dir,info) + + return + end subroutine psb_z_insert_heap + + subroutine psb_z_heap_get_first(key,heap,info) + implicit none + + class(psb_z_heap), intent(inout) :: heap + integer(psb_ipk_), intent(out) :: info + complex(psb_dpk_), intent(out) :: key + + + info = psb_success_ + + call psi_z_heap_get_first(key,& + & heap%last,heap%keys,heap%dir,info) + + return + end subroutine psb_z_heap_get_first + + subroutine psb_z_dump_heap(iout,heap,info) + + implicit none + class(psb_z_heap), intent(in) :: heap + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_), intent(in) :: iout + + info = psb_success_ + if (iout < 0) then + write(psb_err_unit,*) 'Invalid file ' + info =-1 + return + end if + + write(iout,*) 'Heap direction ',heap%dir + write(iout,*) 'Heap size ',heap%last + if ((heap%last > 0).and.((.not.allocated(heap%keys)).or.& + & (size(heap%keys) 0).and.((.not.allocated(heap%keys)).or.& + & (size(heap%keys) 0).and.((.not.allocated(heap%idxs)).or.& + & (size(heap%idxs)