diff --git a/base/serial/sort/psb_c_qsort_impl.f90 b/base/serial/sort/psb_c_qsort_impl.f90 index 1f81743f..cf5b7b0e 100644 --- a/base/serial/sort/psb_c_qsort_impl.f90 +++ b/base/serial/sort/psb_c_qsort_impl.f90 @@ -137,584 +137,6 @@ 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 @@ -730,7 +152,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=32 + integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=16 integer(psb_ipk_) :: istack(nparms,maxstack) if (n > ithrs) then @@ -886,7 +308,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=32 + integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=16 integer(psb_ipk_) :: istack(nparms,maxstack) if (n > ithrs) then @@ -1041,7 +463,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=32 + integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=16 integer(psb_ipk_) :: istack(nparms,maxstack) @@ -1182,7 +604,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=32 + integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=16 integer(psb_ipk_) :: istack(nparms,maxstack) @@ -1324,7 +746,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=32 + integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=16 integer(psb_ipk_) :: istack(nparms,maxstack) if (n > ithrs) then @@ -1479,7 +901,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=32 + integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=16 integer(psb_ipk_) :: istack(nparms,maxstack) if (n > ithrs) then @@ -1634,7 +1056,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=32 + integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=16 integer(psb_ipk_) :: istack(nparms,maxstack) @@ -1774,7 +1196,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=32 + integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=16 integer(psb_ipk_) :: istack(nparms,maxstack) @@ -1915,7 +1337,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=32 + integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=16 integer(psb_ipk_) :: istack(nparms,maxstack) if (n > ithrs) then @@ -2071,7 +1493,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=32 + integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=16 integer(psb_ipk_) :: istack(nparms,maxstack) if (n > ithrs) then ! @@ -2224,7 +1646,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=32 + integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=16 integer(psb_ipk_) :: istack(nparms,maxstack) if (n > ithrs) then @@ -2364,7 +1786,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=32 + integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=16 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 f3a40c06..c0eda8ed 100644 --- a/base/serial/sort/psb_d_qsort_impl.f90 +++ b/base/serial/sort/psb_d_qsort_impl.f90 @@ -141,7 +141,8 @@ 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=60 + + integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=16 integer(psb_ipk_) :: istack(nparms,maxstack) if (n > ithrs) then @@ -295,7 +296,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=60 + integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=16 integer(psb_ipk_) :: istack(nparms,maxstack) if (n > ithrs) then @@ -450,7 +451,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=60 + integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=16 integer(psb_ipk_) :: istack(nparms,maxstack) @@ -591,7 +592,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=60 + integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=16 integer(psb_ipk_) :: istack(nparms,maxstack) @@ -719,1177 +720,6 @@ 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 @@ -1904,7 +734,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=60 + integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=16 integer(psb_ipk_) :: istack(nparms,maxstack) if (n > ithrs) then @@ -2060,7 +890,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=60 + integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=16 integer(psb_ipk_) :: istack(nparms,maxstack) if (n > ithrs) then ! @@ -2213,7 +1043,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=60 + integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=16 integer(psb_ipk_) :: istack(nparms,maxstack) if (n > ithrs) then @@ -2353,7 +1183,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=60 + integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=16 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 83aa43bd..6ea51aba 100644 --- a/base/serial/sort/psb_i_qsort_impl.f90 +++ b/base/serial/sort/psb_i_qsort_impl.f90 @@ -141,6 +141,910 @@ 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_) :: 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_i_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_) :: ixt, n1, n2 + + integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=16 + 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_i_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=16 + 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_ilisrx_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_ilisrx_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_ilisrx_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_ilisrx_up(n1,x(ilx:i-1),idx(ilx:i-1)) + endif + endif + enddo + else + call psi_ilisrx_up(n,x,idx) + endif + +end subroutine psi_ilqsrx_up + +subroutine psi_ilqsrx_dw(n,x,idx) + use psb_i_sort_mod, psb_protect_name => psi_ilqsrx_dw + use psb_error_mod + use psi_lcx_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_) :: 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_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_ilisrx_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_ilisrx_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_ilisrx_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_ilisrx_dw(n1,x(ilx:i-1),idx(ilx:i-1)) + endif + endif + enddo + else + call psi_ilisrx_dw(n,x,idx) + endif +end subroutine psi_ilqsrx_dw + +subroutine psi_ilqsr_up(n,x) + use psb_i_sort_mod, psb_protect_name => psi_ilqsr_up + use psb_error_mod + use psi_lcx_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_ilisr_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_ilisr_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_ilisr_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_ilisr_up(n1,x(ilx:i-1)) + endif + endif + enddo + else + call psi_ilisr_up(n,x) + endif + +end subroutine psi_ilqsr_up + +subroutine psi_ilqsr_dw(n,x) + use psb_i_sort_mod, psb_protect_name => psi_ilqsr_dw + use psb_error_mod + use psi_lcx_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_ilisr_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_ilisr_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_ilisr_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_ilisr_dw(n1,x(ilx:i-1)) + endif + endif + enddo + else + call psi_ilisr_dw(n,x) + endif + +end subroutine psi_ilqsr_dw + +subroutine psi_ialqsrx_up(n,x,idx) + use psb_i_sort_mod, psb_protect_name => psi_ialqsrx_up + use psb_error_mod + use psi_alcx_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_) :: ixt, n1, n2 + integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=72 integer(psb_ipk_) :: istack(nparms,maxstack) @@ -251,14 +1155,14 @@ subroutine psi_iqsrx_up(n,x,idx) istack(1,istp) = ilx istack(2,istp) = i-1 else - call psi_iisrx_up(n1,x(ilx:i-1),idx(ilx:i-1)) + call psi_ialisrx_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)) + call psi_ialisrx_up(n2,x(i:iux),idx(i:iux)) endif else if (n2 > ithrs) then @@ -266,25 +1170,26 @@ subroutine psi_iqsrx_up(n,x,idx) istack(1,istp) = i istack(2,istp) = iux else - call psi_iisrx_up(n2,x(i:iux),idx(i:iux)) + call psi_ialisrx_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)) + call psi_ialisrx_up(n1,x(ilx:i-1),idx(ilx:i-1)) endif endif enddo else - call psi_iisrx_up(n,x,idx) + call psi_ialisrx_up(n,x,idx) endif -end subroutine psi_iqsrx_up +end subroutine psi_ialqsrx_up -subroutine psi_iqsrx_dw(n,x,idx) - use psb_i_sort_mod, psb_protect_name => psi_iqsrx_dw +subroutine psi_ialqsrx_dw(n,x,idx) + use psb_i_sort_mod, psb_protect_name => psi_ialqsrx_dw use psb_error_mod + use psi_alcx_mod implicit none integer(psb_ipk_), intent(inout) :: x(:) @@ -405,14 +1310,14 @@ subroutine psi_iqsrx_dw(n,x,idx) istack(1,istp) = ilx istack(2,istp) = i-1 else - call psi_iisrx_dw(n1,x(ilx:i-1),idx(ilx:i-1)) + call psi_ialisrx_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)) + call psi_ialisrx_dw(n2,x(i:iux),idx(i:iux)) endif else if (n2 > ithrs) then @@ -420,26 +1325,26 @@ subroutine psi_iqsrx_dw(n,x,idx) istack(1,istp) = i istack(2,istp) = iux else - call psi_iisrx_dw(n2,x(i:iux),idx(i:iux)) + call psi_ialisrx_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)) + call psi_ialisrx_dw(n1,x(ilx:i-1),idx(ilx:i-1)) endif endif enddo else - call psi_iisrx_dw(n,x,idx) + call psi_ialisrx_dw(n,x,idx) endif +end subroutine psi_ialqsrx_dw -end subroutine psi_iqsrx_dw - -subroutine psi_iqsr_up(n,x) - use psb_i_sort_mod, psb_protect_name => psi_iqsr_up +subroutine psi_ialqsr_up(n,x) + use psb_i_sort_mod, psb_protect_name => psi_ialqsr_up use psb_error_mod + use psi_alcx_mod implicit none integer(psb_ipk_), intent(inout) :: x(:) @@ -450,7 +1355,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=72 + integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=16 integer(psb_ipk_) :: istack(nparms,maxstack) @@ -546,14 +1451,14 @@ subroutine psi_iqsr_up(n,x) istack(1,istp) = ilx istack(2,istp) = i-1 else - call psi_iisr_up(n1,x(ilx:i-1)) + call psi_ialisr_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)) + call psi_ialisr_up(n2,x(i:iux)) endif else if (n2 > ithrs) then @@ -561,31 +1466,30 @@ subroutine psi_iqsr_up(n,x) istack(1,istp) = i istack(2,istp) = iux else - call psi_iisr_up(n2,x(i:iux)) + call psi_ialisr_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)) + call psi_ialisr_up(n1,x(ilx:i-1)) endif endif enddo else - call psi_iisr_up(n,x) + call psi_ialisr_up(n,x) endif +end subroutine psi_ialqsr_up -end subroutine psi_iqsr_up - -subroutine psi_iqsr_dw(n,x) - use psb_i_sort_mod, psb_protect_name => psi_iqsr_dw +subroutine psi_ialqsr_dw(n,x) + use psb_i_sort_mod, psb_protect_name => psi_ialqsr_dw use psb_error_mod + use psi_alcx_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 @@ -687,14 +1591,14 @@ subroutine psi_iqsr_dw(n,x) istack(1,istp) = ilx istack(2,istp) = i-1 else - call psi_iisr_dw(n1,x(ilx:i-1)) + call psi_ialisr_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)) + call psi_ialisr_dw(n2,x(i:iux)) endif else if (n2 > ithrs) then @@ -702,23 +1606,23 @@ subroutine psi_iqsr_dw(n,x) istack(1,istp) = i istack(2,istp) = iux else - call psi_iisr_dw(n2,x(i:iux)) + call psi_ialisr_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)) + call psi_ialisr_dw(n1,x(ilx:i-1)) endif endif enddo else - call psi_iisr_dw(n,x) + call psi_ialisr_dw(n,x) endif +end subroutine psi_ialqsr_dw -end subroutine psi_iqsr_dw - +@CPLXE@ subroutine psi_iaqsrx_up(n,x,idx) use psb_i_sort_mod, psb_protect_name => psi_iaqsrx_up use psb_error_mod @@ -733,7 +1637,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=72 + integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=16 integer(psb_ipk_) :: istack(nparms,maxstack) if (n > ithrs) then @@ -889,7 +1793,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=72 + integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=16 integer(psb_ipk_) :: istack(nparms,maxstack) if (n > ithrs) then ! @@ -1042,7 +1946,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=72 + integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=16 integer(psb_ipk_) :: istack(nparms,maxstack) if (n > ithrs) then @@ -1182,7 +2086,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=72 + integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=16 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 f8903551..a07d3e2b 100644 --- a/base/serial/sort/psb_s_qsort_impl.f90 +++ b/base/serial/sort/psb_s_qsort_impl.f90 @@ -141,7 +141,8 @@ 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=72 + + integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=16 integer(psb_ipk_) :: istack(nparms,maxstack) if (n > ithrs) then @@ -295,7 +296,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=72 + integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=16 integer(psb_ipk_) :: istack(nparms,maxstack) if (n > ithrs) then @@ -450,7 +451,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=72 + integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=16 integer(psb_ipk_) :: istack(nparms,maxstack) @@ -591,7 +592,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=72 + integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=16 integer(psb_ipk_) :: istack(nparms,maxstack) @@ -719,1177 +720,6 @@ 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 @@ -1904,7 +734,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=72 + integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=16 integer(psb_ipk_) :: istack(nparms,maxstack) if (n > ithrs) then @@ -2060,7 +890,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=72 + integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=16 integer(psb_ipk_) :: istack(nparms,maxstack) if (n > ithrs) then ! @@ -2213,7 +1043,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=72 + integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=16 integer(psb_ipk_) :: istack(nparms,maxstack) if (n > ithrs) then @@ -2353,7 +1183,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=72 + integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=16 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 e1a4c0b7..bf09e9bb 100644 --- a/base/serial/sort/psb_z_qsort_impl.f90 +++ b/base/serial/sort/psb_z_qsort_impl.f90 @@ -137,584 +137,6 @@ 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 @@ -730,7 +152,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=24 + integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=16 integer(psb_ipk_) :: istack(nparms,maxstack) if (n > ithrs) then @@ -886,7 +308,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=24 + integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=16 integer(psb_ipk_) :: istack(nparms,maxstack) if (n > ithrs) then @@ -1041,7 +463,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=24 + integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=16 integer(psb_ipk_) :: istack(nparms,maxstack) @@ -1182,7 +604,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=24 + integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=16 integer(psb_ipk_) :: istack(nparms,maxstack) @@ -1324,7 +746,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=24 + integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=16 integer(psb_ipk_) :: istack(nparms,maxstack) if (n > ithrs) then @@ -1479,7 +901,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=24 + integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=16 integer(psb_ipk_) :: istack(nparms,maxstack) if (n > ithrs) then @@ -1634,7 +1056,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=24 + integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=16 integer(psb_ipk_) :: istack(nparms,maxstack) @@ -1774,7 +1196,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=24 + integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=16 integer(psb_ipk_) :: istack(nparms,maxstack) @@ -1915,7 +1337,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=24 + integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=16 integer(psb_ipk_) :: istack(nparms,maxstack) if (n > ithrs) then @@ -2071,7 +1493,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=24 + integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=16 integer(psb_ipk_) :: istack(nparms,maxstack) if (n > ithrs) then ! @@ -2224,7 +1646,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=24 + integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=16 integer(psb_ipk_) :: istack(nparms,maxstack) if (n > ithrs) then @@ -2364,7 +1786,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=24 + integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=16 integer(psb_ipk_) :: istack(nparms,maxstack) if (n > ithrs) then