fdurastante
/
psblas3
mirror of https://github.com/sfilippone/psblas3.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

Reworked IP_REORD.

Browse Source
ILmat
Salvatore Filippone 8 years ago
parent 8723dd8112
commit 744731b4ed
16 changed files with 3901 additions and 23 deletions
Show Stats Download Patch File Download Diff File

6
base/modules/Makefile
Unescape Escape View File

@ -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

7
base/modules/auxil/psb_c_ip_reord_mod.f90
Unescape Escape View File

@ -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

608
base/modules/auxil/psb_c_sort_mod.f90
Unescape Escape View File

@ -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)<heap%last))) then
write(iout,*) 'Inconsistent size/allocation status!!'
else
write(iout,*) heap%keys(1:heap%last)
end if
end subroutine psb_c_dump_heap
subroutine psb_c_free_heap(heap,info)
implicit none
class(psb_c_heap), intent(inout) :: heap
integer(psb_ipk_), intent(out) :: info
info=psb_success_
if (allocated(heap%keys)) deallocate(heap%keys,stat=info)
end subroutine psb_c_free_heap
subroutine psb_c_idx_init_heap(heap,info,dir)
use psb_realloc_mod, only : psb_ensure_size
implicit none
class(psb_c_idx_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)
call psb_ensure_size(psb_heap_resize,heap%idxs,info)
return
end subroutine psb_c_idx_init_heap
function psb_c_idx_howmany(heap) result(res)
implicit none
class(psb_c_idx_heap), intent(in) :: heap
integer(psb_ipk_) :: res
res = heap%last
end function psb_c_idx_howmany
subroutine psb_c_idx_insert_heap(key,index,heap,info)
use psb_realloc_mod, only : psb_ensure_size
implicit none
complex(psb_spk_), intent(in) :: key
integer(psb_ipk_), intent(in) :: index
class(psb_c_idx_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_) &
& call psb_ensure_size(heap%last+1,heap%idxs,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_idx_insert_heap(key,index,&
& heap%last,heap%keys,heap%idxs,heap%dir,info)
return
end subroutine psb_c_idx_insert_heap
subroutine psb_c_idx_heap_get_first(key,index,heap,info)
implicit none
class(psb_c_idx_heap), intent(inout) :: heap
integer(psb_ipk_), intent(out) :: index
integer(psb_ipk_), intent(out) :: info
complex(psb_spk_), intent(out) :: key
info = psb_success_
call psi_c_idx_heap_get_first(key,index,&
& heap%last,heap%keys,heap%idxs,heap%dir,info)
return
end subroutine psb_c_idx_heap_get_first
subroutine psb_c_idx_dump_heap(iout,heap,info)
implicit none
class(psb_c_idx_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)<heap%last))) then
write(iout,*) 'Inconsistent size/allocation status!!'
else if ((heap%last > 0).and.((.not.allocated(heap%idxs)).or.&
& (size(heap%idxs)<heap%last))) then
write(iout,*) 'Inconsistent size/allocation status!!'
else
write(iout,*) heap%keys(1:heap%last)
write(iout,*) heap%idxs(1:heap%last)
end if
end subroutine psb_c_idx_dump_heap
subroutine psb_c_idx_free_heap(heap,info)
implicit none
class(psb_c_idx_heap), intent(inout) :: heap
integer(psb_ipk_), intent(out) :: info
info=psb_success_
if (allocated(heap%keys)) deallocate(heap%keys,stat=info)
if ((info == psb_success_).and.(allocated(heap%idxs))) deallocate(heap%idxs,stat=info)
end subroutine psb_c_idx_free_heap
end module psb_c_sort_mod

7
base/modules/auxil/psb_d_ip_reord_mod.f90
Unescape Escape View File

@ -53,10 +53,9 @@ contains
integer(psb_ipk_), intent(in) :: n
integer(psb_ipk_) :: iaux(0:*)
real(psb_dpk_) :: x(*)
integer(psb_ipk_) :: lswap, lp, k
real(psb_dpk_) :: swap
lp = iaux(0)
k = 1
do
@ -82,10 +81,11 @@ contains
integer(psb_ipk_) :: iaux(0:*)
real(psb_dpk_) :: x(*)
integer(psb_ipk_) :: indx(*)
integer(psb_ipk_) :: lswap, lp, k, ixswap
real(psb_dpk_) :: swap
lp = iaux(0)
k = 1
do
@ -150,7 +150,6 @@ contains
integer(psb_ipk_) :: iaux(0:*)
real(psb_dpk_) :: x(*)
integer(psb_ipk_) :: i1(*), i2(*), i3(*)
integer(psb_ipk_) :: lswap, lp, k, isw1, isw2, isw3
real(psb_dpk_) :: swap

570
base/modules/auxil/psb_d_sort_mod.f90
Unescape Escape View File

@ -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_d_sort_mod
use psb_const_mod
interface psb_msort_unique
subroutine psb_dmsort_u(x,nout,dir)
import
real(psb_dpk_), intent(inout) :: x(:)
integer(psb_ipk_), intent(out) :: nout
integer(psb_ipk_), optional, intent(in) :: dir
end subroutine psb_dmsort_u
end interface psb_msort_unique
type psb_d_heap
integer(psb_ipk_) :: last, dir
real(psb_dpk_), allocatable :: keys(:)
contains
procedure, pass(heap) :: init => 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)<heap%last))) then
write(iout,*) 'Inconsistent size/allocation status!!'
else
write(iout,*) heap%keys(1:heap%last)
end if
end subroutine psb_d_dump_heap
subroutine psb_d_free_heap(heap,info)
implicit none
class(psb_d_heap), intent(inout) :: heap
integer(psb_ipk_), intent(out) :: info
info=psb_success_
if (allocated(heap%keys)) deallocate(heap%keys,stat=info)
end subroutine psb_d_free_heap
subroutine psb_d_idx_init_heap(heap,info,dir)
use psb_realloc_mod, only : psb_ensure_size
implicit none
class(psb_d_idx_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)
call psb_ensure_size(psb_heap_resize,heap%idxs,info)
return
end subroutine psb_d_idx_init_heap
function psb_d_idx_howmany(heap) result(res)
implicit none
class(psb_d_idx_heap), intent(in) :: heap
integer(psb_ipk_) :: res
res = heap%last
end function psb_d_idx_howmany
subroutine psb_d_idx_insert_heap(key,index,heap,info)
use psb_realloc_mod, only : psb_ensure_size
implicit none
real(psb_dpk_), intent(in) :: key
integer(psb_ipk_), intent(in) :: index
class(psb_d_idx_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_) &
& call psb_ensure_size(heap%last+1,heap%idxs,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_idx_insert_heap(key,index,&
& heap%last,heap%keys,heap%idxs,heap%dir,info)
return
end subroutine psb_d_idx_insert_heap
subroutine psb_d_idx_heap_get_first(key,index,heap,info)
implicit none
class(psb_d_idx_heap), intent(inout) :: heap
integer(psb_ipk_), intent(out) :: index
integer(psb_ipk_), intent(out) :: info
real(psb_dpk_), intent(out) :: key
info = psb_success_
call psi_d_idx_heap_get_first(key,index,&
& heap%last,heap%keys,heap%idxs,heap%dir,info)
return
end subroutine psb_d_idx_heap_get_first
subroutine psb_d_idx_dump_heap(iout,heap,info)
implicit none
class(psb_d_idx_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)<heap%last))) then
write(iout,*) 'Inconsistent size/allocation status!!'
else if ((heap%last > 0).and.((.not.allocated(heap%idxs)).or.&
& (size(heap%idxs)<heap%last))) then
write(iout,*) 'Inconsistent size/allocation status!!'
else
write(iout,*) heap%keys(1:heap%last)
write(iout,*) heap%idxs(1:heap%last)
end if
end subroutine psb_d_idx_dump_heap
subroutine psb_d_idx_free_heap(heap,info)
implicit none
class(psb_d_idx_heap), intent(inout) :: heap
integer(psb_ipk_), intent(out) :: info
info=psb_success_
if (allocated(heap%keys)) deallocate(heap%keys,stat=info)
if ((info == psb_success_).and.(allocated(heap%idxs))) deallocate(heap%idxs,stat=info)
end subroutine psb_d_idx_free_heap
end module psb_d_sort_mod

186
base/modules/auxil/psb_e_ip_reord_mod.f90
Unescape Escape View File

@ -0,0 +1,186 @@
!
! 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.
!
!
!
! Reorder (an) input vector(s) based on a list sort output.
! Based on: D. E. Knuth: The Art of Computer Programming
! vol. 3: Sorting and Searching, Addison Wesley, 1973
! ex. 5.2.12
!
!
module psb_e_ip_reord_mod
use psb_const_mod
interface psb_ip_reord
module procedure psb_ip_reord_e1,&
& psb_ip_reord_e1i1, psb_ip_reord_e1i2,&
& psb_ip_reord_e1i3
end interface
contains
subroutine psb_ip_reord_e1(n,x,iaux)
integer(psb_ipk_), intent(in) :: n
integer(psb_epk_) :: iaux(0:*)
integer(psb_epk_) :: x(*)
integer(psb_epk_) :: lswap, lp, k
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
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

577
base/modules/auxil/psb_i_sort_mod.f90
Unescape Escape View File

@ -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)<heap%last))) then
write(iout,*) 'Inconsistent size/allocation status!!'
else
write(iout,*) heap%keys(1:heap%last)
end if
end subroutine psb_i_dump_heap
subroutine psb_i_free_heap(heap,info)
implicit none
class(psb_i_heap), intent(inout) :: heap
integer(psb_ipk_), intent(out) :: info
info=psb_success_
if (allocated(heap%keys)) deallocate(heap%keys,stat=info)
end subroutine psb_i_free_heap
subroutine psb_i_idx_init_heap(heap,info,dir)
use psb_realloc_mod, only : psb_ensure_size
implicit none
class(psb_i_idx_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)
call psb_ensure_size(psb_heap_resize,heap%idxs,info)
return
end subroutine psb_i_idx_init_heap
function psb_i_idx_howmany(heap) result(res)
implicit none
class(psb_i_idx_heap), intent(in) :: heap
integer(psb_ipk_) :: res
res = heap%last
end function psb_i_idx_howmany
subroutine psb_i_idx_insert_heap(key,index,heap,info)
use psb_realloc_mod, only : psb_ensure_size
implicit none
integer(psb_ipk_), intent(in) :: key
integer(psb_ipk_), intent(in) :: index
class(psb_i_idx_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_) &
& call psb_ensure_size(heap%last+1,heap%idxs,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_idx_insert_heap(key,index,&
& heap%last,heap%keys,heap%idxs,heap%dir,info)
return
end subroutine psb_i_idx_insert_heap
subroutine psb_i_idx_heap_get_first(key,index,heap,info)
implicit none
class(psb_i_idx_heap), intent(inout) :: heap
integer(psb_ipk_), intent(out) :: index
integer(psb_ipk_), intent(out) :: info
integer(psb_ipk_), intent(out) :: key
info = psb_success_
call psi_i_idx_heap_get_first(key,index,&
& heap%last,heap%keys,heap%idxs,heap%dir,info)
return
end subroutine psb_i_idx_heap_get_first
subroutine psb_i_idx_dump_heap(iout,heap,info)
implicit none
class(psb_i_idx_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)<heap%last))) then
write(iout,*) 'Inconsistent size/allocation status!!'
else if ((heap%last > 0).and.((.not.allocated(heap%idxs)).or.&
& (size(heap%idxs)<heap%last))) then
write(iout,*) 'Inconsistent size/allocation status!!'
else
write(iout,*) heap%keys(1:heap%last)
write(iout,*) heap%idxs(1:heap%last)
end if
end subroutine psb_i_idx_dump_heap
subroutine psb_i_idx_free_heap(heap,info)
implicit none
class(psb_i_idx_heap), intent(inout) :: heap
integer(psb_ipk_), intent(out) :: info
info=psb_success_
if (allocated(heap%keys)) deallocate(heap%keys,stat=info)
if ((info == psb_success_).and.(allocated(heap%idxs))) deallocate(heap%idxs,stat=info)
end subroutine psb_i_idx_free_heap
end module psb_i_sort_mod

4
base/modules/auxil/psb_ip_reord_mod.f90
Unescape Escape View File

@ -37,8 +37,8 @@
!
!
module psb_ip_reord_mod
use psb_i_ip_reord_mod
use psb_l_ip_reord_mod
use psb_m_ip_reord_mod
use psb_e_ip_reord_mod
use psb_s_ip_reord_mod
use psb_c_ip_reord_mod
use psb_d_ip_reord_mod

577
base/modules/auxil/psb_l_sort_mod.f90
Unescape Escape View File

@ -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_l_sort_mod
use psb_const_mod
interface psb_isaperm
logical function psb_lisaperm(n,eip)
import
integer(psb_lpk_), intent(in) :: n
integer(psb_lpk_), intent(in) :: eip(n)
end function psb_lisaperm
end interface psb_isaperm
interface psb_msort_unique
subroutine psb_lmsort_u(x,nout,dir)
import
integer(psb_lpk_), intent(inout) :: x(:)
integer(psb_lpk_), intent(out) :: nout
integer(psb_ipk_), optional, intent(in) :: dir
end subroutine psb_lmsort_u
end interface psb_msort_unique
type psb_l_heap
integer(psb_ipk_) :: last, dir
integer(psb_lpk_), allocatable :: keys(:)
contains
procedure, pass(heap) :: init => 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)<heap%last))) then
write(iout,*) 'Inconsistent size/allocation status!!'
else
write(iout,*) heap%keys(1:heap%last)
end if
end subroutine psb_l_dump_heap
subroutine psb_l_free_heap(heap,info)
implicit none
class(psb_l_heap), intent(inout) :: heap
integer(psb_ipk_), intent(out) :: info
info=psb_success_
if (allocated(heap%keys)) deallocate(heap%keys,stat=info)
end subroutine psb_l_free_heap
subroutine psb_l_idx_init_heap(heap,info,dir)
use psb_realloc_mod, only : psb_ensure_size
implicit none
class(psb_l_idx_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)
call psb_ensure_size(psb_heap_resize,heap%idxs,info)
return
end subroutine psb_l_idx_init_heap
function psb_l_idx_howmany(heap) result(res)
implicit none
class(psb_l_idx_heap), intent(in) :: heap
integer(psb_ipk_) :: res
res = heap%last
end function psb_l_idx_howmany
subroutine psb_l_idx_insert_heap(key,index,heap,info)
use psb_realloc_mod, only : psb_ensure_size
implicit none
integer(psb_lpk_), intent(in) :: key
integer(psb_lpk_), intent(in) :: index
class(psb_l_idx_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_) &
& call psb_ensure_size(heap%last+1,heap%idxs,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_idx_insert_heap(key,index,&
& heap%last,heap%keys,heap%idxs,heap%dir,info)
return
end subroutine psb_l_idx_insert_heap
subroutine psb_l_idx_heap_get_first(key,index,heap,info)
implicit none
class(psb_l_idx_heap), intent(inout) :: heap
integer(psb_lpk_), intent(out) :: index
integer(psb_ipk_), intent(out) :: info
integer(psb_lpk_), intent(out) :: key
info = psb_success_
call psi_l_idx_heap_get_first(key,index,&
& heap%last,heap%keys,heap%idxs,heap%dir,info)
return
end subroutine psb_l_idx_heap_get_first
subroutine psb_l_idx_dump_heap(iout,heap,info)
implicit none
class(psb_l_idx_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)<heap%last))) then
write(iout,*) 'Inconsistent size/allocation status!!'
else if ((heap%last > 0).and.((.not.allocated(heap%idxs)).or.&
& (size(heap%idxs)<heap%last))) then
write(iout,*) 'Inconsistent size/allocation status!!'
else
write(iout,*) heap%keys(1:heap%last)
write(iout,*) heap%idxs(1:heap%last)
end if
end subroutine psb_l_idx_dump_heap
subroutine psb_l_idx_free_heap(heap,info)
implicit none
class(psb_l_idx_heap), intent(inout) :: heap
integer(psb_ipk_), intent(out) :: info
info=psb_success_
if (allocated(heap%keys)) deallocate(heap%keys,stat=info)
if ((info == psb_success_).and.(allocated(heap%idxs))) deallocate(heap%idxs,stat=info)
end subroutine psb_l_idx_free_heap
end module psb_l_sort_mod

186
base/modules/auxil/psb_m_ip_reord_mod.f90
Unescape Escape View File

@ -0,0 +1,186 @@
!
! 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.
!
!
!
! Reorder (an) input vector(s) based on a list sort output.
! Based on: D. E. Knuth: The Art of Computer Programming
! vol. 3: Sorting and Searching, Addison Wesley, 1973
! ex. 5.2.12
!
!
module psb_m_ip_reord_mod
use psb_const_mod
interface psb_ip_reord
module procedure psb_ip_reord_m1,&
& psb_ip_reord_m1i1, psb_ip_reord_m1i2,&
& psb_ip_reord_m1i3
end interface
contains
subroutine psb_ip_reord_m1(n,x,iaux)
integer(psb_ipk_), intent(in) :: n
integer(psb_mpk_) :: iaux(0:*)
integer(psb_mpk_) :: x(*)
integer(psb_mpk_) :: lswap, lp, k
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
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

7
base/modules/auxil/psb_s_ip_reord_mod.f90
Unescape Escape View File

@ -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

570
base/modules/auxil/psb_s_sort_mod.f90
Unescape Escape View File

@ -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)<heap%last))) then
write(iout,*) 'Inconsistent size/allocation status!!'
else
write(iout,*) heap%keys(1:heap%last)
end if
end subroutine psb_s_dump_heap
subroutine psb_s_free_heap(heap,info)
implicit none
class(psb_s_heap), intent(inout) :: heap
integer(psb_ipk_), intent(out) :: info
info=psb_success_
if (allocated(heap%keys)) deallocate(heap%keys,stat=info)
end subroutine psb_s_free_heap
subroutine psb_s_idx_init_heap(heap,info,dir)
use psb_realloc_mod, only : psb_ensure_size
implicit none
class(psb_s_idx_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)
call psb_ensure_size(psb_heap_resize,heap%idxs,info)
return
end subroutine psb_s_idx_init_heap
function psb_s_idx_howmany(heap) result(res)
implicit none
class(psb_s_idx_heap), intent(in) :: heap
integer(psb_ipk_) :: res
res = heap%last
end function psb_s_idx_howmany
subroutine psb_s_idx_insert_heap(key,index,heap,info)
use psb_realloc_mod, only : psb_ensure_size
implicit none
real(psb_spk_), intent(in) :: key
integer(psb_ipk_), intent(in) :: index
class(psb_s_idx_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_) &
& call psb_ensure_size(heap%last+1,heap%idxs,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_idx_insert_heap(key,index,&
& heap%last,heap%keys,heap%idxs,heap%dir,info)
return
end subroutine psb_s_idx_insert_heap
subroutine psb_s_idx_heap_get_first(key,index,heap,info)
implicit none
class(psb_s_idx_heap), intent(inout) :: heap
integer(psb_ipk_), intent(out) :: index
integer(psb_ipk_), intent(out) :: info
real(psb_spk_), intent(out) :: key
info = psb_success_
call psi_s_idx_heap_get_first(key,index,&
& heap%last,heap%keys,heap%idxs,heap%dir,info)
return
end subroutine psb_s_idx_heap_get_first
subroutine psb_s_idx_dump_heap(iout,heap,info)
implicit none
class(psb_s_idx_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)<heap%last))) then
write(iout,*) 'Inconsistent size/allocation status!!'
else if ((heap%last > 0).and.((.not.allocated(heap%idxs)).or.&
& (size(heap%idxs)<heap%last))) then
write(iout,*) 'Inconsistent size/allocation status!!'
else
write(iout,*) heap%keys(1:heap%last)
write(iout,*) heap%idxs(1:heap%last)
end if
end subroutine psb_s_idx_dump_heap
subroutine psb_s_idx_free_heap(heap,info)
implicit none
class(psb_s_idx_heap), intent(inout) :: heap
integer(psb_ipk_), intent(out) :: info
info=psb_success_
if (allocated(heap%keys)) deallocate(heap%keys,stat=info)
if ((info == psb_success_).and.(allocated(heap%idxs))) deallocate(heap%idxs,stat=info)
end subroutine psb_s_idx_free_heap
end module psb_s_sort_mod

7
base/modules/auxil/psb_z_ip_reord_mod.f90
Unescape Escape View File

@ -53,10 +53,9 @@ contains
integer(psb_ipk_), intent(in) :: n
integer(psb_ipk_) :: iaux(0:*)
complex(psb_dpk_) :: x(*)
integer(psb_ipk_) :: lswap, lp, k
complex(psb_dpk_) :: swap
lp = iaux(0)
k = 1
do
@ -82,10 +81,11 @@ contains
integer(psb_ipk_) :: iaux(0:*)
complex(psb_dpk_) :: x(*)
integer(psb_ipk_) :: indx(*)
integer(psb_ipk_) :: lswap, lp, k, ixswap
complex(psb_dpk_) :: swap
lp = iaux(0)
k = 1
do
@ -150,7 +150,6 @@ contains
integer(psb_ipk_) :: iaux(0:*)
complex(psb_dpk_) :: x(*)
integer(psb_ipk_) :: i1(*), i2(*), i3(*)
integer(psb_ipk_) :: lswap, lp, k, isw1, isw2, isw3
complex(psb_dpk_) :: swap

608
base/modules/auxil/psb_z_sort_mod.f90
Unescape Escape View File

@ -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_z_sort_mod
use psb_const_mod
interface psb_msort_unique
subroutine psb_zmsort_u(x,nout,dir)
import
complex(psb_dpk_), intent(inout) :: x(:)
integer(psb_ipk_), intent(out) :: nout
integer(psb_ipk_), optional, intent(in) :: dir
end subroutine psb_zmsort_u
end interface psb_msort_unique
type psb_z_heap
integer(psb_ipk_) :: last, dir
complex(psb_dpk_), allocatable :: keys(:)
contains
procedure, pass(heap) :: init => 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)<heap%last))) then
write(iout,*) 'Inconsistent size/allocation status!!'
else
write(iout,*) heap%keys(1:heap%last)
end if
end subroutine psb_z_dump_heap
subroutine psb_z_free_heap(heap,info)
implicit none
class(psb_z_heap), intent(inout) :: heap
integer(psb_ipk_), intent(out) :: info
info=psb_success_
if (allocated(heap%keys)) deallocate(heap%keys,stat=info)
end subroutine psb_z_free_heap
subroutine psb_z_idx_init_heap(heap,info,dir)
use psb_realloc_mod, only : psb_ensure_size
implicit none
class(psb_z_idx_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)
call psb_ensure_size(psb_heap_resize,heap%idxs,info)
return
end subroutine psb_z_idx_init_heap
function psb_z_idx_howmany(heap) result(res)
implicit none
class(psb_z_idx_heap), intent(in) :: heap
integer(psb_ipk_) :: res
res = heap%last
end function psb_z_idx_howmany
subroutine psb_z_idx_insert_heap(key,index,heap,info)
use psb_realloc_mod, only : psb_ensure_size
implicit none
complex(psb_dpk_), intent(in) :: key
integer(psb_ipk_), intent(in) :: index
class(psb_z_idx_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_) &
& call psb_ensure_size(heap%last+1,heap%idxs,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_idx_insert_heap(key,index,&
& heap%last,heap%keys,heap%idxs,heap%dir,info)
return
end subroutine psb_z_idx_insert_heap
subroutine psb_z_idx_heap_get_first(key,index,heap,info)
implicit none
class(psb_z_idx_heap), intent(inout) :: heap
integer(psb_ipk_), intent(out) :: index
integer(psb_ipk_), intent(out) :: info
complex(psb_dpk_), intent(out) :: key
info = psb_success_
call psi_z_idx_heap_get_first(key,index,&
& heap%last,heap%keys,heap%idxs,heap%dir,info)
return
end subroutine psb_z_idx_heap_get_first
subroutine psb_z_idx_dump_heap(iout,heap,info)
implicit none
class(psb_z_idx_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)<heap%last))) then
write(iout,*) 'Inconsistent size/allocation status!!'
else if ((heap%last > 0).and.((.not.allocated(heap%idxs)).or.&
& (size(heap%idxs)<heap%last))) then
write(iout,*) 'Inconsistent size/allocation status!!'
else
write(iout,*) heap%keys(1:heap%last)
write(iout,*) heap%idxs(1:heap%last)
end if
end subroutine psb_z_idx_dump_heap
subroutine psb_z_idx_free_heap(heap,info)
implicit none
class(psb_z_idx_heap), intent(inout) :: heap
integer(psb_ipk_), intent(out) :: info
info=psb_success_
if (allocated(heap%keys)) deallocate(heap%keys,stat=info)
if ((info == psb_success_).and.(allocated(heap%idxs))) deallocate(heap%idxs,stat=info)
end subroutine psb_z_idx_free_heap
end module psb_z_sort_mod

3
base/serial/sort/psb_e_msort_impl.f90
Unescape Escape View File

@ -285,8 +285,7 @@
return
end subroutine psb_emsort
end subroutine psb_emsort
subroutine psi_e_msort_up(n,k,l,iret)
use psb_const_mod

1
base/serial/sort/psb_m_msort_impl.f90
Unescape Escape View File

@ -285,7 +285,6 @@
return
end subroutine psb_mmsort
subroutine psi_m_msort_up(n,k,l,iret)

Write Preview
Loading…
Cancel
Save