Merge branch 'development' into maint-3.7.0
Salvatore Filippone
3 years ago
commit
0c99e85343
@ -1,610 +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_c_sort_mod
|
||||
use psb_const_mod
|
||||
|
||||
|
||||
@INTE@
|
||||
|
||||
interface psb_msort_unique
|
||||
subroutine psb_cmsort_u(x,nout,dir)
|
||||
import
|
||||
complex(psb_spk_), intent(inout) :: x(:)
|
||||
integer(psb_ipk_), intent(out) :: nout
|
||||
integer(psb_ipk_), optional, intent(in) :: dir
|
||||
end subroutine psb_cmsort_u
|
||||
end interface psb_msort_unique
|
||||
|
||||
type psb_c_heap
|
||||
integer(psb_ipk_) :: last, dir
|
||||
complex(psb_spk_), allocatable :: keys(:)
|
||||
contains
|
||||
procedure, pass(heap) :: init => psb_c_init_heap
|
||||
procedure, pass(heap) :: howmany => psb_c_howmany
|
||||
procedure, pass(heap) :: insert => psb_c_insert_heap
|
||||
procedure, pass(heap) :: get_first => psb_c_heap_get_first
|
||||
procedure, pass(heap) :: dump => psb_c_dump_heap
|
||||
procedure, pass(heap) :: free => psb_c_free_heap
|
||||
end type psb_c_heap
|
||||
|
||||
type psb_c_idx_heap
|
||||
integer(psb_ipk_) :: last, dir
|
||||
complex(psb_spk_), allocatable :: keys(:)
|
||||
integer(psb_ipk_), allocatable :: idxs(:)
|
||||
contains
|
||||
procedure, pass(heap) :: init => psb_c_idx_init_heap
|
||||
procedure, pass(heap) :: howmany => psb_c_idx_howmany
|
||||
procedure, pass(heap) :: insert => psb_c_idx_insert_heap
|
||||
procedure, pass(heap) :: get_first => psb_c_idx_heap_get_first
|
||||
procedure, pass(heap) :: dump => psb_c_idx_dump_heap
|
||||
procedure, pass(heap) :: free => psb_c_idx_free_heap
|
||||
end type psb_c_idx_heap
|
||||
|
||||
|
||||
interface psb_msort
|
||||
subroutine psb_cmsort(x,ix,dir,flag)
|
||||
import
|
||||
complex(psb_spk_), intent(inout) :: x(:)
|
||||
integer(psb_ipk_), optional, intent(in) :: dir, flag
|
||||
integer(psb_ipk_), optional, intent(inout) :: ix(:)
|
||||
end subroutine psb_cmsort
|
||||
end interface psb_msort
|
||||
|
||||
interface
|
||||
subroutine psi_c_lmsort_up(n,k,l,iret)
|
||||
import
|
||||
implicit none
|
||||
integer(psb_ipk_) :: n, iret
|
||||
complex(psb_spk_) :: k(n)
|
||||
integer(psb_ipk_) :: l(0:n+1)
|
||||
end subroutine psi_c_lmsort_up
|
||||
subroutine psi_c_lmsort_dw(n,k,l,iret)
|
||||
import
|
||||
implicit none
|
||||
integer(psb_ipk_) :: n, iret
|
||||
complex(psb_spk_) :: k(n)
|
||||
integer(psb_ipk_) :: l(0:n+1)
|
||||
end subroutine psi_c_lmsort_dw
|
||||
subroutine psi_c_almsort_up(n,k,l,iret)
|
||||
import
|
||||
implicit none
|
||||
integer(psb_ipk_) :: n, iret
|
||||
complex(psb_spk_) :: k(n)
|
||||
integer(psb_ipk_) :: l(0:n+1)
|
||||
end subroutine psi_c_almsort_up
|
||||
subroutine psi_c_almsort_dw(n,k,l,iret)
|
||||
import
|
||||
implicit none
|
||||
integer(psb_ipk_) :: n, iret
|
||||
complex(psb_spk_) :: k(n)
|
||||
integer(psb_ipk_) :: l(0:n+1)
|
||||
end subroutine psi_c_almsort_dw
|
||||
end interface
|
||||
interface
|
||||
subroutine psi_c_amsort_up(n,k,l,iret)
|
||||
import
|
||||
implicit none
|
||||
integer(psb_ipk_) :: n, iret
|
||||
complex(psb_spk_) :: k(n)
|
||||
integer(psb_ipk_) :: l(0:n+1)
|
||||
end subroutine psi_c_amsort_up
|
||||
subroutine psi_c_amsort_dw(n,k,l,iret)
|
||||
import
|
||||
implicit none
|
||||
integer(psb_ipk_) :: n, iret
|
||||
complex(psb_spk_) :: k(n)
|
||||
integer(psb_ipk_) :: l(0:n+1)
|
||||
end subroutine psi_c_amsort_dw
|
||||
end interface
|
||||
|
||||
|
||||
interface psb_qsort
|
||||
subroutine psb_cqsort(x,ix,dir,flag)
|
||||
import
|
||||
complex(psb_spk_), intent(inout) :: x(:)
|
||||
integer(psb_ipk_), optional, intent(in) :: dir, flag
|
||||
integer(psb_ipk_), optional, intent(inout) :: ix(:)
|
||||
end subroutine psb_cqsort
|
||||
end interface psb_qsort
|
||||
|
||||
interface psb_isort
|
||||
subroutine psb_cisort(x,ix,dir,flag)
|
||||
import
|
||||
complex(psb_spk_), intent(inout) :: x(:)
|
||||
integer(psb_ipk_), optional, intent(in) :: dir, flag
|
||||
integer(psb_ipk_), optional, intent(inout) :: ix(:)
|
||||
end subroutine psb_cisort
|
||||
end interface psb_isort
|
||||
|
||||
|
||||
interface psb_hsort
|
||||
subroutine psb_chsort(x,ix,dir,flag)
|
||||
import
|
||||
complex(psb_spk_), intent(inout) :: x(:)
|
||||
integer(psb_ipk_), optional, intent(in) :: dir, flag
|
||||
integer(psb_ipk_), optional, intent(inout) :: ix(:)
|
||||
end subroutine psb_chsort
|
||||
end interface psb_hsort
|
||||
|
||||
|
||||
interface
|
||||
subroutine psi_c_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
|
||||
|
||||
complex(psb_spk_), intent(in) :: key
|
||||
complex(psb_spk_), intent(inout) :: heap(:)
|
||||
integer(psb_ipk_), intent(in) :: dir
|
||||
integer(psb_ipk_), intent(inout) :: last
|
||||
integer(psb_ipk_), intent(out) :: info
|
||||
end subroutine psi_c_insert_heap
|
||||
end interface
|
||||
|
||||
interface
|
||||
subroutine psi_c_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
|
||||
|
||||
complex(psb_spk_), intent(in) :: key
|
||||
complex(psb_spk_), 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_c_idx_insert_heap
|
||||
end interface
|
||||
|
||||
|
||||
interface
|
||||
subroutine psi_c_heap_get_first(key,last,heap,dir,info)
|
||||
import
|
||||
implicit none
|
||||
complex(psb_spk_), intent(inout) :: key
|
||||
integer(psb_ipk_), intent(inout) :: last
|
||||
integer(psb_ipk_), intent(in) :: dir
|
||||
complex(psb_spk_), intent(inout) :: heap(:)
|
||||
integer(psb_ipk_), intent(out) :: info
|
||||
end subroutine psi_c_heap_get_first
|
||||
end interface
|
||||
|
||||
interface
|
||||
subroutine psi_c_idx_heap_get_first(key,index,last,heap,idxs,dir,info)
|
||||
import
|
||||
complex(psb_spk_), intent(inout) :: key
|
||||
integer(psb_ipk_), intent(out) :: index
|
||||
complex(psb_spk_), 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_c_idx_heap_get_first
|
||||
end interface
|
||||
|
||||
interface
|
||||
subroutine psi_clisrx_up(n,x,ix)
|
||||
import
|
||||
complex(psb_spk_), intent(inout) :: x(:)
|
||||
integer(psb_ipk_), intent(inout) :: ix(:)
|
||||
integer(psb_ipk_), intent(in) :: n
|
||||
end subroutine psi_clisrx_up
|
||||
subroutine psi_clisrx_dw(n,x,ix)
|
||||
import
|
||||
complex(psb_spk_), intent(inout) :: x(:)
|
||||
integer(psb_ipk_), intent(inout) :: ix(:)
|
||||
integer(psb_ipk_), intent(in) :: n
|
||||
end subroutine psi_clisrx_dw
|
||||
subroutine psi_clisr_up(n,x)
|
||||
import
|
||||
complex(psb_spk_), intent(inout) :: x(:)
|
||||
integer(psb_ipk_), intent(in) :: n
|
||||
end subroutine psi_clisr_up
|
||||
subroutine psi_clisr_dw(n,x)
|
||||
import
|
||||
complex(psb_spk_), intent(inout) :: x(:)
|
||||
integer(psb_ipk_), intent(in) :: n
|
||||
end subroutine psi_clisr_dw
|
||||
subroutine psi_calisrx_up(n,x,ix)
|
||||
import
|
||||
complex(psb_spk_), intent(inout) :: x(:)
|
||||
integer(psb_ipk_), intent(inout) :: ix(:)
|
||||
integer(psb_ipk_), intent(in) :: n
|
||||
end subroutine psi_calisrx_up
|
||||
subroutine psi_calisrx_dw(n,x,ix)
|
||||
import
|
||||
complex(psb_spk_), intent(inout) :: x(:)
|
||||
integer(psb_ipk_), intent(inout) :: ix(:)
|
||||
integer(psb_ipk_), intent(in) :: n
|
||||
end subroutine psi_calisrx_dw
|
||||
subroutine psi_calisr_up(n,x)
|
||||
import
|
||||
complex(psb_spk_), intent(inout) :: x(:)
|
||||
integer(psb_ipk_), intent(in) :: n
|
||||
end subroutine psi_calisr_up
|
||||
subroutine psi_calisr_dw(n,x)
|
||||
import
|
||||
complex(psb_spk_), intent(inout) :: x(:)
|
||||
integer(psb_ipk_), intent(in) :: n
|
||||
end subroutine psi_calisr_dw
|
||||
subroutine psi_caisrx_up(n,x,ix)
|
||||
import
|
||||
complex(psb_spk_), intent(inout) :: x(:)
|
||||
integer(psb_ipk_), intent(inout) :: ix(:)
|
||||
integer(psb_ipk_), intent(in) :: n
|
||||
end subroutine psi_caisrx_up
|
||||
subroutine psi_caisrx_dw(n,x,ix)
|
||||
import
|
||||
complex(psb_spk_), intent(inout) :: x(:)
|
||||
integer(psb_ipk_), intent(inout) :: ix(:)
|
||||
integer(psb_ipk_), intent(in) :: n
|
||||
end subroutine psi_caisrx_dw
|
||||
subroutine psi_caisr_up(n,x)
|
||||
import
|
||||
complex(psb_spk_), intent(inout) :: x(:)
|
||||
integer(psb_ipk_), intent(in) :: n
|
||||
end subroutine psi_caisr_up
|
||||
subroutine psi_caisr_dw(n,x)
|
||||
import
|
||||
complex(psb_spk_), intent(inout) :: x(:)
|
||||
integer(psb_ipk_), intent(in) :: n
|
||||
end subroutine psi_caisr_dw
|
||||
end interface
|
||||
|
||||
interface
|
||||
subroutine psi_clqsrx_up(n,x,ix)
|
||||
import
|
||||
complex(psb_spk_), intent(inout) :: x(:)
|
||||
integer(psb_ipk_), intent(inout) :: ix(:)
|
||||
integer(psb_ipk_), intent(in) :: n
|
||||
end subroutine psi_clqsrx_up
|
||||
subroutine psi_clqsrx_dw(n,x,ix)
|
||||
import
|
||||
complex(psb_spk_), intent(inout) :: x(:)
|
||||
integer(psb_ipk_), intent(inout) :: ix(:)
|
||||
integer(psb_ipk_), intent(in) :: n
|
||||
end subroutine psi_clqsrx_dw
|
||||
subroutine psi_clqsr_up(n,x)
|
||||
import
|
||||
complex(psb_spk_), intent(inout) :: x(:)
|
||||
integer(psb_ipk_), intent(in) :: n
|
||||
end subroutine psi_clqsr_up
|
||||
subroutine psi_clqsr_dw(n,x)
|
||||
import
|
||||
complex(psb_spk_), intent(inout) :: x(:)
|
||||
integer(psb_ipk_), intent(in) :: n
|
||||
end subroutine psi_clqsr_dw
|
||||
subroutine psi_calqsrx_up(n,x,ix)
|
||||
import
|
||||
complex(psb_spk_), intent(inout) :: x(:)
|
||||
integer(psb_ipk_), intent(inout) :: ix(:)
|
||||
integer(psb_ipk_), intent(in) :: n
|
||||
end subroutine psi_calqsrx_up
|
||||
subroutine psi_calqsrx_dw(n,x,ix)
|
||||
import
|
||||
complex(psb_spk_), intent(inout) :: x(:)
|
||||
integer(psb_ipk_), intent(inout) :: ix(:)
|
||||
integer(psb_ipk_), intent(in) :: n
|
||||
end subroutine psi_calqsrx_dw
|
||||
subroutine psi_calqsr_up(n,x)
|
||||
import
|
||||
complex(psb_spk_), intent(inout) :: x(:)
|
||||
integer(psb_ipk_), intent(in) :: n
|
||||
end subroutine psi_calqsr_up
|
||||
subroutine psi_calqsr_dw(n,x)
|
||||
import
|
||||
complex(psb_spk_), intent(inout) :: x(:)
|
||||
integer(psb_ipk_), intent(in) :: n
|
||||
end subroutine psi_calqsr_dw
|
||||
subroutine psi_caqsrx_up(n,x,ix)
|
||||
import
|
||||
complex(psb_spk_), intent(inout) :: x(:)
|
||||
integer(psb_ipk_), intent(inout) :: ix(:)
|
||||
integer(psb_ipk_), intent(in) :: n
|
||||
end subroutine psi_caqsrx_up
|
||||
subroutine psi_caqsrx_dw(n,x,ix)
|
||||
import
|
||||
complex(psb_spk_), intent(inout) :: x(:)
|
||||
integer(psb_ipk_), intent(inout) :: ix(:)
|
||||
integer(psb_ipk_), intent(in) :: n
|
||||
end subroutine psi_caqsrx_dw
|
||||
subroutine psi_caqsr_up(n,x)
|
||||
import
|
||||
complex(psb_spk_), intent(inout) :: x(:)
|
||||
integer(psb_ipk_), intent(in) :: n
|
||||
end subroutine psi_caqsr_up
|
||||
subroutine psi_caqsr_dw(n,x)
|
||||
import
|
||||
complex(psb_spk_), intent(inout) :: x(:)
|
||||
integer(psb_ipk_), intent(in) :: n
|
||||
end subroutine psi_caqsr_dw
|
||||
end interface
|
||||
|
||||
contains
|
||||
|
||||
subroutine psb_c_init_heap(heap,info,dir)
|
||||
use psb_realloc_mod, only : psb_ensure_size
|
||||
implicit none
|
||||
class(psb_c_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_asort_up_
|
||||
endif
|
||||
select case(heap%dir)
|
||||
case (psb_asort_up_,psb_asort_down_)
|
||||
! ok, do nothing
|
||||
case default
|
||||
write(psb_err_unit,*) 'Invalid direction, defaulting to psb_asort_up_'
|
||||
heap%dir = psb_asort_up_
|
||||
end select
|
||||
call psb_ensure_size(psb_heap_resize,heap%keys,info)
|
||||
|
||||
return
|
||||
end subroutine psb_c_init_heap
|
||||
|
||||
|
||||
function psb_c_howmany(heap) result(res)
|
||||
implicit none
|
||||
class(psb_c_heap), intent(in) :: heap
|
||||
integer(psb_ipk_) :: res
|
||||
res = heap%last
|
||||
end function psb_c_howmany
|
||||
|
||||
subroutine psb_c_insert_heap(key,heap,info)
|
||||
use psb_realloc_mod, only : psb_ensure_size
|
||||
implicit none
|
||||
|
||||
complex(psb_spk_), intent(in) :: key
|
||||
class(psb_c_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_c_insert_heap(key,&
|
||||
& heap%last,heap%keys,heap%dir,info)
|
||||
|
||||
return
|
||||
end subroutine psb_c_insert_heap
|
||||
|
||||
subroutine psb_c_heap_get_first(key,heap,info)
|
||||
implicit none
|
||||
|
||||
class(psb_c_heap), intent(inout) :: heap
|
||||
integer(psb_ipk_), intent(out) :: info
|
||||
complex(psb_spk_), intent(out) :: key
|
||||
|
||||
|
||||
info = psb_success_
|
||||
|
||||
call psi_c_heap_get_first(key,&
|
||||
& heap%last,heap%keys,heap%dir,info)
|
||||
|
||||
return
|
||||
end subroutine psb_c_heap_get_first
|
||||
|
||||
subroutine psb_c_dump_heap(iout,heap,info)
|
||||
|
||||
implicit none
|
||||
class(psb_c_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_c_dump_heap
|
||||
|
||||
subroutine psb_c_free_heap(heap,info)
|
||||
implicit none
|
||||
class(psb_c_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_c_free_heap
|
||||
|
||||
subroutine psb_c_idx_init_heap(heap,info,dir)
|
||||
use psb_realloc_mod, only : psb_ensure_size
|
||||
implicit none
|
||||
class(psb_c_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_asort_up_
|
||||
endif
|
||||
select case(heap%dir)
|
||||
case (psb_asort_up_,psb_asort_down_)
|
||||
! ok, do nothing
|
||||
case default
|
||||
write(psb_err_unit,*) 'Invalid direction, defaulting to psb_asort_up_'
|
||||
heap%dir = psb_asort_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_c_idx_init_heap
|
||||
|
||||
|
||||
function psb_c_idx_howmany(heap) result(res)
|
||||
implicit none
|
||||
class(psb_c_idx_heap), intent(in) :: heap
|
||||
integer(psb_ipk_) :: res
|
||||
res = heap%last
|
||||
end function psb_c_idx_howmany
|
||||
|
||||
subroutine psb_c_idx_insert_heap(key,index,heap,info)
|
||||
use psb_realloc_mod, only : psb_ensure_size
|
||||
implicit none
|
||||
|
||||
complex(psb_spk_), intent(in) :: key
|
||||
integer(psb_ipk_), intent(in) :: index
|
||||
class(psb_c_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_c_idx_insert_heap(key,index,&
|
||||
& heap%last,heap%keys,heap%idxs,heap%dir,info)
|
||||
|
||||
return
|
||||
end subroutine psb_c_idx_insert_heap
|
||||
|
||||
subroutine psb_c_idx_heap_get_first(key,index,heap,info)
|
||||
implicit none
|
||||
|
||||
class(psb_c_idx_heap), intent(inout) :: heap
|
||||
integer(psb_ipk_), intent(out) :: index
|
||||
integer(psb_ipk_), intent(out) :: info
|
||||
complex(psb_spk_), intent(out) :: key
|
||||
|
||||
|
||||
info = psb_success_
|
||||
|
||||
call psi_c_idx_heap_get_first(key,index,&
|
||||
& heap%last,heap%keys,heap%idxs,heap%dir,info)
|
||||
|
||||
return
|
||||
end subroutine psb_c_idx_heap_get_first
|
||||
|
||||
subroutine psb_c_idx_dump_heap(iout,heap,info)
|
||||
|
||||
implicit none
|
||||
class(psb_c_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_c_idx_dump_heap
|
||||
|
||||
subroutine psb_c_idx_free_heap(heap,info)
|
||||
implicit none
|
||||
class(psb_c_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_c_idx_free_heap
|
||||
|
||||
end module psb_c_sort_mod
|
||||
@ -1,572 +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_d_sort_mod
|
||||
use psb_const_mod
|
||||
|
||||
|
||||
@INTE@
|
||||
|
||||
interface psb_msort_unique
|
||||
subroutine psb_dmsort_u(x,nout,dir)
|
||||
import
|
||||
real(psb_dpk_), intent(inout) :: x(:)
|
||||
integer(psb_ipk_), intent(out) :: nout
|
||||
integer(psb_ipk_), optional, intent(in) :: dir
|
||||
end subroutine psb_dmsort_u
|
||||
end interface psb_msort_unique
|
||||
|
||||
type psb_d_heap
|
||||
integer(psb_ipk_) :: last, dir
|
||||
real(psb_dpk_), allocatable :: keys(:)
|
||||
contains
|
||||
procedure, pass(heap) :: init => psb_d_init_heap
|
||||
procedure, pass(heap) :: howmany => psb_d_howmany
|
||||
procedure, pass(heap) :: insert => psb_d_insert_heap
|
||||
procedure, pass(heap) :: get_first => psb_d_heap_get_first
|
||||
procedure, pass(heap) :: dump => psb_d_dump_heap
|
||||
procedure, pass(heap) :: free => psb_d_free_heap
|
||||
end type psb_d_heap
|
||||
|
||||
type psb_d_idx_heap
|
||||
integer(psb_ipk_) :: last, dir
|
||||
real(psb_dpk_), allocatable :: keys(:)
|
||||
integer(psb_ipk_), allocatable :: idxs(:)
|
||||
contains
|
||||
procedure, pass(heap) :: init => psb_d_idx_init_heap
|
||||
procedure, pass(heap) :: howmany => psb_d_idx_howmany
|
||||
procedure, pass(heap) :: insert => psb_d_idx_insert_heap
|
||||
procedure, pass(heap) :: get_first => psb_d_idx_heap_get_first
|
||||
procedure, pass(heap) :: dump => psb_d_idx_dump_heap
|
||||
procedure, pass(heap) :: free => psb_d_idx_free_heap
|
||||
end type psb_d_idx_heap
|
||||
|
||||
|
||||
interface psb_msort
|
||||
subroutine psb_dmsort(x,ix,dir,flag)
|
||||
import
|
||||
real(psb_dpk_), intent(inout) :: x(:)
|
||||
integer(psb_ipk_), optional, intent(in) :: dir, flag
|
||||
integer(psb_ipk_), optional, intent(inout) :: ix(:)
|
||||
end subroutine psb_dmsort
|
||||
end interface psb_msort
|
||||
|
||||
|
||||
interface psb_bsrch
|
||||
function psb_dbsrch(key,n,v) result(ipos)
|
||||
import
|
||||
integer(psb_ipk_) :: ipos, n
|
||||
real(psb_dpk_) :: key
|
||||
real(psb_dpk_) :: v(:)
|
||||
end function psb_dbsrch
|
||||
end interface psb_bsrch
|
||||
|
||||
interface psb_ssrch
|
||||
function psb_dssrch(key,n,v) result(ipos)
|
||||
import
|
||||
implicit none
|
||||
integer(psb_ipk_) :: ipos, n
|
||||
real(psb_dpk_) :: key
|
||||
real(psb_dpk_) :: v(:)
|
||||
end function psb_dssrch
|
||||
end interface psb_ssrch
|
||||
|
||||
interface
|
||||
subroutine psi_d_msort_up(n,k,l,iret)
|
||||
import
|
||||
implicit none
|
||||
integer(psb_ipk_) :: n, iret
|
||||
real(psb_dpk_) :: k(n)
|
||||
integer(psb_ipk_) :: l(0:n+1)
|
||||
end subroutine psi_d_msort_up
|
||||
subroutine psi_d_msort_dw(n,k,l,iret)
|
||||
import
|
||||
implicit none
|
||||
integer(psb_ipk_) :: n, iret
|
||||
real(psb_dpk_) :: k(n)
|
||||
integer(psb_ipk_) :: l(0:n+1)
|
||||
end subroutine psi_d_msort_dw
|
||||
end interface
|
||||
interface
|
||||
subroutine psi_d_amsort_up(n,k,l,iret)
|
||||
import
|
||||
implicit none
|
||||
integer(psb_ipk_) :: n, iret
|
||||
real(psb_dpk_) :: k(n)
|
||||
integer(psb_ipk_) :: l(0:n+1)
|
||||
end subroutine psi_d_amsort_up
|
||||
subroutine psi_d_amsort_dw(n,k,l,iret)
|
||||
import
|
||||
implicit none
|
||||
integer(psb_ipk_) :: n, iret
|
||||
real(psb_dpk_) :: k(n)
|
||||
integer(psb_ipk_) :: l(0:n+1)
|
||||
end subroutine psi_d_amsort_dw
|
||||
end interface
|
||||
|
||||
|
||||
interface psb_qsort
|
||||
subroutine psb_dqsort(x,ix,dir,flag)
|
||||
import
|
||||
real(psb_dpk_), intent(inout) :: x(:)
|
||||
integer(psb_ipk_), optional, intent(in) :: dir, flag
|
||||
integer(psb_ipk_), optional, intent(inout) :: ix(:)
|
||||
end subroutine psb_dqsort
|
||||
end interface psb_qsort
|
||||
|
||||
interface psb_isort
|
||||
subroutine psb_disort(x,ix,dir,flag)
|
||||
import
|
||||
real(psb_dpk_), intent(inout) :: x(:)
|
||||
integer(psb_ipk_), optional, intent(in) :: dir, flag
|
||||
integer(psb_ipk_), optional, intent(inout) :: ix(:)
|
||||
end subroutine psb_disort
|
||||
end interface psb_isort
|
||||
|
||||
|
||||
interface psb_hsort
|
||||
subroutine psb_dhsort(x,ix,dir,flag)
|
||||
import
|
||||
real(psb_dpk_), intent(inout) :: x(:)
|
||||
integer(psb_ipk_), optional, intent(in) :: dir, flag
|
||||
integer(psb_ipk_), optional, intent(inout) :: ix(:)
|
||||
end subroutine psb_dhsort
|
||||
end interface psb_hsort
|
||||
|
||||
|
||||
interface
|
||||
subroutine psi_d_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
|
||||
|
||||
real(psb_dpk_), intent(in) :: key
|
||||
real(psb_dpk_), intent(inout) :: heap(:)
|
||||
integer(psb_ipk_), intent(in) :: dir
|
||||
integer(psb_ipk_), intent(inout) :: last
|
||||
integer(psb_ipk_), intent(out) :: info
|
||||
end subroutine psi_d_insert_heap
|
||||
end interface
|
||||
|
||||
interface
|
||||
subroutine psi_d_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
|
||||
|
||||
real(psb_dpk_), intent(in) :: key
|
||||
real(psb_dpk_), 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_d_idx_insert_heap
|
||||
end interface
|
||||
|
||||
|
||||
interface
|
||||
subroutine psi_d_heap_get_first(key,last,heap,dir,info)
|
||||
import
|
||||
implicit none
|
||||
real(psb_dpk_), intent(inout) :: key
|
||||
integer(psb_ipk_), intent(inout) :: last
|
||||
integer(psb_ipk_), intent(in) :: dir
|
||||
real(psb_dpk_), intent(inout) :: heap(:)
|
||||
integer(psb_ipk_), intent(out) :: info
|
||||
end subroutine psi_d_heap_get_first
|
||||
end interface
|
||||
|
||||
interface
|
||||
subroutine psi_d_idx_heap_get_first(key,index,last,heap,idxs,dir,info)
|
||||
import
|
||||
real(psb_dpk_), intent(inout) :: key
|
||||
integer(psb_ipk_), intent(out) :: index
|
||||
real(psb_dpk_), 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_d_idx_heap_get_first
|
||||
end interface
|
||||
|
||||
interface
|
||||
subroutine psi_disrx_up(n,x,ix)
|
||||
import
|
||||
real(psb_dpk_), intent(inout) :: x(:)
|
||||
integer(psb_ipk_), intent(inout) :: ix(:)
|
||||
integer(psb_ipk_), intent(in) :: n
|
||||
end subroutine psi_disrx_up
|
||||
subroutine psi_disrx_dw(n,x,ix)
|
||||
import
|
||||
real(psb_dpk_), intent(inout) :: x(:)
|
||||
integer(psb_ipk_), intent(inout) :: ix(:)
|
||||
integer(psb_ipk_), intent(in) :: n
|
||||
end subroutine psi_disrx_dw
|
||||
subroutine psi_disr_up(n,x)
|
||||
import
|
||||
real(psb_dpk_), intent(inout) :: x(:)
|
||||
integer(psb_ipk_), intent(in) :: n
|
||||
end subroutine psi_disr_up
|
||||
subroutine psi_disr_dw(n,x)
|
||||
import
|
||||
real(psb_dpk_), intent(inout) :: x(:)
|
||||
integer(psb_ipk_), intent(in) :: n
|
||||
end subroutine psi_disr_dw
|
||||
subroutine psi_daisrx_up(n,x,ix)
|
||||
import
|
||||
real(psb_dpk_), intent(inout) :: x(:)
|
||||
integer(psb_ipk_), intent(inout) :: ix(:)
|
||||
integer(psb_ipk_), intent(in) :: n
|
||||
end subroutine psi_daisrx_up
|
||||
subroutine psi_daisrx_dw(n,x,ix)
|
||||
import
|
||||
real(psb_dpk_), intent(inout) :: x(:)
|
||||
integer(psb_ipk_), intent(inout) :: ix(:)
|
||||
integer(psb_ipk_), intent(in) :: n
|
||||
end subroutine psi_daisrx_dw
|
||||
subroutine psi_daisr_up(n,x)
|
||||
import
|
||||
real(psb_dpk_), intent(inout) :: x(:)
|
||||
integer(psb_ipk_), intent(in) :: n
|
||||
end subroutine psi_daisr_up
|
||||
subroutine psi_daisr_dw(n,x)
|
||||
import
|
||||
real(psb_dpk_), intent(inout) :: x(:)
|
||||
integer(psb_ipk_), intent(in) :: n
|
||||
end subroutine psi_daisr_dw
|
||||
end interface
|
||||
|
||||
interface
|
||||
subroutine psi_dqsrx_up(n,x,ix)
|
||||
import
|
||||
real(psb_dpk_), intent(inout) :: x(:)
|
||||
integer(psb_ipk_), intent(inout) :: ix(:)
|
||||
integer(psb_ipk_), intent(in) :: n
|
||||
end subroutine psi_dqsrx_up
|
||||
subroutine psi_dqsrx_dw(n,x,ix)
|
||||
import
|
||||
real(psb_dpk_), intent(inout) :: x(:)
|
||||
integer(psb_ipk_), intent(inout) :: ix(:)
|
||||
integer(psb_ipk_), intent(in) :: n
|
||||
end subroutine psi_dqsrx_dw
|
||||
subroutine psi_dqsr_up(n,x)
|
||||
import
|
||||
real(psb_dpk_), intent(inout) :: x(:)
|
||||
integer(psb_ipk_), intent(in) :: n
|
||||
end subroutine psi_dqsr_up
|
||||
subroutine psi_dqsr_dw(n,x)
|
||||
import
|
||||
real(psb_dpk_), intent(inout) :: x(:)
|
||||
integer(psb_ipk_), intent(in) :: n
|
||||
end subroutine psi_dqsr_dw
|
||||
subroutine psi_daqsrx_up(n,x,ix)
|
||||
import
|
||||
real(psb_dpk_), intent(inout) :: x(:)
|
||||
integer(psb_ipk_), intent(inout) :: ix(:)
|
||||
integer(psb_ipk_), intent(in) :: n
|
||||
end subroutine psi_daqsrx_up
|
||||
subroutine psi_daqsrx_dw(n,x,ix)
|
||||
import
|
||||
real(psb_dpk_), intent(inout) :: x(:)
|
||||
integer(psb_ipk_), intent(inout) :: ix(:)
|
||||
integer(psb_ipk_), intent(in) :: n
|
||||
end subroutine psi_daqsrx_dw
|
||||
subroutine psi_daqsr_up(n,x)
|
||||
import
|
||||
real(psb_dpk_), intent(inout) :: x(:)
|
||||
integer(psb_ipk_), intent(in) :: n
|
||||
end subroutine psi_daqsr_up
|
||||
subroutine psi_daqsr_dw(n,x)
|
||||
import
|
||||
real(psb_dpk_), intent(inout) :: x(:)
|
||||
integer(psb_ipk_), intent(in) :: n
|
||||
end subroutine psi_daqsr_dw
|
||||
end interface
|
||||
|
||||
contains
|
||||
|
||||
subroutine psb_d_init_heap(heap,info,dir)
|
||||
use psb_realloc_mod, only : psb_ensure_size
|
||||
implicit none
|
||||
class(psb_d_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_d_init_heap
|
||||
|
||||
|
||||
function psb_d_howmany(heap) result(res)
|
||||
implicit none
|
||||
class(psb_d_heap), intent(in) :: heap
|
||||
integer(psb_ipk_) :: res
|
||||
res = heap%last
|
||||
end function psb_d_howmany
|
||||
|
||||
subroutine psb_d_insert_heap(key,heap,info)
|
||||
use psb_realloc_mod, only : psb_ensure_size
|
||||
implicit none
|
||||
|
||||
real(psb_dpk_), intent(in) :: key
|
||||
class(psb_d_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_d_insert_heap(key,&
|
||||
& heap%last,heap%keys,heap%dir,info)
|
||||
|
||||
return
|
||||
end subroutine psb_d_insert_heap
|
||||
|
||||
subroutine psb_d_heap_get_first(key,heap,info)
|
||||
implicit none
|
||||
|
||||
class(psb_d_heap), intent(inout) :: heap
|
||||
integer(psb_ipk_), intent(out) :: info
|
||||
real(psb_dpk_), intent(out) :: key
|
||||
|
||||
|
||||
info = psb_success_
|
||||
|
||||
call psi_d_heap_get_first(key,&
|
||||
& heap%last,heap%keys,heap%dir,info)
|
||||
|
||||
return
|
||||
end subroutine psb_d_heap_get_first
|
||||
|
||||
subroutine psb_d_dump_heap(iout,heap,info)
|
||||
|
||||
implicit none
|
||||
class(psb_d_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_d_dump_heap
|
||||
|
||||
subroutine psb_d_free_heap(heap,info)
|
||||
implicit none
|
||||
class(psb_d_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_d_free_heap
|
||||
|
||||
subroutine psb_d_idx_init_heap(heap,info,dir)
|
||||
use psb_realloc_mod, only : psb_ensure_size
|
||||
implicit none
|
||||
class(psb_d_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_d_idx_init_heap
|
||||
|
||||
|
||||
function psb_d_idx_howmany(heap) result(res)
|
||||
implicit none
|
||||
class(psb_d_idx_heap), intent(in) :: heap
|
||||
integer(psb_ipk_) :: res
|
||||
res = heap%last
|
||||
end function psb_d_idx_howmany
|
||||
|
||||
subroutine psb_d_idx_insert_heap(key,index,heap,info)
|
||||
use psb_realloc_mod, only : psb_ensure_size
|
||||
implicit none
|
||||
|
||||
real(psb_dpk_), intent(in) :: key
|
||||
integer(psb_ipk_), intent(in) :: index
|
||||
class(psb_d_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_d_idx_insert_heap(key,index,&
|
||||
& heap%last,heap%keys,heap%idxs,heap%dir,info)
|
||||
|
||||
return
|
||||
end subroutine psb_d_idx_insert_heap
|
||||
|
||||
subroutine psb_d_idx_heap_get_first(key,index,heap,info)
|
||||
implicit none
|
||||
|
||||
class(psb_d_idx_heap), intent(inout) :: heap
|
||||
integer(psb_ipk_), intent(out) :: index
|
||||
integer(psb_ipk_), intent(out) :: info
|
||||
real(psb_dpk_), intent(out) :: key
|
||||
|
||||
|
||||
info = psb_success_
|
||||
|
||||
call psi_d_idx_heap_get_first(key,index,&
|
||||
& heap%last,heap%keys,heap%idxs,heap%dir,info)
|
||||
|
||||
return
|
||||
end subroutine psb_d_idx_heap_get_first
|
||||
|
||||
subroutine psb_d_idx_dump_heap(iout,heap,info)
|
||||
|
||||
implicit none
|
||||
class(psb_d_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_d_idx_dump_heap
|
||||
|
||||
subroutine psb_d_idx_free_heap(heap,info)
|
||||
implicit none
|
||||
class(psb_d_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_d_idx_free_heap
|
||||
|
||||
end module psb_d_sort_mod
|
||||
@ -1,572 +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_s_sort_mod
|
||||
use psb_const_mod
|
||||
|
||||
|
||||
@INTE@
|
||||
|
||||
interface psb_msort_unique
|
||||
subroutine psb_smsort_u(x,nout,dir)
|
||||
import
|
||||
real(psb_spk_), intent(inout) :: x(:)
|
||||
integer(psb_ipk_), intent(out) :: nout
|
||||
integer(psb_ipk_), optional, intent(in) :: dir
|
||||
end subroutine psb_smsort_u
|
||||
end interface psb_msort_unique
|
||||
|
||||
type psb_s_heap
|
||||
integer(psb_ipk_) :: last, dir
|
||||
real(psb_spk_), allocatable :: keys(:)
|
||||
contains
|
||||
procedure, pass(heap) :: init => psb_s_init_heap
|
||||
procedure, pass(heap) :: howmany => psb_s_howmany
|
||||
procedure, pass(heap) :: insert => psb_s_insert_heap
|
||||
procedure, pass(heap) :: get_first => psb_s_heap_get_first
|
||||
procedure, pass(heap) :: dump => psb_s_dump_heap
|
||||
procedure, pass(heap) :: free => psb_s_free_heap
|
||||
end type psb_s_heap
|
||||
|
||||
type psb_s_idx_heap
|
||||
integer(psb_ipk_) :: last, dir
|
||||
real(psb_spk_), allocatable :: keys(:)
|
||||
integer(psb_ipk_), allocatable :: idxs(:)
|
||||
contains
|
||||
procedure, pass(heap) :: init => psb_s_idx_init_heap
|
||||
procedure, pass(heap) :: howmany => psb_s_idx_howmany
|
||||
procedure, pass(heap) :: insert => psb_s_idx_insert_heap
|
||||
procedure, pass(heap) :: get_first => psb_s_idx_heap_get_first
|
||||
procedure, pass(heap) :: dump => psb_s_idx_dump_heap
|
||||
procedure, pass(heap) :: free => psb_s_idx_free_heap
|
||||
end type psb_s_idx_heap
|
||||
|
||||
|
||||
interface psb_msort
|
||||
subroutine psb_smsort(x,ix,dir,flag)
|
||||
import
|
||||
real(psb_spk_), intent(inout) :: x(:)
|
||||
integer(psb_ipk_), optional, intent(in) :: dir, flag
|
||||
integer(psb_ipk_), optional, intent(inout) :: ix(:)
|
||||
end subroutine psb_smsort
|
||||
end interface psb_msort
|
||||
|
||||
|
||||
interface psb_bsrch
|
||||
function psb_sbsrch(key,n,v) result(ipos)
|
||||
import
|
||||
integer(psb_ipk_) :: ipos, n
|
||||
real(psb_spk_) :: key
|
||||
real(psb_spk_) :: v(:)
|
||||
end function psb_sbsrch
|
||||
end interface psb_bsrch
|
||||
|
||||
interface psb_ssrch
|
||||
function psb_sssrch(key,n,v) result(ipos)
|
||||
import
|
||||
implicit none
|
||||
integer(psb_ipk_) :: ipos, n
|
||||
real(psb_spk_) :: key
|
||||
real(psb_spk_) :: v(:)
|
||||
end function psb_sssrch
|
||||
end interface psb_ssrch
|
||||
|
||||
interface
|
||||
subroutine psi_s_msort_up(n,k,l,iret)
|
||||
import
|
||||
implicit none
|
||||
integer(psb_ipk_) :: n, iret
|
||||
real(psb_spk_) :: k(n)
|
||||
integer(psb_ipk_) :: l(0:n+1)
|
||||
end subroutine psi_s_msort_up
|
||||
subroutine psi_s_msort_dw(n,k,l,iret)
|
||||
import
|
||||
implicit none
|
||||
integer(psb_ipk_) :: n, iret
|
||||
real(psb_spk_) :: k(n)
|
||||
integer(psb_ipk_) :: l(0:n+1)
|
||||
end subroutine psi_s_msort_dw
|
||||
end interface
|
||||
interface
|
||||
subroutine psi_s_amsort_up(n,k,l,iret)
|
||||
import
|
||||
implicit none
|
||||
integer(psb_ipk_) :: n, iret
|
||||
real(psb_spk_) :: k(n)
|
||||
integer(psb_ipk_) :: l(0:n+1)
|
||||
end subroutine psi_s_amsort_up
|
||||
subroutine psi_s_amsort_dw(n,k,l,iret)
|
||||
import
|
||||
implicit none
|
||||
integer(psb_ipk_) :: n, iret
|
||||
real(psb_spk_) :: k(n)
|
||||
integer(psb_ipk_) :: l(0:n+1)
|
||||
end subroutine psi_s_amsort_dw
|
||||
end interface
|
||||
|
||||
|
||||
interface psb_qsort
|
||||
subroutine psb_sqsort(x,ix,dir,flag)
|
||||
import
|
||||
real(psb_spk_), intent(inout) :: x(:)
|
||||
integer(psb_ipk_), optional, intent(in) :: dir, flag
|
||||
integer(psb_ipk_), optional, intent(inout) :: ix(:)
|
||||
end subroutine psb_sqsort
|
||||
end interface psb_qsort
|
||||
|
||||
interface psb_isort
|
||||
subroutine psb_sisort(x,ix,dir,flag)
|
||||
import
|
||||
real(psb_spk_), intent(inout) :: x(:)
|
||||
integer(psb_ipk_), optional, intent(in) :: dir, flag
|
||||
integer(psb_ipk_), optional, intent(inout) :: ix(:)
|
||||
end subroutine psb_sisort
|
||||
end interface psb_isort
|
||||
|
||||
|
||||
interface psb_hsort
|
||||
subroutine psb_shsort(x,ix,dir,flag)
|
||||
import
|
||||
real(psb_spk_), intent(inout) :: x(:)
|
||||
integer(psb_ipk_), optional, intent(in) :: dir, flag
|
||||
integer(psb_ipk_), optional, intent(inout) :: ix(:)
|
||||
end subroutine psb_shsort
|
||||
end interface psb_hsort
|
||||
|
||||
|
||||
interface
|
||||
subroutine psi_s_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
|
||||
|
||||
real(psb_spk_), intent(in) :: key
|
||||
real(psb_spk_), intent(inout) :: heap(:)
|
||||
integer(psb_ipk_), intent(in) :: dir
|
||||
integer(psb_ipk_), intent(inout) :: last
|
||||
integer(psb_ipk_), intent(out) :: info
|
||||
end subroutine psi_s_insert_heap
|
||||
end interface
|
||||
|
||||
interface
|
||||
subroutine psi_s_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
|
||||
|
||||
real(psb_spk_), intent(in) :: key
|
||||
real(psb_spk_), 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_s_idx_insert_heap
|
||||
end interface
|
||||
|
||||
|
||||
interface
|
||||
subroutine psi_s_heap_get_first(key,last,heap,dir,info)
|
||||
import
|
||||
implicit none
|
||||
real(psb_spk_), intent(inout) :: key
|
||||
integer(psb_ipk_), intent(inout) :: last
|
||||
integer(psb_ipk_), intent(in) :: dir
|
||||
real(psb_spk_), intent(inout) :: heap(:)
|
||||
integer(psb_ipk_), intent(out) :: info
|
||||
end subroutine psi_s_heap_get_first
|
||||
end interface
|
||||
|
||||
interface
|
||||
subroutine psi_s_idx_heap_get_first(key,index,last,heap,idxs,dir,info)
|
||||
import
|
||||
real(psb_spk_), intent(inout) :: key
|
||||
integer(psb_ipk_), intent(out) :: index
|
||||
real(psb_spk_), 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_s_idx_heap_get_first
|
||||
end interface
|
||||
|
||||
interface
|
||||
subroutine psi_sisrx_up(n,x,ix)
|
||||
import
|
||||
real(psb_spk_), intent(inout) :: x(:)
|
||||
integer(psb_ipk_), intent(inout) :: ix(:)
|
||||
integer(psb_ipk_), intent(in) :: n
|
||||
end subroutine psi_sisrx_up
|
||||
subroutine psi_sisrx_dw(n,x,ix)
|
||||
import
|
||||
real(psb_spk_), intent(inout) :: x(:)
|
||||
integer(psb_ipk_), intent(inout) :: ix(:)
|
||||
integer(psb_ipk_), intent(in) :: n
|
||||
end subroutine psi_sisrx_dw
|
||||
subroutine psi_sisr_up(n,x)
|
||||
import
|
||||
real(psb_spk_), intent(inout) :: x(:)
|
||||
integer(psb_ipk_), intent(in) :: n
|
||||
end subroutine psi_sisr_up
|
||||
subroutine psi_sisr_dw(n,x)
|
||||
import
|
||||
real(psb_spk_), intent(inout) :: x(:)
|
||||
integer(psb_ipk_), intent(in) :: n
|
||||
end subroutine psi_sisr_dw
|
||||
subroutine psi_saisrx_up(n,x,ix)
|
||||
import
|
||||
real(psb_spk_), intent(inout) :: x(:)
|
||||
integer(psb_ipk_), intent(inout) :: ix(:)
|
||||
integer(psb_ipk_), intent(in) :: n
|
||||
end subroutine psi_saisrx_up
|
||||
subroutine psi_saisrx_dw(n,x,ix)
|
||||
import
|
||||
real(psb_spk_), intent(inout) :: x(:)
|
||||
integer(psb_ipk_), intent(inout) :: ix(:)
|
||||
integer(psb_ipk_), intent(in) :: n
|
||||
end subroutine psi_saisrx_dw
|
||||
subroutine psi_saisr_up(n,x)
|
||||
import
|
||||
real(psb_spk_), intent(inout) :: x(:)
|
||||
integer(psb_ipk_), intent(in) :: n
|
||||
end subroutine psi_saisr_up
|
||||
subroutine psi_saisr_dw(n,x)
|
||||
import
|
||||
real(psb_spk_), intent(inout) :: x(:)
|
||||
integer(psb_ipk_), intent(in) :: n
|
||||
end subroutine psi_saisr_dw
|
||||
end interface
|
||||
|
||||
interface
|
||||
subroutine psi_sqsrx_up(n,x,ix)
|
||||
import
|
||||
real(psb_spk_), intent(inout) :: x(:)
|
||||
integer(psb_ipk_), intent(inout) :: ix(:)
|
||||
integer(psb_ipk_), intent(in) :: n
|
||||
end subroutine psi_sqsrx_up
|
||||
subroutine psi_sqsrx_dw(n,x,ix)
|
||||
import
|
||||
real(psb_spk_), intent(inout) :: x(:)
|
||||
integer(psb_ipk_), intent(inout) :: ix(:)
|
||||
integer(psb_ipk_), intent(in) :: n
|
||||
end subroutine psi_sqsrx_dw
|
||||
subroutine psi_sqsr_up(n,x)
|
||||
import
|
||||
real(psb_spk_), intent(inout) :: x(:)
|
||||
integer(psb_ipk_), intent(in) :: n
|
||||
end subroutine psi_sqsr_up
|
||||
subroutine psi_sqsr_dw(n,x)
|
||||
import
|
||||
real(psb_spk_), intent(inout) :: x(:)
|
||||
integer(psb_ipk_), intent(in) :: n
|
||||
end subroutine psi_sqsr_dw
|
||||
subroutine psi_saqsrx_up(n,x,ix)
|
||||
import
|
||||
real(psb_spk_), intent(inout) :: x(:)
|
||||
integer(psb_ipk_), intent(inout) :: ix(:)
|
||||
integer(psb_ipk_), intent(in) :: n
|
||||
end subroutine psi_saqsrx_up
|
||||
subroutine psi_saqsrx_dw(n,x,ix)
|
||||
import
|
||||
real(psb_spk_), intent(inout) :: x(:)
|
||||
integer(psb_ipk_), intent(inout) :: ix(:)
|
||||
integer(psb_ipk_), intent(in) :: n
|
||||
end subroutine psi_saqsrx_dw
|
||||
subroutine psi_saqsr_up(n,x)
|
||||
import
|
||||
real(psb_spk_), intent(inout) :: x(:)
|
||||
integer(psb_ipk_), intent(in) :: n
|
||||
end subroutine psi_saqsr_up
|
||||
subroutine psi_saqsr_dw(n,x)
|
||||
import
|
||||
real(psb_spk_), intent(inout) :: x(:)
|
||||
integer(psb_ipk_), intent(in) :: n
|
||||
end subroutine psi_saqsr_dw
|
||||
end interface
|
||||
|
||||
contains
|
||||
|
||||
subroutine psb_s_init_heap(heap,info,dir)
|
||||
use psb_realloc_mod, only : psb_ensure_size
|
||||
implicit none
|
||||
class(psb_s_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_s_init_heap
|
||||
|
||||
|
||||
function psb_s_howmany(heap) result(res)
|
||||
implicit none
|
||||
class(psb_s_heap), intent(in) :: heap
|
||||
integer(psb_ipk_) :: res
|
||||
res = heap%last
|
||||
end function psb_s_howmany
|
||||
|
||||
subroutine psb_s_insert_heap(key,heap,info)
|
||||
use psb_realloc_mod, only : psb_ensure_size
|
||||
implicit none
|
||||
|
||||
real(psb_spk_), intent(in) :: key
|
||||
class(psb_s_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_s_insert_heap(key,&
|
||||
& heap%last,heap%keys,heap%dir,info)
|
||||
|
||||
return
|
||||
end subroutine psb_s_insert_heap
|
||||
|
||||
subroutine psb_s_heap_get_first(key,heap,info)
|
||||
implicit none
|
||||
|
||||
class(psb_s_heap), intent(inout) :: heap
|
||||
integer(psb_ipk_), intent(out) :: info
|
||||
real(psb_spk_), intent(out) :: key
|
||||
|
||||
|
||||
info = psb_success_
|
||||
|
||||
call psi_s_heap_get_first(key,&
|
||||
& heap%last,heap%keys,heap%dir,info)
|
||||
|
||||
return
|
||||
end subroutine psb_s_heap_get_first
|
||||
|
||||
subroutine psb_s_dump_heap(iout,heap,info)
|
||||
|
||||
implicit none
|
||||
class(psb_s_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_s_dump_heap
|
||||
|
||||
subroutine psb_s_free_heap(heap,info)
|
||||
implicit none
|
||||
class(psb_s_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_s_free_heap
|
||||
|
||||
subroutine psb_s_idx_init_heap(heap,info,dir)
|
||||
use psb_realloc_mod, only : psb_ensure_size
|
||||
implicit none
|
||||
class(psb_s_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_s_idx_init_heap
|
||||
|
||||
|
||||
function psb_s_idx_howmany(heap) result(res)
|
||||
implicit none
|
||||
class(psb_s_idx_heap), intent(in) :: heap
|
||||
integer(psb_ipk_) :: res
|
||||
res = heap%last
|
||||
end function psb_s_idx_howmany
|
||||
|
||||
subroutine psb_s_idx_insert_heap(key,index,heap,info)
|
||||
use psb_realloc_mod, only : psb_ensure_size
|
||||
implicit none
|
||||
|
||||
real(psb_spk_), intent(in) :: key
|
||||
integer(psb_ipk_), intent(in) :: index
|
||||
class(psb_s_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_s_idx_insert_heap(key,index,&
|
||||
& heap%last,heap%keys,heap%idxs,heap%dir,info)
|
||||
|
||||
return
|
||||
end subroutine psb_s_idx_insert_heap
|
||||
|
||||
subroutine psb_s_idx_heap_get_first(key,index,heap,info)
|
||||
implicit none
|
||||
|
||||
class(psb_s_idx_heap), intent(inout) :: heap
|
||||
integer(psb_ipk_), intent(out) :: index
|
||||
integer(psb_ipk_), intent(out) :: info
|
||||
real(psb_spk_), intent(out) :: key
|
||||
|
||||
|
||||
info = psb_success_
|
||||
|
||||
call psi_s_idx_heap_get_first(key,index,&
|
||||
& heap%last,heap%keys,heap%idxs,heap%dir,info)
|
||||
|
||||
return
|
||||
end subroutine psb_s_idx_heap_get_first
|
||||
|
||||
subroutine psb_s_idx_dump_heap(iout,heap,info)
|
||||
|
||||
implicit none
|
||||
class(psb_s_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_s_idx_dump_heap
|
||||
|
||||
subroutine psb_s_idx_free_heap(heap,info)
|
||||
implicit none
|
||||
class(psb_s_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_s_idx_free_heap
|
||||
|
||||
end module psb_s_sort_mod
|
||||
@ -1,610 +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_z_sort_mod
|
||||
use psb_const_mod
|
||||
|
||||
|
||||
@INTE@
|
||||
|
||||
interface psb_msort_unique
|
||||
subroutine psb_zmsort_u(x,nout,dir)
|
||||
import
|
||||
complex(psb_dpk_), intent(inout) :: x(:)
|
||||
integer(psb_ipk_), intent(out) :: nout
|
||||
integer(psb_ipk_), optional, intent(in) :: dir
|
||||
end subroutine psb_zmsort_u
|
||||
end interface psb_msort_unique
|
||||
|
||||
type psb_z_heap
|
||||
integer(psb_ipk_) :: last, dir
|
||||
complex(psb_dpk_), allocatable :: keys(:)
|
||||
contains
|
||||
procedure, pass(heap) :: init => psb_z_init_heap
|
||||
procedure, pass(heap) :: howmany => psb_z_howmany
|
||||
procedure, pass(heap) :: insert => psb_z_insert_heap
|
||||
procedure, pass(heap) :: get_first => psb_z_heap_get_first
|
||||
procedure, pass(heap) :: dump => psb_z_dump_heap
|
||||
procedure, pass(heap) :: free => psb_z_free_heap
|
||||
end type psb_z_heap
|
||||
|
||||
type psb_z_idx_heap
|
||||
integer(psb_ipk_) :: last, dir
|
||||
complex(psb_dpk_), allocatable :: keys(:)
|
||||
integer(psb_ipk_), allocatable :: idxs(:)
|
||||
contains
|
||||
procedure, pass(heap) :: init => psb_z_idx_init_heap
|
||||
procedure, pass(heap) :: howmany => psb_z_idx_howmany
|
||||
procedure, pass(heap) :: insert => psb_z_idx_insert_heap
|
||||
procedure, pass(heap) :: get_first => psb_z_idx_heap_get_first
|
||||
procedure, pass(heap) :: dump => psb_z_idx_dump_heap
|
||||
procedure, pass(heap) :: free => psb_z_idx_free_heap
|
||||
end type psb_z_idx_heap
|
||||
|
||||
|
||||
interface psb_msort
|
||||
subroutine psb_zmsort(x,ix,dir,flag)
|
||||
import
|
||||
complex(psb_dpk_), intent(inout) :: x(:)
|
||||
integer(psb_ipk_), optional, intent(in) :: dir, flag
|
||||
integer(psb_ipk_), optional, intent(inout) :: ix(:)
|
||||
end subroutine psb_zmsort
|
||||
end interface psb_msort
|
||||
|
||||
interface
|
||||
subroutine psi_z_lmsort_up(n,k,l,iret)
|
||||
import
|
||||
implicit none
|
||||
integer(psb_ipk_) :: n, iret
|
||||
complex(psb_dpk_) :: k(n)
|
||||
integer(psb_ipk_) :: l(0:n+1)
|
||||
end subroutine psi_z_lmsort_up
|
||||
subroutine psi_z_lmsort_dw(n,k,l,iret)
|
||||
import
|
||||
implicit none
|
||||
integer(psb_ipk_) :: n, iret
|
||||
complex(psb_dpk_) :: k(n)
|
||||
integer(psb_ipk_) :: l(0:n+1)
|
||||
end subroutine psi_z_lmsort_dw
|
||||
subroutine psi_z_almsort_up(n,k,l,iret)
|
||||
import
|
||||
implicit none
|
||||
integer(psb_ipk_) :: n, iret
|
||||
complex(psb_dpk_) :: k(n)
|
||||
integer(psb_ipk_) :: l(0:n+1)
|
||||
end subroutine psi_z_almsort_up
|
||||
subroutine psi_z_almsort_dw(n,k,l,iret)
|
||||
import
|
||||
implicit none
|
||||
integer(psb_ipk_) :: n, iret
|
||||
complex(psb_dpk_) :: k(n)
|
||||
integer(psb_ipk_) :: l(0:n+1)
|
||||
end subroutine psi_z_almsort_dw
|
||||
end interface
|
||||
interface
|
||||
subroutine psi_z_amsort_up(n,k,l,iret)
|
||||
import
|
||||
implicit none
|
||||
integer(psb_ipk_) :: n, iret
|
||||
complex(psb_dpk_) :: k(n)
|
||||
integer(psb_ipk_) :: l(0:n+1)
|
||||
end subroutine psi_z_amsort_up
|
||||
subroutine psi_z_amsort_dw(n,k,l,iret)
|
||||
import
|
||||
implicit none
|
||||
integer(psb_ipk_) :: n, iret
|
||||
complex(psb_dpk_) :: k(n)
|
||||
integer(psb_ipk_) :: l(0:n+1)
|
||||
end subroutine psi_z_amsort_dw
|
||||
end interface
|
||||
|
||||
|
||||
interface psb_qsort
|
||||
subroutine psb_zqsort(x,ix,dir,flag)
|
||||
import
|
||||
complex(psb_dpk_), intent(inout) :: x(:)
|
||||
integer(psb_ipk_), optional, intent(in) :: dir, flag
|
||||
integer(psb_ipk_), optional, intent(inout) :: ix(:)
|
||||
end subroutine psb_zqsort
|
||||
end interface psb_qsort
|
||||
|
||||
interface psb_isort
|
||||
subroutine psb_zisort(x,ix,dir,flag)
|
||||
import
|
||||
complex(psb_dpk_), intent(inout) :: x(:)
|
||||
integer(psb_ipk_), optional, intent(in) :: dir, flag
|
||||
integer(psb_ipk_), optional, intent(inout) :: ix(:)
|
||||
end subroutine psb_zisort
|
||||
end interface psb_isort
|
||||
|
||||
|
||||
interface psb_hsort
|
||||
subroutine psb_zhsort(x,ix,dir,flag)
|
||||
import
|
||||
complex(psb_dpk_), intent(inout) :: x(:)
|
||||
integer(psb_ipk_), optional, intent(in) :: dir, flag
|
||||
integer(psb_ipk_), optional, intent(inout) :: ix(:)
|
||||
end subroutine psb_zhsort
|
||||
end interface psb_hsort
|
||||
|
||||
|
||||
interface
|
||||
subroutine psi_z_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
|
||||
|
||||
complex(psb_dpk_), intent(in) :: key
|
||||
complex(psb_dpk_), intent(inout) :: heap(:)
|
||||
integer(psb_ipk_), intent(in) :: dir
|
||||
integer(psb_ipk_), intent(inout) :: last
|
||||
integer(psb_ipk_), intent(out) :: info
|
||||
end subroutine psi_z_insert_heap
|
||||
end interface
|
||||
|
||||
interface
|
||||
subroutine psi_z_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
|
||||
|
||||
complex(psb_dpk_), intent(in) :: key
|
||||
complex(psb_dpk_), 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_z_idx_insert_heap
|
||||
end interface
|
||||
|
||||
|
||||
interface
|
||||
subroutine psi_z_heap_get_first(key,last,heap,dir,info)
|
||||
import
|
||||
implicit none
|
||||
complex(psb_dpk_), intent(inout) :: key
|
||||
integer(psb_ipk_), intent(inout) :: last
|
||||
integer(psb_ipk_), intent(in) :: dir
|
||||
complex(psb_dpk_), intent(inout) :: heap(:)
|
||||
integer(psb_ipk_), intent(out) :: info
|
||||
end subroutine psi_z_heap_get_first
|
||||
end interface
|
||||
|
||||
interface
|
||||
subroutine psi_z_idx_heap_get_first(key,index,last,heap,idxs,dir,info)
|
||||
import
|
||||
complex(psb_dpk_), intent(inout) :: key
|
||||
integer(psb_ipk_), intent(out) :: index
|
||||
complex(psb_dpk_), 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_z_idx_heap_get_first
|
||||
end interface
|
||||
|
||||
interface
|
||||
subroutine psi_zlisrx_up(n,x,ix)
|
||||
import
|
||||
complex(psb_dpk_), intent(inout) :: x(:)
|
||||
integer(psb_ipk_), intent(inout) :: ix(:)
|
||||
integer(psb_ipk_), intent(in) :: n
|
||||
end subroutine psi_zlisrx_up
|
||||
subroutine psi_zlisrx_dw(n,x,ix)
|
||||
import
|
||||
complex(psb_dpk_), intent(inout) :: x(:)
|
||||
integer(psb_ipk_), intent(inout) :: ix(:)
|
||||
integer(psb_ipk_), intent(in) :: n
|
||||
end subroutine psi_zlisrx_dw
|
||||
subroutine psi_zlisr_up(n,x)
|
||||
import
|
||||
complex(psb_dpk_), intent(inout) :: x(:)
|
||||
integer(psb_ipk_), intent(in) :: n
|
||||
end subroutine psi_zlisr_up
|
||||
subroutine psi_zlisr_dw(n,x)
|
||||
import
|
||||
complex(psb_dpk_), intent(inout) :: x(:)
|
||||
integer(psb_ipk_), intent(in) :: n
|
||||
end subroutine psi_zlisr_dw
|
||||
subroutine psi_zalisrx_up(n,x,ix)
|
||||
import
|
||||
complex(psb_dpk_), intent(inout) :: x(:)
|
||||
integer(psb_ipk_), intent(inout) :: ix(:)
|
||||
integer(psb_ipk_), intent(in) :: n
|
||||
end subroutine psi_zalisrx_up
|
||||
subroutine psi_zalisrx_dw(n,x,ix)
|
||||
import
|
||||
complex(psb_dpk_), intent(inout) :: x(:)
|
||||
integer(psb_ipk_), intent(inout) :: ix(:)
|
||||
integer(psb_ipk_), intent(in) :: n
|
||||
end subroutine psi_zalisrx_dw
|
||||
subroutine psi_zalisr_up(n,x)
|
||||
import
|
||||
complex(psb_dpk_), intent(inout) :: x(:)
|
||||
integer(psb_ipk_), intent(in) :: n
|
||||
end subroutine psi_zalisr_up
|
||||
subroutine psi_zalisr_dw(n,x)
|
||||
import
|
||||
complex(psb_dpk_), intent(inout) :: x(:)
|
||||
integer(psb_ipk_), intent(in) :: n
|
||||
end subroutine psi_zalisr_dw
|
||||
subroutine psi_zaisrx_up(n,x,ix)
|
||||
import
|
||||
complex(psb_dpk_), intent(inout) :: x(:)
|
||||
integer(psb_ipk_), intent(inout) :: ix(:)
|
||||
integer(psb_ipk_), intent(in) :: n
|
||||
end subroutine psi_zaisrx_up
|
||||
subroutine psi_zaisrx_dw(n,x,ix)
|
||||
import
|
||||
complex(psb_dpk_), intent(inout) :: x(:)
|
||||
integer(psb_ipk_), intent(inout) :: ix(:)
|
||||
integer(psb_ipk_), intent(in) :: n
|
||||
end subroutine psi_zaisrx_dw
|
||||
subroutine psi_zaisr_up(n,x)
|
||||
import
|
||||
complex(psb_dpk_), intent(inout) :: x(:)
|
||||
integer(psb_ipk_), intent(in) :: n
|
||||
end subroutine psi_zaisr_up
|
||||
subroutine psi_zaisr_dw(n,x)
|
||||
import
|
||||
complex(psb_dpk_), intent(inout) :: x(:)
|
||||
integer(psb_ipk_), intent(in) :: n
|
||||
end subroutine psi_zaisr_dw
|
||||
end interface
|
||||
|
||||
interface
|
||||
subroutine psi_zlqsrx_up(n,x,ix)
|
||||
import
|
||||
complex(psb_dpk_), intent(inout) :: x(:)
|
||||
integer(psb_ipk_), intent(inout) :: ix(:)
|
||||
integer(psb_ipk_), intent(in) :: n
|
||||
end subroutine psi_zlqsrx_up
|
||||
subroutine psi_zlqsrx_dw(n,x,ix)
|
||||
import
|
||||
complex(psb_dpk_), intent(inout) :: x(:)
|
||||
integer(psb_ipk_), intent(inout) :: ix(:)
|
||||
integer(psb_ipk_), intent(in) :: n
|
||||
end subroutine psi_zlqsrx_dw
|
||||
subroutine psi_zlqsr_up(n,x)
|
||||
import
|
||||
complex(psb_dpk_), intent(inout) :: x(:)
|
||||
integer(psb_ipk_), intent(in) :: n
|
||||
end subroutine psi_zlqsr_up
|
||||
subroutine psi_zlqsr_dw(n,x)
|
||||
import
|
||||
complex(psb_dpk_), intent(inout) :: x(:)
|
||||
integer(psb_ipk_), intent(in) :: n
|
||||
end subroutine psi_zlqsr_dw
|
||||
subroutine psi_zalqsrx_up(n,x,ix)
|
||||
import
|
||||
complex(psb_dpk_), intent(inout) :: x(:)
|
||||
integer(psb_ipk_), intent(inout) :: ix(:)
|
||||
integer(psb_ipk_), intent(in) :: n
|
||||
end subroutine psi_zalqsrx_up
|
||||
subroutine psi_zalqsrx_dw(n,x,ix)
|
||||
import
|
||||
complex(psb_dpk_), intent(inout) :: x(:)
|
||||
integer(psb_ipk_), intent(inout) :: ix(:)
|
||||
integer(psb_ipk_), intent(in) :: n
|
||||
end subroutine psi_zalqsrx_dw
|
||||
subroutine psi_zalqsr_up(n,x)
|
||||
import
|
||||
complex(psb_dpk_), intent(inout) :: x(:)
|
||||
integer(psb_ipk_), intent(in) :: n
|
||||
end subroutine psi_zalqsr_up
|
||||
subroutine psi_zalqsr_dw(n,x)
|
||||
import
|
||||
complex(psb_dpk_), intent(inout) :: x(:)
|
||||
integer(psb_ipk_), intent(in) :: n
|
||||
end subroutine psi_zalqsr_dw
|
||||
subroutine psi_zaqsrx_up(n,x,ix)
|
||||
import
|
||||
complex(psb_dpk_), intent(inout) :: x(:)
|
||||
integer(psb_ipk_), intent(inout) :: ix(:)
|
||||
integer(psb_ipk_), intent(in) :: n
|
||||
end subroutine psi_zaqsrx_up
|
||||
subroutine psi_zaqsrx_dw(n,x,ix)
|
||||
import
|
||||
complex(psb_dpk_), intent(inout) :: x(:)
|
||||
integer(psb_ipk_), intent(inout) :: ix(:)
|
||||
integer(psb_ipk_), intent(in) :: n
|
||||
end subroutine psi_zaqsrx_dw
|
||||
subroutine psi_zaqsr_up(n,x)
|
||||
import
|
||||
complex(psb_dpk_), intent(inout) :: x(:)
|
||||
integer(psb_ipk_), intent(in) :: n
|
||||
end subroutine psi_zaqsr_up
|
||||
subroutine psi_zaqsr_dw(n,x)
|
||||
import
|
||||
complex(psb_dpk_), intent(inout) :: x(:)
|
||||
integer(psb_ipk_), intent(in) :: n
|
||||
end subroutine psi_zaqsr_dw
|
||||
end interface
|
||||
|
||||
contains
|
||||
|
||||
subroutine psb_z_init_heap(heap,info,dir)
|
||||
use psb_realloc_mod, only : psb_ensure_size
|
||||
implicit none
|
||||
class(psb_z_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_asort_up_
|
||||
endif
|
||||
select case(heap%dir)
|
||||
case (psb_asort_up_,psb_asort_down_)
|
||||
! ok, do nothing
|
||||
case default
|
||||
write(psb_err_unit,*) 'Invalid direction, defaulting to psb_asort_up_'
|
||||
heap%dir = psb_asort_up_
|
||||
end select
|
||||
call psb_ensure_size(psb_heap_resize,heap%keys,info)
|
||||
|
||||
return
|
||||
end subroutine psb_z_init_heap
|
||||
|
||||
|
||||
function psb_z_howmany(heap) result(res)
|
||||
implicit none
|
||||
class(psb_z_heap), intent(in) :: heap
|
||||
integer(psb_ipk_) :: res
|
||||
res = heap%last
|
||||
end function psb_z_howmany
|
||||
|
||||
subroutine psb_z_insert_heap(key,heap,info)
|
||||
use psb_realloc_mod, only : psb_ensure_size
|
||||
implicit none
|
||||
|
||||
complex(psb_dpk_), intent(in) :: key
|
||||
class(psb_z_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_z_insert_heap(key,&
|
||||
& heap%last,heap%keys,heap%dir,info)
|
||||
|
||||
return
|
||||
end subroutine psb_z_insert_heap
|
||||
|
||||
subroutine psb_z_heap_get_first(key,heap,info)
|
||||
implicit none
|
||||
|
||||
class(psb_z_heap), intent(inout) :: heap
|
||||
integer(psb_ipk_), intent(out) :: info
|
||||
complex(psb_dpk_), intent(out) :: key
|
||||
|
||||
|
||||
info = psb_success_
|
||||
|
||||
call psi_z_heap_get_first(key,&
|
||||
& heap%last,heap%keys,heap%dir,info)
|
||||
|
||||
return
|
||||
end subroutine psb_z_heap_get_first
|
||||
|
||||
subroutine psb_z_dump_heap(iout,heap,info)
|
||||
|
||||
implicit none
|
||||
class(psb_z_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_z_dump_heap
|
||||
|
||||
subroutine psb_z_free_heap(heap,info)
|
||||
implicit none
|
||||
class(psb_z_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_z_free_heap
|
||||
|
||||
subroutine psb_z_idx_init_heap(heap,info,dir)
|
||||
use psb_realloc_mod, only : psb_ensure_size
|
||||
implicit none
|
||||
class(psb_z_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_asort_up_
|
||||
endif
|
||||
select case(heap%dir)
|
||||
case (psb_asort_up_,psb_asort_down_)
|
||||
! ok, do nothing
|
||||
case default
|
||||
write(psb_err_unit,*) 'Invalid direction, defaulting to psb_asort_up_'
|
||||
heap%dir = psb_asort_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_z_idx_init_heap
|
||||
|
||||
|
||||
function psb_z_idx_howmany(heap) result(res)
|
||||
implicit none
|
||||
class(psb_z_idx_heap), intent(in) :: heap
|
||||
integer(psb_ipk_) :: res
|
||||
res = heap%last
|
||||
end function psb_z_idx_howmany
|
||||
|
||||
subroutine psb_z_idx_insert_heap(key,index,heap,info)
|
||||
use psb_realloc_mod, only : psb_ensure_size
|
||||
implicit none
|
||||
|
||||
complex(psb_dpk_), intent(in) :: key
|
||||
integer(psb_ipk_), intent(in) :: index
|
||||
class(psb_z_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_z_idx_insert_heap(key,index,&
|
||||
& heap%last,heap%keys,heap%idxs,heap%dir,info)
|
||||
|
||||
return
|
||||
end subroutine psb_z_idx_insert_heap
|
||||
|
||||
subroutine psb_z_idx_heap_get_first(key,index,heap,info)
|
||||
implicit none
|
||||
|
||||
class(psb_z_idx_heap), intent(inout) :: heap
|
||||
integer(psb_ipk_), intent(out) :: index
|
||||
integer(psb_ipk_), intent(out) :: info
|
||||
complex(psb_dpk_), intent(out) :: key
|
||||
|
||||
|
||||
info = psb_success_
|
||||
|
||||
call psi_z_idx_heap_get_first(key,index,&
|
||||
& heap%last,heap%keys,heap%idxs,heap%dir,info)
|
||||
|
||||
return
|
||||
end subroutine psb_z_idx_heap_get_first
|
||||
|
||||
subroutine psb_z_idx_dump_heap(iout,heap,info)
|
||||
|
||||
implicit none
|
||||
class(psb_z_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_z_idx_dump_heap
|
||||
|
||||
subroutine psb_z_idx_free_heap(heap,info)
|
||||
implicit none
|
||||
class(psb_z_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_z_idx_free_heap
|
||||
|
||||
end module psb_z_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
@ -0,0 +1,49 @@
|
||||
#
|
||||
# Libraries used
|
||||
#
|
||||
INSTALLDIR=../..
|
||||
INCDIR=$(INSTALLDIR)/include/
|
||||
MODDIR=$(INSTALLDIR)/modules/
|
||||
include $(INCDIR)/Make.inc.psblas
|
||||
LIBDIR=$(INSTALLDIR)/lib/
|
||||
PSBLAS_LIB= -L$(LIBDIR) -lpsb_util -lpsb_krylov -lpsb_prec -lpsb_base
|
||||
LDLIBS=$(PSBLDLIBS)
|
||||
|
||||
FINCLUDES=$(FMFLAG)$(MODDIR) $(FMFLAG).
|
||||
|
||||
|
||||
DTOBJS=d_file_spmv.o
|
||||
STOBJS=s_file_spmv.o
|
||||
DPGOBJS=pdgenspmv.o
|
||||
DVECOBJS=vecoperation.o
|
||||
EXEDIR=./runs
|
||||
|
||||
all: runsd d_file_spmv s_file_spmv pdgenspmv vecoperation
|
||||
|
||||
runsd:
|
||||
(if test ! -d runs ; then mkdir runs; fi)
|
||||
|
||||
d_file_spmv: $(DTOBJS)
|
||||
$(FLINK) $(LOPT) $(DTOBJS) -o d_file_spmv $(PSBLAS_LIB) $(LDLIBS)
|
||||
/bin/mv d_file_spmv $(EXEDIR)
|
||||
|
||||
pdgenspmv: $(DPGOBJS)
|
||||
$(FLINK) $(LOPT) $(DPGOBJS) -o pdgenspmv $(PSBLAS_LIB) $(LDLIBS)
|
||||
/bin/mv pdgenspmv $(EXEDIR)
|
||||
|
||||
|
||||
s_file_spmv: $(STOBJS)
|
||||
$(FLINK) $(LOPT) $(STOBJS) -o s_file_spmv $(PSBLAS_LIB) $(LDLIBS)
|
||||
/bin/mv s_file_spmv $(EXEDIR)
|
||||
|
||||
vecoperation: $(DVECOBJS)
|
||||
$(FLINK) $(LOPT) $(DVECOBJS) -o vecoperation $(PSBLAS_LIB) $(LDLIBS)
|
||||
/bin/mv vecoperation $(EXEDIR)
|
||||
|
||||
clean:
|
||||
/bin/rm -f $(DBOBJSS) $(DBOBJS) $(DTOBJS) $(STOBJS) $(DVECOBJS)
|
||||
|
||||
lib:
|
||||
(cd ../../; make library)
|
||||
verycleanlib:
|
||||
(cd ../../; make veryclean)
|
||||
@ -0,0 +1,297 @@
|
||||
!
|
||||
! 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.
|
||||
!
|
||||
!
|
||||
program d_file_spmv
|
||||
use psb_base_mod
|
||||
use psb_util_mod
|
||||
implicit none
|
||||
|
||||
! input parameters
|
||||
character(len=40) :: kmethd, ptype
|
||||
character(len=512) :: mtrx_file, rhs_file
|
||||
|
||||
! sparse matrices
|
||||
type(psb_dspmat_type) :: a
|
||||
type(psb_ldspmat_type) :: aux_a
|
||||
|
||||
! dense matrices
|
||||
real(psb_dpk_), allocatable, target :: aux_b(:,:), d(:)
|
||||
real(psb_dpk_), allocatable , save :: x_col_glob(:), r_col_glob(:)
|
||||
real(psb_dpk_), pointer :: b_col_glob(:)
|
||||
type(psb_d_vect_type) :: b_col, x_col, r_col
|
||||
|
||||
|
||||
! communications data structure
|
||||
type(psb_desc_type):: desc_a
|
||||
|
||||
type(psb_ctxt_type) :: ctxt
|
||||
integer(psb_ipk_) :: iam, np
|
||||
|
||||
! solver paramters
|
||||
integer(psb_ipk_) :: iter, itmax, ierr, itrace, ircode, ipart,&
|
||||
& methd, istopc, irst, nr
|
||||
integer(psb_lpk_) :: lnp
|
||||
integer(psb_epk_) :: amatsize, descsize, annz, nbytes
|
||||
real(psb_dpk_) :: err, eps,cond
|
||||
|
||||
character(len=5) :: afmt
|
||||
character(len=20) :: name
|
||||
character(len=2) :: filefmt
|
||||
integer(psb_ipk_), parameter :: iunit=12
|
||||
integer(psb_ipk_), parameter :: times=20
|
||||
integer(psb_ipk_) :: iparm(20)
|
||||
|
||||
! other variables
|
||||
integer(psb_lpk_) :: i,j,m_problem
|
||||
integer(psb_ipk_) :: internal, m,ii,nnzero, info
|
||||
real(psb_dpk_) :: t1, t2, r_amax, b_amax,&
|
||||
&scale,resmx,resmxp, flops, bdwdth
|
||||
real(psb_dpk_) :: tt1, tt2, tflops
|
||||
integer(psb_ipk_) :: nrhs, nrow, n_row, dim, nv, ne
|
||||
integer(psb_ipk_), allocatable :: ivg(:), ipv(:)
|
||||
|
||||
|
||||
call psb_init(ctxt)
|
||||
call psb_info(ctxt,iam,np)
|
||||
lnp = np
|
||||
if (iam < 0) then
|
||||
! This should not happen, but just in case
|
||||
call psb_exit(ctxt)
|
||||
stop
|
||||
endif
|
||||
|
||||
|
||||
name='d_file_spmv'
|
||||
if(psb_get_errstatus() /= 0) goto 9999
|
||||
info=psb_success_
|
||||
call psb_set_errverbosity(2)
|
||||
!
|
||||
! Hello world
|
||||
!
|
||||
if (iam == psb_root_) then
|
||||
write(*,*) 'Welcome to PSBLAS version: ',psb_version_string_
|
||||
write(*,*) 'This is the ',trim(name),' sample program'
|
||||
read(psb_inp_unit,*) mtrx_file
|
||||
read(psb_inp_unit,*) filefmt
|
||||
read(psb_inp_unit,*) ipart
|
||||
end if
|
||||
call psb_bcast(ctxt,mtrx_file)
|
||||
call psb_bcast(ctxt,filefmt)
|
||||
call psb_bcast(ctxt,ipart)
|
||||
rhs_file = 'NONE'
|
||||
afmt = 'CSR'
|
||||
call psb_barrier(ctxt)
|
||||
t1 = psb_wtime()
|
||||
! read the input matrix to be processed and (possibly) the rhs
|
||||
nrhs = 1
|
||||
|
||||
if (iam==psb_root_) then
|
||||
select case(psb_toupper(filefmt))
|
||||
case('MM')
|
||||
! For Matrix Market we have an input file for the matrix
|
||||
! and an (optional) second file for the RHS.
|
||||
call mm_mat_read(aux_a,info,iunit=iunit,filename=mtrx_file)
|
||||
if (info == psb_success_) then
|
||||
if (rhs_file /= 'NONE') then
|
||||
call mm_array_read(aux_b,info,iunit=iunit,filename=rhs_file)
|
||||
end if
|
||||
end if
|
||||
|
||||
case ('HB')
|
||||
! For Harwell-Boeing we have a single file which may or may not
|
||||
! contain an RHS.
|
||||
call hb_read(aux_a,info,iunit=iunit,b=aux_b,filename=mtrx_file)
|
||||
|
||||
case default
|
||||
info = -1
|
||||
write(psb_err_unit,*) 'Wrong choice for fileformat ', filefmt
|
||||
end select
|
||||
if (info /= psb_success_) then
|
||||
write(psb_err_unit,*) 'Error while reading input matrix '
|
||||
call psb_abort(ctxt)
|
||||
end if
|
||||
|
||||
m_problem = aux_a%get_nrows()
|
||||
call psb_bcast(ctxt,m_problem)
|
||||
|
||||
! At this point aux_b may still be unallocated
|
||||
if (psb_size(aux_b,dim=1)==m_problem) then
|
||||
! if any rhs were present, broadcast the first one
|
||||
write(psb_err_unit,'("Ok, got an rhs ")')
|
||||
b_col_glob =>aux_b(:,1)
|
||||
else
|
||||
write(psb_out_unit,'("Generating an rhs...")')
|
||||
write(psb_out_unit,'(" ")')
|
||||
call psb_realloc(m_problem,1,aux_b,ircode)
|
||||
if (ircode /= 0) then
|
||||
call psb_errpush(psb_err_alloc_dealloc_,name)
|
||||
goto 9999
|
||||
endif
|
||||
|
||||
b_col_glob => aux_b(:,1)
|
||||
do i=1, m_problem
|
||||
b_col_glob(i) = 1.d0
|
||||
enddo
|
||||
endif
|
||||
|
||||
else
|
||||
|
||||
call psb_bcast(ctxt,m_problem)
|
||||
b_col_glob =>aux_b(:,1)
|
||||
|
||||
end if
|
||||
|
||||
! switch over different partition types
|
||||
write(psb_out_unit,'("Number of processors : ",i0)')np
|
||||
if (ipart == 0) then
|
||||
call psb_barrier(ctxt)
|
||||
if (iam==psb_root_) write(psb_out_unit,'("Partition type: block")')
|
||||
allocate(ivg(m_problem),ipv(np))
|
||||
do i=1,m_problem
|
||||
call part_block(i,m_problem,np,ipv,nv)
|
||||
ivg(i) = ipv(1)
|
||||
enddo
|
||||
call psb_matdist(aux_a, a, ctxt,desc_a,info,fmt=afmt,vg=ivg)
|
||||
|
||||
else if (ipart == 2) then
|
||||
if (iam==psb_root_) then
|
||||
write(psb_out_unit,'("Partition type: graph")')
|
||||
write(psb_out_unit,'(" ")')
|
||||
! write(psb_err_unit,'("Build type: graph")')
|
||||
call build_mtpart(aux_a,lnp)
|
||||
|
||||
endif
|
||||
call psb_barrier(ctxt)
|
||||
call distr_mtpart(psb_root_,ctxt)
|
||||
call getv_mtpart(ivg)
|
||||
call psb_matdist(aux_a, a, ctxt, desc_a,info,fmt=afmt,vg=ivg)
|
||||
|
||||
else
|
||||
if (iam==psb_root_) write(psb_out_unit,'("Partition type: default block")')
|
||||
call psb_matdist(aux_a, a, ctxt, desc_a,info,fmt=afmt,parts=part_block)
|
||||
end if
|
||||
|
||||
|
||||
call psb_geall(x_col,desc_a,info)
|
||||
call x_col%set(done)
|
||||
call psb_geasb(x_col,desc_a,info)
|
||||
call psb_geall(b_col,desc_a,info)
|
||||
call x_col%zero()
|
||||
call psb_geasb(b_col,desc_a,info)
|
||||
t2 = psb_wtime() - t1
|
||||
|
||||
|
||||
call psb_amx(ctxt, t2)
|
||||
|
||||
if (iam==psb_root_) then
|
||||
write(psb_out_unit,'(" ")')
|
||||
write(psb_out_unit,'("Time to read and partition matrix : ",es12.5)')t2
|
||||
write(psb_out_unit,'(" ")')
|
||||
end if
|
||||
|
||||
|
||||
call psb_barrier(ctxt)
|
||||
t1 = psb_wtime()
|
||||
do i=1,times
|
||||
call psb_spmm(done,a,x_col,dzero,b_col,desc_a,info,'n')
|
||||
end do
|
||||
call psb_barrier(ctxt)
|
||||
t2 = psb_wtime() - t1
|
||||
call psb_amx(ctxt,t2)
|
||||
|
||||
! FIXME: cache flush needed here
|
||||
call psb_barrier(ctxt)
|
||||
tt1 = psb_wtime()
|
||||
do i=1,times
|
||||
call psb_spmm(done,a,x_col,dzero,b_col,desc_a,info,'t')
|
||||
end do
|
||||
call psb_barrier(ctxt)
|
||||
tt2 = psb_wtime() - tt1
|
||||
call psb_amx(ctxt,tt2)
|
||||
|
||||
nr = desc_a%get_global_rows()
|
||||
annz = a%get_nzeros()
|
||||
amatsize = psb_sizeof(a)
|
||||
descsize = psb_sizeof(desc_a)
|
||||
call psb_sum(ctxt,annz)
|
||||
call psb_sum(ctxt,amatsize)
|
||||
call psb_sum(ctxt,descsize)
|
||||
|
||||
if (iam==psb_root_) then
|
||||
flops = 2.d0*times*annz
|
||||
tflops=flops
|
||||
write(psb_out_unit,'("Matrix: ",a)') mtrx_file
|
||||
write(psb_out_unit,'("Test on : ",i20," processors")') np
|
||||
write(psb_out_unit,'("Size of matrix : ",i20," ")') nr
|
||||
write(psb_out_unit,'("Number of nonzeros : ",i20," ")') annz
|
||||
write(psb_out_unit,'("Memory occupation : ",i20," ")') amatsize
|
||||
write(psb_out_unit,'("Number of flops (",i0," prod) : ",F20.0," ")') times,flops
|
||||
flops = flops / (t2)
|
||||
tflops = tflops / (tt2)
|
||||
write(psb_out_unit,'("Time for ",i0," products (s) : ",F20.3)')times, t2
|
||||
write(psb_out_unit,'("Time per product (ms) : ",F20.3)') t2*1.d3/(1.d0*times)
|
||||
write(psb_out_unit,'("MFLOPS : ",F20.3)') flops/1.d6
|
||||
|
||||
write(psb_out_unit,'("Time for ",i0," products (s) (trans.): ",F20.3)') times,tt2
|
||||
write(psb_out_unit,'("Time per product (ms) (trans.): ",F20.3)') tt2*1.d3/(1.d0*times)
|
||||
write(psb_out_unit,'("MFLOPS (trans.): ",F20.3)') tflops/1.d6
|
||||
|
||||
!
|
||||
! This computation is valid for CSR
|
||||
!
|
||||
nbytes = nr*(2*psb_sizeof_dp + psb_sizeof_ip)+&
|
||||
& annz*(psb_sizeof_dp + psb_sizeof_ip)
|
||||
bdwdth = times*nbytes/(t2*1.d6)
|
||||
write(psb_out_unit,*)
|
||||
write(psb_out_unit,'("MBYTES/S : ",F20.3)') bdwdth
|
||||
bdwdth = times*nbytes/(tt2*1.d6)
|
||||
write(psb_out_unit,'("MBYTES/S (trans): ",F20.3)') bdwdth
|
||||
write(psb_out_unit,'("Storage type for DESC_A: ",a)') desc_a%get_fmt()
|
||||
|
||||
end if
|
||||
|
||||
call psb_gefree(b_col, desc_a,info)
|
||||
call psb_gefree(x_col, desc_a,info)
|
||||
call psb_spfree(a, desc_a,info)
|
||||
call psb_cdfree(desc_a,info)
|
||||
call psb_exit(ctxt)
|
||||
stop
|
||||
|
||||
9999 call psb_error(ctxt)
|
||||
|
||||
stop
|
||||
|
||||
end program d_file_spmv
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
@ -0,0 +1,770 @@
|
||||
!
|
||||
! Parallel Sparse BLAS version 3.5.1
|
||||
! (C) Copyright 2015
|
||||
! 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: ppde.f90
|
||||
!
|
||||
module psb_d_pde3d_mod
|
||||
|
||||
|
||||
use psb_base_mod, only : psb_dpk_, psb_ipk_, psb_lpk_, psb_desc_type,&
|
||||
& psb_dspmat_type, psb_d_vect_type, dzero,&
|
||||
& psb_d_base_sparse_mat, psb_d_base_vect_type, &
|
||||
& psb_i_base_vect_type, psb_l_base_vect_type
|
||||
|
||||
interface
|
||||
function d_func_3d(x,y,z) result(val)
|
||||
import :: psb_dpk_
|
||||
real(psb_dpk_), intent(in) :: x,y,z
|
||||
real(psb_dpk_) :: val
|
||||
end function d_func_3d
|
||||
end interface
|
||||
|
||||
interface psb_gen_pde3d
|
||||
module procedure psb_d_gen_pde3d
|
||||
end interface psb_gen_pde3d
|
||||
|
||||
contains
|
||||
|
||||
function d_null_func_3d(x,y,z) result(val)
|
||||
|
||||
real(psb_dpk_), intent(in) :: x,y,z
|
||||
real(psb_dpk_) :: val
|
||||
|
||||
val = dzero
|
||||
|
||||
end function d_null_func_3d
|
||||
!
|
||||
! functions parametrizing the differential equation
|
||||
!
|
||||
|
||||
!
|
||||
! Note: b1, b2 and b3 are the coefficients of the first
|
||||
! derivative of the unknown function. The default
|
||||
! we apply here is to have them zero, so that the resulting
|
||||
! matrix is symmetric/hermitian and suitable for
|
||||
! testing with CG and FCG.
|
||||
! When testing methods for non-hermitian matrices you can
|
||||
! change the B1/B2/B3 functions to e.g. done/sqrt((3*done))
|
||||
!
|
||||
function b1(x,y,z)
|
||||
use psb_base_mod, only : psb_dpk_, done, dzero
|
||||
implicit none
|
||||
real(psb_dpk_) :: b1
|
||||
real(psb_dpk_), intent(in) :: x,y,z
|
||||
b1=dzero
|
||||
end function b1
|
||||
function b2(x,y,z)
|
||||
use psb_base_mod, only : psb_dpk_, done, dzero
|
||||
implicit none
|
||||
real(psb_dpk_) :: b2
|
||||
real(psb_dpk_), intent(in) :: x,y,z
|
||||
b2=dzero
|
||||
end function b2
|
||||
function b3(x,y,z)
|
||||
use psb_base_mod, only : psb_dpk_, done, dzero
|
||||
implicit none
|
||||
real(psb_dpk_) :: b3
|
||||
real(psb_dpk_), intent(in) :: x,y,z
|
||||
b3=dzero
|
||||
end function b3
|
||||
function c(x,y,z)
|
||||
use psb_base_mod, only : psb_dpk_, done, dzero
|
||||
implicit none
|
||||
real(psb_dpk_) :: c
|
||||
real(psb_dpk_), intent(in) :: x,y,z
|
||||
c=dzero
|
||||
end function c
|
||||
function a1(x,y,z)
|
||||
use psb_base_mod, only : psb_dpk_, done, dzero
|
||||
implicit none
|
||||
real(psb_dpk_) :: a1
|
||||
real(psb_dpk_), intent(in) :: x,y,z
|
||||
a1=done/80
|
||||
end function a1
|
||||
function a2(x,y,z)
|
||||
use psb_base_mod, only : psb_dpk_, done, dzero
|
||||
implicit none
|
||||
real(psb_dpk_) :: a2
|
||||
real(psb_dpk_), intent(in) :: x,y,z
|
||||
a2=done/80
|
||||
end function a2
|
||||
function a3(x,y,z)
|
||||
use psb_base_mod, only : psb_dpk_, done, dzero
|
||||
implicit none
|
||||
real(psb_dpk_) :: a3
|
||||
real(psb_dpk_), intent(in) :: x,y,z
|
||||
a3=done/80
|
||||
end function a3
|
||||
function g(x,y,z)
|
||||
use psb_base_mod, only : psb_dpk_, done, dzero
|
||||
implicit none
|
||||
real(psb_dpk_) :: g
|
||||
real(psb_dpk_), intent(in) :: x,y,z
|
||||
g = dzero
|
||||
if (x == done) then
|
||||
g = done
|
||||
else if (x == dzero) then
|
||||
g = exp(y**2-z**2)
|
||||
end if
|
||||
end function g
|
||||
|
||||
|
||||
!
|
||||
! subroutine to allocate and fill in the coefficient matrix and
|
||||
! the rhs.
|
||||
!
|
||||
subroutine psb_d_gen_pde3d(ctxt,idim,a,bv,xv,desc_a,afmt,info,&
|
||||
& f,amold,vmold,imold,partition,nrl,iv)
|
||||
use psb_base_mod
|
||||
use psb_util_mod
|
||||
!
|
||||
! Discretizes the partial differential equation
|
||||
!
|
||||
! a1 dd(u) a2 dd(u) a3 dd(u) b1 d(u) b2 d(u) b3 d(u)
|
||||
! - ------ - ------ - ------ + ----- + ------ + ------ + c u = f
|
||||
! dxdx dydy dzdz dx dy dz
|
||||
!
|
||||
! with Dirichlet boundary conditions
|
||||
! u = g
|
||||
!
|
||||
! on the unit cube 0<=x,y,z<=1.
|
||||
!
|
||||
!
|
||||
! Note that if b1=b2=b3=c=0., the PDE is the Laplace equation.
|
||||
!
|
||||
implicit none
|
||||
integer(psb_ipk_) :: idim
|
||||
type(psb_dspmat_type) :: a
|
||||
type(psb_d_vect_type) :: xv,bv
|
||||
type(psb_desc_type) :: desc_a
|
||||
type(psb_ctxt_type) :: ctxt
|
||||
integer(psb_ipk_) :: info
|
||||
character(len=*) :: afmt
|
||||
procedure(d_func_3d), optional :: f
|
||||
class(psb_d_base_sparse_mat), optional :: amold
|
||||
class(psb_d_base_vect_type), optional :: vmold
|
||||
class(psb_i_base_vect_type), optional :: imold
|
||||
integer(psb_ipk_), optional :: partition, nrl,iv(:)
|
||||
|
||||
! Local variables.
|
||||
|
||||
integer(psb_ipk_), parameter :: nb=20
|
||||
type(psb_d_csc_sparse_mat) :: acsc
|
||||
type(psb_d_coo_sparse_mat) :: acoo
|
||||
type(psb_d_csr_sparse_mat) :: acsr
|
||||
real(psb_dpk_) :: zt(nb),x,y,z
|
||||
integer(psb_ipk_) :: nnz,nr,nlr,i,j,ii,ib,k, partition_
|
||||
integer(psb_lpk_) :: m,n,glob_row,nt
|
||||
integer(psb_ipk_) :: ix,iy,iz,ia,indx_owner
|
||||
! For 3D partition
|
||||
! Note: integer control variables going directly into an MPI call
|
||||
! must be 4 bytes, i.e. psb_mpk_
|
||||
integer(psb_mpk_) :: npdims(3), npp, minfo
|
||||
integer(psb_ipk_) :: npx,npy,npz, iamx,iamy,iamz,mynx,myny,mynz
|
||||
integer(psb_ipk_), allocatable :: bndx(:),bndy(:),bndz(:)
|
||||
! Process grid
|
||||
integer(psb_ipk_) :: np, iam
|
||||
integer(psb_ipk_) :: icoeff
|
||||
integer(psb_lpk_), allocatable :: irow(:),icol(:),myidx(:)
|
||||
real(psb_dpk_), allocatable :: val(:)
|
||||
! deltah dimension of each grid cell
|
||||
! deltat discretization time
|
||||
real(psb_dpk_) :: deltah, sqdeltah, deltah2
|
||||
real(psb_dpk_), parameter :: rhs=dzero,one=done,zero=dzero
|
||||
real(psb_dpk_) :: t0, t1, t2, t3, tasb, talc, ttot, tgen, tcdasb
|
||||
integer(psb_ipk_) :: err_act
|
||||
procedure(d_func_3d), pointer :: f_
|
||||
character(len=20) :: name, ch_err,tmpfmt
|
||||
|
||||
info = psb_success_
|
||||
name = 'create_matrix'
|
||||
call psb_erractionsave(err_act)
|
||||
|
||||
call psb_info(ctxt, iam, np)
|
||||
|
||||
|
||||
if (present(f)) then
|
||||
f_ => f
|
||||
else
|
||||
f_ => d_null_func_3d
|
||||
end if
|
||||
|
||||
deltah = done/(idim+1)
|
||||
sqdeltah = deltah*deltah
|
||||
deltah2 = (2*done)* deltah
|
||||
|
||||
if (present(partition)) then
|
||||
if ((1<= partition).and.(partition <= 3)) then
|
||||
partition_ = partition
|
||||
else
|
||||
write(*,*) 'Invalid partition choice ',partition,' defaulting to 3'
|
||||
partition_ = 3
|
||||
end if
|
||||
else
|
||||
partition_ = 3
|
||||
end if
|
||||
|
||||
! initialize array descriptor and sparse matrix storage. provide an
|
||||
! estimate of the number of non zeroes
|
||||
|
||||
m = (1_psb_lpk_*idim)*idim*idim
|
||||
n = m
|
||||
nnz = ((n*7)/(np))
|
||||
if(iam == psb_root_) write(psb_out_unit,'("Generating Matrix (size=",i0,")...")')n
|
||||
t0 = psb_wtime()
|
||||
select case(partition_)
|
||||
case(1)
|
||||
! A BLOCK partition
|
||||
if (present(nrl)) then
|
||||
nr = nrl
|
||||
else
|
||||
!
|
||||
! Using a simple BLOCK distribution.
|
||||
!
|
||||
nt = (m+np-1)/np
|
||||
nr = max(0,min(nt,m-(iam*nt)))
|
||||
end if
|
||||
|
||||
nt = nr
|
||||
call psb_sum(ctxt,nt)
|
||||
if (nt /= m) then
|
||||
write(psb_err_unit,*) iam, 'Initialization error ',nr,nt,m
|
||||
info = -1
|
||||
call psb_barrier(ctxt)
|
||||
call psb_abort(ctxt)
|
||||
return
|
||||
end if
|
||||
|
||||
!
|
||||
! First example of use of CDALL: specify for each process a number of
|
||||
! contiguous rows
|
||||
!
|
||||
call psb_cdall(ctxt,desc_a,info,nl=nr)
|
||||
myidx = desc_a%get_global_indices()
|
||||
nlr = size(myidx)
|
||||
|
||||
case(2)
|
||||
! A partition defined by the user through IV
|
||||
|
||||
if (present(iv)) then
|
||||
if (size(iv) /= m) then
|
||||
write(psb_err_unit,*) iam, 'Initialization error: wrong IV size',size(iv),m
|
||||
info = -1
|
||||
call psb_barrier(ctxt)
|
||||
call psb_abort(ctxt)
|
||||
return
|
||||
end if
|
||||
else
|
||||
write(psb_err_unit,*) iam, 'Initialization error: IV not present'
|
||||
info = -1
|
||||
call psb_barrier(ctxt)
|
||||
call psb_abort(ctxt)
|
||||
return
|
||||
end if
|
||||
|
||||
!
|
||||
! Second example of use of CDALL: specify for each row the
|
||||
! process that owns it
|
||||
!
|
||||
call psb_cdall(ctxt,desc_a,info,vg=iv)
|
||||
myidx = desc_a%get_global_indices()
|
||||
nlr = size(myidx)
|
||||
|
||||
case(3)
|
||||
! A 3-dimensional partition
|
||||
|
||||
! A nifty MPI function will split the process list
|
||||
npdims = 0
|
||||
call mpi_dims_create(np,3,npdims,info)
|
||||
npx = npdims(1)
|
||||
npy = npdims(2)
|
||||
npz = npdims(3)
|
||||
|
||||
allocate(bndx(0:npx),bndy(0:npy),bndz(0:npz))
|
||||
! We can reuse idx2ijk for process indices as well.
|
||||
call idx2ijk(iamx,iamy,iamz,iam,npx,npy,npz,base=0)
|
||||
! Now let's split the 3D cube in hexahedra
|
||||
call dist1Didx(bndx,idim,npx)
|
||||
mynx = bndx(iamx+1)-bndx(iamx)
|
||||
call dist1Didx(bndy,idim,npy)
|
||||
myny = bndy(iamy+1)-bndy(iamy)
|
||||
call dist1Didx(bndz,idim,npz)
|
||||
mynz = bndz(iamz+1)-bndz(iamz)
|
||||
|
||||
! How many indices do I own?
|
||||
nlr = mynx*myny*mynz
|
||||
allocate(myidx(nlr))
|
||||
! Now, let's generate the list of indices I own
|
||||
nr = 0
|
||||
do i=bndx(iamx),bndx(iamx+1)-1
|
||||
do j=bndy(iamy),bndy(iamy+1)-1
|
||||
do k=bndz(iamz),bndz(iamz+1)-1
|
||||
nr = nr + 1
|
||||
call ijk2idx(myidx(nr),i,j,k,idim,idim,idim)
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
if (nr /= nlr) then
|
||||
write(psb_err_unit,*) iam,iamx,iamy,iamz, 'Initialization error: NR vs NLR ',&
|
||||
& nr,nlr,mynx,myny,mynz
|
||||
info = -1
|
||||
call psb_barrier(ctxt)
|
||||
call psb_abort(ctxt)
|
||||
end if
|
||||
|
||||
!
|
||||
! Third example of use of CDALL: specify for each process
|
||||
! the set of global indices it owns.
|
||||
!
|
||||
call psb_cdall(ctxt,desc_a,info,vl=myidx)
|
||||
|
||||
case default
|
||||
write(psb_err_unit,*) iam, 'Initialization error: should not get here'
|
||||
info = -1
|
||||
call psb_barrier(ctxt)
|
||||
call psb_abort(ctxt)
|
||||
return
|
||||
end select
|
||||
|
||||
|
||||
if (info == psb_success_) call psb_spall(a,desc_a,info,nnz=nnz)
|
||||
! define rhs from boundary conditions; also build initial guess
|
||||
if (info == psb_success_) call psb_geall(xv,desc_a,info)
|
||||
if (info == psb_success_) call psb_geall(bv,desc_a,info)
|
||||
|
||||
call psb_barrier(ctxt)
|
||||
talc = psb_wtime()-t0
|
||||
|
||||
if (info /= psb_success_) then
|
||||
info=psb_err_from_subroutine_
|
||||
ch_err='allocation rout.'
|
||||
call psb_errpush(info,name,a_err=ch_err)
|
||||
goto 9999
|
||||
end if
|
||||
|
||||
! we build an auxiliary matrix consisting of one row at a
|
||||
! time; just a small matrix. might be extended to generate
|
||||
! a bunch of rows per call.
|
||||
!
|
||||
allocate(val(20*nb),irow(20*nb),&
|
||||
&icol(20*nb),stat=info)
|
||||
if (info /= psb_success_ ) then
|
||||
info=psb_err_alloc_dealloc_
|
||||
call psb_errpush(info,name)
|
||||
goto 9999
|
||||
endif
|
||||
|
||||
|
||||
! loop over rows belonging to current process in a block
|
||||
! distribution.
|
||||
|
||||
call psb_barrier(ctxt)
|
||||
t1 = psb_wtime()
|
||||
do ii=1, nlr,nb
|
||||
ib = min(nb,nlr-ii+1)
|
||||
icoeff = 1
|
||||
do k=1,ib
|
||||
i=ii+k-1
|
||||
! local matrix pointer
|
||||
glob_row=myidx(i)
|
||||
! compute gridpoint coordinates
|
||||
call idx2ijk(ix,iy,iz,glob_row,idim,idim,idim)
|
||||
! x, y, z coordinates
|
||||
x = (ix-1)*deltah
|
||||
y = (iy-1)*deltah
|
||||
z = (iz-1)*deltah
|
||||
zt(k) = f_(x,y,z)
|
||||
! internal point: build discretization
|
||||
!
|
||||
! term depending on (x-1,y,z)
|
||||
!
|
||||
val(icoeff) = -a1(x,y,z)/sqdeltah-b1(x,y,z)/deltah2
|
||||
if (ix == 1) then
|
||||
zt(k) = g(dzero,y,z)*(-val(icoeff)) + zt(k)
|
||||
else
|
||||
call ijk2idx(icol(icoeff),ix-1,iy,iz,idim,idim,idim)
|
||||
irow(icoeff) = glob_row
|
||||
icoeff = icoeff+1
|
||||
endif
|
||||
! term depending on (x,y-1,z)
|
||||
val(icoeff) = -a2(x,y,z)/sqdeltah-b2(x,y,z)/deltah2
|
||||
if (iy == 1) then
|
||||
zt(k) = g(x,dzero,z)*(-val(icoeff)) + zt(k)
|
||||
else
|
||||
call ijk2idx(icol(icoeff),ix,iy-1,iz,idim,idim,idim)
|
||||
irow(icoeff) = glob_row
|
||||
icoeff = icoeff+1
|
||||
endif
|
||||
! term depending on (x,y,z-1)
|
||||
val(icoeff)=-a3(x,y,z)/sqdeltah-b3(x,y,z)/deltah2
|
||||
if (iz == 1) then
|
||||
zt(k) = g(x,y,dzero)*(-val(icoeff)) + zt(k)
|
||||
else
|
||||
call ijk2idx(icol(icoeff),ix,iy,iz-1,idim,idim,idim)
|
||||
irow(icoeff) = glob_row
|
||||
icoeff = icoeff+1
|
||||
endif
|
||||
|
||||
! term depending on (x,y,z)
|
||||
val(icoeff)=(2*done)*(a1(x,y,z)+a2(x,y,z)+a3(x,y,z))/sqdeltah &
|
||||
& + c(x,y,z)
|
||||
call ijk2idx(icol(icoeff),ix,iy,iz,idim,idim,idim)
|
||||
irow(icoeff) = glob_row
|
||||
icoeff = icoeff+1
|
||||
! term depending on (x,y,z+1)
|
||||
val(icoeff)=-a3(x,y,z)/sqdeltah+b3(x,y,z)/deltah2
|
||||
if (iz == idim) then
|
||||
zt(k) = g(x,y,done)*(-val(icoeff)) + zt(k)
|
||||
else
|
||||
call ijk2idx(icol(icoeff),ix,iy,iz+1,idim,idim,idim)
|
||||
irow(icoeff) = glob_row
|
||||
icoeff = icoeff+1
|
||||
endif
|
||||
! term depending on (x,y+1,z)
|
||||
val(icoeff)=-a2(x,y,z)/sqdeltah+b2(x,y,z)/deltah2
|
||||
if (iy == idim) then
|
||||
zt(k) = g(x,done,z)*(-val(icoeff)) + zt(k)
|
||||
else
|
||||
call ijk2idx(icol(icoeff),ix,iy+1,iz,idim,idim,idim)
|
||||
irow(icoeff) = glob_row
|
||||
icoeff = icoeff+1
|
||||
endif
|
||||
! term depending on (x+1,y,z)
|
||||
val(icoeff)=-a1(x,y,z)/sqdeltah+b1(x,y,z)/deltah2
|
||||
if (ix==idim) then
|
||||
zt(k) = g(done,y,z)*(-val(icoeff)) + zt(k)
|
||||
else
|
||||
call ijk2idx(icol(icoeff),ix+1,iy,iz,idim,idim,idim)
|
||||
irow(icoeff) = glob_row
|
||||
icoeff = icoeff+1
|
||||
endif
|
||||
|
||||
end do
|
||||
call psb_spins(icoeff-1,irow,icol,val,a,desc_a,info)
|
||||
if(info /= psb_success_) exit
|
||||
call psb_geins(ib,myidx(ii:ii+ib-1),zt(1:ib),bv,desc_a,info)
|
||||
if(info /= psb_success_) exit
|
||||
zt(:)=dzero
|
||||
call psb_geins(ib,myidx(ii:ii+ib-1),zt(1:ib),xv,desc_a,info)
|
||||
if(info /= psb_success_) exit
|
||||
end do
|
||||
|
||||
tgen = psb_wtime()-t1
|
||||
if(info /= psb_success_) then
|
||||
info=psb_err_from_subroutine_
|
||||
ch_err='insert rout.'
|
||||
call psb_errpush(info,name,a_err=ch_err)
|
||||
goto 9999
|
||||
end if
|
||||
|
||||
deallocate(val,irow,icol)
|
||||
|
||||
call psb_barrier(ctxt)
|
||||
t1 = psb_wtime()
|
||||
call psb_cdasb(desc_a,info,mold=imold)
|
||||
tcdasb = psb_wtime()-t1
|
||||
call psb_barrier(ctxt)
|
||||
t1 = psb_wtime()
|
||||
if (info == psb_success_) then
|
||||
if (present(amold)) then
|
||||
call psb_spasb(a,desc_a,info,dupl=psb_dupl_err_,mold=amold)
|
||||
else
|
||||
call psb_spasb(a,desc_a,info,dupl=psb_dupl_err_,afmt=afmt)
|
||||
end if
|
||||
end if
|
||||
call psb_barrier(ctxt)
|
||||
if(info /= psb_success_) then
|
||||
info=psb_err_from_subroutine_
|
||||
ch_err='asb rout.'
|
||||
call psb_errpush(info,name,a_err=ch_err)
|
||||
goto 9999
|
||||
end if
|
||||
if (info == psb_success_) call psb_geasb(xv,desc_a,info,mold=vmold)
|
||||
if (info == psb_success_) call psb_geasb(bv,desc_a,info,mold=vmold)
|
||||
if(info /= psb_success_) then
|
||||
info=psb_err_from_subroutine_
|
||||
ch_err='asb rout.'
|
||||
call psb_errpush(info,name,a_err=ch_err)
|
||||
goto 9999
|
||||
end if
|
||||
tasb = psb_wtime()-t1
|
||||
call psb_barrier(ctxt)
|
||||
ttot = psb_wtime() - t0
|
||||
|
||||
call psb_amx(ctxt,talc)
|
||||
call psb_amx(ctxt,tgen)
|
||||
call psb_amx(ctxt,tasb)
|
||||
call psb_amx(ctxt,ttot)
|
||||
if(iam == psb_root_) then
|
||||
tmpfmt = a%get_fmt()
|
||||
write(psb_out_unit,'("The matrix has been generated and assembled in ",a3," format.")')&
|
||||
& tmpfmt
|
||||
write(psb_out_unit,'("-allocation time : ",es12.5)') talc
|
||||
write(psb_out_unit,'("-coeff. gen. time : ",es12.5)') tgen
|
||||
write(psb_out_unit,'("-desc asbly time : ",es12.5)') tcdasb
|
||||
write(psb_out_unit,'("- mat asbly time : ",es12.5)') tasb
|
||||
write(psb_out_unit,'("-total time : ",es12.5)') ttot
|
||||
|
||||
end if
|
||||
call psb_erractionrestore(err_act)
|
||||
return
|
||||
|
||||
9999 call psb_error_handler(ctxt,err_act)
|
||||
|
||||
return
|
||||
end subroutine psb_d_gen_pde3d
|
||||
|
||||
|
||||
end module psb_d_pde3d_mod
|
||||
|
||||
program pdgenspmv
|
||||
use psb_base_mod
|
||||
use psb_util_mod
|
||||
use psb_d_pde3d_mod
|
||||
implicit none
|
||||
|
||||
! input parameters
|
||||
character(len=20) :: kmethd, ptype
|
||||
character(len=5) :: afmt
|
||||
integer(psb_ipk_) :: idim
|
||||
|
||||
! miscellaneous
|
||||
real(psb_dpk_), parameter :: one = done
|
||||
real(psb_dpk_) :: t1, t2, tprec, flops, tflops, tt1, tt2, bdwdth
|
||||
|
||||
! sparse matrix and preconditioner
|
||||
type(psb_dspmat_type) :: a
|
||||
! descriptor
|
||||
type(psb_desc_type) :: desc_a
|
||||
! dense matrices
|
||||
type(psb_d_vect_type) :: xv,bv, vtst
|
||||
real(psb_dpk_), allocatable :: tst(:)
|
||||
! blacs parameters
|
||||
type(psb_ctxt_type) :: ctxt
|
||||
integer(psb_ipk_) :: iam, np
|
||||
|
||||
! solver parameters
|
||||
integer(psb_ipk_) :: iter, itmax,itrace, istopc, irst, nr, ipart
|
||||
integer(psb_epk_) :: amatsize, precsize, descsize, d2size, annz, nbytes
|
||||
real(psb_dpk_) :: err, eps
|
||||
integer(psb_ipk_), parameter :: times=10
|
||||
|
||||
! other variables
|
||||
integer(psb_ipk_) :: info, i
|
||||
character(len=20) :: name,ch_err
|
||||
character(len=40) :: fname
|
||||
|
||||
info=psb_success_
|
||||
|
||||
|
||||
call psb_init(ctxt)
|
||||
call psb_info(ctxt,iam,np)
|
||||
|
||||
if (iam < 0) then
|
||||
! This should not happen, but just in case
|
||||
call psb_exit(ctxt)
|
||||
stop
|
||||
endif
|
||||
if(psb_get_errstatus() /= 0) goto 9999
|
||||
name='pde90'
|
||||
call psb_set_errverbosity(itwo)
|
||||
!
|
||||
! Hello world
|
||||
!
|
||||
if (iam == psb_root_) then
|
||||
write(*,*) 'Welcome to PSBLAS version: ',psb_version_string_
|
||||
write(*,*) 'This is the ',trim(name),' sample program'
|
||||
end if
|
||||
!
|
||||
! get parameters
|
||||
!
|
||||
call get_parms(ctxt,afmt,idim)
|
||||
|
||||
!
|
||||
! allocate and fill in the coefficient matrix, rhs and initial guess
|
||||
!
|
||||
call psb_barrier(ctxt)
|
||||
t1 = psb_wtime()
|
||||
call psb_gen_pde3d(ctxt,idim,a,bv,xv,desc_a,afmt,info)
|
||||
call psb_barrier(ctxt)
|
||||
t2 = psb_wtime() - t1
|
||||
if(info /= psb_success_) then
|
||||
info=psb_err_from_subroutine_
|
||||
ch_err='psb_gen_pde3d'
|
||||
call psb_errpush(info,name,a_err=ch_err)
|
||||
goto 9999
|
||||
end if
|
||||
if (iam == psb_root_) write(psb_out_unit,'("Overall matrix creation time : ",es12.5)')t2
|
||||
if (iam == psb_root_) write(psb_out_unit,'(" ")')
|
||||
|
||||
call xv%set(done)
|
||||
|
||||
call psb_barrier(ctxt)
|
||||
t1 = psb_wtime()
|
||||
!
|
||||
! Perform Ax multiple times to compute average performance
|
||||
!
|
||||
do i=1,times
|
||||
call psb_spmm(done,a,xv,dzero,bv,desc_a,info,'n')
|
||||
end do
|
||||
call psb_barrier(ctxt)
|
||||
t2 = psb_wtime() - t1
|
||||
call psb_amx(ctxt,t2)
|
||||
|
||||
! FIXME: cache flush needed here
|
||||
call psb_barrier(ctxt)
|
||||
tt1 = psb_wtime()
|
||||
!
|
||||
! Perform A^Tx multiple times to compute average performance
|
||||
!
|
||||
do i=1,times
|
||||
call psb_spmm(done,a,xv,dzero,bv,desc_a,info,'t')
|
||||
end do
|
||||
call psb_barrier(ctxt)
|
||||
tt2 = psb_wtime() - tt1
|
||||
call psb_amx(ctxt,tt2)
|
||||
|
||||
call psb_amx(ctxt,t2)
|
||||
nr = desc_a%get_global_rows()
|
||||
annz = a%get_nzeros()
|
||||
amatsize = a%sizeof()
|
||||
descsize = psb_sizeof(desc_a)
|
||||
call psb_sum(ctxt,annz)
|
||||
call psb_sum(ctxt,amatsize)
|
||||
call psb_sum(ctxt,descsize)
|
||||
|
||||
if (iam == psb_root_) then
|
||||
flops = 2.d0*times*annz
|
||||
tflops=flops
|
||||
write(psb_out_unit,'("Matrix: ell1 ",i0)') idim
|
||||
write(psb_out_unit,'("Test on : ",i20," processors")') np
|
||||
write(psb_out_unit,'("Size of matrix : ",i20," ")') nr
|
||||
write(psb_out_unit,'("Number of nonzeros : ",i20," ")') annz
|
||||
write(psb_out_unit,'("Memory occupation : ",i20," ")') amatsize
|
||||
write(psb_out_unit,'("Number of flops (",i0," prod) : ",F20.0," ")') times,flops
|
||||
flops = flops / (t2)
|
||||
tflops = tflops / (tt2)
|
||||
write(psb_out_unit,'("Time for ",i0," products (s) : ",F20.3)')times, t2
|
||||
write(psb_out_unit,'("Time per product (ms) : ",F20.3)') t2*1.d3/(1.d0*times)
|
||||
write(psb_out_unit,'("MFLOPS : ",F20.3)') flops/1.d6
|
||||
|
||||
write(psb_out_unit,'("Time for ",i0," products (s) (trans.): ",F20.3)') times,tt2
|
||||
write(psb_out_unit,'("Time per product (ms) (trans.): ",F20.3)') tt2*1.d3/(1.d0*times)
|
||||
write(psb_out_unit,'("MFLOPS (trans.): ",F20.3)') tflops/1.d6
|
||||
|
||||
!
|
||||
! This computation is valid for CSR
|
||||
!
|
||||
nbytes = nr*(2*psb_sizeof_dp + psb_sizeof_ip)+&
|
||||
& annz*(psb_sizeof_dp + psb_sizeof_ip)
|
||||
bdwdth = times*nbytes/(t2*1.d6)
|
||||
write(psb_out_unit,*)
|
||||
write(psb_out_unit,'("MBYTES/S : ",F20.3)') bdwdth
|
||||
bdwdth = times*nbytes/(tt2*1.d6)
|
||||
write(psb_out_unit,'("MBYTES/S (trans): ",F20.3)') bdwdth
|
||||
write(psb_out_unit,'("Storage type for DESC_A: ",a)') desc_a%get_fmt()
|
||||
write(psb_out_unit,'("Total memory occupation for DESC_A: ",i12)')descsize
|
||||
|
||||
end if
|
||||
|
||||
|
||||
!
|
||||
! cleanup storage and exit
|
||||
!
|
||||
call psb_gefree(bv,desc_a,info)
|
||||
call psb_gefree(xv,desc_a,info)
|
||||
call psb_spfree(a,desc_a,info)
|
||||
call psb_cdfree(desc_a,info)
|
||||
if(info /= psb_success_) then
|
||||
info=psb_err_from_subroutine_
|
||||
ch_err='free routine'
|
||||
call psb_errpush(info,name,a_err=ch_err)
|
||||
goto 9999
|
||||
end if
|
||||
|
||||
call psb_exit(ctxt)
|
||||
stop
|
||||
|
||||
9999 call psb_error(ctxt)
|
||||
|
||||
stop
|
||||
|
||||
contains
|
||||
!
|
||||
! get iteration parameters from standard input
|
||||
!
|
||||
subroutine get_parms(ctxt,afmt,idim)
|
||||
type(psb_ctxt_type) :: ctxt
|
||||
character(len=*) :: afmt
|
||||
integer(psb_ipk_) :: idim
|
||||
integer(psb_ipk_) :: np, iam
|
||||
integer(psb_ipk_) :: intbuf(10), ip
|
||||
|
||||
call psb_info(ctxt, iam, np)
|
||||
|
||||
if (iam == 0) then
|
||||
read(psb_inp_unit,*) afmt
|
||||
read(psb_inp_unit,*) idim
|
||||
endif
|
||||
call psb_bcast(ctxt,afmt)
|
||||
call psb_bcast(ctxt,idim)
|
||||
|
||||
if (iam == 0) then
|
||||
write(psb_out_unit,'("Testing matrix : ell1")')
|
||||
write(psb_out_unit,'("Grid dimensions : ",i4,"x",i4,"x",i4)')idim,idim,idim
|
||||
write(psb_out_unit,'("Number of processors : ",i0)')np
|
||||
write(psb_out_unit,'("Data distribution : BLOCK")')
|
||||
write(psb_out_unit,'(" ")')
|
||||
end if
|
||||
return
|
||||
|
||||
end subroutine get_parms
|
||||
!
|
||||
! print an error message
|
||||
!
|
||||
subroutine pr_usage(iout)
|
||||
integer(psb_ipk_) :: iout
|
||||
write(iout,*)'incorrect parameter(s) found'
|
||||
write(iout,*)' usage: pde90 methd prec dim &
|
||||
&[istop itmax itrace]'
|
||||
write(iout,*)' where:'
|
||||
write(iout,*)' methd: cgstab cgs rgmres bicgstabl'
|
||||
write(iout,*)' prec : bjac diag none'
|
||||
write(iout,*)' dim number of points along each axis'
|
||||
write(iout,*)' the size of the resulting linear '
|
||||
write(iout,*)' system is dim**3'
|
||||
write(iout,*)' istop stopping criterion 1, 2 '
|
||||
write(iout,*)' itmax maximum number of iterations [500] '
|
||||
write(iout,*)' itrace <=0 (no tracing, default) or '
|
||||
write(iout,*)' >= 1 do tracing every itrace'
|
||||
write(iout,*)' iterations '
|
||||
end subroutine pr_usage
|
||||
end program pdgenspmv
|
||||
@ -0,0 +1,5 @@
|
||||
pde100.mtx
|
||||
MM
|
||||
0
|
||||
|
||||
|
||||
@ -0,0 +1,3 @@
|
||||
CSR
|
||||
50
|
||||
|
||||
@ -0,0 +1,295 @@
|
||||
!
|
||||
! 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.
|
||||
!
|
||||
!
|
||||
program s_file_spmv
|
||||
use psb_base_mod
|
||||
use psb_util_mod
|
||||
implicit none
|
||||
|
||||
! input parameters
|
||||
character(len=40) :: kmethd, ptype, mtrx_file, rhs_file
|
||||
|
||||
! sparse matrices
|
||||
type(psb_sspmat_type) :: a
|
||||
type(psb_lsspmat_type) :: aux_a
|
||||
|
||||
! dense matrices
|
||||
real(psb_spk_), allocatable, target :: aux_b(:,:), d(:)
|
||||
real(psb_spk_), allocatable , save :: x_col_glob(:), r_col_glob(:)
|
||||
real(psb_spk_), pointer :: b_col_glob(:)
|
||||
type(psb_s_vect_type) :: b_col, x_col, r_col
|
||||
|
||||
|
||||
! communications data structure
|
||||
type(psb_desc_type):: desc_a
|
||||
|
||||
type(psb_ctxt_type) :: ctxt
|
||||
integer(psb_ipk_) :: iam, np
|
||||
|
||||
! solver paramters
|
||||
integer(psb_ipk_) :: iter, itmax, ierr, itrace, ircode, ipart,&
|
||||
& methd, istopc, irst, nr
|
||||
integer(psb_lpk_) :: lnp
|
||||
integer(psb_epk_) :: amatsize, descsize, annz, nbytes
|
||||
real(psb_spk_) :: err, eps,cond
|
||||
|
||||
character(len=5) :: afmt
|
||||
character(len=20) :: name
|
||||
character(len=2) :: filefmt
|
||||
integer(psb_ipk_), parameter :: iunit=12
|
||||
integer(psb_ipk_), parameter :: times=20
|
||||
integer(psb_ipk_) :: iparm(20)
|
||||
|
||||
! other variables
|
||||
integer(psb_lpk_) :: i,j,m_problem
|
||||
integer(psb_ipk_) :: internal, m,ii,nnzero,info
|
||||
real(psb_dpk_) :: t1, t2, r_amax, b_amax,&
|
||||
&scale,resmx,resmxp, flops, bdwdth
|
||||
real(psb_dpk_) :: tt1, tt2, tflops
|
||||
integer(psb_ipk_) :: nrhs, nrow, n_row, dim, nv, ne
|
||||
integer(psb_ipk_), allocatable :: ivg(:), ipv(:)
|
||||
|
||||
|
||||
call psb_init(ctxt)
|
||||
call psb_info(ctxt,iam,np)
|
||||
|
||||
if (iam < 0) then
|
||||
! This should not happen, but just in case
|
||||
call psb_exit(ctxt)
|
||||
stop
|
||||
endif
|
||||
|
||||
|
||||
name='s_file_spmv'
|
||||
if(psb_get_errstatus() /= 0) goto 9999
|
||||
info=psb_success_
|
||||
call psb_set_errverbosity(2)
|
||||
!
|
||||
! Hello world
|
||||
!
|
||||
if (iam == psb_root_) then
|
||||
write(*,*) 'Welcome to PSBLAS version: ',psb_version_string_
|
||||
write(*,*) 'This is the ',trim(name),' sample program'
|
||||
read(psb_inp_unit,*) mtrx_file
|
||||
read(psb_inp_unit,*) filefmt
|
||||
read(psb_inp_unit,*) ipart
|
||||
end if
|
||||
call psb_bcast(ctxt,mtrx_file)
|
||||
call psb_bcast(ctxt,filefmt)
|
||||
call psb_bcast(ctxt,ipart)
|
||||
rhs_file = 'NONE'
|
||||
afmt = 'CSR'
|
||||
call psb_barrier(ctxt)
|
||||
t1 = psb_wtime()
|
||||
! read the input matrix to be processed and (possibly) the rhs
|
||||
nrhs = 1
|
||||
|
||||
if (iam==psb_root_) then
|
||||
select case(psb_toupper(filefmt))
|
||||
case('MM')
|
||||
! For Matrix Market we have an input file for the matrix
|
||||
! and an (optional) second file for the RHS.
|
||||
call mm_mat_read(aux_a,info,iunit=iunit,filename=mtrx_file)
|
||||
if (info == psb_success_) then
|
||||
if (rhs_file /= 'NONE') then
|
||||
call mm_array_read(aux_b,info,iunit=iunit,filename=rhs_file)
|
||||
end if
|
||||
end if
|
||||
|
||||
case ('HB')
|
||||
! For Harwell-Boeing we have a single file which may or may not
|
||||
! contain an RHS.
|
||||
call hb_read(aux_a,info,iunit=iunit,b=aux_b,filename=mtrx_file)
|
||||
|
||||
case default
|
||||
info = -1
|
||||
write(psb_err_unit,*) 'Wrong choice for fileformat ', filefmt
|
||||
end select
|
||||
if (info /= psb_success_) then
|
||||
write(psb_err_unit,*) 'Error while reading input matrix '
|
||||
call psb_abort(ctxt)
|
||||
end if
|
||||
|
||||
m_problem = aux_a%get_nrows()
|
||||
call psb_bcast(ctxt,m_problem)
|
||||
|
||||
! At this point aux_b may still be unallocated
|
||||
if (psb_size(aux_b,dim=1)==m_problem) then
|
||||
! if any rhs were present, broadcast the first one
|
||||
write(psb_err_unit,'("Ok, got an rhs ")')
|
||||
b_col_glob =>aux_b(:,1)
|
||||
else
|
||||
write(psb_out_unit,'("Generating an rhs...")')
|
||||
write(psb_out_unit,'(" ")')
|
||||
call psb_realloc(m_problem,1,aux_b,ircode)
|
||||
if (ircode /= 0) then
|
||||
call psb_errpush(psb_err_alloc_dealloc_,name)
|
||||
goto 9999
|
||||
endif
|
||||
|
||||
b_col_glob => aux_b(:,1)
|
||||
do i=1, m_problem
|
||||
b_col_glob(i) = 1.d0
|
||||
enddo
|
||||
endif
|
||||
|
||||
else
|
||||
|
||||
call psb_bcast(ctxt,m_problem)
|
||||
b_col_glob =>aux_b(:,1)
|
||||
|
||||
end if
|
||||
|
||||
! switch over different partition types
|
||||
write(psb_out_unit,'("Number of processors : ",i0)')np
|
||||
if (ipart == 0) then
|
||||
call psb_barrier(ctxt)
|
||||
if (iam==psb_root_) write(psb_out_unit,'("Partition type: block")')
|
||||
allocate(ivg(m_problem),ipv(np))
|
||||
do i=1,m_problem
|
||||
call part_block(i,m_problem,np,ipv,nv)
|
||||
ivg(i) = ipv(1)
|
||||
enddo
|
||||
call psb_matdist(aux_a, a, ctxt,desc_a,info,fmt=afmt,vg=ivg)
|
||||
|
||||
else if (ipart == 2) then
|
||||
if (iam==psb_root_) then
|
||||
write(psb_out_unit,'("Partition type: graph")')
|
||||
write(psb_out_unit,'(" ")')
|
||||
! write(psb_err_unit,'("Build type: graph")')
|
||||
call build_mtpart(aux_a,lnp)
|
||||
|
||||
endif
|
||||
call psb_barrier(ctxt)
|
||||
call distr_mtpart(psb_root_,ctxt)
|
||||
call getv_mtpart(ivg)
|
||||
call psb_matdist(aux_a, a, ctxt, desc_a,info,fmt=afmt,vg=ivg)
|
||||
|
||||
else
|
||||
if (iam==psb_root_) write(psb_out_unit,'("Partition type: default block")')
|
||||
call psb_matdist(aux_a, a, ctxt, desc_a,info,fmt=afmt,parts=part_block)
|
||||
end if
|
||||
|
||||
call psb_geall(x_col,desc_a,info)
|
||||
call x_col%set(sone)
|
||||
call psb_geasb(x_col,desc_a,info)
|
||||
call psb_geall(b_col,desc_a,info)
|
||||
call x_col%zero()
|
||||
call psb_geasb(b_col,desc_a,info)
|
||||
t2 = psb_wtime() - t1
|
||||
|
||||
|
||||
call psb_amx(ctxt, t2)
|
||||
|
||||
if (iam==psb_root_) then
|
||||
write(psb_out_unit,'(" ")')
|
||||
write(psb_out_unit,'("Time to read and partition matrix : ",es12.5)')t2
|
||||
write(psb_out_unit,'(" ")')
|
||||
end if
|
||||
|
||||
|
||||
call psb_barrier(ctxt)
|
||||
t1 = psb_wtime()
|
||||
do i=1,times
|
||||
call psb_spmm(sone,a,x_col,szero,b_col,desc_a,info,'n')
|
||||
end do
|
||||
call psb_barrier(ctxt)
|
||||
t2 = psb_wtime() - t1
|
||||
call psb_amx(ctxt,t2)
|
||||
|
||||
! FIXME: cache flush needed here
|
||||
call psb_barrier(ctxt)
|
||||
tt1 = psb_wtime()
|
||||
do i=1,times
|
||||
call psb_spmm(sone,a,x_col,szero,b_col,desc_a,info,'t')
|
||||
end do
|
||||
call psb_barrier(ctxt)
|
||||
tt2 = psb_wtime() - tt1
|
||||
call psb_amx(ctxt,tt2)
|
||||
|
||||
nr = desc_a%get_global_rows()
|
||||
annz = a%get_nzeros()
|
||||
amatsize = psb_sizeof(a)
|
||||
descsize = psb_sizeof(desc_a)
|
||||
call psb_sum(ctxt,annz)
|
||||
call psb_sum(ctxt,amatsize)
|
||||
call psb_sum(ctxt,descsize)
|
||||
|
||||
if (iam==psb_root_) then
|
||||
flops = 2.d0*times*annz
|
||||
tflops=flops
|
||||
write(psb_out_unit,'("Matrix: ",a)') mtrx_file
|
||||
write(psb_out_unit,'("Test on : ",i20," processors")') np
|
||||
write(psb_out_unit,'("Size of matrix : ",i20," ")') nr
|
||||
write(psb_out_unit,'("Number of nonzeros : ",i20," ")') annz
|
||||
write(psb_out_unit,'("Memory occupation : ",i20," ")') amatsize
|
||||
write(psb_out_unit,'("Number of flops (",i0," prod) : ",F20.0," ")') times,flops
|
||||
flops = flops / (t2)
|
||||
tflops = tflops / (tt2)
|
||||
write(psb_out_unit,'("Time for ",i0," products (s) : ",F20.3)')times, t2
|
||||
write(psb_out_unit,'("Time per product (ms) : ",F20.3)') t2*1.d3/(1.d0*times)
|
||||
write(psb_out_unit,'("MFLOPS : ",F20.3)') flops/1.d6
|
||||
|
||||
write(psb_out_unit,'("Time for ",i0," products (s) (trans.): ",F20.3)') times,tt2
|
||||
write(psb_out_unit,'("Time per product (ms) (trans.): ",F20.3)') tt2*1.d3/(1.d0*times)
|
||||
write(psb_out_unit,'("MFLOPS (trans.): ",F20.3)') tflops/1.d6
|
||||
|
||||
!
|
||||
! This computation is valid for CSR
|
||||
!
|
||||
nbytes = nr*(2*psb_sizeof_sp + psb_sizeof_ip)+ &
|
||||
& annz*(psb_sizeof_sp + psb_sizeof_ip)
|
||||
bdwdth = times*nbytes/(t2*1.d6)
|
||||
write(psb_out_unit,*)
|
||||
write(psb_out_unit,'("MBYTES/S : ",F20.3)') bdwdth
|
||||
bdwdth = times*nbytes/(tt2*1.d6)
|
||||
write(psb_out_unit,'("MBYTES/S (trans): ",F20.3)') bdwdth
|
||||
|
||||
end if
|
||||
|
||||
call psb_gefree(b_col, desc_a,info)
|
||||
call psb_gefree(x_col, desc_a,info)
|
||||
call psb_spfree(a, desc_a,info)
|
||||
call psb_cdfree(desc_a,info)
|
||||
|
||||
call psb_exit(ctxt)
|
||||
stop
|
||||
|
||||
9999 call psb_error(ctxt)
|
||||
|
||||
stop
|
||||
|
||||
end program s_file_spmv
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
Some files were not shown because too many files have changed in this diff Show More
Loading…
Reference in New Issue