diff --git a/base/modules/Makefile b/base/modules/Makefile index ec53f942c..19b43558a 100644 --- a/base/modules/Makefile +++ b/base/modules/Makefile @@ -27,7 +27,11 @@ UTIL_MODS = aux/psb_string_mod.o desc/psb_desc_const_mod.o desc/psb_indx_map_mod aux/psi_s_serial_mod.o aux/psi_d_serial_mod.o aux/psi_c_serial_mod.o aux/psi_z_serial_mod.o \ psi_mod.o psi_i_mod.o psi_s_mod.o psi_d_mod.o psi_c_mod.o psi_z_mod.o\ aux/psb_ip_reord_mod.o\ - aux/psb_i_sort_mod.o aux/psb_s_sort_mod.o aux/psb_d_sort_mod.o \ + aux/psb_i_ip_reord_mod.o aux/psb_l_ip_reord_mod.o \ + aux/psb_s_ip_reord_mod.o aux/psb_d_ip_reord_mod.o \ + aux/psb_c_ip_reord_mod.o aux/psb_z_ip_reord_mod.o \ + aux/psb_i_sort_mod.o aux/psb_l_sort_mod.o \ + aux/psb_s_sort_mod.o aux/psb_d_sort_mod.o \ aux/psb_c_sort_mod.o aux/psb_z_sort_mod.o \ psb_check_mod.o aux/psb_hash_mod.o\ serial/psb_base_mat_mod.o serial/psb_mat_mod.o\ @@ -64,12 +68,18 @@ aux/psb_string_mod.o desc/psb_desc_const_mod.o psi_comm_buffers_mod.o: psb_const aux/psb_hash_mod.o: psb_realloc_mod.o psb_const_mod.o aux/psb_i_sort_mod.o aux/psb_s_sort_mod.o aux/psb_d_sort_mod.o aux/psb_c_sort_mod.o aux/psb_z_sort_mod.o \ aux/psb_ip_reord_mod.o aux/psi_serial_mod.o aux/psb_sort_mod.o: $(BASIC_MODS) -aux/psb_sort_mod.o: aux/psb_i_sort_mod.o aux/psb_s_sort_mod.o aux/psb_d_sort_mod.o \ - aux/psb_c_sort_mod.o aux/psb_z_sort_mod.o aux/psb_ip_reord_mod.o aux/psi_serial_mod.o +aux/psb_sort_mod.o: aux/psb_i_sort_mod.o aux/psb_l_sort_mod.o \ + aux/psb_s_sort_mod.o aux/psb_d_sort_mod.o \ + aux/psb_c_sort_mod.o aux/psb_z_sort_mod.o \ + aux/psb_ip_reord_mod.o aux/psi_serial_mod.o aux/psi_serial_mod.o: aux/psi_i_serial_mod.o \ aux/psi_s_serial_mod.o aux/psi_d_serial_mod.o aux/psi_c_serial_mod.o aux/psi_z_serial_mod.o aux/psi_i_serial_mod.o aux/psi_s_serial_mod.o aux/psi_d_serial_mod.o aux/psi_c_serial_mod.o aux/psi_z_serial_mod.o: psb_const_mod.o +aux/psb_ip_reord_mod.o: aux/psb_i_ip_reord_mod.o aux/psb_l_ip_reord_mod.o \ + aux/psb_s_ip_reord_mod.o aux/psb_d_ip_reord_mod.o \ + aux/psb_c_ip_reord_mod.o aux/psb_z_ip_reord_mod.o + serial/psb_base_mat_mod.o: aux/psi_serial_mod.o serial/psb_s_base_mat_mod.o serial/psb_d_base_mat_mod.o serial/psb_c_base_mat_mod.o serial/psb_z_base_mat_mod.o: serial/psb_base_mat_mod.o serial/psb_s_base_mat_mod.o: serial/psb_s_base_vect_mod.o diff --git a/base/modules/aux/psb_c_ip_reord_mod.f90 b/base/modules/aux/psb_c_ip_reord_mod.f90 new file mode 100644 index 000000000..c31643fbf --- /dev/null +++ b/base/modules/aux/psb_c_ip_reord_mod.f90 @@ -0,0 +1,187 @@ +! +! Parallel Sparse BLAS version 3.5 +! (C) Copyright 2006-2018 +! Salvatore Filippone +! Alfredo Buttari +! +! Redistribution and use in source and binary forms, with or without +! modification, are permitted provided that the following conditions +! are met: +! 1. Redistributions of source code must retain the above copyright +! notice, this list of conditions and the following disclaimer. +! 2. Redistributions in binary form must reproduce the above copyright +! notice, this list of conditions, and the following disclaimer in the +! documentation and/or other materials provided with the distribution. +! 3. The name of the PSBLAS group or the names of its contributors may +! not be used to endorse or promote products derived from this +! software without specific written permission. +! +! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +! ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED +! TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR +! PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE PSBLAS GROUP OR ITS CONTRIBUTORS +! BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +! CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +! SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +! INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +! CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +! ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +! POSSIBILITY OF SUCH DAMAGE. +! +! +! +! Reorder (an) input vector(s) based on a list sort output. +! Based on: D. E. Knuth: The Art of Computer Programming +! vol. 3: Sorting and Searching, Addison Wesley, 1973 +! ex. 5.2.12 +! +! +module psb_c_ip_reord_mod + use psb_const_mod + + interface psb_ip_reord + module procedure psb_ip_reord_c1,& + & psb_ip_reord_c1i1, psb_ip_reord_c1i2,& + & psb_ip_reord_c1i3 + + + end interface + +contains + + subroutine psb_ip_reord_c1(n,x,iaux) + integer(psb_ipk_), intent(in) :: n + integer(psb_ipk_) :: iaux(0:*) + complex(psb_spk_) :: x(*) + + integer(psb_ipk_) :: lswap, lp, k + complex(psb_spk_) :: swap + + lp = iaux(0) + k = 1 + do + if ((lp == 0).or.(k>n)) exit + do + if (lp >= k) exit + lp = iaux(lp) + end do + swap = x(lp) + x(lp) = x(k) + x(k) = swap + lswap = iaux(lp) + iaux(lp) = iaux(k) + iaux(k) = lp + lp = lswap + k = k + 1 + enddo + return + end subroutine psb_ip_reord_c1 + + subroutine psb_ip_reord_c1i1(n,x,indx,iaux) + integer(psb_ipk_), intent(in) :: n + integer(psb_ipk_) :: iaux(0:*) + complex(psb_spk_) :: x(*) + integer(psb_ipk_) :: indx(*) + + integer(psb_ipk_) :: lswap, lp, k, ixswap + complex(psb_spk_) :: swap + + lp = iaux(0) + k = 1 + do + if ((lp == 0).or.(k>n)) exit + do + if (lp >= k) exit + lp = iaux(lp) + end do + swap = x(lp) + x(lp) = x(k) + x(k) = swap + ixswap = indx(lp) + indx(lp) = indx(k) + indx(k) = ixswap + lswap = iaux(lp) + iaux(lp) = iaux(k) + iaux(k) = lp + lp = lswap + k = k + 1 + enddo + return + end subroutine psb_ip_reord_c1i1 + + subroutine psb_ip_reord_c1i2(n,x,i1,i2,iaux) + integer(psb_ipk_), intent(in) :: n + integer(psb_ipk_) :: iaux(0:*) + complex(psb_spk_) :: x(*) + integer(psb_ipk_) :: i1(*), i2(*) + + + integer(psb_ipk_) :: lswap, lp, k, isw1, isw2 + complex(psb_spk_) :: swap + + lp = iaux(0) + k = 1 + do + if ((lp == 0).or.(k>n)) exit + do + if (lp >= k) exit + lp = iaux(lp) + end do + swap = x(lp) + x(lp) = x(k) + x(k) = swap + isw1 = i1(lp) + i1(lp) = i1(k) + i1(k) = isw1 + isw2 = i2(lp) + i2(lp) = i2(k) + i2(k) = isw2 + lswap = iaux(lp) + iaux(lp) = iaux(k) + iaux(k) = lp + lp = lswap + k = k + 1 + enddo + return + end subroutine psb_ip_reord_c1i2 + + subroutine psb_ip_reord_c1i3(n,x,i1,i2,i3,iaux) + integer(psb_ipk_), intent(in) :: n + integer(psb_ipk_) :: iaux(0:*) + complex(psb_spk_) :: x(*) + integer(psb_ipk_) :: i1(*), i2(*), i3(*) + + + integer(psb_ipk_) :: lswap, lp, k, isw1, isw2, isw3 + complex(psb_spk_) :: swap + + lp = iaux(0) + k = 1 + do + if ((lp == 0).or.(k>n)) exit + do + if (lp >= k) exit + lp = iaux(lp) + end do + swap = x(lp) + x(lp) = x(k) + x(k) = swap + isw1 = i1(lp) + i1(lp) = i1(k) + i1(k) = isw1 + isw2 = i2(lp) + i2(lp) = i2(k) + i2(k) = isw2 + isw3 = i3(lp) + i3(lp) = i3(k) + i3(k) = isw3 + lswap = iaux(lp) + iaux(lp) = iaux(k) + iaux(k) = lp + lp = lswap + k = k + 1 + enddo + return + end subroutine psb_ip_reord_c1i3 + +end module psb_c_ip_reord_mod diff --git a/base/modules/aux/psb_c_sort_mod.f90 b/base/modules/aux/psb_c_sort_mod.f90 index 53e98a843..a2f0f3c6b 100644 --- a/base/modules/aux/psb_c_sort_mod.f90 +++ b/base/modules/aux/psb_c_sort_mod.f90 @@ -44,6 +44,15 @@ module psb_c_sort_mod use psb_const_mod + interface psb_msort_unique + subroutine psb_cmsort_u(x,nout,dir) + import + complex(psb_spk_), intent(inout) :: x(:) + integer(psb_ipk_), intent(out) :: nout + integer(psb_ipk_), optional, intent(in) :: dir + end subroutine psb_cmsort_u + end interface psb_msort_unique + type psb_c_heap integer(psb_ipk_) :: last, dir complex(psb_spk_), allocatable :: keys(:) @@ -544,7 +553,8 @@ contains implicit none class(psb_c_idx_heap), intent(inout) :: heap - integer(psb_ipk_), intent(out) :: index,info + integer(psb_ipk_), intent(out) :: index + integer(psb_ipk_), intent(out) :: info complex(psb_spk_), intent(out) :: key diff --git a/base/modules/aux/psb_d_ip_reord_mod.f90 b/base/modules/aux/psb_d_ip_reord_mod.f90 new file mode 100644 index 000000000..87be4333a --- /dev/null +++ b/base/modules/aux/psb_d_ip_reord_mod.f90 @@ -0,0 +1,187 @@ +! +! Parallel Sparse BLAS version 3.5 +! (C) Copyright 2006-2018 +! Salvatore Filippone +! Alfredo Buttari +! +! Redistribution and use in source and binary forms, with or without +! modification, are permitted provided that the following conditions +! are met: +! 1. Redistributions of source code must retain the above copyright +! notice, this list of conditions and the following disclaimer. +! 2. Redistributions in binary form must reproduce the above copyright +! notice, this list of conditions, and the following disclaimer in the +! documentation and/or other materials provided with the distribution. +! 3. The name of the PSBLAS group or the names of its contributors may +! not be used to endorse or promote products derived from this +! software without specific written permission. +! +! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +! ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED +! TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR +! PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE PSBLAS GROUP OR ITS CONTRIBUTORS +! BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +! CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +! SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +! INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +! CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +! ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +! POSSIBILITY OF SUCH DAMAGE. +! +! +! +! Reorder (an) input vector(s) based on a list sort output. +! Based on: D. E. Knuth: The Art of Computer Programming +! vol. 3: Sorting and Searching, Addison Wesley, 1973 +! ex. 5.2.12 +! +! +module psb_d_ip_reord_mod + use psb_const_mod + + interface psb_ip_reord + module procedure psb_ip_reord_d1,& + & psb_ip_reord_d1i1, psb_ip_reord_d1i2,& + & psb_ip_reord_d1i3 + + + end interface + +contains + + subroutine psb_ip_reord_d1(n,x,iaux) + integer(psb_ipk_), intent(in) :: n + integer(psb_ipk_) :: iaux(0:*) + real(psb_dpk_) :: x(*) + + integer(psb_ipk_) :: lswap, lp, k + real(psb_dpk_) :: swap + + lp = iaux(0) + k = 1 + do + if ((lp == 0).or.(k>n)) exit + do + if (lp >= k) exit + lp = iaux(lp) + end do + swap = x(lp) + x(lp) = x(k) + x(k) = swap + lswap = iaux(lp) + iaux(lp) = iaux(k) + iaux(k) = lp + lp = lswap + k = k + 1 + enddo + return + end subroutine psb_ip_reord_d1 + + subroutine psb_ip_reord_d1i1(n,x,indx,iaux) + integer(psb_ipk_), intent(in) :: n + integer(psb_ipk_) :: iaux(0:*) + real(psb_dpk_) :: x(*) + integer(psb_ipk_) :: indx(*) + + integer(psb_ipk_) :: lswap, lp, k, ixswap + real(psb_dpk_) :: swap + + lp = iaux(0) + k = 1 + do + if ((lp == 0).or.(k>n)) exit + do + if (lp >= k) exit + lp = iaux(lp) + end do + swap = x(lp) + x(lp) = x(k) + x(k) = swap + ixswap = indx(lp) + indx(lp) = indx(k) + indx(k) = ixswap + lswap = iaux(lp) + iaux(lp) = iaux(k) + iaux(k) = lp + lp = lswap + k = k + 1 + enddo + return + end subroutine psb_ip_reord_d1i1 + + subroutine psb_ip_reord_d1i2(n,x,i1,i2,iaux) + integer(psb_ipk_), intent(in) :: n + integer(psb_ipk_) :: iaux(0:*) + real(psb_dpk_) :: x(*) + integer(psb_ipk_) :: i1(*), i2(*) + + + integer(psb_ipk_) :: lswap, lp, k, isw1, isw2 + real(psb_dpk_) :: swap + + lp = iaux(0) + k = 1 + do + if ((lp == 0).or.(k>n)) exit + do + if (lp >= k) exit + lp = iaux(lp) + end do + swap = x(lp) + x(lp) = x(k) + x(k) = swap + isw1 = i1(lp) + i1(lp) = i1(k) + i1(k) = isw1 + isw2 = i2(lp) + i2(lp) = i2(k) + i2(k) = isw2 + lswap = iaux(lp) + iaux(lp) = iaux(k) + iaux(k) = lp + lp = lswap + k = k + 1 + enddo + return + end subroutine psb_ip_reord_d1i2 + + subroutine psb_ip_reord_d1i3(n,x,i1,i2,i3,iaux) + integer(psb_ipk_), intent(in) :: n + integer(psb_ipk_) :: iaux(0:*) + real(psb_dpk_) :: x(*) + integer(psb_ipk_) :: i1(*), i2(*), i3(*) + + + integer(psb_ipk_) :: lswap, lp, k, isw1, isw2, isw3 + real(psb_dpk_) :: swap + + lp = iaux(0) + k = 1 + do + if ((lp == 0).or.(k>n)) exit + do + if (lp >= k) exit + lp = iaux(lp) + end do + swap = x(lp) + x(lp) = x(k) + x(k) = swap + isw1 = i1(lp) + i1(lp) = i1(k) + i1(k) = isw1 + isw2 = i2(lp) + i2(lp) = i2(k) + i2(k) = isw2 + isw3 = i3(lp) + i3(lp) = i3(k) + i3(k) = isw3 + lswap = iaux(lp) + iaux(lp) = iaux(k) + iaux(k) = lp + lp = lswap + k = k + 1 + enddo + return + end subroutine psb_ip_reord_d1i3 + +end module psb_d_ip_reord_mod diff --git a/base/modules/aux/psb_d_sort_mod.f90 b/base/modules/aux/psb_d_sort_mod.f90 index 72a83355b..4e9daee0c 100644 --- a/base/modules/aux/psb_d_sort_mod.f90 +++ b/base/modules/aux/psb_d_sort_mod.f90 @@ -44,6 +44,15 @@ module psb_d_sort_mod use psb_const_mod + interface psb_msort_unique + subroutine psb_dmsort_u(x,nout,dir) + import + real(psb_dpk_), intent(inout) :: x(:) + integer(psb_ipk_), intent(out) :: nout + integer(psb_ipk_), optional, intent(in) :: dir + end subroutine psb_dmsort_u + end interface psb_msort_unique + type psb_d_heap integer(psb_ipk_) :: last, dir real(psb_dpk_), allocatable :: keys(:) @@ -79,6 +88,26 @@ module psb_d_sort_mod end subroutine psb_dmsort end interface psb_msort + + interface psb_bsrch + function psb_dbsrch(key,n,v) result(ipos) + import + integer(psb_ipk_) :: ipos, n + real(psb_dpk_) :: key + real(psb_dpk_) :: v(:) + end function psb_dbsrch + end interface psb_bsrch + + interface psb_ssrch + function psb_dssrch(key,n,v) result(ipos) + import + implicit none + integer(psb_ipk_) :: ipos, n + real(psb_dpk_) :: key + real(psb_dpk_) :: v(:) + end function psb_dssrch + end interface psb_ssrch + interface subroutine psi_d_msort_up(n,k,l,iret) import @@ -486,7 +515,8 @@ contains implicit none class(psb_d_idx_heap), intent(inout) :: heap - integer(psb_ipk_), intent(out) :: index,info + integer(psb_ipk_), intent(out) :: index + integer(psb_ipk_), intent(out) :: info real(psb_dpk_), intent(out) :: key diff --git a/base/modules/aux/psb_i_ip_reord_mod.f90 b/base/modules/aux/psb_i_ip_reord_mod.f90 new file mode 100644 index 000000000..20ac6f643 --- /dev/null +++ b/base/modules/aux/psb_i_ip_reord_mod.f90 @@ -0,0 +1,187 @@ +! +! Parallel Sparse BLAS version 3.5 +! (C) Copyright 2006-2018 +! Salvatore Filippone +! Alfredo Buttari +! +! Redistribution and use in source and binary forms, with or without +! modification, are permitted provided that the following conditions +! are met: +! 1. Redistributions of source code must retain the above copyright +! notice, this list of conditions and the following disclaimer. +! 2. Redistributions in binary form must reproduce the above copyright +! notice, this list of conditions, and the following disclaimer in the +! documentation and/or other materials provided with the distribution. +! 3. The name of the PSBLAS group or the names of its contributors may +! not be used to endorse or promote products derived from this +! software without specific written permission. +! +! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +! ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED +! TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR +! PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE PSBLAS GROUP OR ITS CONTRIBUTORS +! BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +! CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +! SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +! INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +! CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +! ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +! POSSIBILITY OF SUCH DAMAGE. +! +! +! +! Reorder (an) input vector(s) based on a list sort output. +! Based on: D. E. Knuth: The Art of Computer Programming +! vol. 3: Sorting and Searching, Addison Wesley, 1973 +! ex. 5.2.12 +! +! +module psb_i_ip_reord_mod + use psb_const_mod + + interface psb_ip_reord + module procedure psb_ip_reord_i1,& + & psb_ip_reord_i1i1, psb_ip_reord_i1i2,& + & psb_ip_reord_i1i3 + + + end interface + +contains + + subroutine psb_ip_reord_i1(n,x,iaux) + integer(psb_ipk_), intent(in) :: n + integer(psb_ipk_) :: iaux(0:*) + integer(psb_ipk_) :: x(*) + + integer(psb_ipk_) :: lswap, lp, k + integer(psb_ipk_) :: swap + + lp = iaux(0) + k = 1 + do + if ((lp == 0).or.(k>n)) exit + do + if (lp >= k) exit + lp = iaux(lp) + end do + swap = x(lp) + x(lp) = x(k) + x(k) = swap + lswap = iaux(lp) + iaux(lp) = iaux(k) + iaux(k) = lp + lp = lswap + k = k + 1 + enddo + return + end subroutine psb_ip_reord_i1 + + subroutine psb_ip_reord_i1i1(n,x,indx,iaux) + integer(psb_ipk_), intent(in) :: n + integer(psb_ipk_) :: iaux(0:*) + integer(psb_ipk_) :: x(*) + integer(psb_ipk_) :: indx(*) + + integer(psb_ipk_) :: lswap, lp, k, ixswap + integer(psb_ipk_) :: swap + + lp = iaux(0) + k = 1 + do + if ((lp == 0).or.(k>n)) exit + do + if (lp >= k) exit + lp = iaux(lp) + end do + swap = x(lp) + x(lp) = x(k) + x(k) = swap + ixswap = indx(lp) + indx(lp) = indx(k) + indx(k) = ixswap + lswap = iaux(lp) + iaux(lp) = iaux(k) + iaux(k) = lp + lp = lswap + k = k + 1 + enddo + return + end subroutine psb_ip_reord_i1i1 + + subroutine psb_ip_reord_i1i2(n,x,i1,i2,iaux) + integer(psb_ipk_), intent(in) :: n + integer(psb_ipk_) :: iaux(0:*) + integer(psb_ipk_) :: x(*) + integer(psb_ipk_) :: i1(*), i2(*) + + + integer(psb_ipk_) :: lswap, lp, k, isw1, isw2 + integer(psb_ipk_) :: swap + + lp = iaux(0) + k = 1 + do + if ((lp == 0).or.(k>n)) exit + do + if (lp >= k) exit + lp = iaux(lp) + end do + swap = x(lp) + x(lp) = x(k) + x(k) = swap + isw1 = i1(lp) + i1(lp) = i1(k) + i1(k) = isw1 + isw2 = i2(lp) + i2(lp) = i2(k) + i2(k) = isw2 + lswap = iaux(lp) + iaux(lp) = iaux(k) + iaux(k) = lp + lp = lswap + k = k + 1 + enddo + return + end subroutine psb_ip_reord_i1i2 + + subroutine psb_ip_reord_i1i3(n,x,i1,i2,i3,iaux) + integer(psb_ipk_), intent(in) :: n + integer(psb_ipk_) :: iaux(0:*) + integer(psb_ipk_) :: x(*) + integer(psb_ipk_) :: i1(*), i2(*), i3(*) + + + integer(psb_ipk_) :: lswap, lp, k, isw1, isw2, isw3 + integer(psb_ipk_) :: swap + + lp = iaux(0) + k = 1 + do + if ((lp == 0).or.(k>n)) exit + do + if (lp >= k) exit + lp = iaux(lp) + end do + swap = x(lp) + x(lp) = x(k) + x(k) = swap + isw1 = i1(lp) + i1(lp) = i1(k) + i1(k) = isw1 + isw2 = i2(lp) + i2(lp) = i2(k) + i2(k) = isw2 + isw3 = i3(lp) + i3(lp) = i3(k) + i3(k) = isw3 + lswap = iaux(lp) + iaux(lp) = iaux(k) + iaux(k) = lp + lp = lswap + k = k + 1 + enddo + return + end subroutine psb_ip_reord_i1i3 + +end module psb_i_ip_reord_mod diff --git a/base/modules/aux/psb_i_sort_mod.f90 b/base/modules/aux/psb_i_sort_mod.f90 index 29c94c8fb..a7b13682a 100644 --- a/base/modules/aux/psb_i_sort_mod.f90 +++ b/base/modules/aux/psb_i_sort_mod.f90 @@ -43,42 +43,17 @@ module psb_i_sort_mod use psb_const_mod - interface psb_iblsrch - function psb_iblsrch(key,n,v) result(ipos) - import :: psb_ipk_ - integer(psb_ipk_) :: ipos, key, n - integer(psb_ipk_) :: v(:) - end function psb_iblsrch - end interface psb_iblsrch - - interface psb_ibsrch - function psb_ibsrch(key,n,v) result(ipos) - import :: psb_ipk_ - integer(psb_ipk_) :: ipos, key, n - integer(psb_ipk_) :: v(:) - end function psb_ibsrch - end interface psb_ibsrch - - interface psb_issrch - function psb_issrch(key,n,v) result(ipos) - import :: psb_ipk_ - implicit none - integer(psb_ipk_) :: ipos, key, n - integer(psb_ipk_) :: v(:) - end function psb_issrch - end interface psb_issrch - interface psb_isaperm - logical function psb_isaperm(n,eip) - import :: psb_ipk_ + logical function psb_iisaperm(n,eip) + import integer(psb_ipk_), intent(in) :: n integer(psb_ipk_), intent(in) :: eip(n) - end function psb_isaperm + end function psb_iisaperm end interface psb_isaperm interface psb_msort_unique subroutine psb_imsort_u(x,nout,dir) - import :: psb_ipk_, psb_spk_, psb_dpk_ + import integer(psb_ipk_), intent(inout) :: x(:) integer(psb_ipk_), intent(out) :: nout integer(psb_ipk_), optional, intent(in) :: dir @@ -120,6 +95,26 @@ module psb_i_sort_mod end subroutine psb_imsort end interface psb_msort + + interface psb_bsrch + function psb_ibsrch(key,n,v) result(ipos) + import + integer(psb_ipk_) :: ipos, n + integer(psb_ipk_) :: key + integer(psb_ipk_) :: v(:) + end function psb_ibsrch + end interface psb_bsrch + + interface psb_ssrch + function psb_issrch(key,n,v) result(ipos) + import + implicit none + integer(psb_ipk_) :: ipos, n + integer(psb_ipk_) :: key + integer(psb_ipk_) :: v(:) + end function psb_issrch + end interface psb_ssrch + interface subroutine psi_i_msort_up(n,k,l,iret) import @@ -527,7 +522,8 @@ contains implicit none class(psb_i_idx_heap), intent(inout) :: heap - integer(psb_ipk_), intent(out) :: index,info + integer(psb_ipk_), intent(out) :: index + integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: key diff --git a/base/modules/aux/psb_ip_reord_mod.f90 b/base/modules/aux/psb_ip_reord_mod.f90 index 325bc1beb..7a756bf45 100644 --- a/base/modules/aux/psb_ip_reord_mod.f90 +++ b/base/modules/aux/psb_ip_reord_mod.f90 @@ -37,634 +37,10 @@ ! ! module psb_ip_reord_mod - use psb_const_mod - - interface psb_ip_reord - module procedure psb_ip_reord_i1,& - & psb_ip_reord_s1, psb_ip_reord_d1,& - & psb_ip_reord_c1, psb_ip_reord_z1,& - & psb_ip_reord_i1i1,& - & psb_ip_reord_s1i1, psb_ip_reord_d1i1,& - & psb_ip_reord_c1i1, psb_ip_reord_z1i1,& - & psb_ip_reord_s1i2, psb_ip_reord_d1i2,& - & psb_ip_reord_c1i2, psb_ip_reord_z1i2,& - & psb_ip_reord_s1i3, psb_ip_reord_d1i3,& - & psb_ip_reord_c1i3, psb_ip_reord_z1i3 - - - end interface - -contains - - subroutine psb_ip_reord_i1(n,x,iaux) - integer(psb_ipk_), intent(in) :: n - integer(psb_ipk_) :: iaux(0:*) - integer(psb_ipk_) :: x(*) - - integer(psb_ipk_) :: lswap, lp, k - integer(psb_ipk_) :: swap - - lp = iaux(0) - k = 1 - do - if ((lp == 0).or.(k>n)) exit - do - if (lp >= k) exit - lp = iaux(lp) - end do - swap = x(lp) - x(lp) = x(k) - x(k) = swap - lswap = iaux(lp) - iaux(lp) = iaux(k) - iaux(k) = lp - lp = lswap - k = k + 1 - enddo - return - end subroutine psb_ip_reord_i1 - - - subroutine psb_ip_reord_s1(n,x,iaux) - integer(psb_ipk_), intent(in) :: n - integer(psb_ipk_) :: iaux(0:*) - real(psb_spk_) :: x(*) - - integer(psb_ipk_) :: lswap, lp, k - real(psb_spk_) :: swap - - lp = iaux(0) - k = 1 - do - if ((lp == 0).or.(k>n)) exit - do - if (lp >= k) exit - lp = iaux(lp) - end do - swap = x(lp) - x(lp) = x(k) - x(k) = swap - lswap = iaux(lp) - iaux(lp) = iaux(k) - iaux(k) = lp - lp = lswap - k = k + 1 - enddo - return - end subroutine psb_ip_reord_s1 - - subroutine psb_ip_reord_d1(n,x,iaux) - integer(psb_ipk_), intent(in) :: n - integer(psb_ipk_) :: iaux(0:*) - real(psb_dpk_) :: x(*) - - integer(psb_ipk_) :: lswap, lp, k - real(psb_dpk_) :: swap - - lp = iaux(0) - k = 1 - do - if ((lp == 0).or.(k>n)) exit - do - if (lp >= k) exit - lp = iaux(lp) - end do - swap = x(lp) - x(lp) = x(k) - x(k) = swap - lswap = iaux(lp) - iaux(lp) = iaux(k) - iaux(k) = lp - lp = lswap - k = k + 1 - enddo - return - end subroutine psb_ip_reord_d1 - - - - subroutine psb_ip_reord_c1(n,x,iaux) - integer(psb_ipk_), intent(in) :: n - integer(psb_ipk_) :: iaux(0:*) - complex(psb_spk_) :: x(*) - - integer(psb_ipk_) :: lswap, lp, k - complex(psb_spk_) :: swap - - lp = iaux(0) - k = 1 - do - if ((lp == 0).or.(k>n)) exit - do - if (lp >= k) exit - lp = iaux(lp) - end do - swap = x(lp) - x(lp) = x(k) - x(k) = swap - lswap = iaux(lp) - iaux(lp) = iaux(k) - iaux(k) = lp - lp = lswap - k = k + 1 - enddo - return - end subroutine psb_ip_reord_c1 - - subroutine psb_ip_reord_z1(n,x,iaux) - integer(psb_ipk_), intent(in) :: n - integer(psb_ipk_) :: iaux(0:*) - complex(psb_dpk_) :: x(*) - - integer(psb_ipk_) :: lswap, lp, k - complex(psb_dpk_) :: swap - - lp = iaux(0) - k = 1 - do - if ((lp == 0).or.(k>n)) exit - do - if (lp >= k) exit - lp = iaux(lp) - end do - swap = x(lp) - x(lp) = x(k) - x(k) = swap - lswap = iaux(lp) - iaux(lp) = iaux(k) - iaux(k) = lp - lp = lswap - k = k + 1 - enddo - return - end subroutine psb_ip_reord_z1 - - - subroutine psb_ip_reord_i1i1(n,x,indx,iaux) - integer(psb_ipk_), intent(in) :: n - integer(psb_ipk_) :: iaux(0:*) - integer(psb_ipk_) :: x(*) - integer(psb_ipk_) :: indx(*) - - integer(psb_ipk_) :: lswap, lp, k, ixswap - integer(psb_ipk_) :: swap - - lp = iaux(0) - k = 1 - do - if ((lp == 0).or.(k>n)) exit - do - if (lp >= k) exit - lp = iaux(lp) - end do - swap = x(lp) - x(lp) = x(k) - x(k) = swap - ixswap = indx(lp) - indx(lp) = indx(k) - indx(k) = ixswap - lswap = iaux(lp) - iaux(lp) = iaux(k) - iaux(k) = lp - lp = lswap - k = k + 1 - enddo - return - end subroutine psb_ip_reord_i1i1 - - - subroutine psb_ip_reord_s1i1(n,x,indx,iaux) - integer(psb_ipk_), intent(in) :: n - integer(psb_ipk_) :: iaux(0:*) - real(psb_spk_) :: x(*) - integer(psb_ipk_) :: indx(*) - - - integer(psb_ipk_) :: lswap, lp, k, ixswap - real(psb_spk_) :: swap - - lp = iaux(0) - k = 1 - do - if ((lp == 0).or.(k>n)) exit - do - if (lp >= k) exit - lp = iaux(lp) - end do - swap = x(lp) - x(lp) = x(k) - x(k) = swap - ixswap = indx(lp) - indx(lp) = indx(k) - indx(k) = ixswap - lswap = iaux(lp) - iaux(lp) = iaux(k) - iaux(k) = lp - lp = lswap - k = k + 1 - enddo - return - end subroutine psb_ip_reord_s1i1 - - subroutine psb_ip_reord_d1i1(n,x,indx,iaux) - integer(psb_ipk_), intent(in) :: n - integer(psb_ipk_) :: iaux(0:*) - real(psb_dpk_) :: x(*) - integer(psb_ipk_) :: indx(*) - - integer(psb_ipk_) :: lswap, lp, k, ixswap - real(psb_dpk_) :: swap - - lp = iaux(0) - k = 1 - do - if ((lp == 0).or.(k>n)) exit - do - if (lp >= k) exit - lp = iaux(lp) - end do - swap = x(lp) - x(lp) = x(k) - x(k) = swap - ixswap = indx(lp) - indx(lp) = indx(k) - indx(k) = ixswap - lswap = iaux(lp) - iaux(lp) = iaux(k) - iaux(k) = lp - lp = lswap - k = k + 1 - enddo - return - end subroutine psb_ip_reord_d1i1 - - - - subroutine psb_ip_reord_c1i1(n,x,indx,iaux) - integer(psb_ipk_), intent(in) :: n - integer(psb_ipk_) :: iaux(0:*) - complex(psb_spk_) :: x(*) - integer(psb_ipk_) :: indx(*) - - integer(psb_ipk_) :: lswap, lp, k, ixswap - complex(psb_spk_) :: swap - - lp = iaux(0) - k = 1 - do - if ((lp == 0).or.(k>n)) exit - do - if (lp >= k) exit - lp = iaux(lp) - end do - swap = x(lp) - x(lp) = x(k) - x(k) = swap - ixswap = indx(lp) - indx(lp) = indx(k) - indx(k) = ixswap - lswap = iaux(lp) - iaux(lp) = iaux(k) - iaux(k) = lp - lp = lswap - k = k + 1 - enddo - return - end subroutine psb_ip_reord_c1i1 - - subroutine psb_ip_reord_z1i1(n,x,indx,iaux) - integer(psb_ipk_), intent(in) :: n - integer(psb_ipk_) :: iaux(0:*) - complex(psb_dpk_) :: x(*) - integer(psb_ipk_) :: indx(*) - - integer(psb_ipk_) :: lswap, lp, k, ixswap - complex(psb_dpk_) :: swap - - lp = iaux(0) - k = 1 - do - if ((lp == 0).or.(k>n)) exit - do - if (lp >= k) exit - lp = iaux(lp) - end do - swap = x(lp) - x(lp) = x(k) - x(k) = swap - ixswap = indx(lp) - indx(lp) = indx(k) - indx(k) = ixswap - lswap = iaux(lp) - iaux(lp) = iaux(k) - iaux(k) = lp - lp = lswap - k = k + 1 - enddo - return - end subroutine psb_ip_reord_z1i1 - - - subroutine psb_ip_reord_s1i2(n,x,i1,i2,iaux) - integer(psb_ipk_), intent(in) :: n - integer(psb_ipk_) :: iaux(0:*) - real(psb_spk_) :: x(*) - integer(psb_ipk_) :: i1(*), i2(*) - - - integer(psb_ipk_) :: lswap, lp, k, isw1, isw2 - real(psb_spk_) :: swap - - lp = iaux(0) - k = 1 - do - if ((lp == 0).or.(k>n)) exit - do - if (lp >= k) exit - lp = iaux(lp) - end do - swap = x(lp) - x(lp) = x(k) - x(k) = swap - isw1 = i1(lp) - i1(lp) = i1(k) - i1(k) = isw1 - isw2 = i2(lp) - i2(lp) = i2(k) - i2(k) = isw2 - lswap = iaux(lp) - iaux(lp) = iaux(k) - iaux(k) = lp - lp = lswap - k = k + 1 - enddo - return - end subroutine psb_ip_reord_s1i2 - - subroutine psb_ip_reord_d1i2(n,x,i1,i2,iaux) - integer(psb_ipk_), intent(in) :: n - integer(psb_ipk_) :: iaux(0:*) - real(psb_dpk_) :: x(*) - integer(psb_ipk_) :: i1(*), i2(*) - - - integer(psb_ipk_) :: lswap, lp, k, isw1, isw2 - real(psb_dpk_) :: swap - - lp = iaux(0) - k = 1 - do - if ((lp == 0).or.(k>n)) exit - do - if (lp >= k) exit - lp = iaux(lp) - end do - swap = x(lp) - x(lp) = x(k) - x(k) = swap - isw1 = i1(lp) - i1(lp) = i1(k) - i1(k) = isw1 - isw2 = i2(lp) - i2(lp) = i2(k) - i2(k) = isw2 - lswap = iaux(lp) - iaux(lp) = iaux(k) - iaux(k) = lp - lp = lswap - k = k + 1 - enddo - return - end subroutine psb_ip_reord_d1i2 - - subroutine psb_ip_reord_c1i2(n,x,i1,i2,iaux) - integer(psb_ipk_), intent(in) :: n - integer(psb_ipk_) :: iaux(0:*) - complex(psb_spk_) :: x(*) - integer(psb_ipk_) :: i1(*), i2(*) - - - integer(psb_ipk_) :: lswap, lp, k, isw1, isw2 - complex(psb_spk_) :: swap - - lp = iaux(0) - k = 1 - do - if ((lp == 0).or.(k>n)) exit - do - if (lp >= k) exit - lp = iaux(lp) - end do - swap = x(lp) - x(lp) = x(k) - x(k) = swap - isw1 = i1(lp) - i1(lp) = i1(k) - i1(k) = isw1 - isw2 = i2(lp) - i2(lp) = i2(k) - i2(k) = isw2 - lswap = iaux(lp) - iaux(lp) = iaux(k) - iaux(k) = lp - lp = lswap - k = k + 1 - enddo - return - end subroutine psb_ip_reord_c1i2 - - subroutine psb_ip_reord_z1i2(n,x,i1,i2,iaux) - integer(psb_ipk_), intent(in) :: n - integer(psb_ipk_) :: iaux(0:*) - complex(psb_dpk_) :: x(*) - integer(psb_ipk_) :: i1(*), i2(*) - - - integer(psb_ipk_) :: lswap, lp, k, isw1, isw2 - complex(psb_dpk_) :: swap - - lp = iaux(0) - k = 1 - do - if ((lp == 0).or.(k>n)) exit - do - if (lp >= k) exit - lp = iaux(lp) - end do - swap = x(lp) - x(lp) = x(k) - x(k) = swap - isw1 = i1(lp) - i1(lp) = i1(k) - i1(k) = isw1 - isw2 = i2(lp) - i2(lp) = i2(k) - i2(k) = isw2 - lswap = iaux(lp) - iaux(lp) = iaux(k) - iaux(k) = lp - lp = lswap - k = k + 1 - enddo - return - end subroutine psb_ip_reord_z1i2 - - - subroutine psb_ip_reord_s1i3(n,x,i1,i2,i3,iaux) - integer(psb_ipk_), intent(in) :: n - integer(psb_ipk_) :: iaux(0:*) - real(psb_spk_) :: x(*) - integer(psb_ipk_) :: i1(*), i2(*), i3(*) - - - integer(psb_ipk_) :: lswap, lp, k, isw1, isw2, isw3 - real(psb_spk_) :: swap - - lp = iaux(0) - k = 1 - do - if ((lp == 0).or.(k>n)) exit - do - if (lp >= k) exit - lp = iaux(lp) - end do - swap = x(lp) - x(lp) = x(k) - x(k) = swap - isw1 = i1(lp) - i1(lp) = i1(k) - i1(k) = isw1 - isw2 = i2(lp) - i2(lp) = i2(k) - i2(k) = isw2 - isw3 = i3(lp) - i3(lp) = i3(k) - i3(k) = isw3 - lswap = iaux(lp) - iaux(lp) = iaux(k) - iaux(k) = lp - lp = lswap - k = k + 1 - enddo - return - end subroutine psb_ip_reord_s1i3 - - subroutine psb_ip_reord_d1i3(n,x,i1,i2,i3,iaux) - integer(psb_ipk_), intent(in) :: n - integer(psb_ipk_) :: iaux(0:*) - real(psb_dpk_) :: x(*) - integer(psb_ipk_) :: i1(*), i2(*),i3(*) - - - integer(psb_ipk_) :: lswap, lp, k, isw1, isw2,isw3 - real(psb_dpk_) :: swap - - lp = iaux(0) - k = 1 - do - if ((lp == 0).or.(k>n)) exit - do - if (lp >= k) exit - lp = iaux(lp) - end do - swap = x(lp) - x(lp) = x(k) - x(k) = swap - isw1 = i1(lp) - i1(lp) = i1(k) - i1(k) = isw1 - isw2 = i2(lp) - i2(lp) = i2(k) - i2(k) = isw2 - isw3 = i3(lp) - i3(lp) = i3(k) - i3(k) = isw3 - lswap = iaux(lp) - iaux(lp) = iaux(k) - iaux(k) = lp - lp = lswap - k = k + 1 - enddo - return - end subroutine psb_ip_reord_d1i3 - - subroutine psb_ip_reord_c1i3(n,x,i1,i2,i3,iaux) - integer(psb_ipk_), intent(in) :: n - integer(psb_ipk_) :: iaux(0:*) - complex(psb_spk_) :: x(*) - integer(psb_ipk_) :: i1(*), i2(*), i3(*) - - - integer(psb_ipk_) :: lswap, lp, k, isw1, isw2, isw3 - complex(psb_spk_) :: swap - - lp = iaux(0) - k = 1 - do - if ((lp == 0).or.(k>n)) exit - do - if (lp >= k) exit - lp = iaux(lp) - end do - swap = x(lp) - x(lp) = x(k) - x(k) = swap - isw1 = i1(lp) - i1(lp) = i1(k) - i1(k) = isw1 - isw2 = i2(lp) - i2(lp) = i2(k) - i2(k) = isw2 - isw3 = i3(lp) - i3(lp) = i3(k) - i3(k) = isw3 - lswap = iaux(lp) - iaux(lp) = iaux(k) - iaux(k) = lp - lp = lswap - k = k + 1 - enddo - return - end subroutine psb_ip_reord_c1i3 - - subroutine psb_ip_reord_z1i3(n,x,i1,i2,i3,iaux) - integer(psb_ipk_), intent(in) :: n - integer(psb_ipk_) :: iaux(0:*) - complex(psb_dpk_) :: x(*) - integer(psb_ipk_) :: i1(*), i2(*), i3(*) - - - integer(psb_ipk_) :: lswap, lp, k, isw1, isw2, isw3 - complex(psb_dpk_) :: swap - - lp = iaux(0) - k = 1 - do - if ((lp == 0).or.(k>n)) exit - do - if (lp >= k) exit - lp = iaux(lp) - end do - swap = x(lp) - x(lp) = x(k) - x(k) = swap - isw1 = i1(lp) - i1(lp) = i1(k) - i1(k) = isw1 - isw2 = i2(lp) - i2(lp) = i2(k) - i2(k) = isw2 - isw3 = i3(lp) - i3(lp) = i3(k) - i3(k) = isw3 - lswap = iaux(lp) - iaux(lp) = iaux(k) - iaux(k) = lp - lp = lswap - k = k + 1 - enddo - return - end subroutine psb_ip_reord_z1i3 - - + use psb_i_ip_reord_mod + use psb_l_ip_reord_mod + use psb_s_ip_reord_mod + use psb_c_ip_reord_mod + use psb_d_ip_reord_mod + use psb_z_ip_reord_mod end module psb_ip_reord_mod diff --git a/base/modules/aux/psb_l_ip_reord_mod.f90 b/base/modules/aux/psb_l_ip_reord_mod.f90 new file mode 100644 index 000000000..b46010747 --- /dev/null +++ b/base/modules/aux/psb_l_ip_reord_mod.f90 @@ -0,0 +1,187 @@ +! +! Parallel Sparse BLAS version 3.5 +! (C) Copyright 2006-2018 +! Salvatore Filippone +! Alfredo Buttari +! +! Redistribution and use in source and binary forms, with or without +! modification, are permitted provided that the following conditions +! are met: +! 1. Redistributions of source code must retain the above copyright +! notice, this list of conditions and the following disclaimer. +! 2. Redistributions in binary form must reproduce the above copyright +! notice, this list of conditions, and the following disclaimer in the +! documentation and/or other materials provided with the distribution. +! 3. The name of the PSBLAS group or the names of its contributors may +! not be used to endorse or promote products derived from this +! software without specific written permission. +! +! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +! ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED +! TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR +! PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE PSBLAS GROUP OR ITS CONTRIBUTORS +! BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +! CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +! SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +! INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +! CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +! ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +! POSSIBILITY OF SUCH DAMAGE. +! +! +! +! Reorder (an) input vector(s) based on a list sort output. +! Based on: D. E. Knuth: The Art of Computer Programming +! vol. 3: Sorting and Searching, Addison Wesley, 1973 +! ex. 5.2.12 +! +! +module psb_l_ip_reord_mod + use psb_const_mod + + interface psb_ip_reord + module procedure psb_ip_reord_l1,& + & psb_ip_reord_l1i1, psb_ip_reord_l1i2,& + & psb_ip_reord_l1i3 + + + end interface + +contains + + subroutine psb_ip_reord_l1(n,x,iaux) + integer(psb_ipk_), intent(in) :: n + integer(psb_lpk_) :: iaux(0:*) + integer(psb_lpk_) :: x(*) + + integer(psb_lpk_) :: lswap, lp, k + integer(psb_lpk_) :: swap + + lp = iaux(0) + k = 1 + do + if ((lp == 0).or.(k>n)) exit + do + if (lp >= k) exit + lp = iaux(lp) + end do + swap = x(lp) + x(lp) = x(k) + x(k) = swap + lswap = iaux(lp) + iaux(lp) = iaux(k) + iaux(k) = lp + lp = lswap + k = k + 1 + enddo + return + end subroutine psb_ip_reord_l1 + + subroutine psb_ip_reord_l1i1(n,x,indx,iaux) + integer(psb_ipk_), intent(in) :: n + integer(psb_lpk_) :: iaux(0:*) + integer(psb_lpk_) :: x(*) + integer(psb_lpk_) :: indx(*) + + integer(psb_lpk_) :: lswap, lp, k, ixswap + integer(psb_lpk_) :: swap + + lp = iaux(0) + k = 1 + do + if ((lp == 0).or.(k>n)) exit + do + if (lp >= k) exit + lp = iaux(lp) + end do + swap = x(lp) + x(lp) = x(k) + x(k) = swap + ixswap = indx(lp) + indx(lp) = indx(k) + indx(k) = ixswap + lswap = iaux(lp) + iaux(lp) = iaux(k) + iaux(k) = lp + lp = lswap + k = k + 1 + enddo + return + end subroutine psb_ip_reord_l1i1 + + subroutine psb_ip_reord_l1i2(n,x,i1,i2,iaux) + integer(psb_ipk_), intent(in) :: n + integer(psb_lpk_) :: iaux(0:*) + integer(psb_lpk_) :: x(*) + integer(psb_lpk_) :: i1(*), i2(*) + + + integer(psb_lpk_) :: lswap, lp, k, isw1, isw2 + integer(psb_lpk_) :: swap + + lp = iaux(0) + k = 1 + do + if ((lp == 0).or.(k>n)) exit + do + if (lp >= k) exit + lp = iaux(lp) + end do + swap = x(lp) + x(lp) = x(k) + x(k) = swap + isw1 = i1(lp) + i1(lp) = i1(k) + i1(k) = isw1 + isw2 = i2(lp) + i2(lp) = i2(k) + i2(k) = isw2 + lswap = iaux(lp) + iaux(lp) = iaux(k) + iaux(k) = lp + lp = lswap + k = k + 1 + enddo + return + end subroutine psb_ip_reord_l1i2 + + subroutine psb_ip_reord_l1i3(n,x,i1,i2,i3,iaux) + integer(psb_ipk_), intent(in) :: n + integer(psb_lpk_) :: iaux(0:*) + integer(psb_lpk_) :: x(*) + integer(psb_lpk_) :: i1(*), i2(*), i3(*) + + + integer(psb_lpk_) :: lswap, lp, k, isw1, isw2, isw3 + integer(psb_lpk_) :: swap + + lp = iaux(0) + k = 1 + do + if ((lp == 0).or.(k>n)) exit + do + if (lp >= k) exit + lp = iaux(lp) + end do + swap = x(lp) + x(lp) = x(k) + x(k) = swap + isw1 = i1(lp) + i1(lp) = i1(k) + i1(k) = isw1 + isw2 = i2(lp) + i2(lp) = i2(k) + i2(k) = isw2 + isw3 = i3(lp) + i3(lp) = i3(k) + i3(k) = isw3 + lswap = iaux(lp) + iaux(lp) = iaux(k) + iaux(k) = lp + lp = lswap + k = k + 1 + enddo + return + end subroutine psb_ip_reord_l1i3 + +end module psb_l_ip_reord_mod diff --git a/base/modules/aux/psb_l_sort_mod.f90 b/base/modules/aux/psb_l_sort_mod.f90 new file mode 100644 index 000000000..3f156ef68 --- /dev/null +++ b/base/modules/aux/psb_l_sort_mod.f90 @@ -0,0 +1,577 @@ +! +! Parallel Sparse BLAS version 3.5 +! (C) Copyright 2006-2018 +! Salvatore Filippone +! Alfredo Buttari +! +! Redistribution and use in source and binary forms, with or without +! modification, are permitted provided that the following conditions +! are met: +! 1. Redistributions of source code must retain the above copyright +! notice, this list of conditions and the following disclaimer. +! 2. Redistributions in binary form must reproduce the above copyright +! notice, this list of conditions, and the following disclaimer in the +! documentation and/or other materials provided with the distribution. +! 3. The name of the PSBLAS group or the names of its contributors may +! not be used to endorse or promote products derived from this +! software without specific written permission. +! +! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +! ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED +! TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR +! PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE PSBLAS GROUP OR ITS CONTRIBUTORS +! BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +! CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +! SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +! INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +! CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +! ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +! POSSIBILITY OF SUCH DAMAGE. +! +! +! +! Sorting routines +! References: +! D. Knuth +! The Art of Computer Programming, vol. 3 +! Addison-Wesley +! +! Aho, Hopcroft, Ullman +! Data Structures and Algorithms +! Addison-Wesley +! +module psb_l_sort_mod + use psb_const_mod + + interface psb_isaperm + logical function psb_lisaperm(n,eip) + import + integer(psb_lpk_), intent(in) :: n + integer(psb_lpk_), intent(in) :: eip(n) + end function psb_lisaperm + end interface psb_isaperm + + interface psb_msort_unique + subroutine psb_lmsort_u(x,nout,dir) + import + integer(psb_lpk_), intent(inout) :: x(:) + integer(psb_lpk_), intent(out) :: nout + integer(psb_ipk_), optional, intent(in) :: dir + end subroutine psb_lmsort_u + end interface psb_msort_unique + + type psb_l_heap + integer(psb_ipk_) :: last, dir + integer(psb_lpk_), allocatable :: keys(:) + contains + procedure, pass(heap) :: init => psb_l_init_heap + procedure, pass(heap) :: howmany => psb_l_howmany + procedure, pass(heap) :: insert => psb_l_insert_heap + procedure, pass(heap) :: get_first => psb_l_heap_get_first + procedure, pass(heap) :: dump => psb_l_dump_heap + procedure, pass(heap) :: free => psb_l_free_heap + end type psb_l_heap + + type psb_l_idx_heap + integer(psb_ipk_) :: last, dir + integer(psb_lpk_), allocatable :: keys(:) + integer(psb_lpk_), allocatable :: idxs(:) + contains + procedure, pass(heap) :: init => psb_l_idx_init_heap + procedure, pass(heap) :: howmany => psb_l_idx_howmany + procedure, pass(heap) :: insert => psb_l_idx_insert_heap + procedure, pass(heap) :: get_first => psb_l_idx_heap_get_first + procedure, pass(heap) :: dump => psb_l_idx_dump_heap + procedure, pass(heap) :: free => psb_l_idx_free_heap + end type psb_l_idx_heap + + + interface psb_msort + subroutine psb_lmsort(x,ix,dir,flag) + import + integer(psb_lpk_), intent(inout) :: x(:) + integer(psb_ipk_), optional, intent(in) :: dir, flag + integer(psb_lpk_), optional, intent(inout) :: ix(:) + end subroutine psb_lmsort + end interface psb_msort + + + interface psb_bsrch + function psb_lbsrch(key,n,v) result(ipos) + import + integer(psb_ipk_) :: ipos, n + integer(psb_lpk_) :: key + integer(psb_lpk_) :: v(:) + end function psb_lbsrch + end interface psb_bsrch + + interface psb_ssrch + function psb_lssrch(key,n,v) result(ipos) + import + implicit none + integer(psb_ipk_) :: ipos, n + integer(psb_lpk_) :: key + integer(psb_lpk_) :: v(:) + end function psb_lssrch + end interface psb_ssrch + + interface + subroutine psi_l_msort_up(n,k,l,iret) + import + implicit none + integer(psb_ipk_) :: n, iret + integer(psb_lpk_) :: k(n) + integer(psb_lpk_) :: l(0:n+1) + end subroutine psi_l_msort_up + subroutine psi_l_msort_dw(n,k,l,iret) + import + implicit none + integer(psb_ipk_) :: n, iret + integer(psb_lpk_) :: k(n) + integer(psb_lpk_) :: l(0:n+1) + end subroutine psi_l_msort_dw + end interface + interface + subroutine psi_l_amsort_up(n,k,l,iret) + import + implicit none + integer(psb_ipk_) :: n, iret + integer(psb_lpk_) :: k(n) + integer(psb_lpk_) :: l(0:n+1) + end subroutine psi_l_amsort_up + subroutine psi_l_amsort_dw(n,k,l,iret) + import + implicit none + integer(psb_ipk_) :: n, iret + integer(psb_lpk_) :: k(n) + integer(psb_lpk_) :: l(0:n+1) + end subroutine psi_l_amsort_dw + end interface + + + interface psb_qsort + subroutine psb_lqsort(x,ix,dir,flag) + import + integer(psb_lpk_), intent(inout) :: x(:) + integer(psb_ipk_), optional, intent(in) :: dir, flag + integer(psb_lpk_), optional, intent(inout) :: ix(:) + end subroutine psb_lqsort + end interface psb_qsort + + interface psb_isort + subroutine psb_lisort(x,ix,dir,flag) + import + integer(psb_lpk_), intent(inout) :: x(:) + integer(psb_ipk_), optional, intent(in) :: dir, flag + integer(psb_lpk_), optional, intent(inout) :: ix(:) + end subroutine psb_lisort + end interface psb_isort + + + interface psb_hsort + subroutine psb_lhsort(x,ix,dir,flag) + import + integer(psb_lpk_), intent(inout) :: x(:) + integer(psb_ipk_), optional, intent(in) :: dir, flag + integer(psb_lpk_), optional, intent(inout) :: ix(:) + end subroutine psb_lhsort + end interface psb_hsort + + + interface + subroutine psi_l_insert_heap(key,last,heap,dir,info) + import + implicit none + + ! + ! Input: + ! key: the new value + ! last: pointer to the last occupied element in heap + ! heap: the heap + ! dir: sorting direction + + integer(psb_lpk_), intent(in) :: key + integer(psb_lpk_), intent(inout) :: heap(:) + integer(psb_ipk_), intent(in) :: dir + integer(psb_ipk_), intent(inout) :: last + integer(psb_ipk_), intent(out) :: info + end subroutine psi_l_insert_heap + end interface + + interface + subroutine psi_l_idx_insert_heap(key,index,last,heap,idxs,dir,info) + import + implicit none + + ! + ! Input: + ! key: the new value + ! last: pointer to the last occupied element in heap + ! heap: the heap + ! dir: sorting direction + + integer(psb_lpk_), intent(in) :: key + integer(psb_lpk_), intent(inout) :: heap(:) + integer(psb_lpk_), intent(in) :: index + integer(psb_ipk_), intent(in) :: dir + integer(psb_lpk_), intent(inout) :: idxs(:) + integer(psb_ipk_), intent(inout) :: last + integer(psb_ipk_), intent(out) :: info + end subroutine psi_l_idx_insert_heap + end interface + + + interface + subroutine psi_l_heap_get_first(key,last,heap,dir,info) + import + implicit none + integer(psb_lpk_), intent(inout) :: key + integer(psb_ipk_), intent(inout) :: last + integer(psb_ipk_), intent(in) :: dir + integer(psb_lpk_), intent(inout) :: heap(:) + integer(psb_ipk_), intent(out) :: info + end subroutine psi_l_heap_get_first + end interface + + interface + subroutine psi_l_idx_heap_get_first(key,index,last,heap,idxs,dir,info) + import + integer(psb_lpk_), intent(inout) :: key + integer(psb_lpk_), intent(out) :: index + integer(psb_lpk_), intent(inout) :: heap(:) + integer(psb_ipk_), intent(in) :: dir + integer(psb_ipk_), intent(inout) :: last + integer(psb_lpk_), intent(inout) :: idxs(:) + integer(psb_ipk_), intent(out) :: info + end subroutine psi_l_idx_heap_get_first + end interface + + interface + subroutine psi_lisrx_up(n,x,ix) + import + integer(psb_lpk_), intent(inout) :: x(:) + integer(psb_lpk_), intent(inout) :: ix(:) + integer(psb_lpk_), intent(in) :: n + end subroutine psi_lisrx_up + subroutine psi_lisrx_dw(n,x,ix) + import + integer(psb_lpk_), intent(inout) :: x(:) + integer(psb_lpk_), intent(inout) :: ix(:) + integer(psb_lpk_), intent(in) :: n + end subroutine psi_lisrx_dw + subroutine psi_lisr_up(n,x) + import + integer(psb_lpk_), intent(inout) :: x(:) + integer(psb_lpk_), intent(in) :: n + end subroutine psi_lisr_up + subroutine psi_lisr_dw(n,x) + import + integer(psb_lpk_), intent(inout) :: x(:) + integer(psb_lpk_), intent(in) :: n + end subroutine psi_lisr_dw + subroutine psi_laisrx_up(n,x,ix) + import + integer(psb_lpk_), intent(inout) :: x(:) + integer(psb_lpk_), intent(inout) :: ix(:) + integer(psb_lpk_), intent(in) :: n + end subroutine psi_laisrx_up + subroutine psi_laisrx_dw(n,x,ix) + import + integer(psb_lpk_), intent(inout) :: x(:) + integer(psb_lpk_), intent(inout) :: ix(:) + integer(psb_lpk_), intent(in) :: n + end subroutine psi_laisrx_dw + subroutine psi_laisr_up(n,x) + import + integer(psb_lpk_), intent(inout) :: x(:) + integer(psb_lpk_), intent(in) :: n + end subroutine psi_laisr_up + subroutine psi_laisr_dw(n,x) + import + integer(psb_lpk_), intent(inout) :: x(:) + integer(psb_lpk_), intent(in) :: n + end subroutine psi_laisr_dw + end interface + + interface + subroutine psi_lqsrx_up(n,x,ix) + import + integer(psb_lpk_), intent(inout) :: x(:) + integer(psb_lpk_), intent(inout) :: ix(:) + integer(psb_lpk_), intent(in) :: n + end subroutine psi_lqsrx_up + subroutine psi_lqsrx_dw(n,x,ix) + import + integer(psb_lpk_), intent(inout) :: x(:) + integer(psb_lpk_), intent(inout) :: ix(:) + integer(psb_lpk_), intent(in) :: n + end subroutine psi_lqsrx_dw + subroutine psi_lqsr_up(n,x) + import + integer(psb_lpk_), intent(inout) :: x(:) + integer(psb_lpk_), intent(in) :: n + end subroutine psi_lqsr_up + subroutine psi_lqsr_dw(n,x) + import + integer(psb_lpk_), intent(inout) :: x(:) + integer(psb_lpk_), intent(in) :: n + end subroutine psi_lqsr_dw + subroutine psi_laqsrx_up(n,x,ix) + import + integer(psb_lpk_), intent(inout) :: x(:) + integer(psb_lpk_), intent(inout) :: ix(:) + integer(psb_lpk_), intent(in) :: n + end subroutine psi_laqsrx_up + subroutine psi_laqsrx_dw(n,x,ix) + import + integer(psb_lpk_), intent(inout) :: x(:) + integer(psb_lpk_), intent(inout) :: ix(:) + integer(psb_lpk_), intent(in) :: n + end subroutine psi_laqsrx_dw + subroutine psi_laqsr_up(n,x) + import + integer(psb_lpk_), intent(inout) :: x(:) + integer(psb_lpk_), intent(in) :: n + end subroutine psi_laqsr_up + subroutine psi_laqsr_dw(n,x) + import + integer(psb_lpk_), intent(inout) :: x(:) + integer(psb_lpk_), intent(in) :: n + end subroutine psi_laqsr_dw + end interface + +contains + + subroutine psb_l_init_heap(heap,info,dir) + use psb_realloc_mod, only : psb_ensure_size + implicit none + class(psb_l_heap), intent(inout) :: heap + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_), intent(in), optional :: dir + + info = psb_success_ + heap%last=0 + if (present(dir)) then + heap%dir = dir + else + heap%dir = psb_sort_up_ + endif + select case(heap%dir) + case (psb_sort_up_,psb_sort_down_,psb_asort_up_,psb_asort_down_) + ! ok, do nothing + case default + write(psb_err_unit,*) 'Invalid direction, defaulting to psb_sort_up_' + heap%dir = psb_sort_up_ + end select + call psb_ensure_size(psb_heap_resize,heap%keys,info) + + return + end subroutine psb_l_init_heap + + + function psb_l_howmany(heap) result(res) + implicit none + class(psb_l_heap), intent(in) :: heap + integer(psb_ipk_) :: res + res = heap%last + end function psb_l_howmany + + subroutine psb_l_insert_heap(key,heap,info) + use psb_realloc_mod, only : psb_ensure_size + implicit none + + integer(psb_lpk_), intent(in) :: key + class(psb_l_heap), intent(inout) :: heap + integer(psb_ipk_), intent(out) :: info + + info = psb_success_ + if (heap%last < 0) then + write(psb_err_unit,*) 'Invalid last in heap ',heap%last + info = heap%last + return + endif + + call psb_ensure_size(heap%last+1,heap%keys,info,addsz=psb_heap_resize) + if (info /= psb_success_) then + write(psb_err_unit,*) 'Memory allocation failure in heap_insert' + info = -5 + return + end if + call psi_l_insert_heap(key,& + & heap%last,heap%keys,heap%dir,info) + + return + end subroutine psb_l_insert_heap + + subroutine psb_l_heap_get_first(key,heap,info) + implicit none + + class(psb_l_heap), intent(inout) :: heap + integer(psb_ipk_), intent(out) :: info + integer(psb_lpk_), intent(out) :: key + + + info = psb_success_ + + call psi_l_heap_get_first(key,& + & heap%last,heap%keys,heap%dir,info) + + return + end subroutine psb_l_heap_get_first + + subroutine psb_l_dump_heap(iout,heap,info) + + implicit none + class(psb_l_heap), intent(in) :: heap + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_), intent(in) :: iout + + info = psb_success_ + if (iout < 0) then + write(psb_err_unit,*) 'Invalid file ' + info =-1 + return + end if + + write(iout,*) 'Heap direction ',heap%dir + write(iout,*) 'Heap size ',heap%last + if ((heap%last > 0).and.((.not.allocated(heap%keys)).or.& + & (size(heap%keys) 0).and.((.not.allocated(heap%keys)).or.& + & (size(heap%keys) 0).and.((.not.allocated(heap%idxs)).or.& + & (size(heap%idxs)n)) exit + do + if (lp >= k) exit + lp = iaux(lp) + end do + swap = x(lp) + x(lp) = x(k) + x(k) = swap + lswap = iaux(lp) + iaux(lp) = iaux(k) + iaux(k) = lp + lp = lswap + k = k + 1 + enddo + return + end subroutine psb_ip_reord_s1 + + subroutine psb_ip_reord_s1i1(n,x,indx,iaux) + integer(psb_ipk_), intent(in) :: n + integer(psb_ipk_) :: iaux(0:*) + real(psb_spk_) :: x(*) + integer(psb_ipk_) :: indx(*) + + integer(psb_ipk_) :: lswap, lp, k, ixswap + real(psb_spk_) :: swap + + lp = iaux(0) + k = 1 + do + if ((lp == 0).or.(k>n)) exit + do + if (lp >= k) exit + lp = iaux(lp) + end do + swap = x(lp) + x(lp) = x(k) + x(k) = swap + ixswap = indx(lp) + indx(lp) = indx(k) + indx(k) = ixswap + lswap = iaux(lp) + iaux(lp) = iaux(k) + iaux(k) = lp + lp = lswap + k = k + 1 + enddo + return + end subroutine psb_ip_reord_s1i1 + + subroutine psb_ip_reord_s1i2(n,x,i1,i2,iaux) + integer(psb_ipk_), intent(in) :: n + integer(psb_ipk_) :: iaux(0:*) + real(psb_spk_) :: x(*) + integer(psb_ipk_) :: i1(*), i2(*) + + + integer(psb_ipk_) :: lswap, lp, k, isw1, isw2 + real(psb_spk_) :: swap + + lp = iaux(0) + k = 1 + do + if ((lp == 0).or.(k>n)) exit + do + if (lp >= k) exit + lp = iaux(lp) + end do + swap = x(lp) + x(lp) = x(k) + x(k) = swap + isw1 = i1(lp) + i1(lp) = i1(k) + i1(k) = isw1 + isw2 = i2(lp) + i2(lp) = i2(k) + i2(k) = isw2 + lswap = iaux(lp) + iaux(lp) = iaux(k) + iaux(k) = lp + lp = lswap + k = k + 1 + enddo + return + end subroutine psb_ip_reord_s1i2 + + subroutine psb_ip_reord_s1i3(n,x,i1,i2,i3,iaux) + integer(psb_ipk_), intent(in) :: n + integer(psb_ipk_) :: iaux(0:*) + real(psb_spk_) :: x(*) + integer(psb_ipk_) :: i1(*), i2(*), i3(*) + + + integer(psb_ipk_) :: lswap, lp, k, isw1, isw2, isw3 + real(psb_spk_) :: swap + + lp = iaux(0) + k = 1 + do + if ((lp == 0).or.(k>n)) exit + do + if (lp >= k) exit + lp = iaux(lp) + end do + swap = x(lp) + x(lp) = x(k) + x(k) = swap + isw1 = i1(lp) + i1(lp) = i1(k) + i1(k) = isw1 + isw2 = i2(lp) + i2(lp) = i2(k) + i2(k) = isw2 + isw3 = i3(lp) + i3(lp) = i3(k) + i3(k) = isw3 + lswap = iaux(lp) + iaux(lp) = iaux(k) + iaux(k) = lp + lp = lswap + k = k + 1 + enddo + return + end subroutine psb_ip_reord_s1i3 + +end module psb_s_ip_reord_mod diff --git a/base/modules/aux/psb_s_sort_mod.f90 b/base/modules/aux/psb_s_sort_mod.f90 index 4272dfe98..bd67ac068 100644 --- a/base/modules/aux/psb_s_sort_mod.f90 +++ b/base/modules/aux/psb_s_sort_mod.f90 @@ -44,6 +44,15 @@ module psb_s_sort_mod use psb_const_mod + interface psb_msort_unique + subroutine psb_smsort_u(x,nout,dir) + import + real(psb_spk_), intent(inout) :: x(:) + integer(psb_ipk_), intent(out) :: nout + integer(psb_ipk_), optional, intent(in) :: dir + end subroutine psb_smsort_u + end interface psb_msort_unique + type psb_s_heap integer(psb_ipk_) :: last, dir real(psb_spk_), allocatable :: keys(:) @@ -79,6 +88,26 @@ module psb_s_sort_mod end subroutine psb_smsort end interface psb_msort + + interface psb_bsrch + function psb_sbsrch(key,n,v) result(ipos) + import + integer(psb_ipk_) :: ipos, n + real(psb_spk_) :: key + real(psb_spk_) :: v(:) + end function psb_sbsrch + end interface psb_bsrch + + interface psb_ssrch + function psb_sssrch(key,n,v) result(ipos) + import + implicit none + integer(psb_ipk_) :: ipos, n + real(psb_spk_) :: key + real(psb_spk_) :: v(:) + end function psb_sssrch + end interface psb_ssrch + interface subroutine psi_s_msort_up(n,k,l,iret) import @@ -486,7 +515,8 @@ contains implicit none class(psb_s_idx_heap), intent(inout) :: heap - integer(psb_ipk_), intent(out) :: index,info + integer(psb_ipk_), intent(out) :: index + integer(psb_ipk_), intent(out) :: info real(psb_spk_), intent(out) :: key diff --git a/base/modules/aux/psb_sort_mod.f90 b/base/modules/aux/psb_sort_mod.f90 index b0ffd9864..99f1a1e4e 100644 --- a/base/modules/aux/psb_sort_mod.f90 +++ b/base/modules/aux/psb_sort_mod.f90 @@ -46,6 +46,7 @@ module psb_sort_mod use psb_const_mod use psb_ip_reord_mod use psb_i_sort_mod + use psb_l_sort_mod use psb_s_sort_mod use psb_c_sort_mod use psb_d_sort_mod diff --git a/base/modules/aux/psb_z_ip_reord_mod.f90 b/base/modules/aux/psb_z_ip_reord_mod.f90 new file mode 100644 index 000000000..5d2d36997 --- /dev/null +++ b/base/modules/aux/psb_z_ip_reord_mod.f90 @@ -0,0 +1,187 @@ +! +! Parallel Sparse BLAS version 3.5 +! (C) Copyright 2006-2018 +! Salvatore Filippone +! Alfredo Buttari +! +! Redistribution and use in source and binary forms, with or without +! modification, are permitted provided that the following conditions +! are met: +! 1. Redistributions of source code must retain the above copyright +! notice, this list of conditions and the following disclaimer. +! 2. Redistributions in binary form must reproduce the above copyright +! notice, this list of conditions, and the following disclaimer in the +! documentation and/or other materials provided with the distribution. +! 3. The name of the PSBLAS group or the names of its contributors may +! not be used to endorse or promote products derived from this +! software without specific written permission. +! +! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +! ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED +! TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR +! PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE PSBLAS GROUP OR ITS CONTRIBUTORS +! BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +! CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +! SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +! INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +! CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +! ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +! POSSIBILITY OF SUCH DAMAGE. +! +! +! +! Reorder (an) input vector(s) based on a list sort output. +! Based on: D. E. Knuth: The Art of Computer Programming +! vol. 3: Sorting and Searching, Addison Wesley, 1973 +! ex. 5.2.12 +! +! +module psb_z_ip_reord_mod + use psb_const_mod + + interface psb_ip_reord + module procedure psb_ip_reord_z1,& + & psb_ip_reord_z1i1, psb_ip_reord_z1i2,& + & psb_ip_reord_z1i3 + + + end interface + +contains + + subroutine psb_ip_reord_z1(n,x,iaux) + integer(psb_ipk_), intent(in) :: n + integer(psb_ipk_) :: iaux(0:*) + complex(psb_dpk_) :: x(*) + + integer(psb_ipk_) :: lswap, lp, k + complex(psb_dpk_) :: swap + + lp = iaux(0) + k = 1 + do + if ((lp == 0).or.(k>n)) exit + do + if (lp >= k) exit + lp = iaux(lp) + end do + swap = x(lp) + x(lp) = x(k) + x(k) = swap + lswap = iaux(lp) + iaux(lp) = iaux(k) + iaux(k) = lp + lp = lswap + k = k + 1 + enddo + return + end subroutine psb_ip_reord_z1 + + subroutine psb_ip_reord_z1i1(n,x,indx,iaux) + integer(psb_ipk_), intent(in) :: n + integer(psb_ipk_) :: iaux(0:*) + complex(psb_dpk_) :: x(*) + integer(psb_ipk_) :: indx(*) + + integer(psb_ipk_) :: lswap, lp, k, ixswap + complex(psb_dpk_) :: swap + + lp = iaux(0) + k = 1 + do + if ((lp == 0).or.(k>n)) exit + do + if (lp >= k) exit + lp = iaux(lp) + end do + swap = x(lp) + x(lp) = x(k) + x(k) = swap + ixswap = indx(lp) + indx(lp) = indx(k) + indx(k) = ixswap + lswap = iaux(lp) + iaux(lp) = iaux(k) + iaux(k) = lp + lp = lswap + k = k + 1 + enddo + return + end subroutine psb_ip_reord_z1i1 + + subroutine psb_ip_reord_z1i2(n,x,i1,i2,iaux) + integer(psb_ipk_), intent(in) :: n + integer(psb_ipk_) :: iaux(0:*) + complex(psb_dpk_) :: x(*) + integer(psb_ipk_) :: i1(*), i2(*) + + + integer(psb_ipk_) :: lswap, lp, k, isw1, isw2 + complex(psb_dpk_) :: swap + + lp = iaux(0) + k = 1 + do + if ((lp == 0).or.(k>n)) exit + do + if (lp >= k) exit + lp = iaux(lp) + end do + swap = x(lp) + x(lp) = x(k) + x(k) = swap + isw1 = i1(lp) + i1(lp) = i1(k) + i1(k) = isw1 + isw2 = i2(lp) + i2(lp) = i2(k) + i2(k) = isw2 + lswap = iaux(lp) + iaux(lp) = iaux(k) + iaux(k) = lp + lp = lswap + k = k + 1 + enddo + return + end subroutine psb_ip_reord_z1i2 + + subroutine psb_ip_reord_z1i3(n,x,i1,i2,i3,iaux) + integer(psb_ipk_), intent(in) :: n + integer(psb_ipk_) :: iaux(0:*) + complex(psb_dpk_) :: x(*) + integer(psb_ipk_) :: i1(*), i2(*), i3(*) + + + integer(psb_ipk_) :: lswap, lp, k, isw1, isw2, isw3 + complex(psb_dpk_) :: swap + + lp = iaux(0) + k = 1 + do + if ((lp == 0).or.(k>n)) exit + do + if (lp >= k) exit + lp = iaux(lp) + end do + swap = x(lp) + x(lp) = x(k) + x(k) = swap + isw1 = i1(lp) + i1(lp) = i1(k) + i1(k) = isw1 + isw2 = i2(lp) + i2(lp) = i2(k) + i2(k) = isw2 + isw3 = i3(lp) + i3(lp) = i3(k) + i3(k) = isw3 + lswap = iaux(lp) + iaux(lp) = iaux(k) + iaux(k) = lp + lp = lswap + k = k + 1 + enddo + return + end subroutine psb_ip_reord_z1i3 + +end module psb_z_ip_reord_mod diff --git a/base/modules/aux/psb_z_sort_mod.f90 b/base/modules/aux/psb_z_sort_mod.f90 index 15294b1c3..9fe7199dd 100644 --- a/base/modules/aux/psb_z_sort_mod.f90 +++ b/base/modules/aux/psb_z_sort_mod.f90 @@ -44,6 +44,15 @@ module psb_z_sort_mod use psb_const_mod + interface psb_msort_unique + subroutine psb_zmsort_u(x,nout,dir) + import + complex(psb_dpk_), intent(inout) :: x(:) + integer(psb_ipk_), intent(out) :: nout + integer(psb_ipk_), optional, intent(in) :: dir + end subroutine psb_zmsort_u + end interface psb_msort_unique + type psb_z_heap integer(psb_ipk_) :: last, dir complex(psb_dpk_), allocatable :: keys(:) @@ -544,7 +553,8 @@ contains implicit none class(psb_z_idx_heap), intent(inout) :: heap - integer(psb_ipk_), intent(out) :: index,info + integer(psb_ipk_), intent(out) :: index + integer(psb_ipk_), intent(out) :: info complex(psb_dpk_), intent(out) :: key diff --git a/base/modules/desc/psb_gen_block_map_mod.f90 b/base/modules/desc/psb_gen_block_map_mod.f90 index ff6aadae1..3afae81b8 100644 --- a/base/modules/desc/psb_gen_block_map_mod.f90 +++ b/base/modules/desc/psb_gen_block_map_mod.f90 @@ -90,7 +90,8 @@ module psb_gen_block_map_mod & block_get_fmt, block_l2gs1, block_l2gs2, block_l2gv1,& & block_l2gv2, block_g2ls1, block_g2ls2, block_g2lv1,& & block_g2lv2, block_g2ls1_ins, block_g2ls2_ins,& - & block_g2lv1_ins, block_g2lv2_ins, block_clone, block_reinit + & block_g2lv1_ins, block_g2lv2_ins, block_clone, block_reinit,& + & gen_block_search integer(psb_ipk_), private :: laddsz=500 @@ -992,7 +993,6 @@ contains subroutine block_fnd_owner(idx,iprc,idxmap,info) use psb_penv_mod - use psb_sort_mod implicit none integer(psb_ipk_), intent(in) :: idx(:) integer(psb_ipk_), allocatable, intent(out) :: iprc(:) @@ -1009,7 +1009,7 @@ contains return end if do i=1, nv - ip = psb_iblsrch(idx(i)-1,np+1,idxmap%vnl) + ip = gen_block_search(idx(i)-1,np+1,idxmap%vnl) iprc(i) = ip - 1 end do @@ -1226,52 +1226,50 @@ contains return end subroutine block_reinit -!!$ -!!$ subroutine block_reinit(idxmap,info) -!!$ use psb_penv_mod -!!$ use psb_error_mod -!!$ use psb_realloc_mod -!!$ implicit none -!!$ class(psb_gen_block_map), intent(inout) :: idxmap -!!$ integer(psb_ipk_), intent(out) :: info -!!$ integer(psb_ipk_) :: err_act, nr,nc,k, nl, ictxt -!!$ integer(psb_ipk_), allocatable :: idx(:),lidx(:) -!!$ character(len=20) :: name='block_reinit' -!!$ logical, parameter :: debug=.false. -!!$ -!!$ info = psb_success_ -!!$ call psb_get_erraction(err_act) -!!$ ictxt = idxmap%get_ctxt() -!!$ nr = idxmap%get_lr() -!!$ nc = idxmap%get_lc() -!!$ if (nc>nr) then -!!$ lidx = (/(k,k=nr+1,nc)/) -!!$ idx = (/(k,k=nr+1,nc)/) -!!$ call idxmap%l2gip(idx,info) -!!$ end if -!!$ if (info /= 0) & -!!$ & write(0,*) 'From l2gip',info -!!$ -!!$ call idxmap%init(ictxt,nr,info) -!!$ if (nc>nr) then -!!$ call idxmap%g2lip_ins(idx,info,lidx=lidx) -!!$ end if -!!$ -!!$ -!!$ if (info /= psb_success_) then -!!$ info = psb_err_from_subroutine_ -!!$ call psb_errpush(info,name) -!!$ goto 9999 -!!$ end if -!!$ call psb_erractionrestore(err_act) -!!$ return -!!$ -!!$ -!!$9999 call psb_error_handler(err_act) -!!$ -!!$ return -!!$ end subroutine block_reinit -!!$ + + ! + ! This is a purely internal version of "binary" search + ! specialized for gen_block usage. + ! + function gen_block_search(key,n,v) result(ipos) + implicit none + integer(psb_ipk_) :: ipos, key, n + integer(psb_ipk_) :: v(:) + + integer(psb_ipk_) :: lb, ub, m + + if (n < 5) then + ! don't bother with binary search for very + ! small vectors + ipos = 0 + do + if (ipos == n) return + if (key < v(ipos+1)) return + ipos = ipos + 1 + end do + else + lb = 1 + ub = n + ipos = -1 + + do while (lb <= ub) + m = (lb+ub)/2 + if (key==v(m)) then + ipos = m + return + else if (key < v(m)) then + ub = m-1 + else + lb = m + 1 + end if + enddo + if (v(ub) > key) then + ub = ub - 1 + end if + ipos = ub + endif + return + end function gen_block_search end module psb_gen_block_map_mod diff --git a/base/serial/sort/Makefile b/base/serial/sort/Makefile index 4dc7435a2..ddf08b96b 100644 --- a/base/serial/sort/Makefile +++ b/base/serial/sort/Makefile @@ -5,12 +5,13 @@ include ../../../Make.inc # BOBJS=psi_lcx_mod.o psi_alcx_mod.o psi_acx_mod.o IOBJS=psb_i_hsort_impl.o psb_i_isort_impl.o psb_i_msort_impl.o psb_i_qsort_impl.o +LOBJS=psb_l_hsort_impl.o psb_l_isort_impl.o psb_l_msort_impl.o psb_l_qsort_impl.o SOBJS=psb_s_hsort_impl.o psb_s_isort_impl.o psb_s_msort_impl.o psb_s_qsort_impl.o DOBJS=psb_d_hsort_impl.o psb_d_isort_impl.o psb_d_msort_impl.o psb_d_qsort_impl.o COBJS=psb_c_hsort_impl.o psb_c_isort_impl.o psb_c_msort_impl.o psb_c_qsort_impl.o ZOBJS=psb_z_hsort_impl.o psb_z_isort_impl.o psb_z_msort_impl.o psb_z_qsort_impl.o -OBJS=$(BOBJS) $(IOBJS) $(SOBJS) $(DOBJS) $(COBJS) $(ZOBJS) +OBJS=$(BOBJS) $(IOBJS) $(LOBJS) $(SOBJS) $(DOBJS) $(COBJS) $(ZOBJS) # # Where the library should go, and how it is called. @@ -34,7 +35,7 @@ lib: $(OBJS) # A bit excessive, but safe $(OBJS): $(MODDIR)/psb_base_mod.o -$(IOBJS) $(SOBJS) $(DOBJS) $(COBJS) $(ZOBJS): $(BOBJS) +$(IOBJS) $(LOBJS) $(SOBJS) $(DOBJS) $(COBJS) $(ZOBJS): $(BOBJS) clean: cleanobjs veryclean: cleanobjs diff --git a/base/serial/sort/psb_c_isort_impl.f90 b/base/serial/sort/psb_c_isort_impl.f90 index f5b7f9740..99353014a 100644 --- a/base/serial/sort/psb_c_isort_impl.f90 +++ b/base/serial/sort/psb_c_isort_impl.f90 @@ -48,7 +48,8 @@ subroutine psb_cisort(x,ix,dir,flag) integer(psb_ipk_), optional, intent(in) :: dir, flag integer(psb_ipk_), optional, intent(inout) :: ix(:) - integer(psb_ipk_) :: dir_, flag_, n, err_act, i + integer(psb_ipk_) :: dir_, flag_, err_act + integer(psb_ipk_) :: n, i integer(psb_ipk_) :: ierr(5) character(len=20) :: name diff --git a/base/serial/sort/psb_c_msort_impl.f90 b/base/serial/sort/psb_c_msort_impl.f90 index db5cbdf9f..c41cb9c82 100644 --- a/base/serial/sort/psb_c_msort_impl.f90 +++ b/base/serial/sort/psb_c_msort_impl.f90 @@ -41,6 +41,41 @@ ! Addison-Wesley ! + subroutine psb_cmsort_u(x,nout,dir) + use psb_c_sort_mod, psb_protect_name => psb_cmsort_u + use psb_error_mod + implicit none + complex(psb_spk_), intent(inout) :: x(:) + integer(psb_ipk_), intent(out) :: nout + integer(psb_ipk_), optional, intent(in) :: dir + + integer(psb_ipk_) :: n, k + integer(psb_ipk_) :: err_act + + integer(psb_ipk_) :: ierr(5) + character(len=20) :: name + + name='psb_msort_u' + call psb_erractionsave(err_act) + + n = size(x) + + call psb_msort(x,dir=dir) + nout = min(1,n) + do k=2,n + if (x(k) /= x(nout)) then + nout = nout + 1 + x(nout) = x(k) + endif + enddo + + return + +9999 call psb_error_handler(err_act) + + return + end subroutine psb_cmsort_u + diff --git a/base/serial/sort/psb_c_qsort_impl.f90 b/base/serial/sort/psb_c_qsort_impl.f90 index f22f3a862..169888cc0 100644 --- a/base/serial/sort/psb_c_qsort_impl.f90 +++ b/base/serial/sort/psb_c_qsort_impl.f90 @@ -48,8 +48,8 @@ subroutine psb_cqsort(x,ix,dir,flag) integer(psb_ipk_), optional, intent(in) :: dir, flag integer(psb_ipk_), optional, intent(inout) :: ix(:) - integer(psb_ipk_) :: dir_, flag_, n, err_act, i - + integer(psb_ipk_) :: dir_, flag_, err_act, i + integer(psb_ipk_) :: n integer(psb_ipk_) :: ierr(5) character(len=20) :: name @@ -150,7 +150,8 @@ subroutine psi_clqsrx_up(n,x,idx) ! .. 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_) :: n1, n2 + integer(psb_ipk_) :: ixt integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=32 integer(psb_ipk_) :: istack(nparms,maxstack) @@ -306,7 +307,8 @@ subroutine psi_clqsrx_dw(n,x,idx) ! .. 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_) :: n1, n2 + integer(psb_ipk_) :: ixt integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=32 integer(psb_ipk_) :: istack(nparms,maxstack) @@ -744,7 +746,8 @@ subroutine psi_calqsrx_up(n,x,idx) ! .. 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_) :: n1, n2 + integer(psb_ipk_) :: ixt integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=32 integer(psb_ipk_) :: istack(nparms,maxstack) @@ -899,7 +902,8 @@ subroutine psi_calqsrx_dw(n,x,idx) ! .. 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_) :: n1, n2 + integer(psb_ipk_) :: ixt integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=32 integer(psb_ipk_) :: istack(nparms,maxstack) @@ -1335,7 +1339,8 @@ subroutine psi_caqsrx_up(n,x,idx) real(psb_spk_) :: piv, xk complex(psb_spk_) :: xt integer(psb_ipk_) :: i, j, ilx, iux, istp, lpiv - integer(psb_ipk_) :: ixt, n1, n2 + integer(psb_ipk_) :: n1, n2 + integer(psb_ipk_) :: ixt integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=32 integer(psb_ipk_) :: istack(nparms,maxstack) @@ -1491,7 +1496,8 @@ subroutine psi_caqsrx_dw(n,x,idx) real(psb_spk_) :: piv, xk complex(psb_spk_) :: xt integer(psb_ipk_) :: i, j, ilx, iux, istp, lpiv - integer(psb_ipk_) :: ixt, n1, n2 + integer(psb_ipk_) :: n1, n2 + integer(psb_ipk_) :: ixt integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=32 integer(psb_ipk_) :: istack(nparms,maxstack) @@ -1644,7 +1650,8 @@ subroutine psi_caqsr_up(n,x) real(psb_spk_) :: piv, xk complex(psb_spk_) :: xt integer(psb_ipk_) :: i, j, ilx, iux, istp, lpiv - integer(psb_ipk_) :: ixt, n1, n2 + integer(psb_ipk_) :: n1, n2 + integer(psb_ipk_) :: ixt integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=32 integer(psb_ipk_) :: istack(nparms,maxstack) @@ -1784,7 +1791,8 @@ subroutine psi_caqsr_dw(n,x) real(psb_spk_) :: piv, xk complex(psb_spk_) :: xt integer(psb_ipk_) :: i, j, ilx, iux, istp, lpiv - integer(psb_ipk_) :: ixt, n1, n2 + integer(psb_ipk_) :: n1, n2 + integer(psb_ipk_) :: ixt integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=32 integer(psb_ipk_) :: istack(nparms,maxstack) diff --git a/base/serial/sort/psb_d_isort_impl.f90 b/base/serial/sort/psb_d_isort_impl.f90 index ddf41a381..262931c36 100644 --- a/base/serial/sort/psb_d_isort_impl.f90 +++ b/base/serial/sort/psb_d_isort_impl.f90 @@ -48,7 +48,8 @@ subroutine psb_disort(x,ix,dir,flag) integer(psb_ipk_), optional, intent(in) :: dir, flag integer(psb_ipk_), optional, intent(inout) :: ix(:) - integer(psb_ipk_) :: dir_, flag_, n, err_act, i + integer(psb_ipk_) :: dir_, flag_, err_act + integer(psb_ipk_) :: n, i integer(psb_ipk_) :: ierr(5) character(len=20) :: name diff --git a/base/serial/sort/psb_d_msort_impl.f90 b/base/serial/sort/psb_d_msort_impl.f90 index 5bda802dd..b0d7f8b54 100644 --- a/base/serial/sort/psb_d_msort_impl.f90 +++ b/base/serial/sort/psb_d_msort_impl.f90 @@ -41,6 +41,99 @@ ! Addison-Wesley ! + subroutine psb_dmsort_u(x,nout,dir) + use psb_d_sort_mod, psb_protect_name => psb_dmsort_u + use psb_error_mod + implicit none + real(psb_dpk_), intent(inout) :: x(:) + integer(psb_ipk_), intent(out) :: nout + integer(psb_ipk_), optional, intent(in) :: dir + + integer(psb_ipk_) :: n, k + integer(psb_ipk_) :: err_act + + integer(psb_ipk_) :: ierr(5) + character(len=20) :: name + + name='psb_msort_u' + call psb_erractionsave(err_act) + + n = size(x) + + call psb_msort(x,dir=dir) + nout = min(1,n) + do k=2,n + if (x(k) /= x(nout)) then + nout = nout + 1 + x(nout) = x(k) + endif + enddo + + return + +9999 call psb_error_handler(err_act) + + return + end subroutine psb_dmsort_u + + + function psb_dbsrch(key,n,v) result(ipos) + use psb_d_sort_mod, psb_protect_name => psb_dbsrch + implicit none + integer(psb_ipk_) :: ipos, n + real(psb_dpk_) :: key + real(psb_dpk_) :: v(:) + + integer(psb_ipk_) :: lb, ub, m, i + + ipos = -1 + if (n<5) then + do i=1,n + if (key.eq.v(i)) then + ipos = i + return + end if + enddo + return + end if + + lb = 1 + ub = n + + do while (lb.le.ub) + m = (lb+ub)/2 + if (key.eq.v(m)) then + ipos = m + lb = ub + 1 + else if (key < v(m)) then + ub = m-1 + else + lb = m + 1 + end if + enddo + return + end function psb_dbsrch + + function psb_dssrch(key,n,v) result(ipos) + use psb_d_sort_mod, psb_protect_name => psb_dssrch + implicit none + integer(psb_ipk_) :: ipos, n + real(psb_dpk_) :: key + real(psb_dpk_) :: v(:) + + integer(psb_ipk_) :: i + + ipos = -1 + do i=1,n + if (key.eq.v(i)) then + ipos = i + return + end if + enddo + + return + end function psb_dssrch + subroutine psb_dmsort(x,ix,dir,flag) use psb_d_sort_mod, psb_protect_name => psb_dmsort use psb_error_mod diff --git a/base/serial/sort/psb_d_qsort_impl.f90 b/base/serial/sort/psb_d_qsort_impl.f90 index 105e94d00..1ccbbf5fa 100644 --- a/base/serial/sort/psb_d_qsort_impl.f90 +++ b/base/serial/sort/psb_d_qsort_impl.f90 @@ -48,8 +48,8 @@ subroutine psb_dqsort(x,ix,dir,flag) integer(psb_ipk_), optional, intent(in) :: dir, flag integer(psb_ipk_), optional, intent(inout) :: ix(:) - integer(psb_ipk_) :: dir_, flag_, n, err_act, i - + integer(psb_ipk_) :: dir_, flag_, err_act, i + integer(psb_ipk_) :: n integer(psb_ipk_) :: ierr(5) character(len=20) :: name @@ -140,7 +140,8 @@ subroutine psi_dqsrx_up(n,x,idx) ! .. 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_) :: n1, n2 + integer(psb_ipk_) :: ixt integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=60 integer(psb_ipk_) :: istack(nparms,maxstack) @@ -293,7 +294,8 @@ subroutine psi_dqsrx_dw(n,x,idx) ! .. 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_) :: n1, n2 + integer(psb_ipk_) :: ixt integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=60 integer(psb_ipk_) :: istack(nparms,maxstack) @@ -731,7 +733,8 @@ subroutine psi_daqsrx_up(n,x,idx) real(psb_dpk_) :: piv, xk real(psb_dpk_) :: xt integer(psb_ipk_) :: i, j, ilx, iux, istp, lpiv - integer(psb_ipk_) :: ixt, n1, n2 + integer(psb_ipk_) :: n1, n2 + integer(psb_ipk_) :: ixt integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=60 integer(psb_ipk_) :: istack(nparms,maxstack) @@ -887,7 +890,8 @@ subroutine psi_daqsrx_dw(n,x,idx) real(psb_dpk_) :: piv, xk real(psb_dpk_) :: xt integer(psb_ipk_) :: i, j, ilx, iux, istp, lpiv - integer(psb_ipk_) :: ixt, n1, n2 + integer(psb_ipk_) :: n1, n2 + integer(psb_ipk_) :: ixt integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=60 integer(psb_ipk_) :: istack(nparms,maxstack) @@ -1040,7 +1044,8 @@ subroutine psi_daqsr_up(n,x) real(psb_dpk_) :: piv, xk real(psb_dpk_) :: xt integer(psb_ipk_) :: i, j, ilx, iux, istp, lpiv - integer(psb_ipk_) :: ixt, n1, n2 + integer(psb_ipk_) :: n1, n2 + integer(psb_ipk_) :: ixt integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=60 integer(psb_ipk_) :: istack(nparms,maxstack) @@ -1180,7 +1185,8 @@ subroutine psi_daqsr_dw(n,x) real(psb_dpk_) :: piv, xk real(psb_dpk_) :: xt integer(psb_ipk_) :: i, j, ilx, iux, istp, lpiv - integer(psb_ipk_) :: ixt, n1, n2 + integer(psb_ipk_) :: n1, n2 + integer(psb_ipk_) :: ixt integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=60 integer(psb_ipk_) :: istack(nparms,maxstack) diff --git a/base/serial/sort/psb_i_isort_impl.f90 b/base/serial/sort/psb_i_isort_impl.f90 index ccaec9b7c..d3a82950e 100644 --- a/base/serial/sort/psb_i_isort_impl.f90 +++ b/base/serial/sort/psb_i_isort_impl.f90 @@ -48,7 +48,8 @@ subroutine psb_iisort(x,ix,dir,flag) integer(psb_ipk_), optional, intent(in) :: dir, flag integer(psb_ipk_), optional, intent(inout) :: ix(:) - integer(psb_ipk_) :: dir_, flag_, n, err_act, i + integer(psb_ipk_) :: dir_, flag_, err_act + integer(psb_ipk_) :: n, i integer(psb_ipk_) :: ierr(5) character(len=20) :: name diff --git a/base/serial/sort/psb_i_msort_impl.f90 b/base/serial/sort/psb_i_msort_impl.f90 index 934e2e4f1..db03c8b9b 100644 --- a/base/serial/sort/psb_i_msort_impl.f90 +++ b/base/serial/sort/psb_i_msort_impl.f90 @@ -40,8 +40,8 @@ ! Data Structures and Algorithms ! Addison-Wesley ! - logical function psb_isaperm(n,eip) - use psb_i_sort_mod, psb_protect_name => psb_isaperm + logical function psb_iisaperm(n,eip) + use psb_i_sort_mod, psb_protect_name => psb_iisaperm implicit none integer(psb_ipk_), intent(in) :: n @@ -50,7 +50,7 @@ integer(psb_ipk_) :: i,j,m, info - psb_isaperm = .true. + psb_iisaperm = .true. if (n <= 0) return allocate(ip(n), stat=info) if (info /= psb_success_) return @@ -61,7 +61,7 @@ ip(i) = eip(i) if ((ip(i) < 1).or.(ip(i) > n)) then write(psb_err_unit,*) 'Out of bounds in isaperm' ,ip(i), n - psb_isaperm = .false. + psb_iisaperm = .false. return endif enddo @@ -85,7 +85,7 @@ enddo ip(m) = abs(ip(m)) if (j /= m) then - psb_isaperm = .false. + psb_iisaperm = .false. goto 9999 endif end if @@ -93,61 +93,67 @@ 9999 continue return - end function psb_isaperm + end function psb_iisaperm - function psb_iblsrch(key,n,v) result(ipos) - use psb_i_sort_mod, psb_protect_name => psb_iblsrch - implicit none - integer(psb_ipk_) :: ipos, key, n - integer(psb_ipk_) :: v(:) - integer(psb_ipk_) :: lb, ub, m - - if (n < 5) then - ! don't bother with binary search for very - ! small vectors - ipos = 0 - do - if (ipos == n) return - if (key < v(ipos+1)) return - ipos = ipos + 1 - end do - else - lb = 1 - ub = n - ipos = -1 - - do while (lb <= ub) - m = (lb+ub)/2 - if (key==v(m)) then - ipos = m - return - else if (key < v(m)) then - ub = m-1 - else - lb = m + 1 - end if - enddo - if (v(ub) > key) then -!!$ write(0,*) 'Check: ',ub,v(ub),key - ub = ub - 1 - end if - ipos = ub - endif + subroutine psb_imsort_u(x,nout,dir) + use psb_i_sort_mod, psb_protect_name => psb_imsort_u + use psb_error_mod + implicit none + integer(psb_ipk_), intent(inout) :: x(:) + integer(psb_ipk_), intent(out) :: nout + integer(psb_ipk_), optional, intent(in) :: dir + + integer(psb_ipk_) :: n, k + integer(psb_ipk_) :: err_act + + integer(psb_ipk_) :: ierr(5) + character(len=20) :: name + + name='psb_msort_u' + call psb_erractionsave(err_act) + + n = size(x) + + call psb_msort(x,dir=dir) + nout = min(1,n) + do k=2,n + if (x(k) /= x(nout)) then + nout = nout + 1 + x(nout) = x(k) + endif + enddo + return - end function psb_iblsrch + +9999 call psb_error_handler(err_act) + + return + end subroutine psb_imsort_u + function psb_ibsrch(key,n,v) result(ipos) use psb_i_sort_mod, psb_protect_name => psb_ibsrch implicit none - integer(psb_ipk_) :: ipos, key, n + integer(psb_ipk_) :: ipos, n + integer(psb_ipk_) :: key integer(psb_ipk_) :: v(:) - integer(psb_ipk_) :: lb, ub, m + integer(psb_ipk_) :: lb, ub, m, i + ipos = -1 + if (n<5) then + do i=1,n + if (key.eq.v(i)) then + ipos = i + return + end if + enddo + return + end if + lb = 1 ub = n - ipos = -1 do while (lb.le.ub) m = (lb+ub)/2 @@ -166,7 +172,8 @@ function psb_issrch(key,n,v) result(ipos) use psb_i_sort_mod, psb_protect_name => psb_issrch implicit none - integer(psb_ipk_) :: ipos, key, n + integer(psb_ipk_) :: ipos, n + integer(psb_ipk_) :: key integer(psb_ipk_) :: v(:) integer(psb_ipk_) :: i @@ -182,56 +189,6 @@ return end function psb_issrch - - subroutine psb_imsort_u(x,nout,dir) - use psb_i_sort_mod, psb_protect_name => psb_imsort_u - use psb_error_mod - implicit none - integer(psb_ipk_), intent(inout) :: x(:) - integer(psb_ipk_), intent(out) :: nout - integer(psb_ipk_), optional, intent(in) :: dir - - integer(psb_ipk_) :: dir_, n, err_act, k - - integer(psb_ipk_) :: ierr(5) - character(len=20) :: name - - name='psb_msort_u' - call psb_erractionsave(err_act) - - if (present(dir)) then - dir_ = dir - else - dir_= psb_sort_up_ - end if - select case(dir_) - case( psb_sort_up_, psb_sort_down_, psb_asort_up_, psb_asort_down_) - ! OK keep going - case default - ierr(1) = 3; ierr(2) = dir_; - call psb_errpush(psb_err_input_value_invalid_i_,name,i_err=ierr) - goto 9999 - end select - - n = size(x) - - call psb_imsort(x,dir=dir_) - nout = min(1,n) - do k=2,n - if (x(k) /= x(nout)) then - nout = nout + 1 - x(nout) = x(k) - endif - enddo - - return - -9999 call psb_error_handler(err_act) - - return - end subroutine psb_imsort_u - - subroutine psb_imsort(x,ix,dir,flag) use psb_i_sort_mod, psb_protect_name => psb_imsort use psb_error_mod diff --git a/base/serial/sort/psb_i_qsort_impl.f90 b/base/serial/sort/psb_i_qsort_impl.f90 index a16ef9108..355bd47cd 100644 --- a/base/serial/sort/psb_i_qsort_impl.f90 +++ b/base/serial/sort/psb_i_qsort_impl.f90 @@ -48,8 +48,8 @@ subroutine psb_iqsort(x,ix,dir,flag) integer(psb_ipk_), optional, intent(in) :: dir, flag integer(psb_ipk_), optional, intent(inout) :: ix(:) - integer(psb_ipk_) :: dir_, flag_, n, err_act, i - + integer(psb_ipk_) :: dir_, flag_, err_act, i + integer(psb_ipk_) :: n integer(psb_ipk_) :: ierr(5) character(len=20) :: name @@ -140,7 +140,8 @@ subroutine psi_iqsrx_up(n,x,idx) ! .. 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_) :: n1, n2 + integer(psb_ipk_) :: ixt integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=72 integer(psb_ipk_) :: istack(nparms,maxstack) @@ -293,7 +294,8 @@ subroutine psi_iqsrx_dw(n,x,idx) ! .. 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_) :: n1, n2 + integer(psb_ipk_) :: ixt integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=72 integer(psb_ipk_) :: istack(nparms,maxstack) @@ -731,7 +733,8 @@ subroutine psi_iaqsrx_up(n,x,idx) integer(psb_ipk_) :: piv, xk integer(psb_ipk_) :: xt integer(psb_ipk_) :: i, j, ilx, iux, istp, lpiv - integer(psb_ipk_) :: ixt, n1, n2 + integer(psb_ipk_) :: n1, n2 + integer(psb_ipk_) :: ixt integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=72 integer(psb_ipk_) :: istack(nparms,maxstack) @@ -887,7 +890,8 @@ subroutine psi_iaqsrx_dw(n,x,idx) integer(psb_ipk_) :: piv, xk integer(psb_ipk_) :: xt integer(psb_ipk_) :: i, j, ilx, iux, istp, lpiv - integer(psb_ipk_) :: ixt, n1, n2 + integer(psb_ipk_) :: n1, n2 + integer(psb_ipk_) :: ixt integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=72 integer(psb_ipk_) :: istack(nparms,maxstack) @@ -1040,7 +1044,8 @@ subroutine psi_iaqsr_up(n,x) integer(psb_ipk_) :: piv, xk integer(psb_ipk_) :: xt integer(psb_ipk_) :: i, j, ilx, iux, istp, lpiv - integer(psb_ipk_) :: ixt, n1, n2 + integer(psb_ipk_) :: n1, n2 + integer(psb_ipk_) :: ixt integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=72 integer(psb_ipk_) :: istack(nparms,maxstack) @@ -1180,7 +1185,8 @@ subroutine psi_iaqsr_dw(n,x) integer(psb_ipk_) :: piv, xk integer(psb_ipk_) :: xt integer(psb_ipk_) :: i, j, ilx, iux, istp, lpiv - integer(psb_ipk_) :: ixt, n1, n2 + integer(psb_ipk_) :: n1, n2 + integer(psb_ipk_) :: ixt integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=72 integer(psb_ipk_) :: istack(nparms,maxstack) diff --git a/base/serial/sort/psb_l_hsort_impl.f90 b/base/serial/sort/psb_l_hsort_impl.f90 new file mode 100644 index 000000000..951499ead --- /dev/null +++ b/base/serial/sort/psb_l_hsort_impl.f90 @@ -0,0 +1,678 @@ +! +! Parallel Sparse BLAS version 3.5 +! (C) Copyright 2006-2018 +! Salvatore Filippone +! Alfredo Buttari +! +! Redistribution and use in source and binary forms, with or without +! modification, are permitted provided that the following conditions +! are met: +! 1. Redistributions of source code must retain the above copyright +! notice, this list of conditions and the following disclaimer. +! 2. Redistributions in binary form must reproduce the above copyright +! notice, this list of conditions, and the following disclaimer in the +! documentation and/or other materials provided with the distribution. +! 3. The name of the PSBLAS group or the names of its contributors may +! not be used to endorse or promote products derived from this +! software without specific written permission. +! +! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +! ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED +! TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR +! PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE PSBLAS GROUP OR ITS CONTRIBUTORS +! BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +! CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +! SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +! INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +! CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +! ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +! POSSIBILITY OF SUCH DAMAGE. +! +! +! +! The merge-sort and quicksort routines are implemented in the +! serial/aux directory +! References: +! D. Knuth +! The Art of Computer Programming, vol. 3 +! Addison-Wesley +! +! Aho, Hopcroft, Ullman +! Data Structures and Algorithms +! Addison-Wesley +! +subroutine psb_lhsort(x,ix,dir,flag) + use psb_l_sort_mod, psb_protect_name => psb_lhsort + use psb_error_mod + implicit none + integer(psb_lpk_), intent(inout) :: x(:) + integer(psb_ipk_), optional, intent(in) :: dir, flag + integer(psb_lpk_), optional, intent(inout) :: ix(:) + + integer(psb_ipk_) :: dir_, flag_, n, i, l, err_act,info + integer(psb_lpk_) :: key + integer(psb_lpk_) :: index + + integer(psb_ipk_) :: ierr(5) + character(len=20) :: name + + name='psb_hsort' + call psb_erractionsave(err_act) + + if (present(flag)) then + flag_ = flag + else + flag_ = psb_sort_ovw_idx_ + end if + select case(flag_) + case( psb_sort_ovw_idx_, psb_sort_keep_idx_) + ! OK keep going + case default + ierr(1) = 4; ierr(2) = flag_; + call psb_errpush(psb_err_input_value_invalid_i_,name,i_err=ierr) + goto 9999 + end select + + if (present(dir)) then + dir_ = dir + else + dir_= psb_sort_up_ + end if + + select case(dir_) + case(psb_sort_up_,psb_sort_down_) + ! OK + case (psb_asort_up_,psb_asort_down_) + ! OK + case default + ierr(1) = 3; ierr(2) = dir_; + call psb_errpush(psb_err_input_value_invalid_i_,name,i_err=ierr) + goto 9999 + end select + + n = size(x) + + ! + ! Dirty trick to sort with heaps: if we want + ! to sort in place upwards, first we set up a heap so that + ! we can easily get the LARGEST element, then we take it out + ! and put it in the last entry, and so on. + ! So, we invert dir_ + ! + dir_ = -dir_ + + if (present(ix)) then + if (size(ix) < n) then + ierr(1) = 2; ierr(2) = size(ix); + call psb_errpush(psb_err_input_asize_invalid_i_,name,i_err=ierr) + goto 9999 + end if + if (flag_ == psb_sort_ovw_idx_) then + do i=1, n + ix(i) = i + end do + end if + l = 0 + do i=1, n + key = x(i) + index = ix(i) + call psi_l_idx_insert_heap(key,index,l,x,ix,dir_,info) + if (l /= i) then + write(psb_err_unit,*) 'Mismatch while heapifying ! ' + end if + end do + do i=n, 2, -1 + call psi_l_idx_heap_get_first(key,index,l,x,ix,dir_,info) + if (l /= i-1) then + write(psb_err_unit,*) 'Mismatch while pulling out of heap ',l,i + end if + x(i) = key + ix(i) = index + end do + else if (.not.present(ix)) then + l = 0 + do i=1, n + key = x(i) + call psi_l_insert_heap(key,l,x,dir_,info) + if (l /= i) then + write(psb_err_unit,*) 'Mismatch while heapifying ! ',l,i + end if + end do + do i=n, 2, -1 + call psi_l_heap_get_first(key,l,x,dir_,info) + if (l /= i-1) then + write(psb_err_unit,*) 'Mismatch while pulling out of heap ',l,i + end if + x(i) = key + end do + end if + + + return + +9999 call psb_error_handler(err_act) + + return +end subroutine psb_lhsort + + + +! +! These are packaged so that they can be used to implement +! a heapsort, should the need arise +! +! +! Programming note: +! In the implementation of the heap_get_first function +! we have code like this +! +! if ( ( heap(2*i) < heap(2*i+1) ) .or.& +! & (2*i == last)) then +! j = 2*i +! else +! j = 2*i + 1 +! end if +! +! It looks like the 2*i+1 could overflow the array, but this +! is not true because there is a guard statement +! if (i>last/2) exit +! and because last has just been reduced by 1 when defining the return value, +! therefore 2*i+1 may be greater than the current value of last, +! but cannot be greater than the value of last when the routine was entered +! hence it is safe. +! +! +! + +subroutine psi_l_insert_heap(key,last,heap,dir,info) + use psb_l_sort_mod, psb_protect_name => psi_l_insert_heap + implicit none + + ! + ! Input: + ! key: the new value + ! last: pointer to the last occupied element in heap + ! heap: the heap + ! dir: sorting direction + + integer(psb_lpk_), intent(in) :: key + integer(psb_ipk_), intent(in) :: dir + integer(psb_lpk_), intent(inout) :: heap(:) + integer(psb_lpk_), intent(inout) :: last + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_) :: i, i2 + integer(psb_lpk_) :: temp + + info = psb_success_ + if (last < 0) then + write(psb_err_unit,*) 'Invalid last in heap ',last + info = last + return + endif + last = last + 1 + if (last > size(heap)) then + write(psb_err_unit,*) 'out of bounds ' + info = -1 + return + end if + i = last + heap(i) = key + + select case(dir) + case (psb_sort_up_) + + do + if (i<=1) exit + i2 = i/2 + if (heap(i) < heap(i2)) then + temp = heap(i) + heap(i) = heap(i2) + heap(i2) = temp + i = i2 + else + exit + end if + end do + + + case (psb_sort_down_) + + do + if (i<=1) exit + i2 = i/2 + if (heap(i) > heap(i2)) then + temp = heap(i) + heap(i) = heap(i2) + heap(i2) = temp + i = i2 + else + exit + end if + end do + + case (psb_asort_up_) + + do + if (i<=1) exit + i2 = i/2 + if (abs(heap(i)) < abs(heap(i2))) then + temp = heap(i) + heap(i) = heap(i2) + heap(i2) = temp + i = i2 + else + exit + end if + end do + + + case (psb_asort_down_) + + do + if (i<=1) exit + i2 = i/2 + if (abs(heap(i)) > abs(heap(i2))) then + temp = heap(i) + heap(i) = heap(i2) + heap(i2) = temp + i = i2 + else + exit + end if + end do + + + case default + write(psb_err_unit,*) 'Invalid direction in heap ',dir + end select + + return +end subroutine psi_l_insert_heap + + +subroutine psi_l_heap_get_first(key,last,heap,dir,info) + use psb_l_sort_mod, psb_protect_name => psi_l_heap_get_first + implicit none + + integer(psb_lpk_), intent(inout) :: key + integer(psb_lpk_), intent(inout) :: last + integer(psb_ipk_), intent(in) :: dir + integer(psb_lpk_), intent(inout) :: heap(:) + integer(psb_ipk_), intent(out) :: info + + integer(psb_ipk_) :: i, j + integer(psb_lpk_) :: temp + + + info = psb_success_ + if (last <= 0) then + key = 0 + info = -1 + return + endif + + key = heap(1) + heap(1) = heap(last) + last = last - 1 + + select case(dir) + case (psb_sort_up_) + + i = 1 + do + if (i > (last/2)) exit + if ( (heap(2*i) < heap(2*i+1)) .or.& + & (2*i == last)) then + j = 2*i + else + j = 2*i + 1 + end if + + if (heap(i) > heap(j)) then + temp = heap(i) + heap(i) = heap(j) + heap(j) = temp + i = j + else + exit + end if + end do + + + case (psb_sort_down_) + + i = 1 + do + if (i > (last/2)) exit + if ( (heap(2*i) > heap(2*i+1)) .or.& + & (2*i == last)) then + j = 2*i + else + j = 2*i + 1 + end if + + if (heap(i) < heap(j)) then + temp = heap(i) + heap(i) = heap(j) + heap(j) = temp + i = j + else + exit + end if + end do + + case (psb_asort_up_) + + i = 1 + do + if (i > (last/2)) exit + if ( (abs(heap(2*i)) < abs(heap(2*i+1))) .or.& + & (2*i == last)) then + j = 2*i + else + j = 2*i + 1 + end if + + if (abs(heap(i)) > abs(heap(j))) then + temp = heap(i) + heap(i) = heap(j) + heap(j) = temp + i = j + else + exit + end if + end do + + + case (psb_asort_down_) + + i = 1 + do + if (i > (last/2)) exit + if ( (abs(heap(2*i)) > abs(heap(2*i+1))) .or.& + & (2*i == last)) then + j = 2*i + else + j = 2*i + 1 + end if + + if (abs(heap(i)) < abs(heap(j))) then + temp = heap(i) + heap(i) = heap(j) + heap(j) = temp + i = j + else + exit + end if + end do + + case default + write(psb_err_unit,*) 'Invalid direction in heap ',dir + end select + + return +end subroutine psi_l_heap_get_first + + +subroutine psi_l_idx_insert_heap(key,index,last,heap,idxs,dir,info) + use psb_l_sort_mod, psb_protect_name => psi_l_idx_insert_heap + + implicit none + ! + ! Input: + ! key: the new value + ! index: the new index + ! last: pointer to the last occupied element in heap + ! heap: the heap + ! idxs: the indices + ! dir: sorting direction + + integer(psb_lpk_), intent(in) :: key + integer(psb_ipk_), intent(in) :: index,dir + integer(psb_lpk_), intent(inout) :: heap(:) + integer(psb_ipk_), intent(inout) :: idxs(:),last + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_) :: i, i2, itemp + integer(psb_lpk_) :: temp + + info = psb_success_ + if (last < 0) then + write(psb_err_unit,*) 'Invalid last in heap ',last + info = last + return + endif + + last = last + 1 + if (last > size(heap)) then + write(psb_err_unit,*) 'out of bounds ' + info = -1 + return + end if + + i = last + heap(i) = key + idxs(i) = index + + select case(dir) + case (psb_sort_up_) + + do + if (i<=1) exit + i2 = i/2 + if (heap(i) < heap(i2)) then + itemp = idxs(i) + idxs(i) = idxs(i2) + idxs(i2) = itemp + temp = heap(i) + heap(i) = heap(i2) + heap(i2) = temp + i = i2 + else + exit + end if + end do + + + case (psb_sort_down_) + + do + if (i<=1) exit + i2 = i/2 + if (heap(i) > heap(i2)) then + itemp = idxs(i) + idxs(i) = idxs(i2) + idxs(i2) = itemp + temp = heap(i) + heap(i) = heap(i2) + heap(i2) = temp + i = i2 + else + exit + end if + end do + + case (psb_asort_up_) + + do + if (i<=1) exit + i2 = i/2 + if (abs(heap(i)) < abs(heap(i2))) then + itemp = idxs(i) + idxs(i) = idxs(i2) + idxs(i2) = itemp + temp = heap(i) + heap(i) = heap(i2) + heap(i2) = temp + i = i2 + else + exit + end if + end do + + + case (psb_asort_down_) + + do + if (i<=1) exit + i2 = i/2 + if (abs(heap(i)) > abs(heap(i2))) then + itemp = idxs(i) + idxs(i) = idxs(i2) + idxs(i2) = itemp + temp = heap(i) + heap(i) = heap(i2) + heap(i2) = temp + i = i2 + else + exit + end if + end do + + + case default + write(psb_err_unit,*) 'Invalid direction in heap ',dir + end select + + return +end subroutine psi_l_idx_insert_heap + +subroutine psi_l_idx_heap_get_first(key,index,last,heap,idxs,dir,info) + use psb_l_sort_mod, psb_protect_name => psi_l_idx_heap_get_first + implicit none + + integer(psb_lpk_), intent(inout) :: heap(:) + integer(psb_ipk_), intent(out) :: index,info + integer(psb_ipk_), intent(inout) :: last,idxs(:) + integer(psb_ipk_), intent(in) :: dir + integer(psb_lpk_), intent(out) :: key + + integer(psb_ipk_) :: i, j,itemp + integer(psb_lpk_) :: temp + + info = psb_success_ + if (last <= 0) then + key = 0 + index = 0 + info = -1 + return + endif + + key = heap(1) + index = idxs(1) + heap(1) = heap(last) + idxs(1) = idxs(last) + last = last - 1 + + select case(dir) + case (psb_sort_up_) + + i = 1 + do + if (i > (last/2)) exit + if ( (heap(2*i) < heap(2*i+1)) .or.& + & (2*i == last)) then + j = 2*i + else + j = 2*i + 1 + end if + + if (heap(i) > heap(j)) then + itemp = idxs(i) + idxs(i) = idxs(j) + idxs(j) = itemp + temp = heap(i) + heap(i) = heap(j) + heap(j) = temp + i = j + else + exit + end if + end do + + + case (psb_sort_down_) + + i = 1 + do + if (i > (last/2)) exit + if ( (heap(2*i) > heap(2*i+1)) .or.& + & (2*i == last)) then + j = 2*i + else + j = 2*i + 1 + end if + + if (heap(i) < heap(j)) then + itemp = idxs(i) + idxs(i) = idxs(j) + idxs(j) = itemp + temp = heap(i) + heap(i) = heap(j) + heap(j) = temp + i = j + else + exit + end if + end do + + case (psb_asort_up_) + + i = 1 + do + if (i > (last/2)) exit + if ( (abs(heap(2*i)) < abs(heap(2*i+1))) .or.& + & (2*i == last)) then + j = 2*i + else + j = 2*i + 1 + end if + + if (abs(heap(i)) > abs(heap(j))) then + itemp = idxs(i) + idxs(i) = idxs(j) + idxs(j) = itemp + temp = heap(i) + heap(i) = heap(j) + heap(j) = temp + i = j + else + exit + end if + end do + + + case (psb_asort_down_) + + i = 1 + do + if (i > (last/2)) exit + if ( (abs(heap(2*i)) > abs(heap(2*i+1))) .or.& + & (2*i == last)) then + j = 2*i + else + j = 2*i + 1 + end if + + if (abs(heap(i)) < abs(heap(j))) then + itemp = idxs(i) + idxs(i) = idxs(j) + idxs(j) = itemp + temp = heap(i) + heap(i) = heap(j) + heap(j) = temp + i = j + else + exit + end if + end do + + case default + write(psb_err_unit,*) 'Invalid direction in heap ',dir + end select + + return +end subroutine psi_l_idx_heap_get_first + + + + diff --git a/base/serial/sort/psb_l_isort_impl.f90 b/base/serial/sort/psb_l_isort_impl.f90 new file mode 100644 index 000000000..be0b98f2e --- /dev/null +++ b/base/serial/sort/psb_l_isort_impl.f90 @@ -0,0 +1,341 @@ +! +! Parallel Sparse BLAS version 3.5 +! (C) Copyright 2006-2018 +! Salvatore Filippone +! Alfredo Buttari +! +! Redistribution and use in source and binary forms, with or without +! modification, are permitted provided that the following conditions +! are met: +! 1. Redistributions of source code must retain the above copyright +! notice, this list of conditions and the following disclaimer. +! 2. Redistributions in binary form must reproduce the above copyright +! notice, this list of conditions, and the following disclaimer in the +! documentation and/or other materials provided with the distribution. +! 3. The name of the PSBLAS group or the names of its contributors may +! not be used to endorse or promote products derived from this +! software without specific written permission. +! +! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +! ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED +! TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR +! PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE PSBLAS GROUP OR ITS CONTRIBUTORS +! BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +! CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +! SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +! INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +! CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +! ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +! POSSIBILITY OF SUCH DAMAGE. +! +! +! +! The insertion sort routines +! References: +! D. Knuth +! The Art of Computer Programming, vol. 3 +! Addison-Wesley +! +! Aho, Hopcroft, Ullman +! Data Structures and Algorithms +! Addison-Wesley +! +subroutine psb_lisort(x,ix,dir,flag) + use psb_l_sort_mod, psb_protect_name => psb_lisort + use psb_error_mod + implicit none + integer(psb_lpk_), intent(inout) :: x(:) + integer(psb_ipk_), optional, intent(in) :: dir, flag + integer(psb_lpk_), optional, intent(inout) :: ix(:) + + integer(psb_ipk_) :: dir_, flag_, err_act + integer(psb_lpk_) :: n, i + + integer(psb_ipk_) :: ierr(5) + character(len=20) :: name + + name='psb_lisort' + call psb_erractionsave(err_act) + + if (present(flag)) then + flag_ = flag + else + flag_ = psb_sort_ovw_idx_ + end if + select case(flag_) + case( psb_sort_ovw_idx_, psb_sort_keep_idx_) + ! OK keep going + case default + ierr(1) = 4; ierr(2) = flag_; + call psb_errpush(psb_err_input_value_invalid_i_,name,i_err=ierr) + goto 9999 + end select + + if (present(dir)) then + dir_ = dir + else + dir_= psb_sort_up_ + end if + + n = size(x) + + if (present(ix)) then + if (size(ix) < n) then + ierr(1) = 2; ierr(2) = size(ix); + call psb_errpush(psb_err_input_asize_invalid_i_,name,i_err=ierr) + goto 9999 + end if + if (flag_==psb_sort_ovw_idx_) then + do i=1,n + ix(i) = i + end do + end if + + select case(dir_) + case (psb_sort_up_) + call psi_lisrx_up(n,x,ix) + case (psb_sort_down_) + call psi_lisrx_dw(n,x,ix) + case (psb_asort_up_) + call psi_laisrx_up(n,x,ix) + case (psb_asort_down_) + call psi_laisrx_dw(n,x,ix) + case default + ierr(1) = 3; ierr(2) = dir_; + call psb_errpush(psb_err_input_value_invalid_i_,name,i_err=ierr) + goto 9999 + end select + else + select case(dir_) + case (psb_sort_up_) + call psi_lisr_up(n,x) + case (psb_sort_down_) + call psi_lisr_dw(n,x) + case (psb_asort_up_) + call psi_laisr_up(n,x) + case (psb_asort_down_) + call psi_laisr_dw(n,x) + case default + ierr(1) = 3; ierr(2) = dir_; + call psb_errpush(psb_err_input_value_invalid_i_,name,i_err=ierr) + goto 9999 + end select + + end if + + return + +9999 call psb_error_handler(err_act) + + return +end subroutine psb_lisort + +subroutine psi_lisrx_up(n,x,idx) + use psb_l_sort_mod, psb_protect_name => psi_lisrx_up + use psb_error_mod + implicit none + integer(psb_lpk_), intent(inout) :: x(:) + integer(psb_lpk_), intent(inout) :: idx(:) + integer(psb_lpk_), intent(in) :: n + integer(psb_lpk_) :: i,j,ix + integer(psb_lpk_) :: xx + + do j=n-1,1,-1 + if (x(j+1) < x(j)) then + xx = x(j) + ix = idx(j) + i=j+1 + do + x(i-1) = x(i) + idx(i-1) = idx(i) + i = i+1 + if (i>n) exit + if (x(i) >= xx) exit + end do + x(i-1) = xx + idx(i-1) = ix + endif + enddo +end subroutine psi_lisrx_up + +subroutine psi_lisrx_dw(n,x,idx) + use psb_l_sort_mod, psb_protect_name => psi_lisrx_dw + use psb_error_mod + implicit none + integer(psb_lpk_), intent(inout) :: x(:) + integer(psb_lpk_), intent(inout) :: idx(:) + integer(psb_lpk_), intent(in) :: n + integer(psb_lpk_) :: i,j,ix + integer(psb_lpk_) :: xx + + do j=n-1,1,-1 + if (x(j+1) > x(j)) then + xx = x(j) + ix = idx(j) + i=j+1 + do + x(i-1) = x(i) + idx(i-1) = idx(i) + i = i+1 + if (i>n) exit + if (x(i) <= xx) exit + end do + x(i-1) = xx + idx(i-1) = ix + endif + enddo +end subroutine psi_lisrx_dw + + +subroutine psi_lisr_up(n,x) + use psb_l_sort_mod, psb_protect_name => psi_lisr_up + use psb_error_mod + implicit none + integer(psb_lpk_), intent(inout) :: x(:) + integer(psb_lpk_), intent(in) :: n + integer(psb_lpk_) :: i,j + integer(psb_lpk_) :: xx + + do j=n-1,1,-1 + if (x(j+1) < x(j)) then + xx = x(j) + i=j+1 + do + x(i-1) = x(i) + i = i+1 + if (i>n) exit + if (x(i) >= xx) exit + end do + x(i-1) = xx + endif + enddo +end subroutine psi_lisr_up + +subroutine psi_lisr_dw(n,x) + use psb_l_sort_mod, psb_protect_name => psi_lisr_dw + use psb_error_mod + implicit none + integer(psb_lpk_), intent(inout) :: x(:) + integer(psb_lpk_), intent(in) :: n + integer(psb_lpk_) :: i,j + integer(psb_lpk_) :: xx + + do j=n-1,1,-1 + if (x(j+1) > x(j)) then + xx = x(j) + i=j+1 + do + x(i-1) = x(i) + i = i+1 + if (i>n) exit + if (x(i) <= xx) exit + end do + x(i-1) = xx + endif + enddo +end subroutine psi_lisr_dw + +subroutine psi_laisrx_up(n,x,idx) + use psb_l_sort_mod, psb_protect_name => psi_laisrx_up + use psb_error_mod + implicit none + integer(psb_lpk_), intent(inout) :: x(:) + integer(psb_lpk_), intent(inout) :: idx(:) + integer(psb_lpk_), intent(in) :: n + integer(psb_lpk_) :: i,j,ix + integer(psb_lpk_) :: xx + + do j=n-1,1,-1 + if (abs(x(j+1)) < abs(x(j))) then + xx = x(j) + ix = idx(j) + i=j+1 + do + x(i-1) = x(i) + idx(i-1) = idx(i) + i = i+1 + if (i>n) exit + if (abs(x(i)) >= abs(xx)) exit + end do + x(i-1) = xx + idx(i-1) = ix + endif + enddo +end subroutine psi_laisrx_up + +subroutine psi_laisrx_dw(n,x,idx) + use psb_l_sort_mod, psb_protect_name => psi_laisrx_dw + use psb_error_mod + implicit none + integer(psb_lpk_), intent(inout) :: x(:) + integer(psb_lpk_), intent(inout) :: idx(:) + integer(psb_lpk_), intent(in) :: n + integer(psb_lpk_) :: i,j,ix + integer(psb_lpk_) :: xx + + do j=n-1,1,-1 + if (abs(x(j+1)) > abs(x(j))) then + xx = x(j) + ix = idx(j) + i=j+1 + do + x(i-1) = x(i) + idx(i-1) = idx(i) + i = i+1 + if (i>n) exit + if (abs(x(i)) <= abs(xx)) exit + end do + x(i-1) = xx + idx(i-1) = ix + endif + enddo +end subroutine psi_laisrx_dw + +subroutine psi_laisr_up(n,x) + use psb_l_sort_mod, psb_protect_name => psi_laisr_up + use psb_error_mod + implicit none + integer(psb_lpk_), intent(inout) :: x(:) + integer(psb_lpk_), intent(in) :: n + integer(psb_lpk_) :: i,j + integer(psb_lpk_) :: xx + + do j=n-1,1,-1 + if (abs(x(j+1)) < abs(x(j))) then + xx = x(j) + i=j+1 + do + x(i-1) = x(i) + i = i+1 + if (i>n) exit + if (abs(x(i)) >= abs(xx)) exit + end do + x(i-1) = xx + endif + enddo +end subroutine psi_laisr_up + +subroutine psi_laisr_dw(n,x) + use psb_l_sort_mod, psb_protect_name => psi_laisr_dw + use psb_error_mod + implicit none + integer(psb_lpk_), intent(inout) :: x(:) + integer(psb_lpk_), intent(in) :: n + integer(psb_lpk_) :: i,j + integer(psb_lpk_) :: xx + + do j=n-1,1,-1 + if (abs(x(j+1)) > abs(x(j))) then + xx = x(j) + i=j+1 + do + x(i-1) = x(i) + i = i+1 + if (i>n) exit + if (abs(x(i)) <= abs(xx)) exit + end do + x(i-1) = xx + endif + enddo +end subroutine psi_laisr_dw + diff --git a/base/serial/sort/psb_l_msort_impl.f90 b/base/serial/sort/psb_l_msort_impl.f90 new file mode 100644 index 000000000..309176962 --- /dev/null +++ b/base/serial/sort/psb_l_msort_impl.f90 @@ -0,0 +1,713 @@ +! +! Parallel Sparse BLAS version 3.5 +! (C) Copyright 2006-2018 +! Salvatore Filippone +! Alfredo Buttari +! +! Redistribution and use in source and binary forms, with or without +! modification, are permitted provided that the following conditions +! are met: +! 1. Redistributions of source code must retain the above copyright +! notice, this list of conditions and the following disclaimer. +! 2. Redistributions in binary form must reproduce the above copyright +! notice, this list of conditions, and the following disclaimer in the +! documentation and/or other materials provided with the distribution. +! 3. The name of the PSBLAS group or the names of its contributors may +! not be used to endorse or promote products derived from this +! software without specific written permission. +! +! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +! ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED +! TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR +! PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE PSBLAS GROUP OR ITS CONTRIBUTORS +! BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +! CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +! SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +! INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +! CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +! ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +! POSSIBILITY OF SUCH DAMAGE. +! +! + ! + ! The merge-sort routines + ! References: + ! D. Knuth + ! The Art of Computer Programming, vol. 3 + ! Addison-Wesley + ! + ! Aho, Hopcroft, Ullman + ! Data Structures and Algorithms + ! Addison-Wesley + ! + logical function psb_lisaperm(n,eip) + use psb_l_sort_mod, psb_protect_name => psb_lisaperm + implicit none + + integer(psb_lpk_), intent(in) :: n + integer(psb_lpk_), intent(in) :: eip(n) + integer(psb_lpk_), allocatable :: ip(:) + integer(psb_lpk_) :: i,j,m, info + + + psb_lisaperm = .true. + if (n <= 0) return + allocate(ip(n), stat=info) + if (info /= psb_success_) return + ! + ! sanity check first + ! + do i=1, n + ip(i) = eip(i) + if ((ip(i) < 1).or.(ip(i) > n)) then + write(psb_err_unit,*) 'Out of bounds in isaperm' ,ip(i), n + psb_lisaperm = .false. + return + endif + enddo + + ! + ! now work through the cycles, by marking each successive item as negative. + ! no cycle should intersect with any other, hence the >= 1 check. + ! + do m = 1, n + i = ip(m) + if (i < 0) then + ip(m) = -i + else if (i /= m) then + j = ip(i) + ip(i) = -j + i = j + do while ((j >= 1).and.(j /= m)) + j = ip(i) + ip(i) = -j + i = j + enddo + ip(m) = abs(ip(m)) + if (j /= m) then + psb_lisaperm = .false. + goto 9999 + endif + end if + enddo +9999 continue + + return + end function psb_lisaperm + + + subroutine psb_lmsort_u(x,nout,dir) + use psb_l_sort_mod, psb_protect_name => psb_lmsort_u + use psb_error_mod + implicit none + integer(psb_lpk_), intent(inout) :: x(:) + integer(psb_lpk_), intent(out) :: nout + integer(psb_ipk_), optional, intent(in) :: dir + + integer(psb_lpk_) :: n, k + integer(psb_ipk_) :: err_act + + integer(psb_ipk_) :: ierr(5) + character(len=20) :: name + + name='psb_msort_u' + call psb_erractionsave(err_act) + + n = size(x) + + call psb_msort(x,dir=dir) + nout = min(1,n) + do k=2,n + if (x(k) /= x(nout)) then + nout = nout + 1 + x(nout) = x(k) + endif + enddo + + return + +9999 call psb_error_handler(err_act) + + return + end subroutine psb_lmsort_u + + + function psb_lbsrch(key,n,v) result(ipos) + use psb_l_sort_mod, psb_protect_name => psb_lbsrch + implicit none + integer(psb_ipk_) :: ipos, n + integer(psb_lpk_) :: key + integer(psb_lpk_) :: v(:) + + integer(psb_ipk_) :: lb, ub, m, i + + ipos = -1 + if (n<5) then + do i=1,n + if (key.eq.v(i)) then + ipos = i + return + end if + enddo + return + end if + + lb = 1 + ub = n + + do while (lb.le.ub) + m = (lb+ub)/2 + if (key.eq.v(m)) then + ipos = m + lb = ub + 1 + else if (key < v(m)) then + ub = m-1 + else + lb = m + 1 + end if + enddo + return + end function psb_lbsrch + + function psb_lssrch(key,n,v) result(ipos) + use psb_l_sort_mod, psb_protect_name => psb_lssrch + implicit none + integer(psb_ipk_) :: ipos, n + integer(psb_lpk_) :: key + integer(psb_lpk_) :: v(:) + + integer(psb_ipk_) :: i + + ipos = -1 + do i=1,n + if (key.eq.v(i)) then + ipos = i + return + end if + enddo + + return + end function psb_lssrch + + subroutine psb_lmsort(x,ix,dir,flag) + use psb_l_sort_mod, psb_protect_name => psb_lmsort + use psb_error_mod + use psb_ip_reord_mod + implicit none + integer(psb_lpk_), intent(inout) :: x(:) + integer(psb_ipk_), optional, intent(in) :: dir, flag + integer(psb_lpk_), optional, intent(inout) :: ix(:) + + integer(psb_ipk_) :: dir_, flag_, n, err_act + + integer(psb_lpk_), allocatable :: iaux(:) + integer(psb_ipk_) :: iret, info, i + integer(psb_ipk_) :: ierr(5) + character(len=20) :: name + + name='psb_lmsort' + call psb_erractionsave(err_act) + + if (present(dir)) then + dir_ = dir + else + dir_= psb_sort_up_ + end if + select case(dir_) + case( psb_sort_up_, psb_sort_down_, psb_asort_up_, psb_asort_down_) + ! OK keep going + case default + ierr(1) = 3; ierr(2) = dir_; + call psb_errpush(psb_err_input_value_invalid_i_,name,i_err=ierr) + goto 9999 + end select + + n = size(x) + + if (present(ix)) then + if (size(ix) < n) then + ierr(1) = 2; ierr(2) = size(ix); + call psb_errpush(psb_err_input_asize_invalid_i_,name,i_err=ierr) + goto 9999 + end if + if (present(flag)) then + flag_ = flag + else + flag_ = psb_sort_ovw_idx_ + end if + select case(flag_) + case(psb_sort_ovw_idx_) + do i=1,n + ix(i) = i + end do + case (psb_sort_keep_idx_) + ! OK keep going + case default + ierr(1) = 4; ierr(2) = flag_; + call psb_errpush(psb_err_input_value_invalid_i_,name,i_err=ierr) + goto 9999 + end select + + end if + + allocate(iaux(0:n+1),stat=info) + if (info /= psb_success_) then + call psb_errpush(psb_err_alloc_dealloc_,r_name='psb_l_msort') + goto 9999 + endif + + select case(dir_) + case (psb_sort_up_) + call psi_l_msort_up(n,x,iaux,iret) + case (psb_sort_down_) + call psi_l_msort_dw(n,x,iaux,iret) + case (psb_asort_up_) + call psi_l_amsort_up(n,x,iaux,iret) + case (psb_asort_down_) + call psi_l_amsort_dw(n,x,iaux,iret) + end select + ! + ! Do the actual reordering, since the inner routines + ! only provide linked pointers. + ! + if (iret == 0 ) then + if (present(ix)) then + call psb_ip_reord(n,x,ix,iaux) + else + call psb_ip_reord(n,x,iaux) + end if + end if + + + return + +9999 call psb_error_handler(err_act) + + return + + + end subroutine psb_lmsort + + subroutine psi_l_msort_up(n,k,l,iret) + use psb_const_mod + implicit none + integer(psb_ipk_) :: n, iret + integer(psb_lpk_) :: k(n) + integer(psb_lpk_) :: l(0:n+1) + ! + integer(psb_lpk_) :: p,q,s,t + ! .. + iret = 0 + ! first step: we are preparing ordered sublists, exploiting + ! what order was already in the input data; negative links + ! mark the end of the sublists + l(0) = 1 + t = n + 1 + do p = 1,n - 1 + if (k(p) <= k(p+1)) then + l(p) = p + 1 + else + l(t) = - (p+1) + t = p + end if + end do + l(t) = 0 + l(n) = 0 + ! see if the input was already sorted + if (l(n+1) == 0) then + iret = 1 + return + else + l(n+1) = abs(l(n+1)) + end if + + mergepass: do + ! otherwise, begin a pass through the list. + ! throughout all the subroutine we have: + ! p, q: pointing to the sublists being merged + ! s: pointing to the most recently processed record + ! t: pointing to the end of previously completed sublist + s = 0 + t = n + 1 + p = l(s) + q = l(t) + if (q == 0) exit mergepass + + outer: do + + if (k(p) > k(q)) then + + l(s) = sign(q,l(s)) + s = q + q = l(q) + if (q > 0) then + do + if (k(p) <= k(q)) cycle outer + s = q + q = l(q) + if (q <= 0) exit + end do + end if + l(s) = p + s = t + do + t = p + p = l(p) + if (p <= 0) exit + end do + + else + + l(s) = sign(p,l(s)) + s = p + p = l(p) + if (p>0) then + do + if (k(p) > k(q)) cycle outer + s = p + p = l(p) + if (p <= 0) exit + end do + end if + ! otherwise, one sublist ended, and we append to it the rest + ! of the other one. + l(s) = q + s = t + do + t = q + q = l(q) + if (q <= 0) exit + end do + end if + + p = -p + q = -q + if (q == 0) then + l(s) = sign(p,l(s)) + l(t) = 0 + exit outer + end if + end do outer + end do mergepass + + end subroutine psi_l_msort_up + + subroutine psi_l_msort_dw(n,k,l,iret) + use psb_const_mod + implicit none + integer(psb_ipk_) :: n, iret + integer(psb_lpk_) :: k(n) + integer(psb_lpk_) :: l(0:n+1) + ! + integer(psb_lpk_) :: p,q,s,t + ! .. + iret = 0 + ! first step: we are preparing ordered sublists, exploiting + ! what order was already in the input data; negative links + ! mark the end of the sublists + l(0) = 1 + t = n + 1 + do p = 1,n - 1 + if (k(p) >= k(p+1)) then + l(p) = p + 1 + else + l(t) = - (p+1) + t = p + end if + end do + l(t) = 0 + l(n) = 0 + ! see if the input was already sorted + if (l(n+1) == 0) then + iret = 1 + return + else + l(n+1) = abs(l(n+1)) + end if + + mergepass: do + ! otherwise, begin a pass through the list. + ! throughout all the subroutine we have: + ! p, q: pointing to the sublists being merged + ! s: pointing to the most recently processed record + ! t: pointing to the end of previously completed sublist + s = 0 + t = n + 1 + p = l(s) + q = l(t) + if (q == 0) exit mergepass + + outer: do + + if (k(p) < k(q)) then + + l(s) = sign(q,l(s)) + s = q + q = l(q) + if (q > 0) then + do + if (k(p) >= k(q)) cycle outer + s = q + q = l(q) + if (q <= 0) exit + end do + end if + l(s) = p + s = t + do + t = p + p = l(p) + if (p <= 0) exit + end do + + else + + l(s) = sign(p,l(s)) + s = p + p = l(p) + if (p>0) then + do + if (k(p) < k(q)) cycle outer + s = p + p = l(p) + if (p <= 0) exit + end do + end if + ! otherwise, one sublist ended, and we append to it the rest + ! of the other one. + l(s) = q + s = t + do + t = q + q = l(q) + if (q <= 0) exit + end do + end if + + p = -p + q = -q + if (q == 0) then + l(s) = sign(p,l(s)) + l(t) = 0 + exit outer + end if + end do outer + end do mergepass + + end subroutine psi_l_msort_dw + + subroutine psi_l_amsort_up(n,k,l,iret) + use psb_const_mod + implicit none + integer(psb_ipk_) :: n, iret + integer(psb_lpk_) :: k(n) + integer(psb_lpk_) :: l(0:n+1) + ! + integer(psb_lpk_) :: p,q,s,t + ! .. + iret = 0 + ! first step: we are preparing ordered sublists, exploiting + ! what order was already in the input data; negative links + ! mark the end of the sublists + l(0) = 1 + t = n + 1 + do p = 1,n - 1 + if (abs(k(p)) <= abs(k(p+1))) then + l(p) = p + 1 + else + l(t) = - (p+1) + t = p + end if + end do + l(t) = 0 + l(n) = 0 + ! see if the input was already sorted + if (l(n+1) == 0) then + iret = 1 + return + else + l(n+1) = abs(l(n+1)) + end if + + mergepass: do + ! otherwise, begin a pass through the list. + ! throughout all the subroutine we have: + ! p, q: pointing to the sublists being merged + ! s: pointing to the most recently processed record + ! t: pointing to the end of previously completed sublist + s = 0 + t = n + 1 + p = l(s) + q = l(t) + if (q == 0) exit mergepass + + outer: do + + if (abs(k(p)) > abs(k(q))) then + + l(s) = sign(q,l(s)) + s = q + q = l(q) + if (q > 0) then + do + if (abs(k(p)) <= abs(k(q))) cycle outer + s = q + q = l(q) + if (q <= 0) exit + end do + end if + l(s) = p + s = t + do + t = p + p = l(p) + if (p <= 0) exit + end do + + else + + l(s) = sign(p,l(s)) + s = p + p = l(p) + if (p>0) then + do + if (abs(k(p)) > abs(k(q))) cycle outer + s = p + p = l(p) + if (p <= 0) exit + end do + end if + ! otherwise, one sublist ended, and we append to it the rest + ! of the other one. + l(s) = q + s = t + do + t = q + q = l(q) + if (q <= 0) exit + end do + end if + + p = -p + q = -q + if (q == 0) then + l(s) = sign(p,l(s)) + l(t) = 0 + exit outer + end if + end do outer + end do mergepass + + end subroutine psi_l_amsort_up + + subroutine psi_l_amsort_dw(n,k,l,iret) + use psb_const_mod + implicit none + integer(psb_ipk_) :: n, iret + integer(psb_lpk_) :: k(n) + integer(psb_lpk_) :: l(0:n+1) + ! + integer(psb_lpk_) :: p,q,s,t + ! .. + iret = 0 + ! first step: we are preparing ordered sublists, exploiting + ! what order was already in the input data; negative links + ! mark the end of the sublists + l(0) = 1 + t = n + 1 + do p = 1,n - 1 + if (abs(k(p)) >= abs(k(p+1))) then + l(p) = p + 1 + else + l(t) = - (p+1) + t = p + end if + end do + l(t) = 0 + l(n) = 0 + ! see if the input was already sorted + if (l(n+1) == 0) then + iret = 1 + return + else + l(n+1) = abs(l(n+1)) + end if + + mergepass: do + ! otherwise, begin a pass through the list. + ! throughout all the subroutine we have: + ! p, q: pointing to the sublists being merged + ! s: pointing to the most recently processed record + ! t: pointing to the end of previously completed sublist + s = 0 + t = n + 1 + p = l(s) + q = l(t) + if (q == 0) exit mergepass + + outer: do + + if (abs(k(p)) < abs(k(q))) then + + l(s) = sign(q,l(s)) + s = q + q = l(q) + if (q > 0) then + do + if (abs(k(p)) >= abs(k(q))) cycle outer + s = q + q = l(q) + if (q <= 0) exit + end do + end if + l(s) = p + s = t + do + t = p + p = l(p) + if (p <= 0) exit + end do + + else + + l(s) = sign(p,l(s)) + s = p + p = l(p) + if (p>0) then + do + if (abs(k(p)) < abs(k(q))) cycle outer + s = p + p = l(p) + if (p <= 0) exit + end do + end if + ! otherwise, one sublist ended, and we append to it the rest + ! of the other one. + l(s) = q + s = t + do + t = q + q = l(q) + if (q <= 0) exit + end do + end if + + p = -p + q = -q + if (q == 0) then + l(s) = sign(p,l(s)) + l(t) = 0 + exit outer + end if + end do outer + end do mergepass + + end subroutine psi_l_amsort_dw + + + + + + + + diff --git a/base/serial/sort/psb_l_qsort_impl.f90 b/base/serial/sort/psb_l_qsort_impl.f90 new file mode 100644 index 000000000..4c1be19b1 --- /dev/null +++ b/base/serial/sort/psb_l_qsort_impl.f90 @@ -0,0 +1,1318 @@ +! +! Parallel Sparse BLAS version 3.5 +! (C) Copyright 2006-2018 +! Salvatore Filippone +! Alfredo Buttari +! +! Redistribution and use in source and binary forms, with or without +! modification, are permitted provided that the following conditions +! are met: +! 1. Redistributions of source code must retain the above copyright +! notice, this list of conditions and the following disclaimer. +! 2. Redistributions in binary form must reproduce the above copyright +! notice, this list of conditions, and the following disclaimer in the +! documentation and/or other materials provided with the distribution. +! 3. The name of the PSBLAS group or the names of its contributors may +! not be used to endorse or promote products derived from this +! software without specific written permission. +! +! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +! ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED +! TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR +! PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE PSBLAS GROUP OR ITS CONTRIBUTORS +! BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +! CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +! SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +! INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +! CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +! ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +! POSSIBILITY OF SUCH DAMAGE. +! +! +! +! The quicksort routines +! References: +! D. Knuth +! The Art of Computer Programming, vol. 3 +! Addison-Wesley +! +! Aho, Hopcroft, Ullman +! Data Structures and Algorithms +! Addison-Wesley +! +subroutine psb_lqsort(x,ix,dir,flag) + use psb_l_sort_mod, psb_protect_name => psb_lqsort + use psb_error_mod + implicit none + integer(psb_lpk_), intent(inout) :: x(:) + integer(psb_ipk_), optional, intent(in) :: dir, flag + integer(psb_lpk_), optional, intent(inout) :: ix(:) + + integer(psb_ipk_) :: dir_, flag_, err_act, i + integer(psb_lpk_) :: n + integer(psb_ipk_) :: ierr(5) + character(len=20) :: name + + name='psb_lqsort' + call psb_erractionsave(err_act) + + if (present(flag)) then + flag_ = flag + else + flag_ = psb_sort_ovw_idx_ + end if + select case(flag_) + case( psb_sort_ovw_idx_, psb_sort_keep_idx_) + ! OK keep going + case default + ierr(1) = 4; ierr(2) = flag_; + call psb_errpush(psb_err_input_value_invalid_i_,name,i_err=ierr) + goto 9999 + end select + + if (present(dir)) then + dir_ = dir + else + dir_= psb_sort_up_ + end if + + n = size(x) + + if (present(ix)) then + if (size(ix) < n) then + ierr(1) = 2; ierr(2) = size(ix); + call psb_errpush(psb_err_input_asize_invalid_i_,name,i_err=ierr) + goto 9999 + end if + if (flag_==psb_sort_ovw_idx_) then + do i=1,n + ix(i) = i + end do + end if + + select case(dir_) + case (psb_sort_up_) + call psi_lqsrx_up(n,x,ix) + case (psb_sort_down_) + call psi_lqsrx_dw(n,x,ix) + case (psb_asort_up_) + call psi_laqsrx_up(n,x,ix) + case (psb_asort_down_) + call psi_laqsrx_dw(n,x,ix) + case default + ierr(1) = 3; ierr(2) = dir_; + call psb_errpush(psb_err_input_value_invalid_i_,name,i_err=ierr) + goto 9999 + end select + else + select case(dir_) + case (psb_sort_up_) + call psi_lqsr_up(n,x) + case (psb_sort_down_) + call psi_lqsr_dw(n,x) + case (psb_asort_up_) + call psi_laqsr_up(n,x) + case (psb_asort_down_) + call psi_laqsr_dw(n,x) + case default + ierr(1) = 3; ierr(2) = dir_; + call psb_errpush(psb_err_input_value_invalid_i_,name,i_err=ierr) + goto 9999 + end select + + end if + + return + +9999 call psb_error_handler(err_act) + + return +end subroutine psb_lqsort + +subroutine psi_lqsrx_up(n,x,idx) + use psb_l_sort_mod, psb_protect_name => psi_lqsrx_up + use psb_error_mod + implicit none + + integer(psb_lpk_), intent(inout) :: x(:) + integer(psb_lpk_), intent(inout) :: idx(:) + integer(psb_lpk_), intent(in) :: n + ! .. Local Scalars .. + integer(psb_lpk_) :: piv, xk, xt + integer(psb_ipk_) :: i, j, ilx, iux, istp, lpiv + integer(psb_lpk_) :: n1, n2 + integer(psb_lpk_) :: ixt + integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=72 + integer(psb_ipk_) :: istack(nparms,maxstack) + + if (n > ithrs) then + ! + ! Init stack pointer + ! + istp = 1 + istack(1,istp) = 1 + istack(2,istp) = n + + do + if (istp <= 0) exit + ilx = istack(1,istp) + iux = istack(2,istp) + istp = istp - 1 + ! + ! Choose a pivot with median-of-three heuristics, leave it + ! in the LPIV location + ! + i = ilx + j = iux + lpiv = (i+j)/2 + piv = x(lpiv) + if (piv < x(i)) then + xt = x(i) + ixt = idx(i) + x(i) = x(lpiv) + idx(i) = idx(lpiv) + x(lpiv) = xt + idx(lpiv) = ixt + piv = x(lpiv) + endif + if (piv > x(j)) then + xt = x(j) + ixt = idx(j) + x(j) = x(lpiv) + idx(j) = idx(lpiv) + x(lpiv) = xt + idx(lpiv) = ixt + piv = x(lpiv) + endif + if (piv < x(i)) then + xt = x(i) + ixt = idx(i) + x(i) = x(lpiv) + idx(i) = idx(lpiv) + x(lpiv) = xt + idx(lpiv) = ixt + piv = x(lpiv) + endif + ! + ! now piv is correct; place it into first location + xt = x(i) + ixt = idx(i) + x(i) = x(lpiv) + idx(i) = idx(lpiv) + x(lpiv) = xt + idx(lpiv) = ixt + piv = x(lpiv) + + i = ilx - 1 + j = iux + 1 + + outer_up: do + in_up1: do + i = i + 1 + xk = x(i) + if (xk >= piv) exit in_up1 + end do in_up1 + ! + ! Ensure finite termination for next loop + ! + xt = xk + x(i) = piv + in_up2:do + j = j - 1 + xk = x(j) + if (xk <= piv) exit in_up2 + end do in_up2 + x(i) = xt + + if (j > i) then + xt = x(i) + ixt = idx(i) + x(i) = x(j) + idx(i) = idx(j) + x(j) = xt + idx(j) = ixt + else + exit outer_up + end if + end do outer_up + if (i == ilx) then + if (x(i) /= piv) then + call psb_errpush(psb_err_internal_error_,& + & r_name='psi_lqsrx',a_err='impossible pivot condition') + call psb_error() + endif + i = i + 1 + endif + + n1 = (i-1)-ilx+1 + n2 = iux-(i)+1 + if (n1 > n2) then + if (n1 > ithrs) then + istp = istp + 1 + istack(1,istp) = ilx + istack(2,istp) = i-1 + else + call psi_lisrx_up(n1,x(ilx:i-1),idx(ilx:i-1)) + endif + if (n2 > ithrs) then + istp = istp + 1 + istack(1,istp) = i + istack(2,istp) = iux + else + call psi_lisrx_up(n2,x(i:iux),idx(i:iux)) + endif + else + if (n2 > ithrs) then + istp = istp + 1 + istack(1,istp) = i + istack(2,istp) = iux + else + call psi_lisrx_up(n2,x(i:iux),idx(i:iux)) + endif + if (n1 > ithrs) then + istp = istp + 1 + istack(1,istp) = ilx + istack(2,istp) = i-1 + else + call psi_lisrx_up(n1,x(ilx:i-1),idx(ilx:i-1)) + endif + endif + enddo + else + call psi_lisrx_up(n,x,idx) + endif +end subroutine psi_lqsrx_up + +subroutine psi_lqsrx_dw(n,x,idx) + use psb_l_sort_mod, psb_protect_name => psi_lqsrx_dw + use psb_error_mod + implicit none + + integer(psb_lpk_), intent(inout) :: x(:) + integer(psb_lpk_), intent(inout) :: idx(:) + integer(psb_lpk_), intent(in) :: n + ! .. Local Scalars .. + integer(psb_lpk_) :: piv, xk, xt + integer(psb_ipk_) :: i, j, ilx, iux, istp, lpiv + integer(psb_lpk_) :: n1, n2 + integer(psb_lpk_) :: ixt + + integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=72 + integer(psb_ipk_) :: istack(nparms,maxstack) + + if (n > ithrs) then + ! + ! Init stack pointer + ! + istp = 1 + istack(1,istp) = 1 + istack(2,istp) = n + + do + if (istp <= 0) exit + ilx = istack(1,istp) + iux = istack(2,istp) + istp = istp - 1 + ! + ! Choose a pivot with median-of-three heuristics, leave it + ! in the LPIV location + ! + i = ilx + j = iux + lpiv = (i+j)/2 + piv = x(lpiv) + if (piv > x(i)) then + xt = x(i) + ixt = idx(i) + x(i) = x(lpiv) + idx(i) = idx(lpiv) + x(lpiv) = xt + idx(lpiv) = ixt + piv = x(lpiv) + endif + if (piv < x(j)) then + xt = x(j) + ixt = idx(j) + x(j) = x(lpiv) + idx(j) = idx(lpiv) + x(lpiv) = xt + idx(lpiv) = ixt + piv = x(lpiv) + endif + if (piv > x(i)) then + xt = x(i) + ixt = idx(i) + x(i) = x(lpiv) + idx(i) = idx(lpiv) + x(lpiv) = xt + idx(lpiv) = ixt + piv = x(lpiv) + endif + ! + ! now piv is correct; place it into first location + xt = x(i) + ixt = idx(i) + x(i) = x(lpiv) + idx(i) = idx(lpiv) + x(lpiv) = xt + idx(lpiv) = ixt + piv = x(lpiv) + + i = ilx - 1 + j = iux + 1 + + outer_dw: do + in_dw1: do + i = i + 1 + xk = x(i) + if (xk <= piv) exit in_dw1 + end do in_dw1 + ! + ! Ensure finite termination for next loop + ! + xt = xk + x(i) = piv + in_dw2:do + j = j - 1 + xk = x(j) + if (xk >= piv) exit in_dw2 + end do in_dw2 + x(i) = xt + + if (j > i) then + xt = x(i) + ixt = idx(i) + x(i) = x(j) + idx(i) = idx(j) + x(j) = xt + idx(j) = ixt + else + exit outer_dw + end if + end do outer_dw + if (i == ilx) then + if (x(i) /= piv) then + call psb_errpush(psb_err_internal_error_,& + & r_name='psi_lqsrx',a_err='impossible pivot condition') + call psb_error() + endif + i = i + 1 + endif + + n1 = (i-1)-ilx+1 + n2 = iux-(i)+1 + if (n1 > n2) then + if (n1 > ithrs) then + istp = istp + 1 + istack(1,istp) = ilx + istack(2,istp) = i-1 + else + call psi_lisrx_dw(n1,x(ilx:i-1),idx(ilx:i-1)) + endif + if (n2 > ithrs) then + istp = istp + 1 + istack(1,istp) = i + istack(2,istp) = iux + else + call psi_lisrx_dw(n2,x(i:iux),idx(i:iux)) + endif + else + if (n2 > ithrs) then + istp = istp + 1 + istack(1,istp) = i + istack(2,istp) = iux + else + call psi_lisrx_dw(n2,x(i:iux),idx(i:iux)) + endif + if (n1 > ithrs) then + istp = istp + 1 + istack(1,istp) = ilx + istack(2,istp) = i-1 + else + call psi_lisrx_dw(n1,x(ilx:i-1),idx(ilx:i-1)) + endif + endif + enddo + else + call psi_lisrx_dw(n,x,idx) + endif + +end subroutine psi_lqsrx_dw + +subroutine psi_lqsr_up(n,x) + use psb_l_sort_mod, psb_protect_name => psi_lqsr_up + use psb_error_mod + implicit none + + integer(psb_lpk_), intent(inout) :: x(:) + integer(psb_lpk_), intent(in) :: n + ! .. + ! .. Local Scalars .. + integer(psb_lpk_) :: piv, xt, xk + integer(psb_ipk_) :: i, j, ilx, iux, istp, lpiv + integer(psb_lpk_) :: n1, n2 + + integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=72 + integer(psb_ipk_) :: istack(nparms,maxstack) + + + if (n > ithrs) then + ! + ! Init stack pointer + ! + istp = 1 + istack(1,istp) = 1 + istack(2,istp) = n + + do + if (istp <= 0) exit + ilx = istack(1,istp) + iux = istack(2,istp) + istp = istp - 1 + ! + ! Choose a pivot with median-of-three heuristics, leave it + ! in the LPIV location + ! + i = ilx + j = iux + lpiv = (i+j)/2 + piv = x(lpiv) + if (piv < x(i)) then + xt = x(i) + x(i) = x(lpiv) + x(lpiv) = xt + piv = x(lpiv) + endif + if (piv > x(j)) then + xt = x(j) + x(j) = x(lpiv) + x(lpiv) = xt + piv = x(lpiv) + endif + if (piv < x(i)) then + xt = x(i) + x(i) = x(lpiv) + x(lpiv) = xt + piv = x(lpiv) + endif + ! + ! now piv is correct; place it into first location + + xt = x(i) + x(i) = x(lpiv) + x(lpiv) = xt + + i = ilx - 1 + j = iux + 1 + + outer_up: do + in_up1: do + i = i + 1 + xk = x(i) + if (xk >= piv) exit in_up1 + end do in_up1 + ! + ! Ensure finite termination for next loop + ! + xt = xk + x(i) = piv + in_up2:do + j = j - 1 + xk = x(j) + if (xk <= piv) exit in_up2 + end do in_up2 + x(i) = xt + + if (j > i) then + xt = x(i) + x(i) = x(j) + x(j) = xt + else + exit outer_up + end if + end do outer_up + if (i == ilx) then + if (x(i) /= piv) then + call psb_errpush(psb_err_internal_error_,& + & r_name='psi_lqsr',a_err='impossible pivot condition') + call psb_error() + endif + i = i + 1 + endif + + n1 = (i-1)-ilx+1 + n2 = iux-(i)+1 + if (n1 > n2) then + if (n1 > ithrs) then + istp = istp + 1 + istack(1,istp) = ilx + istack(2,istp) = i-1 + else + call psi_lisr_up(n1,x(ilx:i-1)) + endif + if (n2 > ithrs) then + istp = istp + 1 + istack(1,istp) = i + istack(2,istp) = iux + else + call psi_lisr_up(n2,x(i:iux)) + endif + else + if (n2 > ithrs) then + istp = istp + 1 + istack(1,istp) = i + istack(2,istp) = iux + else + call psi_lisr_up(n2,x(i:iux)) + endif + if (n1 > ithrs) then + istp = istp + 1 + istack(1,istp) = ilx + istack(2,istp) = i-1 + else + call psi_lisr_up(n1,x(ilx:i-1)) + endif + endif + enddo + else + call psi_lisr_up(n,x) + endif + +end subroutine psi_lqsr_up + +subroutine psi_lqsr_dw(n,x) + use psb_l_sort_mod, psb_protect_name => psi_lqsr_dw + use psb_error_mod + implicit none + + integer(psb_lpk_), intent(inout) :: x(:) + integer(psb_lpk_), intent(in) :: n + ! .. + ! .. Local Scalars .. + integer(psb_lpk_) :: piv, xt, xk + integer(psb_ipk_) :: i, j, ilx, iux, istp, lpiv + integer(psb_lpk_) :: n1, n2 + + integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=72 + integer(psb_ipk_) :: istack(nparms,maxstack) + + + if (n > ithrs) then + ! + ! Init stack pointer + ! + istp = 1 + istack(1,istp) = 1 + istack(2,istp) = n + + do + if (istp <= 0) exit + ilx = istack(1,istp) + iux = istack(2,istp) + istp = istp - 1 + ! + ! Choose a pivot with median-of-three heuristics, leave it + ! in the LPIV location + ! + i = ilx + j = iux + lpiv = (i+j)/2 + piv = x(lpiv) + if (piv > x(i)) then + xt = x(i) + x(i) = x(lpiv) + x(lpiv) = xt + piv = x(lpiv) + endif + if (piv < x(j)) then + xt = x(j) + x(j) = x(lpiv) + x(lpiv) = xt + piv = x(lpiv) + endif + if (piv > x(i)) then + xt = x(i) + x(i) = x(lpiv) + x(lpiv) = xt + piv = x(lpiv) + endif + ! + ! now piv is correct; place it into first location + + xt = x(i) + x(i) = x(lpiv) + x(lpiv) = xt + + i = ilx - 1 + j = iux + 1 + + outer_dw: do + in_dw1: do + i = i + 1 + xk = x(i) + if (xk <= piv) exit in_dw1 + end do in_dw1 + ! + ! Ensure finite termination for next loop + ! + xt = xk + x(i) = piv + in_dw2:do + j = j - 1 + xk = x(j) + if (xk >= piv) exit in_dw2 + end do in_dw2 + x(i) = xt + + if (j > i) then + xt = x(i) + x(i) = x(j) + x(j) = xt + else + exit outer_dw + end if + end do outer_dw + if (i == ilx) then + if (x(i) /= piv) then + call psb_errpush(psb_err_internal_error_, & + & r_name='psi_lqsr',a_err='impossible pivot condition') + call psb_error() + endif + i = i + 1 + endif + + n1 = (i-1)-ilx+1 + n2 = iux-(i)+1 + if (n1 > n2) then + if (n1 > ithrs) then + istp = istp + 1 + istack(1,istp) = ilx + istack(2,istp) = i-1 + else + call psi_lisr_dw(n1,x(ilx:i-1)) + endif + if (n2 > ithrs) then + istp = istp + 1 + istack(1,istp) = i + istack(2,istp) = iux + else + call psi_lisr_dw(n2,x(i:iux)) + endif + else + if (n2 > ithrs) then + istp = istp + 1 + istack(1,istp) = i + istack(2,istp) = iux + else + call psi_lisr_dw(n2,x(i:iux)) + endif + if (n1 > ithrs) then + istp = istp + 1 + istack(1,istp) = ilx + istack(2,istp) = i-1 + else + call psi_lisr_dw(n1,x(ilx:i-1)) + endif + endif + enddo + else + call psi_lisr_dw(n,x) + endif + +end subroutine psi_lqsr_dw + +subroutine psi_laqsrx_up(n,x,idx) + use psb_l_sort_mod, psb_protect_name => psi_laqsrx_up + use psb_error_mod + implicit none + + integer(psb_lpk_), intent(inout) :: x(:) + integer(psb_lpk_), intent(inout) :: idx(:) + integer(psb_lpk_), intent(in) :: n + ! .. Local Scalars .. + integer(psb_lpk_) :: piv, xk + integer(psb_lpk_) :: xt + integer(psb_ipk_) :: i, j, ilx, iux, istp, lpiv + integer(psb_lpk_) :: n1, n2 + integer(psb_lpk_) :: ixt + + integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=72 + integer(psb_ipk_) :: istack(nparms,maxstack) + + if (n > ithrs) then + ! + ! Init stack pointer + ! + istp = 1 + istack(1,istp) = 1 + istack(2,istp) = n + + do + if (istp <= 0) exit + ilx = istack(1,istp) + iux = istack(2,istp) + istp = istp - 1 + ! + ! Choose a pivot with median-of-three heuristics, leave it + ! in the LPIV location + ! + i = ilx + j = iux + lpiv = (i+j)/2 + piv = abs(x(lpiv)) + if (piv < abs(x(i))) then + xt = x(i) + ixt = idx(i) + x(i) = x(lpiv) + idx(i) = idx(lpiv) + x(lpiv) = xt + idx(lpiv) = ixt + piv = abs(x(lpiv)) + endif + if (piv > abs(x(j))) then + xt = x(j) + ixt = idx(j) + x(j) = x(lpiv) + idx(j) = idx(lpiv) + x(lpiv) = xt + idx(lpiv) = ixt + piv = abs(x(lpiv)) + endif + if (piv < abs(x(i))) then + xt = x(i) + ixt = idx(i) + x(i) = x(lpiv) + idx(i) = idx(lpiv) + x(lpiv) = xt + idx(lpiv) = ixt + piv = abs(x(lpiv)) + endif + ! + ! now piv is correct; place it into first location + xt = x(i) + ixt = idx(i) + x(i) = x(lpiv) + idx(i) = idx(lpiv) + x(lpiv) = xt + idx(lpiv) = ixt + + i = ilx - 1 + j = iux + 1 + + outer_up: do + in_up1: do + i = i + 1 + xk = abs(x(i)) + if (xk >= piv) exit in_up1 + end do in_up1 + ! + ! Ensure finite termination for next loop + ! + xt = x(i) + x(i) = piv + in_up2:do + j = j - 1 + xk = abs(x(j)) + if (xk <= piv) exit in_up2 + end do in_up2 + x(i) = xt + + if (j > i) then + xt = x(i) + ixt = idx(i) + x(i) = x(j) + idx(i) = idx(j) + x(j) = xt + idx(j) = ixt + else + exit outer_up + end if + end do outer_up + if (i == ilx) then + if (x(i) /= piv) then + call psb_errpush(psb_err_internal_error_, & + & r_name='psi_laqsrx',a_err='impossible pivot condition') + call psb_error() + endif + i = i + 1 + endif + + n1 = (i-1)-ilx+1 + n2 = iux-(i)+1 + if (n1 > n2) then + if (n1 > ithrs) then + istp = istp + 1 + istack(1,istp) = ilx + istack(2,istp) = i-1 + else + call psi_laisrx_up(n1,x(ilx:i-1),idx(ilx:i-1)) + endif + if (n2 > ithrs) then + istp = istp + 1 + istack(1,istp) = i + istack(2,istp) = iux + else + call psi_laisrx_up(n2,x(i:iux),idx(i:iux)) + endif + else + if (n2 > ithrs) then + istp = istp + 1 + istack(1,istp) = i + istack(2,istp) = iux + else + call psi_laisrx_up(n2,x(i:iux),idx(i:iux)) + endif + if (n1 > ithrs) then + istp = istp + 1 + istack(1,istp) = ilx + istack(2,istp) = i-1 + else + call psi_laisrx_up(n1,x(ilx:i-1),idx(ilx:i-1)) + endif + endif + enddo + else + call psi_laisrx_up(n,x,idx) + endif + + +end subroutine psi_laqsrx_up + +subroutine psi_laqsrx_dw(n,x,idx) + use psb_l_sort_mod, psb_protect_name => psi_laqsrx_dw + use psb_error_mod + implicit none + + integer(psb_lpk_), intent(inout) :: x(:) + integer(psb_lpk_), intent(inout) :: idx(:) + integer(psb_lpk_), intent(in) :: n + ! .. Local Scalars .. + integer(psb_lpk_) :: piv, xk + integer(psb_lpk_) :: xt + integer(psb_ipk_) :: i, j, ilx, iux, istp, lpiv + integer(psb_lpk_) :: n1, n2 + integer(psb_lpk_) :: ixt + + integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=72 + integer(psb_ipk_) :: istack(nparms,maxstack) + if (n > ithrs) then + ! + ! Init stack pointer + ! + istp = 1 + istack(1,istp) = 1 + istack(2,istp) = n + + do + if (istp <= 0) exit + ilx = istack(1,istp) + iux = istack(2,istp) + istp = istp - 1 + ! + ! Choose a pivot with median-of-three heuristics, leave it + ! in the LPIV location + ! + i = ilx + j = iux + lpiv = (i+j)/2 + piv = abs(x(lpiv)) + if (piv > abs(x(i))) then + xt = x(i) + ixt = idx(i) + x(i) = x(lpiv) + idx(i) = idx(lpiv) + x(lpiv) = xt + idx(lpiv) = ixt + piv = abs(x(lpiv)) + endif + if (piv < abs(x(j))) then + xt = x(j) + ixt = idx(j) + x(j) = x(lpiv) + idx(j) = idx(lpiv) + x(lpiv) = xt + idx(lpiv) = ixt + piv = abs(x(lpiv)) + endif + if (piv > abs(x(i))) then + xt = x(i) + ixt = idx(i) + x(i) = x(lpiv) + idx(i) = idx(lpiv) + x(lpiv) = xt + idx(lpiv) = ixt + piv = abs(x(lpiv)) + endif + ! + ! now piv is correct; place it into first location + xt = x(i) + ixt = idx(i) + x(i) = x(lpiv) + idx(i) = idx(lpiv) + x(lpiv) = xt + idx(lpiv) = ixt + + i = ilx - 1 + j = iux + 1 + + outer_dw: do + in_dw1: do + i = i + 1 + xk = abs(x(i)) + if (xk <= piv) exit in_dw1 + end do in_dw1 + ! + ! Ensure finite termination for next loop + ! + xt = x(i) + x(i) = piv + in_dw2:do + j = j - 1 + xk = abs(x(j)) + if (xk >= piv) exit in_dw2 + end do in_dw2 + x(i) = xt + + if (j > i) then + xt = x(i) + ixt = idx(i) + x(i) = x(j) + idx(i) = idx(j) + x(j) = xt + idx(j) = ixt + else + exit outer_dw + end if + end do outer_dw + if (i == ilx) then + if (x(i) /= piv) then + call psb_errpush(psb_err_internal_error_,& + & r_name='psi_laqsrx',a_err='impossible pivot condition') + call psb_error() + endif + i = i + 1 + endif + + n1 = (i-1)-ilx+1 + n2 = iux-(i)+1 + if (n1 > n2) then + if (n1 > ithrs) then + istp = istp + 1 + istack(1,istp) = ilx + istack(2,istp) = i-1 + else + call psi_laisrx_dw(n1,x(ilx:i-1),idx(ilx:i-1)) + endif + if (n2 > ithrs) then + istp = istp + 1 + istack(1,istp) = i + istack(2,istp) = iux + else + call psi_laisrx_dw(n2,x(i:iux),idx(i:iux)) + endif + else + if (n2 > ithrs) then + istp = istp + 1 + istack(1,istp) = i + istack(2,istp) = iux + else + call psi_laisrx_dw(n2,x(i:iux),idx(i:iux)) + endif + if (n1 > ithrs) then + istp = istp + 1 + istack(1,istp) = ilx + istack(2,istp) = i-1 + else + call psi_laisrx_dw(n1,x(ilx:i-1),idx(ilx:i-1)) + endif + endif + enddo + else + call psi_laisrx_dw(n,x,idx) + endif + +end subroutine psi_laqsrx_dw + +subroutine psi_laqsr_up(n,x) + use psb_l_sort_mod, psb_protect_name => psi_laqsr_up + use psb_error_mod + implicit none + + integer(psb_lpk_), intent(inout) :: x(:) + integer(psb_lpk_), intent(in) :: n + ! .. Local Scalars .. + integer(psb_lpk_) :: piv, xk + integer(psb_lpk_) :: xt + integer(psb_ipk_) :: i, j, ilx, iux, istp, lpiv + integer(psb_lpk_) :: n1, n2 + integer(psb_lpk_) :: ixt + + integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=72 + integer(psb_ipk_) :: istack(nparms,maxstack) + + if (n > ithrs) then + ! + ! Init stack pointer + ! + istp = 1 + istack(1,istp) = 1 + istack(2,istp) = n + + do + if (istp <= 0) exit + ilx = istack(1,istp) + iux = istack(2,istp) + istp = istp - 1 + ! + ! Choose a pivot with median-of-three heuristics, leave it + ! in the LPIV location + ! + i = ilx + j = iux + lpiv = (i+j)/2 + piv = abs(x(lpiv)) + if (piv < abs(x(i))) then + xt = x(i) + x(i) = x(lpiv) + x(lpiv) = xt + piv = abs(x(lpiv)) + endif + if (piv > abs(x(j))) then + xt = x(j) + x(j) = x(lpiv) + x(lpiv) = xt + piv = abs(x(lpiv)) + endif + if (piv < abs(x(i))) then + xt = x(i) + x(i) = x(lpiv) + x(lpiv) = xt + piv = abs(x(lpiv)) + endif + ! + ! now piv is correct; place it into first location + + xt = x(i) + x(i) = x(lpiv) + x(lpiv) = xt + + i = ilx - 1 + j = iux + 1 + + outer_up: do + in_up1: do + i = i + 1 + xk = abs(x(i)) + if (xk >= piv) exit in_up1 + end do in_up1 + ! + ! Ensure finite termination for next loop + ! + xt = x(i) + x(i) = piv + in_up2:do + j = j - 1 + xk = abs(x(j)) + if (xk <= piv) exit in_up2 + end do in_up2 + x(i) = xt + + if (j > i) then + xt = x(i) + x(i) = x(j) + x(j) = xt + else + exit outer_up + end if + end do outer_up + if (i == ilx) then + if (x(i) /= piv) then + call psb_errpush(psb_err_internal_error_, & + & r_name='psi_lqasr',a_err='impossible pivot condition') + call psb_error() + endif + i = i + 1 + endif + + n1 = (i-1)-ilx+1 + n2 = iux-(i)+1 + if (n1 > n2) then + if (n1 > ithrs) then + istp = istp + 1 + istack(1,istp) = ilx + istack(2,istp) = i-1 + else + call psi_laisr_up(n1,x(ilx:i-1)) + endif + if (n2 > ithrs) then + istp = istp + 1 + istack(1,istp) = i + istack(2,istp) = iux + else + call psi_laisr_up(n2,x(i:iux)) + endif + else + if (n2 > ithrs) then + istp = istp + 1 + istack(1,istp) = i + istack(2,istp) = iux + else + call psi_laisr_up(n2,x(i:iux)) + endif + if (n1 > ithrs) then + istp = istp + 1 + istack(1,istp) = ilx + istack(2,istp) = i-1 + else + call psi_laisr_up(n1,x(ilx:i-1)) + endif + endif + enddo + else + call psi_laisr_up(n,x) + endif + +end subroutine psi_laqsr_up + +subroutine psi_laqsr_dw(n,x) + use psb_l_sort_mod, psb_protect_name => psi_laqsr_dw + use psb_error_mod + implicit none + + integer(psb_lpk_), intent(inout) :: x(:) + integer(psb_lpk_), intent(in) :: n + ! .. Local Scalars .. + integer(psb_lpk_) :: piv, xk + integer(psb_lpk_) :: xt + integer(psb_ipk_) :: i, j, ilx, iux, istp, lpiv + integer(psb_lpk_) :: n1, n2 + integer(psb_lpk_) :: ixt + + integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=72 + integer(psb_ipk_) :: istack(nparms,maxstack) + + if (n > ithrs) then + ! + ! Init stack pointer + ! + istp = 1 + istack(1,istp) = 1 + istack(2,istp) = n + + do + if (istp <= 0) exit + ilx = istack(1,istp) + iux = istack(2,istp) + istp = istp - 1 + ! + ! Choose a pivot with median-of-three heuristics, leave it + ! in the LPIV location + ! + i = ilx + j = iux + lpiv = (i+j)/2 + piv = abs(x(lpiv)) + if (piv > abs(x(i))) then + xt = x(i) + x(i) = x(lpiv) + x(lpiv) = xt + piv = abs(x(lpiv)) + endif + if (piv < abs(x(j))) then + xt = x(j) + x(j) = x(lpiv) + x(lpiv) = xt + piv = abs(x(lpiv)) + endif + if (piv > abs(x(i))) then + xt = x(i) + x(i) = x(lpiv) + x(lpiv) = xt + piv = abs(x(lpiv)) + endif + ! + ! now piv is correct; place it into first location + + xt = x(i) + x(i) = x(lpiv) + x(lpiv) = xt + + i = ilx - 1 + j = iux + 1 + + outer_dw: do + in_dw1: do + i = i + 1 + xk = abs(x(i)) + if (xk <= piv) exit in_dw1 + end do in_dw1 + ! + ! Ensure finite termination for next loop + ! + xt = x(i) + x(i) = piv + in_dw2:do + j = j - 1 + xk = abs(x(j)) + if (xk >= piv) exit in_dw2 + end do in_dw2 + x(i) = xt + + if (j > i) then + xt = x(i) + x(i) = x(j) + x(j) = xt + else + exit outer_dw + end if + end do outer_dw + if (i == ilx) then + if (x(i) /= piv) then + call psb_errpush(psb_err_internal_error_,& + & r_name='psi_lqasr',a_err='impossible pivot condition') + call psb_error() + endif + i = i + 1 + endif + + n1 = (i-1)-ilx+1 + n2 = iux-(i)+1 + if (n1 > n2) then + if (n1 > ithrs) then + istp = istp + 1 + istack(1,istp) = ilx + istack(2,istp) = i-1 + else + call psi_laisr_dw(n1,x(ilx:i-1)) + endif + if (n2 > ithrs) then + istp = istp + 1 + istack(1,istp) = i + istack(2,istp) = iux + else + call psi_laisr_dw(n2,x(i:iux)) + endif + else + if (n2 > ithrs) then + istp = istp + 1 + istack(1,istp) = i + istack(2,istp) = iux + else + call psi_laisr_dw(n2,x(i:iux)) + endif + if (n1 > ithrs) then + istp = istp + 1 + istack(1,istp) = ilx + istack(2,istp) = i-1 + else + call psi_laisr_dw(n1,x(ilx:i-1)) + endif + endif + enddo + else + call psi_laisr_dw(n,x) + endif + +end subroutine psi_laqsr_dw + + diff --git a/base/serial/sort/psb_s_isort_impl.f90 b/base/serial/sort/psb_s_isort_impl.f90 index fbac121b8..4d45539c0 100644 --- a/base/serial/sort/psb_s_isort_impl.f90 +++ b/base/serial/sort/psb_s_isort_impl.f90 @@ -48,7 +48,8 @@ subroutine psb_sisort(x,ix,dir,flag) integer(psb_ipk_), optional, intent(in) :: dir, flag integer(psb_ipk_), optional, intent(inout) :: ix(:) - integer(psb_ipk_) :: dir_, flag_, n, err_act, i + integer(psb_ipk_) :: dir_, flag_, err_act + integer(psb_ipk_) :: n, i integer(psb_ipk_) :: ierr(5) character(len=20) :: name diff --git a/base/serial/sort/psb_s_msort_impl.f90 b/base/serial/sort/psb_s_msort_impl.f90 index 53b557123..3b3732912 100644 --- a/base/serial/sort/psb_s_msort_impl.f90 +++ b/base/serial/sort/psb_s_msort_impl.f90 @@ -41,6 +41,99 @@ ! Addison-Wesley ! + subroutine psb_smsort_u(x,nout,dir) + use psb_s_sort_mod, psb_protect_name => psb_smsort_u + use psb_error_mod + implicit none + real(psb_spk_), intent(inout) :: x(:) + integer(psb_ipk_), intent(out) :: nout + integer(psb_ipk_), optional, intent(in) :: dir + + integer(psb_ipk_) :: n, k + integer(psb_ipk_) :: err_act + + integer(psb_ipk_) :: ierr(5) + character(len=20) :: name + + name='psb_msort_u' + call psb_erractionsave(err_act) + + n = size(x) + + call psb_msort(x,dir=dir) + nout = min(1,n) + do k=2,n + if (x(k) /= x(nout)) then + nout = nout + 1 + x(nout) = x(k) + endif + enddo + + return + +9999 call psb_error_handler(err_act) + + return + end subroutine psb_smsort_u + + + function psb_sbsrch(key,n,v) result(ipos) + use psb_s_sort_mod, psb_protect_name => psb_sbsrch + implicit none + integer(psb_ipk_) :: ipos, n + real(psb_spk_) :: key + real(psb_spk_) :: v(:) + + integer(psb_ipk_) :: lb, ub, m, i + + ipos = -1 + if (n<5) then + do i=1,n + if (key.eq.v(i)) then + ipos = i + return + end if + enddo + return + end if + + lb = 1 + ub = n + + do while (lb.le.ub) + m = (lb+ub)/2 + if (key.eq.v(m)) then + ipos = m + lb = ub + 1 + else if (key < v(m)) then + ub = m-1 + else + lb = m + 1 + end if + enddo + return + end function psb_sbsrch + + function psb_sssrch(key,n,v) result(ipos) + use psb_s_sort_mod, psb_protect_name => psb_sssrch + implicit none + integer(psb_ipk_) :: ipos, n + real(psb_spk_) :: key + real(psb_spk_) :: v(:) + + integer(psb_ipk_) :: i + + ipos = -1 + do i=1,n + if (key.eq.v(i)) then + ipos = i + return + end if + enddo + + return + end function psb_sssrch + subroutine psb_smsort(x,ix,dir,flag) use psb_s_sort_mod, psb_protect_name => psb_smsort use psb_error_mod diff --git a/base/serial/sort/psb_s_qsort_impl.f90 b/base/serial/sort/psb_s_qsort_impl.f90 index 7110fa4fc..3bc0d06f2 100644 --- a/base/serial/sort/psb_s_qsort_impl.f90 +++ b/base/serial/sort/psb_s_qsort_impl.f90 @@ -48,8 +48,8 @@ subroutine psb_sqsort(x,ix,dir,flag) integer(psb_ipk_), optional, intent(in) :: dir, flag integer(psb_ipk_), optional, intent(inout) :: ix(:) - integer(psb_ipk_) :: dir_, flag_, n, err_act, i - + integer(psb_ipk_) :: dir_, flag_, err_act, i + integer(psb_ipk_) :: n integer(psb_ipk_) :: ierr(5) character(len=20) :: name @@ -140,7 +140,8 @@ subroutine psi_sqsrx_up(n,x,idx) ! .. 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_) :: n1, n2 + integer(psb_ipk_) :: ixt integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=72 integer(psb_ipk_) :: istack(nparms,maxstack) @@ -293,7 +294,8 @@ subroutine psi_sqsrx_dw(n,x,idx) ! .. 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_) :: n1, n2 + integer(psb_ipk_) :: ixt integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=72 integer(psb_ipk_) :: istack(nparms,maxstack) @@ -731,7 +733,8 @@ subroutine psi_saqsrx_up(n,x,idx) real(psb_spk_) :: piv, xk real(psb_spk_) :: xt integer(psb_ipk_) :: i, j, ilx, iux, istp, lpiv - integer(psb_ipk_) :: ixt, n1, n2 + integer(psb_ipk_) :: n1, n2 + integer(psb_ipk_) :: ixt integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=72 integer(psb_ipk_) :: istack(nparms,maxstack) @@ -887,7 +890,8 @@ subroutine psi_saqsrx_dw(n,x,idx) real(psb_spk_) :: piv, xk real(psb_spk_) :: xt integer(psb_ipk_) :: i, j, ilx, iux, istp, lpiv - integer(psb_ipk_) :: ixt, n1, n2 + integer(psb_ipk_) :: n1, n2 + integer(psb_ipk_) :: ixt integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=72 integer(psb_ipk_) :: istack(nparms,maxstack) @@ -1040,7 +1044,8 @@ subroutine psi_saqsr_up(n,x) real(psb_spk_) :: piv, xk real(psb_spk_) :: xt integer(psb_ipk_) :: i, j, ilx, iux, istp, lpiv - integer(psb_ipk_) :: ixt, n1, n2 + integer(psb_ipk_) :: n1, n2 + integer(psb_ipk_) :: ixt integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=72 integer(psb_ipk_) :: istack(nparms,maxstack) @@ -1180,7 +1185,8 @@ subroutine psi_saqsr_dw(n,x) real(psb_spk_) :: piv, xk real(psb_spk_) :: xt integer(psb_ipk_) :: i, j, ilx, iux, istp, lpiv - integer(psb_ipk_) :: ixt, n1, n2 + integer(psb_ipk_) :: n1, n2 + integer(psb_ipk_) :: ixt integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=72 integer(psb_ipk_) :: istack(nparms,maxstack) diff --git a/base/serial/sort/psb_z_isort_impl.f90 b/base/serial/sort/psb_z_isort_impl.f90 index 98cc055ee..e51dd941a 100644 --- a/base/serial/sort/psb_z_isort_impl.f90 +++ b/base/serial/sort/psb_z_isort_impl.f90 @@ -48,7 +48,8 @@ subroutine psb_zisort(x,ix,dir,flag) integer(psb_ipk_), optional, intent(in) :: dir, flag integer(psb_ipk_), optional, intent(inout) :: ix(:) - integer(psb_ipk_) :: dir_, flag_, n, err_act, i + integer(psb_ipk_) :: dir_, flag_, err_act + integer(psb_ipk_) :: n, i integer(psb_ipk_) :: ierr(5) character(len=20) :: name diff --git a/base/serial/sort/psb_z_msort_impl.f90 b/base/serial/sort/psb_z_msort_impl.f90 index d0f3823bb..201556f96 100644 --- a/base/serial/sort/psb_z_msort_impl.f90 +++ b/base/serial/sort/psb_z_msort_impl.f90 @@ -41,6 +41,41 @@ ! Addison-Wesley ! + subroutine psb_zmsort_u(x,nout,dir) + use psb_z_sort_mod, psb_protect_name => psb_zmsort_u + use psb_error_mod + implicit none + complex(psb_dpk_), intent(inout) :: x(:) + integer(psb_ipk_), intent(out) :: nout + integer(psb_ipk_), optional, intent(in) :: dir + + integer(psb_ipk_) :: n, k + integer(psb_ipk_) :: err_act + + integer(psb_ipk_) :: ierr(5) + character(len=20) :: name + + name='psb_msort_u' + call psb_erractionsave(err_act) + + n = size(x) + + call psb_msort(x,dir=dir) + nout = min(1,n) + do k=2,n + if (x(k) /= x(nout)) then + nout = nout + 1 + x(nout) = x(k) + endif + enddo + + return + +9999 call psb_error_handler(err_act) + + return + end subroutine psb_zmsort_u + diff --git a/base/serial/sort/psb_z_qsort_impl.f90 b/base/serial/sort/psb_z_qsort_impl.f90 index cdb58dccd..5acef4fd4 100644 --- a/base/serial/sort/psb_z_qsort_impl.f90 +++ b/base/serial/sort/psb_z_qsort_impl.f90 @@ -48,8 +48,8 @@ subroutine psb_zqsort(x,ix,dir,flag) integer(psb_ipk_), optional, intent(in) :: dir, flag integer(psb_ipk_), optional, intent(inout) :: ix(:) - integer(psb_ipk_) :: dir_, flag_, n, err_act, i - + integer(psb_ipk_) :: dir_, flag_, err_act, i + integer(psb_ipk_) :: n integer(psb_ipk_) :: ierr(5) character(len=20) :: name @@ -150,7 +150,8 @@ subroutine psi_zlqsrx_up(n,x,idx) ! .. 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_) :: n1, n2 + integer(psb_ipk_) :: ixt integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=24 integer(psb_ipk_) :: istack(nparms,maxstack) @@ -306,7 +307,8 @@ subroutine psi_zlqsrx_dw(n,x,idx) ! .. 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_) :: n1, n2 + integer(psb_ipk_) :: ixt integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=24 integer(psb_ipk_) :: istack(nparms,maxstack) @@ -744,7 +746,8 @@ subroutine psi_zalqsrx_up(n,x,idx) ! .. 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_) :: n1, n2 + integer(psb_ipk_) :: ixt integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=24 integer(psb_ipk_) :: istack(nparms,maxstack) @@ -899,7 +902,8 @@ subroutine psi_zalqsrx_dw(n,x,idx) ! .. 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_) :: n1, n2 + integer(psb_ipk_) :: ixt integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=24 integer(psb_ipk_) :: istack(nparms,maxstack) @@ -1335,7 +1339,8 @@ subroutine psi_zaqsrx_up(n,x,idx) real(psb_dpk_) :: piv, xk complex(psb_dpk_) :: xt integer(psb_ipk_) :: i, j, ilx, iux, istp, lpiv - integer(psb_ipk_) :: ixt, n1, n2 + integer(psb_ipk_) :: n1, n2 + integer(psb_ipk_) :: ixt integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=24 integer(psb_ipk_) :: istack(nparms,maxstack) @@ -1491,7 +1496,8 @@ subroutine psi_zaqsrx_dw(n,x,idx) real(psb_dpk_) :: piv, xk complex(psb_dpk_) :: xt integer(psb_ipk_) :: i, j, ilx, iux, istp, lpiv - integer(psb_ipk_) :: ixt, n1, n2 + integer(psb_ipk_) :: n1, n2 + integer(psb_ipk_) :: ixt integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=24 integer(psb_ipk_) :: istack(nparms,maxstack) @@ -1644,7 +1650,8 @@ subroutine psi_zaqsr_up(n,x) real(psb_dpk_) :: piv, xk complex(psb_dpk_) :: xt integer(psb_ipk_) :: i, j, ilx, iux, istp, lpiv - integer(psb_ipk_) :: ixt, n1, n2 + integer(psb_ipk_) :: n1, n2 + integer(psb_ipk_) :: ixt integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=24 integer(psb_ipk_) :: istack(nparms,maxstack) @@ -1784,7 +1791,8 @@ subroutine psi_zaqsr_dw(n,x) real(psb_dpk_) :: piv, xk complex(psb_dpk_) :: xt integer(psb_ipk_) :: i, j, ilx, iux, istp, lpiv - integer(psb_ipk_) :: ixt, n1, n2 + integer(psb_ipk_) :: n1, n2 + integer(psb_ipk_) :: ixt integer(psb_ipk_), parameter :: maxstack=64,nparms=3,ithrs=24 integer(psb_ipk_) :: istack(nparms,maxstack)