diff --git a/base/internals/psi_adjcncy_fnd_owner.F90 b/base/internals/psi_adjcncy_fnd_owner.F90 index 03aa0f09..639cdb5d 100644 --- a/base/internals/psi_adjcncy_fnd_owner.F90 +++ b/base/internals/psi_adjcncy_fnd_owner.F90 @@ -71,7 +71,7 @@ subroutine psi_adjcncy_fnd_owner(idx,iprc,adj,idxmap,info) #endif integer(psb_lpk_), intent(in) :: idx(:) integer(psb_ipk_), allocatable, intent(out) :: iprc(:) - integer(psb_ipk_), intent(in) :: adj(:) + integer(psb_ipk_), intent(inout) :: adj(:) class(psb_indx_map), intent(in) :: idxmap integer(psb_ipk_), intent(out) :: info diff --git a/base/internals/psi_sort_dl.f90 b/base/internals/psi_sort_dl.f90 index 55cc280f..ef3ac74d 100644 --- a/base/internals/psi_sort_dl.f90 +++ b/base/internals/psi_sort_dl.f90 @@ -84,9 +84,10 @@ subroutine psi_i_csr_sort_dl(dl_ptr,c_dep_list,l_dep_list,ctxt,info) use psb_sort_mod implicit none - integer(psb_ipk_), intent(inout) :: c_dep_list(:), dl_ptr(0:), l_dep_list(0:) + integer(psb_ipk_), intent(in) :: dl_ptr(0:) + integer(psb_ipk_), intent(inout) :: c_dep_list(:), l_dep_list(0:) type(psb_ctxt_type), intent(in) :: ctxt - integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_), intent(out) :: info ! Local variables integer(psb_ipk_), allocatable :: dg(:), dgp(:),& & idx(:), upd(:), edges(:,:), ich(:) diff --git a/base/modules/auxil/psb_c_hsort_mod.f90 b/base/modules/auxil/psb_c_hsort_mod.f90 index 43b22ec4..e7eb2fbf 100644 --- a/base/modules/auxil/psb_c_hsort_mod.f90 +++ b/base/modules/auxil/psb_c_hsort_mod.f90 @@ -100,8 +100,8 @@ module psb_c_hsort_mod 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 + 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 diff --git a/base/modules/auxil/psb_c_hsort_x_mod.f90 b/base/modules/auxil/psb_c_hsort_x_mod.f90 index 4331c567..c0e39411 100644 --- a/base/modules/auxil/psb_c_hsort_x_mod.f90 +++ b/base/modules/auxil/psb_c_hsort_x_mod.f90 @@ -45,7 +45,8 @@ module psb_c_hsort_x_mod use psb_c_hsort_mod type psb_c_heap - integer(psb_ipk_) :: last, dir + integer(psb_ipk_) :: dir + integer(psb_ipk_) :: last complex(psb_spk_), allocatable :: keys(:) contains procedure, pass(heap) :: init => psb_c_init_heap @@ -57,7 +58,8 @@ module psb_c_hsort_x_mod end type psb_c_heap type psb_c_idx_heap - integer(psb_ipk_) :: last, dir + integer(psb_ipk_) :: dir + integer(psb_ipk_) :: last complex(psb_spk_), allocatable :: keys(:) integer(psb_ipk_), allocatable :: idxs(:) contains @@ -121,7 +123,7 @@ contains return endif - call psb_ensure_size(heap%last+1,heap%keys,info,addsz=psb_heap_resize) + call psb_ensure_size(heap%last+1,heap%keys,info,addsz=(1_psb_ipk_)*psb_heap_resize) if (info /= psb_success_) then write(psb_err_unit,*) 'Memory allocation failure in heap_insert' info = -5 @@ -234,9 +236,9 @@ contains return endif - call psb_ensure_size(heap%last+1,heap%keys,info,addsz=psb_heap_resize) + call psb_ensure_size(heap%last+1,heap%keys,info,addsz=(1_psb_ipk_)*psb_heap_resize) if (info == psb_success_) & - & call psb_ensure_size(heap%last+1,heap%idxs,info,addsz=psb_heap_resize) + & call psb_ensure_size(heap%last+1,heap%idxs,info,addsz=(1_psb_ipk_)*psb_heap_resize) if (info /= psb_success_) then write(psb_err_unit,*) 'Memory allocation failure in heap_insert' info = -5 @@ -254,7 +256,7 @@ contains 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 + complex(psb_spk_), intent(inout) :: key info = psb_success_ diff --git a/base/modules/auxil/psb_d_hsort_mod.f90 b/base/modules/auxil/psb_d_hsort_mod.f90 index c1a7523c..570f27b1 100644 --- a/base/modules/auxil/psb_d_hsort_mod.f90 +++ b/base/modules/auxil/psb_d_hsort_mod.f90 @@ -100,8 +100,8 @@ module psb_d_hsort_mod 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 + 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 diff --git a/base/modules/auxil/psb_d_hsort_x_mod.f90 b/base/modules/auxil/psb_d_hsort_x_mod.f90 index df290e38..7273e972 100644 --- a/base/modules/auxil/psb_d_hsort_x_mod.f90 +++ b/base/modules/auxil/psb_d_hsort_x_mod.f90 @@ -45,7 +45,8 @@ module psb_d_hsort_x_mod use psb_d_hsort_mod type psb_d_heap - integer(psb_ipk_) :: last, dir + integer(psb_ipk_) :: dir + integer(psb_ipk_) :: last real(psb_dpk_), allocatable :: keys(:) contains procedure, pass(heap) :: init => psb_d_init_heap @@ -57,7 +58,8 @@ module psb_d_hsort_x_mod end type psb_d_heap type psb_d_idx_heap - integer(psb_ipk_) :: last, dir + integer(psb_ipk_) :: dir + integer(psb_ipk_) :: last real(psb_dpk_), allocatable :: keys(:) integer(psb_ipk_), allocatable :: idxs(:) contains @@ -121,7 +123,7 @@ contains return endif - call psb_ensure_size(heap%last+1,heap%keys,info,addsz=psb_heap_resize) + call psb_ensure_size(heap%last+1,heap%keys,info,addsz=(1_psb_ipk_)*psb_heap_resize) if (info /= psb_success_) then write(psb_err_unit,*) 'Memory allocation failure in heap_insert' info = -5 @@ -234,9 +236,9 @@ contains return endif - call psb_ensure_size(heap%last+1,heap%keys,info,addsz=psb_heap_resize) + call psb_ensure_size(heap%last+1,heap%keys,info,addsz=(1_psb_ipk_)*psb_heap_resize) if (info == psb_success_) & - & call psb_ensure_size(heap%last+1,heap%idxs,info,addsz=psb_heap_resize) + & call psb_ensure_size(heap%last+1,heap%idxs,info,addsz=(1_psb_ipk_)*psb_heap_resize) if (info /= psb_success_) then write(psb_err_unit,*) 'Memory allocation failure in heap_insert' info = -5 @@ -254,7 +256,7 @@ contains 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 + real(psb_dpk_), intent(inout) :: key info = psb_success_ diff --git a/base/modules/auxil/psb_e_hsort_mod.f90 b/base/modules/auxil/psb_e_hsort_mod.f90 index 3ce7ac45..5433cec4 100644 --- a/base/modules/auxil/psb_e_hsort_mod.f90 +++ b/base/modules/auxil/psb_e_hsort_mod.f90 @@ -67,8 +67,8 @@ module psb_e_hsort_mod integer(psb_epk_), intent(in) :: key integer(psb_epk_), intent(inout) :: heap(:) - integer(psb_ipk_), intent(in) :: dir - integer(psb_ipk_), intent(inout) :: last + integer(psb_epk_), intent(in) :: dir + integer(psb_epk_), intent(inout) :: last integer(psb_ipk_), intent(out) :: info end subroutine psi_e_insert_heap end interface psi_insert_heap @@ -88,9 +88,9 @@ module psb_e_hsort_mod integer(psb_epk_), intent(in) :: key integer(psb_epk_), intent(inout) :: heap(:) integer(psb_epk_), intent(in) :: index - integer(psb_ipk_), intent(in) :: dir + integer(psb_epk_), intent(in) :: dir integer(psb_epk_), intent(inout) :: idxs(:) - integer(psb_ipk_), intent(inout) :: last + integer(psb_epk_), intent(inout) :: last integer(psb_ipk_), intent(out) :: info end subroutine psi_e_idx_insert_heap end interface psi_idx_insert_heap @@ -100,9 +100,9 @@ module psb_e_hsort_mod subroutine psi_e_heap_get_first(key,last,heap,dir,info) import implicit none - integer(psb_epk_), intent(inout) :: key - integer(psb_ipk_), intent(inout) :: last - integer(psb_ipk_), intent(in) :: dir + integer(psb_epk_), intent(inout) :: key + integer(psb_epk_), intent(inout) :: last + integer(psb_epk_), intent(in) :: dir integer(psb_epk_), intent(inout) :: heap(:) integer(psb_ipk_), intent(out) :: info end subroutine psi_e_heap_get_first @@ -114,8 +114,8 @@ module psb_e_hsort_mod integer(psb_epk_), intent(inout) :: key integer(psb_epk_), intent(out) :: index integer(psb_epk_), intent(inout) :: heap(:) - integer(psb_ipk_), intent(in) :: dir - integer(psb_ipk_), intent(inout) :: last + integer(psb_epk_), intent(in) :: dir + integer(psb_epk_), intent(inout) :: last integer(psb_epk_), intent(inout) :: idxs(:) integer(psb_ipk_), intent(out) :: info end subroutine psi_e_idx_heap_get_first diff --git a/base/modules/auxil/psb_i2_hsort_mod.f90 b/base/modules/auxil/psb_i2_hsort_mod.f90 index 0878f86e..6a6c96de 100644 --- a/base/modules/auxil/psb_i2_hsort_mod.f90 +++ b/base/modules/auxil/psb_i2_hsort_mod.f90 @@ -100,8 +100,8 @@ module psb_i2_hsort_mod subroutine psi_i2_heap_get_first(key,last,heap,dir,info) import implicit none - integer(psb_i2pk_), intent(inout) :: key - integer(psb_ipk_), intent(inout) :: last + integer(psb_i2pk_), intent(inout) :: key + integer(psb_ipk_), intent(inout) :: last integer(psb_ipk_), intent(in) :: dir integer(psb_i2pk_), intent(inout) :: heap(:) integer(psb_ipk_), intent(out) :: info diff --git a/base/modules/auxil/psb_i_hsort_x_mod.f90 b/base/modules/auxil/psb_i_hsort_x_mod.f90 index cf148d20..0d1288a6 100644 --- a/base/modules/auxil/psb_i_hsort_x_mod.f90 +++ b/base/modules/auxil/psb_i_hsort_x_mod.f90 @@ -46,7 +46,8 @@ module psb_i_hsort_x_mod use psb_m_hsort_mod type psb_i_heap - integer(psb_ipk_) :: last, dir + integer(psb_ipk_) :: dir + integer(psb_ipk_) :: last integer(psb_ipk_), allocatable :: keys(:) contains procedure, pass(heap) :: init => psb_i_init_heap @@ -58,7 +59,8 @@ module psb_i_hsort_x_mod end type psb_i_heap type psb_i_idx_heap - integer(psb_ipk_) :: last, dir + integer(psb_ipk_) :: dir + integer(psb_ipk_) :: last integer(psb_ipk_), allocatable :: keys(:) integer(psb_ipk_), allocatable :: idxs(:) contains @@ -122,7 +124,7 @@ contains return endif - call psb_ensure_size(heap%last+1,heap%keys,info,addsz=psb_heap_resize) + call psb_ensure_size(heap%last+1,heap%keys,info,addsz=(1_psb_ipk_)*psb_heap_resize) if (info /= psb_success_) then write(psb_err_unit,*) 'Memory allocation failure in heap_insert' info = -5 @@ -235,9 +237,9 @@ contains return endif - call psb_ensure_size(heap%last+1,heap%keys,info,addsz=psb_heap_resize) + call psb_ensure_size(heap%last+1,heap%keys,info,addsz=(1_psb_ipk_)*psb_heap_resize) if (info == psb_success_) & - & call psb_ensure_size(heap%last+1,heap%idxs,info,addsz=psb_heap_resize) + & call psb_ensure_size(heap%last+1,heap%idxs,info,addsz=(1_psb_ipk_)*psb_heap_resize) if (info /= psb_success_) then write(psb_err_unit,*) 'Memory allocation failure in heap_insert' info = -5 @@ -255,7 +257,7 @@ contains 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 + integer(psb_ipk_), intent(inout) :: key info = psb_success_ diff --git a/base/modules/auxil/psb_l_hsort_x_mod.f90 b/base/modules/auxil/psb_l_hsort_x_mod.f90 index bf9bd435..487e8ce9 100644 --- a/base/modules/auxil/psb_l_hsort_x_mod.f90 +++ b/base/modules/auxil/psb_l_hsort_x_mod.f90 @@ -46,7 +46,8 @@ module psb_l_hsort_x_mod use psb_m_hsort_mod type psb_l_heap - integer(psb_ipk_) :: last, dir + integer(psb_lpk_) :: dir + integer(psb_lpk_) :: last integer(psb_lpk_), allocatable :: keys(:) contains procedure, pass(heap) :: init => psb_l_init_heap @@ -58,7 +59,8 @@ module psb_l_hsort_x_mod end type psb_l_heap type psb_l_idx_heap - integer(psb_ipk_) :: last, dir + integer(psb_lpk_) :: dir + integer(psb_lpk_) :: last integer(psb_lpk_), allocatable :: keys(:) integer(psb_lpk_), allocatable :: idxs(:) contains @@ -122,7 +124,7 @@ contains return endif - call psb_ensure_size(heap%last+1,heap%keys,info,addsz=psb_heap_resize) + call psb_ensure_size(heap%last+1,heap%keys,info,addsz=(1_psb_lpk_)*psb_heap_resize) if (info /= psb_success_) then write(psb_err_unit,*) 'Memory allocation failure in heap_insert' info = -5 @@ -235,9 +237,9 @@ contains return endif - call psb_ensure_size(heap%last+1,heap%keys,info,addsz=psb_heap_resize) + call psb_ensure_size(heap%last+1,heap%keys,info,addsz=(1_psb_lpk_)*psb_heap_resize) if (info == psb_success_) & - & call psb_ensure_size(heap%last+1,heap%idxs,info,addsz=psb_heap_resize) + & call psb_ensure_size(heap%last+1,heap%idxs,info,addsz=(1_psb_lpk_)*psb_heap_resize) if (info /= psb_success_) then write(psb_err_unit,*) 'Memory allocation failure in heap_insert' info = -5 @@ -255,7 +257,7 @@ contains 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 + integer(psb_lpk_), intent(inout) :: key info = psb_success_ diff --git a/base/modules/auxil/psb_m_hsort_mod.f90 b/base/modules/auxil/psb_m_hsort_mod.f90 index e2d7aef6..22ffdf34 100644 --- a/base/modules/auxil/psb_m_hsort_mod.f90 +++ b/base/modules/auxil/psb_m_hsort_mod.f90 @@ -100,8 +100,8 @@ module psb_m_hsort_mod subroutine psi_m_heap_get_first(key,last,heap,dir,info) import implicit none - integer(psb_mpk_), intent(inout) :: key - integer(psb_ipk_), intent(inout) :: last + integer(psb_mpk_), intent(inout) :: key + integer(psb_ipk_), intent(inout) :: last integer(psb_ipk_), intent(in) :: dir integer(psb_mpk_), intent(inout) :: heap(:) integer(psb_ipk_), intent(out) :: info diff --git a/base/modules/auxil/psb_s_hsort_mod.f90 b/base/modules/auxil/psb_s_hsort_mod.f90 index 3d30bafd..4cee5508 100644 --- a/base/modules/auxil/psb_s_hsort_mod.f90 +++ b/base/modules/auxil/psb_s_hsort_mod.f90 @@ -100,8 +100,8 @@ module psb_s_hsort_mod 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 + 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 diff --git a/base/modules/auxil/psb_s_hsort_x_mod.f90 b/base/modules/auxil/psb_s_hsort_x_mod.f90 index 3b395fb1..34f69ea4 100644 --- a/base/modules/auxil/psb_s_hsort_x_mod.f90 +++ b/base/modules/auxil/psb_s_hsort_x_mod.f90 @@ -45,7 +45,8 @@ module psb_s_hsort_x_mod use psb_s_hsort_mod type psb_s_heap - integer(psb_ipk_) :: last, dir + integer(psb_ipk_) :: dir + integer(psb_ipk_) :: last real(psb_spk_), allocatable :: keys(:) contains procedure, pass(heap) :: init => psb_s_init_heap @@ -57,7 +58,8 @@ module psb_s_hsort_x_mod end type psb_s_heap type psb_s_idx_heap - integer(psb_ipk_) :: last, dir + integer(psb_ipk_) :: dir + integer(psb_ipk_) :: last real(psb_spk_), allocatable :: keys(:) integer(psb_ipk_), allocatable :: idxs(:) contains @@ -121,7 +123,7 @@ contains return endif - call psb_ensure_size(heap%last+1,heap%keys,info,addsz=psb_heap_resize) + call psb_ensure_size(heap%last+1,heap%keys,info,addsz=(1_psb_ipk_)*psb_heap_resize) if (info /= psb_success_) then write(psb_err_unit,*) 'Memory allocation failure in heap_insert' info = -5 @@ -234,9 +236,9 @@ contains return endif - call psb_ensure_size(heap%last+1,heap%keys,info,addsz=psb_heap_resize) + call psb_ensure_size(heap%last+1,heap%keys,info,addsz=(1_psb_ipk_)*psb_heap_resize) if (info == psb_success_) & - & call psb_ensure_size(heap%last+1,heap%idxs,info,addsz=psb_heap_resize) + & call psb_ensure_size(heap%last+1,heap%idxs,info,addsz=(1_psb_ipk_)*psb_heap_resize) if (info /= psb_success_) then write(psb_err_unit,*) 'Memory allocation failure in heap_insert' info = -5 @@ -254,7 +256,7 @@ contains 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 + real(psb_spk_), intent(inout) :: key info = psb_success_ diff --git a/base/modules/auxil/psb_z_hsort_mod.f90 b/base/modules/auxil/psb_z_hsort_mod.f90 index 573eef6b..98e47da2 100644 --- a/base/modules/auxil/psb_z_hsort_mod.f90 +++ b/base/modules/auxil/psb_z_hsort_mod.f90 @@ -100,8 +100,8 @@ module psb_z_hsort_mod 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 + 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 diff --git a/base/modules/auxil/psb_z_hsort_x_mod.f90 b/base/modules/auxil/psb_z_hsort_x_mod.f90 index fc0edd9b..39f52e4f 100644 --- a/base/modules/auxil/psb_z_hsort_x_mod.f90 +++ b/base/modules/auxil/psb_z_hsort_x_mod.f90 @@ -45,7 +45,8 @@ module psb_z_hsort_x_mod use psb_z_hsort_mod type psb_z_heap - integer(psb_ipk_) :: last, dir + integer(psb_ipk_) :: dir + integer(psb_ipk_) :: last complex(psb_dpk_), allocatable :: keys(:) contains procedure, pass(heap) :: init => psb_z_init_heap @@ -57,7 +58,8 @@ module psb_z_hsort_x_mod end type psb_z_heap type psb_z_idx_heap - integer(psb_ipk_) :: last, dir + integer(psb_ipk_) :: dir + integer(psb_ipk_) :: last complex(psb_dpk_), allocatable :: keys(:) integer(psb_ipk_), allocatable :: idxs(:) contains @@ -121,7 +123,7 @@ contains return endif - call psb_ensure_size(heap%last+1,heap%keys,info,addsz=psb_heap_resize) + call psb_ensure_size(heap%last+1,heap%keys,info,addsz=(1_psb_ipk_)*psb_heap_resize) if (info /= psb_success_) then write(psb_err_unit,*) 'Memory allocation failure in heap_insert' info = -5 @@ -234,9 +236,9 @@ contains return endif - call psb_ensure_size(heap%last+1,heap%keys,info,addsz=psb_heap_resize) + call psb_ensure_size(heap%last+1,heap%keys,info,addsz=(1_psb_ipk_)*psb_heap_resize) if (info == psb_success_) & - & call psb_ensure_size(heap%last+1,heap%idxs,info,addsz=psb_heap_resize) + & call psb_ensure_size(heap%last+1,heap%idxs,info,addsz=(1_psb_ipk_)*psb_heap_resize) if (info /= psb_success_) then write(psb_err_unit,*) 'Memory allocation failure in heap_insert' info = -5 @@ -254,7 +256,7 @@ contains 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 + complex(psb_dpk_), intent(inout) :: key info = psb_success_ diff --git a/base/modules/auxil/psi_c_serial_mod.f90 b/base/modules/auxil/psi_c_serial_mod.f90 index f017f350..191e7ef3 100644 --- a/base/modules/auxil/psi_c_serial_mod.f90 +++ b/base/modules/auxil/psi_c_serial_mod.f90 @@ -93,7 +93,7 @@ module psi_c_serial_mod integer(psb_ipk_), intent(in) :: m complex(psb_spk_), intent (in) :: x(:) complex(psb_spk_), intent (in) :: y(:) - complex(psb_spk_), intent (in) :: z(:) + complex(psb_spk_), intent (inout) :: z(:) complex(psb_spk_), intent (in) :: alpha, beta integer(psb_ipk_), intent(out) :: info end subroutine psi_caxpbyv2 diff --git a/base/modules/auxil/psi_d_serial_mod.f90 b/base/modules/auxil/psi_d_serial_mod.f90 index c27aa600..f1dbc16c 100644 --- a/base/modules/auxil/psi_d_serial_mod.f90 +++ b/base/modules/auxil/psi_d_serial_mod.f90 @@ -93,7 +93,7 @@ module psi_d_serial_mod integer(psb_ipk_), intent(in) :: m real(psb_dpk_), intent (in) :: x(:) real(psb_dpk_), intent (in) :: y(:) - real(psb_dpk_), intent (in) :: z(:) + real(psb_dpk_), intent (inout) :: z(:) real(psb_dpk_), intent (in) :: alpha, beta integer(psb_ipk_), intent(out) :: info end subroutine psi_daxpbyv2 diff --git a/base/modules/auxil/psi_e_serial_mod.f90 b/base/modules/auxil/psi_e_serial_mod.f90 index 99a91985..909b025c 100644 --- a/base/modules/auxil/psi_e_serial_mod.f90 +++ b/base/modules/auxil/psi_e_serial_mod.f90 @@ -93,7 +93,7 @@ module psi_e_serial_mod integer(psb_ipk_), intent(in) :: m integer(psb_epk_), intent (in) :: x(:) integer(psb_epk_), intent (in) :: y(:) - integer(psb_epk_), intent (in) :: z(:) + integer(psb_epk_), intent (inout) :: z(:) integer(psb_epk_), intent (in) :: alpha, beta integer(psb_ipk_), intent(out) :: info end subroutine psi_eaxpbyv2 diff --git a/base/modules/auxil/psi_i2_serial_mod.f90 b/base/modules/auxil/psi_i2_serial_mod.f90 index 565955e7..38ef9c38 100644 --- a/base/modules/auxil/psi_i2_serial_mod.f90 +++ b/base/modules/auxil/psi_i2_serial_mod.f90 @@ -93,7 +93,7 @@ module psi_i2_serial_mod integer(psb_ipk_), intent(in) :: m integer(psb_i2pk_), intent (in) :: x(:) integer(psb_i2pk_), intent (in) :: y(:) - integer(psb_i2pk_), intent (in) :: z(:) + integer(psb_i2pk_), intent (inout) :: z(:) integer(psb_i2pk_), intent (in) :: alpha, beta integer(psb_ipk_), intent(out) :: info end subroutine psi_i2axpbyv2 diff --git a/base/modules/auxil/psi_m_serial_mod.f90 b/base/modules/auxil/psi_m_serial_mod.f90 index 17ea8dc4..a80d0ffe 100644 --- a/base/modules/auxil/psi_m_serial_mod.f90 +++ b/base/modules/auxil/psi_m_serial_mod.f90 @@ -93,7 +93,7 @@ module psi_m_serial_mod integer(psb_ipk_), intent(in) :: m integer(psb_mpk_), intent (in) :: x(:) integer(psb_mpk_), intent (in) :: y(:) - integer(psb_mpk_), intent (in) :: z(:) + integer(psb_mpk_), intent (inout) :: z(:) integer(psb_mpk_), intent (in) :: alpha, beta integer(psb_ipk_), intent(out) :: info end subroutine psi_maxpbyv2 diff --git a/base/modules/auxil/psi_s_serial_mod.f90 b/base/modules/auxil/psi_s_serial_mod.f90 index ed7b5d9f..3c0a3bdc 100644 --- a/base/modules/auxil/psi_s_serial_mod.f90 +++ b/base/modules/auxil/psi_s_serial_mod.f90 @@ -93,7 +93,7 @@ module psi_s_serial_mod integer(psb_ipk_), intent(in) :: m real(psb_spk_), intent (in) :: x(:) real(psb_spk_), intent (in) :: y(:) - real(psb_spk_), intent (in) :: z(:) + real(psb_spk_), intent (inout) :: z(:) real(psb_spk_), intent (in) :: alpha, beta integer(psb_ipk_), intent(out) :: info end subroutine psi_saxpbyv2 diff --git a/base/modules/auxil/psi_z_serial_mod.f90 b/base/modules/auxil/psi_z_serial_mod.f90 index 9de8451b..7bc9728e 100644 --- a/base/modules/auxil/psi_z_serial_mod.f90 +++ b/base/modules/auxil/psi_z_serial_mod.f90 @@ -93,7 +93,7 @@ module psi_z_serial_mod integer(psb_ipk_), intent(in) :: m complex(psb_dpk_), intent (in) :: x(:) complex(psb_dpk_), intent (in) :: y(:) - complex(psb_dpk_), intent (in) :: z(:) + complex(psb_dpk_), intent (inout) :: z(:) complex(psb_dpk_), intent (in) :: alpha, beta integer(psb_ipk_), intent(out) :: info end subroutine psi_zaxpbyv2 diff --git a/base/modules/desc/psb_indx_map_mod.f90 b/base/modules/desc/psb_indx_map_mod.f90 index 0c0d8199..79526d59 100644 --- a/base/modules/desc/psb_indx_map_mod.f90 +++ b/base/modules/desc/psb_indx_map_mod.f90 @@ -296,7 +296,7 @@ module psb_indx_map_mod implicit none integer(psb_lpk_), intent(in) :: idx(:) integer(psb_ipk_), allocatable, intent(out) :: iprc(:) - integer(psb_ipk_), allocatable, intent(inout) :: adj(:) + integer(psb_ipk_), intent(inout) :: adj(:) class(psb_indx_map), intent(in) :: idxmap integer(psb_ipk_), intent(out) :: info end subroutine psi_adjcncy_fnd_owner diff --git a/base/modules/psi_i_mod.F90 b/base/modules/psi_i_mod.F90 index 7310065a..e852dd99 100644 --- a/base/modules/psi_i_mod.F90 +++ b/base/modules/psi_i_mod.F90 @@ -87,10 +87,10 @@ module psi_i_mod subroutine psi_i_csr_sort_dl(dl_ptr,c_dep_list,l_dep_list,ctxt,info) import implicit none - integer(psb_ipk_), intent(in) :: c_dep_list(:), dl_ptr(0:) - integer(psb_ipk_), intent(inout) :: l_dep_list(0:) - type(psb_ctxt_type) :: ctxt - integer(psb_ipk_) :: info + integer(psb_ipk_), intent(in) :: dl_ptr(0:) + integer(psb_ipk_), intent(inout) :: c_dep_list(:), l_dep_list(0:) + type(psb_ctxt_type), intent(in) :: ctxt + integer(psb_ipk_), intent(out) :: info end subroutine psi_i_csr_sort_dl end interface diff --git a/base/psblas/psb_caxpby.f90 b/base/psblas/psb_caxpby.f90 index 9ac48603..da3dd93b 100644 --- a/base/psblas/psb_caxpby.f90 +++ b/base/psblas/psb_caxpby.f90 @@ -642,7 +642,7 @@ subroutine psb_caxpbyvout(alpha, x, beta,y, z, desc_a,info) end if if(desc_a%get_local_rows() > 0) then - call caxpby(desc_a%get_local_cols(),ione,& + call caxpbyv2(desc_a%get_local_cols(),ione,& & alpha,x,lldx,beta,& & y,lldy,z,lldz,info) end if diff --git a/base/psblas/psb_daxpby.f90 b/base/psblas/psb_daxpby.f90 index f2768789..c386f8f2 100644 --- a/base/psblas/psb_daxpby.f90 +++ b/base/psblas/psb_daxpby.f90 @@ -642,7 +642,7 @@ subroutine psb_daxpbyvout(alpha, x, beta,y, z, desc_a,info) end if if(desc_a%get_local_rows() > 0) then - call daxpby(desc_a%get_local_cols(),ione,& + call daxpbyv2(desc_a%get_local_cols(),ione,& & alpha,x,lldx,beta,& & y,lldy,z,lldz,info) end if diff --git a/base/psblas/psb_saxpby.f90 b/base/psblas/psb_saxpby.f90 index 774c1ad7..78f4d01a 100644 --- a/base/psblas/psb_saxpby.f90 +++ b/base/psblas/psb_saxpby.f90 @@ -642,7 +642,7 @@ subroutine psb_saxpbyvout(alpha, x, beta,y, z, desc_a,info) end if if(desc_a%get_local_rows() > 0) then - call saxpby(desc_a%get_local_cols(),ione,& + call saxpbyv2(desc_a%get_local_cols(),ione,& & alpha,x,lldx,beta,& & y,lldy,z,lldz,info) end if diff --git a/base/psblas/psb_zaxpby.f90 b/base/psblas/psb_zaxpby.f90 index 1165ea8a..2258f38f 100644 --- a/base/psblas/psb_zaxpby.f90 +++ b/base/psblas/psb_zaxpby.f90 @@ -642,7 +642,7 @@ subroutine psb_zaxpbyvout(alpha, x, beta,y, z, desc_a,info) end if if(desc_a%get_local_rows() > 0) then - call zaxpby(desc_a%get_local_cols(),ione,& + call zaxpbyv2(desc_a%get_local_cols(),ione,& & alpha,x,lldx,beta,& & y,lldy,z,lldz,info) end if diff --git a/base/serial/impl/psb_c_coo_impl.F90 b/base/serial/impl/psb_c_coo_impl.F90 index 442b2a48..bd03b1b8 100644 --- a/base/serial/impl/psb_c_coo_impl.F90 +++ b/base/serial/impl/psb_c_coo_impl.F90 @@ -4616,7 +4616,7 @@ function psb_lc_coo_csnm1(a) result(res) use psb_c_base_mat_mod, psb_protect_name => psb_lc_coo_csnm1 implicit none - class(psb_c_coo_sparse_mat), intent(in) :: a + class(psb_lc_coo_sparse_mat), intent(in) :: a real(psb_spk_) :: res integer(psb_lpk_) :: i,j,k,m,n, nnz, ir, jc, nc, info diff --git a/base/serial/impl/psb_d_coo_impl.F90 b/base/serial/impl/psb_d_coo_impl.F90 index 88c1ef16..d4d88027 100644 --- a/base/serial/impl/psb_d_coo_impl.F90 +++ b/base/serial/impl/psb_d_coo_impl.F90 @@ -4616,7 +4616,7 @@ function psb_ld_coo_csnm1(a) result(res) use psb_d_base_mat_mod, psb_protect_name => psb_ld_coo_csnm1 implicit none - class(psb_d_coo_sparse_mat), intent(in) :: a + class(psb_ld_coo_sparse_mat), intent(in) :: a real(psb_dpk_) :: res integer(psb_lpk_) :: i,j,k,m,n, nnz, ir, jc, nc, info diff --git a/base/serial/impl/psb_s_coo_impl.F90 b/base/serial/impl/psb_s_coo_impl.F90 index 5a0a3279..5ae9dd96 100644 --- a/base/serial/impl/psb_s_coo_impl.F90 +++ b/base/serial/impl/psb_s_coo_impl.F90 @@ -4616,7 +4616,7 @@ function psb_ls_coo_csnm1(a) result(res) use psb_s_base_mat_mod, psb_protect_name => psb_ls_coo_csnm1 implicit none - class(psb_s_coo_sparse_mat), intent(in) :: a + class(psb_ls_coo_sparse_mat), intent(in) :: a real(psb_spk_) :: res integer(psb_lpk_) :: i,j,k,m,n, nnz, ir, jc, nc, info diff --git a/base/serial/impl/psb_z_coo_impl.F90 b/base/serial/impl/psb_z_coo_impl.F90 index 7ea439d0..49e6c7fe 100644 --- a/base/serial/impl/psb_z_coo_impl.F90 +++ b/base/serial/impl/psb_z_coo_impl.F90 @@ -4616,7 +4616,7 @@ function psb_lz_coo_csnm1(a) result(res) use psb_z_base_mat_mod, psb_protect_name => psb_lz_coo_csnm1 implicit none - class(psb_z_coo_sparse_mat), intent(in) :: a + class(psb_lz_coo_sparse_mat), intent(in) :: a real(psb_dpk_) :: res integer(psb_lpk_) :: i,j,k,m,n, nnz, ir, jc, nc, info diff --git a/base/serial/sort/psb_c_hsort_impl.f90 b/base/serial/sort/psb_c_hsort_impl.f90 index 8a5fe3c7..c68ec73b 100644 --- a/base/serial/sort/psb_c_hsort_impl.f90 +++ b/base/serial/sort/psb_c_hsort_impl.f90 @@ -49,7 +49,8 @@ subroutine psb_chsort(x,ix,dir,flag) integer(psb_ipk_), optional, intent(in) :: dir, flag integer(psb_ipk_), optional, intent(inout) :: ix(:) - integer(psb_ipk_) :: dir_, flag_, n, i, l, err_act,info + integer(psb_ipk_) :: flag_, n, i, err_act,info + integer(psb_ipk_) :: dir_, l complex(psb_spk_) :: key integer(psb_ipk_) :: index @@ -159,7 +160,7 @@ end subroutine psb_chsort ! ! These are packaged so that they can be used to implement -! a heapsort, should the need arise +! a heapsort. ! ! ! Programming note: @@ -646,7 +647,8 @@ subroutine psi_c_idx_insert_heap(key,index,last,heap,idxs,dir,info) ! dir: sorting direction complex(psb_spk_), intent(in) :: key - integer(psb_ipk_), intent(in) :: index,dir + integer(psb_ipk_), intent(in) :: index + integer(psb_ipk_), intent(in) :: dir complex(psb_spk_), intent(inout) :: heap(:) integer(psb_ipk_), intent(inout) :: idxs(:) integer(psb_ipk_), intent(inout) :: last diff --git a/base/serial/sort/psb_c_msort_impl.f90 b/base/serial/sort/psb_c_msort_impl.f90 index fa87a425..751f6098 100644 --- a/base/serial/sort/psb_c_msort_impl.f90 +++ b/base/serial/sort/psb_c_msort_impl.f90 @@ -1,34 +1,34 @@ -! -! 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. -! -! + ! + ! 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. + ! + ! ! ! The merge-sort routines ! References: @@ -41,777 +41,771 @@ ! Addison-Wesley ! - subroutine psb_cmsort_u(x,nout,dir) - use psb_sort_mod, psb_protect_name => psb_cmsort_u - use psb_error_mod - implicit none - complex(psb_spk_), intent(inout) :: x(:) - integer(psb_ipk_), intent(out) :: nout - integer(psb_ipk_), optional, intent(in) :: dir +subroutine psb_cmsort_u(x,nout,dir) + use psb_sort_mod, psb_protect_name => psb_cmsort_u + use psb_error_mod + implicit none + complex(psb_spk_), intent(inout) :: x(:) + integer(psb_ipk_), intent(out) :: nout + integer(psb_ipk_), optional, intent(in) :: dir - integer(psb_ipk_) :: n, k - integer(psb_ipk_) :: err_act + integer(psb_ipk_) :: n, k + integer(psb_ipk_) :: err_act - integer(psb_ipk_) :: ierr(5) - character(len=20) :: name + integer(psb_ipk_) :: ierr(5) + character(len=20) :: name - name='psb_msort_u' - call psb_erractionsave(err_act) + name='psb_msort_u' + call psb_erractionsave(err_act) - n = size(x) + n = size(x) - call psb_msort(x,dir=dir) - nout = min(1,n) - do k=2,n - if (x(k) /= x(nout)) then - nout = nout + 1 - x(nout) = x(k) - endif - enddo + call psb_msort(x,dir=dir) + nout = min(1,n) + do k=2,n + if (x(k) /= x(nout)) then + nout = nout + 1 + x(nout) = x(k) + endif + enddo - return + return 9999 call psb_error_handler(err_act) - return - end subroutine psb_cmsort_u - - - - - - - - - subroutine psb_cmsort(x,ix,dir,flag) - use psb_sort_mod, psb_protect_name => psb_cmsort - use psb_error_mod - use psb_ip_reord_mod - implicit none - complex(psb_spk_), intent(inout) :: x(:) - integer(psb_ipk_), optional, intent(in) :: dir, flag - integer(psb_ipk_), optional, intent(inout) :: ix(:) - - integer(psb_ipk_) :: dir_, flag_, n, err_act - - integer(psb_ipk_), allocatable :: iaux(:) - integer(psb_ipk_) :: iret, info, i - integer(psb_ipk_) :: ierr(5) - character(len=20) :: name - - name='psb_cmsort' - call psb_erractionsave(err_act) - - if (present(dir)) then - dir_ = dir - else - dir_= psb_asort_up_ + return +end subroutine psb_cmsort_u + + +subroutine psb_cmsort(x,ix,dir,flag) + use psb_sort_mod, psb_protect_name => psb_cmsort + use psb_error_mod + use psb_ip_reord_mod + implicit none + complex(psb_spk_), intent(inout) :: x(:) + integer(psb_ipk_), optional, intent(in) :: dir, flag + integer(psb_ipk_), optional, intent(inout) :: ix(:) + + integer(psb_ipk_) :: dir_, flag_, n, err_act + + integer(psb_ipk_), allocatable :: iaux(:) + integer(psb_ipk_) :: iret, info, i + integer(psb_ipk_) :: ierr(5) + character(len=20) :: name + + name='psb_cmsort' + call psb_erractionsave(err_act) + + if (present(dir)) then + dir_ = dir + else + dir_= psb_asort_up_ + end if + select case(dir_) + case( psb_lsort_up_, psb_lsort_down_, psb_alsort_up_, psb_alsort_down_,& + & psb_asort_up_, psb_asort_down_) + ! OK keep going + case default + ierr(1) = 3; ierr(2) = dir_; + call psb_errpush(psb_err_input_value_invalid_i_,name,i_err=ierr) + goto 9999 + end select + + n = size(x) + + if (present(ix)) then + if (size(ix) < n) then + ierr(1) = 2; ierr(2) = size(ix); + call psb_errpush(psb_err_input_asize_invalid_i_,name,i_err=ierr) + goto 9999 + end if + if (present(flag)) then + flag_ = flag + else + flag_ = psb_sort_ovw_idx_ end if - select case(dir_) - case( psb_lsort_up_, psb_lsort_down_, psb_alsort_up_, psb_alsort_down_,& - & psb_asort_up_, psb_asort_down_) + select case(flag_) + case(psb_sort_ovw_idx_) + do i=1,n + ix(i) = i + end do + case (psb_sort_keep_idx_) ! OK keep going case default - ierr(1) = 3; ierr(2) = dir_; + ierr(1) = 4; ierr(2) = flag_; call psb_errpush(psb_err_input_value_invalid_i_,name,i_err=ierr) goto 9999 end select - - n = size(x) - - if (present(ix)) then - if (size(ix) < n) then - ierr(1) = 2; ierr(2) = size(ix); - call psb_errpush(psb_err_input_asize_invalid_i_,name,i_err=ierr) - goto 9999 - end if - if (present(flag)) then - flag_ = flag - else - flag_ = psb_sort_ovw_idx_ - end if - select case(flag_) - case(psb_sort_ovw_idx_) - do i=1,n - ix(i) = i - end do - case (psb_sort_keep_idx_) - ! OK keep going - case default - ierr(1) = 4; ierr(2) = flag_; - call psb_errpush(psb_err_input_value_invalid_i_,name,i_err=ierr) - goto 9999 - end select + end if + + allocate(iaux(0:n+1),stat=info) + if (info /= psb_success_) then + call psb_errpush(psb_err_alloc_dealloc_,r_name='psb_c_msort') + goto 9999 + endif + + select case(dir_) + case (psb_lsort_up_) + call psi_c_lmsort_up(n,x,iaux,iret) + case (psb_lsort_down_) + call psi_c_lmsort_dw(n,x,iaux,iret) + case (psb_alsort_up_) + call psi_c_almsort_up(n,x,iaux,iret) + case (psb_alsort_down_) + call psi_c_almsort_dw(n,x,iaux,iret) + case (psb_asort_up_) + call psi_c_amsort_up(n,x,iaux,iret) + case (psb_asort_down_) + call psi_c_amsort_dw(n,x,iaux,iret) + end select + ! + ! Do the actual reordering, since the inner routines + ! only provide linked pointers. + ! + if (iret == 0 ) then + if (present(ix)) then + call psb_ip_reord(n,x,ix,iaux) + else + call psb_ip_reord(n,x,iaux) end if + end if - allocate(iaux(0:n+1),stat=info) - if (info /= psb_success_) then - call psb_errpush(psb_err_alloc_dealloc_,r_name='psb_c_msort') - goto 9999 - endif + return - select case(dir_) - case (psb_lsort_up_) - call psi_c_lmsort_up(n,x,iaux,iret) - case (psb_lsort_down_) - call psi_c_lmsort_dw(n,x,iaux,iret) - case (psb_alsort_up_) - call psi_c_almsort_up(n,x,iaux,iret) - case (psb_alsort_down_) - call psi_c_almsort_dw(n,x,iaux,iret) - case (psb_asort_up_) - call psi_c_amsort_up(n,x,iaux,iret) - case (psb_asort_down_) - call psi_c_amsort_dw(n,x,iaux,iret) - end select - ! - ! Do the actual reordering, since the inner routines - ! only provide linked pointers. - ! - if (iret == 0 ) then - if (present(ix)) then - call psb_ip_reord(n,x,ix,iaux) - else - call psb_ip_reord(n,x,iaux) - end if - end if +9999 call psb_error_handler(err_act) - return + return -9999 call psb_error_handler(err_act) - return - - - end subroutine psb_cmsort - - subroutine psi_c_lmsort_up(n,k,l,iret) - use psb_const_mod - use psi_lcx_mod - implicit none - integer(psb_ipk_) :: n, iret - complex(psb_spk_) :: k(n) - integer(psb_ipk_) :: l(0:n+1) - ! - integer(psb_ipk_) :: p,q,s,t - ! .. - iret = 0 - ! first step: we are preparing ordered sublists, exploiting - ! what order was already in the input data; negative links - ! mark the end of the sublists - l(0) = 1 - t = n + 1 - do p = 1,n - 1 - if (k(p) <= k(p+1)) then - l(p) = p + 1 - else - l(t) = - (p+1) - t = p - end if - end do - l(t) = 0 - l(n) = 0 - ! see if the input was already sorted - if (l(n+1) == 0) then - iret = 1 - return +end subroutine psb_cmsort + +subroutine psi_c_lmsort_up(n,k,l,iret) + use psb_const_mod + use psi_lcx_mod + implicit none + integer(psb_ipk_) :: n, iret + complex(psb_spk_) :: k(n) + integer(psb_ipk_) :: l(0:n+1) + ! + integer(psb_ipk_) :: p,q,s,t + ! .. + iret = 0 + ! first step: we are preparing ordered sublists, exploiting + ! what order was already in the input data; negative links + ! mark the end of the sublists + l(0) = 1 + t = n + 1 + do p = 1,n - 1 + if (k(p) <= k(p+1)) then + l(p) = p + 1 else - l(n+1) = abs(l(n+1)) + l(t) = - (p+1) + t = p end if + end do + l(t) = 0 + l(n) = 0 + ! see if the input was already sorted + if (l(n+1) == 0) then + iret = 1 + return + else + l(n+1) = abs(l(n+1)) + end if + + mergepass: do + ! otherwise, begin a pass through the list. + ! throughout all the subroutine we have: + ! p, q: pointing to the sublists being merged + ! s: pointing to the most recently processed record + ! t: pointing to the end of previously completed sublist + s = 0 + t = n + 1 + p = l(s) + q = l(t) + if (q == 0) exit mergepass - mergepass: do - ! otherwise, begin a pass through the list. - ! throughout all the subroutine we have: - ! p, q: pointing to the sublists being merged - ! s: pointing to the most recently processed record - ! t: pointing to the end of previously completed sublist - s = 0 - t = n + 1 - p = l(s) - q = l(t) - if (q == 0) exit mergepass - - outer: do + outer: do - if (k(p) > k(q)) then + if (k(p) > k(q)) then - l(s) = sign(q,l(s)) - s = q - q = l(q) - if (q > 0) then - do - if (k(p) <= k(q)) cycle outer - s = q - q = l(q) - if (q <= 0) exit - end do - end if - l(s) = p - s = t + l(s) = sign(q,l(s)) + s = q + q = l(q) + if (q > 0) then do - t = p - p = l(p) - if (p <= 0) exit + if (k(p) <= k(q)) cycle outer + s = q + q = l(q) + if (q <= 0) exit end do + end if + l(s) = p + s = t + do + t = p + p = l(p) + if (p <= 0) exit + end do - else + else - l(s) = sign(p,l(s)) - s = p - p = l(p) - if (p>0) then - do - if (k(p) > k(q)) cycle outer - s = p - p = l(p) - if (p <= 0) exit - end do - end if - ! otherwise, one sublist ended, and we append to it the rest - ! of the other one. - l(s) = q - s = t + l(s) = sign(p,l(s)) + s = p + p = l(p) + if (p>0) then do - t = q - q = l(q) - if (q <= 0) exit + if (k(p) > k(q)) cycle outer + s = p + p = l(p) + if (p <= 0) exit end do end if + ! otherwise, one sublist ended, and we append to it the rest + ! of the other one. + l(s) = q + s = t + do + t = q + q = l(q) + if (q <= 0) exit + end do + end if - p = -p - q = -q - if (q == 0) then - l(s) = sign(p,l(s)) - l(t) = 0 - exit outer - end if - end do outer - end do mergepass - - end subroutine psi_c_lmsort_up - - subroutine psi_c_lmsort_dw(n,k,l,iret) - use psb_const_mod - use psi_lcx_mod - implicit none - integer(psb_ipk_) :: n, iret - complex(psb_spk_) :: k(n) - integer(psb_ipk_) :: l(0:n+1) - ! - integer(psb_ipk_) :: p,q,s,t - ! .. - iret = 0 - ! first step: we are preparing ordered sublists, exploiting - ! what order was already in the input data; negative links - ! mark the end of the sublists - l(0) = 1 - t = n + 1 - do p = 1,n - 1 - if (k(p) >= k(p+1)) then - l(p) = p + 1 - else - l(t) = - (p+1) - t = p + p = -p + q = -q + if (q == 0) then + l(s) = sign(p,l(s)) + l(t) = 0 + exit outer end if - end do - l(t) = 0 - l(n) = 0 - ! see if the input was already sorted - if (l(n+1) == 0) then - iret = 1 - return + end do outer + end do mergepass + +end subroutine psi_c_lmsort_up + +subroutine psi_c_lmsort_dw(n,k,l,iret) + use psb_const_mod + use psi_lcx_mod + implicit none + integer(psb_ipk_) :: n, iret + complex(psb_spk_) :: k(n) + integer(psb_ipk_) :: l(0:n+1) + ! + integer(psb_ipk_) :: p,q,s,t + ! .. + iret = 0 + ! first step: we are preparing ordered sublists, exploiting + ! what order was already in the input data; negative links + ! mark the end of the sublists + l(0) = 1 + t = n + 1 + do p = 1,n - 1 + if (k(p) >= k(p+1)) then + l(p) = p + 1 else - l(n+1) = abs(l(n+1)) + l(t) = - (p+1) + t = p end if + end do + l(t) = 0 + l(n) = 0 + ! see if the input was already sorted + if (l(n+1) == 0) then + iret = 1 + return + else + l(n+1) = abs(l(n+1)) + end if + + mergepass: do + ! otherwise, begin a pass through the list. + ! throughout all the subroutine we have: + ! p, q: pointing to the sublists being merged + ! s: pointing to the most recently processed record + ! t: pointing to the end of previously completed sublist + s = 0 + t = n + 1 + p = l(s) + q = l(t) + if (q == 0) exit mergepass - mergepass: do - ! otherwise, begin a pass through the list. - ! throughout all the subroutine we have: - ! p, q: pointing to the sublists being merged - ! s: pointing to the most recently processed record - ! t: pointing to the end of previously completed sublist - s = 0 - t = n + 1 - p = l(s) - q = l(t) - if (q == 0) exit mergepass - - outer: do + outer: do - if (k(p) < k(q)) then + if (k(p) < k(q)) then - l(s) = sign(q,l(s)) - s = q - q = l(q) - if (q > 0) then - do - if (k(p) >= k(q)) cycle outer - s = q - q = l(q) - if (q <= 0) exit - end do - end if - l(s) = p - s = t + l(s) = sign(q,l(s)) + s = q + q = l(q) + if (q > 0) then do - t = p - p = l(p) - if (p <= 0) exit + if (k(p) >= k(q)) cycle outer + s = q + q = l(q) + if (q <= 0) exit end do + end if + l(s) = p + s = t + do + t = p + p = l(p) + if (p <= 0) exit + end do - else + else - l(s) = sign(p,l(s)) - s = p - p = l(p) - if (p>0) then - do - if (k(p) < k(q)) cycle outer - s = p - p = l(p) - if (p <= 0) exit - end do - end if - ! otherwise, one sublist ended, and we append to it the rest - ! of the other one. - l(s) = q - s = t + l(s) = sign(p,l(s)) + s = p + p = l(p) + if (p>0) then do - t = q - q = l(q) - if (q <= 0) exit + if (k(p) < k(q)) cycle outer + s = p + p = l(p) + if (p <= 0) exit end do end if + ! otherwise, one sublist ended, and we append to it the rest + ! of the other one. + l(s) = q + s = t + do + t = q + q = l(q) + if (q <= 0) exit + end do + end if - p = -p - q = -q - if (q == 0) then - l(s) = sign(p,l(s)) - l(t) = 0 - exit outer - end if - end do outer - end do mergepass - - end subroutine psi_c_lmsort_dw - - subroutine psi_c_amsort_up(n,k,l,iret) - use psb_const_mod - use psi_acx_mod - implicit none - integer(psb_ipk_) :: n, iret - complex(psb_spk_) :: k(n) - integer(psb_ipk_) :: l(0:n+1) - ! - integer(psb_ipk_) :: p,q,s,t - ! .. - iret = 0 - ! first step: we are preparing ordered sublists, exploiting - ! what order was already in the input data; negative links - ! mark the end of the sublists - l(0) = 1 - t = n + 1 - do p = 1,n - 1 - if (k(p) <= k(p+1)) then - l(p) = p + 1 - else - l(t) = - (p+1) - t = p + p = -p + q = -q + if (q == 0) then + l(s) = sign(p,l(s)) + l(t) = 0 + exit outer end if - end do - l(t) = 0 - l(n) = 0 - ! see if the input was already sorted - if (l(n+1) == 0) then - iret = 1 - return + end do outer + end do mergepass + +end subroutine psi_c_lmsort_dw + +subroutine psi_c_amsort_up(n,k,l,iret) + use psb_const_mod + use psi_acx_mod + implicit none + integer(psb_ipk_) :: n, iret + complex(psb_spk_) :: k(n) + integer(psb_ipk_) :: l(0:n+1) + ! + integer(psb_ipk_) :: p,q,s,t + ! .. + iret = 0 + ! first step: we are preparing ordered sublists, exploiting + ! what order was already in the input data; negative links + ! mark the end of the sublists + l(0) = 1 + t = n + 1 + do p = 1,n - 1 + if (k(p) <= k(p+1)) then + l(p) = p + 1 else - l(n+1) = abs(l(n+1)) + l(t) = - (p+1) + t = p end if + end do + l(t) = 0 + l(n) = 0 + ! see if the input was already sorted + if (l(n+1) == 0) then + iret = 1 + return + else + l(n+1) = abs(l(n+1)) + end if + + mergepass: do + ! otherwise, begin a pass through the list. + ! throughout all the subroutine we have: + ! p, q: pointing to the sublists being merged + ! s: pointing to the most recently processed record + ! t: pointing to the end of previously completed sublist + s = 0 + t = n + 1 + p = l(s) + q = l(t) + if (q == 0) exit mergepass - mergepass: do - ! otherwise, begin a pass through the list. - ! throughout all the subroutine we have: - ! p, q: pointing to the sublists being merged - ! s: pointing to the most recently processed record - ! t: pointing to the end of previously completed sublist - s = 0 - t = n + 1 - p = l(s) - q = l(t) - if (q == 0) exit mergepass - - outer: do + outer: do - if (k(p) > k(q)) then + if (k(p) > k(q)) then - l(s) = sign(q,l(s)) - s = q - q = l(q) - if (q > 0) then - do - if (k(p) <= k(q)) cycle outer - s = q - q = l(q) - if (q <= 0) exit - end do - end if - l(s) = p - s = t + l(s) = sign(q,l(s)) + s = q + q = l(q) + if (q > 0) then do - t = p - p = l(p) - if (p <= 0) exit + if (k(p) <= k(q)) cycle outer + s = q + q = l(q) + if (q <= 0) exit end do + end if + l(s) = p + s = t + do + t = p + p = l(p) + if (p <= 0) exit + end do - else + else - l(s) = sign(p,l(s)) - s = p - p = l(p) - if (p>0) then - do - if (k(p) > k(q)) cycle outer - s = p - p = l(p) - if (p <= 0) exit - end do - end if - ! otherwise, one sublist ended, and we append to it the rest - ! of the other one. - l(s) = q - s = t + l(s) = sign(p,l(s)) + s = p + p = l(p) + if (p>0) then do - t = q - q = l(q) - if (q <= 0) exit + if (k(p) > k(q)) cycle outer + s = p + p = l(p) + if (p <= 0) exit end do end if + ! otherwise, one sublist ended, and we append to it the rest + ! of the other one. + l(s) = q + s = t + do + t = q + q = l(q) + if (q <= 0) exit + end do + end if - p = -p - q = -q - if (q == 0) then - l(s) = sign(p,l(s)) - l(t) = 0 - exit outer - end if - end do outer - end do mergepass - - end subroutine psi_c_amsort_up - - subroutine psi_c_amsort_dw(n,k,l,iret) - use psb_const_mod - use psi_acx_mod - implicit none - integer(psb_ipk_) :: n, iret - complex(psb_spk_) :: k(n) - integer(psb_ipk_) :: l(0:n+1) - ! - integer(psb_ipk_) :: p,q,s,t - ! .. - iret = 0 - ! first step: we are preparing ordered sublists, exploiting - ! what order was already in the input data; negative links - ! mark the end of the sublists - l(0) = 1 - t = n + 1 - do p = 1,n - 1 - if (k(p) >= k(p+1)) then - l(p) = p + 1 - else - l(t) = - (p+1) - t = p + p = -p + q = -q + if (q == 0) then + l(s) = sign(p,l(s)) + l(t) = 0 + exit outer end if - end do - l(t) = 0 - l(n) = 0 - ! see if the input was already sorted - if (l(n+1) == 0) then - iret = 1 - return + end do outer + end do mergepass + +end subroutine psi_c_amsort_up + +subroutine psi_c_amsort_dw(n,k,l,iret) + use psb_const_mod + use psi_acx_mod + implicit none + integer(psb_ipk_) :: n, iret + complex(psb_spk_) :: k(n) + integer(psb_ipk_) :: l(0:n+1) + ! + integer(psb_ipk_) :: p,q,s,t + ! .. + iret = 0 + ! first step: we are preparing ordered sublists, exploiting + ! what order was already in the input data; negative links + ! mark the end of the sublists + l(0) = 1 + t = n + 1 + do p = 1,n - 1 + if (k(p) >= k(p+1)) then + l(p) = p + 1 else - l(n+1) = abs(l(n+1)) + l(t) = - (p+1) + t = p end if + end do + l(t) = 0 + l(n) = 0 + ! see if the input was already sorted + if (l(n+1) == 0) then + iret = 1 + return + else + l(n+1) = abs(l(n+1)) + end if + + mergepass: do + ! otherwise, begin a pass through the list. + ! throughout all the subroutine we have: + ! p, q: pointing to the sublists being merged + ! s: pointing to the most recently processed record + ! t: pointing to the end of previously completed sublist + s = 0 + t = n + 1 + p = l(s) + q = l(t) + if (q == 0) exit mergepass - mergepass: do - ! otherwise, begin a pass through the list. - ! throughout all the subroutine we have: - ! p, q: pointing to the sublists being merged - ! s: pointing to the most recently processed record - ! t: pointing to the end of previously completed sublist - s = 0 - t = n + 1 - p = l(s) - q = l(t) - if (q == 0) exit mergepass - - outer: do + outer: do - if (k(p) < k(q)) then + if (k(p) < k(q)) then - l(s) = sign(q,l(s)) - s = q - q = l(q) - if (q > 0) then - do - if (k(p) >= k(q)) cycle outer - s = q - q = l(q) - if (q <= 0) exit - end do - end if - l(s) = p - s = t + l(s) = sign(q,l(s)) + s = q + q = l(q) + if (q > 0) then do - t = p - p = l(p) - if (p <= 0) exit + if (k(p) >= k(q)) cycle outer + s = q + q = l(q) + if (q <= 0) exit end do + end if + l(s) = p + s = t + do + t = p + p = l(p) + if (p <= 0) exit + end do - else + else - l(s) = sign(p,l(s)) - s = p - p = l(p) - if (p>0) then - do - if (k(p) < k(q)) cycle outer - s = p - p = l(p) - if (p <= 0) exit - end do - end if - ! otherwise, one sublist ended, and we append to it the rest - ! of the other one. - l(s) = q - s = t + l(s) = sign(p,l(s)) + s = p + p = l(p) + if (p>0) then do - t = q - q = l(q) - if (q <= 0) exit + if (k(p) < k(q)) cycle outer + s = p + p = l(p) + if (p <= 0) exit end do end if + ! otherwise, one sublist ended, and we append to it the rest + ! of the other one. + l(s) = q + s = t + do + t = q + q = l(q) + if (q <= 0) exit + end do + end if - p = -p - q = -q - if (q == 0) then - l(s) = sign(p,l(s)) - l(t) = 0 - exit outer - end if - end do outer - end do mergepass - - end subroutine psi_c_amsort_dw - - subroutine psi_c_almsort_up(n,k,l,iret) - use psb_const_mod - use psi_alcx_mod - implicit none - integer(psb_ipk_) :: n, iret - complex(psb_spk_) :: k(n) - integer(psb_ipk_) :: l(0:n+1) - ! - integer(psb_ipk_) :: p,q,s,t - ! .. - iret = 0 - ! first step: we are preparing ordered sublists, exploiting - ! what order was already in the input data; negative links - ! mark the end of the sublists - l(0) = 1 - t = n + 1 - do p = 1,n - 1 - if (k(p) <= k(p+1)) then - l(p) = p + 1 - else - l(t) = - (p+1) - t = p + p = -p + q = -q + if (q == 0) then + l(s) = sign(p,l(s)) + l(t) = 0 + exit outer end if - end do - l(t) = 0 - l(n) = 0 - ! see if the input was already sorted - if (l(n+1) == 0) then - iret = 1 - return + end do outer + end do mergepass + +end subroutine psi_c_amsort_dw + +subroutine psi_c_almsort_up(n,k,l,iret) + use psb_const_mod + use psi_alcx_mod + implicit none + integer(psb_ipk_) :: n, iret + complex(psb_spk_) :: k(n) + integer(psb_ipk_) :: l(0:n+1) + ! + integer(psb_ipk_) :: p,q,s,t + ! .. + iret = 0 + ! first step: we are preparing ordered sublists, exploiting + ! what order was already in the input data; negative links + ! mark the end of the sublists + l(0) = 1 + t = n + 1 + do p = 1,n - 1 + if (k(p) <= k(p+1)) then + l(p) = p + 1 else - l(n+1) = abs(l(n+1)) + l(t) = - (p+1) + t = p end if + end do + l(t) = 0 + l(n) = 0 + ! see if the input was already sorted + if (l(n+1) == 0) then + iret = 1 + return + else + l(n+1) = abs(l(n+1)) + end if + + mergepass: do + ! otherwise, begin a pass through the list. + ! throughout all the subroutine we have: + ! p, q: pointing to the sublists being merged + ! s: pointing to the most recently processed record + ! t: pointing to the end of previously completed sublist + s = 0 + t = n + 1 + p = l(s) + q = l(t) + if (q == 0) exit mergepass - mergepass: do - ! otherwise, begin a pass through the list. - ! throughout all the subroutine we have: - ! p, q: pointing to the sublists being merged - ! s: pointing to the most recently processed record - ! t: pointing to the end of previously completed sublist - s = 0 - t = n + 1 - p = l(s) - q = l(t) - if (q == 0) exit mergepass - - outer: do + outer: do - if (k(p) > k(q)) then + if (k(p) > k(q)) then - l(s) = sign(q,l(s)) - s = q - q = l(q) - if (q > 0) then - do - if (k(p) <= k(q)) cycle outer - s = q - q = l(q) - if (q <= 0) exit - end do - end if - l(s) = p - s = t + l(s) = sign(q,l(s)) + s = q + q = l(q) + if (q > 0) then do - t = p - p = l(p) - if (p <= 0) exit + if (k(p) <= k(q)) cycle outer + s = q + q = l(q) + if (q <= 0) exit end do + end if + l(s) = p + s = t + do + t = p + p = l(p) + if (p <= 0) exit + end do - else + else - l(s) = sign(p,l(s)) - s = p - p = l(p) - if (p>0) then - do - if (k(p) > k(q)) cycle outer - s = p - p = l(p) - if (p <= 0) exit - end do - end if - ! otherwise, one sublist ended, and we append to it the rest - ! of the other one. - l(s) = q - s = t + l(s) = sign(p,l(s)) + s = p + p = l(p) + if (p>0) then do - t = q - q = l(q) - if (q <= 0) exit + if (k(p) > k(q)) cycle outer + s = p + p = l(p) + if (p <= 0) exit end do end if + ! otherwise, one sublist ended, and we append to it the rest + ! of the other one. + l(s) = q + s = t + do + t = q + q = l(q) + if (q <= 0) exit + end do + end if - p = -p - q = -q - if (q == 0) then - l(s) = sign(p,l(s)) - l(t) = 0 - exit outer - end if - end do outer - end do mergepass - - end subroutine psi_c_almsort_up - - subroutine psi_c_almsort_dw(n,k,l,iret) - use psb_const_mod - use psi_alcx_mod - implicit none - integer(psb_ipk_) :: n, iret - complex(psb_spk_) :: k(n) - integer(psb_ipk_) :: l(0:n+1) - ! - integer(psb_ipk_) :: p,q,s,t - ! .. - iret = 0 - ! first step: we are preparing ordered sublists, exploiting - ! what order was already in the input data; negative links - ! mark the end of the sublists - l(0) = 1 - t = n + 1 - do p = 1,n - 1 - if (k(p) >= k(p+1)) then - l(p) = p + 1 - else - l(t) = - (p+1) - t = p + p = -p + q = -q + if (q == 0) then + l(s) = sign(p,l(s)) + l(t) = 0 + exit outer end if - end do - l(t) = 0 - l(n) = 0 - ! see if the input was already sorted - if (l(n+1) == 0) then - iret = 1 - return + end do outer + end do mergepass + +end subroutine psi_c_almsort_up + +subroutine psi_c_almsort_dw(n,k,l,iret) + use psb_const_mod + use psi_alcx_mod + implicit none + integer(psb_ipk_) :: n, iret + complex(psb_spk_) :: k(n) + integer(psb_ipk_) :: l(0:n+1) + ! + integer(psb_ipk_) :: p,q,s,t + ! .. + iret = 0 + ! first step: we are preparing ordered sublists, exploiting + ! what order was already in the input data; negative links + ! mark the end of the sublists + l(0) = 1 + t = n + 1 + do p = 1,n - 1 + if (k(p) >= k(p+1)) then + l(p) = p + 1 else - l(n+1) = abs(l(n+1)) + l(t) = - (p+1) + t = p end if + end do + l(t) = 0 + l(n) = 0 + ! see if the input was already sorted + if (l(n+1) == 0) then + iret = 1 + return + else + l(n+1) = abs(l(n+1)) + end if + + mergepass: do + ! otherwise, begin a pass through the list. + ! throughout all the subroutine we have: + ! p, q: pointing to the sublists being merged + ! s: pointing to the most recently processed record + ! t: pointing to the end of previously completed sublist + s = 0 + t = n + 1 + p = l(s) + q = l(t) + if (q == 0) exit mergepass - mergepass: do - ! otherwise, begin a pass through the list. - ! throughout all the subroutine we have: - ! p, q: pointing to the sublists being merged - ! s: pointing to the most recently processed record - ! t: pointing to the end of previously completed sublist - s = 0 - t = n + 1 - p = l(s) - q = l(t) - if (q == 0) exit mergepass - - outer: do + outer: do - if (k(p) < k(q)) then + if (k(p) < k(q)) then - l(s) = sign(q,l(s)) - s = q - q = l(q) - if (q > 0) then - do - if (k(p) >= k(q)) cycle outer - s = q - q = l(q) - if (q <= 0) exit - end do - end if - l(s) = p - s = t + l(s) = sign(q,l(s)) + s = q + q = l(q) + if (q > 0) then do - t = p - p = l(p) - if (p <= 0) exit + if (k(p) >= k(q)) cycle outer + s = q + q = l(q) + if (q <= 0) exit end do + end if + l(s) = p + s = t + do + t = p + p = l(p) + if (p <= 0) exit + end do - else + else - l(s) = sign(p,l(s)) - s = p - p = l(p) - if (p>0) then - do - if (k(p) < k(q)) cycle outer - s = p - p = l(p) - if (p <= 0) exit - end do - end if - ! otherwise, one sublist ended, and we append to it the rest - ! of the other one. - l(s) = q - s = t + l(s) = sign(p,l(s)) + s = p + p = l(p) + if (p>0) then do - t = q - q = l(q) - if (q <= 0) exit + if (k(p) < k(q)) cycle outer + s = p + p = l(p) + if (p <= 0) exit end do end if + ! otherwise, one sublist ended, and we append to it the rest + ! of the other one. + l(s) = q + s = t + do + t = q + q = l(q) + if (q <= 0) exit + end do + end if - p = -p - q = -q - if (q == 0) then - l(s) = sign(p,l(s)) - l(t) = 0 - exit outer - end if - end do outer - end do mergepass + p = -p + q = -q + if (q == 0) then + l(s) = sign(p,l(s)) + l(t) = 0 + exit outer + end if + end do outer + end do mergepass - end subroutine psi_c_almsort_dw +end subroutine psi_c_almsort_dw diff --git a/base/serial/sort/psb_d_hsort_impl.f90 b/base/serial/sort/psb_d_hsort_impl.f90 index b83e93cc..ffb952d5 100644 --- a/base/serial/sort/psb_d_hsort_impl.f90 +++ b/base/serial/sort/psb_d_hsort_impl.f90 @@ -49,7 +49,8 @@ subroutine psb_dhsort(x,ix,dir,flag) integer(psb_ipk_), optional, intent(in) :: dir, flag integer(psb_ipk_), optional, intent(inout) :: ix(:) - integer(psb_ipk_) :: dir_, flag_, n, i, l, err_act,info + integer(psb_ipk_) :: flag_, n, i, err_act,info + integer(psb_ipk_) :: dir_, l real(psb_dpk_) :: key integer(psb_ipk_) :: index @@ -159,7 +160,7 @@ end subroutine psb_dhsort ! ! These are packaged so that they can be used to implement -! a heapsort, should the need arise +! a heapsort. ! ! ! Programming note: @@ -540,11 +541,12 @@ subroutine psi_d_idx_heap_get_first(key,index,last,heap,idxs,dir,info) use psb_sort_mod, psb_protect_name => psi_d_idx_heap_get_first implicit none + real(psb_dpk_), intent(inout) :: key real(psb_dpk_), intent(inout) :: heap(:) - integer(psb_ipk_), intent(out) :: index,info + integer(psb_ipk_), intent(out) :: index + integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(inout) :: last,idxs(:) integer(psb_ipk_), intent(in) :: dir - real(psb_dpk_), intent(out) :: key integer(psb_ipk_) :: i, j,itemp real(psb_dpk_) :: temp diff --git a/base/serial/sort/psb_d_msort_impl.f90 b/base/serial/sort/psb_d_msort_impl.f90 index 01491b0b..11029818 100644 --- a/base/serial/sort/psb_d_msort_impl.f90 +++ b/base/serial/sort/psb_d_msort_impl.f90 @@ -1,34 +1,34 @@ -! -! 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. -! -! + ! + ! 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. + ! + ! ! ! The merge-sort routines ! References: @@ -41,618 +41,612 @@ ! Addison-Wesley ! - subroutine psb_dmsort_u(x,nout,dir) - use psb_sort_mod, psb_protect_name => psb_dmsort_u - use psb_error_mod - implicit none - real(psb_dpk_), intent(inout) :: x(:) - integer(psb_ipk_), intent(out) :: nout - integer(psb_ipk_), optional, intent(in) :: dir +subroutine psb_dmsort_u(x,nout,dir) + use psb_sort_mod, psb_protect_name => psb_dmsort_u + use psb_error_mod + implicit none + real(psb_dpk_), intent(inout) :: x(:) + integer(psb_ipk_), intent(out) :: nout + integer(psb_ipk_), optional, intent(in) :: dir - integer(psb_ipk_) :: n, k - integer(psb_ipk_) :: err_act + integer(psb_ipk_) :: n, k + integer(psb_ipk_) :: err_act - integer(psb_ipk_) :: ierr(5) - character(len=20) :: name + integer(psb_ipk_) :: ierr(5) + character(len=20) :: name - name='psb_msort_u' - call psb_erractionsave(err_act) + name='psb_msort_u' + call psb_erractionsave(err_act) - n = size(x) + n = size(x) - call psb_msort(x,dir=dir) - nout = min(1,n) - do k=2,n - if (x(k) /= x(nout)) then - nout = nout + 1 - x(nout) = x(k) - endif - enddo + call psb_msort(x,dir=dir) + nout = min(1,n) + do k=2,n + if (x(k) /= x(nout)) then + nout = nout + 1 + x(nout) = x(k) + endif + enddo - return + return 9999 call psb_error_handler(err_act) - return - end subroutine psb_dmsort_u + return +end subroutine psb_dmsort_u - function psb_dbsrch(key,n,v) result(ipos) - use psb_sort_mod, psb_protect_name => psb_dbsrch - implicit none - integer(psb_ipk_) :: ipos, n - real(psb_dpk_) :: key - real(psb_dpk_) :: v(:) +function psb_dbsrch(key,n,v) result(ipos) + use psb_sort_mod, psb_protect_name => psb_dbsrch + implicit none + integer(psb_ipk_) :: ipos, n + real(psb_dpk_) :: key + real(psb_dpk_) :: v(:) - integer(psb_ipk_) :: lb, ub, m, i + integer(psb_ipk_) :: lb, ub, m, i - ipos = -1 - if (n<5) then - do i=1,n - if (key.eq.v(i)) then - ipos = i - return - end if - enddo - return - end if - - lb = 1 - ub = n - - do while (lb.le.ub) - m = (lb+ub)/2 - if (key.eq.v(m)) then - ipos = m - lb = ub + 1 - else if (key < v(m)) then - ub = m-1 - else - lb = m + 1 - end if - enddo - return - end function psb_dbsrch - - function psb_dssrch(key,n,v) result(ipos) - use psb_sort_mod, psb_protect_name => psb_dssrch - implicit none - integer(psb_ipk_) :: ipos, n - real(psb_dpk_) :: key - real(psb_dpk_) :: v(:) - - integer(psb_ipk_) :: i - - ipos = -1 + ipos = -1 + if (n<5) then do i=1,n if (key.eq.v(i)) then ipos = i return end if enddo - return - end function psb_dssrch - - subroutine psb_dmsort(x,ix,dir,flag) - use psb_sort_mod, psb_protect_name => psb_dmsort - use psb_error_mod - use psb_ip_reord_mod - implicit none - real(psb_dpk_), intent(inout) :: x(:) - integer(psb_ipk_), optional, intent(in) :: dir, flag - integer(psb_ipk_), optional, intent(inout) :: ix(:) - - integer(psb_ipk_) :: dir_, flag_, n, err_act - - integer(psb_ipk_), allocatable :: iaux(:) - integer(psb_ipk_) :: iret, info, i - integer(psb_ipk_) :: ierr(5) - character(len=20) :: name - - name='psb_dmsort' - call psb_erractionsave(err_act) - - if (present(dir)) then - dir_ = dir - else - dir_= psb_sort_up_ + end if + + lb = 1 + ub = n + + do while (lb.le.ub) + m = (lb+ub)/2 + if (key.eq.v(m)) then + ipos = m + lb = ub + 1 + else if (key < v(m)) then + ub = m-1 + else + lb = m + 1 end if - select case(dir_) - case( psb_sort_up_, psb_sort_down_, psb_asort_up_, psb_asort_down_) + enddo + return +end function psb_dbsrch + +function psb_dssrch(key,n,v) result(ipos) + use psb_sort_mod, psb_protect_name => psb_dssrch + implicit none + integer(psb_ipk_) :: ipos, n + real(psb_dpk_) :: key + real(psb_dpk_) :: v(:) + + integer(psb_ipk_) :: i + + ipos = -1 + do i=1,n + if (key.eq.v(i)) then + ipos = i + return + end if + enddo + + return +end function psb_dssrch + +subroutine psb_dmsort(x,ix,dir,flag) + use psb_sort_mod, psb_protect_name => psb_dmsort + use psb_error_mod + use psb_ip_reord_mod + implicit none + real(psb_dpk_), intent(inout) :: x(:) + integer(psb_ipk_), optional, intent(in) :: dir, flag + integer(psb_ipk_), optional, intent(inout) :: ix(:) + + integer(psb_ipk_) :: dir_, flag_, n, err_act + + integer(psb_ipk_), allocatable :: iaux(:) + integer(psb_ipk_) :: iret, info, i + integer(psb_ipk_) :: ierr(5) + character(len=20) :: name + + name='psb_dmsort' + call psb_erractionsave(err_act) + + if (present(dir)) then + dir_ = dir + else + dir_= psb_sort_up_ + end if + select case(dir_) + case( psb_sort_up_, psb_sort_down_, psb_asort_up_, psb_asort_down_) + ! OK keep going + case default + ierr(1) = 3; ierr(2) = dir_; + call psb_errpush(psb_err_input_value_invalid_i_,name,i_err=ierr) + goto 9999 + end select + + n = size(x) + + if (present(ix)) then + if (size(ix) < n) then + ierr(1) = 2; ierr(2) = size(ix); + call psb_errpush(psb_err_input_asize_invalid_i_,name,i_err=ierr) + goto 9999 + end if + if (present(flag)) then + flag_ = flag + else + flag_ = psb_sort_ovw_idx_ + end if + select case(flag_) + case(psb_sort_ovw_idx_) + do i=1,n + ix(i) = i + end do + case (psb_sort_keep_idx_) ! OK keep going case default - ierr(1) = 3; ierr(2) = dir_; + ierr(1) = 4; ierr(2) = flag_; call psb_errpush(psb_err_input_value_invalid_i_,name,i_err=ierr) goto 9999 end select - n = size(x) - - if (present(ix)) then - if (size(ix) < n) then - ierr(1) = 2; ierr(2) = size(ix); - call psb_errpush(psb_err_input_asize_invalid_i_,name,i_err=ierr) - goto 9999 - end if - if (present(flag)) then - flag_ = flag - else - flag_ = psb_sort_ovw_idx_ - end if - select case(flag_) - case(psb_sort_ovw_idx_) - do i=1,n - ix(i) = i - end do - case (psb_sort_keep_idx_) - ! OK keep going - case default - ierr(1) = 4; ierr(2) = flag_; - call psb_errpush(psb_err_input_value_invalid_i_,name,i_err=ierr) - goto 9999 - end select - - end if - - allocate(iaux(0:n+1),stat=info) - if (info /= psb_success_) then - call psb_errpush(psb_err_alloc_dealloc_,r_name='psb_d_msort') - goto 9999 - endif - - select case(dir_) - case (psb_sort_up_) - call psi_d_msort_up(n,x,iaux,iret) - case (psb_sort_down_) - call psi_d_msort_dw(n,x,iaux,iret) - case (psb_asort_up_) - call psi_d_amsort_up(n,x,iaux,iret) - case (psb_asort_down_) - call psi_d_amsort_dw(n,x,iaux,iret) - end select - ! - ! Do the actual reordering, since the inner routines - ! only provide linked pointers. - ! - if (iret == 0 ) then - if (present(ix)) then - call psb_ip_reord(n,x,ix,iaux) - else - call psb_ip_reord(n,x,iaux) - end if + end if + + allocate(iaux(0:n+1),stat=info) + if (info /= psb_success_) then + call psb_errpush(psb_err_alloc_dealloc_,r_name='psb_d_msort') + goto 9999 + endif + + select case(dir_) + case (psb_sort_up_) + call psi_d_msort_up(n,x,iaux,iret) + case (psb_sort_down_) + call psi_d_msort_dw(n,x,iaux,iret) + case (psb_asort_up_) + call psi_d_amsort_up(n,x,iaux,iret) + case (psb_asort_down_) + call psi_d_amsort_dw(n,x,iaux,iret) + end select + ! + ! Do the actual reordering, since the inner routines + ! only provide linked pointers. + ! + if (iret == 0 ) then + if (present(ix)) then + call psb_ip_reord(n,x,ix,iaux) + else + call psb_ip_reord(n,x,iaux) end if + end if - return + return 9999 call psb_error_handler(err_act) - return + return - end subroutine psb_dmsort - - subroutine psi_d_msort_up(n,k,l,iret) - use psb_const_mod - implicit none - integer(psb_ipk_) :: n, iret - real(psb_dpk_) :: k(n) - integer(psb_ipk_) :: l(0:n+1) - ! - integer(psb_ipk_) :: p,q,s,t - ! .. - iret = 0 - ! first step: we are preparing ordered sublists, exploiting - ! what order was already in the input data; negative links - ! mark the end of the sublists - l(0) = 1 - t = n + 1 - do p = 1,n - 1 - if (k(p) <= k(p+1)) then - l(p) = p + 1 - else - l(t) = - (p+1) - t = p - end if - end do - l(t) = 0 - l(n) = 0 - ! see if the input was already sorted - if (l(n+1) == 0) then - iret = 1 - return +end subroutine psb_dmsort + +subroutine psi_d_msort_up(n,k,l,iret) + use psb_const_mod + implicit none + integer(psb_ipk_) :: n, iret + real(psb_dpk_) :: k(n) + integer(psb_ipk_) :: l(0:n+1) + ! + integer(psb_ipk_) :: p,q,s,t + ! .. + iret = 0 + ! first step: we are preparing ordered sublists, exploiting + ! what order was already in the input data; negative links + ! mark the end of the sublists + l(0) = 1 + t = n + 1 + do p = 1,n - 1 + if (k(p) <= k(p+1)) then + l(p) = p + 1 else - l(n+1) = abs(l(n+1)) + l(t) = - (p+1) + t = p end if + end do + l(t) = 0 + l(n) = 0 + ! see if the input was already sorted + if (l(n+1) == 0) then + iret = 1 + return + else + l(n+1) = abs(l(n+1)) + end if + + mergepass: do + ! otherwise, begin a pass through the list. + ! throughout all the subroutine we have: + ! p, q: pointing to the sublists being merged + ! s: pointing to the most recently processed record + ! t: pointing to the end of previously completed sublist + s = 0 + t = n + 1 + p = l(s) + q = l(t) + if (q == 0) exit mergepass - mergepass: do - ! otherwise, begin a pass through the list. - ! throughout all the subroutine we have: - ! p, q: pointing to the sublists being merged - ! s: pointing to the most recently processed record - ! t: pointing to the end of previously completed sublist - s = 0 - t = n + 1 - p = l(s) - q = l(t) - if (q == 0) exit mergepass - - outer: do - - if (k(p) > k(q)) then - - l(s) = sign(q,l(s)) - s = q - q = l(q) - if (q > 0) then - do - if (k(p) <= k(q)) cycle outer - s = q - q = l(q) - if (q <= 0) exit - end do - end if - l(s) = p - s = t - do - t = p - p = l(p) - if (p <= 0) exit - end do + outer: do - else + if (k(p) > k(q)) then - l(s) = sign(p,l(s)) - s = p - p = l(p) - if (p>0) then - do - if (k(p) > k(q)) cycle outer - s = p - p = l(p) - if (p <= 0) exit - end do - end if - ! otherwise, one sublist ended, and we append to it the rest - ! of the other one. - l(s) = q - s = t + l(s) = sign(q,l(s)) + s = q + q = l(q) + if (q > 0) then do - t = q + if (k(p) <= k(q)) cycle outer + s = q q = l(q) if (q <= 0) exit end do end if + l(s) = p + s = t + do + t = p + p = l(p) + if (p <= 0) exit + end do - p = -p - q = -q - if (q == 0) then - l(s) = sign(p,l(s)) - l(t) = 0 - exit outer + else + + l(s) = sign(p,l(s)) + s = p + p = l(p) + if (p>0) then + do + if (k(p) > k(q)) cycle outer + s = p + p = l(p) + if (p <= 0) exit + end do end if - end do outer - end do mergepass - - end subroutine psi_d_msort_up - - subroutine psi_d_msort_dw(n,k,l,iret) - use psb_const_mod - implicit none - integer(psb_ipk_) :: n, iret - real(psb_dpk_) :: k(n) - integer(psb_ipk_) :: l(0:n+1) - ! - integer(psb_ipk_) :: p,q,s,t - ! .. - iret = 0 - ! first step: we are preparing ordered sublists, exploiting - ! what order was already in the input data; negative links - ! mark the end of the sublists - l(0) = 1 - t = n + 1 - do p = 1,n - 1 - if (k(p) >= k(p+1)) then - l(p) = p + 1 - else - l(t) = - (p+1) - t = p + ! otherwise, one sublist ended, and we append to it the rest + ! of the other one. + l(s) = q + s = t + do + t = q + q = l(q) + if (q <= 0) exit + end do end if - end do - l(t) = 0 - l(n) = 0 - ! see if the input was already sorted - if (l(n+1) == 0) then - iret = 1 - return - else - l(n+1) = abs(l(n+1)) - end if - mergepass: do - ! otherwise, begin a pass through the list. - ! throughout all the subroutine we have: - ! p, q: pointing to the sublists being merged - ! s: pointing to the most recently processed record - ! t: pointing to the end of previously completed sublist - s = 0 - t = n + 1 - p = l(s) - q = l(t) - if (q == 0) exit mergepass + p = -p + q = -q + if (q == 0) then + l(s) = sign(p,l(s)) + l(t) = 0 + exit outer + end if + end do outer + end do mergepass - outer: do +end subroutine psi_d_msort_up - if (k(p) < k(q)) then +subroutine psi_d_msort_dw(n,k,l,iret) + use psb_const_mod + implicit none + integer(psb_ipk_) :: n, iret + real(psb_dpk_) :: k(n) + integer(psb_ipk_) :: l(0:n+1) + ! + integer(psb_ipk_) :: p,q,s,t + ! .. + iret = 0 + ! first step: we are preparing ordered sublists, exploiting + ! what order was already in the input data; negative links + ! mark the end of the sublists + l(0) = 1 + t = n + 1 + do p = 1,n - 1 + if (k(p) >= k(p+1)) then + l(p) = p + 1 + else + l(t) = - (p+1) + t = p + end if + end do + l(t) = 0 + l(n) = 0 + ! see if the input was already sorted + if (l(n+1) == 0) then + iret = 1 + return + else + l(n+1) = abs(l(n+1)) + end if + + mergepass: do + ! otherwise, begin a pass through the list. + ! throughout all the subroutine we have: + ! p, q: pointing to the sublists being merged + ! s: pointing to the most recently processed record + ! t: pointing to the end of previously completed sublist + s = 0 + t = n + 1 + p = l(s) + q = l(t) + if (q == 0) exit mergepass - l(s) = sign(q,l(s)) - s = q - q = l(q) - if (q > 0) then - do - if (k(p) >= k(q)) cycle outer - s = q - q = l(q) - if (q <= 0) exit - end do - end if - l(s) = p - s = t - do - t = p - p = l(p) - if (p <= 0) exit - end do + outer: do - else + if (k(p) < k(q)) then - l(s) = sign(p,l(s)) - s = p - p = l(p) - if (p>0) then - do - if (k(p) < k(q)) cycle outer - s = p - p = l(p) - if (p <= 0) exit - end do - end if - ! otherwise, one sublist ended, and we append to it the rest - ! of the other one. - l(s) = q - s = t + l(s) = sign(q,l(s)) + s = q + q = l(q) + if (q > 0) then do - t = q + if (k(p) >= k(q)) cycle outer + s = q q = l(q) if (q <= 0) exit end do end if + l(s) = p + s = t + do + t = p + p = l(p) + if (p <= 0) exit + end do - p = -p - q = -q - if (q == 0) then - l(s) = sign(p,l(s)) - l(t) = 0 - exit outer + else + + l(s) = sign(p,l(s)) + s = p + p = l(p) + if (p>0) then + do + if (k(p) < k(q)) cycle outer + s = p + p = l(p) + if (p <= 0) exit + end do end if - end do outer - end do mergepass - - end subroutine psi_d_msort_dw - - subroutine psi_d_amsort_up(n,k,l,iret) - use psb_const_mod - implicit none - integer(psb_ipk_) :: n, iret - real(psb_dpk_) :: k(n) - integer(psb_ipk_) :: l(0:n+1) - ! - integer(psb_ipk_) :: p,q,s,t - ! .. - iret = 0 - ! first step: we are preparing ordered sublists, exploiting - ! what order was already in the input data; negative links - ! mark the end of the sublists - l(0) = 1 - t = n + 1 - do p = 1,n - 1 - if (abs(k(p)) <= abs(k(p+1))) then - l(p) = p + 1 - else - l(t) = - (p+1) - t = p + ! otherwise, one sublist ended, and we append to it the rest + ! of the other one. + l(s) = q + s = t + do + t = q + q = l(q) + if (q <= 0) exit + end do end if - end do - l(t) = 0 - l(n) = 0 - ! see if the input was already sorted - if (l(n+1) == 0) then - iret = 1 - return - else - l(n+1) = abs(l(n+1)) - end if - mergepass: do - ! otherwise, begin a pass through the list. - ! throughout all the subroutine we have: - ! p, q: pointing to the sublists being merged - ! s: pointing to the most recently processed record - ! t: pointing to the end of previously completed sublist - s = 0 - t = n + 1 - p = l(s) - q = l(t) - if (q == 0) exit mergepass + p = -p + q = -q + if (q == 0) then + l(s) = sign(p,l(s)) + l(t) = 0 + exit outer + end if + end do outer + end do mergepass - outer: do +end subroutine psi_d_msort_dw - if (abs(k(p)) > abs(k(q))) then +subroutine psi_d_amsort_up(n,k,l,iret) + use psb_const_mod + implicit none + integer(psb_ipk_) :: n, iret + real(psb_dpk_) :: k(n) + integer(psb_ipk_) :: l(0:n+1) + ! + integer(psb_ipk_) :: p,q,s,t + ! .. + iret = 0 + ! first step: we are preparing ordered sublists, exploiting + ! what order was already in the input data; negative links + ! mark the end of the sublists + l(0) = 1 + t = n + 1 + do p = 1,n - 1 + if (abs(k(p)) <= abs(k(p+1))) then + l(p) = p + 1 + else + l(t) = - (p+1) + t = p + end if + end do + l(t) = 0 + l(n) = 0 + ! see if the input was already sorted + if (l(n+1) == 0) then + iret = 1 + return + else + l(n+1) = abs(l(n+1)) + end if + + mergepass: do + ! otherwise, begin a pass through the list. + ! throughout all the subroutine we have: + ! p, q: pointing to the sublists being merged + ! s: pointing to the most recently processed record + ! t: pointing to the end of previously completed sublist + s = 0 + t = n + 1 + p = l(s) + q = l(t) + if (q == 0) exit mergepass - l(s) = sign(q,l(s)) - s = q - q = l(q) - if (q > 0) then - do - if (abs(k(p)) <= abs(k(q))) cycle outer - s = q - q = l(q) - if (q <= 0) exit - end do - end if - l(s) = p - s = t - do - t = p - p = l(p) - if (p <= 0) exit - end do + outer: do - else + if (abs(k(p)) > abs(k(q))) then - l(s) = sign(p,l(s)) - s = p - p = l(p) - if (p>0) then - do - if (abs(k(p)) > abs(k(q))) cycle outer - s = p - p = l(p) - if (p <= 0) exit - end do - end if - ! otherwise, one sublist ended, and we append to it the rest - ! of the other one. - l(s) = q - s = t + l(s) = sign(q,l(s)) + s = q + q = l(q) + if (q > 0) then do - t = q + if (abs(k(p)) <= abs(k(q))) cycle outer + s = q q = l(q) if (q <= 0) exit end do end if + l(s) = p + s = t + do + t = p + p = l(p) + if (p <= 0) exit + end do + + else - p = -p - q = -q - if (q == 0) then - l(s) = sign(p,l(s)) - l(t) = 0 - exit outer + l(s) = sign(p,l(s)) + s = p + p = l(p) + if (p>0) then + do + if (abs(k(p)) > abs(k(q))) cycle outer + s = p + p = l(p) + if (p <= 0) exit + end do end if - end do outer - end do mergepass - - end subroutine psi_d_amsort_up - - subroutine psi_d_amsort_dw(n,k,l,iret) - use psb_const_mod - implicit none - integer(psb_ipk_) :: n, iret - real(psb_dpk_) :: k(n) - integer(psb_ipk_) :: l(0:n+1) - ! - integer(psb_ipk_) :: p,q,s,t - ! .. - iret = 0 - ! first step: we are preparing ordered sublists, exploiting - ! what order was already in the input data; negative links - ! mark the end of the sublists - l(0) = 1 - t = n + 1 - do p = 1,n - 1 - if (abs(k(p)) >= abs(k(p+1))) then - l(p) = p + 1 - else - l(t) = - (p+1) - t = p + ! otherwise, one sublist ended, and we append to it the rest + ! of the other one. + l(s) = q + s = t + do + t = q + q = l(q) + if (q <= 0) exit + end do end if - end do - l(t) = 0 - l(n) = 0 - ! see if the input was already sorted - if (l(n+1) == 0) then - iret = 1 - return - else - l(n+1) = abs(l(n+1)) - end if - mergepass: do - ! otherwise, begin a pass through the list. - ! throughout all the subroutine we have: - ! p, q: pointing to the sublists being merged - ! s: pointing to the most recently processed record - ! t: pointing to the end of previously completed sublist - s = 0 - t = n + 1 - p = l(s) - q = l(t) - if (q == 0) exit mergepass + p = -p + q = -q + if (q == 0) then + l(s) = sign(p,l(s)) + l(t) = 0 + exit outer + end if + end do outer + end do mergepass - outer: do +end subroutine psi_d_amsort_up - if (abs(k(p)) < abs(k(q))) then +subroutine psi_d_amsort_dw(n,k,l,iret) + use psb_const_mod + implicit none + integer(psb_ipk_) :: n, iret + real(psb_dpk_) :: k(n) + integer(psb_ipk_) :: l(0:n+1) + ! + integer(psb_ipk_) :: p,q,s,t + ! .. + iret = 0 + ! first step: we are preparing ordered sublists, exploiting + ! what order was already in the input data; negative links + ! mark the end of the sublists + l(0) = 1 + t = n + 1 + do p = 1,n - 1 + if (abs(k(p)) >= abs(k(p+1))) then + l(p) = p + 1 + else + l(t) = - (p+1) + t = p + end if + end do + l(t) = 0 + l(n) = 0 + ! see if the input was already sorted + if (l(n+1) == 0) then + iret = 1 + return + else + l(n+1) = abs(l(n+1)) + end if + + mergepass: do + ! otherwise, begin a pass through the list. + ! throughout all the subroutine we have: + ! p, q: pointing to the sublists being merged + ! s: pointing to the most recently processed record + ! t: pointing to the end of previously completed sublist + s = 0 + t = n + 1 + p = l(s) + q = l(t) + if (q == 0) exit mergepass - l(s) = sign(q,l(s)) - s = q - q = l(q) - if (q > 0) then - do - if (abs(k(p)) >= abs(k(q))) cycle outer - s = q - q = l(q) - if (q <= 0) exit - end do - end if - l(s) = p - s = t - do - t = p - p = l(p) - if (p <= 0) exit - end do + outer: do - else + if (abs(k(p)) < abs(k(q))) then - l(s) = sign(p,l(s)) - s = p - p = l(p) - if (p>0) then - do - if (abs(k(p)) < abs(k(q))) cycle outer - s = p - p = l(p) - if (p <= 0) exit - end do - end if - ! otherwise, one sublist ended, and we append to it the rest - ! of the other one. - l(s) = q - s = t + l(s) = sign(q,l(s)) + s = q + q = l(q) + if (q > 0) then do - t = q + if (abs(k(p)) >= abs(k(q))) cycle outer + s = q q = l(q) if (q <= 0) exit end do end if + l(s) = p + s = t + do + t = p + p = l(p) + if (p <= 0) exit + end do - p = -p - q = -q - if (q == 0) then - l(s) = sign(p,l(s)) - l(t) = 0 - exit outer - end if - end do outer - end do mergepass - - end subroutine psi_d_amsort_dw - - - + else + l(s) = sign(p,l(s)) + s = p + p = l(p) + if (p>0) then + do + if (abs(k(p)) < abs(k(q))) cycle outer + s = p + p = l(p) + if (p <= 0) exit + end do + end if + ! otherwise, one sublist ended, and we append to it the rest + ! of the other one. + l(s) = q + s = t + do + t = q + q = l(q) + if (q <= 0) exit + end do + end if + p = -p + q = -q + if (q == 0) then + l(s) = sign(p,l(s)) + l(t) = 0 + exit outer + end if + end do outer + end do mergepass +end subroutine psi_d_amsort_dw diff --git a/base/serial/sort/psb_e_hsort_impl.f90 b/base/serial/sort/psb_e_hsort_impl.f90 index d6990f36..f1a1a78f 100644 --- a/base/serial/sort/psb_e_hsort_impl.f90 +++ b/base/serial/sort/psb_e_hsort_impl.f90 @@ -49,7 +49,8 @@ subroutine psb_ehsort(x,ix,dir,flag) integer(psb_ipk_), optional, intent(in) :: dir, flag integer(psb_epk_), optional, intent(inout) :: ix(:) - integer(psb_ipk_) :: dir_, flag_, n, i, l, err_act,info + integer(psb_ipk_) :: flag_, n, i, err_act,info + integer(psb_epk_) :: dir_, l integer(psb_epk_) :: key integer(psb_epk_) :: index @@ -159,7 +160,7 @@ end subroutine psb_ehsort ! ! These are packaged so that they can be used to implement -! a heapsort, should the need arise +! a heapsort. ! ! ! Programming note: @@ -196,7 +197,7 @@ subroutine psi_e_insert_heap(key,last,heap,dir,info) ! dir: sorting direction integer(psb_epk_), intent(in) :: key - integer(psb_ipk_), intent(in) :: dir + integer(psb_epk_), intent(in) :: dir integer(psb_epk_), intent(inout) :: heap(:) integer(psb_epk_), intent(inout) :: last integer(psb_ipk_), intent(out) :: info @@ -296,7 +297,7 @@ subroutine psi_e_heap_get_first(key,last,heap,dir,info) integer(psb_epk_), intent(inout) :: key integer(psb_epk_), intent(inout) :: last - integer(psb_ipk_), intent(in) :: dir + integer(psb_epk_), intent(in) :: dir integer(psb_epk_), intent(inout) :: heap(:) integer(psb_ipk_), intent(out) :: info @@ -428,9 +429,9 @@ subroutine psi_e_idx_insert_heap(key,index,last,heap,idxs,dir,info) ! dir: sorting direction integer(psb_epk_), intent(in) :: key - integer(psb_ipk_), intent(in) :: index,dir + integer(psb_epk_), intent(in) :: index,dir integer(psb_epk_), intent(inout) :: heap(:) - integer(psb_ipk_), intent(inout) :: idxs(:),last + integer(psb_epk_), intent(inout) :: idxs(:),last integer(psb_ipk_), intent(out) :: info integer(psb_ipk_) :: i, i2, itemp integer(psb_epk_) :: temp @@ -540,11 +541,12 @@ subroutine psi_e_idx_heap_get_first(key,index,last,heap,idxs,dir,info) use psb_sort_mod, psb_protect_name => psi_e_idx_heap_get_first implicit none + integer(psb_epk_), intent(inout) :: key integer(psb_epk_), intent(inout) :: heap(:) - integer(psb_ipk_), intent(out) :: index,info - integer(psb_ipk_), intent(inout) :: last,idxs(:) - integer(psb_ipk_), intent(in) :: dir - integer(psb_epk_), intent(out) :: key + integer(psb_epk_), intent(out) :: index + integer(psb_ipk_), intent(out) :: info + integer(psb_epk_), intent(inout) :: last,idxs(:) + integer(psb_epk_), intent(in) :: dir integer(psb_ipk_) :: i, j,itemp integer(psb_epk_) :: temp diff --git a/base/serial/sort/psb_e_msort_impl.f90 b/base/serial/sort/psb_e_msort_impl.f90 index 2950d7bd..d8cd6404 100644 --- a/base/serial/sort/psb_e_msort_impl.f90 +++ b/base/serial/sort/psb_e_msort_impl.f90 @@ -1,34 +1,34 @@ -! -! 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. -! -! + ! + ! 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. + ! + ! ! ! The merge-sort routines ! References: @@ -40,674 +40,668 @@ ! Data Structures and Algorithms ! Addison-Wesley ! - logical function psb_eisaperm(n,eip) - use psb_sort_mod, psb_protect_name => psb_eisaperm - implicit none - - integer(psb_epk_), intent(in) :: n - integer(psb_epk_), intent(in) :: eip(n) - integer(psb_epk_), allocatable :: ip(:) - integer(psb_epk_) :: i,j,m, info - - - psb_eisaperm = .true. - if (n <= 0) return - allocate(ip(n), stat=info) - if (info /= psb_success_) return - ! - ! sanity check first - ! - do i=1, n - ip(i) = eip(i) - if ((ip(i) < 1).or.(ip(i) > n)) then - write(psb_err_unit,*) 'Out of bounds in isaperm' ,ip(i), n - psb_eisaperm = .false. - return - endif - enddo +logical function psb_eisaperm(n,eip) + use psb_sort_mod, psb_protect_name => psb_eisaperm + implicit none - ! - ! now work through the cycles, by marking each successive item as negative. - ! no cycle should intersect with any other, hence the >= 1 check. - ! - do m = 1, n - i = ip(m) - if (i < 0) then - ip(m) = -i - else if (i /= m) then - j = ip(i) - ip(i) = -j - i = j - do while ((j >= 1).and.(j /= m)) - j = ip(i) - ip(i) = -j - i = j - enddo - ip(m) = abs(ip(m)) - if (j /= m) then - psb_eisaperm = .false. - goto 9999 - endif - end if - enddo -9999 continue + integer(psb_epk_), intent(in) :: n + integer(psb_epk_), intent(in) :: eip(n) + integer(psb_epk_), allocatable :: ip(:) + integer(psb_epk_) :: i,j,m, info - return - end function psb_eisaperm + psb_eisaperm = .true. + if (n <= 0) return + allocate(ip(n), stat=info) + if (info /= psb_success_) return + ! + ! sanity check first + ! + do i=1, n + ip(i) = eip(i) + if ((ip(i) < 1).or.(ip(i) > n)) then + write(psb_err_unit,*) 'Out of bounds in isaperm' ,ip(i), n + psb_eisaperm = .false. + return + endif + enddo - subroutine psb_emsort_u(x,nout,dir) - use psb_sort_mod, psb_protect_name => psb_emsort_u - use psb_error_mod - implicit none - integer(psb_epk_), intent(inout) :: x(:) - integer(psb_epk_), intent(out) :: nout - integer(psb_ipk_), optional, intent(in) :: dir + ! + ! now work through the cycles, by marking each successive item as negative. + ! no cycle should intersect with any other, hence the >= 1 check. + ! + do m = 1, n + i = ip(m) + if (i < 0) then + ip(m) = -i + else if (i /= m) then + j = ip(i) + ip(i) = -j + i = j + do while ((j >= 1).and.(j /= m)) + j = ip(i) + ip(i) = -j + i = j + enddo + ip(m) = abs(ip(m)) + if (j /= m) then + psb_eisaperm = .false. + goto 9999 + endif + end if + enddo +9999 continue - integer(psb_epk_) :: n, k - integer(psb_ipk_) :: err_act + return +end function psb_eisaperm - integer(psb_ipk_) :: ierr(5) - character(len=20) :: name - name='psb_msort_u' - call psb_erractionsave(err_act) +subroutine psb_emsort_u(x,nout,dir) + use psb_sort_mod, psb_protect_name => psb_emsort_u + use psb_error_mod + implicit none + integer(psb_epk_), intent(inout) :: x(:) + integer(psb_epk_), intent(out) :: nout + integer(psb_ipk_), optional, intent(in) :: dir - n = size(x) + integer(psb_epk_) :: n, k + integer(psb_ipk_) :: err_act - call psb_msort(x,dir=dir) - nout = min(1,n) - do k=2,n - if (x(k) /= x(nout)) then - nout = nout + 1 - x(nout) = x(k) - endif - enddo + integer(psb_ipk_) :: ierr(5) + character(len=20) :: name - return + name='psb_msort_u' + call psb_erractionsave(err_act) -9999 call psb_error_handler(err_act) + n = size(x) - return - end subroutine psb_emsort_u + call psb_msort(x,dir=dir) + nout = min(1,n) + do k=2,n + if (x(k) /= x(nout)) then + nout = nout + 1 + x(nout) = x(k) + endif + enddo + return - function psb_ebsrch(key,n,v) result(ipos) - use psb_sort_mod, psb_protect_name => psb_ebsrch - implicit none - integer(psb_ipk_) :: ipos, n - integer(psb_epk_) :: key - integer(psb_epk_) :: v(:) +9999 call psb_error_handler(err_act) - integer(psb_ipk_) :: lb, ub, m, i + return +end subroutine psb_emsort_u - ipos = -1 - if (n<5) then - do i=1,n - if (key.eq.v(i)) then - ipos = i - return - end if - enddo - return - end if - - lb = 1 - ub = n - - do while (lb.le.ub) - m = (lb+ub)/2 - if (key.eq.v(m)) then - ipos = m - lb = ub + 1 - else if (key < v(m)) then - ub = m-1 - else - lb = m + 1 - end if - enddo - return - end function psb_ebsrch - function psb_essrch(key,n,v) result(ipos) - use psb_sort_mod, psb_protect_name => psb_essrch - implicit none - integer(psb_ipk_) :: ipos, n - integer(psb_epk_) :: key - integer(psb_epk_) :: v(:) +function psb_ebsrch(key,n,v) result(ipos) + use psb_sort_mod, psb_protect_name => psb_ebsrch + implicit none + integer(psb_ipk_) :: ipos, n + integer(psb_epk_) :: key + integer(psb_epk_) :: v(:) - integer(psb_ipk_) :: i + integer(psb_ipk_) :: lb, ub, m, i - ipos = -1 + ipos = -1 + if (n<5) then do i=1,n if (key.eq.v(i)) then ipos = i return end if enddo - return - end function psb_essrch - - subroutine psb_emsort(x,ix,dir,flag) - use psb_sort_mod, psb_protect_name => psb_emsort - use psb_error_mod - use psb_ip_reord_mod - implicit none - integer(psb_epk_), intent(inout) :: x(:) - integer(psb_ipk_), optional, intent(in) :: dir, flag - integer(psb_epk_), optional, intent(inout) :: ix(:) - - integer(psb_ipk_) :: dir_, flag_, n, err_act - - integer(psb_epk_), allocatable :: iaux(:) - integer(psb_ipk_) :: iret, info, i - integer(psb_ipk_) :: ierr(5) - character(len=20) :: name - - name='psb_emsort' - call psb_erractionsave(err_act) - - if (present(dir)) then - dir_ = dir - else - dir_= psb_sort_up_ + end if + + lb = 1 + ub = n + + do while (lb.le.ub) + m = (lb+ub)/2 + if (key.eq.v(m)) then + ipos = m + lb = ub + 1 + else if (key < v(m)) then + ub = m-1 + else + lb = m + 1 end if - select case(dir_) - case( psb_sort_up_, psb_sort_down_, psb_asort_up_, psb_asort_down_) + enddo + return +end function psb_ebsrch + +function psb_essrch(key,n,v) result(ipos) + use psb_sort_mod, psb_protect_name => psb_essrch + implicit none + integer(psb_ipk_) :: ipos, n + integer(psb_epk_) :: key + integer(psb_epk_) :: v(:) + + integer(psb_ipk_) :: i + + ipos = -1 + do i=1,n + if (key.eq.v(i)) then + ipos = i + return + end if + enddo + + return +end function psb_essrch + +subroutine psb_emsort(x,ix,dir,flag) + use psb_sort_mod, psb_protect_name => psb_emsort + use psb_error_mod + use psb_ip_reord_mod + implicit none + integer(psb_epk_), intent(inout) :: x(:) + integer(psb_ipk_), optional, intent(in) :: dir, flag + integer(psb_epk_), optional, intent(inout) :: ix(:) + + integer(psb_ipk_) :: dir_, flag_, n, err_act + + integer(psb_epk_), allocatable :: iaux(:) + integer(psb_ipk_) :: iret, info, i + integer(psb_ipk_) :: ierr(5) + character(len=20) :: name + + name='psb_emsort' + call psb_erractionsave(err_act) + + if (present(dir)) then + dir_ = dir + else + dir_= psb_sort_up_ + end if + select case(dir_) + case( psb_sort_up_, psb_sort_down_, psb_asort_up_, psb_asort_down_) + ! OK keep going + case default + ierr(1) = 3; ierr(2) = dir_; + call psb_errpush(psb_err_input_value_invalid_i_,name,i_err=ierr) + goto 9999 + end select + + n = size(x) + + if (present(ix)) then + if (size(ix) < n) then + ierr(1) = 2; ierr(2) = size(ix); + call psb_errpush(psb_err_input_asize_invalid_i_,name,i_err=ierr) + goto 9999 + end if + if (present(flag)) then + flag_ = flag + else + flag_ = psb_sort_ovw_idx_ + end if + select case(flag_) + case(psb_sort_ovw_idx_) + do i=1,n + ix(i) = i + end do + case (psb_sort_keep_idx_) ! OK keep going case default - ierr(1) = 3; ierr(2) = dir_; + ierr(1) = 4; ierr(2) = flag_; call psb_errpush(psb_err_input_value_invalid_i_,name,i_err=ierr) goto 9999 end select - n = size(x) - - if (present(ix)) then - if (size(ix) < n) then - ierr(1) = 2; ierr(2) = size(ix); - call psb_errpush(psb_err_input_asize_invalid_i_,name,i_err=ierr) - goto 9999 - end if - if (present(flag)) then - flag_ = flag - else - flag_ = psb_sort_ovw_idx_ - end if - select case(flag_) - case(psb_sort_ovw_idx_) - do i=1,n - ix(i) = i - end do - case (psb_sort_keep_idx_) - ! OK keep going - case default - ierr(1) = 4; ierr(2) = flag_; - call psb_errpush(psb_err_input_value_invalid_i_,name,i_err=ierr) - goto 9999 - end select - - end if - - allocate(iaux(0:n+1),stat=info) - if (info /= psb_success_) then - call psb_errpush(psb_err_alloc_dealloc_,r_name='psb_e_msort') - goto 9999 - endif - - select case(dir_) - case (psb_sort_up_) - call psi_e_msort_up(n,x,iaux,iret) - case (psb_sort_down_) - call psi_e_msort_dw(n,x,iaux,iret) - case (psb_asort_up_) - call psi_e_amsort_up(n,x,iaux,iret) - case (psb_asort_down_) - call psi_e_amsort_dw(n,x,iaux,iret) - end select - ! - ! Do the actual reordering, since the inner routines - ! only provide linked pointers. - ! - if (iret == 0 ) then - if (present(ix)) then - call psb_ip_reord(n,x,ix,iaux) - else - call psb_ip_reord(n,x,iaux) - end if + end if + + allocate(iaux(0:n+1),stat=info) + if (info /= psb_success_) then + call psb_errpush(psb_err_alloc_dealloc_,r_name='psb_e_msort') + goto 9999 + endif + + select case(dir_) + case (psb_sort_up_) + call psi_e_msort_up(n,x,iaux,iret) + case (psb_sort_down_) + call psi_e_msort_dw(n,x,iaux,iret) + case (psb_asort_up_) + call psi_e_amsort_up(n,x,iaux,iret) + case (psb_asort_down_) + call psi_e_amsort_dw(n,x,iaux,iret) + end select + ! + ! Do the actual reordering, since the inner routines + ! only provide linked pointers. + ! + if (iret == 0 ) then + if (present(ix)) then + call psb_ip_reord(n,x,ix,iaux) + else + call psb_ip_reord(n,x,iaux) end if + end if - return + return 9999 call psb_error_handler(err_act) - return + return - end subroutine psb_emsort - - subroutine psi_e_msort_up(n,k,l,iret) - use psb_const_mod - implicit none - integer(psb_ipk_) :: n, iret - integer(psb_epk_) :: k(n) - integer(psb_epk_) :: l(0:n+1) - ! - integer(psb_epk_) :: p,q,s,t - ! .. - iret = 0 - ! first step: we are preparing ordered sublists, exploiting - ! what order was already in the input data; negative links - ! mark the end of the sublists - l(0) = 1 - t = n + 1 - do p = 1,n - 1 - if (k(p) <= k(p+1)) then - l(p) = p + 1 - else - l(t) = - (p+1) - t = p - end if - end do - l(t) = 0 - l(n) = 0 - ! see if the input was already sorted - if (l(n+1) == 0) then - iret = 1 - return +end subroutine psb_emsort + +subroutine psi_e_msort_up(n,k,l,iret) + use psb_const_mod + implicit none + integer(psb_ipk_) :: n, iret + integer(psb_epk_) :: k(n) + integer(psb_epk_) :: l(0:n+1) + ! + integer(psb_epk_) :: p,q,s,t + ! .. + iret = 0 + ! first step: we are preparing ordered sublists, exploiting + ! what order was already in the input data; negative links + ! mark the end of the sublists + l(0) = 1 + t = n + 1 + do p = 1,n - 1 + if (k(p) <= k(p+1)) then + l(p) = p + 1 else - l(n+1) = abs(l(n+1)) + l(t) = - (p+1) + t = p end if + end do + l(t) = 0 + l(n) = 0 + ! see if the input was already sorted + if (l(n+1) == 0) then + iret = 1 + return + else + l(n+1) = abs(l(n+1)) + end if + + mergepass: do + ! otherwise, begin a pass through the list. + ! throughout all the subroutine we have: + ! p, q: pointing to the sublists being merged + ! s: pointing to the most recently processed record + ! t: pointing to the end of previously completed sublist + s = 0 + t = n + 1 + p = l(s) + q = l(t) + if (q == 0) exit mergepass - mergepass: do - ! otherwise, begin a pass through the list. - ! throughout all the subroutine we have: - ! p, q: pointing to the sublists being merged - ! s: pointing to the most recently processed record - ! t: pointing to the end of previously completed sublist - s = 0 - t = n + 1 - p = l(s) - q = l(t) - if (q == 0) exit mergepass - - outer: do + outer: do - if (k(p) > k(q)) then + if (k(p) > k(q)) then - l(s) = sign(q,l(s)) - s = q - q = l(q) - if (q > 0) then - do - if (k(p) <= k(q)) cycle outer - s = q - q = l(q) - if (q <= 0) exit - end do - end if - l(s) = p - s = t + l(s) = sign(q,l(s)) + s = q + q = l(q) + if (q > 0) then do - t = p - p = l(p) - if (p <= 0) exit - end do - - else - - l(s) = sign(p,l(s)) - s = p - p = l(p) - if (p>0) then - do - if (k(p) > k(q)) cycle outer - s = p - p = l(p) - if (p <= 0) exit - end do - end if - ! otherwise, one sublist ended, and we append to it the rest - ! of the other one. - l(s) = q - s = t - do - t = q + if (k(p) <= k(q)) cycle outer + s = q q = l(q) if (q <= 0) exit end do end if + l(s) = p + s = t + do + t = p + p = l(p) + if (p <= 0) exit + end do - p = -p - q = -q - if (q == 0) then - l(s) = sign(p,l(s)) - l(t) = 0 - exit outer + else + + l(s) = sign(p,l(s)) + s = p + p = l(p) + if (p>0) then + do + if (k(p) > k(q)) cycle outer + s = p + p = l(p) + if (p <= 0) exit + end do end if - end do outer - end do mergepass - - end subroutine psi_e_msort_up - - subroutine psi_e_msort_dw(n,k,l,iret) - use psb_const_mod - implicit none - integer(psb_ipk_) :: n, iret - integer(psb_epk_) :: k(n) - integer(psb_epk_) :: l(0:n+1) - ! - integer(psb_epk_) :: p,q,s,t - ! .. - iret = 0 - ! first step: we are preparing ordered sublists, exploiting - ! what order was already in the input data; negative links - ! mark the end of the sublists - l(0) = 1 - t = n + 1 - do p = 1,n - 1 - if (k(p) >= k(p+1)) then - l(p) = p + 1 - else - l(t) = - (p+1) - t = p + ! otherwise, one sublist ended, and we append to it the rest + ! of the other one. + l(s) = q + s = t + do + t = q + q = l(q) + if (q <= 0) exit + end do end if - end do - l(t) = 0 - l(n) = 0 - ! see if the input was already sorted - if (l(n+1) == 0) then - iret = 1 - return - else - l(n+1) = abs(l(n+1)) - end if - mergepass: do - ! otherwise, begin a pass through the list. - ! throughout all the subroutine we have: - ! p, q: pointing to the sublists being merged - ! s: pointing to the most recently processed record - ! t: pointing to the end of previously completed sublist - s = 0 - t = n + 1 - p = l(s) - q = l(t) - if (q == 0) exit mergepass + p = -p + q = -q + if (q == 0) then + l(s) = sign(p,l(s)) + l(t) = 0 + exit outer + end if + end do outer + end do mergepass - outer: do +end subroutine psi_e_msort_up - if (k(p) < k(q)) then +subroutine psi_e_msort_dw(n,k,l,iret) + use psb_const_mod + implicit none + integer(psb_ipk_) :: n, iret + integer(psb_epk_) :: k(n) + integer(psb_epk_) :: l(0:n+1) + ! + integer(psb_epk_) :: p,q,s,t + ! .. + iret = 0 + ! first step: we are preparing ordered sublists, exploiting + ! what order was already in the input data; negative links + ! mark the end of the sublists + l(0) = 1 + t = n + 1 + do p = 1,n - 1 + if (k(p) >= k(p+1)) then + l(p) = p + 1 + else + l(t) = - (p+1) + t = p + end if + end do + l(t) = 0 + l(n) = 0 + ! see if the input was already sorted + if (l(n+1) == 0) then + iret = 1 + return + else + l(n+1) = abs(l(n+1)) + end if + + mergepass: do + ! otherwise, begin a pass through the list. + ! throughout all the subroutine we have: + ! p, q: pointing to the sublists being merged + ! s: pointing to the most recently processed record + ! t: pointing to the end of previously completed sublist + s = 0 + t = n + 1 + p = l(s) + q = l(t) + if (q == 0) exit mergepass - l(s) = sign(q,l(s)) - s = q - q = l(q) - if (q > 0) then - do - if (k(p) >= k(q)) cycle outer - s = q - q = l(q) - if (q <= 0) exit - end do - end if - l(s) = p - s = t - do - t = p - p = l(p) - if (p <= 0) exit - end do + outer: do - else + if (k(p) < k(q)) then - l(s) = sign(p,l(s)) - s = p - p = l(p) - if (p>0) then - do - if (k(p) < k(q)) cycle outer - s = p - p = l(p) - if (p <= 0) exit - end do - end if - ! otherwise, one sublist ended, and we append to it the rest - ! of the other one. - l(s) = q - s = t + l(s) = sign(q,l(s)) + s = q + q = l(q) + if (q > 0) then do - t = q + if (k(p) >= k(q)) cycle outer + s = q q = l(q) if (q <= 0) exit end do end if + l(s) = p + s = t + do + t = p + p = l(p) + if (p <= 0) exit + end do - p = -p - q = -q - if (q == 0) then - l(s) = sign(p,l(s)) - l(t) = 0 - exit outer + else + + l(s) = sign(p,l(s)) + s = p + p = l(p) + if (p>0) then + do + if (k(p) < k(q)) cycle outer + s = p + p = l(p) + if (p <= 0) exit + end do end if - end do outer - end do mergepass - - end subroutine psi_e_msort_dw - - subroutine psi_e_amsort_up(n,k,l,iret) - use psb_const_mod - implicit none - integer(psb_ipk_) :: n, iret - integer(psb_epk_) :: k(n) - integer(psb_epk_) :: l(0:n+1) - ! - integer(psb_epk_) :: p,q,s,t - ! .. - iret = 0 - ! first step: we are preparing ordered sublists, exploiting - ! what order was already in the input data; negative links - ! mark the end of the sublists - l(0) = 1 - t = n + 1 - do p = 1,n - 1 - if (abs(k(p)) <= abs(k(p+1))) then - l(p) = p + 1 - else - l(t) = - (p+1) - t = p + ! otherwise, one sublist ended, and we append to it the rest + ! of the other one. + l(s) = q + s = t + do + t = q + q = l(q) + if (q <= 0) exit + end do end if - end do - l(t) = 0 - l(n) = 0 - ! see if the input was already sorted - if (l(n+1) == 0) then - iret = 1 - return - else - l(n+1) = abs(l(n+1)) - end if - mergepass: do - ! otherwise, begin a pass through the list. - ! throughout all the subroutine we have: - ! p, q: pointing to the sublists being merged - ! s: pointing to the most recently processed record - ! t: pointing to the end of previously completed sublist - s = 0 - t = n + 1 - p = l(s) - q = l(t) - if (q == 0) exit mergepass + p = -p + q = -q + if (q == 0) then + l(s) = sign(p,l(s)) + l(t) = 0 + exit outer + end if + end do outer + end do mergepass - outer: do +end subroutine psi_e_msort_dw - if (abs(k(p)) > abs(k(q))) then +subroutine psi_e_amsort_up(n,k,l,iret) + use psb_const_mod + implicit none + integer(psb_ipk_) :: n, iret + integer(psb_epk_) :: k(n) + integer(psb_epk_) :: l(0:n+1) + ! + integer(psb_epk_) :: p,q,s,t + ! .. + iret = 0 + ! first step: we are preparing ordered sublists, exploiting + ! what order was already in the input data; negative links + ! mark the end of the sublists + l(0) = 1 + t = n + 1 + do p = 1,n - 1 + if (abs(k(p)) <= abs(k(p+1))) then + l(p) = p + 1 + else + l(t) = - (p+1) + t = p + end if + end do + l(t) = 0 + l(n) = 0 + ! see if the input was already sorted + if (l(n+1) == 0) then + iret = 1 + return + else + l(n+1) = abs(l(n+1)) + end if + + mergepass: do + ! otherwise, begin a pass through the list. + ! throughout all the subroutine we have: + ! p, q: pointing to the sublists being merged + ! s: pointing to the most recently processed record + ! t: pointing to the end of previously completed sublist + s = 0 + t = n + 1 + p = l(s) + q = l(t) + if (q == 0) exit mergepass - l(s) = sign(q,l(s)) - s = q - q = l(q) - if (q > 0) then - do - if (abs(k(p)) <= abs(k(q))) cycle outer - s = q - q = l(q) - if (q <= 0) exit - end do - end if - l(s) = p - s = t - do - t = p - p = l(p) - if (p <= 0) exit - end do + outer: do - else + if (abs(k(p)) > abs(k(q))) then - l(s) = sign(p,l(s)) - s = p - p = l(p) - if (p>0) then - do - if (abs(k(p)) > abs(k(q))) cycle outer - s = p - p = l(p) - if (p <= 0) exit - end do - end if - ! otherwise, one sublist ended, and we append to it the rest - ! of the other one. - l(s) = q - s = t + l(s) = sign(q,l(s)) + s = q + q = l(q) + if (q > 0) then do - t = q + if (abs(k(p)) <= abs(k(q))) cycle outer + s = q q = l(q) if (q <= 0) exit end do end if + l(s) = p + s = t + do + t = p + p = l(p) + if (p <= 0) exit + end do - p = -p - q = -q - if (q == 0) then - l(s) = sign(p,l(s)) - l(t) = 0 - exit outer + else + + l(s) = sign(p,l(s)) + s = p + p = l(p) + if (p>0) then + do + if (abs(k(p)) > abs(k(q))) cycle outer + s = p + p = l(p) + if (p <= 0) exit + end do end if - end do outer - end do mergepass - - end subroutine psi_e_amsort_up - - subroutine psi_e_amsort_dw(n,k,l,iret) - use psb_const_mod - implicit none - integer(psb_ipk_) :: n, iret - integer(psb_epk_) :: k(n) - integer(psb_epk_) :: l(0:n+1) - ! - integer(psb_epk_) :: p,q,s,t - ! .. - iret = 0 - ! first step: we are preparing ordered sublists, exploiting - ! what order was already in the input data; negative links - ! mark the end of the sublists - l(0) = 1 - t = n + 1 - do p = 1,n - 1 - if (abs(k(p)) >= abs(k(p+1))) then - l(p) = p + 1 - else - l(t) = - (p+1) - t = p + ! otherwise, one sublist ended, and we append to it the rest + ! of the other one. + l(s) = q + s = t + do + t = q + q = l(q) + if (q <= 0) exit + end do end if - end do - l(t) = 0 - l(n) = 0 - ! see if the input was already sorted - if (l(n+1) == 0) then - iret = 1 - return - else - l(n+1) = abs(l(n+1)) - end if - mergepass: do - ! otherwise, begin a pass through the list. - ! throughout all the subroutine we have: - ! p, q: pointing to the sublists being merged - ! s: pointing to the most recently processed record - ! t: pointing to the end of previously completed sublist - s = 0 - t = n + 1 - p = l(s) - q = l(t) - if (q == 0) exit mergepass + p = -p + q = -q + if (q == 0) then + l(s) = sign(p,l(s)) + l(t) = 0 + exit outer + end if + end do outer + end do mergepass - outer: do +end subroutine psi_e_amsort_up - if (abs(k(p)) < abs(k(q))) then +subroutine psi_e_amsort_dw(n,k,l,iret) + use psb_const_mod + implicit none + integer(psb_ipk_) :: n, iret + integer(psb_epk_) :: k(n) + integer(psb_epk_) :: l(0:n+1) + ! + integer(psb_epk_) :: p,q,s,t + ! .. + iret = 0 + ! first step: we are preparing ordered sublists, exploiting + ! what order was already in the input data; negative links + ! mark the end of the sublists + l(0) = 1 + t = n + 1 + do p = 1,n - 1 + if (abs(k(p)) >= abs(k(p+1))) then + l(p) = p + 1 + else + l(t) = - (p+1) + t = p + end if + end do + l(t) = 0 + l(n) = 0 + ! see if the input was already sorted + if (l(n+1) == 0) then + iret = 1 + return + else + l(n+1) = abs(l(n+1)) + end if + + mergepass: do + ! otherwise, begin a pass through the list. + ! throughout all the subroutine we have: + ! p, q: pointing to the sublists being merged + ! s: pointing to the most recently processed record + ! t: pointing to the end of previously completed sublist + s = 0 + t = n + 1 + p = l(s) + q = l(t) + if (q == 0) exit mergepass - l(s) = sign(q,l(s)) - s = q - q = l(q) - if (q > 0) then - do - if (abs(k(p)) >= abs(k(q))) cycle outer - s = q - q = l(q) - if (q <= 0) exit - end do - end if - l(s) = p - s = t - do - t = p - p = l(p) - if (p <= 0) exit - end do + outer: do - else + if (abs(k(p)) < abs(k(q))) then - l(s) = sign(p,l(s)) - s = p - p = l(p) - if (p>0) then - do - if (abs(k(p)) < abs(k(q))) cycle outer - s = p - p = l(p) - if (p <= 0) exit - end do - end if - ! otherwise, one sublist ended, and we append to it the rest - ! of the other one. - l(s) = q - s = t + l(s) = sign(q,l(s)) + s = q + q = l(q) + if (q > 0) then do - t = q + if (abs(k(p)) >= abs(k(q))) cycle outer + s = q q = l(q) if (q <= 0) exit end do end if + l(s) = p + s = t + do + t = p + p = l(p) + if (p <= 0) exit + end do - p = -p - q = -q - if (q == 0) then - l(s) = sign(p,l(s)) - l(t) = 0 - exit outer - end if - end do outer - end do mergepass - - end subroutine psi_e_amsort_dw - - - + else + l(s) = sign(p,l(s)) + s = p + p = l(p) + if (p>0) then + do + if (abs(k(p)) < abs(k(q))) cycle outer + s = p + p = l(p) + if (p <= 0) exit + end do + end if + ! otherwise, one sublist ended, and we append to it the rest + ! of the other one. + l(s) = q + s = t + do + t = q + q = l(q) + if (q <= 0) exit + end do + end if + p = -p + q = -q + if (q == 0) then + l(s) = sign(p,l(s)) + l(t) = 0 + exit outer + end if + end do outer + end do mergepass +end subroutine psi_e_amsort_dw diff --git a/base/serial/sort/psb_m_hsort_impl.f90 b/base/serial/sort/psb_m_hsort_impl.f90 index dad77210..5dc92082 100644 --- a/base/serial/sort/psb_m_hsort_impl.f90 +++ b/base/serial/sort/psb_m_hsort_impl.f90 @@ -49,7 +49,8 @@ subroutine psb_mhsort(x,ix,dir,flag) integer(psb_ipk_), optional, intent(in) :: dir, flag integer(psb_ipk_), optional, intent(inout) :: ix(:) - integer(psb_ipk_) :: dir_, flag_, n, i, l, err_act,info + integer(psb_ipk_) :: flag_, n, i, err_act,info + integer(psb_ipk_) :: dir_, l integer(psb_mpk_) :: key integer(psb_ipk_) :: index @@ -159,7 +160,7 @@ end subroutine psb_mhsort ! ! These are packaged so that they can be used to implement -! a heapsort, should the need arise +! a heapsort. ! ! ! Programming note: @@ -540,11 +541,12 @@ subroutine psi_m_idx_heap_get_first(key,index,last,heap,idxs,dir,info) use psb_sort_mod, psb_protect_name => psi_m_idx_heap_get_first implicit none + integer(psb_mpk_), intent(inout) :: key integer(psb_mpk_), intent(inout) :: heap(:) - integer(psb_ipk_), intent(out) :: index,info + integer(psb_ipk_), intent(out) :: index + integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(inout) :: last,idxs(:) integer(psb_ipk_), intent(in) :: dir - integer(psb_mpk_), intent(out) :: key integer(psb_ipk_) :: i, j,itemp integer(psb_mpk_) :: temp diff --git a/base/serial/sort/psb_m_msort_impl.f90 b/base/serial/sort/psb_m_msort_impl.f90 index abbe4049..cd99a3c5 100644 --- a/base/serial/sort/psb_m_msort_impl.f90 +++ b/base/serial/sort/psb_m_msort_impl.f90 @@ -1,34 +1,34 @@ -! -! 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. -! -! + ! + ! 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. + ! + ! ! ! The merge-sort routines ! References: @@ -40,674 +40,668 @@ ! Data Structures and Algorithms ! Addison-Wesley ! - logical function psb_misaperm(n,eip) - use psb_sort_mod, psb_protect_name => psb_misaperm - implicit none - - integer(psb_mpk_), intent(in) :: n - integer(psb_mpk_), intent(in) :: eip(n) - integer(psb_mpk_), allocatable :: ip(:) - integer(psb_mpk_) :: i,j,m, info - - - psb_misaperm = .true. - if (n <= 0) return - allocate(ip(n), stat=info) - if (info /= psb_success_) return - ! - ! sanity check first - ! - do i=1, n - ip(i) = eip(i) - if ((ip(i) < 1).or.(ip(i) > n)) then - write(psb_err_unit,*) 'Out of bounds in isaperm' ,ip(i), n - psb_misaperm = .false. - return - endif - enddo +logical function psb_misaperm(n,eip) + use psb_sort_mod, psb_protect_name => psb_misaperm + implicit none - ! - ! now work through the cycles, by marking each successive item as negative. - ! no cycle should intersect with any other, hence the >= 1 check. - ! - do m = 1, n - i = ip(m) - if (i < 0) then - ip(m) = -i - else if (i /= m) then - j = ip(i) - ip(i) = -j - i = j - do while ((j >= 1).and.(j /= m)) - j = ip(i) - ip(i) = -j - i = j - enddo - ip(m) = abs(ip(m)) - if (j /= m) then - psb_misaperm = .false. - goto 9999 - endif - end if - enddo -9999 continue + integer(psb_mpk_), intent(in) :: n + integer(psb_mpk_), intent(in) :: eip(n) + integer(psb_mpk_), allocatable :: ip(:) + integer(psb_mpk_) :: i,j,m, info - return - end function psb_misaperm + psb_misaperm = .true. + if (n <= 0) return + allocate(ip(n), stat=info) + if (info /= psb_success_) return + ! + ! sanity check first + ! + do i=1, n + ip(i) = eip(i) + if ((ip(i) < 1).or.(ip(i) > n)) then + write(psb_err_unit,*) 'Out of bounds in isaperm' ,ip(i), n + psb_misaperm = .false. + return + endif + enddo - subroutine psb_mmsort_u(x,nout,dir) - use psb_sort_mod, psb_protect_name => psb_mmsort_u - use psb_error_mod - implicit none - integer(psb_mpk_), intent(inout) :: x(:) - integer(psb_ipk_), intent(out) :: nout - integer(psb_ipk_), optional, intent(in) :: dir + ! + ! now work through the cycles, by marking each successive item as negative. + ! no cycle should intersect with any other, hence the >= 1 check. + ! + do m = 1, n + i = ip(m) + if (i < 0) then + ip(m) = -i + else if (i /= m) then + j = ip(i) + ip(i) = -j + i = j + do while ((j >= 1).and.(j /= m)) + j = ip(i) + ip(i) = -j + i = j + enddo + ip(m) = abs(ip(m)) + if (j /= m) then + psb_misaperm = .false. + goto 9999 + endif + end if + enddo +9999 continue - integer(psb_ipk_) :: n, k - integer(psb_ipk_) :: err_act + return +end function psb_misaperm - integer(psb_ipk_) :: ierr(5) - character(len=20) :: name - name='psb_msort_u' - call psb_erractionsave(err_act) +subroutine psb_mmsort_u(x,nout,dir) + use psb_sort_mod, psb_protect_name => psb_mmsort_u + use psb_error_mod + implicit none + integer(psb_mpk_), intent(inout) :: x(:) + integer(psb_ipk_), intent(out) :: nout + integer(psb_ipk_), optional, intent(in) :: dir - n = size(x) + integer(psb_ipk_) :: n, k + integer(psb_ipk_) :: err_act - call psb_msort(x,dir=dir) - nout = min(1,n) - do k=2,n - if (x(k) /= x(nout)) then - nout = nout + 1 - x(nout) = x(k) - endif - enddo + integer(psb_ipk_) :: ierr(5) + character(len=20) :: name - return + name='psb_msort_u' + call psb_erractionsave(err_act) -9999 call psb_error_handler(err_act) + n = size(x) - return - end subroutine psb_mmsort_u + call psb_msort(x,dir=dir) + nout = min(1,n) + do k=2,n + if (x(k) /= x(nout)) then + nout = nout + 1 + x(nout) = x(k) + endif + enddo + return - function psb_mbsrch(key,n,v) result(ipos) - use psb_sort_mod, psb_protect_name => psb_mbsrch - implicit none - integer(psb_ipk_) :: ipos, n - integer(psb_mpk_) :: key - integer(psb_mpk_) :: v(:) +9999 call psb_error_handler(err_act) - integer(psb_ipk_) :: lb, ub, m, i + return +end subroutine psb_mmsort_u - ipos = -1 - if (n<5) then - do i=1,n - if (key.eq.v(i)) then - ipos = i - return - end if - enddo - return - end if - - lb = 1 - ub = n - - do while (lb.le.ub) - m = (lb+ub)/2 - if (key.eq.v(m)) then - ipos = m - lb = ub + 1 - else if (key < v(m)) then - ub = m-1 - else - lb = m + 1 - end if - enddo - return - end function psb_mbsrch - function psb_mssrch(key,n,v) result(ipos) - use psb_sort_mod, psb_protect_name => psb_mssrch - implicit none - integer(psb_ipk_) :: ipos, n - integer(psb_mpk_) :: key - integer(psb_mpk_) :: v(:) +function psb_mbsrch(key,n,v) result(ipos) + use psb_sort_mod, psb_protect_name => psb_mbsrch + implicit none + integer(psb_ipk_) :: ipos, n + integer(psb_mpk_) :: key + integer(psb_mpk_) :: v(:) - integer(psb_ipk_) :: i + integer(psb_ipk_) :: lb, ub, m, i - ipos = -1 + ipos = -1 + if (n<5) then do i=1,n if (key.eq.v(i)) then ipos = i return end if enddo - return - end function psb_mssrch - - subroutine psb_mmsort(x,ix,dir,flag) - use psb_sort_mod, psb_protect_name => psb_mmsort - use psb_error_mod - use psb_ip_reord_mod - implicit none - integer(psb_mpk_), intent(inout) :: x(:) - integer(psb_ipk_), optional, intent(in) :: dir, flag - integer(psb_ipk_), optional, intent(inout) :: ix(:) - - integer(psb_ipk_) :: dir_, flag_, n, err_act - - integer(psb_ipk_), allocatable :: iaux(:) - integer(psb_ipk_) :: iret, info, i - integer(psb_ipk_) :: ierr(5) - character(len=20) :: name - - name='psb_mmsort' - call psb_erractionsave(err_act) - - if (present(dir)) then - dir_ = dir - else - dir_= psb_sort_up_ + end if + + lb = 1 + ub = n + + do while (lb.le.ub) + m = (lb+ub)/2 + if (key.eq.v(m)) then + ipos = m + lb = ub + 1 + else if (key < v(m)) then + ub = m-1 + else + lb = m + 1 end if - select case(dir_) - case( psb_sort_up_, psb_sort_down_, psb_asort_up_, psb_asort_down_) + enddo + return +end function psb_mbsrch + +function psb_mssrch(key,n,v) result(ipos) + use psb_sort_mod, psb_protect_name => psb_mssrch + implicit none + integer(psb_ipk_) :: ipos, n + integer(psb_mpk_) :: key + integer(psb_mpk_) :: v(:) + + integer(psb_ipk_) :: i + + ipos = -1 + do i=1,n + if (key.eq.v(i)) then + ipos = i + return + end if + enddo + + return +end function psb_mssrch + +subroutine psb_mmsort(x,ix,dir,flag) + use psb_sort_mod, psb_protect_name => psb_mmsort + use psb_error_mod + use psb_ip_reord_mod + implicit none + integer(psb_mpk_), intent(inout) :: x(:) + integer(psb_ipk_), optional, intent(in) :: dir, flag + integer(psb_ipk_), optional, intent(inout) :: ix(:) + + integer(psb_ipk_) :: dir_, flag_, n, err_act + + integer(psb_ipk_), allocatable :: iaux(:) + integer(psb_ipk_) :: iret, info, i + integer(psb_ipk_) :: ierr(5) + character(len=20) :: name + + name='psb_mmsort' + call psb_erractionsave(err_act) + + if (present(dir)) then + dir_ = dir + else + dir_= psb_sort_up_ + end if + select case(dir_) + case( psb_sort_up_, psb_sort_down_, psb_asort_up_, psb_asort_down_) + ! OK keep going + case default + ierr(1) = 3; ierr(2) = dir_; + call psb_errpush(psb_err_input_value_invalid_i_,name,i_err=ierr) + goto 9999 + end select + + n = size(x) + + if (present(ix)) then + if (size(ix) < n) then + ierr(1) = 2; ierr(2) = size(ix); + call psb_errpush(psb_err_input_asize_invalid_i_,name,i_err=ierr) + goto 9999 + end if + if (present(flag)) then + flag_ = flag + else + flag_ = psb_sort_ovw_idx_ + end if + select case(flag_) + case(psb_sort_ovw_idx_) + do i=1,n + ix(i) = i + end do + case (psb_sort_keep_idx_) ! OK keep going case default - ierr(1) = 3; ierr(2) = dir_; + ierr(1) = 4; ierr(2) = flag_; call psb_errpush(psb_err_input_value_invalid_i_,name,i_err=ierr) goto 9999 end select - n = size(x) - - if (present(ix)) then - if (size(ix) < n) then - ierr(1) = 2; ierr(2) = size(ix); - call psb_errpush(psb_err_input_asize_invalid_i_,name,i_err=ierr) - goto 9999 - end if - if (present(flag)) then - flag_ = flag - else - flag_ = psb_sort_ovw_idx_ - end if - select case(flag_) - case(psb_sort_ovw_idx_) - do i=1,n - ix(i) = i - end do - case (psb_sort_keep_idx_) - ! OK keep going - case default - ierr(1) = 4; ierr(2) = flag_; - call psb_errpush(psb_err_input_value_invalid_i_,name,i_err=ierr) - goto 9999 - end select - - end if - - allocate(iaux(0:n+1),stat=info) - if (info /= psb_success_) then - call psb_errpush(psb_err_alloc_dealloc_,r_name='psb_m_msort') - goto 9999 - endif - - select case(dir_) - case (psb_sort_up_) - call psi_m_msort_up(n,x,iaux,iret) - case (psb_sort_down_) - call psi_m_msort_dw(n,x,iaux,iret) - case (psb_asort_up_) - call psi_m_amsort_up(n,x,iaux,iret) - case (psb_asort_down_) - call psi_m_amsort_dw(n,x,iaux,iret) - end select - ! - ! Do the actual reordering, since the inner routines - ! only provide linked pointers. - ! - if (iret == 0 ) then - if (present(ix)) then - call psb_ip_reord(n,x,ix,iaux) - else - call psb_ip_reord(n,x,iaux) - end if + end if + + allocate(iaux(0:n+1),stat=info) + if (info /= psb_success_) then + call psb_errpush(psb_err_alloc_dealloc_,r_name='psb_m_msort') + goto 9999 + endif + + select case(dir_) + case (psb_sort_up_) + call psi_m_msort_up(n,x,iaux,iret) + case (psb_sort_down_) + call psi_m_msort_dw(n,x,iaux,iret) + case (psb_asort_up_) + call psi_m_amsort_up(n,x,iaux,iret) + case (psb_asort_down_) + call psi_m_amsort_dw(n,x,iaux,iret) + end select + ! + ! Do the actual reordering, since the inner routines + ! only provide linked pointers. + ! + if (iret == 0 ) then + if (present(ix)) then + call psb_ip_reord(n,x,ix,iaux) + else + call psb_ip_reord(n,x,iaux) end if + end if - return + return 9999 call psb_error_handler(err_act) - return + return - end subroutine psb_mmsort - - subroutine psi_m_msort_up(n,k,l,iret) - use psb_const_mod - implicit none - integer(psb_ipk_) :: n, iret - integer(psb_mpk_) :: k(n) - integer(psb_ipk_) :: l(0:n+1) - ! - integer(psb_ipk_) :: p,q,s,t - ! .. - iret = 0 - ! first step: we are preparing ordered sublists, exploiting - ! what order was already in the input data; negative links - ! mark the end of the sublists - l(0) = 1 - t = n + 1 - do p = 1,n - 1 - if (k(p) <= k(p+1)) then - l(p) = p + 1 - else - l(t) = - (p+1) - t = p - end if - end do - l(t) = 0 - l(n) = 0 - ! see if the input was already sorted - if (l(n+1) == 0) then - iret = 1 - return +end subroutine psb_mmsort + +subroutine psi_m_msort_up(n,k,l,iret) + use psb_const_mod + implicit none + integer(psb_ipk_) :: n, iret + integer(psb_mpk_) :: k(n) + integer(psb_ipk_) :: l(0:n+1) + ! + integer(psb_ipk_) :: p,q,s,t + ! .. + iret = 0 + ! first step: we are preparing ordered sublists, exploiting + ! what order was already in the input data; negative links + ! mark the end of the sublists + l(0) = 1 + t = n + 1 + do p = 1,n - 1 + if (k(p) <= k(p+1)) then + l(p) = p + 1 else - l(n+1) = abs(l(n+1)) + l(t) = - (p+1) + t = p end if + end do + l(t) = 0 + l(n) = 0 + ! see if the input was already sorted + if (l(n+1) == 0) then + iret = 1 + return + else + l(n+1) = abs(l(n+1)) + end if + + mergepass: do + ! otherwise, begin a pass through the list. + ! throughout all the subroutine we have: + ! p, q: pointing to the sublists being merged + ! s: pointing to the most recently processed record + ! t: pointing to the end of previously completed sublist + s = 0 + t = n + 1 + p = l(s) + q = l(t) + if (q == 0) exit mergepass - mergepass: do - ! otherwise, begin a pass through the list. - ! throughout all the subroutine we have: - ! p, q: pointing to the sublists being merged - ! s: pointing to the most recently processed record - ! t: pointing to the end of previously completed sublist - s = 0 - t = n + 1 - p = l(s) - q = l(t) - if (q == 0) exit mergepass - - outer: do + outer: do - if (k(p) > k(q)) then + if (k(p) > k(q)) then - l(s) = sign(q,l(s)) - s = q - q = l(q) - if (q > 0) then - do - if (k(p) <= k(q)) cycle outer - s = q - q = l(q) - if (q <= 0) exit - end do - end if - l(s) = p - s = t + l(s) = sign(q,l(s)) + s = q + q = l(q) + if (q > 0) then do - t = p - p = l(p) - if (p <= 0) exit - end do - - else - - l(s) = sign(p,l(s)) - s = p - p = l(p) - if (p>0) then - do - if (k(p) > k(q)) cycle outer - s = p - p = l(p) - if (p <= 0) exit - end do - end if - ! otherwise, one sublist ended, and we append to it the rest - ! of the other one. - l(s) = q - s = t - do - t = q + if (k(p) <= k(q)) cycle outer + s = q q = l(q) if (q <= 0) exit end do end if + l(s) = p + s = t + do + t = p + p = l(p) + if (p <= 0) exit + end do - p = -p - q = -q - if (q == 0) then - l(s) = sign(p,l(s)) - l(t) = 0 - exit outer + else + + l(s) = sign(p,l(s)) + s = p + p = l(p) + if (p>0) then + do + if (k(p) > k(q)) cycle outer + s = p + p = l(p) + if (p <= 0) exit + end do end if - end do outer - end do mergepass - - end subroutine psi_m_msort_up - - subroutine psi_m_msort_dw(n,k,l,iret) - use psb_const_mod - implicit none - integer(psb_ipk_) :: n, iret - integer(psb_mpk_) :: k(n) - integer(psb_ipk_) :: l(0:n+1) - ! - integer(psb_ipk_) :: p,q,s,t - ! .. - iret = 0 - ! first step: we are preparing ordered sublists, exploiting - ! what order was already in the input data; negative links - ! mark the end of the sublists - l(0) = 1 - t = n + 1 - do p = 1,n - 1 - if (k(p) >= k(p+1)) then - l(p) = p + 1 - else - l(t) = - (p+1) - t = p + ! otherwise, one sublist ended, and we append to it the rest + ! of the other one. + l(s) = q + s = t + do + t = q + q = l(q) + if (q <= 0) exit + end do end if - end do - l(t) = 0 - l(n) = 0 - ! see if the input was already sorted - if (l(n+1) == 0) then - iret = 1 - return - else - l(n+1) = abs(l(n+1)) - end if - mergepass: do - ! otherwise, begin a pass through the list. - ! throughout all the subroutine we have: - ! p, q: pointing to the sublists being merged - ! s: pointing to the most recently processed record - ! t: pointing to the end of previously completed sublist - s = 0 - t = n + 1 - p = l(s) - q = l(t) - if (q == 0) exit mergepass + p = -p + q = -q + if (q == 0) then + l(s) = sign(p,l(s)) + l(t) = 0 + exit outer + end if + end do outer + end do mergepass - outer: do +end subroutine psi_m_msort_up - if (k(p) < k(q)) then +subroutine psi_m_msort_dw(n,k,l,iret) + use psb_const_mod + implicit none + integer(psb_ipk_) :: n, iret + integer(psb_mpk_) :: k(n) + integer(psb_ipk_) :: l(0:n+1) + ! + integer(psb_ipk_) :: p,q,s,t + ! .. + iret = 0 + ! first step: we are preparing ordered sublists, exploiting + ! what order was already in the input data; negative links + ! mark the end of the sublists + l(0) = 1 + t = n + 1 + do p = 1,n - 1 + if (k(p) >= k(p+1)) then + l(p) = p + 1 + else + l(t) = - (p+1) + t = p + end if + end do + l(t) = 0 + l(n) = 0 + ! see if the input was already sorted + if (l(n+1) == 0) then + iret = 1 + return + else + l(n+1) = abs(l(n+1)) + end if + + mergepass: do + ! otherwise, begin a pass through the list. + ! throughout all the subroutine we have: + ! p, q: pointing to the sublists being merged + ! s: pointing to the most recently processed record + ! t: pointing to the end of previously completed sublist + s = 0 + t = n + 1 + p = l(s) + q = l(t) + if (q == 0) exit mergepass - l(s) = sign(q,l(s)) - s = q - q = l(q) - if (q > 0) then - do - if (k(p) >= k(q)) cycle outer - s = q - q = l(q) - if (q <= 0) exit - end do - end if - l(s) = p - s = t - do - t = p - p = l(p) - if (p <= 0) exit - end do + outer: do - else + if (k(p) < k(q)) then - l(s) = sign(p,l(s)) - s = p - p = l(p) - if (p>0) then - do - if (k(p) < k(q)) cycle outer - s = p - p = l(p) - if (p <= 0) exit - end do - end if - ! otherwise, one sublist ended, and we append to it the rest - ! of the other one. - l(s) = q - s = t + l(s) = sign(q,l(s)) + s = q + q = l(q) + if (q > 0) then do - t = q + if (k(p) >= k(q)) cycle outer + s = q q = l(q) if (q <= 0) exit end do end if + l(s) = p + s = t + do + t = p + p = l(p) + if (p <= 0) exit + end do - p = -p - q = -q - if (q == 0) then - l(s) = sign(p,l(s)) - l(t) = 0 - exit outer + else + + l(s) = sign(p,l(s)) + s = p + p = l(p) + if (p>0) then + do + if (k(p) < k(q)) cycle outer + s = p + p = l(p) + if (p <= 0) exit + end do end if - end do outer - end do mergepass - - end subroutine psi_m_msort_dw - - subroutine psi_m_amsort_up(n,k,l,iret) - use psb_const_mod - implicit none - integer(psb_ipk_) :: n, iret - integer(psb_mpk_) :: k(n) - integer(psb_ipk_) :: l(0:n+1) - ! - integer(psb_ipk_) :: p,q,s,t - ! .. - iret = 0 - ! first step: we are preparing ordered sublists, exploiting - ! what order was already in the input data; negative links - ! mark the end of the sublists - l(0) = 1 - t = n + 1 - do p = 1,n - 1 - if (abs(k(p)) <= abs(k(p+1))) then - l(p) = p + 1 - else - l(t) = - (p+1) - t = p + ! otherwise, one sublist ended, and we append to it the rest + ! of the other one. + l(s) = q + s = t + do + t = q + q = l(q) + if (q <= 0) exit + end do end if - end do - l(t) = 0 - l(n) = 0 - ! see if the input was already sorted - if (l(n+1) == 0) then - iret = 1 - return - else - l(n+1) = abs(l(n+1)) - end if - mergepass: do - ! otherwise, begin a pass through the list. - ! throughout all the subroutine we have: - ! p, q: pointing to the sublists being merged - ! s: pointing to the most recently processed record - ! t: pointing to the end of previously completed sublist - s = 0 - t = n + 1 - p = l(s) - q = l(t) - if (q == 0) exit mergepass + p = -p + q = -q + if (q == 0) then + l(s) = sign(p,l(s)) + l(t) = 0 + exit outer + end if + end do outer + end do mergepass - outer: do +end subroutine psi_m_msort_dw - if (abs(k(p)) > abs(k(q))) then +subroutine psi_m_amsort_up(n,k,l,iret) + use psb_const_mod + implicit none + integer(psb_ipk_) :: n, iret + integer(psb_mpk_) :: k(n) + integer(psb_ipk_) :: l(0:n+1) + ! + integer(psb_ipk_) :: p,q,s,t + ! .. + iret = 0 + ! first step: we are preparing ordered sublists, exploiting + ! what order was already in the input data; negative links + ! mark the end of the sublists + l(0) = 1 + t = n + 1 + do p = 1,n - 1 + if (abs(k(p)) <= abs(k(p+1))) then + l(p) = p + 1 + else + l(t) = - (p+1) + t = p + end if + end do + l(t) = 0 + l(n) = 0 + ! see if the input was already sorted + if (l(n+1) == 0) then + iret = 1 + return + else + l(n+1) = abs(l(n+1)) + end if + + mergepass: do + ! otherwise, begin a pass through the list. + ! throughout all the subroutine we have: + ! p, q: pointing to the sublists being merged + ! s: pointing to the most recently processed record + ! t: pointing to the end of previously completed sublist + s = 0 + t = n + 1 + p = l(s) + q = l(t) + if (q == 0) exit mergepass - l(s) = sign(q,l(s)) - s = q - q = l(q) - if (q > 0) then - do - if (abs(k(p)) <= abs(k(q))) cycle outer - s = q - q = l(q) - if (q <= 0) exit - end do - end if - l(s) = p - s = t - do - t = p - p = l(p) - if (p <= 0) exit - end do + outer: do - else + if (abs(k(p)) > abs(k(q))) then - l(s) = sign(p,l(s)) - s = p - p = l(p) - if (p>0) then - do - if (abs(k(p)) > abs(k(q))) cycle outer - s = p - p = l(p) - if (p <= 0) exit - end do - end if - ! otherwise, one sublist ended, and we append to it the rest - ! of the other one. - l(s) = q - s = t + l(s) = sign(q,l(s)) + s = q + q = l(q) + if (q > 0) then do - t = q + if (abs(k(p)) <= abs(k(q))) cycle outer + s = q q = l(q) if (q <= 0) exit end do end if + l(s) = p + s = t + do + t = p + p = l(p) + if (p <= 0) exit + end do - p = -p - q = -q - if (q == 0) then - l(s) = sign(p,l(s)) - l(t) = 0 - exit outer + else + + l(s) = sign(p,l(s)) + s = p + p = l(p) + if (p>0) then + do + if (abs(k(p)) > abs(k(q))) cycle outer + s = p + p = l(p) + if (p <= 0) exit + end do end if - end do outer - end do mergepass - - end subroutine psi_m_amsort_up - - subroutine psi_m_amsort_dw(n,k,l,iret) - use psb_const_mod - implicit none - integer(psb_ipk_) :: n, iret - integer(psb_mpk_) :: k(n) - integer(psb_ipk_) :: l(0:n+1) - ! - integer(psb_ipk_) :: p,q,s,t - ! .. - iret = 0 - ! first step: we are preparing ordered sublists, exploiting - ! what order was already in the input data; negative links - ! mark the end of the sublists - l(0) = 1 - t = n + 1 - do p = 1,n - 1 - if (abs(k(p)) >= abs(k(p+1))) then - l(p) = p + 1 - else - l(t) = - (p+1) - t = p + ! otherwise, one sublist ended, and we append to it the rest + ! of the other one. + l(s) = q + s = t + do + t = q + q = l(q) + if (q <= 0) exit + end do end if - end do - l(t) = 0 - l(n) = 0 - ! see if the input was already sorted - if (l(n+1) == 0) then - iret = 1 - return - else - l(n+1) = abs(l(n+1)) - end if - mergepass: do - ! otherwise, begin a pass through the list. - ! throughout all the subroutine we have: - ! p, q: pointing to the sublists being merged - ! s: pointing to the most recently processed record - ! t: pointing to the end of previously completed sublist - s = 0 - t = n + 1 - p = l(s) - q = l(t) - if (q == 0) exit mergepass + p = -p + q = -q + if (q == 0) then + l(s) = sign(p,l(s)) + l(t) = 0 + exit outer + end if + end do outer + end do mergepass - outer: do +end subroutine psi_m_amsort_up - if (abs(k(p)) < abs(k(q))) then +subroutine psi_m_amsort_dw(n,k,l,iret) + use psb_const_mod + implicit none + integer(psb_ipk_) :: n, iret + integer(psb_mpk_) :: k(n) + integer(psb_ipk_) :: l(0:n+1) + ! + integer(psb_ipk_) :: p,q,s,t + ! .. + iret = 0 + ! first step: we are preparing ordered sublists, exploiting + ! what order was already in the input data; negative links + ! mark the end of the sublists + l(0) = 1 + t = n + 1 + do p = 1,n - 1 + if (abs(k(p)) >= abs(k(p+1))) then + l(p) = p + 1 + else + l(t) = - (p+1) + t = p + end if + end do + l(t) = 0 + l(n) = 0 + ! see if the input was already sorted + if (l(n+1) == 0) then + iret = 1 + return + else + l(n+1) = abs(l(n+1)) + end if + + mergepass: do + ! otherwise, begin a pass through the list. + ! throughout all the subroutine we have: + ! p, q: pointing to the sublists being merged + ! s: pointing to the most recently processed record + ! t: pointing to the end of previously completed sublist + s = 0 + t = n + 1 + p = l(s) + q = l(t) + if (q == 0) exit mergepass - l(s) = sign(q,l(s)) - s = q - q = l(q) - if (q > 0) then - do - if (abs(k(p)) >= abs(k(q))) cycle outer - s = q - q = l(q) - if (q <= 0) exit - end do - end if - l(s) = p - s = t - do - t = p - p = l(p) - if (p <= 0) exit - end do + outer: do - else + if (abs(k(p)) < abs(k(q))) then - l(s) = sign(p,l(s)) - s = p - p = l(p) - if (p>0) then - do - if (abs(k(p)) < abs(k(q))) cycle outer - s = p - p = l(p) - if (p <= 0) exit - end do - end if - ! otherwise, one sublist ended, and we append to it the rest - ! of the other one. - l(s) = q - s = t + l(s) = sign(q,l(s)) + s = q + q = l(q) + if (q > 0) then do - t = q + if (abs(k(p)) >= abs(k(q))) cycle outer + s = q q = l(q) if (q <= 0) exit end do end if + l(s) = p + s = t + do + t = p + p = l(p) + if (p <= 0) exit + end do - p = -p - q = -q - if (q == 0) then - l(s) = sign(p,l(s)) - l(t) = 0 - exit outer - end if - end do outer - end do mergepass - - end subroutine psi_m_amsort_dw - - - + else + l(s) = sign(p,l(s)) + s = p + p = l(p) + if (p>0) then + do + if (abs(k(p)) < abs(k(q))) cycle outer + s = p + p = l(p) + if (p <= 0) exit + end do + end if + ! otherwise, one sublist ended, and we append to it the rest + ! of the other one. + l(s) = q + s = t + do + t = q + q = l(q) + if (q <= 0) exit + end do + end if + p = -p + q = -q + if (q == 0) then + l(s) = sign(p,l(s)) + l(t) = 0 + exit outer + end if + end do outer + end do mergepass +end subroutine psi_m_amsort_dw diff --git a/base/serial/sort/psb_s_hsort_impl.f90 b/base/serial/sort/psb_s_hsort_impl.f90 index 4737d159..77fefe14 100644 --- a/base/serial/sort/psb_s_hsort_impl.f90 +++ b/base/serial/sort/psb_s_hsort_impl.f90 @@ -49,7 +49,8 @@ subroutine psb_shsort(x,ix,dir,flag) integer(psb_ipk_), optional, intent(in) :: dir, flag integer(psb_ipk_), optional, intent(inout) :: ix(:) - integer(psb_ipk_) :: dir_, flag_, n, i, l, err_act,info + integer(psb_ipk_) :: flag_, n, i, err_act,info + integer(psb_ipk_) :: dir_, l real(psb_spk_) :: key integer(psb_ipk_) :: index @@ -159,7 +160,7 @@ end subroutine psb_shsort ! ! These are packaged so that they can be used to implement -! a heapsort, should the need arise +! a heapsort. ! ! ! Programming note: @@ -540,11 +541,12 @@ subroutine psi_s_idx_heap_get_first(key,index,last,heap,idxs,dir,info) use psb_sort_mod, psb_protect_name => psi_s_idx_heap_get_first implicit none + real(psb_spk_), intent(inout) :: key real(psb_spk_), intent(inout) :: heap(:) - integer(psb_ipk_), intent(out) :: index,info + integer(psb_ipk_), intent(out) :: index + integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(inout) :: last,idxs(:) integer(psb_ipk_), intent(in) :: dir - real(psb_spk_), intent(out) :: key integer(psb_ipk_) :: i, j,itemp real(psb_spk_) :: temp diff --git a/base/serial/sort/psb_s_msort_impl.f90 b/base/serial/sort/psb_s_msort_impl.f90 index a1af2f56..dfd7508c 100644 --- a/base/serial/sort/psb_s_msort_impl.f90 +++ b/base/serial/sort/psb_s_msort_impl.f90 @@ -1,34 +1,34 @@ -! -! 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. -! -! + ! + ! 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. + ! + ! ! ! The merge-sort routines ! References: @@ -41,618 +41,612 @@ ! Addison-Wesley ! - subroutine psb_smsort_u(x,nout,dir) - use psb_sort_mod, psb_protect_name => psb_smsort_u - use psb_error_mod - implicit none - real(psb_spk_), intent(inout) :: x(:) - integer(psb_ipk_), intent(out) :: nout - integer(psb_ipk_), optional, intent(in) :: dir +subroutine psb_smsort_u(x,nout,dir) + use psb_sort_mod, psb_protect_name => psb_smsort_u + use psb_error_mod + implicit none + real(psb_spk_), intent(inout) :: x(:) + integer(psb_ipk_), intent(out) :: nout + integer(psb_ipk_), optional, intent(in) :: dir - integer(psb_ipk_) :: n, k - integer(psb_ipk_) :: err_act + integer(psb_ipk_) :: n, k + integer(psb_ipk_) :: err_act - integer(psb_ipk_) :: ierr(5) - character(len=20) :: name + integer(psb_ipk_) :: ierr(5) + character(len=20) :: name - name='psb_msort_u' - call psb_erractionsave(err_act) + name='psb_msort_u' + call psb_erractionsave(err_act) - n = size(x) + n = size(x) - call psb_msort(x,dir=dir) - nout = min(1,n) - do k=2,n - if (x(k) /= x(nout)) then - nout = nout + 1 - x(nout) = x(k) - endif - enddo + call psb_msort(x,dir=dir) + nout = min(1,n) + do k=2,n + if (x(k) /= x(nout)) then + nout = nout + 1 + x(nout) = x(k) + endif + enddo - return + return 9999 call psb_error_handler(err_act) - return - end subroutine psb_smsort_u + return +end subroutine psb_smsort_u - function psb_sbsrch(key,n,v) result(ipos) - use psb_sort_mod, psb_protect_name => psb_sbsrch - implicit none - integer(psb_ipk_) :: ipos, n - real(psb_spk_) :: key - real(psb_spk_) :: v(:) +function psb_sbsrch(key,n,v) result(ipos) + use psb_sort_mod, psb_protect_name => psb_sbsrch + implicit none + integer(psb_ipk_) :: ipos, n + real(psb_spk_) :: key + real(psb_spk_) :: v(:) - integer(psb_ipk_) :: lb, ub, m, i + integer(psb_ipk_) :: lb, ub, m, i - ipos = -1 - if (n<5) then - do i=1,n - if (key.eq.v(i)) then - ipos = i - return - end if - enddo - return - end if - - lb = 1 - ub = n - - do while (lb.le.ub) - m = (lb+ub)/2 - if (key.eq.v(m)) then - ipos = m - lb = ub + 1 - else if (key < v(m)) then - ub = m-1 - else - lb = m + 1 - end if - enddo - return - end function psb_sbsrch - - function psb_sssrch(key,n,v) result(ipos) - use psb_sort_mod, psb_protect_name => psb_sssrch - implicit none - integer(psb_ipk_) :: ipos, n - real(psb_spk_) :: key - real(psb_spk_) :: v(:) - - integer(psb_ipk_) :: i - - ipos = -1 + ipos = -1 + if (n<5) then do i=1,n if (key.eq.v(i)) then ipos = i return end if enddo - return - end function psb_sssrch - - subroutine psb_smsort(x,ix,dir,flag) - use psb_sort_mod, psb_protect_name => psb_smsort - use psb_error_mod - use psb_ip_reord_mod - implicit none - real(psb_spk_), intent(inout) :: x(:) - integer(psb_ipk_), optional, intent(in) :: dir, flag - integer(psb_ipk_), optional, intent(inout) :: ix(:) - - integer(psb_ipk_) :: dir_, flag_, n, err_act - - integer(psb_ipk_), allocatable :: iaux(:) - integer(psb_ipk_) :: iret, info, i - integer(psb_ipk_) :: ierr(5) - character(len=20) :: name - - name='psb_smsort' - call psb_erractionsave(err_act) - - if (present(dir)) then - dir_ = dir - else - dir_= psb_sort_up_ + end if + + lb = 1 + ub = n + + do while (lb.le.ub) + m = (lb+ub)/2 + if (key.eq.v(m)) then + ipos = m + lb = ub + 1 + else if (key < v(m)) then + ub = m-1 + else + lb = m + 1 end if - select case(dir_) - case( psb_sort_up_, psb_sort_down_, psb_asort_up_, psb_asort_down_) + enddo + return +end function psb_sbsrch + +function psb_sssrch(key,n,v) result(ipos) + use psb_sort_mod, psb_protect_name => psb_sssrch + implicit none + integer(psb_ipk_) :: ipos, n + real(psb_spk_) :: key + real(psb_spk_) :: v(:) + + integer(psb_ipk_) :: i + + ipos = -1 + do i=1,n + if (key.eq.v(i)) then + ipos = i + return + end if + enddo + + return +end function psb_sssrch + +subroutine psb_smsort(x,ix,dir,flag) + use psb_sort_mod, psb_protect_name => psb_smsort + use psb_error_mod + use psb_ip_reord_mod + implicit none + real(psb_spk_), intent(inout) :: x(:) + integer(psb_ipk_), optional, intent(in) :: dir, flag + integer(psb_ipk_), optional, intent(inout) :: ix(:) + + integer(psb_ipk_) :: dir_, flag_, n, err_act + + integer(psb_ipk_), allocatable :: iaux(:) + integer(psb_ipk_) :: iret, info, i + integer(psb_ipk_) :: ierr(5) + character(len=20) :: name + + name='psb_smsort' + call psb_erractionsave(err_act) + + if (present(dir)) then + dir_ = dir + else + dir_= psb_sort_up_ + end if + select case(dir_) + case( psb_sort_up_, psb_sort_down_, psb_asort_up_, psb_asort_down_) + ! OK keep going + case default + ierr(1) = 3; ierr(2) = dir_; + call psb_errpush(psb_err_input_value_invalid_i_,name,i_err=ierr) + goto 9999 + end select + + n = size(x) + + if (present(ix)) then + if (size(ix) < n) then + ierr(1) = 2; ierr(2) = size(ix); + call psb_errpush(psb_err_input_asize_invalid_i_,name,i_err=ierr) + goto 9999 + end if + if (present(flag)) then + flag_ = flag + else + flag_ = psb_sort_ovw_idx_ + end if + select case(flag_) + case(psb_sort_ovw_idx_) + do i=1,n + ix(i) = i + end do + case (psb_sort_keep_idx_) ! OK keep going case default - ierr(1) = 3; ierr(2) = dir_; + ierr(1) = 4; ierr(2) = flag_; call psb_errpush(psb_err_input_value_invalid_i_,name,i_err=ierr) goto 9999 end select - n = size(x) - - if (present(ix)) then - if (size(ix) < n) then - ierr(1) = 2; ierr(2) = size(ix); - call psb_errpush(psb_err_input_asize_invalid_i_,name,i_err=ierr) - goto 9999 - end if - if (present(flag)) then - flag_ = flag - else - flag_ = psb_sort_ovw_idx_ - end if - select case(flag_) - case(psb_sort_ovw_idx_) - do i=1,n - ix(i) = i - end do - case (psb_sort_keep_idx_) - ! OK keep going - case default - ierr(1) = 4; ierr(2) = flag_; - call psb_errpush(psb_err_input_value_invalid_i_,name,i_err=ierr) - goto 9999 - end select - - end if - - allocate(iaux(0:n+1),stat=info) - if (info /= psb_success_) then - call psb_errpush(psb_err_alloc_dealloc_,r_name='psb_s_msort') - goto 9999 - endif - - select case(dir_) - case (psb_sort_up_) - call psi_s_msort_up(n,x,iaux,iret) - case (psb_sort_down_) - call psi_s_msort_dw(n,x,iaux,iret) - case (psb_asort_up_) - call psi_s_amsort_up(n,x,iaux,iret) - case (psb_asort_down_) - call psi_s_amsort_dw(n,x,iaux,iret) - end select - ! - ! Do the actual reordering, since the inner routines - ! only provide linked pointers. - ! - if (iret == 0 ) then - if (present(ix)) then - call psb_ip_reord(n,x,ix,iaux) - else - call psb_ip_reord(n,x,iaux) - end if + end if + + allocate(iaux(0:n+1),stat=info) + if (info /= psb_success_) then + call psb_errpush(psb_err_alloc_dealloc_,r_name='psb_s_msort') + goto 9999 + endif + + select case(dir_) + case (psb_sort_up_) + call psi_s_msort_up(n,x,iaux,iret) + case (psb_sort_down_) + call psi_s_msort_dw(n,x,iaux,iret) + case (psb_asort_up_) + call psi_s_amsort_up(n,x,iaux,iret) + case (psb_asort_down_) + call psi_s_amsort_dw(n,x,iaux,iret) + end select + ! + ! Do the actual reordering, since the inner routines + ! only provide linked pointers. + ! + if (iret == 0 ) then + if (present(ix)) then + call psb_ip_reord(n,x,ix,iaux) + else + call psb_ip_reord(n,x,iaux) end if + end if - return + return 9999 call psb_error_handler(err_act) - return + return - end subroutine psb_smsort - - subroutine psi_s_msort_up(n,k,l,iret) - use psb_const_mod - implicit none - integer(psb_ipk_) :: n, iret - real(psb_spk_) :: k(n) - integer(psb_ipk_) :: l(0:n+1) - ! - integer(psb_ipk_) :: p,q,s,t - ! .. - iret = 0 - ! first step: we are preparing ordered sublists, exploiting - ! what order was already in the input data; negative links - ! mark the end of the sublists - l(0) = 1 - t = n + 1 - do p = 1,n - 1 - if (k(p) <= k(p+1)) then - l(p) = p + 1 - else - l(t) = - (p+1) - t = p - end if - end do - l(t) = 0 - l(n) = 0 - ! see if the input was already sorted - if (l(n+1) == 0) then - iret = 1 - return +end subroutine psb_smsort + +subroutine psi_s_msort_up(n,k,l,iret) + use psb_const_mod + implicit none + integer(psb_ipk_) :: n, iret + real(psb_spk_) :: k(n) + integer(psb_ipk_) :: l(0:n+1) + ! + integer(psb_ipk_) :: p,q,s,t + ! .. + iret = 0 + ! first step: we are preparing ordered sublists, exploiting + ! what order was already in the input data; negative links + ! mark the end of the sublists + l(0) = 1 + t = n + 1 + do p = 1,n - 1 + if (k(p) <= k(p+1)) then + l(p) = p + 1 else - l(n+1) = abs(l(n+1)) + l(t) = - (p+1) + t = p end if + end do + l(t) = 0 + l(n) = 0 + ! see if the input was already sorted + if (l(n+1) == 0) then + iret = 1 + return + else + l(n+1) = abs(l(n+1)) + end if + + mergepass: do + ! otherwise, begin a pass through the list. + ! throughout all the subroutine we have: + ! p, q: pointing to the sublists being merged + ! s: pointing to the most recently processed record + ! t: pointing to the end of previously completed sublist + s = 0 + t = n + 1 + p = l(s) + q = l(t) + if (q == 0) exit mergepass - mergepass: do - ! otherwise, begin a pass through the list. - ! throughout all the subroutine we have: - ! p, q: pointing to the sublists being merged - ! s: pointing to the most recently processed record - ! t: pointing to the end of previously completed sublist - s = 0 - t = n + 1 - p = l(s) - q = l(t) - if (q == 0) exit mergepass - - outer: do - - if (k(p) > k(q)) then - - l(s) = sign(q,l(s)) - s = q - q = l(q) - if (q > 0) then - do - if (k(p) <= k(q)) cycle outer - s = q - q = l(q) - if (q <= 0) exit - end do - end if - l(s) = p - s = t - do - t = p - p = l(p) - if (p <= 0) exit - end do + outer: do - else + if (k(p) > k(q)) then - l(s) = sign(p,l(s)) - s = p - p = l(p) - if (p>0) then - do - if (k(p) > k(q)) cycle outer - s = p - p = l(p) - if (p <= 0) exit - end do - end if - ! otherwise, one sublist ended, and we append to it the rest - ! of the other one. - l(s) = q - s = t + l(s) = sign(q,l(s)) + s = q + q = l(q) + if (q > 0) then do - t = q + if (k(p) <= k(q)) cycle outer + s = q q = l(q) if (q <= 0) exit end do end if + l(s) = p + s = t + do + t = p + p = l(p) + if (p <= 0) exit + end do - p = -p - q = -q - if (q == 0) then - l(s) = sign(p,l(s)) - l(t) = 0 - exit outer + else + + l(s) = sign(p,l(s)) + s = p + p = l(p) + if (p>0) then + do + if (k(p) > k(q)) cycle outer + s = p + p = l(p) + if (p <= 0) exit + end do end if - end do outer - end do mergepass - - end subroutine psi_s_msort_up - - subroutine psi_s_msort_dw(n,k,l,iret) - use psb_const_mod - implicit none - integer(psb_ipk_) :: n, iret - real(psb_spk_) :: k(n) - integer(psb_ipk_) :: l(0:n+1) - ! - integer(psb_ipk_) :: p,q,s,t - ! .. - iret = 0 - ! first step: we are preparing ordered sublists, exploiting - ! what order was already in the input data; negative links - ! mark the end of the sublists - l(0) = 1 - t = n + 1 - do p = 1,n - 1 - if (k(p) >= k(p+1)) then - l(p) = p + 1 - else - l(t) = - (p+1) - t = p + ! otherwise, one sublist ended, and we append to it the rest + ! of the other one. + l(s) = q + s = t + do + t = q + q = l(q) + if (q <= 0) exit + end do end if - end do - l(t) = 0 - l(n) = 0 - ! see if the input was already sorted - if (l(n+1) == 0) then - iret = 1 - return - else - l(n+1) = abs(l(n+1)) - end if - mergepass: do - ! otherwise, begin a pass through the list. - ! throughout all the subroutine we have: - ! p, q: pointing to the sublists being merged - ! s: pointing to the most recently processed record - ! t: pointing to the end of previously completed sublist - s = 0 - t = n + 1 - p = l(s) - q = l(t) - if (q == 0) exit mergepass + p = -p + q = -q + if (q == 0) then + l(s) = sign(p,l(s)) + l(t) = 0 + exit outer + end if + end do outer + end do mergepass - outer: do +end subroutine psi_s_msort_up - if (k(p) < k(q)) then +subroutine psi_s_msort_dw(n,k,l,iret) + use psb_const_mod + implicit none + integer(psb_ipk_) :: n, iret + real(psb_spk_) :: k(n) + integer(psb_ipk_) :: l(0:n+1) + ! + integer(psb_ipk_) :: p,q,s,t + ! .. + iret = 0 + ! first step: we are preparing ordered sublists, exploiting + ! what order was already in the input data; negative links + ! mark the end of the sublists + l(0) = 1 + t = n + 1 + do p = 1,n - 1 + if (k(p) >= k(p+1)) then + l(p) = p + 1 + else + l(t) = - (p+1) + t = p + end if + end do + l(t) = 0 + l(n) = 0 + ! see if the input was already sorted + if (l(n+1) == 0) then + iret = 1 + return + else + l(n+1) = abs(l(n+1)) + end if + + mergepass: do + ! otherwise, begin a pass through the list. + ! throughout all the subroutine we have: + ! p, q: pointing to the sublists being merged + ! s: pointing to the most recently processed record + ! t: pointing to the end of previously completed sublist + s = 0 + t = n + 1 + p = l(s) + q = l(t) + if (q == 0) exit mergepass - l(s) = sign(q,l(s)) - s = q - q = l(q) - if (q > 0) then - do - if (k(p) >= k(q)) cycle outer - s = q - q = l(q) - if (q <= 0) exit - end do - end if - l(s) = p - s = t - do - t = p - p = l(p) - if (p <= 0) exit - end do + outer: do - else + if (k(p) < k(q)) then - l(s) = sign(p,l(s)) - s = p - p = l(p) - if (p>0) then - do - if (k(p) < k(q)) cycle outer - s = p - p = l(p) - if (p <= 0) exit - end do - end if - ! otherwise, one sublist ended, and we append to it the rest - ! of the other one. - l(s) = q - s = t + l(s) = sign(q,l(s)) + s = q + q = l(q) + if (q > 0) then do - t = q + if (k(p) >= k(q)) cycle outer + s = q q = l(q) if (q <= 0) exit end do end if + l(s) = p + s = t + do + t = p + p = l(p) + if (p <= 0) exit + end do - p = -p - q = -q - if (q == 0) then - l(s) = sign(p,l(s)) - l(t) = 0 - exit outer + else + + l(s) = sign(p,l(s)) + s = p + p = l(p) + if (p>0) then + do + if (k(p) < k(q)) cycle outer + s = p + p = l(p) + if (p <= 0) exit + end do end if - end do outer - end do mergepass - - end subroutine psi_s_msort_dw - - subroutine psi_s_amsort_up(n,k,l,iret) - use psb_const_mod - implicit none - integer(psb_ipk_) :: n, iret - real(psb_spk_) :: k(n) - integer(psb_ipk_) :: l(0:n+1) - ! - integer(psb_ipk_) :: p,q,s,t - ! .. - iret = 0 - ! first step: we are preparing ordered sublists, exploiting - ! what order was already in the input data; negative links - ! mark the end of the sublists - l(0) = 1 - t = n + 1 - do p = 1,n - 1 - if (abs(k(p)) <= abs(k(p+1))) then - l(p) = p + 1 - else - l(t) = - (p+1) - t = p + ! otherwise, one sublist ended, and we append to it the rest + ! of the other one. + l(s) = q + s = t + do + t = q + q = l(q) + if (q <= 0) exit + end do end if - end do - l(t) = 0 - l(n) = 0 - ! see if the input was already sorted - if (l(n+1) == 0) then - iret = 1 - return - else - l(n+1) = abs(l(n+1)) - end if - mergepass: do - ! otherwise, begin a pass through the list. - ! throughout all the subroutine we have: - ! p, q: pointing to the sublists being merged - ! s: pointing to the most recently processed record - ! t: pointing to the end of previously completed sublist - s = 0 - t = n + 1 - p = l(s) - q = l(t) - if (q == 0) exit mergepass + p = -p + q = -q + if (q == 0) then + l(s) = sign(p,l(s)) + l(t) = 0 + exit outer + end if + end do outer + end do mergepass - outer: do +end subroutine psi_s_msort_dw - if (abs(k(p)) > abs(k(q))) then +subroutine psi_s_amsort_up(n,k,l,iret) + use psb_const_mod + implicit none + integer(psb_ipk_) :: n, iret + real(psb_spk_) :: k(n) + integer(psb_ipk_) :: l(0:n+1) + ! + integer(psb_ipk_) :: p,q,s,t + ! .. + iret = 0 + ! first step: we are preparing ordered sublists, exploiting + ! what order was already in the input data; negative links + ! mark the end of the sublists + l(0) = 1 + t = n + 1 + do p = 1,n - 1 + if (abs(k(p)) <= abs(k(p+1))) then + l(p) = p + 1 + else + l(t) = - (p+1) + t = p + end if + end do + l(t) = 0 + l(n) = 0 + ! see if the input was already sorted + if (l(n+1) == 0) then + iret = 1 + return + else + l(n+1) = abs(l(n+1)) + end if + + mergepass: do + ! otherwise, begin a pass through the list. + ! throughout all the subroutine we have: + ! p, q: pointing to the sublists being merged + ! s: pointing to the most recently processed record + ! t: pointing to the end of previously completed sublist + s = 0 + t = n + 1 + p = l(s) + q = l(t) + if (q == 0) exit mergepass - l(s) = sign(q,l(s)) - s = q - q = l(q) - if (q > 0) then - do - if (abs(k(p)) <= abs(k(q))) cycle outer - s = q - q = l(q) - if (q <= 0) exit - end do - end if - l(s) = p - s = t - do - t = p - p = l(p) - if (p <= 0) exit - end do + outer: do - else + if (abs(k(p)) > abs(k(q))) then - l(s) = sign(p,l(s)) - s = p - p = l(p) - if (p>0) then - do - if (abs(k(p)) > abs(k(q))) cycle outer - s = p - p = l(p) - if (p <= 0) exit - end do - end if - ! otherwise, one sublist ended, and we append to it the rest - ! of the other one. - l(s) = q - s = t + l(s) = sign(q,l(s)) + s = q + q = l(q) + if (q > 0) then do - t = q + if (abs(k(p)) <= abs(k(q))) cycle outer + s = q q = l(q) if (q <= 0) exit end do end if + l(s) = p + s = t + do + t = p + p = l(p) + if (p <= 0) exit + end do + + else - p = -p - q = -q - if (q == 0) then - l(s) = sign(p,l(s)) - l(t) = 0 - exit outer + l(s) = sign(p,l(s)) + s = p + p = l(p) + if (p>0) then + do + if (abs(k(p)) > abs(k(q))) cycle outer + s = p + p = l(p) + if (p <= 0) exit + end do end if - end do outer - end do mergepass - - end subroutine psi_s_amsort_up - - subroutine psi_s_amsort_dw(n,k,l,iret) - use psb_const_mod - implicit none - integer(psb_ipk_) :: n, iret - real(psb_spk_) :: k(n) - integer(psb_ipk_) :: l(0:n+1) - ! - integer(psb_ipk_) :: p,q,s,t - ! .. - iret = 0 - ! first step: we are preparing ordered sublists, exploiting - ! what order was already in the input data; negative links - ! mark the end of the sublists - l(0) = 1 - t = n + 1 - do p = 1,n - 1 - if (abs(k(p)) >= abs(k(p+1))) then - l(p) = p + 1 - else - l(t) = - (p+1) - t = p + ! otherwise, one sublist ended, and we append to it the rest + ! of the other one. + l(s) = q + s = t + do + t = q + q = l(q) + if (q <= 0) exit + end do end if - end do - l(t) = 0 - l(n) = 0 - ! see if the input was already sorted - if (l(n+1) == 0) then - iret = 1 - return - else - l(n+1) = abs(l(n+1)) - end if - mergepass: do - ! otherwise, begin a pass through the list. - ! throughout all the subroutine we have: - ! p, q: pointing to the sublists being merged - ! s: pointing to the most recently processed record - ! t: pointing to the end of previously completed sublist - s = 0 - t = n + 1 - p = l(s) - q = l(t) - if (q == 0) exit mergepass + p = -p + q = -q + if (q == 0) then + l(s) = sign(p,l(s)) + l(t) = 0 + exit outer + end if + end do outer + end do mergepass - outer: do +end subroutine psi_s_amsort_up - if (abs(k(p)) < abs(k(q))) then +subroutine psi_s_amsort_dw(n,k,l,iret) + use psb_const_mod + implicit none + integer(psb_ipk_) :: n, iret + real(psb_spk_) :: k(n) + integer(psb_ipk_) :: l(0:n+1) + ! + integer(psb_ipk_) :: p,q,s,t + ! .. + iret = 0 + ! first step: we are preparing ordered sublists, exploiting + ! what order was already in the input data; negative links + ! mark the end of the sublists + l(0) = 1 + t = n + 1 + do p = 1,n - 1 + if (abs(k(p)) >= abs(k(p+1))) then + l(p) = p + 1 + else + l(t) = - (p+1) + t = p + end if + end do + l(t) = 0 + l(n) = 0 + ! see if the input was already sorted + if (l(n+1) == 0) then + iret = 1 + return + else + l(n+1) = abs(l(n+1)) + end if + + mergepass: do + ! otherwise, begin a pass through the list. + ! throughout all the subroutine we have: + ! p, q: pointing to the sublists being merged + ! s: pointing to the most recently processed record + ! t: pointing to the end of previously completed sublist + s = 0 + t = n + 1 + p = l(s) + q = l(t) + if (q == 0) exit mergepass - l(s) = sign(q,l(s)) - s = q - q = l(q) - if (q > 0) then - do - if (abs(k(p)) >= abs(k(q))) cycle outer - s = q - q = l(q) - if (q <= 0) exit - end do - end if - l(s) = p - s = t - do - t = p - p = l(p) - if (p <= 0) exit - end do + outer: do - else + if (abs(k(p)) < abs(k(q))) then - l(s) = sign(p,l(s)) - s = p - p = l(p) - if (p>0) then - do - if (abs(k(p)) < abs(k(q))) cycle outer - s = p - p = l(p) - if (p <= 0) exit - end do - end if - ! otherwise, one sublist ended, and we append to it the rest - ! of the other one. - l(s) = q - s = t + l(s) = sign(q,l(s)) + s = q + q = l(q) + if (q > 0) then do - t = q + if (abs(k(p)) >= abs(k(q))) cycle outer + s = q q = l(q) if (q <= 0) exit end do end if + l(s) = p + s = t + do + t = p + p = l(p) + if (p <= 0) exit + end do - p = -p - q = -q - if (q == 0) then - l(s) = sign(p,l(s)) - l(t) = 0 - exit outer - end if - end do outer - end do mergepass - - end subroutine psi_s_amsort_dw - - - + else + l(s) = sign(p,l(s)) + s = p + p = l(p) + if (p>0) then + do + if (abs(k(p)) < abs(k(q))) cycle outer + s = p + p = l(p) + if (p <= 0) exit + end do + end if + ! otherwise, one sublist ended, and we append to it the rest + ! of the other one. + l(s) = q + s = t + do + t = q + q = l(q) + if (q <= 0) exit + end do + end if + p = -p + q = -q + if (q == 0) then + l(s) = sign(p,l(s)) + l(t) = 0 + exit outer + end if + end do outer + end do mergepass +end subroutine psi_s_amsort_dw diff --git a/base/serial/sort/psb_z_hsort_impl.f90 b/base/serial/sort/psb_z_hsort_impl.f90 index 7223e211..199f5663 100644 --- a/base/serial/sort/psb_z_hsort_impl.f90 +++ b/base/serial/sort/psb_z_hsort_impl.f90 @@ -49,7 +49,8 @@ subroutine psb_zhsort(x,ix,dir,flag) integer(psb_ipk_), optional, intent(in) :: dir, flag integer(psb_ipk_), optional, intent(inout) :: ix(:) - integer(psb_ipk_) :: dir_, flag_, n, i, l, err_act,info + integer(psb_ipk_) :: flag_, n, i, err_act,info + integer(psb_ipk_) :: dir_, l complex(psb_dpk_) :: key integer(psb_ipk_) :: index @@ -159,7 +160,7 @@ end subroutine psb_zhsort ! ! These are packaged so that they can be used to implement -! a heapsort, should the need arise +! a heapsort. ! ! ! Programming note: @@ -646,7 +647,8 @@ subroutine psi_z_idx_insert_heap(key,index,last,heap,idxs,dir,info) ! dir: sorting direction complex(psb_dpk_), intent(in) :: key - integer(psb_ipk_), intent(in) :: index,dir + integer(psb_ipk_), intent(in) :: index + integer(psb_ipk_), intent(in) :: dir complex(psb_dpk_), intent(inout) :: heap(:) integer(psb_ipk_), intent(inout) :: idxs(:) integer(psb_ipk_), intent(inout) :: last diff --git a/base/serial/sort/psb_z_msort_impl.f90 b/base/serial/sort/psb_z_msort_impl.f90 index e885176e..525ed572 100644 --- a/base/serial/sort/psb_z_msort_impl.f90 +++ b/base/serial/sort/psb_z_msort_impl.f90 @@ -1,34 +1,34 @@ -! -! 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. -! -! + ! + ! 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. + ! + ! ! ! The merge-sort routines ! References: @@ -41,777 +41,771 @@ ! Addison-Wesley ! - subroutine psb_zmsort_u(x,nout,dir) - use psb_sort_mod, psb_protect_name => psb_zmsort_u - use psb_error_mod - implicit none - complex(psb_dpk_), intent(inout) :: x(:) - integer(psb_ipk_), intent(out) :: nout - integer(psb_ipk_), optional, intent(in) :: dir +subroutine psb_zmsort_u(x,nout,dir) + use psb_sort_mod, psb_protect_name => psb_zmsort_u + use psb_error_mod + implicit none + complex(psb_dpk_), intent(inout) :: x(:) + integer(psb_ipk_), intent(out) :: nout + integer(psb_ipk_), optional, intent(in) :: dir - integer(psb_ipk_) :: n, k - integer(psb_ipk_) :: err_act + integer(psb_ipk_) :: n, k + integer(psb_ipk_) :: err_act - integer(psb_ipk_) :: ierr(5) - character(len=20) :: name + integer(psb_ipk_) :: ierr(5) + character(len=20) :: name - name='psb_msort_u' - call psb_erractionsave(err_act) + name='psb_msort_u' + call psb_erractionsave(err_act) - n = size(x) + n = size(x) - call psb_msort(x,dir=dir) - nout = min(1,n) - do k=2,n - if (x(k) /= x(nout)) then - nout = nout + 1 - x(nout) = x(k) - endif - enddo + call psb_msort(x,dir=dir) + nout = min(1,n) + do k=2,n + if (x(k) /= x(nout)) then + nout = nout + 1 + x(nout) = x(k) + endif + enddo - return + return 9999 call psb_error_handler(err_act) - return - end subroutine psb_zmsort_u - - - - - - - - - subroutine psb_zmsort(x,ix,dir,flag) - use psb_sort_mod, psb_protect_name => psb_zmsort - use psb_error_mod - use psb_ip_reord_mod - implicit none - complex(psb_dpk_), intent(inout) :: x(:) - integer(psb_ipk_), optional, intent(in) :: dir, flag - integer(psb_ipk_), optional, intent(inout) :: ix(:) - - integer(psb_ipk_) :: dir_, flag_, n, err_act - - integer(psb_ipk_), allocatable :: iaux(:) - integer(psb_ipk_) :: iret, info, i - integer(psb_ipk_) :: ierr(5) - character(len=20) :: name - - name='psb_zmsort' - call psb_erractionsave(err_act) - - if (present(dir)) then - dir_ = dir - else - dir_= psb_asort_up_ + return +end subroutine psb_zmsort_u + + +subroutine psb_zmsort(x,ix,dir,flag) + use psb_sort_mod, psb_protect_name => psb_zmsort + use psb_error_mod + use psb_ip_reord_mod + implicit none + complex(psb_dpk_), intent(inout) :: x(:) + integer(psb_ipk_), optional, intent(in) :: dir, flag + integer(psb_ipk_), optional, intent(inout) :: ix(:) + + integer(psb_ipk_) :: dir_, flag_, n, err_act + + integer(psb_ipk_), allocatable :: iaux(:) + integer(psb_ipk_) :: iret, info, i + integer(psb_ipk_) :: ierr(5) + character(len=20) :: name + + name='psb_zmsort' + call psb_erractionsave(err_act) + + if (present(dir)) then + dir_ = dir + else + dir_= psb_asort_up_ + end if + select case(dir_) + case( psb_lsort_up_, psb_lsort_down_, psb_alsort_up_, psb_alsort_down_,& + & psb_asort_up_, psb_asort_down_) + ! OK keep going + case default + ierr(1) = 3; ierr(2) = dir_; + call psb_errpush(psb_err_input_value_invalid_i_,name,i_err=ierr) + goto 9999 + end select + + n = size(x) + + if (present(ix)) then + if (size(ix) < n) then + ierr(1) = 2; ierr(2) = size(ix); + call psb_errpush(psb_err_input_asize_invalid_i_,name,i_err=ierr) + goto 9999 + end if + if (present(flag)) then + flag_ = flag + else + flag_ = psb_sort_ovw_idx_ end if - select case(dir_) - case( psb_lsort_up_, psb_lsort_down_, psb_alsort_up_, psb_alsort_down_,& - & psb_asort_up_, psb_asort_down_) + select case(flag_) + case(psb_sort_ovw_idx_) + do i=1,n + ix(i) = i + end do + case (psb_sort_keep_idx_) ! OK keep going case default - ierr(1) = 3; ierr(2) = dir_; + ierr(1) = 4; ierr(2) = flag_; call psb_errpush(psb_err_input_value_invalid_i_,name,i_err=ierr) goto 9999 end select - - n = size(x) - - if (present(ix)) then - if (size(ix) < n) then - ierr(1) = 2; ierr(2) = size(ix); - call psb_errpush(psb_err_input_asize_invalid_i_,name,i_err=ierr) - goto 9999 - end if - if (present(flag)) then - flag_ = flag - else - flag_ = psb_sort_ovw_idx_ - end if - select case(flag_) - case(psb_sort_ovw_idx_) - do i=1,n - ix(i) = i - end do - case (psb_sort_keep_idx_) - ! OK keep going - case default - ierr(1) = 4; ierr(2) = flag_; - call psb_errpush(psb_err_input_value_invalid_i_,name,i_err=ierr) - goto 9999 - end select + end if + + allocate(iaux(0:n+1),stat=info) + if (info /= psb_success_) then + call psb_errpush(psb_err_alloc_dealloc_,r_name='psb_z_msort') + goto 9999 + endif + + select case(dir_) + case (psb_lsort_up_) + call psi_z_lmsort_up(n,x,iaux,iret) + case (psb_lsort_down_) + call psi_z_lmsort_dw(n,x,iaux,iret) + case (psb_alsort_up_) + call psi_z_almsort_up(n,x,iaux,iret) + case (psb_alsort_down_) + call psi_z_almsort_dw(n,x,iaux,iret) + case (psb_asort_up_) + call psi_z_amsort_up(n,x,iaux,iret) + case (psb_asort_down_) + call psi_z_amsort_dw(n,x,iaux,iret) + end select + ! + ! Do the actual reordering, since the inner routines + ! only provide linked pointers. + ! + if (iret == 0 ) then + if (present(ix)) then + call psb_ip_reord(n,x,ix,iaux) + else + call psb_ip_reord(n,x,iaux) end if + end if - allocate(iaux(0:n+1),stat=info) - if (info /= psb_success_) then - call psb_errpush(psb_err_alloc_dealloc_,r_name='psb_z_msort') - goto 9999 - endif + return - select case(dir_) - case (psb_lsort_up_) - call psi_z_lmsort_up(n,x,iaux,iret) - case (psb_lsort_down_) - call psi_z_lmsort_dw(n,x,iaux,iret) - case (psb_alsort_up_) - call psi_z_almsort_up(n,x,iaux,iret) - case (psb_alsort_down_) - call psi_z_almsort_dw(n,x,iaux,iret) - case (psb_asort_up_) - call psi_z_amsort_up(n,x,iaux,iret) - case (psb_asort_down_) - call psi_z_amsort_dw(n,x,iaux,iret) - end select - ! - ! Do the actual reordering, since the inner routines - ! only provide linked pointers. - ! - if (iret == 0 ) then - if (present(ix)) then - call psb_ip_reord(n,x,ix,iaux) - else - call psb_ip_reord(n,x,iaux) - end if - end if +9999 call psb_error_handler(err_act) - return + return -9999 call psb_error_handler(err_act) - return - - - end subroutine psb_zmsort - - subroutine psi_z_lmsort_up(n,k,l,iret) - use psb_const_mod - use psi_lcx_mod - implicit none - integer(psb_ipk_) :: n, iret - complex(psb_dpk_) :: k(n) - integer(psb_ipk_) :: l(0:n+1) - ! - integer(psb_ipk_) :: p,q,s,t - ! .. - iret = 0 - ! first step: we are preparing ordered sublists, exploiting - ! what order was already in the input data; negative links - ! mark the end of the sublists - l(0) = 1 - t = n + 1 - do p = 1,n - 1 - if (k(p) <= k(p+1)) then - l(p) = p + 1 - else - l(t) = - (p+1) - t = p - end if - end do - l(t) = 0 - l(n) = 0 - ! see if the input was already sorted - if (l(n+1) == 0) then - iret = 1 - return +end subroutine psb_zmsort + +subroutine psi_z_lmsort_up(n,k,l,iret) + use psb_const_mod + use psi_lcx_mod + implicit none + integer(psb_ipk_) :: n, iret + complex(psb_dpk_) :: k(n) + integer(psb_ipk_) :: l(0:n+1) + ! + integer(psb_ipk_) :: p,q,s,t + ! .. + iret = 0 + ! first step: we are preparing ordered sublists, exploiting + ! what order was already in the input data; negative links + ! mark the end of the sublists + l(0) = 1 + t = n + 1 + do p = 1,n - 1 + if (k(p) <= k(p+1)) then + l(p) = p + 1 else - l(n+1) = abs(l(n+1)) + l(t) = - (p+1) + t = p end if + end do + l(t) = 0 + l(n) = 0 + ! see if the input was already sorted + if (l(n+1) == 0) then + iret = 1 + return + else + l(n+1) = abs(l(n+1)) + end if + + mergepass: do + ! otherwise, begin a pass through the list. + ! throughout all the subroutine we have: + ! p, q: pointing to the sublists being merged + ! s: pointing to the most recently processed record + ! t: pointing to the end of previously completed sublist + s = 0 + t = n + 1 + p = l(s) + q = l(t) + if (q == 0) exit mergepass - mergepass: do - ! otherwise, begin a pass through the list. - ! throughout all the subroutine we have: - ! p, q: pointing to the sublists being merged - ! s: pointing to the most recently processed record - ! t: pointing to the end of previously completed sublist - s = 0 - t = n + 1 - p = l(s) - q = l(t) - if (q == 0) exit mergepass - - outer: do + outer: do - if (k(p) > k(q)) then + if (k(p) > k(q)) then - l(s) = sign(q,l(s)) - s = q - q = l(q) - if (q > 0) then - do - if (k(p) <= k(q)) cycle outer - s = q - q = l(q) - if (q <= 0) exit - end do - end if - l(s) = p - s = t + l(s) = sign(q,l(s)) + s = q + q = l(q) + if (q > 0) then do - t = p - p = l(p) - if (p <= 0) exit + if (k(p) <= k(q)) cycle outer + s = q + q = l(q) + if (q <= 0) exit end do + end if + l(s) = p + s = t + do + t = p + p = l(p) + if (p <= 0) exit + end do - else + else - l(s) = sign(p,l(s)) - s = p - p = l(p) - if (p>0) then - do - if (k(p) > k(q)) cycle outer - s = p - p = l(p) - if (p <= 0) exit - end do - end if - ! otherwise, one sublist ended, and we append to it the rest - ! of the other one. - l(s) = q - s = t + l(s) = sign(p,l(s)) + s = p + p = l(p) + if (p>0) then do - t = q - q = l(q) - if (q <= 0) exit + if (k(p) > k(q)) cycle outer + s = p + p = l(p) + if (p <= 0) exit end do end if + ! otherwise, one sublist ended, and we append to it the rest + ! of the other one. + l(s) = q + s = t + do + t = q + q = l(q) + if (q <= 0) exit + end do + end if - p = -p - q = -q - if (q == 0) then - l(s) = sign(p,l(s)) - l(t) = 0 - exit outer - end if - end do outer - end do mergepass - - end subroutine psi_z_lmsort_up - - subroutine psi_z_lmsort_dw(n,k,l,iret) - use psb_const_mod - use psi_lcx_mod - implicit none - integer(psb_ipk_) :: n, iret - complex(psb_dpk_) :: k(n) - integer(psb_ipk_) :: l(0:n+1) - ! - integer(psb_ipk_) :: p,q,s,t - ! .. - iret = 0 - ! first step: we are preparing ordered sublists, exploiting - ! what order was already in the input data; negative links - ! mark the end of the sublists - l(0) = 1 - t = n + 1 - do p = 1,n - 1 - if (k(p) >= k(p+1)) then - l(p) = p + 1 - else - l(t) = - (p+1) - t = p + p = -p + q = -q + if (q == 0) then + l(s) = sign(p,l(s)) + l(t) = 0 + exit outer end if - end do - l(t) = 0 - l(n) = 0 - ! see if the input was already sorted - if (l(n+1) == 0) then - iret = 1 - return + end do outer + end do mergepass + +end subroutine psi_z_lmsort_up + +subroutine psi_z_lmsort_dw(n,k,l,iret) + use psb_const_mod + use psi_lcx_mod + implicit none + integer(psb_ipk_) :: n, iret + complex(psb_dpk_) :: k(n) + integer(psb_ipk_) :: l(0:n+1) + ! + integer(psb_ipk_) :: p,q,s,t + ! .. + iret = 0 + ! first step: we are preparing ordered sublists, exploiting + ! what order was already in the input data; negative links + ! mark the end of the sublists + l(0) = 1 + t = n + 1 + do p = 1,n - 1 + if (k(p) >= k(p+1)) then + l(p) = p + 1 else - l(n+1) = abs(l(n+1)) + l(t) = - (p+1) + t = p end if + end do + l(t) = 0 + l(n) = 0 + ! see if the input was already sorted + if (l(n+1) == 0) then + iret = 1 + return + else + l(n+1) = abs(l(n+1)) + end if + + mergepass: do + ! otherwise, begin a pass through the list. + ! throughout all the subroutine we have: + ! p, q: pointing to the sublists being merged + ! s: pointing to the most recently processed record + ! t: pointing to the end of previously completed sublist + s = 0 + t = n + 1 + p = l(s) + q = l(t) + if (q == 0) exit mergepass - mergepass: do - ! otherwise, begin a pass through the list. - ! throughout all the subroutine we have: - ! p, q: pointing to the sublists being merged - ! s: pointing to the most recently processed record - ! t: pointing to the end of previously completed sublist - s = 0 - t = n + 1 - p = l(s) - q = l(t) - if (q == 0) exit mergepass - - outer: do + outer: do - if (k(p) < k(q)) then + if (k(p) < k(q)) then - l(s) = sign(q,l(s)) - s = q - q = l(q) - if (q > 0) then - do - if (k(p) >= k(q)) cycle outer - s = q - q = l(q) - if (q <= 0) exit - end do - end if - l(s) = p - s = t + l(s) = sign(q,l(s)) + s = q + q = l(q) + if (q > 0) then do - t = p - p = l(p) - if (p <= 0) exit + if (k(p) >= k(q)) cycle outer + s = q + q = l(q) + if (q <= 0) exit end do + end if + l(s) = p + s = t + do + t = p + p = l(p) + if (p <= 0) exit + end do - else + else - l(s) = sign(p,l(s)) - s = p - p = l(p) - if (p>0) then - do - if (k(p) < k(q)) cycle outer - s = p - p = l(p) - if (p <= 0) exit - end do - end if - ! otherwise, one sublist ended, and we append to it the rest - ! of the other one. - l(s) = q - s = t + l(s) = sign(p,l(s)) + s = p + p = l(p) + if (p>0) then do - t = q - q = l(q) - if (q <= 0) exit + if (k(p) < k(q)) cycle outer + s = p + p = l(p) + if (p <= 0) exit end do end if + ! otherwise, one sublist ended, and we append to it the rest + ! of the other one. + l(s) = q + s = t + do + t = q + q = l(q) + if (q <= 0) exit + end do + end if - p = -p - q = -q - if (q == 0) then - l(s) = sign(p,l(s)) - l(t) = 0 - exit outer - end if - end do outer - end do mergepass - - end subroutine psi_z_lmsort_dw - - subroutine psi_z_amsort_up(n,k,l,iret) - use psb_const_mod - use psi_acx_mod - implicit none - integer(psb_ipk_) :: n, iret - complex(psb_dpk_) :: k(n) - integer(psb_ipk_) :: l(0:n+1) - ! - integer(psb_ipk_) :: p,q,s,t - ! .. - iret = 0 - ! first step: we are preparing ordered sublists, exploiting - ! what order was already in the input data; negative links - ! mark the end of the sublists - l(0) = 1 - t = n + 1 - do p = 1,n - 1 - if (k(p) <= k(p+1)) then - l(p) = p + 1 - else - l(t) = - (p+1) - t = p + p = -p + q = -q + if (q == 0) then + l(s) = sign(p,l(s)) + l(t) = 0 + exit outer end if - end do - l(t) = 0 - l(n) = 0 - ! see if the input was already sorted - if (l(n+1) == 0) then - iret = 1 - return + end do outer + end do mergepass + +end subroutine psi_z_lmsort_dw + +subroutine psi_z_amsort_up(n,k,l,iret) + use psb_const_mod + use psi_acx_mod + implicit none + integer(psb_ipk_) :: n, iret + complex(psb_dpk_) :: k(n) + integer(psb_ipk_) :: l(0:n+1) + ! + integer(psb_ipk_) :: p,q,s,t + ! .. + iret = 0 + ! first step: we are preparing ordered sublists, exploiting + ! what order was already in the input data; negative links + ! mark the end of the sublists + l(0) = 1 + t = n + 1 + do p = 1,n - 1 + if (k(p) <= k(p+1)) then + l(p) = p + 1 else - l(n+1) = abs(l(n+1)) + l(t) = - (p+1) + t = p end if + end do + l(t) = 0 + l(n) = 0 + ! see if the input was already sorted + if (l(n+1) == 0) then + iret = 1 + return + else + l(n+1) = abs(l(n+1)) + end if + + mergepass: do + ! otherwise, begin a pass through the list. + ! throughout all the subroutine we have: + ! p, q: pointing to the sublists being merged + ! s: pointing to the most recently processed record + ! t: pointing to the end of previously completed sublist + s = 0 + t = n + 1 + p = l(s) + q = l(t) + if (q == 0) exit mergepass - mergepass: do - ! otherwise, begin a pass through the list. - ! throughout all the subroutine we have: - ! p, q: pointing to the sublists being merged - ! s: pointing to the most recently processed record - ! t: pointing to the end of previously completed sublist - s = 0 - t = n + 1 - p = l(s) - q = l(t) - if (q == 0) exit mergepass - - outer: do + outer: do - if (k(p) > k(q)) then + if (k(p) > k(q)) then - l(s) = sign(q,l(s)) - s = q - q = l(q) - if (q > 0) then - do - if (k(p) <= k(q)) cycle outer - s = q - q = l(q) - if (q <= 0) exit - end do - end if - l(s) = p - s = t + l(s) = sign(q,l(s)) + s = q + q = l(q) + if (q > 0) then do - t = p - p = l(p) - if (p <= 0) exit + if (k(p) <= k(q)) cycle outer + s = q + q = l(q) + if (q <= 0) exit end do + end if + l(s) = p + s = t + do + t = p + p = l(p) + if (p <= 0) exit + end do - else + else - l(s) = sign(p,l(s)) - s = p - p = l(p) - if (p>0) then - do - if (k(p) > k(q)) cycle outer - s = p - p = l(p) - if (p <= 0) exit - end do - end if - ! otherwise, one sublist ended, and we append to it the rest - ! of the other one. - l(s) = q - s = t + l(s) = sign(p,l(s)) + s = p + p = l(p) + if (p>0) then do - t = q - q = l(q) - if (q <= 0) exit + if (k(p) > k(q)) cycle outer + s = p + p = l(p) + if (p <= 0) exit end do end if + ! otherwise, one sublist ended, and we append to it the rest + ! of the other one. + l(s) = q + s = t + do + t = q + q = l(q) + if (q <= 0) exit + end do + end if - p = -p - q = -q - if (q == 0) then - l(s) = sign(p,l(s)) - l(t) = 0 - exit outer - end if - end do outer - end do mergepass - - end subroutine psi_z_amsort_up - - subroutine psi_z_amsort_dw(n,k,l,iret) - use psb_const_mod - use psi_acx_mod - implicit none - integer(psb_ipk_) :: n, iret - complex(psb_dpk_) :: k(n) - integer(psb_ipk_) :: l(0:n+1) - ! - integer(psb_ipk_) :: p,q,s,t - ! .. - iret = 0 - ! first step: we are preparing ordered sublists, exploiting - ! what order was already in the input data; negative links - ! mark the end of the sublists - l(0) = 1 - t = n + 1 - do p = 1,n - 1 - if (k(p) >= k(p+1)) then - l(p) = p + 1 - else - l(t) = - (p+1) - t = p + p = -p + q = -q + if (q == 0) then + l(s) = sign(p,l(s)) + l(t) = 0 + exit outer end if - end do - l(t) = 0 - l(n) = 0 - ! see if the input was already sorted - if (l(n+1) == 0) then - iret = 1 - return + end do outer + end do mergepass + +end subroutine psi_z_amsort_up + +subroutine psi_z_amsort_dw(n,k,l,iret) + use psb_const_mod + use psi_acx_mod + implicit none + integer(psb_ipk_) :: n, iret + complex(psb_dpk_) :: k(n) + integer(psb_ipk_) :: l(0:n+1) + ! + integer(psb_ipk_) :: p,q,s,t + ! .. + iret = 0 + ! first step: we are preparing ordered sublists, exploiting + ! what order was already in the input data; negative links + ! mark the end of the sublists + l(0) = 1 + t = n + 1 + do p = 1,n - 1 + if (k(p) >= k(p+1)) then + l(p) = p + 1 else - l(n+1) = abs(l(n+1)) + l(t) = - (p+1) + t = p end if + end do + l(t) = 0 + l(n) = 0 + ! see if the input was already sorted + if (l(n+1) == 0) then + iret = 1 + return + else + l(n+1) = abs(l(n+1)) + end if + + mergepass: do + ! otherwise, begin a pass through the list. + ! throughout all the subroutine we have: + ! p, q: pointing to the sublists being merged + ! s: pointing to the most recently processed record + ! t: pointing to the end of previously completed sublist + s = 0 + t = n + 1 + p = l(s) + q = l(t) + if (q == 0) exit mergepass - mergepass: do - ! otherwise, begin a pass through the list. - ! throughout all the subroutine we have: - ! p, q: pointing to the sublists being merged - ! s: pointing to the most recently processed record - ! t: pointing to the end of previously completed sublist - s = 0 - t = n + 1 - p = l(s) - q = l(t) - if (q == 0) exit mergepass - - outer: do + outer: do - if (k(p) < k(q)) then + if (k(p) < k(q)) then - l(s) = sign(q,l(s)) - s = q - q = l(q) - if (q > 0) then - do - if (k(p) >= k(q)) cycle outer - s = q - q = l(q) - if (q <= 0) exit - end do - end if - l(s) = p - s = t + l(s) = sign(q,l(s)) + s = q + q = l(q) + if (q > 0) then do - t = p - p = l(p) - if (p <= 0) exit + if (k(p) >= k(q)) cycle outer + s = q + q = l(q) + if (q <= 0) exit end do + end if + l(s) = p + s = t + do + t = p + p = l(p) + if (p <= 0) exit + end do - else + else - l(s) = sign(p,l(s)) - s = p - p = l(p) - if (p>0) then - do - if (k(p) < k(q)) cycle outer - s = p - p = l(p) - if (p <= 0) exit - end do - end if - ! otherwise, one sublist ended, and we append to it the rest - ! of the other one. - l(s) = q - s = t + l(s) = sign(p,l(s)) + s = p + p = l(p) + if (p>0) then do - t = q - q = l(q) - if (q <= 0) exit + if (k(p) < k(q)) cycle outer + s = p + p = l(p) + if (p <= 0) exit end do end if + ! otherwise, one sublist ended, and we append to it the rest + ! of the other one. + l(s) = q + s = t + do + t = q + q = l(q) + if (q <= 0) exit + end do + end if - p = -p - q = -q - if (q == 0) then - l(s) = sign(p,l(s)) - l(t) = 0 - exit outer - end if - end do outer - end do mergepass - - end subroutine psi_z_amsort_dw - - subroutine psi_z_almsort_up(n,k,l,iret) - use psb_const_mod - use psi_alcx_mod - implicit none - integer(psb_ipk_) :: n, iret - complex(psb_dpk_) :: k(n) - integer(psb_ipk_) :: l(0:n+1) - ! - integer(psb_ipk_) :: p,q,s,t - ! .. - iret = 0 - ! first step: we are preparing ordered sublists, exploiting - ! what order was already in the input data; negative links - ! mark the end of the sublists - l(0) = 1 - t = n + 1 - do p = 1,n - 1 - if (k(p) <= k(p+1)) then - l(p) = p + 1 - else - l(t) = - (p+1) - t = p + p = -p + q = -q + if (q == 0) then + l(s) = sign(p,l(s)) + l(t) = 0 + exit outer end if - end do - l(t) = 0 - l(n) = 0 - ! see if the input was already sorted - if (l(n+1) == 0) then - iret = 1 - return + end do outer + end do mergepass + +end subroutine psi_z_amsort_dw + +subroutine psi_z_almsort_up(n,k,l,iret) + use psb_const_mod + use psi_alcx_mod + implicit none + integer(psb_ipk_) :: n, iret + complex(psb_dpk_) :: k(n) + integer(psb_ipk_) :: l(0:n+1) + ! + integer(psb_ipk_) :: p,q,s,t + ! .. + iret = 0 + ! first step: we are preparing ordered sublists, exploiting + ! what order was already in the input data; negative links + ! mark the end of the sublists + l(0) = 1 + t = n + 1 + do p = 1,n - 1 + if (k(p) <= k(p+1)) then + l(p) = p + 1 else - l(n+1) = abs(l(n+1)) + l(t) = - (p+1) + t = p end if + end do + l(t) = 0 + l(n) = 0 + ! see if the input was already sorted + if (l(n+1) == 0) then + iret = 1 + return + else + l(n+1) = abs(l(n+1)) + end if + + mergepass: do + ! otherwise, begin a pass through the list. + ! throughout all the subroutine we have: + ! p, q: pointing to the sublists being merged + ! s: pointing to the most recently processed record + ! t: pointing to the end of previously completed sublist + s = 0 + t = n + 1 + p = l(s) + q = l(t) + if (q == 0) exit mergepass - mergepass: do - ! otherwise, begin a pass through the list. - ! throughout all the subroutine we have: - ! p, q: pointing to the sublists being merged - ! s: pointing to the most recently processed record - ! t: pointing to the end of previously completed sublist - s = 0 - t = n + 1 - p = l(s) - q = l(t) - if (q == 0) exit mergepass - - outer: do + outer: do - if (k(p) > k(q)) then + if (k(p) > k(q)) then - l(s) = sign(q,l(s)) - s = q - q = l(q) - if (q > 0) then - do - if (k(p) <= k(q)) cycle outer - s = q - q = l(q) - if (q <= 0) exit - end do - end if - l(s) = p - s = t + l(s) = sign(q,l(s)) + s = q + q = l(q) + if (q > 0) then do - t = p - p = l(p) - if (p <= 0) exit + if (k(p) <= k(q)) cycle outer + s = q + q = l(q) + if (q <= 0) exit end do + end if + l(s) = p + s = t + do + t = p + p = l(p) + if (p <= 0) exit + end do - else + else - l(s) = sign(p,l(s)) - s = p - p = l(p) - if (p>0) then - do - if (k(p) > k(q)) cycle outer - s = p - p = l(p) - if (p <= 0) exit - end do - end if - ! otherwise, one sublist ended, and we append to it the rest - ! of the other one. - l(s) = q - s = t + l(s) = sign(p,l(s)) + s = p + p = l(p) + if (p>0) then do - t = q - q = l(q) - if (q <= 0) exit + if (k(p) > k(q)) cycle outer + s = p + p = l(p) + if (p <= 0) exit end do end if + ! otherwise, one sublist ended, and we append to it the rest + ! of the other one. + l(s) = q + s = t + do + t = q + q = l(q) + if (q <= 0) exit + end do + end if - p = -p - q = -q - if (q == 0) then - l(s) = sign(p,l(s)) - l(t) = 0 - exit outer - end if - end do outer - end do mergepass - - end subroutine psi_z_almsort_up - - subroutine psi_z_almsort_dw(n,k,l,iret) - use psb_const_mod - use psi_alcx_mod - implicit none - integer(psb_ipk_) :: n, iret - complex(psb_dpk_) :: k(n) - integer(psb_ipk_) :: l(0:n+1) - ! - integer(psb_ipk_) :: p,q,s,t - ! .. - iret = 0 - ! first step: we are preparing ordered sublists, exploiting - ! what order was already in the input data; negative links - ! mark the end of the sublists - l(0) = 1 - t = n + 1 - do p = 1,n - 1 - if (k(p) >= k(p+1)) then - l(p) = p + 1 - else - l(t) = - (p+1) - t = p + p = -p + q = -q + if (q == 0) then + l(s) = sign(p,l(s)) + l(t) = 0 + exit outer end if - end do - l(t) = 0 - l(n) = 0 - ! see if the input was already sorted - if (l(n+1) == 0) then - iret = 1 - return + end do outer + end do mergepass + +end subroutine psi_z_almsort_up + +subroutine psi_z_almsort_dw(n,k,l,iret) + use psb_const_mod + use psi_alcx_mod + implicit none + integer(psb_ipk_) :: n, iret + complex(psb_dpk_) :: k(n) + integer(psb_ipk_) :: l(0:n+1) + ! + integer(psb_ipk_) :: p,q,s,t + ! .. + iret = 0 + ! first step: we are preparing ordered sublists, exploiting + ! what order was already in the input data; negative links + ! mark the end of the sublists + l(0) = 1 + t = n + 1 + do p = 1,n - 1 + if (k(p) >= k(p+1)) then + l(p) = p + 1 else - l(n+1) = abs(l(n+1)) + l(t) = - (p+1) + t = p end if + end do + l(t) = 0 + l(n) = 0 + ! see if the input was already sorted + if (l(n+1) == 0) then + iret = 1 + return + else + l(n+1) = abs(l(n+1)) + end if + + mergepass: do + ! otherwise, begin a pass through the list. + ! throughout all the subroutine we have: + ! p, q: pointing to the sublists being merged + ! s: pointing to the most recently processed record + ! t: pointing to the end of previously completed sublist + s = 0 + t = n + 1 + p = l(s) + q = l(t) + if (q == 0) exit mergepass - mergepass: do - ! otherwise, begin a pass through the list. - ! throughout all the subroutine we have: - ! p, q: pointing to the sublists being merged - ! s: pointing to the most recently processed record - ! t: pointing to the end of previously completed sublist - s = 0 - t = n + 1 - p = l(s) - q = l(t) - if (q == 0) exit mergepass - - outer: do + outer: do - if (k(p) < k(q)) then + if (k(p) < k(q)) then - l(s) = sign(q,l(s)) - s = q - q = l(q) - if (q > 0) then - do - if (k(p) >= k(q)) cycle outer - s = q - q = l(q) - if (q <= 0) exit - end do - end if - l(s) = p - s = t + l(s) = sign(q,l(s)) + s = q + q = l(q) + if (q > 0) then do - t = p - p = l(p) - if (p <= 0) exit + if (k(p) >= k(q)) cycle outer + s = q + q = l(q) + if (q <= 0) exit end do + end if + l(s) = p + s = t + do + t = p + p = l(p) + if (p <= 0) exit + end do - else + else - l(s) = sign(p,l(s)) - s = p - p = l(p) - if (p>0) then - do - if (k(p) < k(q)) cycle outer - s = p - p = l(p) - if (p <= 0) exit - end do - end if - ! otherwise, one sublist ended, and we append to it the rest - ! of the other one. - l(s) = q - s = t + l(s) = sign(p,l(s)) + s = p + p = l(p) + if (p>0) then do - t = q - q = l(q) - if (q <= 0) exit + if (k(p) < k(q)) cycle outer + s = p + p = l(p) + if (p <= 0) exit end do end if + ! otherwise, one sublist ended, and we append to it the rest + ! of the other one. + l(s) = q + s = t + do + t = q + q = l(q) + if (q <= 0) exit + end do + end if - p = -p - q = -q - if (q == 0) then - l(s) = sign(p,l(s)) - l(t) = 0 - exit outer - end if - end do outer - end do mergepass + p = -p + q = -q + if (q == 0) then + l(s) = sign(p,l(s)) + l(t) = 0 + exit outer + end if + end do outer + end do mergepass - end subroutine psi_z_almsort_dw +end subroutine psi_z_almsort_dw