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

psblas3:

Browse Source 
base/serial/sort/psb_c_qsort_impl.f90
 base/serial/sort/psb_d_qsort_impl.f90
 base/serial/sort/psb_i_qsort_impl.f90
 base/serial/sort/psb_s_qsort_impl.f90
 base/serial/sort/psb_z_qsort_impl.f90

Revert qsort changes: need to sort out the correct template.
psblas3-dense
Salvatore Filippone 9 years ago
parent adddf964d3
commit 206a765adf
5 changed files with 985 additions and 3577 deletions
Show Stats Download Patch File Download Diff File

602
base/serial/sort/psb_c_qsort_impl.f90
Unescape Escape View File

@ -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

1188
base/serial/sort/psb_d_qsort_impl.f90
View File

File diff suppressed because it is too large Load Diff

982
base/serial/sort/psb_i_qsort_impl.f90
View File

File diff suppressed because it is too large Load Diff

1188
base/serial/sort/psb_s_qsort_impl.f90
View File

File diff suppressed because it is too large Load Diff

602
base/serial/sort/psb_z_qsort_impl.f90
Unescape Escape View File

@ -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

Write Preview
Loading…
Cancel
Save