diff --git a/base/serial/sort/psb_i_hsort_impl.f90 b/base/serial/sort/psb_i_hsort_impl.f90 deleted file mode 100644 index 1c53650a..00000000 --- a/base/serial/sort/psb_i_hsort_impl.f90 +++ /dev/null @@ -1,678 +0,0 @@ -! -! 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_ihsort(x,ix,dir,flag) - use psb_sort_mod, psb_protect_name => psb_ihsort - use psb_error_mod - implicit none - integer(psb_ipk_), 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_ipk_) :: 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_i_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_ihsort - - - -! -! 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_i_insert_heap(key,last,heap,dir,info) - use psb_sort_mod, psb_protect_name => psi_i_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_ipk_), intent(in) :: key - integer(psb_ipk_), intent(in) :: dir - integer(psb_ipk_), intent(inout) :: heap(:) - integer(psb_ipk_), intent(inout) :: last - integer(psb_ipk_), intent(out) :: info - integer(psb_ipk_) :: i, i2 - integer(psb_ipk_) :: 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_i_insert_heap - - -subroutine psi_i_heap_get_first(key,last,heap,dir,info) - use psb_sort_mod, psb_protect_name => psi_i_heap_get_first - implicit none - - integer(psb_ipk_), intent(inout) :: key - integer(psb_ipk_), intent(inout) :: last - integer(psb_ipk_), intent(in) :: dir - integer(psb_ipk_), intent(inout) :: heap(:) - integer(psb_ipk_), intent(out) :: info - - integer(psb_ipk_) :: i, j - integer(psb_ipk_) :: 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_i_heap_get_first - - -subroutine psi_i_idx_insert_heap(key,index,last,heap,idxs,dir,info) - use psb_sort_mod, psb_protect_name => psi_i_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_ipk_), intent(in) :: key - integer(psb_ipk_), intent(in) :: index,dir - integer(psb_ipk_), intent(inout) :: heap(:) - integer(psb_ipk_), intent(inout) :: idxs(:),last - integer(psb_ipk_), intent(out) :: info - integer(psb_ipk_) :: i, i2, itemp - integer(psb_ipk_) :: 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_i_idx_insert_heap - -subroutine psi_i_idx_heap_get_first(key,index,last,heap,idxs,dir,info) - use psb_sort_mod, psb_protect_name => psi_i_idx_heap_get_first - implicit none - - integer(psb_ipk_), 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_ipk_), intent(out) :: key - - integer(psb_ipk_) :: i, j,itemp - integer(psb_ipk_) :: 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_i_idx_heap_get_first - - - - diff --git a/base/serial/sort/psb_i_isort_impl.f90 b/base/serial/sort/psb_i_isort_impl.f90 deleted file mode 100644 index 41c1381f..00000000 --- a/base/serial/sort/psb_i_isort_impl.f90 +++ /dev/null @@ -1,341 +0,0 @@ -! -! 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_iisort(x,ix,dir,flag) - use psb_sort_mod, psb_protect_name => psb_iisort - use psb_error_mod - implicit none - integer(psb_ipk_), 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_iisort' - 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_iisrx_up(n,x,ix) - case (psb_sort_down_) - call psi_iisrx_dw(n,x,ix) - case (psb_asort_up_) - call psi_iaisrx_up(n,x,ix) - case (psb_asort_down_) - call psi_iaisrx_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_iisr_up(n,x) - case (psb_sort_down_) - call psi_iisr_dw(n,x) - case (psb_asort_up_) - call psi_iaisr_up(n,x) - case (psb_asort_down_) - call psi_iaisr_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_iisort - -subroutine psi_iisrx_up(n,x,idx) - use psb_sort_mod, psb_protect_name => psi_iisrx_up - use psb_error_mod - implicit none - integer(psb_ipk_), intent(inout) :: x(:) - integer(psb_ipk_), intent(inout) :: idx(:) - integer(psb_ipk_), intent(in) :: n - integer(psb_ipk_) :: i,j,ix - integer(psb_ipk_) :: 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_iisrx_up - -subroutine psi_iisrx_dw(n,x,idx) - use psb_sort_mod, psb_protect_name => psi_iisrx_dw - use psb_error_mod - implicit none - integer(psb_ipk_), intent(inout) :: x(:) - integer(psb_ipk_), intent(inout) :: idx(:) - integer(psb_ipk_), intent(in) :: n - integer(psb_ipk_) :: i,j,ix - integer(psb_ipk_) :: 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_iisrx_dw - - -subroutine psi_iisr_up(n,x) - use psb_sort_mod, psb_protect_name => psi_iisr_up - use psb_error_mod - implicit none - integer(psb_ipk_), intent(inout) :: x(:) - integer(psb_ipk_), intent(in) :: n - integer(psb_ipk_) :: i,j - integer(psb_ipk_) :: 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_iisr_up - -subroutine psi_iisr_dw(n,x) - use psb_sort_mod, psb_protect_name => psi_iisr_dw - use psb_error_mod - implicit none - integer(psb_ipk_), intent(inout) :: x(:) - integer(psb_ipk_), intent(in) :: n - integer(psb_ipk_) :: i,j - integer(psb_ipk_) :: 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_iisr_dw - -subroutine psi_iaisrx_up(n,x,idx) - use psb_sort_mod, psb_protect_name => psi_iaisrx_up - use psb_error_mod - implicit none - integer(psb_ipk_), intent(inout) :: x(:) - integer(psb_ipk_), intent(inout) :: idx(:) - integer(psb_ipk_), intent(in) :: n - integer(psb_ipk_) :: i,j,ix - integer(psb_ipk_) :: 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_iaisrx_up - -subroutine psi_iaisrx_dw(n,x,idx) - use psb_sort_mod, psb_protect_name => psi_iaisrx_dw - use psb_error_mod - implicit none - integer(psb_ipk_), intent(inout) :: x(:) - integer(psb_ipk_), intent(inout) :: idx(:) - integer(psb_ipk_), intent(in) :: n - integer(psb_ipk_) :: i,j,ix - integer(psb_ipk_) :: 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_iaisrx_dw - -subroutine psi_iaisr_up(n,x) - use psb_sort_mod, psb_protect_name => psi_iaisr_up - use psb_error_mod - implicit none - integer(psb_ipk_), intent(inout) :: x(:) - integer(psb_ipk_), intent(in) :: n - integer(psb_ipk_) :: i,j - integer(psb_ipk_) :: 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_iaisr_up - -subroutine psi_iaisr_dw(n,x) - use psb_sort_mod, psb_protect_name => psi_iaisr_dw - use psb_error_mod - implicit none - integer(psb_ipk_), intent(inout) :: x(:) - integer(psb_ipk_), intent(in) :: n - integer(psb_ipk_) :: i,j - integer(psb_ipk_) :: 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_iaisr_dw - diff --git a/base/serial/sort/psb_i_msort_impl.f90 b/base/serial/sort/psb_i_msort_impl.f90 deleted file mode 100644 index 1e9ad9ec..00000000 --- a/base/serial/sort/psb_i_msort_impl.f90 +++ /dev/null @@ -1,713 +0,0 @@ -! -! 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_iisaperm(n,eip) - use psb_sort_mod, psb_protect_name => psb_iisaperm - implicit none - - integer(psb_ipk_), intent(in) :: n - integer(psb_ipk_), intent(in) :: eip(n) - integer(psb_ipk_), allocatable :: ip(:) - integer(psb_ipk_) :: i,j,m, info - - - psb_iisaperm = .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_iisaperm = .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_iisaperm = .false. - goto 9999 - endif - end if - enddo -9999 continue - - return - end function psb_iisaperm - - - subroutine psb_imsort_u(x,nout,dir) - use psb_sort_mod, psb_protect_name => psb_imsort_u - use psb_error_mod - implicit none - integer(psb_ipk_), 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_imsort_u - - - function psb_ibsrch(key,n,v) result(ipos) - use psb_sort_mod, psb_protect_name => psb_ibsrch - implicit none - integer(psb_ipk_) :: ipos, n - integer(psb_ipk_) :: key - integer(psb_ipk_) :: 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_ibsrch - - function psb_issrch(key,n,v) result(ipos) - use psb_sort_mod, psb_protect_name => psb_issrch - implicit none - integer(psb_ipk_) :: ipos, n - integer(psb_ipk_) :: key - integer(psb_ipk_) :: 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_issrch - - subroutine psb_imsort(x,ix,dir,flag) - use psb_sort_mod, psb_protect_name => psb_imsort - use psb_error_mod - use psb_ip_reord_mod - implicit none - integer(psb_ipk_), 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_imsort' - 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_i_msort') - goto 9999 - endif - - select case(dir_) - case (psb_sort_up_) - call psi_i_msort_up(n,x,iaux,iret) - case (psb_sort_down_) - call psi_i_msort_dw(n,x,iaux,iret) - case (psb_asort_up_) - call psi_i_amsort_up(n,x,iaux,iret) - case (psb_asort_down_) - call psi_i_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_imsort - - subroutine psi_i_msort_up(n,k,l,iret) - use psb_const_mod - implicit none - integer(psb_ipk_) :: n, iret - integer(psb_ipk_) :: 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_i_msort_up - - subroutine psi_i_msort_dw(n,k,l,iret) - use psb_const_mod - implicit none - integer(psb_ipk_) :: n, iret - integer(psb_ipk_) :: 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_i_msort_dw - - subroutine psi_i_amsort_up(n,k,l,iret) - use psb_const_mod - implicit none - integer(psb_ipk_) :: n, iret - integer(psb_ipk_) :: 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_i_amsort_up - - subroutine psi_i_amsort_dw(n,k,l,iret) - use psb_const_mod - implicit none - integer(psb_ipk_) :: n, iret - integer(psb_ipk_) :: 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_i_amsort_dw - - - - - - - - diff --git a/base/serial/sort/psb_i_qsort_impl.f90 b/base/serial/sort/psb_i_qsort_impl.f90 deleted file mode 100644 index e5f0ac92..00000000 --- a/base/serial/sort/psb_i_qsort_impl.f90 +++ /dev/null @@ -1,1318 +0,0 @@ -! -! 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 quicksort routines -! References: -! D. Knuth -! The Art of Computer Programming, vol. 3 -! Addison-Wesley -! -! Aho, Hopcroft, Ullman -! Data Structures and Algorithms -! Addison-Wesley -! -subroutine psb_iqsort(x,ix,dir,flag) - use psb_sort_mod, psb_protect_name => psb_iqsort - use psb_error_mod - implicit none - integer(psb_ipk_), 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, i - integer(psb_ipk_) :: n - integer(psb_ipk_) :: ierr(5) - character(len=20) :: name - - name='psb_iqsort' - 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_iqsrx_up(n,x,ix) - case (psb_sort_down_) - call psi_iqsrx_dw(n,x,ix) - case (psb_asort_up_) - call psi_iaqsrx_up(n,x,ix) - case (psb_asort_down_) - call psi_iaqsrx_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_iqsr_up(n,x) - case (psb_sort_down_) - call psi_iqsr_dw(n,x) - case (psb_asort_up_) - call psi_iaqsr_up(n,x) - case (psb_asort_down_) - call psi_iaqsr_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_iqsort - -subroutine psi_iqsrx_up(n,x,idx) - use psb_sort_mod, psb_protect_name => psi_iqsrx_up - use psb_error_mod - implicit none - - integer(psb_ipk_), intent(inout) :: x(:) - integer(psb_ipk_), intent(inout) :: idx(:) - integer(psb_ipk_), intent(in) :: n - ! .. Local Scalars .. - integer(psb_ipk_) :: piv, xk, xt - integer(psb_ipk_) :: i, j, ilx, iux, istp, lpiv - integer(psb_ipk_) :: n1, n2 - integer(psb_ipk_) :: ixt - integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=72 - integer(psb_ipk_) :: istack(nparms,maxstack) - - if (n > ithrs) then - ! - ! Init stack pointer - ! - istp = 1 - istack(1,istp) = 1 - istack(2,istp) = n - - do - if (istp <= 0) exit - ilx = istack(1,istp) - iux = istack(2,istp) - istp = istp - 1 - ! - ! Choose a pivot with median-of-three heuristics, leave it - ! in the LPIV location - ! - i = ilx - j = iux - lpiv = (i+j)/2 - piv = x(lpiv) - if (piv < x(i)) then - xt = x(i) - ixt = idx(i) - x(i) = x(lpiv) - idx(i) = idx(lpiv) - x(lpiv) = xt - idx(lpiv) = ixt - piv = x(lpiv) - endif - if (piv > x(j)) then - xt = x(j) - ixt = idx(j) - x(j) = x(lpiv) - idx(j) = idx(lpiv) - x(lpiv) = xt - idx(lpiv) = ixt - piv = x(lpiv) - endif - if (piv < x(i)) then - xt = x(i) - ixt = idx(i) - x(i) = x(lpiv) - idx(i) = idx(lpiv) - x(lpiv) = xt - idx(lpiv) = ixt - piv = x(lpiv) - endif - ! - ! now piv is correct; place it into first location - xt = x(i) - ixt = idx(i) - x(i) = x(lpiv) - idx(i) = idx(lpiv) - x(lpiv) = xt - idx(lpiv) = ixt - piv = x(lpiv) - - i = ilx - 1 - j = iux + 1 - - outer_up: do - in_up1: do - i = i + 1 - xk = x(i) - if (xk >= piv) exit in_up1 - end do in_up1 - ! - ! Ensure finite termination for next loop - ! - xt = xk - x(i) = piv - in_up2:do - j = j - 1 - xk = x(j) - if (xk <= piv) exit in_up2 - end do in_up2 - x(i) = xt - - if (j > i) then - xt = x(i) - ixt = idx(i) - x(i) = x(j) - idx(i) = idx(j) - x(j) = xt - idx(j) = ixt - else - exit outer_up - end if - end do outer_up - if (i == ilx) then - if (x(i) /= piv) then - call psb_errpush(psb_err_internal_error_,& - & r_name='psi_iqsrx',a_err='impossible pivot condition') - call psb_error() - endif - i = i + 1 - endif - - n1 = (i-1)-ilx+1 - n2 = iux-(i)+1 - if (n1 > n2) then - if (n1 > ithrs) then - istp = istp + 1 - istack(1,istp) = ilx - istack(2,istp) = i-1 - else - call psi_iisrx_up(n1,x(ilx:i-1),idx(ilx:i-1)) - endif - if (n2 > ithrs) then - istp = istp + 1 - istack(1,istp) = i - istack(2,istp) = iux - else - call psi_iisrx_up(n2,x(i:iux),idx(i:iux)) - endif - else - if (n2 > ithrs) then - istp = istp + 1 - istack(1,istp) = i - istack(2,istp) = iux - else - call psi_iisrx_up(n2,x(i:iux),idx(i:iux)) - endif - if (n1 > ithrs) then - istp = istp + 1 - istack(1,istp) = ilx - istack(2,istp) = i-1 - else - call psi_iisrx_up(n1,x(ilx:i-1),idx(ilx:i-1)) - endif - endif - enddo - else - call psi_iisrx_up(n,x,idx) - endif -end subroutine psi_iqsrx_up - -subroutine psi_iqsrx_dw(n,x,idx) - use psb_sort_mod, psb_protect_name => psi_iqsrx_dw - use psb_error_mod - implicit none - - integer(psb_ipk_), intent(inout) :: x(:) - integer(psb_ipk_), intent(inout) :: idx(:) - integer(psb_ipk_), intent(in) :: n - ! .. Local Scalars .. - integer(psb_ipk_) :: piv, xk, xt - integer(psb_ipk_) :: i, j, ilx, iux, istp, lpiv - integer(psb_ipk_) :: n1, n2 - integer(psb_ipk_) :: ixt - - integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=72 - integer(psb_ipk_) :: istack(nparms,maxstack) - - if (n > ithrs) then - ! - ! Init stack pointer - ! - istp = 1 - istack(1,istp) = 1 - istack(2,istp) = n - - do - if (istp <= 0) exit - ilx = istack(1,istp) - iux = istack(2,istp) - istp = istp - 1 - ! - ! Choose a pivot with median-of-three heuristics, leave it - ! in the LPIV location - ! - i = ilx - j = iux - lpiv = (i+j)/2 - piv = x(lpiv) - if (piv > x(i)) then - xt = x(i) - ixt = idx(i) - x(i) = x(lpiv) - idx(i) = idx(lpiv) - x(lpiv) = xt - idx(lpiv) = ixt - piv = x(lpiv) - endif - if (piv < x(j)) then - xt = x(j) - ixt = idx(j) - x(j) = x(lpiv) - idx(j) = idx(lpiv) - x(lpiv) = xt - idx(lpiv) = ixt - piv = x(lpiv) - endif - if (piv > x(i)) then - xt = x(i) - ixt = idx(i) - x(i) = x(lpiv) - idx(i) = idx(lpiv) - x(lpiv) = xt - idx(lpiv) = ixt - piv = x(lpiv) - endif - ! - ! now piv is correct; place it into first location - xt = x(i) - ixt = idx(i) - x(i) = x(lpiv) - idx(i) = idx(lpiv) - x(lpiv) = xt - idx(lpiv) = ixt - piv = x(lpiv) - - i = ilx - 1 - j = iux + 1 - - outer_dw: do - in_dw1: do - i = i + 1 - xk = x(i) - if (xk <= piv) exit in_dw1 - end do in_dw1 - ! - ! Ensure finite termination for next loop - ! - xt = xk - x(i) = piv - in_dw2:do - j = j - 1 - xk = x(j) - if (xk >= piv) exit in_dw2 - end do in_dw2 - x(i) = xt - - if (j > i) then - xt = x(i) - ixt = idx(i) - x(i) = x(j) - idx(i) = idx(j) - x(j) = xt - idx(j) = ixt - else - exit outer_dw - end if - end do outer_dw - if (i == ilx) then - if (x(i) /= piv) then - call psb_errpush(psb_err_internal_error_,& - & r_name='psi_iqsrx',a_err='impossible pivot condition') - call psb_error() - endif - i = i + 1 - endif - - n1 = (i-1)-ilx+1 - n2 = iux-(i)+1 - if (n1 > n2) then - if (n1 > ithrs) then - istp = istp + 1 - istack(1,istp) = ilx - istack(2,istp) = i-1 - else - call psi_iisrx_dw(n1,x(ilx:i-1),idx(ilx:i-1)) - endif - if (n2 > ithrs) then - istp = istp + 1 - istack(1,istp) = i - istack(2,istp) = iux - else - call psi_iisrx_dw(n2,x(i:iux),idx(i:iux)) - endif - else - if (n2 > ithrs) then - istp = istp + 1 - istack(1,istp) = i - istack(2,istp) = iux - else - call psi_iisrx_dw(n2,x(i:iux),idx(i:iux)) - endif - if (n1 > ithrs) then - istp = istp + 1 - istack(1,istp) = ilx - istack(2,istp) = i-1 - else - call psi_iisrx_dw(n1,x(ilx:i-1),idx(ilx:i-1)) - endif - endif - enddo - else - call psi_iisrx_dw(n,x,idx) - endif - -end subroutine psi_iqsrx_dw - -subroutine psi_iqsr_up(n,x) - use psb_sort_mod, psb_protect_name => psi_iqsr_up - use psb_error_mod - implicit none - - integer(psb_ipk_), intent(inout) :: x(:) - integer(psb_ipk_), intent(in) :: n - ! .. - ! .. Local Scalars .. - integer(psb_ipk_) :: piv, xt, xk - integer(psb_ipk_) :: i, j, ilx, iux, istp, lpiv - integer(psb_ipk_) :: n1, n2 - - integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=72 - integer(psb_ipk_) :: istack(nparms,maxstack) - - - if (n > ithrs) then - ! - ! Init stack pointer - ! - istp = 1 - istack(1,istp) = 1 - istack(2,istp) = n - - do - if (istp <= 0) exit - ilx = istack(1,istp) - iux = istack(2,istp) - istp = istp - 1 - ! - ! Choose a pivot with median-of-three heuristics, leave it - ! in the LPIV location - ! - i = ilx - j = iux - lpiv = (i+j)/2 - piv = x(lpiv) - if (piv < x(i)) then - xt = x(i) - x(i) = x(lpiv) - x(lpiv) = xt - piv = x(lpiv) - endif - if (piv > x(j)) then - xt = x(j) - x(j) = x(lpiv) - x(lpiv) = xt - piv = x(lpiv) - endif - if (piv < x(i)) then - xt = x(i) - x(i) = x(lpiv) - x(lpiv) = xt - piv = x(lpiv) - endif - ! - ! now piv is correct; place it into first location - - xt = x(i) - x(i) = x(lpiv) - x(lpiv) = xt - - i = ilx - 1 - j = iux + 1 - - outer_up: do - in_up1: do - i = i + 1 - xk = x(i) - if (xk >= piv) exit in_up1 - end do in_up1 - ! - ! Ensure finite termination for next loop - ! - xt = xk - x(i) = piv - in_up2:do - j = j - 1 - xk = x(j) - if (xk <= piv) exit in_up2 - end do in_up2 - x(i) = xt - - if (j > i) then - xt = x(i) - x(i) = x(j) - x(j) = xt - else - exit outer_up - end if - end do outer_up - if (i == ilx) then - if (x(i) /= piv) then - call psb_errpush(psb_err_internal_error_,& - & r_name='psi_iqsr',a_err='impossible pivot condition') - call psb_error() - endif - i = i + 1 - endif - - n1 = (i-1)-ilx+1 - n2 = iux-(i)+1 - if (n1 > n2) then - if (n1 > ithrs) then - istp = istp + 1 - istack(1,istp) = ilx - istack(2,istp) = i-1 - else - call psi_iisr_up(n1,x(ilx:i-1)) - endif - if (n2 > ithrs) then - istp = istp + 1 - istack(1,istp) = i - istack(2,istp) = iux - else - call psi_iisr_up(n2,x(i:iux)) - endif - else - if (n2 > ithrs) then - istp = istp + 1 - istack(1,istp) = i - istack(2,istp) = iux - else - call psi_iisr_up(n2,x(i:iux)) - endif - if (n1 > ithrs) then - istp = istp + 1 - istack(1,istp) = ilx - istack(2,istp) = i-1 - else - call psi_iisr_up(n1,x(ilx:i-1)) - endif - endif - enddo - else - call psi_iisr_up(n,x) - endif - -end subroutine psi_iqsr_up - -subroutine psi_iqsr_dw(n,x) - use psb_sort_mod, psb_protect_name => psi_iqsr_dw - use psb_error_mod - implicit none - - integer(psb_ipk_), intent(inout) :: x(:) - integer(psb_ipk_), intent(in) :: n - ! .. - ! .. Local Scalars .. - integer(psb_ipk_) :: piv, xt, xk - integer(psb_ipk_) :: i, j, ilx, iux, istp, lpiv - integer(psb_ipk_) :: n1, n2 - - integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=72 - integer(psb_ipk_) :: istack(nparms,maxstack) - - - if (n > ithrs) then - ! - ! Init stack pointer - ! - istp = 1 - istack(1,istp) = 1 - istack(2,istp) = n - - do - if (istp <= 0) exit - ilx = istack(1,istp) - iux = istack(2,istp) - istp = istp - 1 - ! - ! Choose a pivot with median-of-three heuristics, leave it - ! in the LPIV location - ! - i = ilx - j = iux - lpiv = (i+j)/2 - piv = x(lpiv) - if (piv > x(i)) then - xt = x(i) - x(i) = x(lpiv) - x(lpiv) = xt - piv = x(lpiv) - endif - if (piv < x(j)) then - xt = x(j) - x(j) = x(lpiv) - x(lpiv) = xt - piv = x(lpiv) - endif - if (piv > x(i)) then - xt = x(i) - x(i) = x(lpiv) - x(lpiv) = xt - piv = x(lpiv) - endif - ! - ! now piv is correct; place it into first location - - xt = x(i) - x(i) = x(lpiv) - x(lpiv) = xt - - i = ilx - 1 - j = iux + 1 - - outer_dw: do - in_dw1: do - i = i + 1 - xk = x(i) - if (xk <= piv) exit in_dw1 - end do in_dw1 - ! - ! Ensure finite termination for next loop - ! - xt = xk - x(i) = piv - in_dw2:do - j = j - 1 - xk = x(j) - if (xk >= piv) exit in_dw2 - end do in_dw2 - x(i) = xt - - if (j > i) then - xt = x(i) - x(i) = x(j) - x(j) = xt - else - exit outer_dw - end if - end do outer_dw - if (i == ilx) then - if (x(i) /= piv) then - call psb_errpush(psb_err_internal_error_, & - & r_name='psi_iqsr',a_err='impossible pivot condition') - call psb_error() - endif - i = i + 1 - endif - - n1 = (i-1)-ilx+1 - n2 = iux-(i)+1 - if (n1 > n2) then - if (n1 > ithrs) then - istp = istp + 1 - istack(1,istp) = ilx - istack(2,istp) = i-1 - else - call psi_iisr_dw(n1,x(ilx:i-1)) - endif - if (n2 > ithrs) then - istp = istp + 1 - istack(1,istp) = i - istack(2,istp) = iux - else - call psi_iisr_dw(n2,x(i:iux)) - endif - else - if (n2 > ithrs) then - istp = istp + 1 - istack(1,istp) = i - istack(2,istp) = iux - else - call psi_iisr_dw(n2,x(i:iux)) - endif - if (n1 > ithrs) then - istp = istp + 1 - istack(1,istp) = ilx - istack(2,istp) = i-1 - else - call psi_iisr_dw(n1,x(ilx:i-1)) - endif - endif - enddo - else - call psi_iisr_dw(n,x) - endif - -end subroutine psi_iqsr_dw - -subroutine psi_iaqsrx_up(n,x,idx) - use psb_sort_mod, psb_protect_name => psi_iaqsrx_up - use psb_error_mod - implicit none - - integer(psb_ipk_), intent(inout) :: x(:) - integer(psb_ipk_), intent(inout) :: idx(:) - integer(psb_ipk_), intent(in) :: n - ! .. Local Scalars .. - integer(psb_ipk_) :: piv, xk - integer(psb_ipk_) :: xt - integer(psb_ipk_) :: i, j, ilx, iux, istp, lpiv - integer(psb_ipk_) :: n1, n2 - integer(psb_ipk_) :: ixt - - integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=72 - integer(psb_ipk_) :: istack(nparms,maxstack) - - if (n > ithrs) then - ! - ! Init stack pointer - ! - istp = 1 - istack(1,istp) = 1 - istack(2,istp) = n - - do - if (istp <= 0) exit - ilx = istack(1,istp) - iux = istack(2,istp) - istp = istp - 1 - ! - ! Choose a pivot with median-of-three heuristics, leave it - ! in the LPIV location - ! - i = ilx - j = iux - lpiv = (i+j)/2 - piv = abs(x(lpiv)) - if (piv < abs(x(i))) then - xt = x(i) - ixt = idx(i) - x(i) = x(lpiv) - idx(i) = idx(lpiv) - x(lpiv) = xt - idx(lpiv) = ixt - piv = abs(x(lpiv)) - endif - if (piv > abs(x(j))) then - xt = x(j) - ixt = idx(j) - x(j) = x(lpiv) - idx(j) = idx(lpiv) - x(lpiv) = xt - idx(lpiv) = ixt - piv = abs(x(lpiv)) - endif - if (piv < abs(x(i))) then - xt = x(i) - ixt = idx(i) - x(i) = x(lpiv) - idx(i) = idx(lpiv) - x(lpiv) = xt - idx(lpiv) = ixt - piv = abs(x(lpiv)) - endif - ! - ! now piv is correct; place it into first location - xt = x(i) - ixt = idx(i) - x(i) = x(lpiv) - idx(i) = idx(lpiv) - x(lpiv) = xt - idx(lpiv) = ixt - - i = ilx - 1 - j = iux + 1 - - outer_up: do - in_up1: do - i = i + 1 - xk = abs(x(i)) - if (xk >= piv) exit in_up1 - end do in_up1 - ! - ! Ensure finite termination for next loop - ! - xt = x(i) - x(i) = piv - in_up2:do - j = j - 1 - xk = abs(x(j)) - if (xk <= piv) exit in_up2 - end do in_up2 - x(i) = xt - - if (j > i) then - xt = x(i) - ixt = idx(i) - x(i) = x(j) - idx(i) = idx(j) - x(j) = xt - idx(j) = ixt - else - exit outer_up - end if - end do outer_up - if (i == ilx) then - if (x(i) /= piv) then - call psb_errpush(psb_err_internal_error_, & - & r_name='psi_iaqsrx',a_err='impossible pivot condition') - call psb_error() - endif - i = i + 1 - endif - - n1 = (i-1)-ilx+1 - n2 = iux-(i)+1 - if (n1 > n2) then - if (n1 > ithrs) then - istp = istp + 1 - istack(1,istp) = ilx - istack(2,istp) = i-1 - else - call psi_iaisrx_up(n1,x(ilx:i-1),idx(ilx:i-1)) - endif - if (n2 > ithrs) then - istp = istp + 1 - istack(1,istp) = i - istack(2,istp) = iux - else - call psi_iaisrx_up(n2,x(i:iux),idx(i:iux)) - endif - else - if (n2 > ithrs) then - istp = istp + 1 - istack(1,istp) = i - istack(2,istp) = iux - else - call psi_iaisrx_up(n2,x(i:iux),idx(i:iux)) - endif - if (n1 > ithrs) then - istp = istp + 1 - istack(1,istp) = ilx - istack(2,istp) = i-1 - else - call psi_iaisrx_up(n1,x(ilx:i-1),idx(ilx:i-1)) - endif - endif - enddo - else - call psi_iaisrx_up(n,x,idx) - endif - - -end subroutine psi_iaqsrx_up - -subroutine psi_iaqsrx_dw(n,x,idx) - use psb_sort_mod, psb_protect_name => psi_iaqsrx_dw - use psb_error_mod - implicit none - - integer(psb_ipk_), intent(inout) :: x(:) - integer(psb_ipk_), intent(inout) :: idx(:) - integer(psb_ipk_), intent(in) :: n - ! .. Local Scalars .. - integer(psb_ipk_) :: piv, xk - integer(psb_ipk_) :: xt - integer(psb_ipk_) :: i, j, ilx, iux, istp, lpiv - integer(psb_ipk_) :: n1, n2 - integer(psb_ipk_) :: ixt - - integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=72 - integer(psb_ipk_) :: istack(nparms,maxstack) - if (n > ithrs) then - ! - ! Init stack pointer - ! - istp = 1 - istack(1,istp) = 1 - istack(2,istp) = n - - do - if (istp <= 0) exit - ilx = istack(1,istp) - iux = istack(2,istp) - istp = istp - 1 - ! - ! Choose a pivot with median-of-three heuristics, leave it - ! in the LPIV location - ! - i = ilx - j = iux - lpiv = (i+j)/2 - piv = abs(x(lpiv)) - if (piv > abs(x(i))) then - xt = x(i) - ixt = idx(i) - x(i) = x(lpiv) - idx(i) = idx(lpiv) - x(lpiv) = xt - idx(lpiv) = ixt - piv = abs(x(lpiv)) - endif - if (piv < abs(x(j))) then - xt = x(j) - ixt = idx(j) - x(j) = x(lpiv) - idx(j) = idx(lpiv) - x(lpiv) = xt - idx(lpiv) = ixt - piv = abs(x(lpiv)) - endif - if (piv > abs(x(i))) then - xt = x(i) - ixt = idx(i) - x(i) = x(lpiv) - idx(i) = idx(lpiv) - x(lpiv) = xt - idx(lpiv) = ixt - piv = abs(x(lpiv)) - endif - ! - ! now piv is correct; place it into first location - xt = x(i) - ixt = idx(i) - x(i) = x(lpiv) - idx(i) = idx(lpiv) - x(lpiv) = xt - idx(lpiv) = ixt - - i = ilx - 1 - j = iux + 1 - - outer_dw: do - in_dw1: do - i = i + 1 - xk = abs(x(i)) - if (xk <= piv) exit in_dw1 - end do in_dw1 - ! - ! Ensure finite termination for next loop - ! - xt = x(i) - x(i) = piv - in_dw2:do - j = j - 1 - xk = abs(x(j)) - if (xk >= piv) exit in_dw2 - end do in_dw2 - x(i) = xt - - if (j > i) then - xt = x(i) - ixt = idx(i) - x(i) = x(j) - idx(i) = idx(j) - x(j) = xt - idx(j) = ixt - else - exit outer_dw - end if - end do outer_dw - if (i == ilx) then - if (x(i) /= piv) then - call psb_errpush(psb_err_internal_error_,& - & r_name='psi_iaqsrx',a_err='impossible pivot condition') - call psb_error() - endif - i = i + 1 - endif - - n1 = (i-1)-ilx+1 - n2 = iux-(i)+1 - if (n1 > n2) then - if (n1 > ithrs) then - istp = istp + 1 - istack(1,istp) = ilx - istack(2,istp) = i-1 - else - call psi_iaisrx_dw(n1,x(ilx:i-1),idx(ilx:i-1)) - endif - if (n2 > ithrs) then - istp = istp + 1 - istack(1,istp) = i - istack(2,istp) = iux - else - call psi_iaisrx_dw(n2,x(i:iux),idx(i:iux)) - endif - else - if (n2 > ithrs) then - istp = istp + 1 - istack(1,istp) = i - istack(2,istp) = iux - else - call psi_iaisrx_dw(n2,x(i:iux),idx(i:iux)) - endif - if (n1 > ithrs) then - istp = istp + 1 - istack(1,istp) = ilx - istack(2,istp) = i-1 - else - call psi_iaisrx_dw(n1,x(ilx:i-1),idx(ilx:i-1)) - endif - endif - enddo - else - call psi_iaisrx_dw(n,x,idx) - endif - -end subroutine psi_iaqsrx_dw - -subroutine psi_iaqsr_up(n,x) - use psb_sort_mod, psb_protect_name => psi_iaqsr_up - use psb_error_mod - implicit none - - integer(psb_ipk_), intent(inout) :: x(:) - integer(psb_ipk_), intent(in) :: n - ! .. Local Scalars .. - integer(psb_ipk_) :: piv, xk - integer(psb_ipk_) :: xt - integer(psb_ipk_) :: i, j, ilx, iux, istp, lpiv - integer(psb_ipk_) :: n1, n2 - integer(psb_ipk_) :: ixt - - integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=72 - integer(psb_ipk_) :: istack(nparms,maxstack) - - if (n > ithrs) then - ! - ! Init stack pointer - ! - istp = 1 - istack(1,istp) = 1 - istack(2,istp) = n - - do - if (istp <= 0) exit - ilx = istack(1,istp) - iux = istack(2,istp) - istp = istp - 1 - ! - ! Choose a pivot with median-of-three heuristics, leave it - ! in the LPIV location - ! - i = ilx - j = iux - lpiv = (i+j)/2 - piv = abs(x(lpiv)) - if (piv < abs(x(i))) then - xt = x(i) - x(i) = x(lpiv) - x(lpiv) = xt - piv = abs(x(lpiv)) - endif - if (piv > abs(x(j))) then - xt = x(j) - x(j) = x(lpiv) - x(lpiv) = xt - piv = abs(x(lpiv)) - endif - if (piv < abs(x(i))) then - xt = x(i) - x(i) = x(lpiv) - x(lpiv) = xt - piv = abs(x(lpiv)) - endif - ! - ! now piv is correct; place it into first location - - xt = x(i) - x(i) = x(lpiv) - x(lpiv) = xt - - i = ilx - 1 - j = iux + 1 - - outer_up: do - in_up1: do - i = i + 1 - xk = abs(x(i)) - if (xk >= piv) exit in_up1 - end do in_up1 - ! - ! Ensure finite termination for next loop - ! - xt = x(i) - x(i) = piv - in_up2:do - j = j - 1 - xk = abs(x(j)) - if (xk <= piv) exit in_up2 - end do in_up2 - x(i) = xt - - if (j > i) then - xt = x(i) - x(i) = x(j) - x(j) = xt - else - exit outer_up - end if - end do outer_up - if (i == ilx) then - if (x(i) /= piv) then - call psb_errpush(psb_err_internal_error_, & - & r_name='psi_iqasr',a_err='impossible pivot condition') - call psb_error() - endif - i = i + 1 - endif - - n1 = (i-1)-ilx+1 - n2 = iux-(i)+1 - if (n1 > n2) then - if (n1 > ithrs) then - istp = istp + 1 - istack(1,istp) = ilx - istack(2,istp) = i-1 - else - call psi_iaisr_up(n1,x(ilx:i-1)) - endif - if (n2 > ithrs) then - istp = istp + 1 - istack(1,istp) = i - istack(2,istp) = iux - else - call psi_iaisr_up(n2,x(i:iux)) - endif - else - if (n2 > ithrs) then - istp = istp + 1 - istack(1,istp) = i - istack(2,istp) = iux - else - call psi_iaisr_up(n2,x(i:iux)) - endif - if (n1 > ithrs) then - istp = istp + 1 - istack(1,istp) = ilx - istack(2,istp) = i-1 - else - call psi_iaisr_up(n1,x(ilx:i-1)) - endif - endif - enddo - else - call psi_iaisr_up(n,x) - endif - -end subroutine psi_iaqsr_up - -subroutine psi_iaqsr_dw(n,x) - use psb_sort_mod, psb_protect_name => psi_iaqsr_dw - use psb_error_mod - implicit none - - integer(psb_ipk_), intent(inout) :: x(:) - integer(psb_ipk_), intent(in) :: n - ! .. Local Scalars .. - integer(psb_ipk_) :: piv, xk - integer(psb_ipk_) :: xt - integer(psb_ipk_) :: i, j, ilx, iux, istp, lpiv - integer(psb_ipk_) :: n1, n2 - integer(psb_ipk_) :: ixt - - integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=72 - integer(psb_ipk_) :: istack(nparms,maxstack) - - if (n > ithrs) then - ! - ! Init stack pointer - ! - istp = 1 - istack(1,istp) = 1 - istack(2,istp) = n - - do - if (istp <= 0) exit - ilx = istack(1,istp) - iux = istack(2,istp) - istp = istp - 1 - ! - ! Choose a pivot with median-of-three heuristics, leave it - ! in the LPIV location - ! - i = ilx - j = iux - lpiv = (i+j)/2 - piv = abs(x(lpiv)) - if (piv > abs(x(i))) then - xt = x(i) - x(i) = x(lpiv) - x(lpiv) = xt - piv = abs(x(lpiv)) - endif - if (piv < abs(x(j))) then - xt = x(j) - x(j) = x(lpiv) - x(lpiv) = xt - piv = abs(x(lpiv)) - endif - if (piv > abs(x(i))) then - xt = x(i) - x(i) = x(lpiv) - x(lpiv) = xt - piv = abs(x(lpiv)) - endif - ! - ! now piv is correct; place it into first location - - xt = x(i) - x(i) = x(lpiv) - x(lpiv) = xt - - i = ilx - 1 - j = iux + 1 - - outer_dw: do - in_dw1: do - i = i + 1 - xk = abs(x(i)) - if (xk <= piv) exit in_dw1 - end do in_dw1 - ! - ! Ensure finite termination for next loop - ! - xt = x(i) - x(i) = piv - in_dw2:do - j = j - 1 - xk = abs(x(j)) - if (xk >= piv) exit in_dw2 - end do in_dw2 - x(i) = xt - - if (j > i) then - xt = x(i) - x(i) = x(j) - x(j) = xt - else - exit outer_dw - end if - end do outer_dw - if (i == ilx) then - if (x(i) /= piv) then - call psb_errpush(psb_err_internal_error_,& - & r_name='psi_iqasr',a_err='impossible pivot condition') - call psb_error() - endif - i = i + 1 - endif - - n1 = (i-1)-ilx+1 - n2 = iux-(i)+1 - if (n1 > n2) then - if (n1 > ithrs) then - istp = istp + 1 - istack(1,istp) = ilx - istack(2,istp) = i-1 - else - call psi_iaisr_dw(n1,x(ilx:i-1)) - endif - if (n2 > ithrs) then - istp = istp + 1 - istack(1,istp) = i - istack(2,istp) = iux - else - call psi_iaisr_dw(n2,x(i:iux)) - endif - else - if (n2 > ithrs) then - istp = istp + 1 - istack(1,istp) = i - istack(2,istp) = iux - else - call psi_iaisr_dw(n2,x(i:iux)) - endif - if (n1 > ithrs) then - istp = istp + 1 - istack(1,istp) = ilx - istack(2,istp) = i-1 - else - call psi_iaisr_dw(n1,x(ilx:i-1)) - endif - endif - enddo - else - call psi_iaisr_dw(n,x) - endif - -end subroutine psi_iaqsr_dw - - diff --git a/base/serial/sort/psb_l_hsort_impl.f90 b/base/serial/sort/psb_l_hsort_impl.f90 deleted file mode 100644 index f8337eb4..00000000 --- a/base/serial/sort/psb_l_hsort_impl.f90 +++ /dev/null @@ -1,678 +0,0 @@ -! -! 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_lhsort(x,ix,dir,flag) - use psb_sort_mod, psb_protect_name => psb_lhsort - use psb_error_mod - implicit none - integer(psb_lpk_), intent(inout) :: x(:) - integer(psb_ipk_), optional, intent(in) :: dir, flag - integer(psb_lpk_), optional, intent(inout) :: ix(:) - - integer(psb_ipk_) :: dir_, flag_, n, i, l, err_act,info - integer(psb_lpk_) :: key - integer(psb_lpk_) :: 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_l_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_lhsort - - - -! -! 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_l_insert_heap(key,last,heap,dir,info) - use psb_sort_mod, psb_protect_name => psi_l_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_lpk_), intent(in) :: key - integer(psb_ipk_), intent(in) :: dir - integer(psb_lpk_), intent(inout) :: heap(:) - integer(psb_lpk_), intent(inout) :: last - integer(psb_ipk_), intent(out) :: info - integer(psb_ipk_) :: i, i2 - integer(psb_lpk_) :: 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_l_insert_heap - - -subroutine psi_l_heap_get_first(key,last,heap,dir,info) - use psb_sort_mod, psb_protect_name => psi_l_heap_get_first - implicit none - - integer(psb_lpk_), intent(inout) :: key - integer(psb_lpk_), intent(inout) :: last - integer(psb_ipk_), intent(in) :: dir - integer(psb_lpk_), intent(inout) :: heap(:) - integer(psb_ipk_), intent(out) :: info - - integer(psb_ipk_) :: i, j - integer(psb_lpk_) :: 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_l_heap_get_first - - -subroutine psi_l_idx_insert_heap(key,index,last,heap,idxs,dir,info) - use psb_sort_mod, psb_protect_name => psi_l_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_lpk_), intent(in) :: key - integer(psb_ipk_), intent(in) :: index,dir - integer(psb_lpk_), intent(inout) :: heap(:) - integer(psb_ipk_), intent(inout) :: idxs(:),last - integer(psb_ipk_), intent(out) :: info - integer(psb_ipk_) :: i, i2, itemp - integer(psb_lpk_) :: 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_l_idx_insert_heap - -subroutine psi_l_idx_heap_get_first(key,index,last,heap,idxs,dir,info) - use psb_sort_mod, psb_protect_name => psi_l_idx_heap_get_first - implicit none - - integer(psb_lpk_), 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_lpk_), intent(out) :: key - - integer(psb_ipk_) :: i, j,itemp - integer(psb_lpk_) :: 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_l_idx_heap_get_first - - - - diff --git a/base/serial/sort/psb_l_isort_impl.f90 b/base/serial/sort/psb_l_isort_impl.f90 deleted file mode 100644 index 8d101759..00000000 --- a/base/serial/sort/psb_l_isort_impl.f90 +++ /dev/null @@ -1,341 +0,0 @@ -! -! 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_lisort(x,ix,dir,flag) - use psb_sort_mod, psb_protect_name => psb_lisort - use psb_error_mod - implicit none - integer(psb_lpk_), intent(inout) :: x(:) - integer(psb_ipk_), optional, intent(in) :: dir, flag - integer(psb_lpk_), optional, intent(inout) :: ix(:) - - integer(psb_ipk_) :: dir_, flag_, err_act - integer(psb_lpk_) :: n, i - - integer(psb_ipk_) :: ierr(5) - character(len=20) :: name - - name='psb_lisort' - 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_lisrx_up(n,x,ix) - case (psb_sort_down_) - call psi_lisrx_dw(n,x,ix) - case (psb_asort_up_) - call psi_laisrx_up(n,x,ix) - case (psb_asort_down_) - call psi_laisrx_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_lisr_up(n,x) - case (psb_sort_down_) - call psi_lisr_dw(n,x) - case (psb_asort_up_) - call psi_laisr_up(n,x) - case (psb_asort_down_) - call psi_laisr_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_lisort - -subroutine psi_lisrx_up(n,x,idx) - use psb_sort_mod, psb_protect_name => psi_lisrx_up - use psb_error_mod - implicit none - integer(psb_lpk_), intent(inout) :: x(:) - integer(psb_lpk_), intent(inout) :: idx(:) - integer(psb_lpk_), intent(in) :: n - integer(psb_lpk_) :: i,j,ix - integer(psb_lpk_) :: 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_lisrx_up - -subroutine psi_lisrx_dw(n,x,idx) - use psb_sort_mod, psb_protect_name => psi_lisrx_dw - use psb_error_mod - implicit none - integer(psb_lpk_), intent(inout) :: x(:) - integer(psb_lpk_), intent(inout) :: idx(:) - integer(psb_lpk_), intent(in) :: n - integer(psb_lpk_) :: i,j,ix - integer(psb_lpk_) :: 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_lisrx_dw - - -subroutine psi_lisr_up(n,x) - use psb_sort_mod, psb_protect_name => psi_lisr_up - use psb_error_mod - implicit none - integer(psb_lpk_), intent(inout) :: x(:) - integer(psb_lpk_), intent(in) :: n - integer(psb_lpk_) :: i,j - integer(psb_lpk_) :: 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_lisr_up - -subroutine psi_lisr_dw(n,x) - use psb_sort_mod, psb_protect_name => psi_lisr_dw - use psb_error_mod - implicit none - integer(psb_lpk_), intent(inout) :: x(:) - integer(psb_lpk_), intent(in) :: n - integer(psb_lpk_) :: i,j - integer(psb_lpk_) :: 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_lisr_dw - -subroutine psi_laisrx_up(n,x,idx) - use psb_sort_mod, psb_protect_name => psi_laisrx_up - use psb_error_mod - implicit none - integer(psb_lpk_), intent(inout) :: x(:) - integer(psb_lpk_), intent(inout) :: idx(:) - integer(psb_lpk_), intent(in) :: n - integer(psb_lpk_) :: i,j,ix - integer(psb_lpk_) :: 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_laisrx_up - -subroutine psi_laisrx_dw(n,x,idx) - use psb_sort_mod, psb_protect_name => psi_laisrx_dw - use psb_error_mod - implicit none - integer(psb_lpk_), intent(inout) :: x(:) - integer(psb_lpk_), intent(inout) :: idx(:) - integer(psb_lpk_), intent(in) :: n - integer(psb_lpk_) :: i,j,ix - integer(psb_lpk_) :: 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_laisrx_dw - -subroutine psi_laisr_up(n,x) - use psb_sort_mod, psb_protect_name => psi_laisr_up - use psb_error_mod - implicit none - integer(psb_lpk_), intent(inout) :: x(:) - integer(psb_lpk_), intent(in) :: n - integer(psb_lpk_) :: i,j - integer(psb_lpk_) :: 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_laisr_up - -subroutine psi_laisr_dw(n,x) - use psb_sort_mod, psb_protect_name => psi_laisr_dw - use psb_error_mod - implicit none - integer(psb_lpk_), intent(inout) :: x(:) - integer(psb_lpk_), intent(in) :: n - integer(psb_lpk_) :: i,j - integer(psb_lpk_) :: 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_laisr_dw - diff --git a/base/serial/sort/psb_l_msort_impl.f90 b/base/serial/sort/psb_l_msort_impl.f90 deleted file mode 100644 index 2508b332..00000000 --- a/base/serial/sort/psb_l_msort_impl.f90 +++ /dev/null @@ -1,713 +0,0 @@ -! -! 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_lisaperm(n,eip) - use psb_sort_mod, psb_protect_name => psb_lisaperm - implicit none - - integer(psb_lpk_), intent(in) :: n - integer(psb_lpk_), intent(in) :: eip(n) - integer(psb_lpk_), allocatable :: ip(:) - integer(psb_lpk_) :: i,j,m, info - - - psb_lisaperm = .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_lisaperm = .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_lisaperm = .false. - goto 9999 - endif - end if - enddo -9999 continue - - return - end function psb_lisaperm - - - subroutine psb_lmsort_u(x,nout,dir) - use psb_sort_mod, psb_protect_name => psb_lmsort_u - use psb_error_mod - implicit none - integer(psb_lpk_), intent(inout) :: x(:) - integer(psb_lpk_), intent(out) :: nout - integer(psb_ipk_), optional, intent(in) :: dir - - integer(psb_lpk_) :: 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_lmsort_u - - - function psb_lbsrch(key,n,v) result(ipos) - use psb_sort_mod, psb_protect_name => psb_lbsrch - implicit none - integer(psb_ipk_) :: ipos, n - integer(psb_lpk_) :: key - integer(psb_lpk_) :: 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_lbsrch - - function psb_lssrch(key,n,v) result(ipos) - use psb_sort_mod, psb_protect_name => psb_lssrch - implicit none - integer(psb_ipk_) :: ipos, n - integer(psb_lpk_) :: key - integer(psb_lpk_) :: 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_lssrch - - subroutine psb_lmsort(x,ix,dir,flag) - use psb_sort_mod, psb_protect_name => psb_lmsort - use psb_error_mod - use psb_ip_reord_mod - implicit none - integer(psb_lpk_), intent(inout) :: x(:) - integer(psb_ipk_), optional, intent(in) :: dir, flag - integer(psb_lpk_), optional, intent(inout) :: ix(:) - - integer(psb_ipk_) :: dir_, flag_, n, err_act - - integer(psb_lpk_), allocatable :: iaux(:) - integer(psb_ipk_) :: iret, info, i - integer(psb_ipk_) :: ierr(5) - character(len=20) :: name - - name='psb_lmsort' - 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_l_msort') - goto 9999 - endif - - select case(dir_) - case (psb_sort_up_) - call psi_l_msort_up(n,x,iaux,iret) - case (psb_sort_down_) - call psi_l_msort_dw(n,x,iaux,iret) - case (psb_asort_up_) - call psi_l_amsort_up(n,x,iaux,iret) - case (psb_asort_down_) - call psi_l_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_lmsort - - subroutine psi_l_msort_up(n,k,l,iret) - use psb_const_mod - implicit none - integer(psb_ipk_) :: n, iret - integer(psb_lpk_) :: k(n) - integer(psb_lpk_) :: l(0:n+1) - ! - integer(psb_lpk_) :: 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_l_msort_up - - subroutine psi_l_msort_dw(n,k,l,iret) - use psb_const_mod - implicit none - integer(psb_ipk_) :: n, iret - integer(psb_lpk_) :: k(n) - integer(psb_lpk_) :: l(0:n+1) - ! - integer(psb_lpk_) :: 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_l_msort_dw - - subroutine psi_l_amsort_up(n,k,l,iret) - use psb_const_mod - implicit none - integer(psb_ipk_) :: n, iret - integer(psb_lpk_) :: k(n) - integer(psb_lpk_) :: l(0:n+1) - ! - integer(psb_lpk_) :: 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_l_amsort_up - - subroutine psi_l_amsort_dw(n,k,l,iret) - use psb_const_mod - implicit none - integer(psb_ipk_) :: n, iret - integer(psb_lpk_) :: k(n) - integer(psb_lpk_) :: l(0:n+1) - ! - integer(psb_lpk_) :: 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_l_amsort_dw - - - - - - - - diff --git a/base/serial/sort/psb_l_qsort_impl.f90 b/base/serial/sort/psb_l_qsort_impl.f90 deleted file mode 100644 index f02428e2..00000000 --- a/base/serial/sort/psb_l_qsort_impl.f90 +++ /dev/null @@ -1,1318 +0,0 @@ -! -! 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 quicksort routines -! References: -! D. Knuth -! The Art of Computer Programming, vol. 3 -! Addison-Wesley -! -! Aho, Hopcroft, Ullman -! Data Structures and Algorithms -! Addison-Wesley -! -subroutine psb_lqsort(x,ix,dir,flag) - use psb_sort_mod, psb_protect_name => psb_lqsort - use psb_error_mod - implicit none - integer(psb_lpk_), intent(inout) :: x(:) - integer(psb_ipk_), optional, intent(in) :: dir, flag - integer(psb_lpk_), optional, intent(inout) :: ix(:) - - integer(psb_ipk_) :: dir_, flag_, err_act, i - integer(psb_lpk_) :: n - integer(psb_ipk_) :: ierr(5) - character(len=20) :: name - - name='psb_lqsort' - 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_lqsrx_up(n,x,ix) - case (psb_sort_down_) - call psi_lqsrx_dw(n,x,ix) - case (psb_asort_up_) - call psi_laqsrx_up(n,x,ix) - case (psb_asort_down_) - call psi_laqsrx_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_lqsr_up(n,x) - case (psb_sort_down_) - call psi_lqsr_dw(n,x) - case (psb_asort_up_) - call psi_laqsr_up(n,x) - case (psb_asort_down_) - call psi_laqsr_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_lqsort - -subroutine psi_lqsrx_up(n,x,idx) - use psb_sort_mod, psb_protect_name => psi_lqsrx_up - use psb_error_mod - implicit none - - integer(psb_lpk_), intent(inout) :: x(:) - integer(psb_lpk_), intent(inout) :: idx(:) - integer(psb_lpk_), intent(in) :: n - ! .. Local Scalars .. - integer(psb_lpk_) :: piv, xk, xt - integer(psb_ipk_) :: i, j, ilx, iux, istp, lpiv - integer(psb_lpk_) :: n1, n2 - integer(psb_lpk_) :: ixt - integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=72 - integer(psb_ipk_) :: istack(nparms,maxstack) - - if (n > ithrs) then - ! - ! Init stack pointer - ! - istp = 1 - istack(1,istp) = 1 - istack(2,istp) = n - - do - if (istp <= 0) exit - ilx = istack(1,istp) - iux = istack(2,istp) - istp = istp - 1 - ! - ! Choose a pivot with median-of-three heuristics, leave it - ! in the LPIV location - ! - i = ilx - j = iux - lpiv = (i+j)/2 - piv = x(lpiv) - if (piv < x(i)) then - xt = x(i) - ixt = idx(i) - x(i) = x(lpiv) - idx(i) = idx(lpiv) - x(lpiv) = xt - idx(lpiv) = ixt - piv = x(lpiv) - endif - if (piv > x(j)) then - xt = x(j) - ixt = idx(j) - x(j) = x(lpiv) - idx(j) = idx(lpiv) - x(lpiv) = xt - idx(lpiv) = ixt - piv = x(lpiv) - endif - if (piv < x(i)) then - xt = x(i) - ixt = idx(i) - x(i) = x(lpiv) - idx(i) = idx(lpiv) - x(lpiv) = xt - idx(lpiv) = ixt - piv = x(lpiv) - endif - ! - ! now piv is correct; place it into first location - xt = x(i) - ixt = idx(i) - x(i) = x(lpiv) - idx(i) = idx(lpiv) - x(lpiv) = xt - idx(lpiv) = ixt - piv = x(lpiv) - - i = ilx - 1 - j = iux + 1 - - outer_up: do - in_up1: do - i = i + 1 - xk = x(i) - if (xk >= piv) exit in_up1 - end do in_up1 - ! - ! Ensure finite termination for next loop - ! - xt = xk - x(i) = piv - in_up2:do - j = j - 1 - xk = x(j) - if (xk <= piv) exit in_up2 - end do in_up2 - x(i) = xt - - if (j > i) then - xt = x(i) - ixt = idx(i) - x(i) = x(j) - idx(i) = idx(j) - x(j) = xt - idx(j) = ixt - else - exit outer_up - end if - end do outer_up - if (i == ilx) then - if (x(i) /= piv) then - call psb_errpush(psb_err_internal_error_,& - & r_name='psi_lqsrx',a_err='impossible pivot condition') - call psb_error() - endif - i = i + 1 - endif - - n1 = (i-1)-ilx+1 - n2 = iux-(i)+1 - if (n1 > n2) then - if (n1 > ithrs) then - istp = istp + 1 - istack(1,istp) = ilx - istack(2,istp) = i-1 - else - call psi_lisrx_up(n1,x(ilx:i-1),idx(ilx:i-1)) - endif - if (n2 > ithrs) then - istp = istp + 1 - istack(1,istp) = i - istack(2,istp) = iux - else - call psi_lisrx_up(n2,x(i:iux),idx(i:iux)) - endif - else - if (n2 > ithrs) then - istp = istp + 1 - istack(1,istp) = i - istack(2,istp) = iux - else - call psi_lisrx_up(n2,x(i:iux),idx(i:iux)) - endif - if (n1 > ithrs) then - istp = istp + 1 - istack(1,istp) = ilx - istack(2,istp) = i-1 - else - call psi_lisrx_up(n1,x(ilx:i-1),idx(ilx:i-1)) - endif - endif - enddo - else - call psi_lisrx_up(n,x,idx) - endif -end subroutine psi_lqsrx_up - -subroutine psi_lqsrx_dw(n,x,idx) - use psb_sort_mod, psb_protect_name => psi_lqsrx_dw - use psb_error_mod - implicit none - - integer(psb_lpk_), intent(inout) :: x(:) - integer(psb_lpk_), intent(inout) :: idx(:) - integer(psb_lpk_), intent(in) :: n - ! .. Local Scalars .. - integer(psb_lpk_) :: piv, xk, xt - integer(psb_ipk_) :: i, j, ilx, iux, istp, lpiv - integer(psb_lpk_) :: n1, n2 - integer(psb_lpk_) :: ixt - - integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=72 - integer(psb_ipk_) :: istack(nparms,maxstack) - - if (n > ithrs) then - ! - ! Init stack pointer - ! - istp = 1 - istack(1,istp) = 1 - istack(2,istp) = n - - do - if (istp <= 0) exit - ilx = istack(1,istp) - iux = istack(2,istp) - istp = istp - 1 - ! - ! Choose a pivot with median-of-three heuristics, leave it - ! in the LPIV location - ! - i = ilx - j = iux - lpiv = (i+j)/2 - piv = x(lpiv) - if (piv > x(i)) then - xt = x(i) - ixt = idx(i) - x(i) = x(lpiv) - idx(i) = idx(lpiv) - x(lpiv) = xt - idx(lpiv) = ixt - piv = x(lpiv) - endif - if (piv < x(j)) then - xt = x(j) - ixt = idx(j) - x(j) = x(lpiv) - idx(j) = idx(lpiv) - x(lpiv) = xt - idx(lpiv) = ixt - piv = x(lpiv) - endif - if (piv > x(i)) then - xt = x(i) - ixt = idx(i) - x(i) = x(lpiv) - idx(i) = idx(lpiv) - x(lpiv) = xt - idx(lpiv) = ixt - piv = x(lpiv) - endif - ! - ! now piv is correct; place it into first location - xt = x(i) - ixt = idx(i) - x(i) = x(lpiv) - idx(i) = idx(lpiv) - x(lpiv) = xt - idx(lpiv) = ixt - piv = x(lpiv) - - i = ilx - 1 - j = iux + 1 - - outer_dw: do - in_dw1: do - i = i + 1 - xk = x(i) - if (xk <= piv) exit in_dw1 - end do in_dw1 - ! - ! Ensure finite termination for next loop - ! - xt = xk - x(i) = piv - in_dw2:do - j = j - 1 - xk = x(j) - if (xk >= piv) exit in_dw2 - end do in_dw2 - x(i) = xt - - if (j > i) then - xt = x(i) - ixt = idx(i) - x(i) = x(j) - idx(i) = idx(j) - x(j) = xt - idx(j) = ixt - else - exit outer_dw - end if - end do outer_dw - if (i == ilx) then - if (x(i) /= piv) then - call psb_errpush(psb_err_internal_error_,& - & r_name='psi_lqsrx',a_err='impossible pivot condition') - call psb_error() - endif - i = i + 1 - endif - - n1 = (i-1)-ilx+1 - n2 = iux-(i)+1 - if (n1 > n2) then - if (n1 > ithrs) then - istp = istp + 1 - istack(1,istp) = ilx - istack(2,istp) = i-1 - else - call psi_lisrx_dw(n1,x(ilx:i-1),idx(ilx:i-1)) - endif - if (n2 > ithrs) then - istp = istp + 1 - istack(1,istp) = i - istack(2,istp) = iux - else - call psi_lisrx_dw(n2,x(i:iux),idx(i:iux)) - endif - else - if (n2 > ithrs) then - istp = istp + 1 - istack(1,istp) = i - istack(2,istp) = iux - else - call psi_lisrx_dw(n2,x(i:iux),idx(i:iux)) - endif - if (n1 > ithrs) then - istp = istp + 1 - istack(1,istp) = ilx - istack(2,istp) = i-1 - else - call psi_lisrx_dw(n1,x(ilx:i-1),idx(ilx:i-1)) - endif - endif - enddo - else - call psi_lisrx_dw(n,x,idx) - endif - -end subroutine psi_lqsrx_dw - -subroutine psi_lqsr_up(n,x) - use psb_sort_mod, psb_protect_name => psi_lqsr_up - use psb_error_mod - implicit none - - integer(psb_lpk_), intent(inout) :: x(:) - integer(psb_lpk_), intent(in) :: n - ! .. - ! .. Local Scalars .. - integer(psb_lpk_) :: piv, xt, xk - integer(psb_ipk_) :: i, j, ilx, iux, istp, lpiv - integer(psb_lpk_) :: n1, n2 - - integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=72 - integer(psb_ipk_) :: istack(nparms,maxstack) - - - if (n > ithrs) then - ! - ! Init stack pointer - ! - istp = 1 - istack(1,istp) = 1 - istack(2,istp) = n - - do - if (istp <= 0) exit - ilx = istack(1,istp) - iux = istack(2,istp) - istp = istp - 1 - ! - ! Choose a pivot with median-of-three heuristics, leave it - ! in the LPIV location - ! - i = ilx - j = iux - lpiv = (i+j)/2 - piv = x(lpiv) - if (piv < x(i)) then - xt = x(i) - x(i) = x(lpiv) - x(lpiv) = xt - piv = x(lpiv) - endif - if (piv > x(j)) then - xt = x(j) - x(j) = x(lpiv) - x(lpiv) = xt - piv = x(lpiv) - endif - if (piv < x(i)) then - xt = x(i) - x(i) = x(lpiv) - x(lpiv) = xt - piv = x(lpiv) - endif - ! - ! now piv is correct; place it into first location - - xt = x(i) - x(i) = x(lpiv) - x(lpiv) = xt - - i = ilx - 1 - j = iux + 1 - - outer_up: do - in_up1: do - i = i + 1 - xk = x(i) - if (xk >= piv) exit in_up1 - end do in_up1 - ! - ! Ensure finite termination for next loop - ! - xt = xk - x(i) = piv - in_up2:do - j = j - 1 - xk = x(j) - if (xk <= piv) exit in_up2 - end do in_up2 - x(i) = xt - - if (j > i) then - xt = x(i) - x(i) = x(j) - x(j) = xt - else - exit outer_up - end if - end do outer_up - if (i == ilx) then - if (x(i) /= piv) then - call psb_errpush(psb_err_internal_error_,& - & r_name='psi_lqsr',a_err='impossible pivot condition') - call psb_error() - endif - i = i + 1 - endif - - n1 = (i-1)-ilx+1 - n2 = iux-(i)+1 - if (n1 > n2) then - if (n1 > ithrs) then - istp = istp + 1 - istack(1,istp) = ilx - istack(2,istp) = i-1 - else - call psi_lisr_up(n1,x(ilx:i-1)) - endif - if (n2 > ithrs) then - istp = istp + 1 - istack(1,istp) = i - istack(2,istp) = iux - else - call psi_lisr_up(n2,x(i:iux)) - endif - else - if (n2 > ithrs) then - istp = istp + 1 - istack(1,istp) = i - istack(2,istp) = iux - else - call psi_lisr_up(n2,x(i:iux)) - endif - if (n1 > ithrs) then - istp = istp + 1 - istack(1,istp) = ilx - istack(2,istp) = i-1 - else - call psi_lisr_up(n1,x(ilx:i-1)) - endif - endif - enddo - else - call psi_lisr_up(n,x) - endif - -end subroutine psi_lqsr_up - -subroutine psi_lqsr_dw(n,x) - use psb_sort_mod, psb_protect_name => psi_lqsr_dw - use psb_error_mod - implicit none - - integer(psb_lpk_), intent(inout) :: x(:) - integer(psb_lpk_), intent(in) :: n - ! .. - ! .. Local Scalars .. - integer(psb_lpk_) :: piv, xt, xk - integer(psb_ipk_) :: i, j, ilx, iux, istp, lpiv - integer(psb_lpk_) :: n1, n2 - - integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=72 - integer(psb_ipk_) :: istack(nparms,maxstack) - - - if (n > ithrs) then - ! - ! Init stack pointer - ! - istp = 1 - istack(1,istp) = 1 - istack(2,istp) = n - - do - if (istp <= 0) exit - ilx = istack(1,istp) - iux = istack(2,istp) - istp = istp - 1 - ! - ! Choose a pivot with median-of-three heuristics, leave it - ! in the LPIV location - ! - i = ilx - j = iux - lpiv = (i+j)/2 - piv = x(lpiv) - if (piv > x(i)) then - xt = x(i) - x(i) = x(lpiv) - x(lpiv) = xt - piv = x(lpiv) - endif - if (piv < x(j)) then - xt = x(j) - x(j) = x(lpiv) - x(lpiv) = xt - piv = x(lpiv) - endif - if (piv > x(i)) then - xt = x(i) - x(i) = x(lpiv) - x(lpiv) = xt - piv = x(lpiv) - endif - ! - ! now piv is correct; place it into first location - - xt = x(i) - x(i) = x(lpiv) - x(lpiv) = xt - - i = ilx - 1 - j = iux + 1 - - outer_dw: do - in_dw1: do - i = i + 1 - xk = x(i) - if (xk <= piv) exit in_dw1 - end do in_dw1 - ! - ! Ensure finite termination for next loop - ! - xt = xk - x(i) = piv - in_dw2:do - j = j - 1 - xk = x(j) - if (xk >= piv) exit in_dw2 - end do in_dw2 - x(i) = xt - - if (j > i) then - xt = x(i) - x(i) = x(j) - x(j) = xt - else - exit outer_dw - end if - end do outer_dw - if (i == ilx) then - if (x(i) /= piv) then - call psb_errpush(psb_err_internal_error_, & - & r_name='psi_lqsr',a_err='impossible pivot condition') - call psb_error() - endif - i = i + 1 - endif - - n1 = (i-1)-ilx+1 - n2 = iux-(i)+1 - if (n1 > n2) then - if (n1 > ithrs) then - istp = istp + 1 - istack(1,istp) = ilx - istack(2,istp) = i-1 - else - call psi_lisr_dw(n1,x(ilx:i-1)) - endif - if (n2 > ithrs) then - istp = istp + 1 - istack(1,istp) = i - istack(2,istp) = iux - else - call psi_lisr_dw(n2,x(i:iux)) - endif - else - if (n2 > ithrs) then - istp = istp + 1 - istack(1,istp) = i - istack(2,istp) = iux - else - call psi_lisr_dw(n2,x(i:iux)) - endif - if (n1 > ithrs) then - istp = istp + 1 - istack(1,istp) = ilx - istack(2,istp) = i-1 - else - call psi_lisr_dw(n1,x(ilx:i-1)) - endif - endif - enddo - else - call psi_lisr_dw(n,x) - endif - -end subroutine psi_lqsr_dw - -subroutine psi_laqsrx_up(n,x,idx) - use psb_sort_mod, psb_protect_name => psi_laqsrx_up - use psb_error_mod - implicit none - - integer(psb_lpk_), intent(inout) :: x(:) - integer(psb_lpk_), intent(inout) :: idx(:) - integer(psb_lpk_), intent(in) :: n - ! .. Local Scalars .. - integer(psb_lpk_) :: piv, xk - integer(psb_lpk_) :: xt - integer(psb_ipk_) :: i, j, ilx, iux, istp, lpiv - integer(psb_lpk_) :: n1, n2 - integer(psb_lpk_) :: ixt - - integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=72 - integer(psb_ipk_) :: istack(nparms,maxstack) - - if (n > ithrs) then - ! - ! Init stack pointer - ! - istp = 1 - istack(1,istp) = 1 - istack(2,istp) = n - - do - if (istp <= 0) exit - ilx = istack(1,istp) - iux = istack(2,istp) - istp = istp - 1 - ! - ! Choose a pivot with median-of-three heuristics, leave it - ! in the LPIV location - ! - i = ilx - j = iux - lpiv = (i+j)/2 - piv = abs(x(lpiv)) - if (piv < abs(x(i))) then - xt = x(i) - ixt = idx(i) - x(i) = x(lpiv) - idx(i) = idx(lpiv) - x(lpiv) = xt - idx(lpiv) = ixt - piv = abs(x(lpiv)) - endif - if (piv > abs(x(j))) then - xt = x(j) - ixt = idx(j) - x(j) = x(lpiv) - idx(j) = idx(lpiv) - x(lpiv) = xt - idx(lpiv) = ixt - piv = abs(x(lpiv)) - endif - if (piv < abs(x(i))) then - xt = x(i) - ixt = idx(i) - x(i) = x(lpiv) - idx(i) = idx(lpiv) - x(lpiv) = xt - idx(lpiv) = ixt - piv = abs(x(lpiv)) - endif - ! - ! now piv is correct; place it into first location - xt = x(i) - ixt = idx(i) - x(i) = x(lpiv) - idx(i) = idx(lpiv) - x(lpiv) = xt - idx(lpiv) = ixt - - i = ilx - 1 - j = iux + 1 - - outer_up: do - in_up1: do - i = i + 1 - xk = abs(x(i)) - if (xk >= piv) exit in_up1 - end do in_up1 - ! - ! Ensure finite termination for next loop - ! - xt = x(i) - x(i) = piv - in_up2:do - j = j - 1 - xk = abs(x(j)) - if (xk <= piv) exit in_up2 - end do in_up2 - x(i) = xt - - if (j > i) then - xt = x(i) - ixt = idx(i) - x(i) = x(j) - idx(i) = idx(j) - x(j) = xt - idx(j) = ixt - else - exit outer_up - end if - end do outer_up - if (i == ilx) then - if (x(i) /= piv) then - call psb_errpush(psb_err_internal_error_, & - & r_name='psi_laqsrx',a_err='impossible pivot condition') - call psb_error() - endif - i = i + 1 - endif - - n1 = (i-1)-ilx+1 - n2 = iux-(i)+1 - if (n1 > n2) then - if (n1 > ithrs) then - istp = istp + 1 - istack(1,istp) = ilx - istack(2,istp) = i-1 - else - call psi_laisrx_up(n1,x(ilx:i-1),idx(ilx:i-1)) - endif - if (n2 > ithrs) then - istp = istp + 1 - istack(1,istp) = i - istack(2,istp) = iux - else - call psi_laisrx_up(n2,x(i:iux),idx(i:iux)) - endif - else - if (n2 > ithrs) then - istp = istp + 1 - istack(1,istp) = i - istack(2,istp) = iux - else - call psi_laisrx_up(n2,x(i:iux),idx(i:iux)) - endif - if (n1 > ithrs) then - istp = istp + 1 - istack(1,istp) = ilx - istack(2,istp) = i-1 - else - call psi_laisrx_up(n1,x(ilx:i-1),idx(ilx:i-1)) - endif - endif - enddo - else - call psi_laisrx_up(n,x,idx) - endif - - -end subroutine psi_laqsrx_up - -subroutine psi_laqsrx_dw(n,x,idx) - use psb_sort_mod, psb_protect_name => psi_laqsrx_dw - use psb_error_mod - implicit none - - integer(psb_lpk_), intent(inout) :: x(:) - integer(psb_lpk_), intent(inout) :: idx(:) - integer(psb_lpk_), intent(in) :: n - ! .. Local Scalars .. - integer(psb_lpk_) :: piv, xk - integer(psb_lpk_) :: xt - integer(psb_ipk_) :: i, j, ilx, iux, istp, lpiv - integer(psb_lpk_) :: n1, n2 - integer(psb_lpk_) :: ixt - - integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=72 - integer(psb_ipk_) :: istack(nparms,maxstack) - if (n > ithrs) then - ! - ! Init stack pointer - ! - istp = 1 - istack(1,istp) = 1 - istack(2,istp) = n - - do - if (istp <= 0) exit - ilx = istack(1,istp) - iux = istack(2,istp) - istp = istp - 1 - ! - ! Choose a pivot with median-of-three heuristics, leave it - ! in the LPIV location - ! - i = ilx - j = iux - lpiv = (i+j)/2 - piv = abs(x(lpiv)) - if (piv > abs(x(i))) then - xt = x(i) - ixt = idx(i) - x(i) = x(lpiv) - idx(i) = idx(lpiv) - x(lpiv) = xt - idx(lpiv) = ixt - piv = abs(x(lpiv)) - endif - if (piv < abs(x(j))) then - xt = x(j) - ixt = idx(j) - x(j) = x(lpiv) - idx(j) = idx(lpiv) - x(lpiv) = xt - idx(lpiv) = ixt - piv = abs(x(lpiv)) - endif - if (piv > abs(x(i))) then - xt = x(i) - ixt = idx(i) - x(i) = x(lpiv) - idx(i) = idx(lpiv) - x(lpiv) = xt - idx(lpiv) = ixt - piv = abs(x(lpiv)) - endif - ! - ! now piv is correct; place it into first location - xt = x(i) - ixt = idx(i) - x(i) = x(lpiv) - idx(i) = idx(lpiv) - x(lpiv) = xt - idx(lpiv) = ixt - - i = ilx - 1 - j = iux + 1 - - outer_dw: do - in_dw1: do - i = i + 1 - xk = abs(x(i)) - if (xk <= piv) exit in_dw1 - end do in_dw1 - ! - ! Ensure finite termination for next loop - ! - xt = x(i) - x(i) = piv - in_dw2:do - j = j - 1 - xk = abs(x(j)) - if (xk >= piv) exit in_dw2 - end do in_dw2 - x(i) = xt - - if (j > i) then - xt = x(i) - ixt = idx(i) - x(i) = x(j) - idx(i) = idx(j) - x(j) = xt - idx(j) = ixt - else - exit outer_dw - end if - end do outer_dw - if (i == ilx) then - if (x(i) /= piv) then - call psb_errpush(psb_err_internal_error_,& - & r_name='psi_laqsrx',a_err='impossible pivot condition') - call psb_error() - endif - i = i + 1 - endif - - n1 = (i-1)-ilx+1 - n2 = iux-(i)+1 - if (n1 > n2) then - if (n1 > ithrs) then - istp = istp + 1 - istack(1,istp) = ilx - istack(2,istp) = i-1 - else - call psi_laisrx_dw(n1,x(ilx:i-1),idx(ilx:i-1)) - endif - if (n2 > ithrs) then - istp = istp + 1 - istack(1,istp) = i - istack(2,istp) = iux - else - call psi_laisrx_dw(n2,x(i:iux),idx(i:iux)) - endif - else - if (n2 > ithrs) then - istp = istp + 1 - istack(1,istp) = i - istack(2,istp) = iux - else - call psi_laisrx_dw(n2,x(i:iux),idx(i:iux)) - endif - if (n1 > ithrs) then - istp = istp + 1 - istack(1,istp) = ilx - istack(2,istp) = i-1 - else - call psi_laisrx_dw(n1,x(ilx:i-1),idx(ilx:i-1)) - endif - endif - enddo - else - call psi_laisrx_dw(n,x,idx) - endif - -end subroutine psi_laqsrx_dw - -subroutine psi_laqsr_up(n,x) - use psb_sort_mod, psb_protect_name => psi_laqsr_up - use psb_error_mod - implicit none - - integer(psb_lpk_), intent(inout) :: x(:) - integer(psb_lpk_), intent(in) :: n - ! .. Local Scalars .. - integer(psb_lpk_) :: piv, xk - integer(psb_lpk_) :: xt - integer(psb_ipk_) :: i, j, ilx, iux, istp, lpiv - integer(psb_lpk_) :: n1, n2 - integer(psb_lpk_) :: ixt - - integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=72 - integer(psb_ipk_) :: istack(nparms,maxstack) - - if (n > ithrs) then - ! - ! Init stack pointer - ! - istp = 1 - istack(1,istp) = 1 - istack(2,istp) = n - - do - if (istp <= 0) exit - ilx = istack(1,istp) - iux = istack(2,istp) - istp = istp - 1 - ! - ! Choose a pivot with median-of-three heuristics, leave it - ! in the LPIV location - ! - i = ilx - j = iux - lpiv = (i+j)/2 - piv = abs(x(lpiv)) - if (piv < abs(x(i))) then - xt = x(i) - x(i) = x(lpiv) - x(lpiv) = xt - piv = abs(x(lpiv)) - endif - if (piv > abs(x(j))) then - xt = x(j) - x(j) = x(lpiv) - x(lpiv) = xt - piv = abs(x(lpiv)) - endif - if (piv < abs(x(i))) then - xt = x(i) - x(i) = x(lpiv) - x(lpiv) = xt - piv = abs(x(lpiv)) - endif - ! - ! now piv is correct; place it into first location - - xt = x(i) - x(i) = x(lpiv) - x(lpiv) = xt - - i = ilx - 1 - j = iux + 1 - - outer_up: do - in_up1: do - i = i + 1 - xk = abs(x(i)) - if (xk >= piv) exit in_up1 - end do in_up1 - ! - ! Ensure finite termination for next loop - ! - xt = x(i) - x(i) = piv - in_up2:do - j = j - 1 - xk = abs(x(j)) - if (xk <= piv) exit in_up2 - end do in_up2 - x(i) = xt - - if (j > i) then - xt = x(i) - x(i) = x(j) - x(j) = xt - else - exit outer_up - end if - end do outer_up - if (i == ilx) then - if (x(i) /= piv) then - call psb_errpush(psb_err_internal_error_, & - & r_name='psi_lqasr',a_err='impossible pivot condition') - call psb_error() - endif - i = i + 1 - endif - - n1 = (i-1)-ilx+1 - n2 = iux-(i)+1 - if (n1 > n2) then - if (n1 > ithrs) then - istp = istp + 1 - istack(1,istp) = ilx - istack(2,istp) = i-1 - else - call psi_laisr_up(n1,x(ilx:i-1)) - endif - if (n2 > ithrs) then - istp = istp + 1 - istack(1,istp) = i - istack(2,istp) = iux - else - call psi_laisr_up(n2,x(i:iux)) - endif - else - if (n2 > ithrs) then - istp = istp + 1 - istack(1,istp) = i - istack(2,istp) = iux - else - call psi_laisr_up(n2,x(i:iux)) - endif - if (n1 > ithrs) then - istp = istp + 1 - istack(1,istp) = ilx - istack(2,istp) = i-1 - else - call psi_laisr_up(n1,x(ilx:i-1)) - endif - endif - enddo - else - call psi_laisr_up(n,x) - endif - -end subroutine psi_laqsr_up - -subroutine psi_laqsr_dw(n,x) - use psb_sort_mod, psb_protect_name => psi_laqsr_dw - use psb_error_mod - implicit none - - integer(psb_lpk_), intent(inout) :: x(:) - integer(psb_lpk_), intent(in) :: n - ! .. Local Scalars .. - integer(psb_lpk_) :: piv, xk - integer(psb_lpk_) :: xt - integer(psb_ipk_) :: i, j, ilx, iux, istp, lpiv - integer(psb_lpk_) :: n1, n2 - integer(psb_lpk_) :: ixt - - integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=72 - integer(psb_ipk_) :: istack(nparms,maxstack) - - if (n > ithrs) then - ! - ! Init stack pointer - ! - istp = 1 - istack(1,istp) = 1 - istack(2,istp) = n - - do - if (istp <= 0) exit - ilx = istack(1,istp) - iux = istack(2,istp) - istp = istp - 1 - ! - ! Choose a pivot with median-of-three heuristics, leave it - ! in the LPIV location - ! - i = ilx - j = iux - lpiv = (i+j)/2 - piv = abs(x(lpiv)) - if (piv > abs(x(i))) then - xt = x(i) - x(i) = x(lpiv) - x(lpiv) = xt - piv = abs(x(lpiv)) - endif - if (piv < abs(x(j))) then - xt = x(j) - x(j) = x(lpiv) - x(lpiv) = xt - piv = abs(x(lpiv)) - endif - if (piv > abs(x(i))) then - xt = x(i) - x(i) = x(lpiv) - x(lpiv) = xt - piv = abs(x(lpiv)) - endif - ! - ! now piv is correct; place it into first location - - xt = x(i) - x(i) = x(lpiv) - x(lpiv) = xt - - i = ilx - 1 - j = iux + 1 - - outer_dw: do - in_dw1: do - i = i + 1 - xk = abs(x(i)) - if (xk <= piv) exit in_dw1 - end do in_dw1 - ! - ! Ensure finite termination for next loop - ! - xt = x(i) - x(i) = piv - in_dw2:do - j = j - 1 - xk = abs(x(j)) - if (xk >= piv) exit in_dw2 - end do in_dw2 - x(i) = xt - - if (j > i) then - xt = x(i) - x(i) = x(j) - x(j) = xt - else - exit outer_dw - end if - end do outer_dw - if (i == ilx) then - if (x(i) /= piv) then - call psb_errpush(psb_err_internal_error_,& - & r_name='psi_lqasr',a_err='impossible pivot condition') - call psb_error() - endif - i = i + 1 - endif - - n1 = (i-1)-ilx+1 - n2 = iux-(i)+1 - if (n1 > n2) then - if (n1 > ithrs) then - istp = istp + 1 - istack(1,istp) = ilx - istack(2,istp) = i-1 - else - call psi_laisr_dw(n1,x(ilx:i-1)) - endif - if (n2 > ithrs) then - istp = istp + 1 - istack(1,istp) = i - istack(2,istp) = iux - else - call psi_laisr_dw(n2,x(i:iux)) - endif - else - if (n2 > ithrs) then - istp = istp + 1 - istack(1,istp) = i - istack(2,istp) = iux - else - call psi_laisr_dw(n2,x(i:iux)) - endif - if (n1 > ithrs) then - istp = istp + 1 - istack(1,istp) = ilx - istack(2,istp) = i-1 - else - call psi_laisr_dw(n1,x(ilx:i-1)) - endif - endif - enddo - else - call psi_laisr_dw(n,x) - endif - -end subroutine psi_laqsr_dw - -