From adddf964d3a55b6f0e2033821000a7fc6cd5f319 Mon Sep 17 00:00:00 2001 From: Salvatore Filippone Date: Fri, 17 Jun 2016 15:22:28 +0000 Subject: [PATCH] psblas3: base/serial/sort/psb_c_qsort_impl.f90 base/serial/sort/psb_d_qsort_impl.f90 base/serial/sort/psb_i_qsort_impl.f90 base/serial/sort/psb_s_qsort_impl.f90 base/serial/sort/psb_z_qsort_impl.f90 Revised quicksort thresholds. --- base/serial/sort/psb_c_qsort_impl.f90 | 602 ++++++++++++- base/serial/sort/psb_d_qsort_impl.f90 | 1188 ++++++++++++++++++++++++- base/serial/sort/psb_i_qsort_impl.f90 | 17 +- base/serial/sort/psb_s_qsort_impl.f90 | 1188 ++++++++++++++++++++++++- base/serial/sort/psb_z_qsort_impl.f90 | 602 ++++++++++++- 5 files changed, 3546 insertions(+), 51 deletions(-) diff --git a/base/serial/sort/psb_c_qsort_impl.f90 b/base/serial/sort/psb_c_qsort_impl.f90 index cf5b7b0e..1f81743f 100644 --- a/base/serial/sort/psb_c_qsort_impl.f90 +++ b/base/serial/sort/psb_c_qsort_impl.f90 @@ -137,6 +137,584 @@ subroutine psb_cqsort(x,ix,dir,flag) return end subroutine psb_cqsort + integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=32 + 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_cqsrx',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_cisrx_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_cisrx_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_cisrx_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_cisrx_up(n1,x(ilx:i-1),idx(ilx:i-1)) + endif + endif + enddo + else + call psi_cisrx_up(n,x,idx) + endif +end subroutine psi_cqsrx_up + +subroutine psi_cqsrx_dw(n,x,idx) + use psb_c_sort_mod, psb_protect_name => psi_cqsrx_dw + use psb_error_mod + implicit none + + complex(psb_spk_), intent(inout) :: x(:) + integer(psb_ipk_), intent(inout) :: idx(:) + integer(psb_ipk_), intent(in) :: n + ! .. Local Scalars .. + complex(psb_spk_) :: piv, xk, xt + integer(psb_ipk_) :: i, j, ilx, iux, istp, lpiv + integer(psb_ipk_) :: ixt, n1, n2 + + integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=32 + 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_cqsrx',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_cisrx_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_cisrx_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_cisrx_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_cisrx_dw(n1,x(ilx:i-1),idx(ilx:i-1)) + endif + endif + enddo + else + call psi_cisrx_dw(n,x,idx) + endif + +end subroutine psi_cqsrx_dw + +subroutine psi_cqsr_up(n,x) + use psb_c_sort_mod, psb_protect_name => psi_cqsr_up + use psb_error_mod + implicit none + + complex(psb_spk_), intent(inout) :: x(:) + integer(psb_ipk_), intent(in) :: n + ! .. + ! .. Local Scalars .. + complex(psb_spk_) :: 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=32 + 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_cqsr',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_cisr_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_cisr_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_cisr_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_cisr_up(n1,x(ilx:i-1)) + endif + endif + enddo + else + call psi_cisr_up(n,x) + endif + +end subroutine psi_cqsr_up + +subroutine psi_cqsr_dw(n,x) + use psb_c_sort_mod, psb_protect_name => psi_cqsr_dw + use psb_error_mod + implicit none + + complex(psb_spk_), intent(inout) :: x(:) + integer(psb_ipk_), intent(in) :: n + ! .. + ! .. Local Scalars .. + complex(psb_spk_) :: 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=32 + 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_cqsr',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_cisr_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_cisr_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_cisr_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_cisr_dw(n1,x(ilx:i-1)) + endif + endif + enddo + else + call psi_cisr_dw(n,x) + endif + +end subroutine psi_cqsr_dw +@REALE@ subroutine psi_clqsrx_up(n,x,idx) use psb_c_sort_mod, psb_protect_name => psi_clqsrx_up @@ -152,7 +730,7 @@ subroutine psi_clqsrx_up(n,x,idx) integer(psb_ipk_) :: i, j, ilx, iux, istp, lpiv integer(psb_ipk_) :: ixt, n1, n2 - integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=16 + integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=32 integer(psb_ipk_) :: istack(nparms,maxstack) if (n > ithrs) then @@ -308,7 +886,7 @@ subroutine psi_clqsrx_dw(n,x,idx) integer(psb_ipk_) :: i, j, ilx, iux, istp, lpiv integer(psb_ipk_) :: ixt, n1, n2 - integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=16 + integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=32 integer(psb_ipk_) :: istack(nparms,maxstack) if (n > ithrs) then @@ -463,7 +1041,7 @@ subroutine psi_clqsr_up(n,x) integer(psb_ipk_) :: i, j, ilx, iux, istp, lpiv integer(psb_ipk_) :: n1, n2 - integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=16 + integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=32 integer(psb_ipk_) :: istack(nparms,maxstack) @@ -604,7 +1182,7 @@ subroutine psi_clqsr_dw(n,x) integer(psb_ipk_) :: i, j, ilx, iux, istp, lpiv integer(psb_ipk_) :: n1, n2 - integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=16 + integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=32 integer(psb_ipk_) :: istack(nparms,maxstack) @@ -746,7 +1324,7 @@ subroutine psi_calqsrx_up(n,x,idx) integer(psb_ipk_) :: i, j, ilx, iux, istp, lpiv integer(psb_ipk_) :: ixt, n1, n2 - integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=16 + integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=32 integer(psb_ipk_) :: istack(nparms,maxstack) if (n > ithrs) then @@ -901,7 +1479,7 @@ subroutine psi_calqsrx_dw(n,x,idx) integer(psb_ipk_) :: i, j, ilx, iux, istp, lpiv integer(psb_ipk_) :: ixt, n1, n2 - integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=16 + integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=32 integer(psb_ipk_) :: istack(nparms,maxstack) if (n > ithrs) then @@ -1056,7 +1634,7 @@ subroutine psi_calqsr_up(n,x) integer(psb_ipk_) :: i, j, ilx, iux, istp, lpiv integer(psb_ipk_) :: n1, n2 - integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=16 + integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=32 integer(psb_ipk_) :: istack(nparms,maxstack) @@ -1196,7 +1774,7 @@ subroutine psi_calqsr_dw(n,x) integer(psb_ipk_) :: i, j, ilx, iux, istp, lpiv integer(psb_ipk_) :: n1, n2 - integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=16 + integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=32 integer(psb_ipk_) :: istack(nparms,maxstack) @@ -1337,7 +1915,7 @@ subroutine psi_caqsrx_up(n,x,idx) integer(psb_ipk_) :: i, j, ilx, iux, istp, lpiv integer(psb_ipk_) :: ixt, n1, n2 - integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=16 + integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=32 integer(psb_ipk_) :: istack(nparms,maxstack) if (n > ithrs) then @@ -1493,7 +2071,7 @@ subroutine psi_caqsrx_dw(n,x,idx) integer(psb_ipk_) :: i, j, ilx, iux, istp, lpiv integer(psb_ipk_) :: ixt, n1, n2 - integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=16 + integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=32 integer(psb_ipk_) :: istack(nparms,maxstack) if (n > ithrs) then ! @@ -1646,7 +2224,7 @@ subroutine psi_caqsr_up(n,x) integer(psb_ipk_) :: i, j, ilx, iux, istp, lpiv integer(psb_ipk_) :: ixt, n1, n2 - integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=16 + integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=32 integer(psb_ipk_) :: istack(nparms,maxstack) if (n > ithrs) then @@ -1786,7 +2364,7 @@ subroutine psi_caqsr_dw(n,x) integer(psb_ipk_) :: i, j, ilx, iux, istp, lpiv integer(psb_ipk_) :: ixt, n1, n2 - integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=16 + integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=32 integer(psb_ipk_) :: istack(nparms,maxstack) if (n > ithrs) then diff --git a/base/serial/sort/psb_d_qsort_impl.f90 b/base/serial/sort/psb_d_qsort_impl.f90 index c0eda8ed..f3a40c06 100644 --- a/base/serial/sort/psb_d_qsort_impl.f90 +++ b/base/serial/sort/psb_d_qsort_impl.f90 @@ -141,8 +141,7 @@ subroutine psi_dqsrx_up(n,x,idx) real(psb_dpk_) :: piv, xk, xt integer(psb_ipk_) :: i, j, ilx, iux, istp, lpiv integer(psb_ipk_) :: ixt, n1, n2 - - integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=16 + integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=60 integer(psb_ipk_) :: istack(nparms,maxstack) if (n > ithrs) then @@ -296,7 +295,7 @@ subroutine psi_dqsrx_dw(n,x,idx) integer(psb_ipk_) :: i, j, ilx, iux, istp, lpiv integer(psb_ipk_) :: ixt, n1, n2 - integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=16 + integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=60 integer(psb_ipk_) :: istack(nparms,maxstack) if (n > ithrs) then @@ -451,7 +450,7 @@ subroutine psi_dqsr_up(n,x) integer(psb_ipk_) :: i, j, ilx, iux, istp, lpiv integer(psb_ipk_) :: n1, n2 - integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=16 + integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=60 integer(psb_ipk_) :: istack(nparms,maxstack) @@ -592,7 +591,7 @@ subroutine psi_dqsr_dw(n,x) integer(psb_ipk_) :: i, j, ilx, iux, istp, lpiv integer(psb_ipk_) :: n1, n2 - integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=16 + integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=60 integer(psb_ipk_) :: istack(nparms,maxstack) @@ -720,6 +719,1177 @@ subroutine psi_dqsr_dw(n,x) end subroutine psi_dqsr_dw + 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_dqsrx',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_dlisrx_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_dlisrx_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_dlisrx_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_dlisrx_up(n1,x(ilx:i-1),idx(ilx:i-1)) + endif + endif + enddo + else + call psi_dlisrx_up(n,x,idx) + endif + +end subroutine psi_dlqsrx_up + +subroutine psi_dlqsrx_dw(n,x,idx) + use psb_d_sort_mod, psb_protect_name => psi_dlqsrx_dw + use psb_error_mod + use psi_lcx_mod + implicit none + + real(psb_dpk_), intent(inout) :: x(:) + integer(psb_ipk_), intent(inout) :: idx(:) + integer(psb_ipk_), intent(in) :: n + ! .. Local Scalars .. + real(psb_dpk_) :: piv, xk, xt + integer(psb_ipk_) :: i, j, ilx, iux, istp, lpiv + integer(psb_ipk_) :: ixt, n1, n2 + + integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=60 + 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_dqsrx',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_dlisrx_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_dlisrx_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_dlisrx_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_dlisrx_dw(n1,x(ilx:i-1),idx(ilx:i-1)) + endif + endif + enddo + else + call psi_dlisrx_dw(n,x,idx) + endif +end subroutine psi_dlqsrx_dw + +subroutine psi_dlqsr_up(n,x) + use psb_d_sort_mod, psb_protect_name => psi_dlqsr_up + use psb_error_mod + use psi_lcx_mod + implicit none + + real(psb_dpk_), intent(inout) :: x(:) + integer(psb_ipk_), intent(in) :: n + ! .. + ! .. Local Scalars .. + real(psb_dpk_) :: 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=60 + 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_dqsr',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_dlisr_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_dlisr_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_dlisr_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_dlisr_up(n1,x(ilx:i-1)) + endif + endif + enddo + else + call psi_dlisr_up(n,x) + endif + +end subroutine psi_dlqsr_up + +subroutine psi_dlqsr_dw(n,x) + use psb_d_sort_mod, psb_protect_name => psi_dlqsr_dw + use psb_error_mod + use psi_lcx_mod + implicit none + + real(psb_dpk_), intent(inout) :: x(:) + integer(psb_ipk_), intent(in) :: n + ! .. Local Scalars .. + real(psb_dpk_) :: 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=60 + 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_dqsr',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_dlisr_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_dlisr_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_dlisr_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_dlisr_dw(n1,x(ilx:i-1)) + endif + endif + enddo + else + call psi_dlisr_dw(n,x) + endif + +end subroutine psi_dlqsr_dw + +subroutine psi_dalqsrx_up(n,x,idx) + use psb_d_sort_mod, psb_protect_name => psi_dalqsrx_up + use psb_error_mod + use psi_alcx_mod + implicit none + + real(psb_dpk_), intent(inout) :: x(:) + integer(psb_ipk_), intent(inout) :: idx(:) + integer(psb_ipk_), intent(in) :: n + ! .. Local Scalars .. + real(psb_dpk_) :: piv, xk, xt + integer(psb_ipk_) :: i, j, ilx, iux, istp, lpiv + integer(psb_ipk_) :: ixt, n1, n2 + + integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=60 + 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_dqsrx',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_dalisrx_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_dalisrx_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_dalisrx_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_dalisrx_up(n1,x(ilx:i-1),idx(ilx:i-1)) + endif + endif + enddo + else + call psi_dalisrx_up(n,x,idx) + endif +end subroutine psi_dalqsrx_up + +subroutine psi_dalqsrx_dw(n,x,idx) + use psb_d_sort_mod, psb_protect_name => psi_dalqsrx_dw + use psb_error_mod + use psi_alcx_mod + implicit none + + real(psb_dpk_), intent(inout) :: x(:) + integer(psb_ipk_), intent(inout) :: idx(:) + integer(psb_ipk_), intent(in) :: n + ! .. Local Scalars .. + real(psb_dpk_) :: piv, xk, xt + integer(psb_ipk_) :: i, j, ilx, iux, istp, lpiv + integer(psb_ipk_) :: ixt, n1, n2 + + integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=60 + 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_dqsrx',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_dalisrx_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_dalisrx_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_dalisrx_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_dalisrx_dw(n1,x(ilx:i-1),idx(ilx:i-1)) + endif + endif + enddo + else + call psi_dalisrx_dw(n,x,idx) + endif +end subroutine psi_dalqsrx_dw + +subroutine psi_dalqsr_up(n,x) + use psb_d_sort_mod, psb_protect_name => psi_dalqsr_up + use psb_error_mod + use psi_alcx_mod + implicit none + + real(psb_dpk_), intent(inout) :: x(:) + integer(psb_ipk_), intent(in) :: n + ! .. + ! .. Local Scalars .. + real(psb_dpk_) :: 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=60 + 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_dqsr',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_dalisr_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_dalisr_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_dalisr_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_dalisr_up(n1,x(ilx:i-1)) + endif + endif + enddo + else + call psi_dalisr_up(n,x) + endif +end subroutine psi_dalqsr_up + +subroutine psi_dalqsr_dw(n,x) + use psb_d_sort_mod, psb_protect_name => psi_dalqsr_dw + use psb_error_mod + use psi_alcx_mod + implicit none + + real(psb_dpk_), intent(inout) :: x(:) + integer(psb_ipk_), intent(in) :: n + ! .. Local Scalars .. + real(psb_dpk_) :: 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=60 + 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_dqsr',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_dalisr_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_dalisr_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_dalisr_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_dalisr_dw(n1,x(ilx:i-1)) + endif + endif + enddo + else + call psi_dalisr_dw(n,x) + endif +end subroutine psi_dalqsr_dw + +@CPLXE@ subroutine psi_daqsrx_up(n,x,idx) use psb_d_sort_mod, psb_protect_name => psi_daqsrx_up use psb_error_mod @@ -734,7 +1904,7 @@ subroutine psi_daqsrx_up(n,x,idx) integer(psb_ipk_) :: i, j, ilx, iux, istp, lpiv integer(psb_ipk_) :: ixt, n1, n2 - integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=16 + integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=60 integer(psb_ipk_) :: istack(nparms,maxstack) if (n > ithrs) then @@ -890,7 +2060,7 @@ subroutine psi_daqsrx_dw(n,x,idx) integer(psb_ipk_) :: i, j, ilx, iux, istp, lpiv integer(psb_ipk_) :: ixt, n1, n2 - integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=16 + integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=60 integer(psb_ipk_) :: istack(nparms,maxstack) if (n > ithrs) then ! @@ -1043,7 +2213,7 @@ subroutine psi_daqsr_up(n,x) integer(psb_ipk_) :: i, j, ilx, iux, istp, lpiv integer(psb_ipk_) :: ixt, n1, n2 - integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=16 + integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=60 integer(psb_ipk_) :: istack(nparms,maxstack) if (n > ithrs) then @@ -1183,7 +2353,7 @@ subroutine psi_daqsr_dw(n,x) integer(psb_ipk_) :: i, j, ilx, iux, istp, lpiv integer(psb_ipk_) :: ixt, n1, n2 - integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=16 + integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=60 integer(psb_ipk_) :: istack(nparms,maxstack) if (n > ithrs) then diff --git a/base/serial/sort/psb_i_qsort_impl.f90 b/base/serial/sort/psb_i_qsort_impl.f90 index f16e0e92..83aa43bd 100644 --- a/base/serial/sort/psb_i_qsort_impl.f90 +++ b/base/serial/sort/psb_i_qsort_impl.f90 @@ -141,8 +141,7 @@ subroutine psi_iqsrx_up(n,x,idx) integer(psb_ipk_) :: piv, xk, xt integer(psb_ipk_) :: i, j, ilx, iux, istp, lpiv integer(psb_ipk_) :: ixt, n1, n2 - - integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=16 + integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=72 integer(psb_ipk_) :: istack(nparms,maxstack) if (n > ithrs) then @@ -296,7 +295,7 @@ subroutine psi_iqsrx_dw(n,x,idx) integer(psb_ipk_) :: i, j, ilx, iux, istp, lpiv integer(psb_ipk_) :: ixt, n1, n2 - integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=16 + integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=72 integer(psb_ipk_) :: istack(nparms,maxstack) if (n > ithrs) then @@ -451,7 +450,7 @@ subroutine psi_iqsr_up(n,x) integer(psb_ipk_) :: i, j, ilx, iux, istp, lpiv integer(psb_ipk_) :: n1, n2 - integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=16 + integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=72 integer(psb_ipk_) :: istack(nparms,maxstack) @@ -592,7 +591,7 @@ subroutine psi_iqsr_dw(n,x) integer(psb_ipk_) :: i, j, ilx, iux, istp, lpiv integer(psb_ipk_) :: n1, n2 - integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=16 + integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=72 integer(psb_ipk_) :: istack(nparms,maxstack) @@ -734,7 +733,7 @@ subroutine psi_iaqsrx_up(n,x,idx) integer(psb_ipk_) :: i, j, ilx, iux, istp, lpiv integer(psb_ipk_) :: ixt, n1, n2 - integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=16 + integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=72 integer(psb_ipk_) :: istack(nparms,maxstack) if (n > ithrs) then @@ -890,7 +889,7 @@ subroutine psi_iaqsrx_dw(n,x,idx) integer(psb_ipk_) :: i, j, ilx, iux, istp, lpiv integer(psb_ipk_) :: ixt, n1, n2 - integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=16 + integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=72 integer(psb_ipk_) :: istack(nparms,maxstack) if (n > ithrs) then ! @@ -1043,7 +1042,7 @@ subroutine psi_iaqsr_up(n,x) integer(psb_ipk_) :: i, j, ilx, iux, istp, lpiv integer(psb_ipk_) :: ixt, n1, n2 - integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=16 + integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=72 integer(psb_ipk_) :: istack(nparms,maxstack) if (n > ithrs) then @@ -1183,7 +1182,7 @@ subroutine psi_iaqsr_dw(n,x) integer(psb_ipk_) :: i, j, ilx, iux, istp, lpiv integer(psb_ipk_) :: ixt, n1, n2 - integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=16 + integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=72 integer(psb_ipk_) :: istack(nparms,maxstack) if (n > ithrs) then diff --git a/base/serial/sort/psb_s_qsort_impl.f90 b/base/serial/sort/psb_s_qsort_impl.f90 index a07d3e2b..f8903551 100644 --- a/base/serial/sort/psb_s_qsort_impl.f90 +++ b/base/serial/sort/psb_s_qsort_impl.f90 @@ -141,8 +141,7 @@ subroutine psi_sqsrx_up(n,x,idx) real(psb_spk_) :: piv, xk, xt integer(psb_ipk_) :: i, j, ilx, iux, istp, lpiv integer(psb_ipk_) :: ixt, n1, n2 - - integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=16 + integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=72 integer(psb_ipk_) :: istack(nparms,maxstack) if (n > ithrs) then @@ -296,7 +295,7 @@ subroutine psi_sqsrx_dw(n,x,idx) integer(psb_ipk_) :: i, j, ilx, iux, istp, lpiv integer(psb_ipk_) :: ixt, n1, n2 - integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=16 + integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=72 integer(psb_ipk_) :: istack(nparms,maxstack) if (n > ithrs) then @@ -451,7 +450,7 @@ subroutine psi_sqsr_up(n,x) integer(psb_ipk_) :: i, j, ilx, iux, istp, lpiv integer(psb_ipk_) :: n1, n2 - integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=16 + integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=72 integer(psb_ipk_) :: istack(nparms,maxstack) @@ -592,7 +591,7 @@ subroutine psi_sqsr_dw(n,x) integer(psb_ipk_) :: i, j, ilx, iux, istp, lpiv integer(psb_ipk_) :: n1, n2 - integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=16 + integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=72 integer(psb_ipk_) :: istack(nparms,maxstack) @@ -720,6 +719,1177 @@ subroutine psi_sqsr_dw(n,x) end subroutine psi_sqsr_dw + 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_sqsrx',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_slisrx_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_slisrx_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_slisrx_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_slisrx_up(n1,x(ilx:i-1),idx(ilx:i-1)) + endif + endif + enddo + else + call psi_slisrx_up(n,x,idx) + endif + +end subroutine psi_slqsrx_up + +subroutine psi_slqsrx_dw(n,x,idx) + use psb_s_sort_mod, psb_protect_name => psi_slqsrx_dw + use psb_error_mod + use psi_lcx_mod + implicit none + + real(psb_spk_), intent(inout) :: x(:) + integer(psb_ipk_), intent(inout) :: idx(:) + integer(psb_ipk_), intent(in) :: n + ! .. Local Scalars .. + real(psb_spk_) :: piv, xk, xt + integer(psb_ipk_) :: i, j, ilx, iux, istp, lpiv + integer(psb_ipk_) :: ixt, 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) + 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_sqsrx',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_slisrx_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_slisrx_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_slisrx_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_slisrx_dw(n1,x(ilx:i-1),idx(ilx:i-1)) + endif + endif + enddo + else + call psi_slisrx_dw(n,x,idx) + endif +end subroutine psi_slqsrx_dw + +subroutine psi_slqsr_up(n,x) + use psb_s_sort_mod, psb_protect_name => psi_slqsr_up + use psb_error_mod + use psi_lcx_mod + implicit none + + real(psb_spk_), intent(inout) :: x(:) + integer(psb_ipk_), intent(in) :: n + ! .. + ! .. Local Scalars .. + real(psb_spk_) :: 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_sqsr',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_slisr_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_slisr_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_slisr_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_slisr_up(n1,x(ilx:i-1)) + endif + endif + enddo + else + call psi_slisr_up(n,x) + endif + +end subroutine psi_slqsr_up + +subroutine psi_slqsr_dw(n,x) + use psb_s_sort_mod, psb_protect_name => psi_slqsr_dw + use psb_error_mod + use psi_lcx_mod + implicit none + + real(psb_spk_), intent(inout) :: x(:) + integer(psb_ipk_), intent(in) :: n + ! .. Local Scalars .. + real(psb_spk_) :: 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_sqsr',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_slisr_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_slisr_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_slisr_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_slisr_dw(n1,x(ilx:i-1)) + endif + endif + enddo + else + call psi_slisr_dw(n,x) + endif + +end subroutine psi_slqsr_dw + +subroutine psi_salqsrx_up(n,x,idx) + use psb_s_sort_mod, psb_protect_name => psi_salqsrx_up + use psb_error_mod + use psi_alcx_mod + implicit none + + real(psb_spk_), intent(inout) :: x(:) + integer(psb_ipk_), intent(inout) :: idx(:) + integer(psb_ipk_), intent(in) :: n + ! .. Local Scalars .. + real(psb_spk_) :: piv, xk, xt + integer(psb_ipk_) :: i, j, ilx, iux, istp, lpiv + integer(psb_ipk_) :: ixt, 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) + 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_sqsrx',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_salisrx_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_salisrx_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_salisrx_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_salisrx_up(n1,x(ilx:i-1),idx(ilx:i-1)) + endif + endif + enddo + else + call psi_salisrx_up(n,x,idx) + endif +end subroutine psi_salqsrx_up + +subroutine psi_salqsrx_dw(n,x,idx) + use psb_s_sort_mod, psb_protect_name => psi_salqsrx_dw + use psb_error_mod + use psi_alcx_mod + implicit none + + real(psb_spk_), intent(inout) :: x(:) + integer(psb_ipk_), intent(inout) :: idx(:) + integer(psb_ipk_), intent(in) :: n + ! .. Local Scalars .. + real(psb_spk_) :: piv, xk, xt + integer(psb_ipk_) :: i, j, ilx, iux, istp, lpiv + integer(psb_ipk_) :: ixt, 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) + 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_sqsrx',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_salisrx_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_salisrx_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_salisrx_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_salisrx_dw(n1,x(ilx:i-1),idx(ilx:i-1)) + endif + endif + enddo + else + call psi_salisrx_dw(n,x,idx) + endif +end subroutine psi_salqsrx_dw + +subroutine psi_salqsr_up(n,x) + use psb_s_sort_mod, psb_protect_name => psi_salqsr_up + use psb_error_mod + use psi_alcx_mod + implicit none + + real(psb_spk_), intent(inout) :: x(:) + integer(psb_ipk_), intent(in) :: n + ! .. + ! .. Local Scalars .. + real(psb_spk_) :: 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_sqsr',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_salisr_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_salisr_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_salisr_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_salisr_up(n1,x(ilx:i-1)) + endif + endif + enddo + else + call psi_salisr_up(n,x) + endif +end subroutine psi_salqsr_up + +subroutine psi_salqsr_dw(n,x) + use psb_s_sort_mod, psb_protect_name => psi_salqsr_dw + use psb_error_mod + use psi_alcx_mod + implicit none + + real(psb_spk_), intent(inout) :: x(:) + integer(psb_ipk_), intent(in) :: n + ! .. Local Scalars .. + real(psb_spk_) :: 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_sqsr',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_salisr_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_salisr_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_salisr_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_salisr_dw(n1,x(ilx:i-1)) + endif + endif + enddo + else + call psi_salisr_dw(n,x) + endif +end subroutine psi_salqsr_dw + +@CPLXE@ subroutine psi_saqsrx_up(n,x,idx) use psb_s_sort_mod, psb_protect_name => psi_saqsrx_up use psb_error_mod @@ -734,7 +1904,7 @@ subroutine psi_saqsrx_up(n,x,idx) integer(psb_ipk_) :: i, j, ilx, iux, istp, lpiv integer(psb_ipk_) :: ixt, n1, n2 - integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=16 + integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=72 integer(psb_ipk_) :: istack(nparms,maxstack) if (n > ithrs) then @@ -890,7 +2060,7 @@ subroutine psi_saqsrx_dw(n,x,idx) integer(psb_ipk_) :: i, j, ilx, iux, istp, lpiv integer(psb_ipk_) :: ixt, n1, n2 - integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=16 + integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=72 integer(psb_ipk_) :: istack(nparms,maxstack) if (n > ithrs) then ! @@ -1043,7 +2213,7 @@ subroutine psi_saqsr_up(n,x) integer(psb_ipk_) :: i, j, ilx, iux, istp, lpiv integer(psb_ipk_) :: ixt, n1, n2 - integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=16 + integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=72 integer(psb_ipk_) :: istack(nparms,maxstack) if (n > ithrs) then @@ -1183,7 +2353,7 @@ subroutine psi_saqsr_dw(n,x) integer(psb_ipk_) :: i, j, ilx, iux, istp, lpiv integer(psb_ipk_) :: ixt, n1, n2 - integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=16 + integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=72 integer(psb_ipk_) :: istack(nparms,maxstack) if (n > ithrs) then diff --git a/base/serial/sort/psb_z_qsort_impl.f90 b/base/serial/sort/psb_z_qsort_impl.f90 index bf09e9bb..e1a4c0b7 100644 --- a/base/serial/sort/psb_z_qsort_impl.f90 +++ b/base/serial/sort/psb_z_qsort_impl.f90 @@ -137,6 +137,584 @@ subroutine psb_zqsort(x,ix,dir,flag) return end subroutine psb_zqsort + integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=24 + 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_zqsrx',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_zisrx_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_zisrx_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_zisrx_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_zisrx_up(n1,x(ilx:i-1),idx(ilx:i-1)) + endif + endif + enddo + else + call psi_zisrx_up(n,x,idx) + endif +end subroutine psi_zqsrx_up + +subroutine psi_zqsrx_dw(n,x,idx) + use psb_z_sort_mod, psb_protect_name => psi_zqsrx_dw + use psb_error_mod + implicit none + + complex(psb_dpk_), intent(inout) :: x(:) + integer(psb_ipk_), intent(inout) :: idx(:) + integer(psb_ipk_), intent(in) :: n + ! .. Local Scalars .. + complex(psb_dpk_) :: piv, xk, xt + integer(psb_ipk_) :: i, j, ilx, iux, istp, lpiv + integer(psb_ipk_) :: ixt, n1, n2 + + integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=24 + 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_zqsrx',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_zisrx_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_zisrx_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_zisrx_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_zisrx_dw(n1,x(ilx:i-1),idx(ilx:i-1)) + endif + endif + enddo + else + call psi_zisrx_dw(n,x,idx) + endif + +end subroutine psi_zqsrx_dw + +subroutine psi_zqsr_up(n,x) + use psb_z_sort_mod, psb_protect_name => psi_zqsr_up + use psb_error_mod + implicit none + + complex(psb_dpk_), intent(inout) :: x(:) + integer(psb_ipk_), intent(in) :: n + ! .. + ! .. Local Scalars .. + complex(psb_dpk_) :: 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=24 + 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_zqsr',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_zisr_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_zisr_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_zisr_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_zisr_up(n1,x(ilx:i-1)) + endif + endif + enddo + else + call psi_zisr_up(n,x) + endif + +end subroutine psi_zqsr_up + +subroutine psi_zqsr_dw(n,x) + use psb_z_sort_mod, psb_protect_name => psi_zqsr_dw + use psb_error_mod + implicit none + + complex(psb_dpk_), intent(inout) :: x(:) + integer(psb_ipk_), intent(in) :: n + ! .. + ! .. Local Scalars .. + complex(psb_dpk_) :: 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=24 + 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_zqsr',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_zisr_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_zisr_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_zisr_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_zisr_dw(n1,x(ilx:i-1)) + endif + endif + enddo + else + call psi_zisr_dw(n,x) + endif + +end subroutine psi_zqsr_dw +@REALE@ subroutine psi_zlqsrx_up(n,x,idx) use psb_z_sort_mod, psb_protect_name => psi_zlqsrx_up @@ -152,7 +730,7 @@ subroutine psi_zlqsrx_up(n,x,idx) integer(psb_ipk_) :: i, j, ilx, iux, istp, lpiv integer(psb_ipk_) :: ixt, n1, n2 - integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=16 + integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=24 integer(psb_ipk_) :: istack(nparms,maxstack) if (n > ithrs) then @@ -308,7 +886,7 @@ subroutine psi_zlqsrx_dw(n,x,idx) integer(psb_ipk_) :: i, j, ilx, iux, istp, lpiv integer(psb_ipk_) :: ixt, n1, n2 - integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=16 + integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=24 integer(psb_ipk_) :: istack(nparms,maxstack) if (n > ithrs) then @@ -463,7 +1041,7 @@ subroutine psi_zlqsr_up(n,x) integer(psb_ipk_) :: i, j, ilx, iux, istp, lpiv integer(psb_ipk_) :: n1, n2 - integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=16 + integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=24 integer(psb_ipk_) :: istack(nparms,maxstack) @@ -604,7 +1182,7 @@ subroutine psi_zlqsr_dw(n,x) integer(psb_ipk_) :: i, j, ilx, iux, istp, lpiv integer(psb_ipk_) :: n1, n2 - integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=16 + integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=24 integer(psb_ipk_) :: istack(nparms,maxstack) @@ -746,7 +1324,7 @@ subroutine psi_zalqsrx_up(n,x,idx) integer(psb_ipk_) :: i, j, ilx, iux, istp, lpiv integer(psb_ipk_) :: ixt, n1, n2 - integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=16 + integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=24 integer(psb_ipk_) :: istack(nparms,maxstack) if (n > ithrs) then @@ -901,7 +1479,7 @@ subroutine psi_zalqsrx_dw(n,x,idx) integer(psb_ipk_) :: i, j, ilx, iux, istp, lpiv integer(psb_ipk_) :: ixt, n1, n2 - integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=16 + integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=24 integer(psb_ipk_) :: istack(nparms,maxstack) if (n > ithrs) then @@ -1056,7 +1634,7 @@ subroutine psi_zalqsr_up(n,x) integer(psb_ipk_) :: i, j, ilx, iux, istp, lpiv integer(psb_ipk_) :: n1, n2 - integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=16 + integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=24 integer(psb_ipk_) :: istack(nparms,maxstack) @@ -1196,7 +1774,7 @@ subroutine psi_zalqsr_dw(n,x) integer(psb_ipk_) :: i, j, ilx, iux, istp, lpiv integer(psb_ipk_) :: n1, n2 - integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=16 + integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=24 integer(psb_ipk_) :: istack(nparms,maxstack) @@ -1337,7 +1915,7 @@ subroutine psi_zaqsrx_up(n,x,idx) integer(psb_ipk_) :: i, j, ilx, iux, istp, lpiv integer(psb_ipk_) :: ixt, n1, n2 - integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=16 + integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=24 integer(psb_ipk_) :: istack(nparms,maxstack) if (n > ithrs) then @@ -1493,7 +2071,7 @@ subroutine psi_zaqsrx_dw(n,x,idx) integer(psb_ipk_) :: i, j, ilx, iux, istp, lpiv integer(psb_ipk_) :: ixt, n1, n2 - integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=16 + integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=24 integer(psb_ipk_) :: istack(nparms,maxstack) if (n > ithrs) then ! @@ -1646,7 +2224,7 @@ subroutine psi_zaqsr_up(n,x) integer(psb_ipk_) :: i, j, ilx, iux, istp, lpiv integer(psb_ipk_) :: ixt, n1, n2 - integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=16 + integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=24 integer(psb_ipk_) :: istack(nparms,maxstack) if (n > ithrs) then @@ -1786,7 +2364,7 @@ subroutine psi_zaqsr_dw(n,x) integer(psb_ipk_) :: i, j, ilx, iux, istp, lpiv integer(psb_ipk_) :: ixt, n1, n2 - integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=16 + integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=24 integer(psb_ipk_) :: istack(nparms,maxstack) if (n > ithrs) then