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

Sort implementation

Browse Source
ILmat
Salvatore Filippone 8 years ago
parent f862fe048a
commit cecd34823e
9 changed files with 6101 additions and 6 deletions
Show Stats Download Patch File Download Diff File

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

@ -45,12 +45,7 @@
module psb_sort_mod module psb_sort_mod
use psb_const_mod use psb_const_mod
use psb_ip_reord_mod use psb_ip_reord_mod
!!$ use psb_i_sort_mod
!!$ use psb_l_sort_mod
!!$ use psb_s_sort_mod
!!$ use psb_c_sort_mod
!!$ use psb_d_sort_mod
!!$ use psb_z_sort_mod
use psb_m_hsort_mod use psb_m_hsort_mod
use psb_m_isort_mod use psb_m_isort_mod
use psb_m_msort_mod use psb_m_msort_mod

678
base/serial/sort/psb_e_hsort_impl.f90
Unescape Escape View File

@ -0,0 +1,678 @@
!
! 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 and quicksort routines are implemented in the
! serial/aux directory
! References:
! D. Knuth
! The Art of Computer Programming, vol. 3
! Addison-Wesley
!
! Aho, Hopcroft, Ullman
! Data Structures and Algorithms
! Addison-Wesley
!
subroutine psb_ehsort(x,ix,dir,flag)
use psb_sort_mod, psb_protect_name => psb_ehsort
use psb_error_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, i, l, err_act,info
integer(psb_epk_) :: key
integer(psb_epk_) :: index
integer(psb_ipk_) :: ierr(5)
character(len=20) :: name
name='psb_hsort'
call psb_erractionsave(err_act)
if (present(flag)) then
flag_ = flag
else
flag_ = psb_sort_ovw_idx_
end if
select case(flag_)
case( psb_sort_ovw_idx_, 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
if (present(dir)) then
dir_ = dir
else
dir_= psb_sort_up_
end if
select case(dir_)
case(psb_sort_up_,psb_sort_down_)
! OK
case (psb_asort_up_,psb_asort_down_)
! OK
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)
!
! Dirty trick to sort with heaps: if we want
! to sort in place upwards, first we set up a heap so that
! we can easily get the LARGEST element, then we take it out
! and put it in the last entry, and so on.
! So, we invert dir_
!
dir_ = -dir_
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 (flag_ == psb_sort_ovw_idx_) then
do i=1, n
ix(i) = i
end do
end if
l = 0
do i=1, n
key = x(i)
index = ix(i)
call psi_idx_insert_heap(key,index,l,x,ix,dir_,info)
if (l /= i) then
write(psb_err_unit,*) 'Mismatch while heapifying ! '
end if
end do
do i=n, 2, -1
call psi_idx_heap_get_first(key,index,l,x,ix,dir_,info)
if (l /= i-1) then
write(psb_err_unit,*) 'Mismatch while pulling out of heap ',l,i
end if
x(i) = key
ix(i) = index
end do
else if (.not.present(ix)) then
l = 0
do i=1, n
key = x(i)
call psi_insert_heap(key,l,x,dir_,info)
if (l /= i) then
write(psb_err_unit,*) 'Mismatch while heapifying ! ',l,i
end if
end do
do i=n, 2, -1
call psi_e_heap_get_first(key,l,x,dir_,info)
if (l /= i-1) then
write(psb_err_unit,*) 'Mismatch while pulling out of heap ',l,i
end if
x(i) = key
end do
end if
return
9999 call psb_error_handler(err_act)
return
end subroutine psb_ehsort
!
! These are packaged so that they can be used to implement
! a heapsort, should the need arise
!
!
! Programming note:
! In the implementation of the heap_get_first function
! we have code like this
!
! if ( ( heap(2*i) < heap(2*i+1) ) .or.&
! & (2*i == last)) then
! j = 2*i
! else
! j = 2*i + 1
! end if
!
! It looks like the 2*i+1 could overflow the array, but this
! is not true because there is a guard statement
! if (i>last/2) exit
! and because last has just been reduced by 1 when defining the return value,
! therefore 2*i+1 may be greater than the current value of last,
! but cannot be greater than the value of last when the routine was entered
! hence it is safe.
!
!
!
subroutine psi_e_insert_heap(key,last,heap,dir,info)
use psb_sort_mod, psb_protect_name => psi_e_insert_heap
implicit none
!
! Input:
! key: the new value
! last: pointer to the last occupied element in heap
! heap: the heap
! dir: sorting direction
integer(psb_epk_), intent(in) :: key
integer(psb_ipk_), intent(in) :: dir
integer(psb_epk_), intent(inout) :: heap(:)
integer(psb_epk_), intent(inout) :: last
integer(psb_ipk_), intent(out) :: info
integer(psb_ipk_) :: i, i2
integer(psb_epk_) :: temp
info = psb_success_
if (last < 0) then
write(psb_err_unit,*) 'Invalid last in heap ',last
info = last
return
endif
last = last + 1
if (last > size(heap)) then
write(psb_err_unit,*) 'out of bounds '
info = -1
return
end if
i = last
heap(i) = key
select case(dir)
case (psb_sort_up_)
do
if (i<=1) exit
i2 = i/2
if (heap(i) < heap(i2)) then
temp = heap(i)
heap(i) = heap(i2)
heap(i2) = temp
i = i2
else
exit
end if
end do
case (psb_sort_down_)
do
if (i<=1) exit
i2 = i/2
if (heap(i) > heap(i2)) then
temp = heap(i)
heap(i) = heap(i2)
heap(i2) = temp
i = i2
else
exit
end if
end do
case (psb_asort_up_)
do
if (i<=1) exit
i2 = i/2
if (abs(heap(i)) < abs(heap(i2))) then
temp = heap(i)
heap(i) = heap(i2)
heap(i2) = temp
i = i2
else
exit
end if
end do
case (psb_asort_down_)
do
if (i<=1) exit
i2 = i/2
if (abs(heap(i)) > abs(heap(i2))) then
temp = heap(i)
heap(i) = heap(i2)
heap(i2) = temp
i = i2
else
exit
end if
end do
case default
write(psb_err_unit,*) 'Invalid direction in heap ',dir
end select
return
end subroutine psi_e_insert_heap
subroutine psi_e_heap_get_first(key,last,heap,dir,info)
use psb_sort_mod, psb_protect_name => psi_e_heap_get_first
implicit none
integer(psb_epk_), intent(inout) :: key
integer(psb_epk_), intent(inout) :: last
integer(psb_ipk_), intent(in) :: dir
integer(psb_epk_), intent(inout) :: heap(:)
integer(psb_ipk_), intent(out) :: info
integer(psb_ipk_) :: i, j
integer(psb_epk_) :: temp
info = psb_success_
if (last <= 0) then
key = 0
info = -1
return
endif
key = heap(1)
heap(1) = heap(last)
last = last - 1
select case(dir)
case (psb_sort_up_)
i = 1
do
if (i > (last/2)) exit
if ( (heap(2*i) < heap(2*i+1)) .or.&
& (2*i == last)) then
j = 2*i
else
j = 2*i + 1
end if
if (heap(i) > heap(j)) then
temp = heap(i)
heap(i) = heap(j)
heap(j) = temp
i = j
else
exit
end if
end do
case (psb_sort_down_)
i = 1
do
if (i > (last/2)) exit
if ( (heap(2*i) > heap(2*i+1)) .or.&
& (2*i == last)) then
j = 2*i
else
j = 2*i + 1
end if
if (heap(i) < heap(j)) then
temp = heap(i)
heap(i) = heap(j)
heap(j) = temp
i = j
else
exit
end if
end do
case (psb_asort_up_)
i = 1
do
if (i > (last/2)) exit
if ( (abs(heap(2*i)) < abs(heap(2*i+1))) .or.&
& (2*i == last)) then
j = 2*i
else
j = 2*i + 1
end if
if (abs(heap(i)) > abs(heap(j))) then
temp = heap(i)
heap(i) = heap(j)
heap(j) = temp
i = j
else
exit
end if
end do
case (psb_asort_down_)
i = 1
do
if (i > (last/2)) exit
if ( (abs(heap(2*i)) > abs(heap(2*i+1))) .or.&
& (2*i == last)) then
j = 2*i
else
j = 2*i + 1
end if
if (abs(heap(i)) < abs(heap(j))) then
temp = heap(i)
heap(i) = heap(j)
heap(j) = temp
i = j
else
exit
end if
end do
case default
write(psb_err_unit,*) 'Invalid direction in heap ',dir
end select
return
end subroutine psi_e_heap_get_first
subroutine psi_e_idx_insert_heap(key,index,last,heap,idxs,dir,info)
use psb_sort_mod, psb_protect_name => psi_e_idx_insert_heap
implicit none
!
! Input:
! key: the new value
! index: the new index
! last: pointer to the last occupied element in heap
! heap: the heap
! idxs: the indices
! dir: sorting direction
integer(psb_epk_), intent(in) :: key
integer(psb_ipk_), intent(in) :: index,dir
integer(psb_epk_), intent(inout) :: heap(:)
integer(psb_ipk_), intent(inout) :: idxs(:),last
integer(psb_ipk_), intent(out) :: info
integer(psb_ipk_) :: i, i2, itemp
integer(psb_epk_) :: temp
info = psb_success_
if (last < 0) then
write(psb_err_unit,*) 'Invalid last in heap ',last
info = last
return
endif
last = last + 1
if (last > size(heap)) then
write(psb_err_unit,*) 'out of bounds '
info = -1
return
end if
i = last
heap(i) = key
idxs(i) = index
select case(dir)
case (psb_sort_up_)
do
if (i<=1) exit
i2 = i/2
if (heap(i) < heap(i2)) then
itemp = idxs(i)
idxs(i) = idxs(i2)
idxs(i2) = itemp
temp = heap(i)
heap(i) = heap(i2)
heap(i2) = temp
i = i2
else
exit
end if
end do
case (psb_sort_down_)
do
if (i<=1) exit
i2 = i/2
if (heap(i) > heap(i2)) then
itemp = idxs(i)
idxs(i) = idxs(i2)
idxs(i2) = itemp
temp = heap(i)
heap(i) = heap(i2)
heap(i2) = temp
i = i2
else
exit
end if
end do
case (psb_asort_up_)
do
if (i<=1) exit
i2 = i/2
if (abs(heap(i)) < abs(heap(i2))) then
itemp = idxs(i)
idxs(i) = idxs(i2)
idxs(i2) = itemp
temp = heap(i)
heap(i) = heap(i2)
heap(i2) = temp
i = i2
else
exit
end if
end do
case (psb_asort_down_)
do
if (i<=1) exit
i2 = i/2
if (abs(heap(i)) > abs(heap(i2))) then
itemp = idxs(i)
idxs(i) = idxs(i2)
idxs(i2) = itemp
temp = heap(i)
heap(i) = heap(i2)
heap(i2) = temp
i = i2
else
exit
end if
end do
case default
write(psb_err_unit,*) 'Invalid direction in heap ',dir
end select
return
end subroutine psi_e_idx_insert_heap
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) :: 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_ipk_) :: i, j,itemp
integer(psb_epk_) :: temp
info = psb_success_
if (last <= 0) then
key = 0
index = 0
info = -1
return
endif
key = heap(1)
index = idxs(1)
heap(1) = heap(last)
idxs(1) = idxs(last)
last = last - 1
select case(dir)
case (psb_sort_up_)
i = 1
do
if (i > (last/2)) exit
if ( (heap(2*i) < heap(2*i+1)) .or.&
& (2*i == last)) then
j = 2*i
else
j = 2*i + 1
end if
if (heap(i) > heap(j)) then
itemp = idxs(i)
idxs(i) = idxs(j)
idxs(j) = itemp
temp = heap(i)
heap(i) = heap(j)
heap(j) = temp
i = j
else
exit
end if
end do
case (psb_sort_down_)
i = 1
do
if (i > (last/2)) exit
if ( (heap(2*i) > heap(2*i+1)) .or.&
& (2*i == last)) then
j = 2*i
else
j = 2*i + 1
end if
if (heap(i) < heap(j)) then
itemp = idxs(i)
idxs(i) = idxs(j)
idxs(j) = itemp
temp = heap(i)
heap(i) = heap(j)
heap(j) = temp
i = j
else
exit
end if
end do
case (psb_asort_up_)
i = 1
do
if (i > (last/2)) exit
if ( (abs(heap(2*i)) < abs(heap(2*i+1))) .or.&
& (2*i == last)) then
j = 2*i
else
j = 2*i + 1
end if
if (abs(heap(i)) > abs(heap(j))) then
itemp = idxs(i)
idxs(i) = idxs(j)
idxs(j) = itemp
temp = heap(i)
heap(i) = heap(j)
heap(j) = temp
i = j
else
exit
end if
end do
case (psb_asort_down_)
i = 1
do
if (i > (last/2)) exit
if ( (abs(heap(2*i)) > abs(heap(2*i+1))) .or.&
& (2*i == last)) then
j = 2*i
else
j = 2*i + 1
end if
if (abs(heap(i)) < abs(heap(j))) then
itemp = idxs(i)
idxs(i) = idxs(j)
idxs(j) = itemp
temp = heap(i)
heap(i) = heap(j)
heap(j) = temp
i = j
else
exit
end if
end do
case default
write(psb_err_unit,*) 'Invalid direction in heap ',dir
end select
return
end subroutine psi_e_idx_heap_get_first

341
base/serial/sort/psb_e_isort_impl.f90
Unescape Escape View File

@ -0,0 +1,341 @@
!
! 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 insertion sort routines
! References:
! D. Knuth
! The Art of Computer Programming, vol. 3
! Addison-Wesley
!
! Aho, Hopcroft, Ullman
! Data Structures and Algorithms
! Addison-Wesley
!
subroutine psb_eisort(x,ix,dir,flag)
use psb_sort_mod, psb_protect_name => psb_eisort
use psb_error_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_, err_act
integer(psb_epk_) :: n, i
integer(psb_ipk_) :: ierr(5)
character(len=20) :: name
name='psb_eisort'
call psb_erractionsave(err_act)
if (present(flag)) then
flag_ = flag
else
flag_ = psb_sort_ovw_idx_
end if
select case(flag_)
case( psb_sort_ovw_idx_, 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
if (present(dir)) then
dir_ = dir
else
dir_= psb_sort_up_
end if
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 (flag_==psb_sort_ovw_idx_) then
do i=1,n
ix(i) = i
end do
end if
select case(dir_)
case (psb_sort_up_)
call psi_eisrx_up(n,x,ix)
case (psb_sort_down_)
call psi_eisrx_dw(n,x,ix)
case (psb_asort_up_)
call psi_eaisrx_up(n,x,ix)
case (psb_asort_down_)
call psi_eaisrx_dw(n,x,ix)
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
else
select case(dir_)
case (psb_sort_up_)
call psi_eisr_up(n,x)
case (psb_sort_down_)
call psi_eisr_dw(n,x)
case (psb_asort_up_)
call psi_eaisr_up(n,x)
case (psb_asort_down_)
call psi_eaisr_dw(n,x)
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
end if
return
9999 call psb_error_handler(err_act)
return
end subroutine psb_eisort
subroutine psi_eisrx_up(n,x,idx)
use psb_sort_mod, psb_protect_name => psi_eisrx_up
use psb_error_mod
implicit none
integer(psb_epk_), intent(inout) :: x(:)
integer(psb_epk_), intent(inout) :: idx(:)
integer(psb_epk_), intent(in) :: n
integer(psb_epk_) :: i,j,ix
integer(psb_epk_) :: xx
do j=n-1,1,-1
if (x(j+1) < x(j)) then
xx = x(j)
ix = idx(j)
i=j+1
do
x(i-1) = x(i)
idx(i-1) = idx(i)
i = i+1
if (i>n) exit
if (x(i) >= xx) exit
end do
x(i-1) = xx
idx(i-1) = ix
endif
enddo
end subroutine psi_eisrx_up
subroutine psi_eisrx_dw(n,x,idx)
use psb_sort_mod, psb_protect_name => psi_eisrx_dw
use psb_error_mod
implicit none
integer(psb_epk_), intent(inout) :: x(:)
integer(psb_epk_), intent(inout) :: idx(:)
integer(psb_epk_), intent(in) :: n
integer(psb_epk_) :: i,j,ix
integer(psb_epk_) :: xx
do j=n-1,1,-1
if (x(j+1) > x(j)) then
xx = x(j)
ix = idx(j)
i=j+1
do
x(i-1) = x(i)
idx(i-1) = idx(i)
i = i+1
if (i>n) exit
if (x(i) <= xx) exit
end do
x(i-1) = xx
idx(i-1) = ix
endif
enddo
end subroutine psi_eisrx_dw
subroutine psi_eisr_up(n,x)
use psb_sort_mod, psb_protect_name => psi_eisr_up
use psb_error_mod
implicit none
integer(psb_epk_), intent(inout) :: x(:)
integer(psb_epk_), intent(in) :: n
integer(psb_epk_) :: i,j
integer(psb_epk_) :: xx
do j=n-1,1,-1
if (x(j+1) < x(j)) then
xx = x(j)
i=j+1
do
x(i-1) = x(i)
i = i+1
if (i>n) exit
if (x(i) >= xx) exit
end do
x(i-1) = xx
endif
enddo
end subroutine psi_eisr_up
subroutine psi_eisr_dw(n,x)
use psb_sort_mod, psb_protect_name => psi_eisr_dw
use psb_error_mod
implicit none
integer(psb_epk_), intent(inout) :: x(:)
integer(psb_epk_), intent(in) :: n
integer(psb_epk_) :: i,j
integer(psb_epk_) :: xx
do j=n-1,1,-1
if (x(j+1) > x(j)) then
xx = x(j)
i=j+1
do
x(i-1) = x(i)
i = i+1
if (i>n) exit
if (x(i) <= xx) exit
end do
x(i-1) = xx
endif
enddo
end subroutine psi_eisr_dw
subroutine psi_eaisrx_up(n,x,idx)
use psb_sort_mod, psb_protect_name => psi_eaisrx_up
use psb_error_mod
implicit none
integer(psb_epk_), intent(inout) :: x(:)
integer(psb_epk_), intent(inout) :: idx(:)
integer(psb_epk_), intent(in) :: n
integer(psb_epk_) :: i,j,ix
integer(psb_epk_) :: xx
do j=n-1,1,-1
if (abs(x(j+1)) < abs(x(j))) then
xx = x(j)
ix = idx(j)
i=j+1
do
x(i-1) = x(i)
idx(i-1) = idx(i)
i = i+1
if (i>n) exit
if (abs(x(i)) >= abs(xx)) exit
end do
x(i-1) = xx
idx(i-1) = ix
endif
enddo
end subroutine psi_eaisrx_up
subroutine psi_eaisrx_dw(n,x,idx)
use psb_sort_mod, psb_protect_name => psi_eaisrx_dw
use psb_error_mod
implicit none
integer(psb_epk_), intent(inout) :: x(:)
integer(psb_epk_), intent(inout) :: idx(:)
integer(psb_epk_), intent(in) :: n
integer(psb_epk_) :: i,j,ix
integer(psb_epk_) :: xx
do j=n-1,1,-1
if (abs(x(j+1)) > abs(x(j))) then
xx = x(j)
ix = idx(j)
i=j+1
do
x(i-1) = x(i)
idx(i-1) = idx(i)
i = i+1
if (i>n) exit
if (abs(x(i)) <= abs(xx)) exit
end do
x(i-1) = xx
idx(i-1) = ix
endif
enddo
end subroutine psi_eaisrx_dw
subroutine psi_eaisr_up(n,x)
use psb_sort_mod, psb_protect_name => psi_eaisr_up
use psb_error_mod
implicit none
integer(psb_epk_), intent(inout) :: x(:)
integer(psb_epk_), intent(in) :: n
integer(psb_epk_) :: i,j
integer(psb_epk_) :: xx
do j=n-1,1,-1
if (abs(x(j+1)) < abs(x(j))) then
xx = x(j)
i=j+1
do
x(i-1) = x(i)
i = i+1
if (i>n) exit
if (abs(x(i)) >= abs(xx)) exit
end do
x(i-1) = xx
endif
enddo
end subroutine psi_eaisr_up
subroutine psi_eaisr_dw(n,x)
use psb_sort_mod, psb_protect_name => psi_eaisr_dw
use psb_error_mod
implicit none
integer(psb_epk_), intent(inout) :: x(:)
integer(psb_epk_), intent(in) :: n
integer(psb_epk_) :: i,j
integer(psb_epk_) :: xx
do j=n-1,1,-1
if (abs(x(j+1)) > abs(x(j))) then
xx = x(j)
i=j+1
do
x(i-1) = x(i)
i = i+1
if (i>n) exit
if (abs(x(i)) <= abs(xx)) exit
end do
x(i-1) = xx
endif
enddo
end subroutine psi_eaisr_dw

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

@ -0,0 +1,713 @@
!
! 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:
! D. Knuth
! The Art of Computer Programming, vol. 3
! Addison-Wesley
!
! Aho, Hopcroft, Ullman
! 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
!
! 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
return
end function psb_eisaperm
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
integer(psb_epk_) :: n, k
integer(psb_ipk_) :: err_act
integer(psb_ipk_) :: ierr(5)
character(len=20) :: name
name='psb_msort_u'
call psb_erractionsave(err_act)
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
return
9999 call psb_error_handler(err_act)
return
end subroutine psb_emsort_u
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_) :: 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_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) = 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
return
9999 call psb_error_handler(err_act)
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
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
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
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
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_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
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
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
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
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_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
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
outer: do
if (abs(k(p)) > abs(k(q))) then
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
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_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
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
outer: do
if (abs(k(p)) < abs(k(q))) then
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
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

1318
base/serial/sort/psb_e_qsort_impl.f90
View File

File diff suppressed because it is too large Load Diff

678
base/serial/sort/psb_m_hsort_impl.f90
Unescape Escape View File

@ -0,0 +1,678 @@
!
! 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 and quicksort routines are implemented in the
! serial/aux directory
! References:
! D. Knuth
! The Art of Computer Programming, vol. 3
! Addison-Wesley
!
! Aho, Hopcroft, Ullman
! Data Structures and Algorithms
! Addison-Wesley
!
subroutine psb_mhsort(x,ix,dir,flag)
use psb_sort_mod, psb_protect_name => psb_mhsort
use psb_error_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, i, l, err_act,info
integer(psb_mpk_) :: key
integer(psb_ipk_) :: index
integer(psb_ipk_) :: ierr(5)
character(len=20) :: name
name='psb_hsort'
call psb_erractionsave(err_act)
if (present(flag)) then
flag_ = flag
else
flag_ = psb_sort_ovw_idx_
end if
select case(flag_)
case( psb_sort_ovw_idx_, 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
if (present(dir)) then
dir_ = dir
else
dir_= psb_sort_up_
end if
select case(dir_)
case(psb_sort_up_,psb_sort_down_)
! OK
case (psb_asort_up_,psb_asort_down_)
! OK
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)
!
! Dirty trick to sort with heaps: if we want
! to sort in place upwards, first we set up a heap so that
! we can easily get the LARGEST element, then we take it out
! and put it in the last entry, and so on.
! So, we invert dir_
!
dir_ = -dir_
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 (flag_ == psb_sort_ovw_idx_) then
do i=1, n
ix(i) = i
end do
end if
l = 0
do i=1, n
key = x(i)
index = ix(i)
call psi_idx_insert_heap(key,index,l,x,ix,dir_,info)
if (l /= i) then
write(psb_err_unit,*) 'Mismatch while heapifying ! '
end if
end do
do i=n, 2, -1
call psi_idx_heap_get_first(key,index,l,x,ix,dir_,info)
if (l /= i-1) then
write(psb_err_unit,*) 'Mismatch while pulling out of heap ',l,i
end if
x(i) = key
ix(i) = index
end do
else if (.not.present(ix)) then
l = 0
do i=1, n
key = x(i)
call psi_insert_heap(key,l,x,dir_,info)
if (l /= i) then
write(psb_err_unit,*) 'Mismatch while heapifying ! ',l,i
end if
end do
do i=n, 2, -1
call psi_m_heap_get_first(key,l,x,dir_,info)
if (l /= i-1) then
write(psb_err_unit,*) 'Mismatch while pulling out of heap ',l,i
end if
x(i) = key
end do
end if
return
9999 call psb_error_handler(err_act)
return
end subroutine psb_mhsort
!
! These are packaged so that they can be used to implement
! a heapsort, should the need arise
!
!
! Programming note:
! In the implementation of the heap_get_first function
! we have code like this
!
! if ( ( heap(2*i) < heap(2*i+1) ) .or.&
! & (2*i == last)) then
! j = 2*i
! else
! j = 2*i + 1
! end if
!
! It looks like the 2*i+1 could overflow the array, but this
! is not true because there is a guard statement
! if (i>last/2) exit
! and because last has just been reduced by 1 when defining the return value,
! therefore 2*i+1 may be greater than the current value of last,
! but cannot be greater than the value of last when the routine was entered
! hence it is safe.
!
!
!
subroutine psi_m_insert_heap(key,last,heap,dir,info)
use psb_sort_mod, psb_protect_name => psi_m_insert_heap
implicit none
!
! Input:
! key: the new value
! last: pointer to the last occupied element in heap
! heap: the heap
! dir: sorting direction
integer(psb_mpk_), intent(in) :: key
integer(psb_ipk_), intent(in) :: dir
integer(psb_mpk_), intent(inout) :: heap(:)
integer(psb_ipk_), intent(inout) :: last
integer(psb_ipk_), intent(out) :: info
integer(psb_ipk_) :: i, i2
integer(psb_mpk_) :: temp
info = psb_success_
if (last < 0) then
write(psb_err_unit,*) 'Invalid last in heap ',last
info = last
return
endif
last = last + 1
if (last > size(heap)) then
write(psb_err_unit,*) 'out of bounds '
info = -1
return
end if
i = last
heap(i) = key
select case(dir)
case (psb_sort_up_)
do
if (i<=1) exit
i2 = i/2
if (heap(i) < heap(i2)) then
temp = heap(i)
heap(i) = heap(i2)
heap(i2) = temp
i = i2
else
exit
end if
end do
case (psb_sort_down_)
do
if (i<=1) exit
i2 = i/2
if (heap(i) > heap(i2)) then
temp = heap(i)
heap(i) = heap(i2)
heap(i2) = temp
i = i2
else
exit
end if
end do
case (psb_asort_up_)
do
if (i<=1) exit
i2 = i/2
if (abs(heap(i)) < abs(heap(i2))) then
temp = heap(i)
heap(i) = heap(i2)
heap(i2) = temp
i = i2
else
exit
end if
end do
case (psb_asort_down_)
do
if (i<=1) exit
i2 = i/2
if (abs(heap(i)) > abs(heap(i2))) then
temp = heap(i)
heap(i) = heap(i2)
heap(i2) = temp
i = i2
else
exit
end if
end do
case default
write(psb_err_unit,*) 'Invalid direction in heap ',dir
end select
return
end subroutine psi_m_insert_heap
subroutine psi_m_heap_get_first(key,last,heap,dir,info)
use psb_sort_mod, psb_protect_name => psi_m_heap_get_first
implicit none
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
integer(psb_ipk_) :: i, j
integer(psb_mpk_) :: temp
info = psb_success_
if (last <= 0) then
key = 0
info = -1
return
endif
key = heap(1)
heap(1) = heap(last)
last = last - 1
select case(dir)
case (psb_sort_up_)
i = 1
do
if (i > (last/2)) exit
if ( (heap(2*i) < heap(2*i+1)) .or.&
& (2*i == last)) then
j = 2*i
else
j = 2*i + 1
end if
if (heap(i) > heap(j)) then
temp = heap(i)
heap(i) = heap(j)
heap(j) = temp
i = j
else
exit
end if
end do
case (psb_sort_down_)
i = 1
do
if (i > (last/2)) exit
if ( (heap(2*i) > heap(2*i+1)) .or.&
& (2*i == last)) then
j = 2*i
else
j = 2*i + 1
end if
if (heap(i) < heap(j)) then
temp = heap(i)
heap(i) = heap(j)
heap(j) = temp
i = j
else
exit
end if
end do
case (psb_asort_up_)
i = 1
do
if (i > (last/2)) exit
if ( (abs(heap(2*i)) < abs(heap(2*i+1))) .or.&
& (2*i == last)) then
j = 2*i
else
j = 2*i + 1
end if
if (abs(heap(i)) > abs(heap(j))) then
temp = heap(i)
heap(i) = heap(j)
heap(j) = temp
i = j
else
exit
end if
end do
case (psb_asort_down_)
i = 1
do
if (i > (last/2)) exit
if ( (abs(heap(2*i)) > abs(heap(2*i+1))) .or.&
& (2*i == last)) then
j = 2*i
else
j = 2*i + 1
end if
if (abs(heap(i)) < abs(heap(j))) then
temp = heap(i)
heap(i) = heap(j)
heap(j) = temp
i = j
else
exit
end if
end do
case default
write(psb_err_unit,*) 'Invalid direction in heap ',dir
end select
return
end subroutine psi_m_heap_get_first
subroutine psi_m_idx_insert_heap(key,index,last,heap,idxs,dir,info)
use psb_sort_mod, psb_protect_name => psi_m_idx_insert_heap
implicit none
!
! Input:
! key: the new value
! index: the new index
! last: pointer to the last occupied element in heap
! heap: the heap
! idxs: the indices
! dir: sorting direction
integer(psb_mpk_), intent(in) :: key
integer(psb_ipk_), intent(in) :: index,dir
integer(psb_mpk_), intent(inout) :: heap(:)
integer(psb_ipk_), intent(inout) :: idxs(:),last
integer(psb_ipk_), intent(out) :: info
integer(psb_ipk_) :: i, i2, itemp
integer(psb_mpk_) :: temp
info = psb_success_
if (last < 0) then
write(psb_err_unit,*) 'Invalid last in heap ',last
info = last
return
endif
last = last + 1
if (last > size(heap)) then
write(psb_err_unit,*) 'out of bounds '
info = -1
return
end if
i = last
heap(i) = key
idxs(i) = index
select case(dir)
case (psb_sort_up_)
do
if (i<=1) exit
i2 = i/2
if (heap(i) < heap(i2)) then
itemp = idxs(i)
idxs(i) = idxs(i2)
idxs(i2) = itemp
temp = heap(i)
heap(i) = heap(i2)
heap(i2) = temp
i = i2
else
exit
end if
end do
case (psb_sort_down_)
do
if (i<=1) exit
i2 = i/2
if (heap(i) > heap(i2)) then
itemp = idxs(i)
idxs(i) = idxs(i2)
idxs(i2) = itemp
temp = heap(i)
heap(i) = heap(i2)
heap(i2) = temp
i = i2
else
exit
end if
end do
case (psb_asort_up_)
do
if (i<=1) exit
i2 = i/2
if (abs(heap(i)) < abs(heap(i2))) then
itemp = idxs(i)
idxs(i) = idxs(i2)
idxs(i2) = itemp
temp = heap(i)
heap(i) = heap(i2)
heap(i2) = temp
i = i2
else
exit
end if
end do
case (psb_asort_down_)
do
if (i<=1) exit
i2 = i/2
if (abs(heap(i)) > abs(heap(i2))) then
itemp = idxs(i)
idxs(i) = idxs(i2)
idxs(i2) = itemp
temp = heap(i)
heap(i) = heap(i2)
heap(i2) = temp
i = i2
else
exit
end if
end do
case default
write(psb_err_unit,*) 'Invalid direction in heap ',dir
end select
return
end subroutine psi_m_idx_insert_heap
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) :: heap(:)
integer(psb_ipk_), intent(out) :: index,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
info = psb_success_
if (last <= 0) then
key = 0
index = 0
info = -1
return
endif
key = heap(1)
index = idxs(1)
heap(1) = heap(last)
idxs(1) = idxs(last)
last = last - 1
select case(dir)
case (psb_sort_up_)
i = 1
do
if (i > (last/2)) exit
if ( (heap(2*i) < heap(2*i+1)) .or.&
& (2*i == last)) then
j = 2*i
else
j = 2*i + 1
end if
if (heap(i) > heap(j)) then
itemp = idxs(i)
idxs(i) = idxs(j)
idxs(j) = itemp
temp = heap(i)
heap(i) = heap(j)
heap(j) = temp
i = j
else
exit
end if
end do
case (psb_sort_down_)
i = 1
do
if (i > (last/2)) exit
if ( (heap(2*i) > heap(2*i+1)) .or.&
& (2*i == last)) then
j = 2*i
else
j = 2*i + 1
end if
if (heap(i) < heap(j)) then
itemp = idxs(i)
idxs(i) = idxs(j)
idxs(j) = itemp
temp = heap(i)
heap(i) = heap(j)
heap(j) = temp
i = j
else
exit
end if
end do
case (psb_asort_up_)
i = 1
do
if (i > (last/2)) exit
if ( (abs(heap(2*i)) < abs(heap(2*i+1))) .or.&
& (2*i == last)) then
j = 2*i
else
j = 2*i + 1
end if
if (abs(heap(i)) > abs(heap(j))) then
itemp = idxs(i)
idxs(i) = idxs(j)
idxs(j) = itemp
temp = heap(i)
heap(i) = heap(j)
heap(j) = temp
i = j
else
exit
end if
end do
case (psb_asort_down_)
i = 1
do
if (i > (last/2)) exit
if ( (abs(heap(2*i)) > abs(heap(2*i+1))) .or.&
& (2*i == last)) then
j = 2*i
else
j = 2*i + 1
end if
if (abs(heap(i)) < abs(heap(j))) then
itemp = idxs(i)
idxs(i) = idxs(j)
idxs(j) = itemp
temp = heap(i)
heap(i) = heap(j)
heap(j) = temp
i = j
else
exit
end if
end do
case default
write(psb_err_unit,*) 'Invalid direction in heap ',dir
end select
return
end subroutine psi_m_idx_heap_get_first

341
base/serial/sort/psb_m_isort_impl.f90
Unescape Escape View File

@ -0,0 +1,341 @@
!
! 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 insertion sort routines
! References:
! D. Knuth
! The Art of Computer Programming, vol. 3
! Addison-Wesley
!
! Aho, Hopcroft, Ullman
! Data Structures and Algorithms
! Addison-Wesley
!
subroutine psb_misort(x,ix,dir,flag)
use psb_sort_mod, psb_protect_name => psb_misort
use psb_error_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_, err_act
integer(psb_ipk_) :: n, i
integer(psb_ipk_) :: ierr(5)
character(len=20) :: name
name='psb_misort'
call psb_erractionsave(err_act)
if (present(flag)) then
flag_ = flag
else
flag_ = psb_sort_ovw_idx_
end if
select case(flag_)
case( psb_sort_ovw_idx_, 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
if (present(dir)) then
dir_ = dir
else
dir_= psb_sort_up_
end if
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 (flag_==psb_sort_ovw_idx_) then
do i=1,n
ix(i) = i
end do
end if
select case(dir_)
case (psb_sort_up_)
call psi_misrx_up(n,x,ix)
case (psb_sort_down_)
call psi_misrx_dw(n,x,ix)
case (psb_asort_up_)
call psi_maisrx_up(n,x,ix)
case (psb_asort_down_)
call psi_maisrx_dw(n,x,ix)
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
else
select case(dir_)
case (psb_sort_up_)
call psi_misr_up(n,x)
case (psb_sort_down_)
call psi_misr_dw(n,x)
case (psb_asort_up_)
call psi_maisr_up(n,x)
case (psb_asort_down_)
call psi_maisr_dw(n,x)
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
end if
return
9999 call psb_error_handler(err_act)
return
end subroutine psb_misort
subroutine psi_misrx_up(n,x,idx)
use psb_sort_mod, psb_protect_name => psi_misrx_up
use psb_error_mod
implicit none
integer(psb_mpk_), intent(inout) :: x(:)
integer(psb_ipk_), intent(inout) :: idx(:)
integer(psb_ipk_), intent(in) :: n
integer(psb_ipk_) :: i,j,ix
integer(psb_mpk_) :: xx
do j=n-1,1,-1
if (x(j+1) < x(j)) then
xx = x(j)
ix = idx(j)
i=j+1
do
x(i-1) = x(i)
idx(i-1) = idx(i)
i = i+1
if (i>n) exit
if (x(i) >= xx) exit
end do
x(i-1) = xx
idx(i-1) = ix
endif
enddo
end subroutine psi_misrx_up
subroutine psi_misrx_dw(n,x,idx)
use psb_sort_mod, psb_protect_name => psi_misrx_dw
use psb_error_mod
implicit none
integer(psb_mpk_), intent(inout) :: x(:)
integer(psb_ipk_), intent(inout) :: idx(:)
integer(psb_ipk_), intent(in) :: n
integer(psb_ipk_) :: i,j,ix
integer(psb_mpk_) :: xx
do j=n-1,1,-1
if (x(j+1) > x(j)) then
xx = x(j)
ix = idx(j)
i=j+1
do
x(i-1) = x(i)
idx(i-1) = idx(i)
i = i+1
if (i>n) exit
if (x(i) <= xx) exit
end do
x(i-1) = xx
idx(i-1) = ix
endif
enddo
end subroutine psi_misrx_dw
subroutine psi_misr_up(n,x)
use psb_sort_mod, psb_protect_name => psi_misr_up
use psb_error_mod
implicit none
integer(psb_mpk_), intent(inout) :: x(:)
integer(psb_ipk_), intent(in) :: n
integer(psb_ipk_) :: i,j
integer(psb_mpk_) :: xx
do j=n-1,1,-1
if (x(j+1) < x(j)) then
xx = x(j)
i=j+1
do
x(i-1) = x(i)
i = i+1
if (i>n) exit
if (x(i) >= xx) exit
end do
x(i-1) = xx
endif
enddo
end subroutine psi_misr_up
subroutine psi_misr_dw(n,x)
use psb_sort_mod, psb_protect_name => psi_misr_dw
use psb_error_mod
implicit none
integer(psb_mpk_), intent(inout) :: x(:)
integer(psb_ipk_), intent(in) :: n
integer(psb_ipk_) :: i,j
integer(psb_mpk_) :: xx
do j=n-1,1,-1
if (x(j+1) > x(j)) then
xx = x(j)
i=j+1
do
x(i-1) = x(i)
i = i+1
if (i>n) exit
if (x(i) <= xx) exit
end do
x(i-1) = xx
endif
enddo
end subroutine psi_misr_dw
subroutine psi_maisrx_up(n,x,idx)
use psb_sort_mod, psb_protect_name => psi_maisrx_up
use psb_error_mod
implicit none
integer(psb_mpk_), intent(inout) :: x(:)
integer(psb_ipk_), intent(inout) :: idx(:)
integer(psb_ipk_), intent(in) :: n
integer(psb_ipk_) :: i,j,ix
integer(psb_mpk_) :: xx
do j=n-1,1,-1
if (abs(x(j+1)) < abs(x(j))) then
xx = x(j)
ix = idx(j)
i=j+1
do
x(i-1) = x(i)
idx(i-1) = idx(i)
i = i+1
if (i>n) exit
if (abs(x(i)) >= abs(xx)) exit
end do
x(i-1) = xx
idx(i-1) = ix
endif
enddo
end subroutine psi_maisrx_up
subroutine psi_maisrx_dw(n,x,idx)
use psb_sort_mod, psb_protect_name => psi_maisrx_dw
use psb_error_mod
implicit none
integer(psb_mpk_), intent(inout) :: x(:)
integer(psb_ipk_), intent(inout) :: idx(:)
integer(psb_ipk_), intent(in) :: n
integer(psb_ipk_) :: i,j,ix
integer(psb_mpk_) :: xx
do j=n-1,1,-1
if (abs(x(j+1)) > abs(x(j))) then
xx = x(j)
ix = idx(j)
i=j+1
do
x(i-1) = x(i)
idx(i-1) = idx(i)
i = i+1
if (i>n) exit
if (abs(x(i)) <= abs(xx)) exit
end do
x(i-1) = xx
idx(i-1) = ix
endif
enddo
end subroutine psi_maisrx_dw
subroutine psi_maisr_up(n,x)
use psb_sort_mod, psb_protect_name => psi_maisr_up
use psb_error_mod
implicit none
integer(psb_mpk_), intent(inout) :: x(:)
integer(psb_ipk_), intent(in) :: n
integer(psb_ipk_) :: i,j
integer(psb_mpk_) :: xx
do j=n-1,1,-1
if (abs(x(j+1)) < abs(x(j))) then
xx = x(j)
i=j+1
do
x(i-1) = x(i)
i = i+1
if (i>n) exit
if (abs(x(i)) >= abs(xx)) exit
end do
x(i-1) = xx
endif
enddo
end subroutine psi_maisr_up
subroutine psi_maisr_dw(n,x)
use psb_sort_mod, psb_protect_name => psi_maisr_dw
use psb_error_mod
implicit none
integer(psb_mpk_), intent(inout) :: x(:)
integer(psb_ipk_), intent(in) :: n
integer(psb_ipk_) :: i,j
integer(psb_mpk_) :: xx
do j=n-1,1,-1
if (abs(x(j+1)) > abs(x(j))) then
xx = x(j)
i=j+1
do
x(i-1) = x(i)
i = i+1
if (i>n) exit
if (abs(x(i)) <= abs(xx)) exit
end do
x(i-1) = xx
endif
enddo
end subroutine psi_maisr_dw

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

@ -0,0 +1,713 @@
!
! 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:
! D. Knuth
! The Art of Computer Programming, vol. 3
! Addison-Wesley
!
! Aho, Hopcroft, Ullman
! 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
!
! 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
return
end function psb_misaperm
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
integer(psb_ipk_) :: n, k
integer(psb_ipk_) :: err_act
integer(psb_ipk_) :: ierr(5)
character(len=20) :: name
name='psb_msort_u'
call psb_erractionsave(err_act)
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
return
9999 call psb_error_handler(err_act)
return
end subroutine psb_mmsort_u
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_) :: 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_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) = 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
return
9999 call psb_error_handler(err_act)
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
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
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
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
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_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
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
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
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
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_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
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
outer: do
if (abs(k(p)) > abs(k(q))) then
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
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_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
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
outer: do
if (abs(k(p)) < abs(k(q))) then
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
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

1318
base/serial/sort/psb_m_qsort_impl.f90
View File

File diff suppressed because it is too large Load Diff
Write Preview
Loading…
Cancel
Save