Merge branch 'dev-openmp' of github.com:sfilippone/psblas3 into dev-openmp
Salvatore Filippone
2 years ago
commit
fd0b1482e5
@ -1,578 +0,0 @@
|
||||
!
|
||||
! 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_i_sort_mod
|
||||
use psb_const_mod
|
||||
|
||||
interface psb_isaperm
|
||||
logical function psb_iisaperm(n,eip)
|
||||
import
|
||||
integer(psb_ipk_), intent(in) :: n
|
||||
integer(psb_ipk_), intent(in) :: eip(n)
|
||||
end function psb_iisaperm
|
||||
end interface psb_isaperm
|
||||
|
||||
|
||||
interface psb_msort_unique
|
||||
subroutine psb_imsort_u(x,nout,dir)
|
||||
import
|
||||
integer(psb_ipk_), intent(inout) :: x(:)
|
||||
integer(psb_ipk_), intent(out) :: nout
|
||||
integer(psb_ipk_), optional, intent(in) :: dir
|
||||
end subroutine psb_imsort_u
|
||||
end interface psb_msort_unique
|
||||
|
||||
type psb_i_heap
|
||||
integer(psb_ipk_) :: last, dir
|
||||
integer(psb_ipk_), allocatable :: keys(:)
|
||||
contains
|
||||
procedure, pass(heap) :: init => psb_i_init_heap
|
||||
procedure, pass(heap) :: howmany => psb_i_howmany
|
||||
procedure, pass(heap) :: insert => psb_i_insert_heap
|
||||
procedure, pass(heap) :: get_first => psb_i_heap_get_first
|
||||
procedure, pass(heap) :: dump => psb_i_dump_heap
|
||||
procedure, pass(heap) :: free => psb_i_free_heap
|
||||
end type psb_i_heap
|
||||
|
||||
type psb_i_idx_heap
|
||||
integer(psb_ipk_) :: last, dir
|
||||
integer(psb_ipk_), allocatable :: keys(:)
|
||||
integer(psb_ipk_), allocatable :: idxs(:)
|
||||
contains
|
||||
procedure, pass(heap) :: init => psb_i_idx_init_heap
|
||||
procedure, pass(heap) :: howmany => psb_i_idx_howmany
|
||||
procedure, pass(heap) :: insert => psb_i_idx_insert_heap
|
||||
procedure, pass(heap) :: get_first => psb_i_idx_heap_get_first
|
||||
procedure, pass(heap) :: dump => psb_i_idx_dump_heap
|
||||
procedure, pass(heap) :: free => psb_i_idx_free_heap
|
||||
end type psb_i_idx_heap
|
||||
|
||||
|
||||
interface psb_msort
|
||||
subroutine psb_imsort(x,ix,dir,flag)
|
||||
import
|
||||
integer(psb_ipk_), intent(inout) :: x(:)
|
||||
integer(psb_ipk_), optional, intent(in) :: dir, flag
|
||||
integer(psb_ipk_), optional, intent(inout) :: ix(:)
|
||||
end subroutine psb_imsort
|
||||
end interface psb_msort
|
||||
|
||||
|
||||
interface psb_bsrch
|
||||
function psb_ibsrch(key,n,v) result(ipos)
|
||||
import
|
||||
integer(psb_ipk_) :: ipos, n
|
||||
integer(psb_ipk_) :: key
|
||||
integer(psb_ipk_) :: v(:)
|
||||
end function psb_ibsrch
|
||||
end interface psb_bsrch
|
||||
|
||||
interface psb_ssrch
|
||||
function psb_issrch(key,n,v) result(ipos)
|
||||
import
|
||||
implicit none
|
||||
integer(psb_ipk_) :: ipos, n
|
||||
integer(psb_ipk_) :: key
|
||||
integer(psb_ipk_) :: v(:)
|
||||
end function psb_issrch
|
||||
end interface psb_ssrch
|
||||
|
||||
interface
|
||||
subroutine psi_i_msort_up(n,k,l,iret)
|
||||
import
|
||||
implicit none
|
||||
integer(psb_ipk_) :: n, iret
|
||||
integer(psb_ipk_) :: k(n)
|
||||
integer(psb_ipk_) :: l(0:n+1)
|
||||
end subroutine psi_i_msort_up
|
||||
subroutine psi_i_msort_dw(n,k,l,iret)
|
||||
import
|
||||
implicit none
|
||||
integer(psb_ipk_) :: n, iret
|
||||
integer(psb_ipk_) :: k(n)
|
||||
integer(psb_ipk_) :: l(0:n+1)
|
||||
end subroutine psi_i_msort_dw
|
||||
end interface
|
||||
interface
|
||||
subroutine psi_i_amsort_up(n,k,l,iret)
|
||||
import
|
||||
implicit none
|
||||
integer(psb_ipk_) :: n, iret
|
||||
integer(psb_ipk_) :: k(n)
|
||||
integer(psb_ipk_) :: l(0:n+1)
|
||||
end subroutine psi_i_amsort_up
|
||||
subroutine psi_i_amsort_dw(n,k,l,iret)
|
||||
import
|
||||
implicit none
|
||||
integer(psb_ipk_) :: n, iret
|
||||
integer(psb_ipk_) :: k(n)
|
||||
integer(psb_ipk_) :: l(0:n+1)
|
||||
end subroutine psi_i_amsort_dw
|
||||
end interface
|
||||
|
||||
|
||||
interface psb_qsort
|
||||
subroutine psb_iqsort(x,ix,dir,flag)
|
||||
import
|
||||
integer(psb_ipk_), intent(inout) :: x(:)
|
||||
integer(psb_ipk_), optional, intent(in) :: dir, flag
|
||||
integer(psb_ipk_), optional, intent(inout) :: ix(:)
|
||||
end subroutine psb_iqsort
|
||||
end interface psb_qsort
|
||||
|
||||
interface psb_isort
|
||||
subroutine psb_iisort(x,ix,dir,flag)
|
||||
import
|
||||
integer(psb_ipk_), intent(inout) :: x(:)
|
||||
integer(psb_ipk_), optional, intent(in) :: dir, flag
|
||||
integer(psb_ipk_), optional, intent(inout) :: ix(:)
|
||||
end subroutine psb_iisort
|
||||
end interface psb_isort
|
||||
|
||||
|
||||
interface psb_hsort
|
||||
subroutine psb_ihsort(x,ix,dir,flag)
|
||||
import
|
||||
integer(psb_ipk_), intent(inout) :: x(:)
|
||||
integer(psb_ipk_), optional, intent(in) :: dir, flag
|
||||
integer(psb_ipk_), optional, intent(inout) :: ix(:)
|
||||
end subroutine psb_ihsort
|
||||
end interface psb_hsort
|
||||
|
||||
|
||||
interface
|
||||
subroutine psi_i_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_ipk_), intent(in) :: key
|
||||
integer(psb_ipk_), intent(inout) :: heap(:)
|
||||
integer(psb_ipk_), intent(in) :: dir
|
||||
integer(psb_ipk_), intent(inout) :: last
|
||||
integer(psb_ipk_), intent(out) :: info
|
||||
end subroutine psi_i_insert_heap
|
||||
end interface
|
||||
|
||||
interface
|
||||
subroutine psi_i_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_ipk_), intent(in) :: key
|
||||
integer(psb_ipk_), intent(inout) :: heap(:)
|
||||
integer(psb_ipk_), intent(in) :: index
|
||||
integer(psb_ipk_), intent(in) :: dir
|
||||
integer(psb_ipk_), intent(inout) :: idxs(:)
|
||||
integer(psb_ipk_), intent(inout) :: last
|
||||
integer(psb_ipk_), intent(out) :: info
|
||||
end subroutine psi_i_idx_insert_heap
|
||||
end interface
|
||||
|
||||
|
||||
interface
|
||||
subroutine psi_i_heap_get_first(key,last,heap,dir,info)
|
||||
import
|
||||
implicit none
|
||||
integer(psb_ipk_), intent(inout) :: key
|
||||
integer(psb_ipk_), intent(inout) :: last
|
||||
integer(psb_ipk_), intent(in) :: dir
|
||||
integer(psb_ipk_), intent(inout) :: heap(:)
|
||||
integer(psb_ipk_), intent(out) :: info
|
||||
end subroutine psi_i_heap_get_first
|
||||
end interface
|
||||
|
||||
interface
|
||||
subroutine psi_i_idx_heap_get_first(key,index,last,heap,idxs,dir,info)
|
||||
import
|
||||
integer(psb_ipk_), intent(inout) :: key
|
||||
integer(psb_ipk_), intent(out) :: index
|
||||
integer(psb_ipk_), intent(inout) :: heap(:)
|
||||
integer(psb_ipk_), intent(in) :: dir
|
||||
integer(psb_ipk_), intent(inout) :: last
|
||||
integer(psb_ipk_), intent(inout) :: idxs(:)
|
||||
integer(psb_ipk_), intent(out) :: info
|
||||
end subroutine psi_i_idx_heap_get_first
|
||||
end interface
|
||||
|
||||
interface
|
||||
subroutine psi_iisrx_up(n,x,ix)
|
||||
import
|
||||
integer(psb_ipk_), intent(inout) :: x(:)
|
||||
integer(psb_ipk_), intent(inout) :: ix(:)
|
||||
integer(psb_ipk_), intent(in) :: n
|
||||
end subroutine psi_iisrx_up
|
||||
subroutine psi_iisrx_dw(n,x,ix)
|
||||
import
|
||||
integer(psb_ipk_), intent(inout) :: x(:)
|
||||
integer(psb_ipk_), intent(inout) :: ix(:)
|
||||
integer(psb_ipk_), intent(in) :: n
|
||||
end subroutine psi_iisrx_dw
|
||||
subroutine psi_iisr_up(n,x)
|
||||
import
|
||||
integer(psb_ipk_), intent(inout) :: x(:)
|
||||
integer(psb_ipk_), intent(in) :: n
|
||||
end subroutine psi_iisr_up
|
||||
subroutine psi_iisr_dw(n,x)
|
||||
import
|
||||
integer(psb_ipk_), intent(inout) :: x(:)
|
||||
integer(psb_ipk_), intent(in) :: n
|
||||
end subroutine psi_iisr_dw
|
||||
subroutine psi_iaisrx_up(n,x,ix)
|
||||
import
|
||||
integer(psb_ipk_), intent(inout) :: x(:)
|
||||
integer(psb_ipk_), intent(inout) :: ix(:)
|
||||
integer(psb_ipk_), intent(in) :: n
|
||||
end subroutine psi_iaisrx_up
|
||||
subroutine psi_iaisrx_dw(n,x,ix)
|
||||
import
|
||||
integer(psb_ipk_), intent(inout) :: x(:)
|
||||
integer(psb_ipk_), intent(inout) :: ix(:)
|
||||
integer(psb_ipk_), intent(in) :: n
|
||||
end subroutine psi_iaisrx_dw
|
||||
subroutine psi_iaisr_up(n,x)
|
||||
import
|
||||
integer(psb_ipk_), intent(inout) :: x(:)
|
||||
integer(psb_ipk_), intent(in) :: n
|
||||
end subroutine psi_iaisr_up
|
||||
subroutine psi_iaisr_dw(n,x)
|
||||
import
|
||||
integer(psb_ipk_), intent(inout) :: x(:)
|
||||
integer(psb_ipk_), intent(in) :: n
|
||||
end subroutine psi_iaisr_dw
|
||||
end interface
|
||||
|
||||
interface
|
||||
subroutine psi_iqsrx_up(n,x,ix)
|
||||
import
|
||||
integer(psb_ipk_), intent(inout) :: x(:)
|
||||
integer(psb_ipk_), intent(inout) :: ix(:)
|
||||
integer(psb_ipk_), intent(in) :: n
|
||||
end subroutine psi_iqsrx_up
|
||||
subroutine psi_iqsrx_dw(n,x,ix)
|
||||
import
|
||||
integer(psb_ipk_), intent(inout) :: x(:)
|
||||
integer(psb_ipk_), intent(inout) :: ix(:)
|
||||
integer(psb_ipk_), intent(in) :: n
|
||||
end subroutine psi_iqsrx_dw
|
||||
subroutine psi_iqsr_up(n,x)
|
||||
import
|
||||
integer(psb_ipk_), intent(inout) :: x(:)
|
||||
integer(psb_ipk_), intent(in) :: n
|
||||
end subroutine psi_iqsr_up
|
||||
subroutine psi_iqsr_dw(n,x)
|
||||
import
|
||||
integer(psb_ipk_), intent(inout) :: x(:)
|
||||
integer(psb_ipk_), intent(in) :: n
|
||||
end subroutine psi_iqsr_dw
|
||||
subroutine psi_iaqsrx_up(n,x,ix)
|
||||
import
|
||||
integer(psb_ipk_), intent(inout) :: x(:)
|
||||
integer(psb_ipk_), intent(inout) :: ix(:)
|
||||
integer(psb_ipk_), intent(in) :: n
|
||||
end subroutine psi_iaqsrx_up
|
||||
subroutine psi_iaqsrx_dw(n,x,ix)
|
||||
import
|
||||
integer(psb_ipk_), intent(inout) :: x(:)
|
||||
integer(psb_ipk_), intent(inout) :: ix(:)
|
||||
integer(psb_ipk_), intent(in) :: n
|
||||
end subroutine psi_iaqsrx_dw
|
||||
subroutine psi_iaqsr_up(n,x)
|
||||
import
|
||||
integer(psb_ipk_), intent(inout) :: x(:)
|
||||
integer(psb_ipk_), intent(in) :: n
|
||||
end subroutine psi_iaqsr_up
|
||||
subroutine psi_iaqsr_dw(n,x)
|
||||
import
|
||||
integer(psb_ipk_), intent(inout) :: x(:)
|
||||
integer(psb_ipk_), intent(in) :: n
|
||||
end subroutine psi_iaqsr_dw
|
||||
end interface
|
||||
|
||||
contains
|
||||
|
||||
subroutine psb_i_init_heap(heap,info,dir)
|
||||
use psb_realloc_mod, only : psb_ensure_size
|
||||
implicit none
|
||||
class(psb_i_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_i_init_heap
|
||||
|
||||
|
||||
function psb_i_howmany(heap) result(res)
|
||||
implicit none
|
||||
class(psb_i_heap), intent(in) :: heap
|
||||
integer(psb_ipk_) :: res
|
||||
res = heap%last
|
||||
end function psb_i_howmany
|
||||
|
||||
subroutine psb_i_insert_heap(key,heap,info)
|
||||
use psb_realloc_mod, only : psb_ensure_size
|
||||
implicit none
|
||||
|
||||
integer(psb_ipk_), intent(in) :: key
|
||||
class(psb_i_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_i_insert_heap(key,&
|
||||
& heap%last,heap%keys,heap%dir,info)
|
||||
|
||||
return
|
||||
end subroutine psb_i_insert_heap
|
||||
|
||||
subroutine psb_i_heap_get_first(key,heap,info)
|
||||
implicit none
|
||||
|
||||
class(psb_i_heap), intent(inout) :: heap
|
||||
integer(psb_ipk_), intent(out) :: info
|
||||
integer(psb_ipk_), intent(out) :: key
|
||||
|
||||
|
||||
info = psb_success_
|
||||
|
||||
call psi_i_heap_get_first(key,&
|
||||
& heap%last,heap%keys,heap%dir,info)
|
||||
|
||||
return
|
||||
end subroutine psb_i_heap_get_first
|
||||
|
||||
subroutine psb_i_dump_heap(iout,heap,info)
|
||||
|
||||
implicit none
|
||||
class(psb_i_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_i_dump_heap
|
||||
|
||||
subroutine psb_i_free_heap(heap,info)
|
||||
implicit none
|
||||
class(psb_i_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_i_free_heap
|
||||
|
||||
subroutine psb_i_idx_init_heap(heap,info,dir)
|
||||
use psb_realloc_mod, only : psb_ensure_size
|
||||
implicit none
|
||||
class(psb_i_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_i_idx_init_heap
|
||||
|
||||
|
||||
function psb_i_idx_howmany(heap) result(res)
|
||||
implicit none
|
||||
class(psb_i_idx_heap), intent(in) :: heap
|
||||
integer(psb_ipk_) :: res
|
||||
res = heap%last
|
||||
end function psb_i_idx_howmany
|
||||
|
||||
subroutine psb_i_idx_insert_heap(key,index,heap,info)
|
||||
use psb_realloc_mod, only : psb_ensure_size
|
||||
implicit none
|
||||
|
||||
integer(psb_ipk_), intent(in) :: key
|
||||
integer(psb_ipk_), intent(in) :: index
|
||||
class(psb_i_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_i_idx_insert_heap(key,index,&
|
||||
& heap%last,heap%keys,heap%idxs,heap%dir,info)
|
||||
|
||||
return
|
||||
end subroutine psb_i_idx_insert_heap
|
||||
|
||||
subroutine psb_i_idx_heap_get_first(key,index,heap,info)
|
||||
implicit none
|
||||
|
||||
class(psb_i_idx_heap), intent(inout) :: heap
|
||||
integer(psb_ipk_), intent(out) :: index
|
||||
integer(psb_ipk_), intent(out) :: info
|
||||
integer(psb_ipk_), intent(out) :: key
|
||||
|
||||
|
||||
info = psb_success_
|
||||
|
||||
call psi_i_idx_heap_get_first(key,index,&
|
||||
& heap%last,heap%keys,heap%idxs,heap%dir,info)
|
||||
|
||||
return
|
||||
end subroutine psb_i_idx_heap_get_first
|
||||
|
||||
subroutine psb_i_idx_dump_heap(iout,heap,info)
|
||||
|
||||
implicit none
|
||||
class(psb_i_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_i_idx_dump_heap
|
||||
|
||||
subroutine psb_i_idx_free_heap(heap,info)
|
||||
implicit none
|
||||
class(psb_i_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_i_idx_free_heap
|
||||
|
||||
end module psb_i_sort_mod
|
||||
@ -1,578 +0,0 @@
|
||||
!
|
||||
! 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
|
||||
File diff suppressed because it is too large
Load Diff
File diff suppressed because it is too large
Load Diff
File diff suppressed because it is too large
Load Diff
File diff suppressed because it is too large
Load Diff
File diff suppressed because it is too large
Load Diff
File diff suppressed because it is too large
Load Diff
File diff suppressed because it is too large
Load Diff
@ -1,678 +0,0 @@
|
||||
!
|
||||
! 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_ihsort(x,ix,dir,flag)
|
||||
use psb_sort_mod, psb_protect_name => psb_ihsort
|
||||
use psb_error_mod
|
||||
implicit none
|
||||
integer(psb_ipk_), intent(inout) :: x(:)
|
||||
integer(psb_ipk_), optional, intent(in) :: dir, flag
|
||||
integer(psb_ipk_), optional, intent(inout) :: ix(:)
|
||||
|
||||
integer(psb_ipk_) :: dir_, flag_, n, i, l, err_act,info
|
||||
integer(psb_ipk_) :: key
|
||||
integer(psb_ipk_) :: 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_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_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_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_i_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_ihsort
|
||||
|
||||
|
||||
|
||||
!
|
||||
! 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_i_insert_heap(key,last,heap,dir,info)
|
||||
use psb_sort_mod, psb_protect_name => psi_i_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_ipk_), intent(in) :: key
|
||||
integer(psb_ipk_), intent(in) :: dir
|
||||
integer(psb_ipk_), intent(inout) :: heap(:)
|
||||
integer(psb_ipk_), intent(inout) :: last
|
||||
integer(psb_ipk_), intent(out) :: info
|
||||
integer(psb_ipk_) :: i, i2
|
||||
integer(psb_ipk_) :: 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_i_insert_heap
|
||||
|
||||
|
||||
subroutine psi_i_heap_get_first(key,last,heap,dir,info)
|
||||
use psb_sort_mod, psb_protect_name => psi_i_heap_get_first
|
||||
implicit none
|
||||
|
||||
integer(psb_ipk_), intent(inout) :: key
|
||||
integer(psb_ipk_), intent(inout) :: last
|
||||
integer(psb_ipk_), intent(in) :: dir
|
||||
integer(psb_ipk_), intent(inout) :: heap(:)
|
||||
integer(psb_ipk_), intent(out) :: info
|
||||
|
||||
integer(psb_ipk_) :: i, j
|
||||
integer(psb_ipk_) :: 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_i_heap_get_first
|
||||
|
||||
|
||||
subroutine psi_i_idx_insert_heap(key,index,last,heap,idxs,dir,info)
|
||||
use psb_sort_mod, psb_protect_name => psi_i_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_ipk_), intent(in) :: key
|
||||
integer(psb_ipk_), intent(in) :: index,dir
|
||||
integer(psb_ipk_), intent(inout) :: heap(:)
|
||||
integer(psb_ipk_), intent(inout) :: idxs(:),last
|
||||
integer(psb_ipk_), intent(out) :: info
|
||||
integer(psb_ipk_) :: i, i2, itemp
|
||||
integer(psb_ipk_) :: 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_i_idx_insert_heap
|
||||
|
||||
subroutine psi_i_idx_heap_get_first(key,index,last,heap,idxs,dir,info)
|
||||
use psb_sort_mod, psb_protect_name => psi_i_idx_heap_get_first
|
||||
implicit none
|
||||
|
||||
integer(psb_ipk_), 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_ipk_), intent(out) :: key
|
||||
|
||||
integer(psb_ipk_) :: i, j,itemp
|
||||
integer(psb_ipk_) :: 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_i_idx_heap_get_first
|
||||
|
||||
|
||||
|
||||
|
||||
@ -1,341 +0,0 @@
|
||||
!
|
||||
! 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_iisort(x,ix,dir,flag)
|
||||
use psb_sort_mod, psb_protect_name => psb_iisort
|
||||
use psb_error_mod
|
||||
implicit none
|
||||
integer(psb_ipk_), intent(inout) :: x(:)
|
||||
integer(psb_ipk_), optional, intent(in) :: dir, flag
|
||||
integer(psb_ipk_), optional, intent(inout) :: ix(:)
|
||||
|
||||
integer(psb_ipk_) :: dir_, flag_, err_act
|
||||
integer(psb_ipk_) :: n, i
|
||||
|
||||
integer(psb_ipk_) :: ierr(5)
|
||||
character(len=20) :: name
|
||||
|
||||
name='psb_iisort'
|
||||
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_iisrx_up(n,x,ix)
|
||||
case (psb_sort_down_)
|
||||
call psi_iisrx_dw(n,x,ix)
|
||||
case (psb_asort_up_)
|
||||
call psi_iaisrx_up(n,x,ix)
|
||||
case (psb_asort_down_)
|
||||
call psi_iaisrx_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_iisr_up(n,x)
|
||||
case (psb_sort_down_)
|
||||
call psi_iisr_dw(n,x)
|
||||
case (psb_asort_up_)
|
||||
call psi_iaisr_up(n,x)
|
||||
case (psb_asort_down_)
|
||||
call psi_iaisr_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_iisort
|
||||
|
||||
subroutine psi_iisrx_up(n,x,idx)
|
||||
use psb_sort_mod, psb_protect_name => psi_iisrx_up
|
||||
use psb_error_mod
|
||||
implicit none
|
||||
integer(psb_ipk_), intent(inout) :: x(:)
|
||||
integer(psb_ipk_), intent(inout) :: idx(:)
|
||||
integer(psb_ipk_), intent(in) :: n
|
||||
integer(psb_ipk_) :: i,j,ix
|
||||
integer(psb_ipk_) :: 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_iisrx_up
|
||||
|
||||
subroutine psi_iisrx_dw(n,x,idx)
|
||||
use psb_sort_mod, psb_protect_name => psi_iisrx_dw
|
||||
use psb_error_mod
|
||||
implicit none
|
||||
integer(psb_ipk_), intent(inout) :: x(:)
|
||||
integer(psb_ipk_), intent(inout) :: idx(:)
|
||||
integer(psb_ipk_), intent(in) :: n
|
||||
integer(psb_ipk_) :: i,j,ix
|
||||
integer(psb_ipk_) :: 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_iisrx_dw
|
||||
|
||||
|
||||
subroutine psi_iisr_up(n,x)
|
||||
use psb_sort_mod, psb_protect_name => psi_iisr_up
|
||||
use psb_error_mod
|
||||
implicit none
|
||||
integer(psb_ipk_), intent(inout) :: x(:)
|
||||
integer(psb_ipk_), intent(in) :: n
|
||||
integer(psb_ipk_) :: i,j
|
||||
integer(psb_ipk_) :: 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_iisr_up
|
||||
|
||||
subroutine psi_iisr_dw(n,x)
|
||||
use psb_sort_mod, psb_protect_name => psi_iisr_dw
|
||||
use psb_error_mod
|
||||
implicit none
|
||||
integer(psb_ipk_), intent(inout) :: x(:)
|
||||
integer(psb_ipk_), intent(in) :: n
|
||||
integer(psb_ipk_) :: i,j
|
||||
integer(psb_ipk_) :: 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_iisr_dw
|
||||
|
||||
subroutine psi_iaisrx_up(n,x,idx)
|
||||
use psb_sort_mod, psb_protect_name => psi_iaisrx_up
|
||||
use psb_error_mod
|
||||
implicit none
|
||||
integer(psb_ipk_), intent(inout) :: x(:)
|
||||
integer(psb_ipk_), intent(inout) :: idx(:)
|
||||
integer(psb_ipk_), intent(in) :: n
|
||||
integer(psb_ipk_) :: i,j,ix
|
||||
integer(psb_ipk_) :: 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_iaisrx_up
|
||||
|
||||
subroutine psi_iaisrx_dw(n,x,idx)
|
||||
use psb_sort_mod, psb_protect_name => psi_iaisrx_dw
|
||||
use psb_error_mod
|
||||
implicit none
|
||||
integer(psb_ipk_), intent(inout) :: x(:)
|
||||
integer(psb_ipk_), intent(inout) :: idx(:)
|
||||
integer(psb_ipk_), intent(in) :: n
|
||||
integer(psb_ipk_) :: i,j,ix
|
||||
integer(psb_ipk_) :: 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_iaisrx_dw
|
||||
|
||||
subroutine psi_iaisr_up(n,x)
|
||||
use psb_sort_mod, psb_protect_name => psi_iaisr_up
|
||||
use psb_error_mod
|
||||
implicit none
|
||||
integer(psb_ipk_), intent(inout) :: x(:)
|
||||
integer(psb_ipk_), intent(in) :: n
|
||||
integer(psb_ipk_) :: i,j
|
||||
integer(psb_ipk_) :: 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_iaisr_up
|
||||
|
||||
subroutine psi_iaisr_dw(n,x)
|
||||
use psb_sort_mod, psb_protect_name => psi_iaisr_dw
|
||||
use psb_error_mod
|
||||
implicit none
|
||||
integer(psb_ipk_), intent(inout) :: x(:)
|
||||
integer(psb_ipk_), intent(in) :: n
|
||||
integer(psb_ipk_) :: i,j
|
||||
integer(psb_ipk_) :: 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_iaisr_dw
|
||||
|
||||
@ -1,713 +0,0 @@
|
||||
!
|
||||
! 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_iisaperm(n,eip)
|
||||
use psb_sort_mod, psb_protect_name => psb_iisaperm
|
||||
implicit none
|
||||
|
||||
integer(psb_ipk_), intent(in) :: n
|
||||
integer(psb_ipk_), intent(in) :: eip(n)
|
||||
integer(psb_ipk_), allocatable :: ip(:)
|
||||
integer(psb_ipk_) :: i,j,m, info
|
||||
|
||||
|
||||
psb_iisaperm = .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_iisaperm = .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_iisaperm = .false.
|
||||
goto 9999
|
||||
endif
|
||||
end if
|
||||
enddo
|
||||
9999 continue
|
||||
|
||||
return
|
||||
end function psb_iisaperm
|
||||
|
||||
|
||||
subroutine psb_imsort_u(x,nout,dir)
|
||||
use psb_sort_mod, psb_protect_name => psb_imsort_u
|
||||
use psb_error_mod
|
||||
implicit none
|
||||
integer(psb_ipk_), intent(inout) :: x(:)
|
||||
integer(psb_ipk_), intent(out) :: nout
|
||||
integer(psb_ipk_), optional, intent(in) :: dir
|
||||
|
||||
integer(psb_ipk_) :: n, k
|
||||
integer(psb_ipk_) :: err_act
|
||||
|
||||
integer(psb_ipk_) :: ierr(5)
|
||||
character(len=20) :: name
|
||||
|
||||
name='psb_msort_u'
|
||||
call psb_erractionsave(err_act)
|
||||
|
||||
n = size(x)
|
||||
|
||||
call psb_msort(x,dir=dir)
|
||||
nout = min(1,n)
|
||||
do k=2,n
|
||||
if (x(k) /= x(nout)) then
|
||||
nout = nout + 1
|
||||
x(nout) = x(k)
|
||||
endif
|
||||
enddo
|
||||
|
||||
return
|
||||
|
||||
9999 call psb_error_handler(err_act)
|
||||
|
||||
return
|
||||
end subroutine psb_imsort_u
|
||||
|
||||
|
||||
function psb_ibsrch(key,n,v) result(ipos)
|
||||
use psb_sort_mod, psb_protect_name => psb_ibsrch
|
||||
implicit none
|
||||
integer(psb_ipk_) :: ipos, n
|
||||
integer(psb_ipk_) :: key
|
||||
integer(psb_ipk_) :: 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_ibsrch
|
||||
|
||||
function psb_issrch(key,n,v) result(ipos)
|
||||
use psb_sort_mod, psb_protect_name => psb_issrch
|
||||
implicit none
|
||||
integer(psb_ipk_) :: ipos, n
|
||||
integer(psb_ipk_) :: key
|
||||
integer(psb_ipk_) :: 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_issrch
|
||||
|
||||
subroutine psb_imsort(x,ix,dir,flag)
|
||||
use psb_sort_mod, psb_protect_name => psb_imsort
|
||||
use psb_error_mod
|
||||
use psb_ip_reord_mod
|
||||
implicit none
|
||||
integer(psb_ipk_), intent(inout) :: x(:)
|
||||
integer(psb_ipk_), optional, intent(in) :: dir, flag
|
||||
integer(psb_ipk_), optional, intent(inout) :: ix(:)
|
||||
|
||||
integer(psb_ipk_) :: dir_, flag_, n, err_act
|
||||
|
||||
integer(psb_ipk_), allocatable :: iaux(:)
|
||||
integer(psb_ipk_) :: iret, info, i
|
||||
integer(psb_ipk_) :: ierr(5)
|
||||
character(len=20) :: name
|
||||
|
||||
name='psb_imsort'
|
||||
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_i_msort')
|
||||
goto 9999
|
||||
endif
|
||||
|
||||
select case(dir_)
|
||||
case (psb_sort_up_)
|
||||
call psi_i_msort_up(n,x,iaux,iret)
|
||||
case (psb_sort_down_)
|
||||
call psi_i_msort_dw(n,x,iaux,iret)
|
||||
case (psb_asort_up_)
|
||||
call psi_i_amsort_up(n,x,iaux,iret)
|
||||
case (psb_asort_down_)
|
||||
call psi_i_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_imsort
|
||||
|
||||
subroutine psi_i_msort_up(n,k,l,iret)
|
||||
use psb_const_mod
|
||||
implicit none
|
||||
integer(psb_ipk_) :: n, iret
|
||||
integer(psb_ipk_) :: k(n)
|
||||
integer(psb_ipk_) :: l(0:n+1)
|
||||
!
|
||||
integer(psb_ipk_) :: 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_i_msort_up
|
||||
|
||||
subroutine psi_i_msort_dw(n,k,l,iret)
|
||||
use psb_const_mod
|
||||
implicit none
|
||||
integer(psb_ipk_) :: n, iret
|
||||
integer(psb_ipk_) :: k(n)
|
||||
integer(psb_ipk_) :: l(0:n+1)
|
||||
!
|
||||
integer(psb_ipk_) :: 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_i_msort_dw
|
||||
|
||||
subroutine psi_i_amsort_up(n,k,l,iret)
|
||||
use psb_const_mod
|
||||
implicit none
|
||||
integer(psb_ipk_) :: n, iret
|
||||
integer(psb_ipk_) :: k(n)
|
||||
integer(psb_ipk_) :: l(0:n+1)
|
||||
!
|
||||
integer(psb_ipk_) :: 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_i_amsort_up
|
||||
|
||||
subroutine psi_i_amsort_dw(n,k,l,iret)
|
||||
use psb_const_mod
|
||||
implicit none
|
||||
integer(psb_ipk_) :: n, iret
|
||||
integer(psb_ipk_) :: k(n)
|
||||
integer(psb_ipk_) :: l(0:n+1)
|
||||
!
|
||||
integer(psb_ipk_) :: 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_i_amsort_dw
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
File diff suppressed because it is too large
Load Diff
@ -1,678 +0,0 @@
|
||||
!
|
||||
! 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_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_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_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_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_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_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_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_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
|
||||
|
||||
|
||||
|
||||
|
||||
@ -1,341 +0,0 @@
|
||||
!
|
||||
! 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_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_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_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_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_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_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_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_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_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
|
||||
|
||||
@ -1,713 +0,0 @@
|
||||
!
|
||||
! 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_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_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_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_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_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
@ -0,0 +1,276 @@
|
||||
!
|
||||
! 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.
|
||||
!
|
||||
!
|
||||
! File: psb_c_remote_mat.f90
|
||||
!
|
||||
! Subroutine:
|
||||
! This routine does the retrieval of remote matrix rows.
|
||||
! Retrieval is done through GETROW, therefore it should work
|
||||
! for any matrix format in A; as for the output, default is CSR.
|
||||
!
|
||||
! There is also a specialized version lc_CSR whose interface
|
||||
! is adapted for the needs of c_par_csr_spspmm.
|
||||
!
|
||||
! There are three possible exchange algorithms:
|
||||
! 1. Use MPI_Alltoallv
|
||||
! 2. Use psb_simple_a2av
|
||||
! 3. Use psb_simple_triad_a2av
|
||||
! Default choice is 3. The MPI variant has proved to be inefficient;
|
||||
! that is because it is not persistent, therefore you pay the initialization price
|
||||
! every time, and it is not optimized for a sparse communication pattern,
|
||||
! most MPI implementations assume that all communications are non-empty.
|
||||
! The PSB_SIMPLE variants reuse the same communicator, and go for a simplistic
|
||||
! sequence of sends/receive that is quite efficient for a sparse communication
|
||||
! pattern. To be refined/reviewed in the future to compare with neighbour
|
||||
! persistent collectives.
|
||||
!
|
||||
! Arguments:
|
||||
! a - type(psb_cspmat_type) The local part of input matrix A
|
||||
! desc_a - type(psb_desc_type). The communication descriptor.
|
||||
! blck - type(psb_cspmat_type) The local part of output matrix BLCK
|
||||
! info - integer. Return code
|
||||
! rowcnv - logical Should row/col indices be converted
|
||||
! colcnv - logical to/from global numbering when sent/received?
|
||||
! default is .TRUE.
|
||||
! rowscale - logical Should row/col indices on output be remapped
|
||||
! colscale - logical from MIN:MAX to 1:(MAX-MIN+1) ?
|
||||
! default is .FALSE.
|
||||
! (commmon use is ROWSCALE=.TRUE., COLSCALE=.FALSE.)
|
||||
! data - integer Which index list in desc_a should be used to retrieve
|
||||
! rows, default psb_comm_halo_
|
||||
! psb_comm_halo_ use halo_index
|
||||
! psb_comm_ext_ use ext_index
|
||||
! psb_comm_ovrl_ DISABLED for this routine.
|
||||
!
|
||||
Subroutine psb_lc_remote_mat(a,desc_a,b,info)
|
||||
use psb_base_mod, psb_protect_name => psb_lc_remote_mat
|
||||
|
||||
#ifdef MPI_MOD
|
||||
use mpi
|
||||
#endif
|
||||
Implicit None
|
||||
#ifdef MPI_H
|
||||
include 'mpif.h'
|
||||
#endif
|
||||
|
||||
Type(psb_lc_coo_sparse_mat),Intent(inout) :: a
|
||||
Type(psb_lc_coo_sparse_mat),Intent(inout) :: b
|
||||
Type(psb_desc_type), Intent(inout) :: desc_a
|
||||
integer(psb_ipk_), intent(out) :: info
|
||||
|
||||
! ...local scalars....
|
||||
type(psb_ctxt_type) :: ctxt
|
||||
integer(psb_ipk_) :: np,me
|
||||
integer(psb_ipk_) :: counter, proc, i, n_el_send,n_el_recv, &
|
||||
& n_elem, j, ipx,idxs,idxr
|
||||
integer(psb_lpk_) :: r, k, irmin, irmax, icmin, icmax, iszs, iszr, &
|
||||
& lidx, l1, lnr, lnc, lnnz, idx, ngtz, tot_elem
|
||||
integer(psb_lpk_) :: nz,nouth
|
||||
integer(psb_ipk_) :: nrcvs, nsnds
|
||||
integer(psb_mpk_) :: icomm, minfo
|
||||
integer(psb_mpk_), allocatable :: brvindx(:), &
|
||||
& rvsz(:), bsdindx(:),sdsz(:), sdsi(:), rvsi(:)
|
||||
integer(psb_lpk_), allocatable :: iasnd(:), jasnd(:)
|
||||
complex(psb_spk_), allocatable :: valsnd(:)
|
||||
type(psb_lc_coo_sparse_mat), allocatable :: acoo
|
||||
class(psb_i_base_vect_type), pointer :: pdxv
|
||||
integer(psb_ipk_), allocatable :: ila(:), iprc(:)
|
||||
logical :: rowcnv_,colcnv_,rowscale_,colscale_
|
||||
character(len=5) :: outfmt_
|
||||
integer(psb_ipk_) :: debug_level, debug_unit, err_act
|
||||
character(len=20) :: name, ch_err
|
||||
|
||||
info=psb_success_
|
||||
name='psb_c_remote_mat'
|
||||
call psb_erractionsave(err_act)
|
||||
if (psb_errstatus_fatal()) then
|
||||
info = psb_err_internal_error_ ; goto 9999
|
||||
end if
|
||||
debug_unit = psb_get_debug_unit()
|
||||
debug_level = psb_get_debug_level()
|
||||
|
||||
ctxt = desc_a%get_context()
|
||||
icomm = desc_a%get_mpic()
|
||||
|
||||
Call psb_info(ctxt, me, np)
|
||||
|
||||
if (debug_level >= psb_debug_outer_) &
|
||||
& write(debug_unit,*) me,' ',trim(name),': Start'
|
||||
|
||||
|
||||
call b%free()
|
||||
|
||||
Allocate(rvsz(np),sdsz(np),sdsi(np),rvsi(np),brvindx(np+1),&
|
||||
& bsdindx(np+1), acoo,stat=info)
|
||||
|
||||
if (info /= psb_success_) then
|
||||
info=psb_err_alloc_dealloc_
|
||||
call psb_errpush(info,name)
|
||||
goto 9999
|
||||
end if
|
||||
|
||||
|
||||
nz = a%get_nzeros()
|
||||
!allocate(ila(nz))
|
||||
!write(0,*) me,name,' size :',nz,size(ila)
|
||||
!call desc_a%g2l(a%ia(1:nz),ila(1:nz),info,owned=.false.)
|
||||
!nouth = count(ila(1:nz)<0)
|
||||
!write(0,*) me,name,' Count out of halo :',nouth
|
||||
!call psb_max(ctxt,nouth)
|
||||
!if ((nouth/=0).and.(me==0)) &
|
||||
! & write(0,*) 'Warning: would require reinit of DESC_A'
|
||||
|
||||
call desc_a%indxmap%fnd_owner(a%ia(1:nz),iprc,info)
|
||||
|
||||
icomm = desc_a%get_mpic()
|
||||
sdsz(:)=0
|
||||
rvsz(:)=0
|
||||
sdsi(:)=0
|
||||
rvsi(:)=0
|
||||
ipx = 1
|
||||
brvindx(:) = 0
|
||||
bsdindx(:) = 0
|
||||
counter=1
|
||||
idx = 0
|
||||
idxs = 0
|
||||
idxr = 0
|
||||
do i=1,nz
|
||||
if (iprc(i) >=0) then
|
||||
sdsz(iprc(i)+1) = sdsz(iprc(i)+1) +1
|
||||
else
|
||||
write(0,*)me,name,' Error from fnd_owner: ',iprc(i)
|
||||
end if
|
||||
end do
|
||||
call mpi_alltoall(sdsz,1,psb_mpi_mpk_,&
|
||||
& rvsz,1,psb_mpi_mpk_,icomm,minfo)
|
||||
if (minfo /= psb_success_) then
|
||||
info=psb_err_from_subroutine_
|
||||
call psb_errpush(info,name,a_err='mpi_alltoall')
|
||||
goto 9999
|
||||
end if
|
||||
!write(0,*)me,name,' sdsz ',sdsz(:),' rvsz:',rvsz(:)
|
||||
nsnds = count(sdsz /= 0)
|
||||
nrcvs = count(rvsz /= 0)
|
||||
idxs = 0
|
||||
idxr = 0
|
||||
counter = 1
|
||||
Do proc=0,np-1
|
||||
bsdindx(proc+1) = idxs
|
||||
idxs = idxs + sdsz(proc+1)
|
||||
brvindx(proc+1) = idxr
|
||||
idxr = idxr + rvsz(proc+1)
|
||||
Enddo
|
||||
|
||||
iszs = sum(sdsz)
|
||||
iszr = sum(rvsz)
|
||||
call acoo%allocate(desc_a%get_global_rows(),desc_a%get_global_cols(),iszr)
|
||||
if (psb_errstatus_fatal()) then
|
||||
write(0,*) 'Error from acoo%allocate '
|
||||
info = 4010
|
||||
goto 9999
|
||||
end if
|
||||
if (debug_level >= psb_debug_outer_)&
|
||||
& write(debug_unit,*) me,' ',trim(name),': Sizes:',acoo%get_size(),&
|
||||
& ' Send:',sdsz(:),' Receive:',rvsz(:)
|
||||
!write(debug_unit,*) me,' ',trim(name),': ',info
|
||||
if (info == psb_success_) call psb_ensure_size(max(iszs,1),iasnd,info)
|
||||
!write(debug_unit,*) me,' ',trim(name),' iasnd: ',info
|
||||
if (info == psb_success_) call psb_ensure_size(max(iszs,1),jasnd,info)
|
||||
!write(debug_unit,*) me,' ',trim(name),' jasnd: ',info
|
||||
if (info == psb_success_) call psb_ensure_size(max(iszs,1),valsnd,info)
|
||||
!write(debug_unit,*) me,' ',trim(name),' valsnd: ',info
|
||||
if (info /= psb_success_) then
|
||||
info=psb_err_from_subroutine_
|
||||
call psb_errpush(info,name,a_err='ensure_size')
|
||||
goto 9999
|
||||
end if
|
||||
do k=1, nz
|
||||
proc = iprc(k)
|
||||
sdsi(proc+1) = sdsi(proc+1) + 1
|
||||
!rvsi(proc) = rvsi(proc) + 1
|
||||
iasnd(bsdindx(proc+1)+sdsi(proc+1)) = a%ia(k)
|
||||
jasnd(bsdindx(proc+1)+sdsi(proc+1)) = a%ja(k)
|
||||
valsnd(bsdindx(proc+1)+sdsi(proc+1)) = a%val(k)
|
||||
end do
|
||||
do proc=0,np-1
|
||||
if (sdsi(proc+1) /= sdsz(proc+1)) &
|
||||
& write(0,*) me,name,'Send mismacth ',sdsi(proc+1),sdsz(proc+1)
|
||||
end do
|
||||
|
||||
select case(psb_get_sp_a2av_alg())
|
||||
case(psb_sp_a2av_smpl_triad_)
|
||||
call psb_simple_triad_a2av(valsnd,iasnd,jasnd,sdsz,bsdindx,&
|
||||
& acoo%val,acoo%ia,acoo%ja,rvsz,brvindx,ctxt,info)
|
||||
case(psb_sp_a2av_smpl_v_)
|
||||
call psb_simple_a2av(valsnd,sdsz,bsdindx,&
|
||||
& acoo%val,rvsz,brvindx,ctxt,info)
|
||||
if (info == psb_success_) call psb_simple_a2av(iasnd,sdsz,bsdindx,&
|
||||
& acoo%ia,rvsz,brvindx,ctxt,info)
|
||||
if (info == psb_success_) call psb_simple_a2av(jasnd,sdsz,bsdindx,&
|
||||
& acoo%ja,rvsz,brvindx,ctxt,info)
|
||||
case(psb_sp_a2av_mpi_)
|
||||
|
||||
call mpi_alltoallv(valsnd,sdsz,bsdindx,psb_mpi_c_spk_,&
|
||||
& acoo%val,rvsz,brvindx,psb_mpi_c_spk_,icomm,minfo)
|
||||
if (minfo == mpi_success) &
|
||||
& call mpi_alltoallv(iasnd,sdsz,bsdindx,psb_mpi_lpk_,&
|
||||
& acoo%ia,rvsz,brvindx,psb_mpi_lpk_,icomm,minfo)
|
||||
if (minfo == mpi_success) &
|
||||
& call mpi_alltoallv(jasnd,sdsz,bsdindx,psb_mpi_lpk_,&
|
||||
& acoo%ja,rvsz,brvindx,psb_mpi_lpk_,icomm,minfo)
|
||||
if (minfo /= mpi_success) info = minfo
|
||||
case default
|
||||
info = psb_err_internal_error_
|
||||
call psb_errpush(info,name,a_err='wrong A2AV alg selector')
|
||||
goto 9999
|
||||
end select
|
||||
|
||||
if (info /= psb_success_) then
|
||||
info=psb_err_from_subroutine_
|
||||
call psb_errpush(info,name,a_err='alltoallv')
|
||||
goto 9999
|
||||
end if
|
||||
call acoo%set_nzeros(iszr)
|
||||
call acoo%mv_to_coo(b,info)
|
||||
|
||||
Deallocate(brvindx,bsdindx,rvsz,sdsz,&
|
||||
& iasnd,jasnd,valsnd,stat=info)
|
||||
if (debug_level >= psb_debug_outer_)&
|
||||
& write(debug_unit,*) me,' ',trim(name),': Done'
|
||||
|
||||
call psb_erractionrestore(err_act)
|
||||
return
|
||||
|
||||
9999 call psb_error_handler(ctxt,err_act)
|
||||
|
||||
return
|
||||
|
||||
End Subroutine psb_lc_remote_mat
|
||||
@ -0,0 +1,223 @@
|
||||
!
|
||||
! 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.
|
||||
!
|
||||
!
|
||||
! File: psb_c_remote_vect.f90
|
||||
!
|
||||
! Subroutine:
|
||||
! This routine does the retrieval of remote vector entries.
|
||||
!
|
||||
! There are three possible exchange algorithms:
|
||||
! 1. Use MPI_Alltoallv
|
||||
! 2. Use psb_simple_a2av
|
||||
! 3. Use psb_simple_triad_a2av
|
||||
! Default choice is 3. The MPI variant has proved to be inefficient;
|
||||
! that is because it is not persistent, therefore you pay the initialization price
|
||||
! every time, and it is not optimized for a sparse communication pattern,
|
||||
! most MPI implementations assume that all communications are non-empty.
|
||||
! The PSB_SIMPLE variants reuse the same communicator, and go for a simplistic
|
||||
! sequence of sends/receive that is quite efficient for a sparse communication
|
||||
! pattern. To be refined/reviewed in the future to compare with neighbour
|
||||
! persistent collectives.
|
||||
!
|
||||
! Arguments:
|
||||
! desc_a - type(psb_desc_type). The communication descriptor.
|
||||
! info - integer. Return code
|
||||
! rowcnv - logical Should row/col indices be converted
|
||||
! colcnv - logical to/from global numbering when sent/received?
|
||||
! default is .TRUE.
|
||||
! rowscale - logical Should row/col indices on output be remapped
|
||||
! colscale - logical from MIN:MAX to 1:(MAX-MIN+1) ?
|
||||
! default is .FALSE.
|
||||
! (commmon use is ROWSCALE=.TRUE., COLSCALE=.FALSE.)
|
||||
! data - integer Which index list in desc_a should be used to retrieve
|
||||
! rows, default psb_comm_halo_
|
||||
! psb_comm_halo_ use halo_index
|
||||
! psb_comm_ext_ use ext_index
|
||||
! psb_comm_ovrl_ DISABLED for this routine.
|
||||
!
|
||||
subroutine psb_c_remote_vect(n,v,iv,desc_a,x,ix, info)
|
||||
use psb_base_mod, psb_protect_name => psb_c_remote_vect
|
||||
|
||||
#ifdef MPI_MOD
|
||||
use mpi
|
||||
#endif
|
||||
Implicit None
|
||||
#ifdef MPI_H
|
||||
include 'mpif.h'
|
||||
#endif
|
||||
integer(psb_ipk_), intent(in) :: n
|
||||
complex(psb_spk_), Intent(in) :: v(:)
|
||||
integer(psb_lpk_), Intent(in) :: iv(:)
|
||||
type(psb_desc_type),intent(in) :: desc_a
|
||||
complex(psb_spk_), allocatable, intent(out) :: x(:)
|
||||
integer(psb_lpk_), allocatable, intent(out) :: ix(:)
|
||||
integer(psb_ipk_), intent(out) :: info
|
||||
! ...local scalars....
|
||||
type(psb_ctxt_type) :: ctxt
|
||||
integer(psb_ipk_) :: np,me
|
||||
integer(psb_ipk_) :: counter, proc, i, &
|
||||
& j, idxs,idxr, k, iszs, iszr
|
||||
integer(psb_ipk_) :: nrcvs, nsnds
|
||||
integer(psb_mpk_) :: icomm, minfo
|
||||
integer(psb_mpk_), allocatable :: brvindx(:), &
|
||||
& rvsz(:), bsdindx(:), sdsz(:), sdsi(:), rvsi(:)
|
||||
integer(psb_lpk_), allocatable :: lsnd(:)
|
||||
complex(psb_spk_), allocatable :: valsnd(:)
|
||||
integer(psb_ipk_), allocatable :: iprc(:)
|
||||
integer(psb_ipk_) :: debug_level, debug_unit, err_act
|
||||
character(len=20) :: name, ch_err
|
||||
|
||||
info=psb_success_
|
||||
name='psb_c_remote_vect'
|
||||
call psb_erractionsave(err_act)
|
||||
if (psb_errstatus_fatal()) then
|
||||
info = psb_err_internal_error_ ; goto 9999
|
||||
end if
|
||||
debug_unit = psb_get_debug_unit()
|
||||
debug_level = psb_get_debug_level()
|
||||
|
||||
ctxt = desc_a%get_context()
|
||||
icomm = desc_a%get_mpic()
|
||||
|
||||
Call psb_info(ctxt, me, np)
|
||||
|
||||
if (debug_level >= psb_debug_outer_) &
|
||||
& write(debug_unit,*) me,' ',trim(name),': Start'
|
||||
|
||||
Allocate(rvsz(np),sdsz(np),sdsi(np),rvsi(np),brvindx(np+1),&
|
||||
& bsdindx(np+1), stat=info)
|
||||
|
||||
if (info /= psb_success_) then
|
||||
info=psb_err_alloc_dealloc_
|
||||
call psb_errpush(info,name)
|
||||
goto 9999
|
||||
end if
|
||||
|
||||
call desc_a%indxmap%fnd_owner(iv(1:n),iprc,info)
|
||||
|
||||
icomm = desc_a%get_mpic()
|
||||
sdsz(:) = 0
|
||||
rvsz(:) = 0
|
||||
sdsi(:) = 0
|
||||
rvsi(:) = 0
|
||||
brvindx(:) = 0
|
||||
bsdindx(:) = 0
|
||||
counter = 1
|
||||
idxs = 0
|
||||
idxr = 0
|
||||
do i=1,n
|
||||
if (iprc(i) >=0) then
|
||||
sdsz(iprc(i)+1) = sdsz(iprc(i)+1) +1
|
||||
else
|
||||
write(0,*)me,name,' Error from fnd_owner: ',iprc(i)
|
||||
end if
|
||||
end do
|
||||
call mpi_alltoall(sdsz,1,psb_mpi_mpk_,&
|
||||
& rvsz,1,psb_mpi_mpk_,icomm,minfo)
|
||||
if (minfo /= psb_success_) then
|
||||
info=psb_err_from_subroutine_
|
||||
call psb_errpush(info,name,a_err='mpi_alltoall')
|
||||
goto 9999
|
||||
end if
|
||||
!write(0,*)me,name,' sdsz ',sdsz(:),' rvsz:',rvsz(:)
|
||||
nsnds = count(sdsz /= 0)
|
||||
nrcvs = count(rvsz /= 0)
|
||||
idxs = 0
|
||||
idxr = 0
|
||||
counter = 1
|
||||
Do proc=0,np-1
|
||||
bsdindx(proc+1) = idxs
|
||||
idxs = idxs + sdsz(proc+1)
|
||||
brvindx(proc+1) = idxr
|
||||
idxr = idxr + rvsz(proc+1)
|
||||
Enddo
|
||||
|
||||
iszs = sum(sdsz)
|
||||
iszr = sum(rvsz)
|
||||
call psb_realloc(iszs,lsnd,info)
|
||||
if (info == 0) call psb_realloc(iszs,valsnd,info)
|
||||
if (info == 0) call psb_realloc(iszr,x,info)
|
||||
if (info == 0) call psb_realloc(iszr,ix,info)
|
||||
if (info /= psb_success_) then
|
||||
info=psb_err_from_subroutine_
|
||||
call psb_errpush(info,name,a_err='realloc')
|
||||
goto 9999
|
||||
end if
|
||||
do k=1, n
|
||||
proc = iprc(k)
|
||||
sdsi(proc+1) = sdsi(proc+1) + 1
|
||||
lsnd(bsdindx(proc+1)+sdsi(proc+1)) = iv(k)
|
||||
valsnd(bsdindx(proc+1)+sdsi(proc+1)) = v(k)
|
||||
end do
|
||||
do proc=0,np-1
|
||||
if (sdsi(proc+1) /= sdsz(proc+1)) &
|
||||
& write(0,*) me,name,'Send mismacth ',sdsi(proc+1),sdsz(proc+1)
|
||||
end do
|
||||
|
||||
select case(psb_get_sp_a2av_alg())
|
||||
case(psb_sp_a2av_smpl_triad_,psb_sp_a2av_smpl_v_)
|
||||
call psb_simple_a2av(valsnd,sdsz,bsdindx,&
|
||||
& x,rvsz,brvindx,ctxt,info)
|
||||
if (info == psb_success_) call psb_simple_a2av(lsnd,sdsz,bsdindx,&
|
||||
& ix,rvsz,brvindx,ctxt,info)
|
||||
case(psb_sp_a2av_mpi_)
|
||||
|
||||
call mpi_alltoallv(valsnd,sdsz,bsdindx,psb_mpi_c_spk_,&
|
||||
& x,rvsz,brvindx,psb_mpi_c_spk_,icomm,minfo)
|
||||
if (minfo == mpi_success) &
|
||||
& call mpi_alltoallv(lsnd,sdsz,bsdindx,psb_mpi_lpk_,&
|
||||
& ix,rvsz,brvindx,psb_mpi_lpk_,icomm,minfo)
|
||||
if (minfo /= mpi_success) info = minfo
|
||||
case default
|
||||
info = psb_err_internal_error_
|
||||
call psb_errpush(info,name,a_err='wrong A2AV alg selector')
|
||||
goto 9999
|
||||
end select
|
||||
|
||||
if (info /= psb_success_) then
|
||||
info=psb_err_from_subroutine_
|
||||
call psb_errpush(info,name,a_err='alltoallv')
|
||||
goto 9999
|
||||
end if
|
||||
|
||||
Deallocate(brvindx,bsdindx,rvsz,sdsz,&
|
||||
& lsnd,valsnd,stat=info)
|
||||
if (debug_level >= psb_debug_outer_)&
|
||||
& write(debug_unit,*) me,' ',trim(name),': Done'
|
||||
|
||||
call psb_erractionrestore(err_act)
|
||||
return
|
||||
|
||||
9999 call psb_error_handler(ctxt,err_act)
|
||||
|
||||
return
|
||||
|
||||
End Subroutine psb_c_remote_vect
|
||||
Some files were not shown because too many files have changed in this diff Show More
Loading…
Reference in New Issue