Changes for L version of sorting routines.
Compiles, but needs thorough testing.ILmat
Salvatore Filippone
8 years ago
parent
bb96330459
commit
51b1d4ade8
@ -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
|
||||
@ -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
|
||||
@ -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
|
||||
@ -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
|
||||
@ -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)<heap%last))) then
|
||||
write(iout,*) 'Inconsistent size/allocation status!!'
|
||||
else
|
||||
write(iout,*) heap%keys(1:heap%last)
|
||||
end if
|
||||
end subroutine psb_l_dump_heap
|
||||
|
||||
subroutine psb_l_free_heap(heap,info)
|
||||
implicit none
|
||||
class(psb_l_heap), intent(inout) :: heap
|
||||
integer(psb_ipk_), intent(out) :: info
|
||||
|
||||
info=psb_success_
|
||||
if (allocated(heap%keys)) deallocate(heap%keys,stat=info)
|
||||
|
||||
end subroutine psb_l_free_heap
|
||||
|
||||
subroutine psb_l_idx_init_heap(heap,info,dir)
|
||||
use psb_realloc_mod, only : psb_ensure_size
|
||||
implicit none
|
||||
class(psb_l_idx_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)
|
||||
call psb_ensure_size(psb_heap_resize,heap%idxs,info)
|
||||
return
|
||||
end subroutine psb_l_idx_init_heap
|
||||
|
||||
|
||||
function psb_l_idx_howmany(heap) result(res)
|
||||
implicit none
|
||||
class(psb_l_idx_heap), intent(in) :: heap
|
||||
integer(psb_ipk_) :: res
|
||||
res = heap%last
|
||||
end function psb_l_idx_howmany
|
||||
|
||||
subroutine psb_l_idx_insert_heap(key,index,heap,info)
|
||||
use psb_realloc_mod, only : psb_ensure_size
|
||||
implicit none
|
||||
|
||||
integer(psb_lpk_), intent(in) :: key
|
||||
integer(psb_lpk_), intent(in) :: index
|
||||
class(psb_l_idx_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_) &
|
||||
& call psb_ensure_size(heap%last+1,heap%idxs,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_idx_insert_heap(key,index,&
|
||||
& heap%last,heap%keys,heap%idxs,heap%dir,info)
|
||||
|
||||
return
|
||||
end subroutine psb_l_idx_insert_heap
|
||||
|
||||
subroutine psb_l_idx_heap_get_first(key,index,heap,info)
|
||||
implicit none
|
||||
|
||||
class(psb_l_idx_heap), intent(inout) :: heap
|
||||
integer(psb_lpk_), intent(out) :: index
|
||||
integer(psb_ipk_), intent(out) :: info
|
||||
integer(psb_lpk_), intent(out) :: key
|
||||
|
||||
|
||||
info = psb_success_
|
||||
|
||||
call psi_l_idx_heap_get_first(key,index,&
|
||||
& heap%last,heap%keys,heap%idxs,heap%dir,info)
|
||||
|
||||
return
|
||||
end subroutine psb_l_idx_heap_get_first
|
||||
|
||||
subroutine psb_l_idx_dump_heap(iout,heap,info)
|
||||
|
||||
implicit none
|
||||
class(psb_l_idx_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)<heap%last))) then
|
||||
write(iout,*) 'Inconsistent size/allocation status!!'
|
||||
else if ((heap%last > 0).and.((.not.allocated(heap%idxs)).or.&
|
||||
& (size(heap%idxs)<heap%last))) then
|
||||
write(iout,*) 'Inconsistent size/allocation status!!'
|
||||
else
|
||||
write(iout,*) heap%keys(1:heap%last)
|
||||
write(iout,*) heap%idxs(1:heap%last)
|
||||
end if
|
||||
end subroutine psb_l_idx_dump_heap
|
||||
|
||||
subroutine psb_l_idx_free_heap(heap,info)
|
||||
implicit none
|
||||
class(psb_l_idx_heap), intent(inout) :: heap
|
||||
integer(psb_ipk_), intent(out) :: info
|
||||
|
||||
info=psb_success_
|
||||
if (allocated(heap%keys)) deallocate(heap%keys,stat=info)
|
||||
if ((info == psb_success_).and.(allocated(heap%idxs))) deallocate(heap%idxs,stat=info)
|
||||
|
||||
end subroutine psb_l_idx_free_heap
|
||||
|
||||
end module psb_l_sort_mod
|
||||
@ -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_s_ip_reord_mod
|
||||
use psb_const_mod
|
||||
|
||||
interface psb_ip_reord
|
||||
module procedure psb_ip_reord_s1,&
|
||||
& psb_ip_reord_s1i1, psb_ip_reord_s1i2,&
|
||||
& psb_ip_reord_s1i3
|
||||
|
||||
|
||||
end interface
|
||||
|
||||
contains
|
||||
|
||||
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_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
|
||||
@ -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
|
||||
@ -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
|
||||
|
||||
|
||||
|
||||
|
||||
@ -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
|
||||
|
||||
@ -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
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
File diff suppressed because it is too large
Load Diff
Loading…
Reference in New Issue