New files for 2-bytes integer support. To be compiled and tested.

This commit is contained in:
Salvatore Filippone
2020-04-02 10:27:49 +02:00
parent db22e832ee
commit 4c7fa23c8a
30 changed files with 9753 additions and 4 deletions
+129
View File
@@ -0,0 +1,129 @@
!
! 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.
!
!
! Subroutine: psi_i2ovrl_restr
! These subroutines restore the overlap region of a vector; they are used
! for the transpose matrix-vector product when there is a nonempty overlap.
!
!
subroutine psi_i2ovrl_restrr1(x,xs,desc_a,info)
use psi_mod, psi_protect_name => psi_i2ovrl_restrr1
implicit none
integer(psb_i2pk_), intent(inout) :: x(:)
integer(psb_i2pk_) :: xs(:)
type(psb_desc_type), intent(in) :: desc_a
integer(psb_ipk_), intent(out) :: info
! locals
integer(psb_ipk_) :: ictxt, np, me, err_act, i, idx, isz
character(len=20) :: name, ch_err
name='psi_i2ovrl_restrr1'
info = psb_success_
call psb_erractionsave(err_act)
if (psb_errstatus_fatal()) then
info = psb_err_internal_error_ ; goto 9999
end if
ictxt = desc_a%get_context()
call psb_info(ictxt, me, np)
if (np == -1) then
info = psb_err_context_error_
call psb_errpush(info,name)
goto 9999
endif
isz = size(desc_a%ovrlap_elem,1)
do i=1, isz
idx = desc_a%ovrlap_elem(i,1)
x(idx) = xs(i)
end do
call psb_erractionrestore(err_act)
return
9999 call psb_error_handler(ictxt,err_act)
return
end subroutine psi_i2ovrl_restrr1
subroutine psi_i2ovrl_restrr2(x,xs,desc_a,info)
use psi_mod, psi_protect_name => psi_i2ovrl_restrr2
implicit none
integer(psb_i2pk_), intent(inout) :: x(:,:)
integer(psb_i2pk_) :: xs(:,:)
type(psb_desc_type), intent(in) :: desc_a
integer(psb_ipk_), intent(out) :: info
! locals
integer(psb_ipk_) :: ictxt, np, me, err_act, i, idx, isz
character(len=20) :: name, ch_err
name='psi_i2ovrl_restrr2'
info = psb_success_
call psb_erractionsave(err_act)
if (psb_errstatus_fatal()) then
info = psb_err_internal_error_ ; goto 9999
end if
ictxt = desc_a%get_context()
call psb_info(ictxt, me, np)
if (np == -1) then
info = psb_err_context_error_
call psb_errpush(info,name)
goto 9999
endif
if (size(x,2) /= size(xs,2)) then
info = psb_err_internal_error_
call psb_errpush(info,name, a_err='Mismacth columns X vs XS')
goto 9999
endif
isz = size(desc_a%ovrlap_elem,1)
do i=1, isz
idx = desc_a%ovrlap_elem(i,1)
x(idx,:) = xs(i,:)
end do
call psb_erractionrestore(err_act)
return
9999 call psb_error_handler(ictxt,err_act)
return
end subroutine psi_i2ovrl_restrr2
+139
View File
@@ -0,0 +1,139 @@
!
! 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.
!
!
!
! Subroutine: psi_i2ovrl_save
! These subroutines save the overlap region of a vector; they are used
! for the transpose matrix-vector product when there is a nonempty overlap.
!
subroutine psi_i2ovrl_saver1(x,xs,desc_a,info)
use psi_mod, psi_protect_name => psi_i2ovrl_saver1
use psb_realloc_mod
implicit none
integer(psb_i2pk_), intent(inout) :: x(:)
integer(psb_i2pk_), allocatable :: xs(:)
type(psb_desc_type), intent(in) :: desc_a
integer(psb_ipk_), intent(out) :: info
! locals
integer(psb_ipk_) :: ictxt, np, me, err_act, i, idx, isz
character(len=20) :: name, ch_err
name='psi_i2ovrl_saver1'
info = psb_success_
call psb_erractionsave(err_act)
if (psb_errstatus_fatal()) then
info = psb_err_internal_error_ ; goto 9999
end if
ictxt = desc_a%get_context()
call psb_info(ictxt, me, np)
if (np == -1) then
info = psb_err_context_error_
call psb_errpush(info,name)
goto 9999
endif
isz = size(desc_a%ovrlap_elem,1)
call psb_realloc(isz,xs,info)
if (info /= psb_success_) then
info = psb_err_alloc_dealloc_
call psb_errpush(info,name)
goto 9999
endif
do i=1, isz
idx = desc_a%ovrlap_elem(i,1)
xs(i) = x(idx)
end do
call psb_erractionrestore(err_act)
return
9999 call psb_error_handler(ictxt,err_act)
return
end subroutine psi_i2ovrl_saver1
subroutine psi_i2ovrl_saver2(x,xs,desc_a,info)
use psi_mod, psi_protect_name => psi_i2ovrl_saver2
use psb_realloc_mod
implicit none
integer(psb_i2pk_), intent(inout) :: x(:,:)
integer(psb_i2pk_), allocatable :: xs(:,:)
type(psb_desc_type), intent(in) :: desc_a
integer(psb_ipk_), intent(out) :: info
! locals
integer(psb_ipk_) :: ictxt, np, me, err_act, i, idx, isz, nc
character(len=20) :: name, ch_err
name='psi_i2ovrl_saver2'
info = psb_success_
call psb_erractionsave(err_act)
if (psb_errstatus_fatal()) then
info = psb_err_internal_error_ ; goto 9999
end if
ictxt = desc_a%get_context()
call psb_info(ictxt, me, np)
if (np == -1) then
info = psb_err_context_error_
call psb_errpush(info,name)
goto 9999
endif
isz = size(desc_a%ovrlap_elem,1)
nc = size(x,2)
call psb_realloc(isz,nc,xs,info)
if (info /= psb_success_) then
info = psb_err_alloc_dealloc_
call psb_errpush(info,name)
goto 9999
endif
do i=1, isz
idx = desc_a%ovrlap_elem(i,1)
xs(i,:) = x(idx,:)
end do
call psb_erractionrestore(err_act)
return
9999 call psb_error_handler(ictxt,err_act)
return
end subroutine psi_i2ovrl_saver2
+173
View File
@@ -0,0 +1,173 @@
!
! 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.
!
!
! Subroutine: psi_i2ovrl_update
! These subroutines update the overlap region of a vector; they are used
! for the transpose matrix-vector product when there is a nonempty overlap,
! or for the application of Additive Schwarz preconditioners.
!
!
subroutine psi_i2ovrl_updr1(x,desc_a,update,info)
use psi_mod, psi_protect_name => psi_i2ovrl_updr1
implicit none
integer(psb_i2pk_), intent(inout), target :: x(:)
type(psb_desc_type), intent(in) :: desc_a
integer(psb_ipk_), intent(in) :: update
integer(psb_ipk_), intent(out) :: info
! locals
integer(psb_ipk_) :: ictxt, np, me, err_act, i, idx, ndm
integer(psb_ipk_) :: ierr(5)
character(len=20) :: name, ch_err
name='psi_i2ovrl_updr1'
info = psb_success_
call psb_erractionsave(err_act)
if (psb_errstatus_fatal()) then
info = psb_err_internal_error_ ; goto 9999
end if
ictxt = desc_a%get_context()
call psb_info(ictxt, me, np)
if (np == -1) then
info = psb_err_context_error_
call psb_errpush(info,name)
goto 9999
endif
! switch on update type
select case (update)
case(psb_square_root_)
do i=1,size(desc_a%ovrlap_elem,1)
idx = desc_a%ovrlap_elem(i,1)
ndm = desc_a%ovrlap_elem(i,2)
x(idx) = x(idx)/sqrt(real(ndm))
end do
case(psb_avg_)
do i=1,size(desc_a%ovrlap_elem,1)
idx = desc_a%ovrlap_elem(i,1)
ndm = desc_a%ovrlap_elem(i,2)
x(idx) = x(idx)/real(ndm)
end do
case(psb_setzero_)
do i=1,size(desc_a%ovrlap_elem,1)
idx = desc_a%ovrlap_elem(i,1)
if (me /= desc_a%ovrlap_elem(i,3))&
& x(idx) = i2zero
end do
case(psb_sum_)
! do nothing
case default
! wrong value for choice argument
info = psb_err_iarg_invalid_value_
ierr(1) = 3; ierr(2)=update;
call psb_errpush(info,name,i_err=ierr)
goto 9999
end select
call psb_erractionrestore(err_act)
return
9999 call psb_error_handler(ictxt,err_act)
return
end subroutine psi_i2ovrl_updr1
subroutine psi_i2ovrl_updr2(x,desc_a,update,info)
use psi_mod, psi_protect_name => psi_i2ovrl_updr2
implicit none
integer(psb_i2pk_), intent(inout), target :: x(:,:)
type(psb_desc_type), intent(in) :: desc_a
integer(psb_ipk_), intent(in) :: update
integer(psb_ipk_), intent(out) :: info
! locals
integer(psb_ipk_) :: ictxt, np, me, err_act, i, idx, ndm
integer(psb_ipk_) :: ierr(5)
character(len=20) :: name, ch_err
name='psi_i2ovrl_updr2'
info = psb_success_
call psb_erractionsave(err_act)
if (psb_errstatus_fatal()) then
info = psb_err_internal_error_ ; goto 9999
end if
ictxt = desc_a%get_context()
call psb_info(ictxt, me, np)
if (np == -1) then
info = psb_err_context_error_
call psb_errpush(info,name)
goto 9999
endif
! switch on update type
select case (update)
case(psb_square_root_)
do i=1,size(desc_a%ovrlap_elem,1)
idx = desc_a%ovrlap_elem(i,1)
ndm = desc_a%ovrlap_elem(i,2)
x(idx,:) = x(idx,:)/sqrt(real(ndm))
end do
case(psb_avg_)
do i=1,size(desc_a%ovrlap_elem,1)
idx = desc_a%ovrlap_elem(i,1)
ndm = desc_a%ovrlap_elem(i,2)
x(idx,:) = x(idx,:)/real(ndm)
end do
case(psb_setzero_)
do i=1,size(desc_a%ovrlap_elem,1)
idx = desc_a%ovrlap_elem(i,1)
if (me /= desc_a%ovrlap_elem(i,3))&
& x(idx,:) = i2zero
end do
case(psb_sum_)
! do nothing
case default
! wrong value for choice argument
info = psb_err_iarg_invalid_value_
ierr(1) = 3; ierr(2)=update;
call psb_errpush(info,name,i_err=ierr)
goto 9999
end select
call psb_erractionrestore(err_act)
return
9999 call psb_error_handler(ictxt,err_act)
return
end subroutine psi_i2ovrl_updr2
+990
View File
@@ -0,0 +1,990 @@
!
! Parallel Sparse BLAS version 3.5
! (C) Copyright 2006-2018
! Salvatore Filippone
! Alfredo Buttari
!
! Redistribution and use in source and binary forms, with or without
! modification, are permitted provided that the following conditions
! are met:
! 1. Redistributions of source code must retain the above copyright
! notice, this list of conditions and the following disclaimer.
! 2. Redistributions in binary form must reproduce the above copyright
! notice, this list of conditions, and the following disclaimer in the
! documentation and/or other materials provided with the distribution.
! 3. The name of the PSBLAS group or the names of its contributors may
! not be used to endorse or promote products derived from this
! software without specific written permission.
!
! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
! ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
! TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
! PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE PSBLAS GROUP OR ITS CONTRIBUTORS
! BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
! CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
! SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
! INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
! CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
! ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
! POSSIBILITY OF SUCH DAMAGE.
!
!
!
! File: psi_i2swapdata.F90
!
! Subroutine: psi_i2swapdatam
! Implements the data exchange among processes. Essentially this is doing
! a variable all-to-all data exchange (ALLTOALLV in MPI parlance), but
! it is capable of pruning empty exchanges, which are very likely in our
! application environment. All the variants have the same structure
! In all these subroutines X may be: I Integer
! S real(psb_spk_)
! D real(psb_dpk_)
! C complex(psb_spk_)
! Z complex(psb_dpk_)
! Basically the operation is as follows: on each process, we identify
! sections SND(Y) and RCV(Y); then we do a send on (PACK(SND(Y)));
! then we receive, and we do an update with Y = UNPACK(RCV(Y)) + BETA * Y
! but only on the elements involved in the UNPACK operation.
! Thus: for halo data exchange, the receive section is confined in the
! halo indices, and BETA=0, whereas for overlap exchange the receive section
! is scattered in the owned indices, and BETA=1.
! The first routine picks the desired exchange index list and passes it to the second.
!
! Arguments:
! flag - integer Choose the algorithm for data exchange:
! this is chosen through bit fields.
! swap_mpi = iand(flag,psb_swap_mpi_) /= 0
! swap_sync = iand(flag,psb_swap_sync_) /= 0
! swap_send = iand(flag,psb_swap_send_) /= 0
! swap_recv = iand(flag,psb_swap_recv_) /= 0
! if (swap_mpi): use underlying MPI_ALLTOALLV.
! if (swap_sync): use PSB_SND and PSB_RCV in
! synchronized pairs
! if (swap_send .and. swap_recv): use mpi_irecv
! and mpi_send
! if (swap_send): use psb_snd (but need another
! call with swap_recv to complete)
! if (swap_recv): use psb_rcv (completing a
! previous call with swap_send)
!
!
! n - integer Number of columns in Y
! beta - integer Choose overwrite or sum.
! y(:,:) - integer The data area
! desc_a - type(psb_desc_type). The communication descriptor.
! work(:) - integer Buffer space. If not sufficient, will do
! our own internal allocation.
! info - integer. return code.
! data - integer which list is to be used to exchange data
! default psb_comm_halo_
! psb_comm_halo_ use halo_index
! psb_comm_ext_ use ext_index
! psb_comm_ovrl_ use ovrl_index
! psb_comm_mov_ use ovr_mst_idx
!
!
subroutine psi_i2swapdatam(flag,n,beta,y,desc_a,work,info,data)
use psi_mod, psb_protect_name => psi_i2swapdatam
use psb_error_mod
use psb_desc_mod
use psb_penv_mod
#ifdef MPI_MOD
use mpi
#endif
implicit none
#ifdef MPI_H
include 'mpif.h'
#endif
integer(psb_ipk_), intent(in) :: flag, n
integer(psb_ipk_), intent(out) :: info
integer(psb_i2pk_) :: y(:,:), beta
integer(psb_i2pk_), target :: work(:)
type(psb_desc_type),target :: desc_a
integer(psb_ipk_), optional :: data
! locals
integer(psb_ipk_) :: ictxt, np, me, icomm, idxs, idxr, totxch, data_, err_act
integer(psb_ipk_), pointer :: d_idx(:)
character(len=20) :: name
info=psb_success_
name='psi_swap_data'
call psb_erractionsave(err_act)
ictxt = desc_a%get_context()
icomm = desc_a%get_mpic()
call psb_info(ictxt,me,np)
if (np == -1) then
info=psb_err_context_error_
call psb_errpush(info,name)
goto 9999
endif
if (.not.psb_is_asb_desc(desc_a)) then
info=psb_err_invalid_cd_state_
call psb_errpush(info,name)
goto 9999
endif
if(present(data)) then
data_ = data
else
data_ = psb_comm_halo_
end if
call desc_a%get_list(data_,d_idx,totxch,idxr,idxs,info)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,a_err='psb_cd_get_list')
goto 9999
end if
call psi_swapdata(ictxt,icomm,flag,n,beta,y,d_idx,totxch,idxs,idxr,work,info)
if (info /= psb_success_) goto 9999
call psb_erractionrestore(err_act)
return
9999 call psb_error_handler(ictxt,err_act)
return
end subroutine psi_i2swapdatam
subroutine psi_i2swapidxm(iictxt,iicomm,flag,n,beta,y,idx, &
& totxch,totsnd,totrcv,work,info)
use psi_mod, psb_protect_name => psi_i2swapidxm
use psb_error_mod
use psb_desc_mod
use psb_penv_mod
#ifdef MPI_MOD
use mpi
#endif
implicit none
#ifdef MPI_H
include 'mpif.h'
#endif
integer(psb_ipk_), intent(in) :: iictxt,iicomm,flag,n
integer(psb_ipk_), intent(out) :: info
integer(psb_i2pk_) :: y(:,:), beta
integer(psb_i2pk_), target :: work(:)
integer(psb_ipk_), intent(in) :: idx(:),totxch,totsnd, totrcv
! locals
integer(psb_mpk_) :: ictxt, icomm, np, me,&
& proc_to_comm, p2ptag, p2pstat(mpi_status_size), iret
integer(psb_mpk_), allocatable, dimension(:) :: bsdidx, brvidx,&
& sdsz, rvsz, prcid, rvhd, sdhd
integer(psb_ipk_) :: nesd, nerv,&
& err_act, i, idx_pt, totsnd_, totrcv_,&
& snd_pt, rcv_pt, pnti
logical :: swap_mpi, swap_sync, swap_send, swap_recv,&
& albf,do_send,do_recv
logical, parameter :: usersend=.false.
integer(psb_i2pk_), pointer, dimension(:) :: sndbuf, rcvbuf
#ifdef HAVE_VOLATILE
volatile :: sndbuf, rcvbuf
#endif
character(len=20) :: name
info=psb_success_
name='psi_swap_data'
call psb_erractionsave(err_act)
ictxt = iictxt
icomm = iicomm
call psb_info(ictxt,me,np)
if (np == -1) then
info=psb_err_context_error_
call psb_errpush(info,name)
goto 9999
endif
swap_mpi = iand(flag,psb_swap_mpi_) /= 0
swap_sync = iand(flag,psb_swap_sync_) /= 0
swap_send = iand(flag,psb_swap_send_) /= 0
swap_recv = iand(flag,psb_swap_recv_) /= 0
do_send = swap_mpi .or. swap_sync .or. swap_send
do_recv = swap_mpi .or. swap_sync .or. swap_recv
totrcv_ = totrcv * n
totsnd_ = totsnd * n
if (swap_mpi) then
allocate(sdsz(0:np-1), rvsz(0:np-1), bsdidx(0:np-1),&
& brvidx(0:np-1), rvhd(0:np-1), sdhd(0:np-1), prcid(0:np-1),&
& stat=info)
if(info /= psb_success_) then
call psb_errpush(psb_err_alloc_dealloc_,name)
goto 9999
end if
rvhd(:) = mpi_request_null
sdsz(:) = 0
rvsz(:) = 0
! prepare info for communications
pnti = 1
snd_pt = 1
rcv_pt = 1
do i=1, totxch
proc_to_comm = idx(pnti+psb_proc_id_)
nerv = idx(pnti+psb_n_elem_recv_)
nesd = idx(pnti+nerv+psb_n_elem_send_)
prcid(proc_to_comm) = psb_get_mpi_rank(ictxt,proc_to_comm)
brvidx(proc_to_comm) = rcv_pt
rvsz(proc_to_comm) = n*nerv
bsdidx(proc_to_comm) = snd_pt
sdsz(proc_to_comm) = n*nesd
rcv_pt = rcv_pt + n*nerv
snd_pt = snd_pt + n*nesd
pnti = pnti + nerv + nesd + 3
end do
else
allocate(rvhd(totxch),prcid(totxch),stat=info)
if(info /= psb_success_) then
call psb_errpush(psb_err_alloc_dealloc_,name)
goto 9999
end if
end if
totrcv_ = max(totrcv_,1)
totsnd_ = max(totsnd_,1)
if((totrcv_+totsnd_) < size(work)) then
sndbuf => work(1:totsnd_)
rcvbuf => work(totsnd_+1:totsnd_+totrcv_)
albf=.false.
else
allocate(sndbuf(totsnd_),rcvbuf(totrcv_), stat=info)
if(info /= psb_success_) then
call psb_errpush(psb_err_alloc_dealloc_,name)
goto 9999
end if
albf=.true.
end if
if (do_send) then
! Pack send buffers
pnti = 1
snd_pt = 1
do i=1, totxch
nerv = idx(pnti+psb_n_elem_recv_)
nesd = idx(pnti+nerv+psb_n_elem_send_)
idx_pt = 1+pnti+nerv+psb_n_elem_send_
call psi_gth(nesd,n,idx(idx_pt:idx_pt+nesd-1),&
& y,sndbuf(snd_pt:snd_pt+n*nesd-1))
snd_pt = snd_pt + n*nesd
pnti = pnti + nerv + nesd + 3
end do
end if
! Case SWAP_MPI
if (swap_mpi) then
! swap elements using mpi_alltoallv
call mpi_alltoallv(sndbuf,sdsz,bsdidx,&
& psb_mpi_i2pk_,rcvbuf,rvsz,&
& brvidx,psb_mpi_i2pk_,icomm,iret)
if(iret /= mpi_success) then
info=psb_err_mpi_error_
call psb_errpush(info,name,m_err=(/iret/))
goto 9999
end if
else if (swap_sync) then
pnti = 1
snd_pt = 1
rcv_pt = 1
do i=1, totxch
proc_to_comm = idx(pnti+psb_proc_id_)
nerv = idx(pnti+psb_n_elem_recv_)
nesd = idx(pnti+nerv+psb_n_elem_send_)
if (proc_to_comm < me) then
if (nesd>0) call psb_snd(ictxt,&
& sndbuf(snd_pt:snd_pt+n*nesd-1), proc_to_comm)
if (nerv>0) call psb_rcv(ictxt,&
& rcvbuf(rcv_pt:rcv_pt+n*nerv-1), proc_to_comm)
else if (proc_to_comm > me) then
if (nerv>0) call psb_rcv(ictxt,&
& rcvbuf(rcv_pt:rcv_pt+n*nerv-1), proc_to_comm)
if (nesd>0) call psb_snd(ictxt,&
& sndbuf(snd_pt:snd_pt+n*nesd-1), proc_to_comm)
else if (proc_to_comm == me) then
if (nesd /= nerv) then
write(psb_err_unit,*) &
& 'Fatal error in swapdata: mismatch on self send',&
& nerv,nesd
end if
rcvbuf(rcv_pt:rcv_pt+n*nerv-1) = sndbuf(snd_pt:snd_pt+n*nesd-1)
end if
rcv_pt = rcv_pt + n*nerv
snd_pt = snd_pt + n*nesd
pnti = pnti + nerv + nesd + 3
end do
else if (swap_send .and. swap_recv) then
! First I post all the non blocking receives
pnti = 1
snd_pt = 1
rcv_pt = 1
do i=1, totxch
proc_to_comm = idx(pnti+psb_proc_id_)
nerv = idx(pnti+psb_n_elem_recv_)
nesd = idx(pnti+nerv+psb_n_elem_send_)
prcid(i) = psb_get_mpi_rank(ictxt,proc_to_comm)
if ((nerv>0).and.(proc_to_comm /= me)) then
p2ptag = psb_int2_swap_tag
call mpi_irecv(rcvbuf(rcv_pt),n*nerv,&
& psb_mpi_i2pk_,prcid(i),&
& p2ptag, icomm,rvhd(i),iret)
end if
rcv_pt = rcv_pt + n*nerv
snd_pt = snd_pt + n*nesd
pnti = pnti + nerv + nesd + 3
end do
! Then I post all the blocking sends
if (usersend) call mpi_barrier(icomm,iret)
pnti = 1
snd_pt = 1
rcv_pt = 1
do i=1, totxch
proc_to_comm = idx(pnti+psb_proc_id_)
nerv = idx(pnti+psb_n_elem_recv_)
nesd = idx(pnti+nerv+psb_n_elem_send_)
p2ptag = psb_int2_swap_tag
if ((nesd>0).and.(proc_to_comm /= me)) then
if (usersend) then
call mpi_rsend(sndbuf(snd_pt),n*nesd,&
& psb_mpi_i2pk_,prcid(i),&
& p2ptag,icomm,iret)
else
call mpi_send(sndbuf(snd_pt),n*nesd,&
& psb_mpi_i2pk_,prcid(i),&
& p2ptag,icomm,iret)
end if
if(iret /= mpi_success) then
info=psb_err_mpi_error_
call psb_errpush(info,name,m_err=(/iret/))
goto 9999
end if
end if
rcv_pt = rcv_pt + n*nerv
snd_pt = snd_pt + n*nesd
pnti = pnti + nerv + nesd + 3
end do
pnti = 1
do i=1, totxch
proc_to_comm = idx(pnti+psb_proc_id_)
nerv = idx(pnti+psb_n_elem_recv_)
nesd = idx(pnti+nerv+psb_n_elem_send_)
p2ptag = psb_int2_swap_tag
if ((proc_to_comm /= me).and.(nerv>0)) then
call mpi_wait(rvhd(i),p2pstat,iret)
if(iret /= mpi_success) then
info=psb_err_mpi_error_
call psb_errpush(info,name,m_err=(/iret/))
goto 9999
end if
else if (proc_to_comm == me) then
if (nesd /= nerv) then
write(psb_err_unit,*)&
& 'Fatal error in swapdata: mismatch on self send', &
& nerv,nesd
end if
rcvbuf(rcv_pt:rcv_pt+n*nerv-1) = sndbuf(snd_pt:snd_pt+n*nesd-1)
end if
pnti = pnti + nerv + nesd + 3
end do
else if (swap_send) then
pnti = 1
snd_pt = 1
rcv_pt = 1
do i=1, totxch
proc_to_comm = idx(pnti+psb_proc_id_)
nerv = idx(pnti+psb_n_elem_recv_)
nesd = idx(pnti+nerv+psb_n_elem_send_)
if (nesd>0) call psb_snd(ictxt,&
& sndbuf(snd_pt:snd_pt+n*nesd-1), proc_to_comm)
rcv_pt = rcv_pt + n*nerv
snd_pt = snd_pt + n*nesd
pnti = pnti + nerv + nesd + 3
end do
else if (swap_recv) then
pnti = 1
snd_pt = 1
rcv_pt = 1
do i=1, totxch
proc_to_comm = idx(pnti+psb_proc_id_)
nerv = idx(pnti+psb_n_elem_recv_)
nesd = idx(pnti+nerv+psb_n_elem_send_)
if (nerv>0) call psb_rcv(ictxt,&
& rcvbuf(rcv_pt:rcv_pt+n*nerv-1), proc_to_comm)
rcv_pt = rcv_pt + n*nerv
snd_pt = snd_pt + n*nesd
pnti = pnti + nerv + nesd + 3
end do
end if
if (do_recv) then
pnti = 1
snd_pt = 1
rcv_pt = 1
do i=1, totxch
proc_to_comm = idx(pnti+psb_proc_id_)
nerv = idx(pnti+psb_n_elem_recv_)
nesd = idx(pnti+nerv+psb_n_elem_send_)
idx_pt = 1+pnti+psb_n_elem_recv_
call psi_sct(nerv,n,idx(idx_pt:idx_pt+nerv-1),&
& rcvbuf(rcv_pt:rcv_pt+n*nerv-1),beta,y)
rcv_pt = rcv_pt + n*nerv
snd_pt = snd_pt + n*nesd
pnti = pnti + nerv + nesd + 3
end do
end if
if (swap_mpi) then
deallocate(sdsz,rvsz,bsdidx,brvidx,rvhd,prcid,sdhd,&
& stat=info)
else
deallocate(rvhd,prcid,stat=info)
end if
if(info /= psb_success_) then
call psb_errpush(psb_err_alloc_dealloc_,name)
goto 9999
end if
if(albf) deallocate(sndbuf,rcvbuf,stat=info)
if(info /= psb_success_) then
call psb_errpush(psb_err_alloc_dealloc_,name)
goto 9999
end if
call psb_erractionrestore(err_act)
return
9999 call psb_error_handler(iictxt,err_act)
return
end subroutine psi_i2swapidxm
!
!
! Subroutine: psi_i2swapdatav
! Implements the data exchange among processes. Essentially this is doing
! a variable all-to-all data exchange (ALLTOALLV in MPI parlance), but
! it is capable of pruning empty exchanges, which are very likely in out
! application environment. All the variants have the same structure
! In all these subroutines X may be: I Integer
! S real(psb_spk_)
! D real(psb_dpk_)
! C complex(psb_spk_)
! Z complex(psb_dpk_)
! Basically the operation is as follows: on each process, we identify
! sections SND(Y) and RCV(Y); then we do a SEND(PACK(SND(Y)));
! then we receive, and we do an update with Y = UNPACK(RCV(Y)) + BETA * Y
! but only on the elements involved in the UNPACK operation.
! Thus: for halo data exchange, the receive section is confined in the
! halo indices, and BETA=0, whereas for overlap exchange the receive section
! is scattered in the owned indices, and BETA=1.
! The first routine picks the desired exchange index list and passes it to the second.
!
! Arguments:
! flag - integer Choose the algorithm for data exchange:
! this is chosen through bit fields.
! swap_mpi = iand(flag,psb_swap_mpi_) /= 0
! swap_sync = iand(flag,psb_swap_sync_) /= 0
! swap_send = iand(flag,psb_swap_send_) /= 0
! swap_recv = iand(flag,psb_swap_recv_) /= 0
! if (swap_mpi): use underlying MPI_ALLTOALLV.
! if (swap_sync): use PSB_SND and PSB_RCV in
! synchronized pairs
! if (swap_send .and. swap_recv): use mpi_irecv
! and mpi_send
! if (swap_send): use psb_snd (but need another
! call with swap_recv to complete)
! if (swap_recv): use psb_rcv (completing a
! previous call with swap_send)
!
!
! n - integer Number of columns in Y
! beta - integer Choose overwrite or sum.
! y(:) - integer The data area
! desc_a - type(psb_desc_type). The communication descriptor.
! work(:) - integer Buffer space. If not sufficient, will do
! our own internal allocation.
! info - integer. return code.
! data - integer which list is to be used to exchange data
! default psb_comm_halo_
! psb_comm_halo_ use halo_index
! psb_comm_ext_ use ext_index
! psb_comm_ovrl_ use ovrl_index
! psb_comm_mov_ use ovr_mst_idx
!
!
subroutine psi_i2swapdatav(flag,beta,y,desc_a,work,info,data)
use psi_mod, psb_protect_name => psi_i2swapdatav
use psb_error_mod
use psb_desc_mod
use psb_penv_mod
#ifdef MPI_MOD
use mpi
#endif
implicit none
#ifdef MPI_H
include 'mpif.h'
#endif
integer(psb_ipk_), intent(in) :: flag
integer(psb_ipk_), intent(out) :: info
integer(psb_i2pk_) :: y(:), beta
integer(psb_i2pk_), target :: work(:)
type(psb_desc_type),target :: desc_a
integer(psb_ipk_), optional :: data
! locals
integer(psb_ipk_) :: ictxt, np, me, icomm, idxs, idxr, totxch, data_, err_act
integer(psb_ipk_), pointer :: d_idx(:)
character(len=20) :: name
info=psb_success_
name='psi_swap_datav'
call psb_erractionsave(err_act)
ictxt = desc_a%get_context()
icomm = desc_a%get_mpic()
call psb_info(ictxt,me,np)
if (np == -1) then
info=psb_err_context_error_
call psb_errpush(info,name)
goto 9999
endif
if (.not.psb_is_asb_desc(desc_a)) then
info=psb_err_invalid_cd_state_
call psb_errpush(info,name)
goto 9999
endif
if (present(data)) then
data_ = data
else
data_ = psb_comm_halo_
end if
call desc_a%get_list(data_,d_idx,totxch,idxr,idxs,info)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,a_err='psb_cd_get_list')
goto 9999
end if
call psi_swapdata(ictxt,icomm,flag,beta,y,d_idx,totxch,idxs,idxr,work,info)
if (info /= psb_success_) goto 9999
call psb_erractionrestore(err_act)
return
9999 call psb_error_handler(ictxt,err_act)
return
end subroutine psi_i2swapdatav
!
!
! Subroutine: psi_i2swapdataidxv
! Does the data exchange among processes.
!
! The real workhorse: the outer routines will only choose the index list
! this one takes the index list and does the actual exchange.
!
!
!
subroutine psi_i2swapidxv(iictxt,iicomm,flag,beta,y,idx, &
& totxch,totsnd,totrcv,work,info)
use psi_mod, psb_protect_name => psi_i2swapidxv
use psb_error_mod
use psb_desc_mod
use psb_penv_mod
#ifdef MPI_MOD
use mpi
#endif
implicit none
#ifdef MPI_H
include 'mpif.h'
#endif
integer(psb_ipk_), intent(in) :: iictxt,iicomm,flag
integer(psb_ipk_), intent(out) :: info
integer(psb_i2pk_) :: y(:), beta
integer(psb_i2pk_), target :: work(:)
integer(psb_ipk_), intent(in) :: idx(:),totxch,totsnd, totrcv
! locals
integer(psb_mpk_) :: ictxt, icomm, np, me,&
& proc_to_comm, p2ptag, p2pstat(mpi_status_size), iret
integer(psb_mpk_), allocatable, dimension(:) :: bsdidx, brvidx,&
& sdsz, rvsz, prcid, rvhd, sdhd
integer(psb_ipk_) :: nesd, nerv,&
& err_act, i, idx_pt, totsnd_, totrcv_,&
& snd_pt, rcv_pt, pnti, n
logical :: swap_mpi, swap_sync, swap_send, swap_recv,&
& albf,do_send,do_recv
logical, parameter :: usersend=.false.
integer(psb_i2pk_), pointer, dimension(:) :: sndbuf, rcvbuf
#ifdef HAVE_VOLATILE
volatile :: sndbuf, rcvbuf
#endif
character(len=20) :: name
info=psb_success_
name='psi_swap_datav'
call psb_erractionsave(err_act)
ictxt = iictxt
icomm = iicomm
call psb_info(ictxt,me,np)
if (np == -1) then
info=psb_err_context_error_
call psb_errpush(info,name)
goto 9999
endif
n=1
swap_mpi = iand(flag,psb_swap_mpi_) /= 0
swap_sync = iand(flag,psb_swap_sync_) /= 0
swap_send = iand(flag,psb_swap_send_) /= 0
swap_recv = iand(flag,psb_swap_recv_) /= 0
do_send = swap_mpi .or. swap_sync .or. swap_send
do_recv = swap_mpi .or. swap_sync .or. swap_recv
totrcv_ = totrcv * n
totsnd_ = totsnd * n
if (swap_mpi) then
allocate(sdsz(0:np-1), rvsz(0:np-1), bsdidx(0:np-1),&
& brvidx(0:np-1), rvhd(0:np-1), sdhd(0:np-1), prcid(0:np-1),&
& stat=info)
if(info /= psb_success_) then
call psb_errpush(psb_err_alloc_dealloc_,name)
goto 9999
end if
rvhd(:) = mpi_request_null
sdsz(:) = 0
rvsz(:) = 0
! prepare info for communications
pnti = 1
snd_pt = 1
rcv_pt = 1
do i=1, totxch
proc_to_comm = idx(pnti+psb_proc_id_)
nerv = idx(pnti+psb_n_elem_recv_)
nesd = idx(pnti+nerv+psb_n_elem_send_)
prcid(proc_to_comm) = psb_get_mpi_rank(ictxt,proc_to_comm)
brvidx(proc_to_comm) = rcv_pt
rvsz(proc_to_comm) = nerv
bsdidx(proc_to_comm) = snd_pt
sdsz(proc_to_comm) = nesd
rcv_pt = rcv_pt + nerv
snd_pt = snd_pt + nesd
pnti = pnti + nerv + nesd + 3
end do
else
allocate(rvhd(totxch),prcid(totxch),stat=info)
if(info /= psb_success_) then
call psb_errpush(psb_err_alloc_dealloc_,name)
goto 9999
end if
end if
totrcv_ = max(totrcv_,1)
totsnd_ = max(totsnd_,1)
if((totrcv_+totsnd_) < size(work)) then
sndbuf => work(1:totsnd_)
rcvbuf => work(totsnd_+1:totsnd_+totrcv_)
albf=.false.
else
allocate(sndbuf(totsnd_),rcvbuf(totrcv_), stat=info)
if(info /= psb_success_) then
call psb_errpush(psb_err_alloc_dealloc_,name)
goto 9999
end if
albf=.true.
end if
if (do_send) then
! Pack send buffers
pnti = 1
snd_pt = 1
do i=1, totxch
nerv = idx(pnti+psb_n_elem_recv_)
nesd = idx(pnti+nerv+psb_n_elem_send_)
idx_pt = 1+pnti+nerv+psb_n_elem_send_
call psi_gth(nesd,idx(idx_pt:idx_pt+nesd-1),&
& y,sndbuf(snd_pt:snd_pt+nesd-1))
snd_pt = snd_pt + nesd
pnti = pnti + nerv + nesd + 3
end do
end if
! Case SWAP_MPI
if (swap_mpi) then
! swap elements using mpi_alltoallv
call mpi_alltoallv(sndbuf,sdsz,bsdidx,&
& psb_mpi_i2pk_,rcvbuf,rvsz,&
& brvidx,psb_mpi_i2pk_,icomm,iret)
if(iret /= mpi_success) then
info=psb_err_mpi_error_
call psb_errpush(info,name,m_err=(/iret/))
goto 9999
end if
else if (swap_sync) then
pnti = 1
snd_pt = 1
rcv_pt = 1
do i=1, totxch
proc_to_comm = idx(pnti+psb_proc_id_)
nerv = idx(pnti+psb_n_elem_recv_)
nesd = idx(pnti+nerv+psb_n_elem_send_)
if (proc_to_comm < me) then
if (nesd>0) call psb_snd(ictxt,&
& sndbuf(snd_pt:snd_pt+nesd-1), proc_to_comm)
if (nerv>0) call psb_rcv(ictxt,&
& rcvbuf(rcv_pt:rcv_pt+nerv-1), proc_to_comm)
else if (proc_to_comm > me) then
if (nerv>0) call psb_rcv(ictxt,&
& rcvbuf(rcv_pt:rcv_pt+nerv-1), proc_to_comm)
if (nesd>0) call psb_snd(ictxt,&
& sndbuf(snd_pt:snd_pt+nesd-1), proc_to_comm)
else if (proc_to_comm == me) then
if (nesd /= nerv) then
write(psb_err_unit,*) &
& 'Fatal error in swapdata: mismatch on self send', &
& nerv,nesd
end if
rcvbuf(rcv_pt:rcv_pt+nerv-1) = sndbuf(snd_pt:snd_pt+nesd-1)
end if
rcv_pt = rcv_pt + nerv
snd_pt = snd_pt + nesd
pnti = pnti + nerv + nesd + 3
end do
else if (swap_send .and. swap_recv) then
! First I post all the non blocking receives
pnti = 1
snd_pt = 1
rcv_pt = 1
do i=1, totxch
proc_to_comm = idx(pnti+psb_proc_id_)
nerv = idx(pnti+psb_n_elem_recv_)
nesd = idx(pnti+nerv+psb_n_elem_send_)
prcid(i) = psb_get_mpi_rank(ictxt,proc_to_comm)
if ((nerv>0).and.(proc_to_comm /= me)) then
p2ptag = psb_int2_swap_tag
call mpi_irecv(rcvbuf(rcv_pt),nerv,&
& psb_mpi_i2pk_,prcid(i),&
& p2ptag, icomm,rvhd(i),iret)
end if
rcv_pt = rcv_pt + nerv
snd_pt = snd_pt + nesd
pnti = pnti + nerv + nesd + 3
end do
! Then I post all the blocking sends
if (usersend) call mpi_barrier(icomm,iret)
pnti = 1
snd_pt = 1
rcv_pt = 1
do i=1, totxch
proc_to_comm = idx(pnti+psb_proc_id_)
nerv = idx(pnti+psb_n_elem_recv_)
nesd = idx(pnti+nerv+psb_n_elem_send_)
p2ptag = psb_int2_swap_tag
if ((nesd>0).and.(proc_to_comm /= me)) then
if (usersend) then
call mpi_rsend(sndbuf(snd_pt),nesd,&
& psb_mpi_i2pk_,prcid(i),&
& p2ptag,icomm,iret)
else
call mpi_send(sndbuf(snd_pt),nesd,&
& psb_mpi_i2pk_,prcid(i),&
& p2ptag,icomm,iret)
end if
if(iret /= mpi_success) then
info=psb_err_mpi_error_
call psb_errpush(info,name,m_err=(/iret/))
goto 9999
end if
end if
rcv_pt = rcv_pt + nerv
snd_pt = snd_pt + nesd
pnti = pnti + nerv + nesd + 3
end do
pnti = 1
do i=1, totxch
proc_to_comm = idx(pnti+psb_proc_id_)
nerv = idx(pnti+psb_n_elem_recv_)
nesd = idx(pnti+nerv+psb_n_elem_send_)
p2ptag = psb_int2_swap_tag
if ((proc_to_comm /= me).and.(nerv>0)) then
call mpi_wait(rvhd(i),p2pstat,iret)
if(iret /= mpi_success) then
info=psb_err_mpi_error_
call psb_errpush(info,name,m_err=(/iret/))
goto 9999
end if
else if (proc_to_comm == me) then
if (nesd /= nerv) then
write(psb_err_unit,*) &
& 'Fatal error in swapdata: mismatch on self send', &
& nerv,nesd
end if
rcvbuf(rcv_pt:rcv_pt+nerv-1) = sndbuf(snd_pt:snd_pt+nesd-1)
end if
pnti = pnti + nerv + nesd + 3
end do
else if (swap_send) then
pnti = 1
snd_pt = 1
rcv_pt = 1
do i=1, totxch
proc_to_comm = idx(pnti+psb_proc_id_)
nerv = idx(pnti+psb_n_elem_recv_)
nesd = idx(pnti+nerv+psb_n_elem_send_)
if (nesd>0) call psb_snd(ictxt,&
& sndbuf(snd_pt:snd_pt+nesd-1), proc_to_comm)
rcv_pt = rcv_pt + nerv
snd_pt = snd_pt + nesd
pnti = pnti + nerv + nesd + 3
end do
else if (swap_recv) then
pnti = 1
snd_pt = 1
rcv_pt = 1
do i=1, totxch
proc_to_comm = idx(pnti+psb_proc_id_)
nerv = idx(pnti+psb_n_elem_recv_)
nesd = idx(pnti+nerv+psb_n_elem_send_)
if (nerv>0) call psb_rcv(ictxt,&
& rcvbuf(rcv_pt:rcv_pt+nerv-1), proc_to_comm)
rcv_pt = rcv_pt + nerv
snd_pt = snd_pt + nesd
pnti = pnti + nerv + nesd + 3
end do
end if
if (do_recv) then
pnti = 1
snd_pt = 1
rcv_pt = 1
do i=1, totxch
proc_to_comm = idx(pnti+psb_proc_id_)
nerv = idx(pnti+psb_n_elem_recv_)
nesd = idx(pnti+nerv+psb_n_elem_send_)
idx_pt = 1+pnti+psb_n_elem_recv_
call psi_sct(nerv,idx(idx_pt:idx_pt+nerv-1),&
& rcvbuf(rcv_pt:rcv_pt+nerv-1),beta,y)
rcv_pt = rcv_pt + nerv
snd_pt = snd_pt + nesd
pnti = pnti + nerv + nesd + 3
end do
end if
if (swap_mpi) then
deallocate(sdsz,rvsz,bsdidx,brvidx,rvhd,prcid,sdhd,&
& stat=info)
else
deallocate(rvhd,prcid,stat=info)
end if
if(info /= psb_success_) then
call psb_errpush(psb_err_alloc_dealloc_,name)
goto 9999
end if
if(albf) deallocate(sndbuf,rcvbuf,stat=info)
if(info /= psb_success_) then
call psb_errpush(psb_err_alloc_dealloc_,name)
goto 9999
end if
call psb_erractionrestore(err_act)
return
9999 call psb_error_handler(iictxt,err_act)
return
end subroutine psi_i2swapidxv
File diff suppressed because it is too large Load Diff
+335
View File
@@ -0,0 +1,335 @@
!
! Parallel Sparse BLAS version 3.5
! (C) Copyright 2006-2018
! Salvatore Filippone
! Alfredo Buttari
!
! Redistribution and use in source and binary forms, with or without
! modification, are permitted provided that the following conditions
! are met:
! 1. Redistributions of source code must retain the above copyright
! notice, this list of conditions and the following disclaimer.
! 2. Redistributions in binary form must reproduce the above copyright
! notice, this list of conditions, and the following disclaimer in the
! documentation and/or other materials provided with the distribution.
! 3. The name of the PSBLAS group or the names of its contributors may
! not be used to endorse or promote products derived from this
! software without specific written permission.
!
! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
! ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
! TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
! PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE PSBLAS GROUP OR ITS CONTRIBUTORS
! BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
! CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
! SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
! INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
! CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
! ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
! POSSIBILITY OF SUCH DAMAGE.
!
!
! File: psb_i2gather.f90
!
! Subroutine: psb_i2gather
! This subroutine gathers pieces of a distributed dense matrix into a local one.
!
! Arguments:
! globx - integer,dimension(:,:). The local matrix into which gather
! the distributed pieces.
! locx - integer,dimension(:,:). The local piece of the distributed
! matrix to be gathered.
! desc_a - type(psb_desc_type). The communication descriptor.
! info - integer. Error code.
! iroot - integer. The process that has to own the
! global matrix. If -1 all
! the processes will have a copy.
!
subroutine psb_i2gatherm(globx, locx, desc_a, info, iroot)
use psb_base_mod, psb_protect_name => psb_i2gatherm
implicit none
integer(psb_i2pk_), intent(in) :: locx(:,:)
integer(psb_i2pk_), intent(out), allocatable :: globx(:,:)
type(psb_desc_type), intent(in) :: desc_a
integer(psb_ipk_), intent(out) :: info
integer(psb_ipk_), intent(in), optional :: iroot
! locals
integer(psb_mpk_) :: ictxt, np, me, root, iiroot, icomm, myrank, rootrank
integer(psb_ipk_) :: ierr(5), err_act, lda_locx, lda_globx, lock, globk,&
& maxk, k, jlx, ilx, i, j
integer(psb_lpk_) :: m, n, ilocx, jlocx, idx, iglobx, jglobx
character(len=20) :: name, ch_err
name='psb_i2gatherm'
info=psb_success_
call psb_erractionsave(err_act)
if (psb_errstatus_fatal()) then
info = psb_err_internal_error_ ; goto 9999
end if
ictxt=desc_a%get_context()
! check on blacs grid
call psb_info(ictxt, me, np)
if (np == -1) then
info = psb_err_context_error_
call psb_errpush(info,name)
goto 9999
endif
if (present(iroot)) then
root = iroot
if((root < -1).or.(root > np)) then
info=psb_err_input_value_invalid_i_
ierr(1) = 5; ierr(2)=root
call psb_errpush(info,name,i_err=ierr)
goto 9999
end if
else
root = -1
end if
if (root == -1) then
iiroot = psb_root_
else
iiroot = root
endif
iglobx = 1
jglobx = 1
ilocx = 1
jlocx = 1
m = desc_a%get_global_rows()
n = desc_a%get_global_cols()
lda_globx = m
lda_locx = size(locx, 1)
lock = size(locx,2)
maxk = lock
k = maxk
call psb_bcast(ictxt,k,root=iiroot)
! there should be a global check on k here!!!
call psb_chkglobvect(m,n,lda_globx,iglobx,jglobx,desc_a,info)
if (info == psb_success_) &
& call psb_chkvect(m,n,lda_locx,ilocx,jlocx,desc_a,info,ilx,jlx)
if(info /= psb_success_) then
info=psb_err_from_subroutine_
ch_err='psb_chk(glob)vect'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
if ((ilx /= 1).or.(iglobx /= 1)) then
info=psb_err_ix_n1_iy_n1_unsupported_
call psb_errpush(info,name)
goto 9999
end if
call psb_realloc(m,k,globx,info)
if (info /= psb_success_) then
info=psb_err_alloc_dealloc_
call psb_errpush(info,name)
goto 9999
end if
globx(:,:)=i2zero
do j=1,k
do i=1,desc_a%get_local_rows()
call psb_loc_to_glob(i,idx,desc_a,info)
globx(idx,j) = locx(i,jlx+j-1)
end do
end do
do j=1,k
! adjust overlapped elements
do i=1, size(desc_a%ovrlap_elem,1)
if (me /= desc_a%ovrlap_elem(i,3)) then
idx = desc_a%ovrlap_elem(i,1)
call psb_loc_to_glob(idx,desc_a,info)
globx(idx,j) = i2zero
end if
end do
end do
call psb_sum(ictxt,globx(1:m,1:k),root=root)
call psb_erractionrestore(err_act)
return
9999 call psb_error_handler(ione*ictxt,err_act)
return
end subroutine psb_i2gatherm
!!$
!!$ Parallel Sparse BLAS version 3.5
!!$ (C) Copyright 2006-2018
!!$ Salvatore Filippone University of Rome Tor Vergata
!!$ 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.
!!$
!!$
! Subroutine: psb_i2gatherv
! This subroutine gathers pieces of a distributed dense vector into a local one.
!
! Arguments:
! globx - integer,dimension(:). The local vector into which gather
! the distributed pieces.
! locx - integer,dimension(:). The local piece of the distributed
! vector to be gathered.
! desc_a - type(psb_desc_type). The communication descriptor.
! info - integer. Error code.
! iroot - integer. The process that has to own the
! global matrix. If -1 all
! the processes will have a copy.
! default: -1
!
subroutine psb_i2gatherv(globx, locx, desc_a, info, iroot)
use psb_base_mod, psb_protect_name => psb_i2gatherv
implicit none
integer(psb_i2pk_), intent(in) :: locx(:)
integer(psb_i2pk_), intent(out), allocatable :: globx(:)
type(psb_desc_type), intent(in) :: desc_a
integer(psb_ipk_), intent(out) :: info
integer(psb_ipk_), intent(in), optional :: iroot
! locals
integer(psb_mpk_) :: ictxt, np, me, root, iiroot, icomm, myrank, rootrank
integer(psb_ipk_) :: ierr(5), err_act, lda_locx, lda_globx, lock, globk,&
& maxk, k, jlx, ilx, i, j
integer(psb_lpk_) :: m, n, ilocx, jlocx, idx, iglobx, jglobx
character(len=20) :: name, ch_err
name='psb_i2gatherv'
info=psb_success_
call psb_erractionsave(err_act)
if (psb_errstatus_fatal()) then
info = psb_err_internal_error_ ; goto 9999
end if
ictxt=desc_a%get_context()
! check on blacs grid
call psb_info(ictxt, me, np)
if (np == -1) then
info = psb_err_context_error_
call psb_errpush(info,name)
goto 9999
endif
if (present(iroot)) then
root = iroot
if((root < -1).or.(root > np)) then
info=psb_err_input_value_invalid_i_
ierr(1)=5; ierr(2)=root
call psb_errpush(info,name,i_err=ierr)
goto 9999
end if
else
root = -1
end if
jglobx=1
iglobx = 1
jlocx=1
ilocx = 1
m = desc_a%get_global_rows()
n = desc_a%get_global_cols()
lda_globx = m
lda_locx = size(locx)
k = 1
! there should be a global check on k here!!!
call psb_chkglobvect(m,n,lda_globx,iglobx,jglobx,desc_a,info)
if (info == psb_success_) &
& call psb_chkvect(m,n,lda_locx,ilocx,jlocx,desc_a,info,ilx,jlx)
if(info /= psb_success_) then
info=psb_err_from_subroutine_
ch_err='psb_chk(glob)vect'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
if ((ilx /= 1).or.(iglobx /= 1)) then
info=psb_err_ix_n1_iy_n1_unsupported_
call psb_errpush(info,name)
goto 9999
end if
call psb_realloc(m,globx,info)
if (info /= psb_success_) then
info=psb_err_alloc_dealloc_
call psb_errpush(info,name)
goto 9999
end if
globx(:)=i2zero
do i=1,desc_a%get_local_rows()
call psb_loc_to_glob(i,idx,desc_a,info)
globx(idx) = locx(i)
end do
! adjust overlapped elements
do i=1, size(desc_a%ovrlap_elem,1)
if (me /= desc_a%ovrlap_elem(i,3)) then
idx = desc_a%ovrlap_elem(i,1)
call psb_loc_to_glob(idx,desc_a,info)
globx(idx) = i2zero
end if
end do
call psb_sum(ictxt,globx(1:m),root=root)
call psb_erractionrestore(err_act)
return
9999 call psb_error_handler(ione*ictxt,err_act)
return
end subroutine psb_i2gatherv
+380
View File
@@ -0,0 +1,380 @@
!
! Parallel Sparse BLAS version 3.5
! (C) Copyright 2006-2018
! Salvatore Filippone
! Alfredo Buttari
!
! Redistribution and use in source and binary forms, with or without
! modification, are permitted provided that the following conditions
! are met:
! 1. Redistributions of source code must retain the above copyright
! notice, this list of conditions and the following disclaimer.
! 2. Redistributions in binary form must reproduce the above copyright
! notice, this list of conditions, and the following disclaimer in the
! documentation and/or other materials provided with the distribution.
! 3. The name of the PSBLAS group or the names of its contributors may
! not be used to endorse or promote products derived from this
! software without specific written permission.
!
! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
! ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
! TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
! PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE PSBLAS GROUP OR ITS CONTRIBUTORS
! BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
! CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
! SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
! INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
! CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
! ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
! POSSIBILITY OF SUCH DAMAGE.
!
!
! File: psb_i2halo.f90
!
! Subroutine: psb_i2halom
! This subroutine performs the exchange of the halo elements in a
! distributed dense matrix between all the processes.
!
! Arguments:
! x - integer,dimension(:,:). The local part of the dense matrix.
! desc_a - type(psb_desc_type). The communication descriptor.
! info - integer. Return code
! jx - integer(optional). The starting column of the global matrix.
! ik - integer(optional). The number of columns to gather.
! work - integer(optional). Work area.
! tran - character(optional). Transpose exchange.
! mode - integer(optional). Communication mode (see Swapdata)
! data - integer Which index list in desc_a should be used
! to retrieve rows, default psb_comm_halo_
! psb_comm_halo_ use halo_index
! psb_comm_ext_ use ext_index
! psb_comm_ovrl_ use ovrl_index
! psb_comm_mov_ use ovr_mst_idx
!
!
subroutine psb_i2halom(x,desc_a,info,jx,ik,work,tran,mode,data)
use psb_base_mod, psb_protect_name => psb_i2halom
use psi_mod
implicit none
integer(psb_i2pk_), intent(inout), target :: x(:,:)
type(psb_desc_type), intent(in) :: desc_a
integer(psb_ipk_), intent(out) :: info
integer(psb_i2pk_), optional, target, intent(inout) :: work(:)
integer(psb_ipk_), intent(in), optional :: mode,jx,ik,data
character, intent(in), optional :: tran
! locals
integer(psb_mpk_) :: ictxt, np, me
integer(psb_ipk_) :: err_act, iix, jjx, k, maxk, nrow, imode, i,&
& liwork,data_, ldx
integer(psb_lpk_) :: m, n, ix, ijx
integer(psb_i2pk_),pointer :: iwork(:), xp(:,:)
character :: tran_
character(len=20) :: name, ch_err
logical :: aliw
name='psb_i2halom'
info=psb_success_
call psb_erractionsave(err_act)
if (psb_errstatus_fatal()) then
info = psb_err_internal_error_ ; goto 9999
end if
ictxt=desc_a%get_context()
! check on blacs grid
call psb_info(ictxt, me, np)
if (np == -1) then
info = psb_err_context_error_
call psb_errpush(info,name)
goto 9999
endif
ix = 1
if (present(jx)) then
ijx = jx
else
ijx = 1
endif
m = desc_a%get_global_rows()
n = desc_a%get_global_cols()
nrow = desc_a%get_local_rows()
maxk=size(x,2)-ijx+1
if(present(ik)) then
if(ik > maxk) then
k=maxk
else
k=ik
end if
else
k = maxk
end if
if (present(tran)) then
tran_ = psb_toupper(tran)
else
tran_ = 'N'
endif
if (present(mode)) then
imode = mode
else
imode = IOR(psb_swap_send_,psb_swap_recv_)
endif
if (present(data)) then
data_ = data
else
data_ = psb_comm_halo_
endif
ldx = size(x,1)
! check vector correctness
call psb_chkvect(m,lone,ldx,ix,ijx,desc_a,info,iix,jjx,check_halo=.true.)
if(info /= psb_success_) then
info=psb_err_from_subroutine_ ; ch_err='psb_chkvect'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
liwork=nrow
if (present(work)) then
if(size(work) >= liwork) then
aliw=.false.
iwork => work
else
aliw=.true.
allocate(iwork(liwork),stat=info)
if(info /= psb_success_) then
info=psb_err_from_subroutine_
ch_err='psb_realloc'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
end if
else
aliw=.true.
allocate(iwork(liwork),stat=info)
if(info /= psb_success_) then
info=psb_err_from_subroutine_
ch_err='psb_realloc'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
end if
! exchange halo elements
xp => x(iix:size(x,1),jjx:jjx+k-1)
if(tran_ == 'N') then
call psi_swapdata(imode,k,i2zero,xp,&
& desc_a,iwork,info,data=data_)
else if((tran_ == 'T').or.(tran_ == 'C')) then
call psi_swaptran(imode,k,i2one,xp,&
&desc_a,iwork,info)
else
info = psb_err_internal_error_
call psb_errpush(info,name,a_err='invalid tran')
goto 9999
end if
if(info /= psb_success_) then
ch_err='PSI_cswapdata'
call psb_errpush(psb_err_from_subroutine_,name,a_err=ch_err)
goto 9999
end if
if (aliw) deallocate(iwork)
nullify(iwork)
call psb_erractionrestore(err_act)
return
9999 call psb_error_handler(ione*ictxt,err_act)
return
end subroutine psb_i2halom
!!$
!!$ Parallel Sparse BLAS version 3.5
!!$ (C) Copyright 2006-2018
!!$ Salvatore Filippone University of Rome Tor Vergata
!!$ 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.
!!$
!!$
!
! Subroutine: psb_i2halov
! This subroutine performs the exchange of the halo elements in a
! distributed dense vector between all the processes.
!
! Arguments:
! x - real,dimension(:). The local part of the dense vector.
! desc_a - type(psb_desc_type). The communication descriptor.
! info - integer. Return code
! jx - integer(optional). The starting column of the global matrix.
! ik - integer(optional). The number of columns to gather.
! work - integer(optional). Work area.
! tran - character(optional). Transpose exchange.
! mode - integer(optional). Communication mode (see Swapdata)
! data - integer Which index list in desc_a should be used
! to retrieve rows, default psb_comm_halo_
! psb_comm_halo_ use halo_index
! psb_comm_ext_ use ext_index
! psb_comm_ovrl_ use ovrl_index
! psb_comm_mov_ use ovr_mst_idx
!
!
subroutine psb_i2halov(x,desc_a,info,work,tran,mode,data)
use psb_base_mod, psb_protect_name => psb_i2halov
use psi_mod
implicit none
integer(psb_i2pk_), intent(inout) :: x(:)
type(psb_desc_type), intent(in) :: desc_a
integer(psb_ipk_), intent(out) :: info
integer(psb_i2pk_), target, optional, intent(inout) :: work(:)
integer(psb_ipk_), intent(in), optional :: mode,data
character, intent(in), optional :: tran
! locals
integer(psb_mpk_) :: ictxt, np, me
integer(psb_ipk_) :: err_act, ldx, iix, jjx, nrow, imode, liwork,data_
integer(psb_lpk_) :: m, n, ix, ijx
integer(psb_i2pk_),pointer :: iwork(:)
character :: tran_
character(len=20) :: name, ch_err
logical :: aliw
name='psb_i2halov'
info=psb_success_
call psb_erractionsave(err_act)
if (psb_errstatus_fatal()) then
info = psb_err_internal_error_ ; goto 9999
end if
ictxt=desc_a%get_context()
! check on blacs grid
call psb_info(ictxt, me, np)
if (np == -1) then
info = psb_err_context_error_
call psb_errpush(info,name)
goto 9999
endif
ix = 1
ijx = 1
m = desc_a%get_global_rows()
n = desc_a%get_global_cols()
nrow = desc_a%get_local_rows()
if (present(tran)) then
tran_ = psb_toupper(tran)
else
tran_ = 'N'
endif
if (present(data)) then
data_ = data
else
data_ = psb_comm_halo_
endif
if (present(mode)) then
imode = mode
else
imode = IOR(psb_swap_send_,psb_swap_recv_)
endif
ldx = size(x,1)
! check vector correctness
call psb_chkvect(m,lone,ldx,ix,ijx,desc_a,info,iix,jjx,check_halo=.true.)
if(info /= psb_success_) then
info=psb_err_from_subroutine_ ; ch_err='psb_chkvect'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
liwork=nrow
if (present(work)) then
if(size(work) >= liwork) then
aliw=.false.
iwork => work
else
aliw=.true.
allocate(iwork(liwork),stat=info)
if(info /= psb_success_) then
info=psb_err_from_subroutine_
ch_err='psb_realloc'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
end if
else
aliw=.true.
allocate(iwork(liwork),stat=info)
if(info /= psb_success_) then
info=psb_err_from_subroutine_
ch_err='psb_realloc'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
end if
! exchange halo elements
if(tran_ == 'N') then
call psi_swapdata(imode,i2zero,x(iix:size(x)),&
& desc_a,iwork,info,data=data_)
else if((tran_ == 'T').or.(tran_ == 'C')) then
call psi_swaptran(imode,i2one,x(iix:size(x)),&
& desc_a,iwork,info)
else
info = psb_err_internal_error_
call psb_errpush(info,name,a_err='invalid tran')
goto 9999
end if
if(info /= psb_success_) then
ch_err='PSI_swapdata'
call psb_errpush(psb_err_from_subroutine_,name,a_err=ch_err)
goto 9999
end if
if (aliw) deallocate(iwork)
nullify(iwork)
call psb_erractionrestore(err_act)
return
9999 call psb_error_handler(ione*ictxt,err_act)
return
end subroutine psb_i2halov
+374
View File
@@ -0,0 +1,374 @@
!
! Parallel Sparse BLAS version 3.5
! (C) Copyright 2006-2018
! Salvatore Filippone
! Alfredo Buttari
!
! Redistribution and use in source and binary forms, with or without
! modification, are permitted provided that the following conditions
! are met:
! 1. Redistributions of source code must retain the above copyright
! notice, this list of conditions and the following disclaimer.
! 2. Redistributions in binary form must reproduce the above copyright
! notice, this list of conditions, and the following disclaimer in the
! documentation and/or other materials provided with the distribution.
! 3. The name of the PSBLAS group or the names of its contributors may
! not be used to endorse or promote products derived from this
! software without specific written permission.
!
! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
! ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
! TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
! PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE PSBLAS GROUP OR ITS CONTRIBUTORS
! BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
! CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
! SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
! INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
! CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
! ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
! POSSIBILITY OF SUCH DAMAGE.
!
!
!
! File: psb_i2ovrl.f90
!
! Subroutine: psb_i2ovrlm
! This subroutine performs the exchange of the overlap elements in a
! distributed dense matrix between all the processes.
!
! Arguments:
! x(:,:) - integer The local part of the dense matrix.
! desc_a - type(psb_desc_type). The communication descriptor.
! info - integer. Return code.
! jx - integer(optional). The starting column of the global matrix
! ik - integer(optional). The number of columns to gather.
! work - integer(optional). A work area.
! update - integer(optional). Type of update:
! psb_none_ do nothing
! psb_sum_ sum of overlaps
! psb_avg_ average of overlaps
! mode - integer(optional). Choose the algorithm for data exchange:
! this is chosen through bit fields.
! - swap_mpi = iand(flag,psb_swap_mpi_) /= 0
! - swap_sync = iand(flag,psb_swap_sync_) /= 0
! - swap_send = iand(flag,psb_swap_send_) /= 0
! - swap_recv = iand(flag,psb_swap_recv_) /= 0
! - if (swap_mpi): use underlying MPI_ALLTOALLV.
! - if (swap_sync): use PSB_SND and PSB_RCV in
! synchronized pairs
! - if (swap_send .and. swap_recv): use mpi_irecv
! and mpi_send
! - if (swap_send): use psb_snd (but need another
! call with swap_recv to complete)
! - if (swap_recv): use psb_rcv (completing a
! previous call with swap_send)
!
!
subroutine psb_i2ovrlm(x,desc_a,info,jx,ik,work,update,mode)
use psb_base_mod, psb_protect_name => psb_i2ovrlm
use psi_mod
implicit none
integer(psb_i2pk_), intent(inout), target :: x(:,:)
type(psb_desc_type), intent(in) :: desc_a
integer(psb_ipk_), intent(out) :: info
integer(psb_i2pk_), optional, target, intent(inout) :: work(:)
integer(psb_ipk_), intent(in), optional :: update,jx,ik,mode
! locals
integer(psb_mpk_) :: ictxt, np, me
integer(psb_ipk_) :: err_act, iix, jjx, nrow, ncol, k, maxk, update_,&
& mode_, liwork, ldx
integer(psb_lpk_) :: m, n, ix, ijx
integer(psb_i2pk_),pointer :: iwork(:), xp(:,:)
logical :: do_swap
character(len=20) :: name, ch_err
logical :: aliw
name='psb_i2ovrlm'
info=psb_success_
call psb_erractionsave(err_act)
if (psb_errstatus_fatal()) then
info = psb_err_internal_error_ ; goto 9999
end if
ictxt=desc_a%get_context()
! check on blacs grid
call psb_info(ictxt, me, np)
if (np == -1) then
info = psb_err_context_error_
call psb_errpush(info,name)
goto 9999
endif
ix = 1
if (present(jx)) then
ijx = jx
else
ijx = 1
endif
m = desc_a%get_global_rows()
n = desc_a%get_global_cols()
nrow = desc_a%get_local_rows()
ncol = desc_a%get_local_cols()
maxk=size(x,2)-ijx+1
if(present(ik)) then
if(ik > maxk) then
k=maxk
else
k=ik
end if
else
k = maxk
end if
if (present(update)) then
update_ = update
else
update_ = psb_avg_
endif
if (present(mode)) then
mode_ = mode
else
mode_ = IOR(psb_swap_send_,psb_swap_recv_)
endif
do_swap = (mode_ /= 0)
ldx = size(x,1)
! check vector correctness
call psb_chkvect(m,lone,ldx,ix,ijx,desc_a,info,iix,jjx)
if(info /= psb_success_) then
info=psb_err_from_subroutine_ ; ch_err='psb_chkvect'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
! check for presence/size of a work area
liwork=ncol
if (present(work)) then
if(size(work) >= liwork) then
aliw=.false.
else
aliw=.true.
end if
else
aliw=.true.
end if
if (aliw) then
allocate(iwork(liwork),stat=info)
if(info /= psb_success_) then
info=psb_err_from_subroutine_
call psb_errpush(info,name,a_err='Allocate')
goto 9999
end if
else
iwork => work
end if
! exchange overlap elements
if(do_swap) then
xp => x(iix:ldx,jjx:jjx+k-1)
call psi_swapdata(mode_,k,i2one,xp,&
& desc_a,iwork,info,data=psb_comm_ovr_)
end if
if (info == psb_success_) call psi_ovrl_upd(xp,desc_a,update_,info)
if (info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='Inner updates')
goto 9999
end if
if (aliw) deallocate(iwork)
nullify(iwork)
call psb_erractionrestore(err_act)
return
9999 call psb_error_handler(ione*ictxt,err_act)
return
end subroutine psb_i2ovrlm
!!$
!!$ Parallel Sparse BLAS version 3.5
!!$ (C) Copyright 2006-2018
!!$ Salvatore Filippone University of Rome Tor Vergata
!!$ 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.
!!$
!!$
!
! Subroutine: psb_i2ovrlv
! This subroutine performs the exchange of the overlap elements in a
! distributed dense vector between all the processes.
!
! Arguments:
! x(:) - integer The local part of the dense vector.
! desc_a - type(psb_desc_type). The communication descriptor.
! info - integer. Return code.
! work - integer(optional). A work area.
! update - integer(optional). Type of update:
! psb_none_ do nothing
! psb_sum_ sum of overlaps
! psb_avg_ average of overlaps
! mode - integer(optional). Choose the algorithm for data exchange:
! this is chosen through bit fields.
! - swap_mpi = iand(flag,psb_swap_mpi_) /= 0
! - swap_sync = iand(flag,psb_swap_sync_) /= 0
! - swap_send = iand(flag,psb_swap_send_) /= 0
! - swap_recv = iand(flag,psb_swap_recv_) /= 0
! - if (swap_mpi): use underlying MPI_ALLTOALLV.
! - if (swap_sync): use PSB_SND and PSB_RCV in
! synchronized pairs
! - if (swap_send .and. swap_recv): use mpi_irecv
! and mpi_send
! - if (swap_send): use psb_snd (but need another
! call with swap_recv to complete)
! - if (swap_recv): use psb_rcv (completing a
! previous call with swap_send)
!
!
subroutine psb_i2ovrlv(x,desc_a,info,work,update,mode)
use psb_base_mod, psb_protect_name => psb_i2ovrlv
use psi_mod
implicit none
integer(psb_i2pk_), intent(inout), target :: x(:)
type(psb_desc_type), intent(in) :: desc_a
integer(psb_ipk_), intent(out) :: info
integer(psb_i2pk_), optional, target, intent(inout) :: work(:)
integer(psb_ipk_), intent(in), optional :: update,mode
! locals
integer(psb_ipk_) :: ictxt, np, me, err_act, iix, jjx, nrow, ncol, &
& k, update_, mode_, liwork, ldx
integer(psb_lpk_) :: m, n, ix, ijx
integer(psb_i2pk_),pointer :: iwork(:)
logical :: do_swap
character(len=20) :: name, ch_err
logical :: aliw
name='psb_i2ovrlv'
info=psb_success_
call psb_erractionsave(err_act)
if (psb_errstatus_fatal()) then
info = psb_err_internal_error_ ; goto 9999
end if
ictxt=desc_a%get_context()
! check on blacs grid
call psb_info(ictxt, me, np)
if (np == -1) then
info = psb_err_context_error_
call psb_errpush(info,name)
goto 9999
endif
ix = 1
ijx = 1
m = desc_a%get_global_rows()
n = desc_a%get_global_cols()
nrow = desc_a%get_local_rows()
ncol = desc_a%get_local_cols()
k = 1
if (present(update)) then
update_ = update
else
update_ = psb_avg_
endif
if (present(mode)) then
mode_ = mode
else
mode_ = IOR(psb_swap_send_,psb_swap_recv_)
endif
do_swap = (mode_ /= 0)
ldx = size(x,1)
! check vector correctness
call psb_chkvect(m,lone,ldx,ix,ijx,desc_a,info,iix,jjx)
if(info /= psb_success_) then
info=psb_err_from_subroutine_
ch_err='psb_chkvect'
call psb_errpush(info,name,a_err=ch_err)
end if
if (iix /= 1) then
info=psb_err_ix_n1_iy_n1_unsupported_
call psb_errpush(info,name)
end if
if(info /= 0) goto 9999
! check for presence/size of a work area
liwork=ncol
if (present(work)) then
if(size(work) >= liwork) then
aliw=.false.
else
aliw=.true.
end if
else
aliw=.true.
end if
if (aliw) then
allocate(iwork(liwork),stat=info)
if(info /= psb_success_) then
info=psb_err_from_subroutine_
call psb_errpush(info,name,a_err='Allocate')
goto 9999
end if
else
iwork => work
end if
! exchange overlap elements
if (do_swap) then
call psi_swapdata(mode_,i2one,x,&
& desc_a,iwork,info,data=psb_comm_ovr_)
end if
if (info == psb_success_) call psi_ovrl_upd(x,desc_a,update_,info)
if (info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='Inner updates')
goto 9999
end if
if (aliw) deallocate(iwork)
nullify(iwork)
call psb_erractionrestore(err_act)
return
9999 call psb_error_handler(ione*ictxt,err_act)
return
end subroutine psb_i2ovrlv
+480
View File
@@ -0,0 +1,480 @@
!
! Parallel Sparse BLAS version 3.5
! (C) Copyright 2006-2018
! Salvatore Filippone
! Alfredo Buttari
!
! Redistribution and use in source and binary forms, with or without
! modification, are permitted provided that the following conditions
! are met:
! 1. Redistributions of source code must retain the above copyright
! notice, this list of conditions and the following disclaimer.
! 2. Redistributions in binary form must reproduce the above copyright
! notice, this list of conditions, and the following disclaimer in the
! documentation and/or other materials provided with the distribution.
! 3. The name of the PSBLAS group or the names of its contributors may
! not be used to endorse or promote products derived from this
! software without specific written permission.
!
! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
! ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
! TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
! PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE PSBLAS GROUP OR ITS CONTRIBUTORS
! BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
! CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
! SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
! INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
! CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
! ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
! POSSIBILITY OF SUCH DAMAGE.
!
!
! File: psb_i2scatter.f90
!
! Subroutine: psb_i2scatterm
! This subroutine scatters a global matrix locally owned by one process
! into pieces that are local to all the processes.
!
! Arguments:
! globx - integer,dimension(:,:). The global matrix to scatter.
! locx - integer,dimension(:,:). The local piece of the distributed matrix.
! desc_a - type(psb_desc_type). The communication descriptor.
! info - integer. Error code.
! iroot - integer(optional). The process that owns the global matrix.
! If -1 all the processes have a copy.
! Default -1
subroutine psb_i2scatterm(globx, locx, desc_a, info, root)
use psb_base_mod, psb_protect_name => psb_i2scatterm
#ifdef MPI_MOD
use mpi
#endif
implicit none
#ifdef MPI_H
include 'mpif.h'
#endif
integer(psb_i2pk_), intent(out), allocatable :: locx(:,:)
integer(psb_i2pk_), intent(in) :: globx(:,:)
type(psb_desc_type), intent(in) :: desc_a
integer(psb_ipk_), intent(out) :: info
integer(psb_ipk_), intent(in), optional :: root
! locals
integer(psb_mpk_) :: ictxt, np, me, iroot, icomm, myrank, rootrank, iam, nlr
integer(psb_ipk_) :: ierr(5), err_act, nrow,&
& ilocx, jlocx, lda_locx, lda_globx, lock, globk, k, maxk, &
& col,pos
integer(psb_lpk_) :: m, n, i, j, idx, iglobx, jglobx
integer(psb_i2pk_),allocatable :: scatterv(:)
integer(psb_mpk_), allocatable :: displ(:), all_dim(:)
integer(psb_lpk_), allocatable :: l_t_g_all(:), ltg(:)
character(len=20) :: name, ch_err
name='psb_scatterm'
info=psb_success_
call psb_erractionsave(err_act)
if (psb_errstatus_fatal()) then
info = psb_err_internal_error_ ; goto 9999
end if
ictxt=desc_a%get_context()
! check on blacs grid
call psb_info(ictxt, iam, np)
if (np == -1) then
info = psb_err_context_error_
call psb_errpush(info,name)
goto 9999
endif
if (present(root)) then
iroot = root
if((iroot < -1).or.(iroot >= np)) then
info=psb_err_input_value_invalid_i_
ierr(1)=5; ierr(2)=iroot
call psb_errpush(info,name,i_err=ierr)
goto 9999
end if
else
iroot = psb_root_
end if
iglobx = 1
jglobx = 1
lda_globx = size(globx,1)
m = desc_a%get_global_rows()
n = desc_a%get_global_cols()
icomm = psb_get_mpi_comm(ictxt)
myrank = psb_get_mpi_rank(ictxt,me)
if (iroot==-1) then
lda_globx = size(globx, 1)
k = size(globx,2)
else
if (iam==iroot) then
k = size(globx,2)
lda_globx = size(globx, 1)
end if
end if
m = desc_a%get_global_rows()
n = desc_a%get_global_cols()
! there should be a global check on k here!!!
if ((iroot==-1).or.(iam==iroot)) &
& call psb_chkglobvect(m,n,lda_globx,iglobx,jglobx,desc_a,info)
if(info /= psb_success_) then
info=psb_err_from_subroutine_
ch_err='psb_chk(glob)vect'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
nrow=desc_a%get_local_rows()
! root has to gather size information
allocate(displ(np),all_dim(np),ltg(nrow),stat=info)
if(info /= psb_success_) then
info=psb_err_from_subroutine_
ch_err='Allocate'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
do i=1, nrow
ltg(i) = i
end do
call psb_loc_to_glob(ltg(1:nrow),desc_a,info)
call psb_geall(locx,desc_a,info,n=k)
if ((iroot == -1).or.(np == 1)) then
! extract my chunk
do j=1,k
do i=1, nrow
locx(i,j)=globx(ltg(i),j)
end do
end do
else
rootrank = psb_get_mpi_rank(ictxt,iroot)
!
! This is potentially unsafe when IPK=8
! But then, IPK=8 is highly experimental anyway.
!
nlr = nrow
call mpi_gather(nlr,1,psb_mpi_mpk_,all_dim,&
& 1,psb_mpi_mpk_,rootrank,icomm,info)
if (iam == iroot) then
displ(1)=0
do i=2,np
displ(i)=displ(i-1)+all_dim(i-1)
end do
! root has to gather loc_glob from each process
allocate(l_t_g_all(sum(all_dim)),scatterv(sum(all_dim)),stat=info)
else
!
! This is to keep debugging compilers from being upset by
! calling an external MPI function with an unallocated array;
! the Fortran side would complain even if the MPI side does
! not use the unallocated stuff.
!
allocate(l_t_g_all(1),scatterv(1),stat=info)
end if
if(info /= psb_success_) then
info=psb_err_from_subroutine_
ch_err='Allocate'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
call mpi_gatherv(ltg,nlr,&
& psb_mpi_lpk_,l_t_g_all,all_dim,&
& displ,psb_mpi_lpk_,rootrank,icomm,info)
do col=1, k
! prepare vector to scatter
if(iam == iroot) then
do i=1,np
pos=displ(i)
do j=1, all_dim(i)
idx=l_t_g_all(pos+j)
scatterv(pos+j)=globx(idx,col)
end do
end do
end if
! scatter
call mpi_scatterv(scatterv,all_dim,displ,&
& psb_mpi_i2pk_,locx(1,col),nrow,&
& psb_mpi_i2pk_,rootrank,icomm,info)
end do
deallocate(l_t_g_all, scatterv,stat=info)
if(info /= psb_success_) then
info=psb_err_from_subroutine_
ch_err='deallocate'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
end if
deallocate(all_dim, displ, ltg,stat=info)
if(info /= psb_success_) then
info=psb_err_from_subroutine_
ch_err='deallocate'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
call psb_erractionrestore(err_act)
return
9999 call psb_error_handler(ione*ictxt,err_act)
return
end subroutine psb_i2scatterm
!!$
!!$ Parallel Sparse BLAS version 3.5
!!$ (C) Copyright 2006-2018
!!$ Salvatore Filippone University of Rome Tor Vergata
!!$ 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.
!!$
!!$
! Subroutine: psb_i2scatterv
! This subroutine scatters a global vector locally owned by one process
! into pieces that are local to all the processes.
!
! Arguments:
! globx - integer,dimension(:). The global vector to scatter.
! locx - integer,dimension(:). The local piece of the ditributed vector.
! desc_a - type(psb_desc_type). The communication descriptor.
! info - integer. Return code
! iroot - integer(optional). The process that owns the global vector. If -1 all
! the processes have a copy.
!
subroutine psb_i2scatterv(globx, locx, desc_a, info, root)
use psb_base_mod, psb_protect_name => psb_i2scatterv
#ifdef MPI_MOD
use mpi
#endif
implicit none
#ifdef MPI_H
include 'mpif.h'
#endif
integer(psb_i2pk_), intent(out), allocatable :: locx(:)
integer(psb_i2pk_), intent(in) :: globx(:)
type(psb_desc_type), intent(in) :: desc_a
integer(psb_ipk_), intent(out) :: info
integer(psb_ipk_), intent(in), optional :: root
! locals
integer(psb_mpk_) :: ictxt, np, iam, iroot, iiroot, icomm, myrank, rootrank, nlr
integer(psb_ipk_) :: ierr(5), err_act, nrow,&
& ilocx, jlocx, lda_locx, lda_globx, k, pos, ilx, jlx
integer(psb_lpk_) :: m, n, i, j, idx, iglobx, jglobx
integer(psb_i2pk_), allocatable :: scatterv(:)
integer(psb_mpk_), allocatable :: displ(:), all_dim(:)
integer(psb_lpk_), allocatable :: l_t_g_all(:), ltg(:)
character(len=20) :: name, ch_err
integer(psb_ipk_) :: debug_level, debug_unit
name='psb_scatterv'
info=psb_success_
call psb_erractionsave(err_act)
if (psb_errstatus_fatal()) then
info = psb_err_internal_error_ ; goto 9999
end if
ictxt=desc_a%get_context()
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
! check on blacs grid
call psb_info(ictxt, iam, np)
if (np == -1) then
info = psb_err_context_error_
call psb_errpush(info,name)
goto 9999
endif
if (present(root)) then
iroot = root
if((iroot < -1).or.(iroot > np)) then
info=psb_err_input_value_invalid_i_
ierr(1) = 5; ierr(2)=iroot
call psb_errpush(info,name,i_err=ierr)
goto 9999
end if
else
iroot = psb_root_
end if
icomm = psb_get_mpi_comm(ictxt)
myrank = psb_get_mpi_rank(ictxt,iam)
iglobx = 1
jglobx = 1
ilocx = 1
jlocx = 1
if ((iroot==-1).or.(iam==iroot))&
& lda_globx = size(globx, 1)
m = desc_a%get_global_rows()
n = desc_a%get_global_cols()
k = 1
! there should be a global check on k here!!!
if ((iroot==-1).or.(iam==iroot)) &
& call psb_chkglobvect(m,n,lda_globx,iglobx,jglobx,desc_a,info)
if(info /= psb_success_) then
info=psb_err_from_subroutine_
ch_err='psb_chk(glob)vect'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
nrow = desc_a%get_local_rows()
allocate(displ(np),all_dim(np),ltg(nrow),stat=info)
if(info /= psb_success_) then
info=psb_err_from_subroutine_
ch_err='Allocate'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
do i=1, nrow
ltg(i) = i
end do
call psb_loc_to_glob(ltg(1:nrow),desc_a,info)
call psb_geall(locx,desc_a,info)
if ((iroot == -1).or.(np == 1)) then
! extract my chunk
do i=1, nrow
locx(i)=globx(ltg(i))
end do
else
rootrank = psb_get_mpi_rank(ictxt,iroot)
!
! This is potentially unsafe when IPK=8
! But then, IPK=8 is highly experimental anyway.
!
nlr = nrow
call mpi_gather(nlr,1,psb_mpi_mpk_,all_dim,&
& 1,psb_mpi_mpk_,rootrank,icomm,info)
if(iam == iroot) then
displ(1)=0
do i=2,np
displ(i)=displ(i-1) + all_dim(i-1)
end do
if (debug_level >= psb_debug_inner_) then
write(debug_unit,*) iam,' ',trim(name),' displ:',displ(1:np), &
&' dim',all_dim(1:np), sum(all_dim)
endif
! root has to gather loc_glob from each process
allocate(l_t_g_all(sum(all_dim)),scatterv(sum(all_dim)),stat=info)
else
!
! This is to keep debugging compilers from being upset by
! calling an external MPI function with an unallocated array;
! the Fortran side would complain even if the MPI side does
! not use the unallocated stuff.
!
allocate(l_t_g_all(1),scatterv(1),stat=info)
end if
if(info /= psb_success_) then
info=psb_err_from_subroutine_
ch_err='Allocate'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
call mpi_gatherv(ltg,nlr,&
& psb_mpi_lpk_,l_t_g_all,all_dim,&
& displ,psb_mpi_lpk_,rootrank,icomm,info)
! prepare vector to scatter
if (iam == iroot) then
do i=1,np
pos=displ(i)
do j=1, all_dim(i)
idx=l_t_g_all(pos+j)
scatterv(pos+j)=globx(idx)
end do
end do
end if
call mpi_scatterv(scatterv,all_dim,displ,&
& psb_mpi_i2pk_,locx,nrow,&
& psb_mpi_i2pk_,rootrank,icomm,info)
deallocate(l_t_g_all, scatterv,stat=info)
if(info /= psb_success_) then
info=psb_err_from_subroutine_
ch_err='deallocate'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
end if
deallocate(all_dim, displ, ltg,stat=info)
if(info /= psb_success_) then
info=psb_err_from_subroutine_
ch_err='deallocate'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
call psb_erractionrestore(err_act)
return
9999 call psb_error_handler(ione*ictxt,err_act)
return
end subroutine psb_i2scatterv