From 9dc5a8bc11f2c58799ba2dab82b130d50184e403 Mon Sep 17 00:00:00 2001 From: Salvatore Filippone Date: Fri, 8 Feb 2008 12:35:30 +0000 Subject: [PATCH] psblas: base/internals/Makefile base/internals/psi_dswapdata.F90 base/internals/psi_dswaptran.F90 base/internals/psi_iswapdata.F90 base/internals/psi_iswaptran.F90 base/internals/psi_zswapdata.F90 base/internals/psi_zswaptran.F90 base/modules/Makefile base/modules/psb_base_mod.f90 base/modules/psb_desc_type.f90 base/modules/psb_inter_desc_type.f90 base/modules/psb_realloc_mod.F90 base/modules/psb_serial_mod.f90 base/modules/psb_spmat_type.f90 base/modules/psb_tools_mod.f90 base/modules/psi_mod.f90 base/modules/psi_serial_mod.f90 base/tools/psb_cdcpy.f90 base/tools/psb_dcdovr.F90 base/tools/psb_glob_to_loc.f90 base/tools/psb_loc_to_glob.f90 base/tools/psb_zcdovr.F90 Merged changes from psblas-intermesh branch up to rev. 2809. --- base/internals/Makefile | 2 +- base/internals/psi_dswapdata.F90 | 400 +++++++++------ base/internals/psi_dswaptran.F90 | 410 +++++++++------ base/internals/psi_iswapdata.F90 | 404 +++++++++------ base/internals/psi_iswaptran.F90 | 407 +++++++++------ base/internals/psi_zswapdata.F90 | 419 ++++++++++------ base/internals/psi_zswaptran.F90 | 409 +++++++++------ base/modules/Makefile | 15 +- base/modules/psb_base_mod.f90 | 2 + base/modules/psb_desc_type.f90 | 147 +++++- base/modules/psb_inter_desc_type.f90 | 655 ++++++++++++++++++++++++ base/modules/psb_realloc_mod.F90 | 323 +++++++++++- base/modules/psb_serial_mod.f90 | 5 + base/modules/psb_spmat_type.f90 | 20 +- base/modules/psb_tools_mod.f90 | 84 +++- base/modules/psi_mod.f90 | 726 +++++++++++++++------------ base/modules/psi_serial_mod.f90 | 518 +++++++++++++++++++ base/tools/psb_cdcpy.f90 | 26 +- base/tools/psb_dcdovr.F90 | 34 +- base/tools/psb_glob_to_loc.f90 | 35 ++ base/tools/psb_loc_to_glob.f90 | 33 ++ base/tools/psb_zcdovr.F90 | 26 +- 22 files changed, 3692 insertions(+), 1408 deletions(-) create mode 100644 base/modules/psb_inter_desc_type.f90 create mode 100644 base/modules/psi_serial_mod.f90 diff --git a/base/internals/Makefile b/base/internals/Makefile index ecad6843..c43b47ed 100644 --- a/base/internals/Makefile +++ b/base/internals/Makefile @@ -21,7 +21,7 @@ lib: mpfobjs $(FOBJS) $(FOBJS2) $(COBJS) $(MPFOBJS2) $(COBJS) $(RANLIB) $(LIBDIR)/$(LIBNAME) - +$(FOBJS) $(FBOJS2): $(MODDIR)/psi_mod$(.mod) mpfobjs: (make $(MPFOBJS) F90="$(MPF90)" FC="$(MPF90)" FCOPT="$(F90COPT)") (make $(FOBJS2) F90="$(MPF77)" FC="$(MPF77)" FCOPT="$(FCOPT)") diff --git a/base/internals/psi_dswapdata.F90 b/base/internals/psi_dswapdata.F90 index 78568524..8930bd28 100644 --- a/base/internals/psi_dswapdata.F90 +++ b/base/internals/psi_dswapdata.F90 @@ -40,7 +40,7 @@ ! D real(kind(1.d0)) ! Z complex(kind(1.d0)) ! 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))); +! 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 @@ -87,7 +87,6 @@ subroutine psi_dswapdatam(flag,n,beta,y,desc_a,work,info,data) use psb_error_mod use psb_descriptor_type use psb_penv_mod -!!$ use psi_gthsct_mod #ifdef MPI_MOD use mpi #endif @@ -104,19 +103,8 @@ subroutine psi_dswapdatam(flag,n,beta,y,desc_a,work,info,data) integer, optional :: data ! locals - integer :: ictxt, np, me, nesd, nerv,& - & proc_to_comm, p2ptag, icomm, p2pstat(mpi_status_size),& - & idxs, idxr, iret, err_act, totxch, i, idx_pt,& - & snd_pt, rcv_pt, pnti, data_ - integer, allocatable, dimension(:) :: bsdidx, brvidx,& - & sdsz, rvsz, prcid, rvhd, sdhd + integer :: ictxt, np, me, icomm, idxs, idxr, totxch, data_, err_act integer, pointer :: d_idx(:) - integer :: int_err(5) - logical :: swap_mpi, swap_sync, swap_send, swap_recv,& - & albf,do_send,do_recv - logical, parameter :: usersend=.false. - - real(kind(1.d0)), pointer, dimension(:) :: sndbuf, rcvbuf character(len=20) :: name info = 0 @@ -138,14 +126,6 @@ subroutine psi_dswapdatam(flag,n,beta,y,desc_a,work,info,data) 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 - if(present(data)) then data_ = data else @@ -158,8 +138,78 @@ subroutine psi_dswapdatam(flag,n,beta,y,desc_a,work,info,data) goto 9999 end if - idxr = idxr * n - idxs = idxs * n + + call psi_swapdata(ictxt,icomm,flag,n,beta,y,d_idx,totxch,idxs,idxr,work,info) + if (info /= 0) goto 9999 + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + if (err_act == psb_act_abort_) then + call psb_error(ictxt) + return + end if + return +end subroutine psi_dswapdatam + +subroutine psi_dswapidxm(ictxt,icomm,flag,n,beta,y,idx,totxch,totsnd,totrcv,work,info) + + use psi_mod, psb_protect_name => psi_dswapidxm + use psb_error_mod + use psb_descriptor_type + use psb_penv_mod +#ifdef MPI_MOD + use mpi +#endif + implicit none +#ifdef MPI_H + include 'mpif.h' +#endif + + integer, intent(in) :: ictxt,icomm,flag,n + integer, intent(out) :: info + real(kind(1.d0)) :: y(:,:), beta + real(kind(1.d0)), target :: work(:) + integer, intent(in) :: idx(:),totxch,totsnd, totrcv + + ! locals + integer :: np, me, nesd, nerv,& + & proc_to_comm, p2ptag, p2pstat(mpi_status_size),& + & iret, err_act, i, idx_pt, totsnd_, totrcv_,& + & snd_pt, rcv_pt, pnti, data_ + integer, allocatable, dimension(:) :: bsdidx, brvidx,& + & sdsz, rvsz, prcid, rvhd, sdhd + integer :: int_err(5) + logical :: swap_mpi, swap_sync, swap_send, swap_recv,& + & albf,do_send,do_recv + logical, parameter :: usersend=.false. + + real(kind(1.d0)), pointer, dimension(:) :: sndbuf, rcvbuf + character(len=20) :: name + + info = 0 + name='psi_swap_data' + call psb_erractionsave(err_act) + + call psb_info(ictxt,me,np) + if (np == -1) then + info = 2010 + 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),& @@ -181,9 +231,9 @@ subroutine psi_dswapdatam(flag,n,beta,y,desc_a,work,info,data) snd_pt = 1 rcv_pt = 1 do i=1, totxch - proc_to_comm = d_idx(pnti+psb_proc_id_) - nerv = d_idx(pnti+psb_n_elem_recv_) - nesd = d_idx(pnti+nerv+psb_n_elem_send_) + proc_to_comm = idx(pnti+psb_proc_id_) + nerv = idx(pnti+psb_n_elem_recv_) + nesd = idx(pnti+nerv+psb_n_elem_send_) call psb_get_rank(prcid(proc_to_comm),ictxt,proc_to_comm) brvidx(proc_to_comm) = rcv_pt @@ -206,14 +256,14 @@ subroutine psi_dswapdatam(flag,n,beta,y,desc_a,work,info,data) end if end if - idxr = max(idxr,1) - idxs = max(idxs,1) - if((idxr+idxs) < size(work)) then - sndbuf => work(1:idxs) - rcvbuf => work(idxs+1:idxs+idxr) + 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(idxs),rcvbuf(idxr), stat=info) + allocate(sndbuf(totsnd_),rcvbuf(totrcv_), stat=info) if(info /= 0) then call psb_errpush(4000,name) goto 9999 @@ -227,10 +277,10 @@ subroutine psi_dswapdatam(flag,n,beta,y,desc_a,work,info,data) pnti = 1 snd_pt = 1 do i=1, totxch - nerv = d_idx(pnti+psb_n_elem_recv_) - nesd = d_idx(pnti+nerv+psb_n_elem_send_) + 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,d_idx(idx_pt:idx_pt+nesd-1),& + 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 @@ -261,9 +311,9 @@ subroutine psi_dswapdatam(flag,n,beta,y,desc_a,work,info,data) snd_pt = 1 rcv_pt = 1 do i=1, totxch - proc_to_comm = d_idx(pnti+psb_proc_id_) - nerv = d_idx(pnti+psb_n_elem_recv_) - nesd = d_idx(pnti+nerv+psb_n_elem_send_) + 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,& @@ -275,6 +325,11 @@ subroutine psi_dswapdatam(flag,n,beta,y,desc_a,work,info,data) & 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(0,*) 'Fatal error in swapdata: mismatch on self sendf',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 @@ -292,11 +347,11 @@ subroutine psi_dswapdatam(flag,n,beta,y,desc_a,work,info,data) snd_pt = 1 rcv_pt = 1 do i=1, totxch - proc_to_comm = d_idx(pnti+psb_proc_id_) - nerv = d_idx(pnti+psb_n_elem_recv_) - nesd = d_idx(pnti+nerv+psb_n_elem_send_) + proc_to_comm = idx(pnti+psb_proc_id_) + nerv = idx(pnti+psb_n_elem_recv_) + nesd = idx(pnti+nerv+psb_n_elem_send_) call psb_get_rank(prcid(i),ictxt,proc_to_comm) - if (nerv>0) then + if ((nerv>0).and.(proc_to_comm/=me)) then p2ptag = krecvid(ictxt,proc_to_comm,me) call mpi_irecv(rcvbuf(rcv_pt),n*nerv,& & mpi_double_precision,prcid(i),& @@ -316,12 +371,12 @@ subroutine psi_dswapdatam(flag,n,beta,y,desc_a,work,info,data) snd_pt = 1 rcv_pt = 1 do i=1, totxch - proc_to_comm = d_idx(pnti+psb_proc_id_) - nerv = d_idx(pnti+psb_n_elem_recv_) - nesd = d_idx(pnti+nerv+psb_n_elem_send_) + proc_to_comm = idx(pnti+psb_proc_id_) + nerv = idx(pnti+psb_n_elem_recv_) + nesd = idx(pnti+nerv+psb_n_elem_send_) p2ptag=ksendid(ictxt,proc_to_comm,me) - if (nesd>0) then + if ((nesd>0).and.(proc_to_comm/=me)) then if (usersend) then call mpi_rsend(sndbuf(snd_pt),n*nesd,& & mpi_double_precision,prcid(i),& @@ -349,9 +404,9 @@ subroutine psi_dswapdatam(flag,n,beta,y,desc_a,work,info,data) pnti = 1 do i=1, totxch - proc_to_comm = d_idx(pnti+psb_proc_id_) - nerv = d_idx(pnti+psb_n_elem_recv_) - nesd = d_idx(pnti+nerv+psb_n_elem_send_) + proc_to_comm = idx(pnti+psb_proc_id_) + nerv = idx(pnti+psb_n_elem_recv_) + nesd = idx(pnti+nerv+psb_n_elem_send_) p2ptag = krecvid(ictxt,proc_to_comm,me) @@ -363,6 +418,11 @@ subroutine psi_dswapdatam(flag,n,beta,y,desc_a,work,info,data) call psb_errpush(info,name,i_err=int_err) goto 9999 end if + else if (proc_to_comm == me) then + if (nesd /= nerv) then + write(0,*) 'Fatal error in swapdata: mismatch on self sendf',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 @@ -375,9 +435,9 @@ subroutine psi_dswapdatam(flag,n,beta,y,desc_a,work,info,data) snd_pt = 1 rcv_pt = 1 do i=1, totxch - proc_to_comm = d_idx(pnti+psb_proc_id_) - nerv = d_idx(pnti+psb_n_elem_recv_) - nesd = d_idx(pnti+nerv+psb_n_elem_send_) + 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) @@ -394,9 +454,9 @@ subroutine psi_dswapdatam(flag,n,beta,y,desc_a,work,info,data) snd_pt = 1 rcv_pt = 1 do i=1, totxch - proc_to_comm = d_idx(pnti+psb_proc_id_) - nerv = d_idx(pnti+psb_n_elem_recv_) - nesd = d_idx(pnti+nerv+psb_n_elem_send_) + 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 @@ -415,11 +475,11 @@ subroutine psi_dswapdatam(flag,n,beta,y,desc_a,work,info,data) snd_pt = 1 rcv_pt = 1 do i=1, totxch - proc_to_comm = d_idx(pnti+psb_proc_id_) - nerv = d_idx(pnti+psb_n_elem_recv_) - nesd = d_idx(pnti+nerv+psb_n_elem_send_) + 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,d_idx(idx_pt:idx_pt+nerv-1),& + 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 @@ -455,39 +515,8 @@ subroutine psi_dswapdatam(flag,n,beta,y,desc_a,work,info,data) return end if return -end subroutine psi_dswapdatam +end subroutine psi_dswapidxm - -!!$ -!!$ Parallel Sparse BLAS v2.0 -!!$ (C) Copyright 2006 Salvatore Filippone University of Rome Tor Vergata -!!$ Alfredo Buttari University of Rome Tor Vergata -!!$ -!!$ 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_dswapdatav @@ -544,7 +573,6 @@ subroutine psi_dswapdatav(flag,beta,y,desc_a,work,info,data) use psb_error_mod use psb_descriptor_type use psb_penv_mod -!!$ use psi_gthsct_mod #ifdef MPI_MOD use mpi #endif @@ -561,20 +589,8 @@ subroutine psi_dswapdatav(flag,beta,y,desc_a,work,info,data) integer, optional :: data ! locals - integer :: ictxt, np, me, nesd, nerv,& - & proc_to_comm, p2ptag, icomm, p2pstat(mpi_status_size),& - & idxs, idxr, iret, err_act, totxch, i, & - & idx_pt, snd_pt, rcv_pt, n, pnti, data_ - - integer, allocatable, dimension(:) :: bsdidx, brvidx,& - & sdsz, rvsz, prcid, rvhd, sdhd + integer :: ictxt, np, me, icomm, idxs, idxr, totxch, data_, err_act integer, pointer :: d_idx(:) - integer :: int_err(5) - logical :: swap_mpi, swap_sync, swap_send, swap_recv,& - & albf,do_send,do_recv - logical, parameter :: usersend=.false. - - real(kind(1.d0)), pointer, dimension(:) :: sndbuf, rcvbuf character(len=20) :: name info = 0 @@ -595,16 +611,7 @@ subroutine psi_dswapdatav(flag,beta,y,desc_a,work,info,data) goto 9999 endif - icomm = desc_a%matrix_data(psb_mpi_c_) - 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 + icomm = psb_cd_get_mpic(desc_a) if(present(data)) then data_ = data @@ -618,8 +625,80 @@ subroutine psi_dswapdatav(flag,beta,y,desc_a,work,info,data) goto 9999 end if - idxr = idxr * n - idxs = idxs * n + call psi_swapdata(ictxt,icomm,flag,beta,y,d_idx,totxch,idxs,idxr,work,info) + if (info /= 0) goto 9999 + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + if (err_act == psb_act_abort_) then + call psb_error(ictxt) + return + end if + return +end subroutine psi_dswapdatav + + +subroutine psi_dswapidxv(ictxt,icomm,flag,beta,y,idx,totxch,totsnd,totrcv,work,info) + + use psi_mod, psb_protect_name => psi_dswapidxv + use psb_error_mod + use psb_descriptor_type + use psb_penv_mod +#ifdef MPI_MOD + use mpi +#endif + implicit none +#ifdef MPI_H + include 'mpif.h' +#endif + + integer, intent(in) :: ictxt,icomm,flag + integer, intent(out) :: info + real(kind(1.d0)) :: y(:), beta + real(kind(1.d0)), target :: work(:) + integer, intent(in) :: idx(:),totxch,totsnd, totrcv + + ! locals + integer :: np, me, nesd, nerv,& + & proc_to_comm, p2ptag, p2pstat(mpi_status_size),& + & iret, err_act, i, totsnd_, totrcv_,& + & idx_pt, snd_pt, rcv_pt, n, pnti, data_ + + integer, allocatable, dimension(:) :: bsdidx, brvidx,& + & sdsz, rvsz, prcid, rvhd, sdhd + integer :: int_err(5) + logical :: swap_mpi, swap_sync, swap_send, swap_recv,& + & albf,do_send,do_recv + logical, parameter :: usersend=.false. + + real(kind(1.d0)), pointer, dimension(:) :: sndbuf, rcvbuf + character(len=20) :: name + + info = 0 + name='psi_swap_datav' + call psb_erractionsave(err_act) + + call psb_info(ictxt,me,np) + if (np == -1) then + info = 2010 + 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),& @@ -640,9 +719,9 @@ subroutine psi_dswapdatav(flag,beta,y,desc_a,work,info,data) snd_pt = 1 rcv_pt = 1 do i=1, totxch - proc_to_comm = d_idx(pnti+psb_proc_id_) - nerv = d_idx(pnti+psb_n_elem_recv_) - nesd = d_idx(pnti+nerv+psb_n_elem_send_) + proc_to_comm = idx(pnti+psb_proc_id_) + nerv = idx(pnti+psb_n_elem_recv_) + nesd = idx(pnti+nerv+psb_n_elem_send_) call psb_get_rank(prcid(proc_to_comm),ictxt,proc_to_comm) brvidx(proc_to_comm) = rcv_pt @@ -666,14 +745,14 @@ subroutine psi_dswapdatav(flag,beta,y,desc_a,work,info,data) end if - idxr = max(idxr,1) - idxs = max(idxs,1) - if((idxr+idxs) < size(work)) then - sndbuf => work(1:idxs) - rcvbuf => work(idxs+1:idxs+idxr) + 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(idxs),rcvbuf(idxr), stat=info) + allocate(sndbuf(totsnd_),rcvbuf(totrcv_), stat=info) if(info /= 0) then call psb_errpush(4000,name) goto 9999 @@ -688,10 +767,10 @@ subroutine psi_dswapdatav(flag,beta,y,desc_a,work,info,data) pnti = 1 snd_pt = 1 do i=1, totxch - nerv = d_idx(pnti+psb_n_elem_recv_) - nesd = d_idx(pnti+nerv+psb_n_elem_send_) + 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,d_idx(idx_pt:idx_pt+nesd-1),& + 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 @@ -719,9 +798,9 @@ subroutine psi_dswapdatav(flag,beta,y,desc_a,work,info,data) snd_pt = 1 rcv_pt = 1 do i=1, totxch - proc_to_comm = d_idx(pnti+psb_proc_id_) - nerv = d_idx(pnti+psb_n_elem_recv_) - nesd = d_idx(pnti+nerv+psb_n_elem_send_) + 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,& @@ -733,6 +812,11 @@ subroutine psi_dswapdatav(flag,beta,y,desc_a,work,info,data) & 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(0,*) 'Fatal error in swapdata: mismatch on self sendf',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 @@ -747,12 +831,12 @@ subroutine psi_dswapdatav(flag,beta,y,desc_a,work,info,data) snd_pt = 1 rcv_pt = 1 do i=1, totxch - proc_to_comm = d_idx(pnti+psb_proc_id_) - nerv = d_idx(pnti+psb_n_elem_recv_) - nesd = d_idx(pnti+nerv+psb_n_elem_send_) + proc_to_comm = idx(pnti+psb_proc_id_) + nerv = idx(pnti+psb_n_elem_recv_) + nesd = idx(pnti+nerv+psb_n_elem_send_) call psb_get_rank(prcid(i),ictxt,proc_to_comm) - if (nerv>0) then + if ((nerv>0).and.(proc_to_comm /= me)) then p2ptag = krecvid(ictxt,proc_to_comm,me) call mpi_irecv(rcvbuf(rcv_pt),nerv,& & mpi_double_precision,prcid(i),& @@ -771,13 +855,13 @@ subroutine psi_dswapdatav(flag,beta,y,desc_a,work,info,data) snd_pt = 1 rcv_pt = 1 do i=1, totxch - proc_to_comm = d_idx(pnti+psb_proc_id_) - nerv = d_idx(pnti+psb_n_elem_recv_) - nesd = d_idx(pnti+nerv+psb_n_elem_send_) + proc_to_comm = idx(pnti+psb_proc_id_) + nerv = idx(pnti+psb_n_elem_recv_) + nesd = idx(pnti+nerv+psb_n_elem_send_) p2ptag=ksendid(ictxt,proc_to_comm,me) - if (nesd>0) then + if ((nesd>0).and.(proc_to_comm /= me)) then if (usersend) then call mpi_rsend(sndbuf(snd_pt),nesd,& & mpi_double_precision,prcid(i),& @@ -803,9 +887,9 @@ subroutine psi_dswapdatav(flag,beta,y,desc_a,work,info,data) pnti = 1 do i=1, totxch - proc_to_comm = d_idx(pnti+psb_proc_id_) - nerv = d_idx(pnti+psb_n_elem_recv_) - nesd = d_idx(pnti+nerv+psb_n_elem_send_) + proc_to_comm = idx(pnti+psb_proc_id_) + nerv = idx(pnti+psb_n_elem_recv_) + nesd = idx(pnti+nerv+psb_n_elem_send_) p2ptag = krecvid(ictxt,proc_to_comm,me) @@ -817,6 +901,11 @@ subroutine psi_dswapdatav(flag,beta,y,desc_a,work,info,data) call psb_errpush(info,name,i_err=int_err) goto 9999 end if + else if (proc_to_comm == me) then + if (nesd /= nerv) then + write(0,*) 'Fatal error in swapdata: mismatch on self sendf',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 @@ -828,9 +917,9 @@ subroutine psi_dswapdatav(flag,beta,y,desc_a,work,info,data) snd_pt = 1 rcv_pt = 1 do i=1, totxch - proc_to_comm = d_idx(pnti+psb_proc_id_) - nerv = d_idx(pnti+psb_n_elem_recv_) - nesd = d_idx(pnti+nerv+psb_n_elem_send_) + 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 @@ -844,9 +933,9 @@ subroutine psi_dswapdatav(flag,beta,y,desc_a,work,info,data) snd_pt = 1 rcv_pt = 1 do i=1, totxch - proc_to_comm = d_idx(pnti+psb_proc_id_) - nerv = d_idx(pnti+psb_n_elem_recv_) - nesd = d_idx(pnti+nerv+psb_n_elem_send_) + 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 @@ -856,18 +945,17 @@ subroutine psi_dswapdatav(flag,beta,y,desc_a,work,info,data) end if - if (do_recv) then pnti = 1 snd_pt = 1 rcv_pt = 1 do i=1, totxch - proc_to_comm = d_idx(pnti+psb_proc_id_) - nerv = d_idx(pnti+psb_n_elem_recv_) - nesd = d_idx(pnti+nerv+psb_n_elem_send_) + 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,d_idx(idx_pt:idx_pt+nerv-1),& + 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 @@ -902,4 +990,4 @@ subroutine psi_dswapdatav(flag,beta,y,desc_a,work,info,data) return end if return -end subroutine psi_dswapdatav +end subroutine psi_dswapidxv diff --git a/base/internals/psi_dswaptran.F90 b/base/internals/psi_dswaptran.F90 index 9fb1a525..48a4ae54 100644 --- a/base/internals/psi_dswaptran.F90 +++ b/base/internals/psi_dswaptran.F90 @@ -90,7 +90,6 @@ subroutine psi_dswaptranm(flag,n,beta,y,desc_a,work,info,data) use psb_error_mod use psb_descriptor_type use psb_penv_mod -!!$ use psi_gthsct_mod #ifdef MPI_MOD use mpi #endif @@ -107,19 +106,9 @@ subroutine psi_dswaptranm(flag,n,beta,y,desc_a,work,info,data) integer, optional :: data ! locals - integer :: ictxt, np, me, nesd, nerv,& - & proc_to_comm, p2ptag, icomm, p2pstat(mpi_status_size),& - & idxs, idxr, iret, err_act, totxch, i, idx_pt,& - & snd_pt, rcv_pt, pnti, data_ - integer, allocatable, dimension(:) :: bsdidx, brvidx,& - & sdsz, rvsz, prcid, rvhd, sdhd + integer :: ictxt, np, me, icomm, idxs, idxr, err_act, totxch, data_ integer, pointer :: d_idx(:) integer :: int_err(5) - logical :: swap_mpi, swap_sync, swap_send, swap_recv,& - & albf,do_send,do_recv - logical, parameter :: usersend=.false. - - real(kind(1.d0)), pointer, dimension(:) :: sndbuf, rcvbuf character(len=20) :: name info = 0 @@ -142,13 +131,6 @@ subroutine psi_dswaptranm(flag,n,beta,y,desc_a,work,info,data) 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 - if(present(data)) then data_ = data else @@ -161,8 +143,77 @@ subroutine psi_dswaptranm(flag,n,beta,y,desc_a,work,info,data) goto 9999 end if - idxr = idxr * n - idxs = idxs * n + call psi_swaptran(ictxt,icomm,flag,n,beta,y,d_idx,totxch,idxs,idxr,work,info) + if (info /= 0) goto 9999 + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + if (err_act == psb_act_abort_) then + call psb_error(ictxt) + return + end if + return +end subroutine psi_dswaptranm + +subroutine psi_dtranidxm(ictxt,icomm,flag,n,beta,y,idx,totxch,totsnd,totrcv,work,info) + + use psi_mod, psb_protect_name => psi_dtranidxm + use psb_error_mod + use psb_descriptor_type + use psb_penv_mod +#ifdef MPI_MOD + use mpi +#endif + implicit none +#ifdef MPI_H + include 'mpif.h' +#endif + + integer, intent(in) :: ictxt,icomm,flag,n + integer, intent(out) :: info + real(kind(1.d0)) :: y(:,:), beta + real(kind(1.d0)), target :: work(:) + integer, intent(in) :: idx(:),totxch,totsnd, totrcv + + ! locals + integer :: np, me, nesd, nerv,& + & proc_to_comm, p2ptag, p2pstat(mpi_status_size),& + & iret, err_act, i, idx_pt, totsnd_, totrcv_,& + & snd_pt, rcv_pt, pnti, data_ + integer, allocatable, dimension(:) :: bsdidx, brvidx,& + & sdsz, rvsz, prcid, rvhd, sdhd + integer :: int_err(5) + logical :: swap_mpi, swap_sync, swap_send, swap_recv,& + & albf,do_send,do_recv + logical, parameter :: usersend=.false. + + real(kind(1.d0)), pointer, dimension(:) :: sndbuf, rcvbuf + character(len=20) :: name + + info = 0 + name='psi_swap_tran' + call psb_erractionsave(err_act) + + call psb_info(ictxt,me,np) + if (np == -1) then + info = 2010 + 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),& @@ -183,9 +234,9 @@ subroutine psi_dswaptranm(flag,n,beta,y,desc_a,work,info,data) snd_pt = 1 rcv_pt = 1 do i=1, totxch - proc_to_comm = d_idx(pnti+psb_proc_id_) - nerv = d_idx(pnti+psb_n_elem_recv_) - nesd = d_idx(pnti+nerv+psb_n_elem_send_) + proc_to_comm = idx(pnti+psb_proc_id_) + nerv = idx(pnti+psb_n_elem_recv_) + nesd = idx(pnti+nerv+psb_n_elem_send_) call psb_get_rank(prcid(proc_to_comm),ictxt,proc_to_comm) @@ -210,14 +261,14 @@ subroutine psi_dswaptranm(flag,n,beta,y,desc_a,work,info,data) end if - idxr = max(idxr,1) - idxs = max(idxs,1) - if((idxr+idxs) < size(work)) then - sndbuf => work(1:idxs) - rcvbuf => work(idxs+1:idxs+idxr) + 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(idxs),rcvbuf(idxr), stat=info) + allocate(sndbuf(totsnd_),rcvbuf(totrcv_), stat=info) if(info /= 0) then call psb_errpush(4000,name) goto 9999 @@ -232,12 +283,12 @@ subroutine psi_dswaptranm(flag,n,beta,y,desc_a,work,info,data) snd_pt = 1 rcv_pt = 1 do i=1, totxch - proc_to_comm = d_idx(pnti+psb_proc_id_) - nerv = d_idx(pnti+psb_n_elem_recv_) - nesd = d_idx(pnti+nerv+psb_n_elem_send_) + 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_gth(nerv,n,d_idx(idx_pt:idx_pt+nerv-1),& + call psi_gth(nerv,n,idx(idx_pt:idx_pt+nerv-1),& & y,rcvbuf(rcv_pt:rcv_pt+n*nerv-1)) rcv_pt = rcv_pt + n*nerv @@ -269,9 +320,9 @@ subroutine psi_dswaptranm(flag,n,beta,y,desc_a,work,info,data) snd_pt = 1 rcv_pt = 1 do i=1, totxch - proc_to_comm = d_idx(pnti+psb_proc_id_) - nerv = d_idx(pnti+psb_n_elem_recv_) - nesd = d_idx(pnti+nerv+psb_n_elem_send_) + 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 (nerv>0) call psb_snd(ictxt,& @@ -283,6 +334,11 @@ subroutine psi_dswaptranm(flag,n,beta,y,desc_a,work,info,data) & sndbuf(snd_pt:snd_pt+n*nesd-1), proc_to_comm) if (nerv>0) call psb_snd(ictxt,& & rcvbuf(rcv_pt:rcv_pt+n*nerv-1), proc_to_comm) + else if (proc_to_comm == me) then + if (nesd /= nerv) then + write(0,*) 'Fatal error in swaptran: mismatch on self sendf',nerv,nesd + end if + sndbuf(snd_pt:snd_pt+n*nesd-1) = rcvbuf(rcv_pt:rcv_pt+n*nerv-1) end if rcv_pt = rcv_pt + n*nerv snd_pt = snd_pt + n*nesd @@ -298,11 +354,11 @@ subroutine psi_dswaptranm(flag,n,beta,y,desc_a,work,info,data) snd_pt = 1 rcv_pt = 1 do i=1, totxch - proc_to_comm = d_idx(pnti+psb_proc_id_) - nerv = d_idx(pnti+psb_n_elem_recv_) - nesd = d_idx(pnti+nerv+psb_n_elem_send_) + proc_to_comm = idx(pnti+psb_proc_id_) + nerv = idx(pnti+psb_n_elem_recv_) + nesd = idx(pnti+nerv+psb_n_elem_send_) call psb_get_rank(prcid(i),ictxt,proc_to_comm) - if (nesd>0) then + if ((nesd>0).and.(proc_to_comm/=me)) then p2ptag = krecvid(ictxt,proc_to_comm,me) call mpi_irecv(sndbuf(snd_pt),n*nesd,& & mpi_double_precision,prcid(i),& @@ -321,11 +377,11 @@ subroutine psi_dswaptranm(flag,n,beta,y,desc_a,work,info,data) snd_pt = 1 rcv_pt = 1 do i=1, totxch - proc_to_comm = d_idx(pnti+psb_proc_id_) - nerv = d_idx(pnti+psb_n_elem_recv_) - nesd = d_idx(pnti+nerv+psb_n_elem_send_) + 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) then + if ((nerv>0).and.(proc_to_comm/=me)) then p2ptag=ksendid(ictxt,proc_to_comm,me) if (usersend) then call mpi_rsend(rcvbuf(rcv_pt),n*nerv,& @@ -334,9 +390,9 @@ subroutine psi_dswaptranm(flag,n,beta,y,desc_a,work,info,data) else call mpi_send(rcvbuf(rcv_pt),n*nerv,& & mpi_double_precision,prcid(i),& - & p2ptag,icomm,iret) + & p2ptag,icomm,iret) end if - + if(iret /= mpi_success) then int_err(1) = iret info=400 @@ -353,9 +409,9 @@ subroutine psi_dswaptranm(flag,n,beta,y,desc_a,work,info,data) pnti = 1 do i=1, totxch - proc_to_comm = d_idx(pnti+psb_proc_id_) - nerv = d_idx(pnti+psb_n_elem_recv_) - nesd = d_idx(pnti+nerv+psb_n_elem_send_) + proc_to_comm = idx(pnti+psb_proc_id_) + nerv = idx(pnti+psb_n_elem_recv_) + nesd = idx(pnti+nerv+psb_n_elem_send_) p2ptag = krecvid(ictxt,proc_to_comm,me) @@ -367,6 +423,11 @@ subroutine psi_dswaptranm(flag,n,beta,y,desc_a,work,info,data) call psb_errpush(info,name,i_err=int_err) goto 9999 end if + else if (proc_to_comm == me) then + if (nesd /= nerv) then + write(0,*) 'Fatal error in swaptran: mismatch on self sendf',nerv,nesd + end if + sndbuf(snd_pt:snd_pt+n*nesd-1) = rcvbuf(rcv_pt:rcv_pt+n*nerv-1) end if pnti = pnti + nerv + nesd + 3 end do @@ -378,9 +439,9 @@ subroutine psi_dswaptranm(flag,n,beta,y,desc_a,work,info,data) snd_pt = 1 rcv_pt = 1 do i=1, totxch - proc_to_comm = d_idx(pnti+psb_proc_id_) - nerv = d_idx(pnti+psb_n_elem_recv_) - nesd = d_idx(pnti+nerv+psb_n_elem_send_) + 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_snd(ictxt,& & rcvbuf(rcv_pt:rcv_pt+n*nerv-1), proc_to_comm) rcv_pt = rcv_pt + n*nerv @@ -395,9 +456,9 @@ subroutine psi_dswaptranm(flag,n,beta,y,desc_a,work,info,data) snd_pt = 1 rcv_pt = 1 do i=1, totxch - proc_to_comm = d_idx(pnti+psb_proc_id_) - nerv = d_idx(pnti+psb_n_elem_recv_) - nesd = d_idx(pnti+nerv+psb_n_elem_send_) + 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_rcv(ictxt,& & sndbuf(snd_pt:snd_pt+n*nesd-1), proc_to_comm) rcv_pt = rcv_pt + n*nerv @@ -413,11 +474,11 @@ subroutine psi_dswaptranm(flag,n,beta,y,desc_a,work,info,data) snd_pt = 1 rcv_pt = 1 do i=1, totxch - proc_to_comm = d_idx(pnti+psb_proc_id_) - nerv = d_idx(pnti+psb_n_elem_recv_) - nesd = d_idx(pnti+nerv+psb_n_elem_send_) + 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+nerv+psb_n_elem_send_ - call psi_sct(nesd,n,d_idx(idx_pt:idx_pt+nesd-1),& + call psi_sct(nesd,n,idx(idx_pt:idx_pt+nesd-1),& & sndbuf(snd_pt:snd_pt+n*nesd-1),beta,y) rcv_pt = rcv_pt + n*nerv snd_pt = snd_pt + n*nesd @@ -453,39 +514,7 @@ subroutine psi_dswaptranm(flag,n,beta,y,desc_a,work,info,data) return end if return -end subroutine psi_dswaptranm - - -!!$ -!!$ Parallel Sparse BLAS v2.0 -!!$ (C) Copyright 2006 Salvatore Filippone University of Rome Tor Vergata -!!$ Alfredo Buttari University of Rome Tor Vergata -!!$ -!!$ 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. -!!$ -!!$ +end subroutine psi_dtranidxm ! ! ! Subroutine: psi_dswaptranv @@ -547,7 +576,6 @@ subroutine psi_dswaptranv(flag,beta,y,desc_a,work,info,data) use psb_error_mod use psb_descriptor_type use psb_penv_mod -!!$ use psi_gthsct_mod #ifdef MPI_MOD use mpi #endif @@ -564,28 +592,17 @@ subroutine psi_dswaptranv(flag,beta,y,desc_a,work,info,data) integer, optional :: data ! locals - integer :: ictxt, np, me, nesd, nerv,& - & proc_to_comm, p2ptag, icomm, p2pstat(mpi_status_size),& - & idxs, idxr, iret, err_act, totxch, i, & - & idx_pt, snd_pt, rcv_pt, n, pnti, data_ - - integer, allocatable, dimension(:) :: bsdidx, brvidx,& - & sdsz, rvsz, prcid, rvhd, sdhd + integer :: ictxt, np, me, icomm, idxs, idxr, totxch, err_act, data_ integer, pointer :: d_idx(:) integer :: int_err(5) - logical :: swap_mpi, swap_sync, swap_send, swap_recv,& - & albf,do_send,do_recv - logical, parameter :: usersend=.false. - - real(kind(1.d0)), pointer, dimension(:) :: sndbuf, rcvbuf character(len=20) :: name info = 0 name='psi_swap_tranv' call psb_erractionsave(err_act) - ictxt=psb_cd_get_context(desc_a) - + ictxt = psb_cd_get_context(desc_a) + icomm = psb_cd_get_mpic(desc_a) call psb_info(ictxt,me,np) if (np == -1) then info = 2010 @@ -599,30 +616,91 @@ subroutine psi_dswaptranv(flag,beta,y,desc_a,work,info,data) goto 9999 endif - icomm = desc_a%matrix_data(psb_mpi_c_) - 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 - if (present(data)) then data_ = data else data_ = psb_comm_halo_ end if - + call psb_cd_get_list(data_,desc_a,d_idx,totxch,idxr,idxs,info) if (info /= 0) then call psb_errpush(4001,name,a_err='psb_cd_get_list') goto 9999 end if - idxr = idxr * n - idxs = idxs * n + call psi_swaptran(ictxt,icomm,flag,beta,y,d_idx,totxch,idxs,idxr,work,info) + if (info /= 0) goto 9999 + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + if (err_act == psb_act_abort_) then + call psb_error(ictxt) + return + end if + return +end subroutine psi_dswaptranv + + + +subroutine psi_dtranidxv(ictxt,icomm,flag,beta,y,idx,totxch,totsnd,totrcv,work,info) + + use psi_mod, psb_protect_name => psi_dtranidxv + use psb_error_mod + use psb_descriptor_type + use psb_penv_mod +#ifdef MPI_MOD + use mpi +#endif + implicit none +#ifdef MPI_H + include 'mpif.h' +#endif + + integer, intent(in) :: ictxt,icomm,flag + integer, intent(out) :: info + real(kind(1.d0)) :: y(:), beta + real(kind(1.d0)), target :: work(:) + integer, intent(in) :: idx(:),totxch,totsnd, totrcv + + ! locals + integer :: np, me, nesd, nerv,& + & proc_to_comm, p2ptag, p2pstat(mpi_status_size),& + & iret, err_act, i, idx_pt, totsnd_, totrcv_,& + & snd_pt, rcv_pt, pnti, data_, n + integer, allocatable, dimension(:) :: bsdidx, brvidx,& + & sdsz, rvsz, prcid, rvhd, sdhd + integer :: int_err(5) + logical :: swap_mpi, swap_sync, swap_send, swap_recv,& + & albf,do_send,do_recv + logical, parameter :: usersend=.false. + + real(kind(1.d0)), pointer, dimension(:) :: sndbuf, rcvbuf + character(len=20) :: name + + info = 0 + name='psi_swap_tran' + call psb_erractionsave(err_act) + + call psb_info(ictxt,me,np) + if (np == -1) then + info = 2010 + 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),& @@ -644,9 +722,9 @@ subroutine psi_dswaptranv(flag,beta,y,desc_a,work,info,data) snd_pt = 1 rcv_pt = 1 do i=1, totxch - proc_to_comm = d_idx(pnti+psb_proc_id_) - nerv = d_idx(pnti+psb_n_elem_recv_) - nesd = d_idx(pnti+nerv+psb_n_elem_send_) + proc_to_comm = idx(pnti+psb_proc_id_) + nerv = idx(pnti+psb_n_elem_recv_) + nesd = idx(pnti+nerv+psb_n_elem_send_) call psb_get_rank(prcid(proc_to_comm),ictxt,proc_to_comm) brvidx(proc_to_comm) = rcv_pt @@ -670,14 +748,14 @@ subroutine psi_dswaptranv(flag,beta,y,desc_a,work,info,data) end if - idxr = max(idxr,1) - idxs = max(idxs,1) - if((idxr+idxs) < size(work)) then - sndbuf => work(1:idxs) - rcvbuf => work(idxs+1:idxs+idxr) + 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(idxs),rcvbuf(idxr), stat=info) + allocate(sndbuf(totsnd_),rcvbuf(totrcv_), stat=info) if(info /= 0) then call psb_errpush(4000,name) goto 9999 @@ -693,12 +771,12 @@ subroutine psi_dswaptranv(flag,beta,y,desc_a,work,info,data) snd_pt = 1 rcv_pt = 1 do i=1, totxch - proc_to_comm = d_idx(pnti+psb_proc_id_) - nerv = d_idx(pnti+psb_n_elem_recv_) - nesd = d_idx(pnti+nerv+psb_n_elem_send_) + 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_gth(nerv,d_idx(idx_pt:idx_pt+nerv-1),& + call psi_gth(nerv,idx(idx_pt:idx_pt+nerv-1),& & y,rcvbuf(rcv_pt:rcv_pt+nerv-1)) rcv_pt = rcv_pt + nerv @@ -728,9 +806,9 @@ subroutine psi_dswaptranv(flag,beta,y,desc_a,work,info,data) snd_pt = 1 rcv_pt = 1 do i=1, totxch - proc_to_comm = d_idx(pnti+psb_proc_id_) - nerv = d_idx(pnti+psb_n_elem_recv_) - nesd = d_idx(pnti+nerv+psb_n_elem_send_) + 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 (nerv>0) call psb_snd(ictxt,& @@ -742,6 +820,11 @@ subroutine psi_dswaptranv(flag,beta,y,desc_a,work,info,data) & sndbuf(snd_pt:snd_pt+nesd-1), proc_to_comm) if (nerv>0) call psb_snd(ictxt,& & rcvbuf(rcv_pt:rcv_pt+nerv-1), proc_to_comm) + else if (proc_to_comm == me) then + if (nesd /= nerv) then + write(0,*) 'Fatal error in swaptran: mismatch on self sendf',nerv,nesd + end if + sndbuf(snd_pt:snd_pt+nesd-1) = rcvbuf(rcv_pt:rcv_pt+nerv-1) end if rcv_pt = rcv_pt + nerv snd_pt = snd_pt + nesd @@ -757,11 +840,11 @@ subroutine psi_dswaptranv(flag,beta,y,desc_a,work,info,data) snd_pt = 1 rcv_pt = 1 do i=1, totxch - proc_to_comm = d_idx(pnti+psb_proc_id_) - nerv = d_idx(pnti+psb_n_elem_recv_) - nesd = d_idx(pnti+nerv+psb_n_elem_send_) + proc_to_comm = idx(pnti+psb_proc_id_) + nerv = idx(pnti+psb_n_elem_recv_) + nesd = idx(pnti+nerv+psb_n_elem_send_) call psb_get_rank(prcid(i),ictxt,proc_to_comm) - if (nesd>0) then + if ((nesd>0).and.(proc_to_comm/=me)) then p2ptag = krecvid(ictxt,proc_to_comm,me) call mpi_irecv(sndbuf(snd_pt),nesd,& & mpi_double_precision,prcid(i),& @@ -780,11 +863,11 @@ subroutine psi_dswaptranv(flag,beta,y,desc_a,work,info,data) snd_pt = 1 rcv_pt = 1 do i=1, totxch - proc_to_comm = d_idx(pnti+psb_proc_id_) - nerv = d_idx(pnti+psb_n_elem_recv_) - nesd = d_idx(pnti+nerv+psb_n_elem_send_) + 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) then + if ((nerv>0).and.(proc_to_comm/=me)) then p2ptag=ksendid(ictxt,proc_to_comm,me) if (usersend) then call mpi_rsend(rcvbuf(rcv_pt),nerv,& @@ -812,11 +895,11 @@ subroutine psi_dswaptranv(flag,beta,y,desc_a,work,info,data) pnti = 1 do i=1, totxch - proc_to_comm = d_idx(pnti+psb_proc_id_) - nerv = d_idx(pnti+psb_n_elem_recv_) - nesd = d_idx(pnti+nerv+psb_n_elem_send_) + proc_to_comm = idx(pnti+psb_proc_id_) + nerv = idx(pnti+psb_n_elem_recv_) + nesd = idx(pnti+nerv+psb_n_elem_send_) p2ptag = krecvid(ictxt,proc_to_comm,me) - + if ((proc_to_comm /= me).and.(nesd>0)) then call mpi_wait(rvhd(i),p2pstat,iret) if(iret /= mpi_success) then @@ -825,10 +908,15 @@ subroutine psi_dswaptranv(flag,beta,y,desc_a,work,info,data) call psb_errpush(info,name,i_err=int_err) goto 9999 end if + else if (proc_to_comm == me) then + if (nesd /= nerv) then + write(0,*) 'Fatal error in swaptran: mismatch on self sendf',nerv,nesd + end if + sndbuf(snd_pt:snd_pt+nesd-1) = rcvbuf(rcv_pt:rcv_pt+nerv-1) end if pnti = pnti + nerv + nesd + 3 end do - + else if (swap_send) then @@ -836,9 +924,9 @@ subroutine psi_dswaptranv(flag,beta,y,desc_a,work,info,data) snd_pt = 1 rcv_pt = 1 do i=1, totxch - proc_to_comm = d_idx(pnti+psb_proc_id_) - nerv = d_idx(pnti+psb_n_elem_recv_) - nesd = d_idx(pnti+nerv+psb_n_elem_send_) + 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_snd(ictxt,& & rcvbuf(rcv_pt:rcv_pt+nerv-1), proc_to_comm) rcv_pt = rcv_pt + nerv @@ -853,9 +941,9 @@ subroutine psi_dswaptranv(flag,beta,y,desc_a,work,info,data) snd_pt = 1 rcv_pt = 1 do i=1, totxch - proc_to_comm = d_idx(pnti+psb_proc_id_) - nerv = d_idx(pnti+psb_n_elem_recv_) - nesd = d_idx(pnti+nerv+psb_n_elem_send_) + 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_rcv(ictxt,& & sndbuf(snd_pt:snd_pt+nesd-1), proc_to_comm) rcv_pt = rcv_pt + nerv @@ -873,11 +961,11 @@ subroutine psi_dswaptranv(flag,beta,y,desc_a,work,info,data) snd_pt = 1 rcv_pt = 1 do i=1, totxch - proc_to_comm = d_idx(pnti+psb_proc_id_) - nerv = d_idx(pnti+psb_n_elem_recv_) - nesd = d_idx(pnti+nerv+psb_n_elem_send_) + 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+nerv+psb_n_elem_send_ - call psi_sct(nesd,d_idx(idx_pt:idx_pt+nesd-1),& + call psi_sct(nesd,idx(idx_pt:idx_pt+nesd-1),& & sndbuf(snd_pt:snd_pt+nesd-1),beta,y) rcv_pt = rcv_pt + nerv snd_pt = snd_pt + nesd @@ -913,4 +1001,4 @@ subroutine psi_dswaptranv(flag,beta,y,desc_a,work,info,data) return end if return -end subroutine psi_dswaptranv +end subroutine psi_dtranidxv diff --git a/base/internals/psi_iswapdata.F90 b/base/internals/psi_iswapdata.F90 index 20f1a8d5..855095db 100644 --- a/base/internals/psi_iswapdata.F90 +++ b/base/internals/psi_iswapdata.F90 @@ -40,7 +40,7 @@ ! D real(kind(1.d0)) ! Z complex(kind(1.d0)) ! 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))); +! 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 @@ -86,7 +86,6 @@ subroutine psi_iswapdatam(flag,n,beta,y,desc_a,work,info,data) use psb_error_mod use psb_descriptor_type use psb_penv_mod -!!$ use psi_gthsct_mod #ifdef MPI_MOD use mpi #endif @@ -103,19 +102,8 @@ subroutine psi_iswapdatam(flag,n,beta,y,desc_a,work,info,data) integer, optional :: data ! locals - integer :: ictxt, np, me, nesd, nerv,& - & proc_to_comm, p2ptag, icomm, p2pstat(mpi_status_size),& - & idxs, idxr, iret, err_act, totxch, i, idx_pt,& - & snd_pt, rcv_pt, pnti, data_ - integer, allocatable, dimension(:) :: bsdidx, brvidx,& - & sdsz, rvsz, prcid, rvhd, sdhd + integer :: ictxt, np, me, icomm, idxs, idxr, totxch, data_, err_act integer, pointer :: d_idx(:) - integer :: int_err(5) - logical :: swap_mpi, swap_sync, swap_send, swap_recv,& - & albf,do_send,do_recv - logical, parameter :: usersend=.false. - - integer, pointer, dimension(:) :: sndbuf, rcvbuf character(len=20) :: name info = 0 @@ -124,7 +112,6 @@ subroutine psi_iswapdatam(flag,n,beta,y,desc_a,work,info,data) ictxt = psb_cd_get_context(desc_a) icomm = psb_cd_get_mpic(desc_a) - call psb_info(ictxt,me,np) if (np == -1) then info = 2010 @@ -138,13 +125,6 @@ subroutine psi_iswapdatam(flag,n,beta,y,desc_a,work,info,data) 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 - if(present(data)) then data_ = data else @@ -157,8 +137,78 @@ subroutine psi_iswapdatam(flag,n,beta,y,desc_a,work,info,data) goto 9999 end if - idxr = idxr * n - idxs = idxs * n + + call psi_swapdata(ictxt,icomm,flag,n,beta,y,d_idx,totxch,idxs,idxr,work,info) + if (info /= 0) goto 9999 + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + if (err_act == psb_act_abort_) then + call psb_error(ictxt) + return + end if + return +end subroutine psi_iswapdatam + +subroutine psi_iswapidxm(ictxt,icomm,flag,n,beta,y,idx,totxch,totsnd,totrcv,work,info) + + use psi_mod, psb_protect_name => psi_iswapidxm + use psb_error_mod + use psb_descriptor_type + use psb_penv_mod +#ifdef MPI_MOD + use mpi +#endif + implicit none +#ifdef MPI_H + include 'mpif.h' +#endif + + integer, intent(in) :: ictxt,icomm,flag,n + integer, intent(out) :: info + integer :: y(:,:), beta + integer, target :: work(:) + integer, intent(in) :: idx(:),totxch,totsnd, totrcv + + ! locals + integer :: np, me, nesd, nerv,& + & proc_to_comm, p2ptag, p2pstat(mpi_status_size),& + & iret, err_act, i, idx_pt, totsnd_, totrcv_,& + & snd_pt, rcv_pt, pnti, data_ + integer, allocatable, dimension(:) :: bsdidx, brvidx,& + & sdsz, rvsz, prcid, rvhd, sdhd + integer :: int_err(5) + logical :: swap_mpi, swap_sync, swap_send, swap_recv,& + & albf,do_send,do_recv + logical, parameter :: usersend=.false. + + integer, pointer, dimension(:) :: sndbuf, rcvbuf + character(len=20) :: name + + info = 0 + name='psi_swap_data' + call psb_erractionsave(err_act) + + call psb_info(ictxt,me,np) + if (np == -1) then + info = 2010 + 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),& @@ -180,9 +230,9 @@ subroutine psi_iswapdatam(flag,n,beta,y,desc_a,work,info,data) snd_pt = 1 rcv_pt = 1 do i=1, totxch - proc_to_comm = d_idx(pnti+psb_proc_id_) - nerv = d_idx(pnti+psb_n_elem_recv_) - nesd = d_idx(pnti+nerv+psb_n_elem_send_) + proc_to_comm = idx(pnti+psb_proc_id_) + nerv = idx(pnti+psb_n_elem_recv_) + nesd = idx(pnti+nerv+psb_n_elem_send_) call psb_get_rank(prcid(proc_to_comm),ictxt,proc_to_comm) brvidx(proc_to_comm) = rcv_pt @@ -205,14 +255,14 @@ subroutine psi_iswapdatam(flag,n,beta,y,desc_a,work,info,data) end if end if - idxr = max(idxr,1) - idxs = max(idxs,1) - if((idxr+idxs) < size(work)) then - sndbuf => work(1:idxs) - rcvbuf => work(idxs+1:idxs+idxr) + 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(idxs),rcvbuf(idxr), stat=info) + allocate(sndbuf(totsnd_),rcvbuf(totrcv_), stat=info) if(info /= 0) then call psb_errpush(4000,name) goto 9999 @@ -226,10 +276,10 @@ subroutine psi_iswapdatam(flag,n,beta,y,desc_a,work,info,data) pnti = 1 snd_pt = 1 do i=1, totxch - nerv = d_idx(pnti+psb_n_elem_recv_) - nesd = d_idx(pnti+nerv+psb_n_elem_send_) + 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,d_idx(idx_pt:idx_pt+nesd-1),& + 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 @@ -260,9 +310,9 @@ subroutine psi_iswapdatam(flag,n,beta,y,desc_a,work,info,data) snd_pt = 1 rcv_pt = 1 do i=1, totxch - proc_to_comm = d_idx(pnti+psb_proc_id_) - nerv = d_idx(pnti+psb_n_elem_recv_) - nesd = d_idx(pnti+nerv+psb_n_elem_send_) + 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,& @@ -274,6 +324,11 @@ subroutine psi_iswapdatam(flag,n,beta,y,desc_a,work,info,data) & 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(0,*) 'Fatal error in swapdata: mismatch on self sendf',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 @@ -291,11 +346,11 @@ subroutine psi_iswapdatam(flag,n,beta,y,desc_a,work,info,data) snd_pt = 1 rcv_pt = 1 do i=1, totxch - proc_to_comm = d_idx(pnti+psb_proc_id_) - nerv = d_idx(pnti+psb_n_elem_recv_) - nesd = d_idx(pnti+nerv+psb_n_elem_send_) + proc_to_comm = idx(pnti+psb_proc_id_) + nerv = idx(pnti+psb_n_elem_recv_) + nesd = idx(pnti+nerv+psb_n_elem_send_) call psb_get_rank(prcid(i),ictxt,proc_to_comm) - if (nerv>0) then + if ((nerv>0).and.(proc_to_comm/=me)) then p2ptag = krecvid(ictxt,proc_to_comm,me) call mpi_irecv(rcvbuf(rcv_pt),n*nerv,& & mpi_integer,prcid(i),& @@ -315,12 +370,12 @@ subroutine psi_iswapdatam(flag,n,beta,y,desc_a,work,info,data) snd_pt = 1 rcv_pt = 1 do i=1, totxch - proc_to_comm = d_idx(pnti+psb_proc_id_) - nerv = d_idx(pnti+psb_n_elem_recv_) - nesd = d_idx(pnti+nerv+psb_n_elem_send_) + proc_to_comm = idx(pnti+psb_proc_id_) + nerv = idx(pnti+psb_n_elem_recv_) + nesd = idx(pnti+nerv+psb_n_elem_send_) p2ptag=ksendid(ictxt,proc_to_comm,me) - if (nesd>0) then + if ((nesd>0).and.(proc_to_comm/=me)) then if (usersend) then call mpi_rsend(sndbuf(snd_pt),n*nesd,& & mpi_integer,prcid(i),& @@ -348,9 +403,9 @@ subroutine psi_iswapdatam(flag,n,beta,y,desc_a,work,info,data) pnti = 1 do i=1, totxch - proc_to_comm = d_idx(pnti+psb_proc_id_) - nerv = d_idx(pnti+psb_n_elem_recv_) - nesd = d_idx(pnti+nerv+psb_n_elem_send_) + proc_to_comm = idx(pnti+psb_proc_id_) + nerv = idx(pnti+psb_n_elem_recv_) + nesd = idx(pnti+nerv+psb_n_elem_send_) p2ptag = krecvid(ictxt,proc_to_comm,me) @@ -362,10 +417,15 @@ subroutine psi_iswapdatam(flag,n,beta,y,desc_a,work,info,data) call psb_errpush(info,name,i_err=int_err) goto 9999 end if + else if (proc_to_comm == me) then + if (nesd /= nerv) then + write(0,*) 'Fatal error in swapdata: mismatch on self sendf',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 @@ -374,9 +434,9 @@ subroutine psi_iswapdatam(flag,n,beta,y,desc_a,work,info,data) snd_pt = 1 rcv_pt = 1 do i=1, totxch - proc_to_comm = d_idx(pnti+psb_proc_id_) - nerv = d_idx(pnti+psb_n_elem_recv_) - nesd = d_idx(pnti+nerv+psb_n_elem_send_) + 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) @@ -393,9 +453,9 @@ subroutine psi_iswapdatam(flag,n,beta,y,desc_a,work,info,data) snd_pt = 1 rcv_pt = 1 do i=1, totxch - proc_to_comm = d_idx(pnti+psb_proc_id_) - nerv = d_idx(pnti+psb_n_elem_recv_) - nesd = d_idx(pnti+nerv+psb_n_elem_send_) + 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 @@ -414,11 +474,11 @@ subroutine psi_iswapdatam(flag,n,beta,y,desc_a,work,info,data) snd_pt = 1 rcv_pt = 1 do i=1, totxch - proc_to_comm = d_idx(pnti+psb_proc_id_) - nerv = d_idx(pnti+psb_n_elem_recv_) - nesd = d_idx(pnti+nerv+psb_n_elem_send_) + 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,d_idx(idx_pt:idx_pt+nerv-1),& + 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 @@ -454,39 +514,8 @@ subroutine psi_iswapdatam(flag,n,beta,y,desc_a,work,info,data) return end if return -end subroutine psi_iswapdatam +end subroutine psi_iswapidxm - -!!$ -!!$ Parallel Sparse BLAS v2.0 -!!$ (C) Copyright 2006 Salvatore Filippone University of Rome Tor Vergata -!!$ Alfredo Buttari University of Rome Tor Vergata -!!$ -!!$ 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_iswapdatav @@ -544,7 +573,6 @@ subroutine psi_iswapdatav(flag,beta,y,desc_a,work,info,data) use psb_error_mod use psb_descriptor_type use psb_penv_mod -!!$ use psi_gthsct_mod #ifdef MPI_MOD use mpi #endif @@ -561,20 +589,8 @@ subroutine psi_iswapdatav(flag,beta,y,desc_a,work,info,data) integer, optional :: data ! locals - integer :: ictxt, np, me, nesd, nerv,& - & proc_to_comm, p2ptag, icomm, p2pstat(mpi_status_size),& - & idxs, idxr, iret, err_act, totxch, i, & - & idx_pt, snd_pt, rcv_pt, n, pnti, data_ - - integer, allocatable, dimension(:) :: bsdidx, brvidx,& - & sdsz, rvsz, prcid, rvhd, sdhd + integer :: ictxt, np, me, icomm, idxs, idxr, totxch, data_, err_act integer, pointer :: d_idx(:) - integer :: int_err(5) - logical :: swap_mpi, swap_sync, swap_send, swap_recv,& - & albf,do_send,do_recv - logical, parameter :: usersend=.false. - - integer, pointer, dimension(:) :: sndbuf, rcvbuf character(len=20) :: name info = 0 @@ -595,16 +611,7 @@ subroutine psi_iswapdatav(flag,beta,y,desc_a,work,info,data) goto 9999 endif - icomm = desc_a%matrix_data(psb_mpi_c_) - 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 + icomm = psb_cd_get_mpic(desc_a) if(present(data)) then data_ = data @@ -618,8 +625,80 @@ subroutine psi_iswapdatav(flag,beta,y,desc_a,work,info,data) goto 9999 end if - idxr = idxr * n - idxs = idxs * n + call psi_swapdata(ictxt,icomm,flag,beta,y,d_idx,totxch,idxs,idxr,work,info) + if (info /= 0) goto 9999 + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + if (err_act == psb_act_abort_) then + call psb_error(ictxt) + return + end if + return +end subroutine psi_iswapdatav + + +subroutine psi_iswapidxv(ictxt,icomm,flag,beta,y,idx,totxch,totsnd,totrcv,work,info) + + use psi_mod, psb_protect_name => psi_iswapidxv + use psb_error_mod + use psb_descriptor_type + use psb_penv_mod +#ifdef MPI_MOD + use mpi +#endif + implicit none +#ifdef MPI_H + include 'mpif.h' +#endif + + integer, intent(in) :: ictxt,icomm,flag + integer, intent(out) :: info + integer :: y(:), beta + integer, target :: work(:) + integer, intent(in) :: idx(:),totxch,totsnd, totrcv + + ! locals + integer :: np, me, nesd, nerv,& + & proc_to_comm, p2ptag, p2pstat(mpi_status_size),& + & iret, err_act, i, totsnd_, totrcv_,& + & idx_pt, snd_pt, rcv_pt, n, pnti, data_ + + integer, allocatable, dimension(:) :: bsdidx, brvidx,& + & sdsz, rvsz, prcid, rvhd, sdhd + integer :: int_err(5) + logical :: swap_mpi, swap_sync, swap_send, swap_recv,& + & albf,do_send,do_recv + logical, parameter :: usersend=.false. + + integer, pointer, dimension(:) :: sndbuf, rcvbuf + character(len=20) :: name + + info = 0 + name='psi_swap_datav' + call psb_erractionsave(err_act) + + call psb_info(ictxt,me,np) + if (np == -1) then + info = 2010 + 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),& @@ -640,9 +719,9 @@ subroutine psi_iswapdatav(flag,beta,y,desc_a,work,info,data) snd_pt = 1 rcv_pt = 1 do i=1, totxch - proc_to_comm = d_idx(pnti+psb_proc_id_) - nerv = d_idx(pnti+psb_n_elem_recv_) - nesd = d_idx(pnti+nerv+psb_n_elem_send_) + proc_to_comm = idx(pnti+psb_proc_id_) + nerv = idx(pnti+psb_n_elem_recv_) + nesd = idx(pnti+nerv+psb_n_elem_send_) call psb_get_rank(prcid(proc_to_comm),ictxt,proc_to_comm) brvidx(proc_to_comm) = rcv_pt @@ -666,14 +745,14 @@ subroutine psi_iswapdatav(flag,beta,y,desc_a,work,info,data) end if - idxr = max(idxr,1) - idxs = max(idxs,1) - if((idxr+idxs) < size(work)) then - sndbuf => work(1:idxs) - rcvbuf => work(idxs+1:idxs+idxr) + 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(idxs),rcvbuf(idxr), stat=info) + allocate(sndbuf(totsnd_),rcvbuf(totrcv_), stat=info) if(info /= 0) then call psb_errpush(4000,name) goto 9999 @@ -688,10 +767,10 @@ subroutine psi_iswapdatav(flag,beta,y,desc_a,work,info,data) pnti = 1 snd_pt = 1 do i=1, totxch - nerv = d_idx(pnti+psb_n_elem_recv_) - nesd = d_idx(pnti+nerv+psb_n_elem_send_) + 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,d_idx(idx_pt:idx_pt+nesd-1),& + 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 @@ -719,9 +798,9 @@ subroutine psi_iswapdatav(flag,beta,y,desc_a,work,info,data) snd_pt = 1 rcv_pt = 1 do i=1, totxch - proc_to_comm = d_idx(pnti+psb_proc_id_) - nerv = d_idx(pnti+psb_n_elem_recv_) - nesd = d_idx(pnti+nerv+psb_n_elem_send_) + 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,& @@ -733,6 +812,11 @@ subroutine psi_iswapdatav(flag,beta,y,desc_a,work,info,data) & 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(0,*) 'Fatal error in swapdata: mismatch on self sendf',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 @@ -747,12 +831,12 @@ subroutine psi_iswapdatav(flag,beta,y,desc_a,work,info,data) snd_pt = 1 rcv_pt = 1 do i=1, totxch - proc_to_comm = d_idx(pnti+psb_proc_id_) - nerv = d_idx(pnti+psb_n_elem_recv_) - nesd = d_idx(pnti+nerv+psb_n_elem_send_) + proc_to_comm = idx(pnti+psb_proc_id_) + nerv = idx(pnti+psb_n_elem_recv_) + nesd = idx(pnti+nerv+psb_n_elem_send_) call psb_get_rank(prcid(i),ictxt,proc_to_comm) - if (nerv>0) then + if ((nerv>0).and.(proc_to_comm /= me)) then p2ptag = krecvid(ictxt,proc_to_comm,me) call mpi_irecv(rcvbuf(rcv_pt),nerv,& & mpi_integer,prcid(i),& @@ -771,13 +855,13 @@ subroutine psi_iswapdatav(flag,beta,y,desc_a,work,info,data) snd_pt = 1 rcv_pt = 1 do i=1, totxch - proc_to_comm = d_idx(pnti+psb_proc_id_) - nerv = d_idx(pnti+psb_n_elem_recv_) - nesd = d_idx(pnti+nerv+psb_n_elem_send_) + proc_to_comm = idx(pnti+psb_proc_id_) + nerv = idx(pnti+psb_n_elem_recv_) + nesd = idx(pnti+nerv+psb_n_elem_send_) p2ptag=ksendid(ictxt,proc_to_comm,me) - if (nesd>0) then + if ((nesd>0).and.(proc_to_comm /= me)) then if (usersend) then call mpi_rsend(sndbuf(snd_pt),nesd,& & mpi_integer,prcid(i),& @@ -803,9 +887,9 @@ subroutine psi_iswapdatav(flag,beta,y,desc_a,work,info,data) pnti = 1 do i=1, totxch - proc_to_comm = d_idx(pnti+psb_proc_id_) - nerv = d_idx(pnti+psb_n_elem_recv_) - nesd = d_idx(pnti+nerv+psb_n_elem_send_) + proc_to_comm = idx(pnti+psb_proc_id_) + nerv = idx(pnti+psb_n_elem_recv_) + nesd = idx(pnti+nerv+psb_n_elem_send_) p2ptag = krecvid(ictxt,proc_to_comm,me) @@ -817,10 +901,15 @@ subroutine psi_iswapdatav(flag,beta,y,desc_a,work,info,data) call psb_errpush(info,name,i_err=int_err) goto 9999 end if + else if (proc_to_comm == me) then + if (nesd /= nerv) then + write(0,*) 'Fatal error in swapdata: mismatch on self sendf',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 @@ -828,9 +917,9 @@ subroutine psi_iswapdatav(flag,beta,y,desc_a,work,info,data) snd_pt = 1 rcv_pt = 1 do i=1, totxch - proc_to_comm = d_idx(pnti+psb_proc_id_) - nerv = d_idx(pnti+psb_n_elem_recv_) - nesd = d_idx(pnti+nerv+psb_n_elem_send_) + 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 @@ -844,9 +933,9 @@ subroutine psi_iswapdatav(flag,beta,y,desc_a,work,info,data) snd_pt = 1 rcv_pt = 1 do i=1, totxch - proc_to_comm = d_idx(pnti+psb_proc_id_) - nerv = d_idx(pnti+psb_n_elem_recv_) - nesd = d_idx(pnti+nerv+psb_n_elem_send_) + 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 @@ -856,18 +945,17 @@ subroutine psi_iswapdatav(flag,beta,y,desc_a,work,info,data) end if - if (do_recv) then pnti = 1 snd_pt = 1 rcv_pt = 1 do i=1, totxch - proc_to_comm = d_idx(pnti+psb_proc_id_) - nerv = d_idx(pnti+psb_n_elem_recv_) - nesd = d_idx(pnti+nerv+psb_n_elem_send_) + 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,d_idx(idx_pt:idx_pt+nerv-1),& + 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 @@ -902,4 +990,4 @@ subroutine psi_iswapdatav(flag,beta,y,desc_a,work,info,data) return end if return -end subroutine psi_iswapdatav +end subroutine psi_iswapidxv diff --git a/base/internals/psi_iswaptran.F90 b/base/internals/psi_iswaptran.F90 index de37f02a..6f16e8d3 100644 --- a/base/internals/psi_iswaptran.F90 +++ b/base/internals/psi_iswaptran.F90 @@ -90,7 +90,6 @@ subroutine psi_iswaptranm(flag,n,beta,y,desc_a,work,info,data) use psb_error_mod use psb_descriptor_type use psb_penv_mod -!!$ use psi_gthsct_mod #ifdef MPI_MOD use mpi #endif @@ -107,19 +106,9 @@ subroutine psi_iswaptranm(flag,n,beta,y,desc_a,work,info,data) integer, optional :: data ! locals - integer :: ictxt, np, me, nesd, nerv,& - & proc_to_comm, p2ptag, icomm, p2pstat(mpi_status_size),& - & idxs, idxr, iret, err_act, totxch, i, idx_pt,& - & snd_pt, rcv_pt, pnti, data_ - integer, allocatable, dimension(:) :: bsdidx, brvidx,& - & sdsz, rvsz, prcid, rvhd, sdhd + integer :: ictxt, np, me, icomm, idxs, idxr, err_act, totxch, data_ integer, pointer :: d_idx(:) integer :: int_err(5) - logical :: swap_mpi, swap_sync, swap_send, swap_recv,& - & albf,do_send,do_recv - logical, parameter :: usersend=.false. - - integer, pointer, dimension(:) :: sndbuf, rcvbuf character(len=20) :: name info = 0 @@ -128,6 +117,7 @@ subroutine psi_iswaptranm(flag,n,beta,y,desc_a,work,info,data) ictxt = psb_cd_get_context(desc_a) icomm = psb_cd_get_mpic(desc_a) + call psb_info(ictxt,me,np) if (np == -1) then info = 2010 @@ -141,13 +131,6 @@ subroutine psi_iswaptranm(flag,n,beta,y,desc_a,work,info,data) 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 - if(present(data)) then data_ = data else @@ -160,8 +143,77 @@ subroutine psi_iswaptranm(flag,n,beta,y,desc_a,work,info,data) goto 9999 end if - idxr = idxr * n - idxs = idxs * n + call psi_swaptran(ictxt,icomm,flag,n,beta,y,d_idx,totxch,idxs,idxr,work,info) + if (info /= 0) goto 9999 + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + if (err_act == psb_act_abort_) then + call psb_error(ictxt) + return + end if + return +end subroutine psi_iswaptranm + +subroutine psi_itranidxm(ictxt,icomm,flag,n,beta,y,idx,totxch,totsnd,totrcv,work,info) + + use psi_mod, psb_protect_name => psi_itranidxm + use psb_error_mod + use psb_descriptor_type + use psb_penv_mod +#ifdef MPI_MOD + use mpi +#endif + implicit none +#ifdef MPI_H + include 'mpif.h' +#endif + + integer, intent(in) :: ictxt,icomm,flag,n + integer, intent(out) :: info + integer :: y(:,:), beta + integer, target :: work(:) + integer, intent(in) :: idx(:),totxch,totsnd, totrcv + + ! locals + integer :: np, me, nesd, nerv,& + & proc_to_comm, p2ptag, p2pstat(mpi_status_size),& + & iret, err_act, i, idx_pt, totsnd_, totrcv_,& + & snd_pt, rcv_pt, pnti, data_ + integer, allocatable, dimension(:) :: bsdidx, brvidx,& + & sdsz, rvsz, prcid, rvhd, sdhd + integer :: int_err(5) + logical :: swap_mpi, swap_sync, swap_send, swap_recv,& + & albf,do_send,do_recv + logical, parameter :: usersend=.false. + + integer, pointer, dimension(:) :: sndbuf, rcvbuf + character(len=20) :: name + + info = 0 + name='psi_swap_tran' + call psb_erractionsave(err_act) + + call psb_info(ictxt,me,np) + if (np == -1) then + info = 2010 + 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),& @@ -182,9 +234,9 @@ subroutine psi_iswaptranm(flag,n,beta,y,desc_a,work,info,data) snd_pt = 1 rcv_pt = 1 do i=1, totxch - proc_to_comm = d_idx(pnti+psb_proc_id_) - nerv = d_idx(pnti+psb_n_elem_recv_) - nesd = d_idx(pnti+nerv+psb_n_elem_send_) + proc_to_comm = idx(pnti+psb_proc_id_) + nerv = idx(pnti+psb_n_elem_recv_) + nesd = idx(pnti+nerv+psb_n_elem_send_) call psb_get_rank(prcid(proc_to_comm),ictxt,proc_to_comm) @@ -209,14 +261,14 @@ subroutine psi_iswaptranm(flag,n,beta,y,desc_a,work,info,data) end if - idxr = max(idxr,1) - idxs = max(idxs,1) - if((idxr+idxs) < size(work)) then - sndbuf => work(1:idxs) - rcvbuf => work(idxs+1:idxs+idxr) + 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(idxs),rcvbuf(idxr), stat=info) + allocate(sndbuf(totsnd_),rcvbuf(totrcv_), stat=info) if(info /= 0) then call psb_errpush(4000,name) goto 9999 @@ -231,12 +283,12 @@ subroutine psi_iswaptranm(flag,n,beta,y,desc_a,work,info,data) snd_pt = 1 rcv_pt = 1 do i=1, totxch - proc_to_comm = d_idx(pnti+psb_proc_id_) - nerv = d_idx(pnti+psb_n_elem_recv_) - nesd = d_idx(pnti+nerv+psb_n_elem_send_) + 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_gth(nerv,n,d_idx(idx_pt:idx_pt+nerv-1),& + call psi_gth(nerv,n,idx(idx_pt:idx_pt+nerv-1),& & y,rcvbuf(rcv_pt:rcv_pt+n*nerv-1)) rcv_pt = rcv_pt + n*nerv @@ -268,9 +320,9 @@ subroutine psi_iswaptranm(flag,n,beta,y,desc_a,work,info,data) snd_pt = 1 rcv_pt = 1 do i=1, totxch - proc_to_comm = d_idx(pnti+psb_proc_id_) - nerv = d_idx(pnti+psb_n_elem_recv_) - nesd = d_idx(pnti+nerv+psb_n_elem_send_) + 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 (nerv>0) call psb_snd(ictxt,& @@ -282,6 +334,11 @@ subroutine psi_iswaptranm(flag,n,beta,y,desc_a,work,info,data) & sndbuf(snd_pt:snd_pt+n*nesd-1), proc_to_comm) if (nerv>0) call psb_snd(ictxt,& & rcvbuf(rcv_pt:rcv_pt+n*nerv-1), proc_to_comm) + else if (proc_to_comm == me) then + if (nesd /= nerv) then + write(0,*) 'Fatal error in swaptran: mismatch on self sendf',nerv,nesd + end if + sndbuf(snd_pt:snd_pt+n*nesd-1) = rcvbuf(rcv_pt:rcv_pt+n*nerv-1) end if rcv_pt = rcv_pt + n*nerv snd_pt = snd_pt + n*nesd @@ -297,11 +354,11 @@ subroutine psi_iswaptranm(flag,n,beta,y,desc_a,work,info,data) snd_pt = 1 rcv_pt = 1 do i=1, totxch - proc_to_comm = d_idx(pnti+psb_proc_id_) - nerv = d_idx(pnti+psb_n_elem_recv_) - nesd = d_idx(pnti+nerv+psb_n_elem_send_) + proc_to_comm = idx(pnti+psb_proc_id_) + nerv = idx(pnti+psb_n_elem_recv_) + nesd = idx(pnti+nerv+psb_n_elem_send_) call psb_get_rank(prcid(i),ictxt,proc_to_comm) - if (nesd>0) then + if ((nesd>0).and.(proc_to_comm/=me)) then p2ptag = krecvid(ictxt,proc_to_comm,me) call mpi_irecv(sndbuf(snd_pt),n*nesd,& & mpi_integer,prcid(i),& @@ -320,11 +377,11 @@ subroutine psi_iswaptranm(flag,n,beta,y,desc_a,work,info,data) snd_pt = 1 rcv_pt = 1 do i=1, totxch - proc_to_comm = d_idx(pnti+psb_proc_id_) - nerv = d_idx(pnti+psb_n_elem_recv_) - nesd = d_idx(pnti+nerv+psb_n_elem_send_) + 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) then + if ((nerv>0).and.(proc_to_comm/=me)) then p2ptag=ksendid(ictxt,proc_to_comm,me) if (usersend) then call mpi_rsend(rcvbuf(rcv_pt),n*nerv,& @@ -352,9 +409,9 @@ subroutine psi_iswaptranm(flag,n,beta,y,desc_a,work,info,data) pnti = 1 do i=1, totxch - proc_to_comm = d_idx(pnti+psb_proc_id_) - nerv = d_idx(pnti+psb_n_elem_recv_) - nesd = d_idx(pnti+nerv+psb_n_elem_send_) + proc_to_comm = idx(pnti+psb_proc_id_) + nerv = idx(pnti+psb_n_elem_recv_) + nesd = idx(pnti+nerv+psb_n_elem_send_) p2ptag = krecvid(ictxt,proc_to_comm,me) @@ -366,6 +423,11 @@ subroutine psi_iswaptranm(flag,n,beta,y,desc_a,work,info,data) call psb_errpush(info,name,i_err=int_err) goto 9999 end if + else if (proc_to_comm == me) then + if (nesd /= nerv) then + write(0,*) 'Fatal error in swaptran: mismatch on self sendf',nerv,nesd + end if + sndbuf(snd_pt:snd_pt+n*nesd-1) = rcvbuf(rcv_pt:rcv_pt+n*nerv-1) end if pnti = pnti + nerv + nesd + 3 end do @@ -377,9 +439,9 @@ subroutine psi_iswaptranm(flag,n,beta,y,desc_a,work,info,data) snd_pt = 1 rcv_pt = 1 do i=1, totxch - proc_to_comm = d_idx(pnti+psb_proc_id_) - nerv = d_idx(pnti+psb_n_elem_recv_) - nesd = d_idx(pnti+nerv+psb_n_elem_send_) + 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_snd(ictxt,& & rcvbuf(rcv_pt:rcv_pt+n*nerv-1), proc_to_comm) rcv_pt = rcv_pt + n*nerv @@ -394,9 +456,9 @@ subroutine psi_iswaptranm(flag,n,beta,y,desc_a,work,info,data) snd_pt = 1 rcv_pt = 1 do i=1, totxch - proc_to_comm = d_idx(pnti+psb_proc_id_) - nerv = d_idx(pnti+psb_n_elem_recv_) - nesd = d_idx(pnti+nerv+psb_n_elem_send_) + 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_rcv(ictxt,& & sndbuf(snd_pt:snd_pt+n*nesd-1), proc_to_comm) rcv_pt = rcv_pt + n*nerv @@ -412,11 +474,11 @@ subroutine psi_iswaptranm(flag,n,beta,y,desc_a,work,info,data) snd_pt = 1 rcv_pt = 1 do i=1, totxch - proc_to_comm = d_idx(pnti+psb_proc_id_) - nerv = d_idx(pnti+psb_n_elem_recv_) - nesd = d_idx(pnti+nerv+psb_n_elem_send_) + 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+nerv+psb_n_elem_send_ - call psi_sct(nesd,n,d_idx(idx_pt:idx_pt+nesd-1),& + call psi_sct(nesd,n,idx(idx_pt:idx_pt+nesd-1),& & sndbuf(snd_pt:snd_pt+n*nesd-1),beta,y) rcv_pt = rcv_pt + n*nerv snd_pt = snd_pt + n*nesd @@ -452,39 +514,7 @@ subroutine psi_iswaptranm(flag,n,beta,y,desc_a,work,info,data) return end if return -end subroutine psi_iswaptranm - - -!!$ -!!$ Parallel Sparse BLAS v2.0 -!!$ (C) Copyright 2006 Salvatore Filippone University of Rome Tor Vergata -!!$ Alfredo Buttari University of Rome Tor Vergata -!!$ -!!$ 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. -!!$ -!!$ +end subroutine psi_itranidxm ! ! ! Subroutine: psi_iswaptranv @@ -546,7 +576,6 @@ subroutine psi_iswaptranv(flag,beta,y,desc_a,work,info,data) use psb_error_mod use psb_descriptor_type use psb_penv_mod -!!$ use psi_gthsct_mod #ifdef MPI_MOD use mpi #endif @@ -563,27 +592,17 @@ subroutine psi_iswaptranv(flag,beta,y,desc_a,work,info,data) integer, optional :: data ! locals - integer :: ictxt, np, me, nesd, nerv,& - & proc_to_comm, p2ptag, icomm, p2pstat(mpi_status_size),& - & idxs, idxr, iret, err_act, totxch, i, & - & idx_pt, snd_pt, rcv_pt, n, pnti, data_ - - integer, allocatable, dimension(:) :: bsdidx, brvidx,& - & sdsz, rvsz, prcid, rvhd, sdhd + integer :: ictxt, np, me, icomm, idxs, idxr, totxch, err_act, data_ integer, pointer :: d_idx(:) integer :: int_err(5) - logical :: swap_mpi, swap_sync, swap_send, swap_recv,& - & albf,do_send,do_recv - logical, parameter :: usersend=.false. - - integer, pointer, dimension(:) :: sndbuf, rcvbuf character(len=20) :: name info = 0 name='psi_swap_tranv' call psb_erractionsave(err_act) - ictxt=psb_cd_get_context(desc_a) + ictxt = psb_cd_get_context(desc_a) + icomm = psb_cd_get_mpic(desc_a) call psb_info(ictxt,me,np) if (np == -1) then info = 2010 @@ -597,30 +616,91 @@ subroutine psi_iswaptranv(flag,beta,y,desc_a,work,info,data) goto 9999 endif - icomm = desc_a%matrix_data(psb_mpi_c_) - 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 - - if(present(data)) then + if (present(data)) then data_ = data else data_ = psb_comm_halo_ end if - + call psb_cd_get_list(data_,desc_a,d_idx,totxch,idxr,idxs,info) if (info /= 0) then call psb_errpush(4001,name,a_err='psb_cd_get_list') goto 9999 end if - idxr = idxr * n - idxs = idxs * n + call psi_swaptran(ictxt,icomm,flag,beta,y,d_idx,totxch,idxs,idxr,work,info) + if (info /= 0) goto 9999 + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + if (err_act == psb_act_abort_) then + call psb_error(ictxt) + return + end if + return +end subroutine psi_iswaptranv + + + +subroutine psi_itranidxv(ictxt,icomm,flag,beta,y,idx,totxch,totsnd,totrcv,work,info) + + use psi_mod, psb_protect_name => psi_itranidxv + use psb_error_mod + use psb_descriptor_type + use psb_penv_mod +#ifdef MPI_MOD + use mpi +#endif + implicit none +#ifdef MPI_H + include 'mpif.h' +#endif + + integer, intent(in) :: ictxt,icomm,flag + integer, intent(out) :: info + integer :: y(:), beta + integer, target :: work(:) + integer, intent(in) :: idx(:),totxch,totsnd, totrcv + + ! locals + integer :: np, me, nesd, nerv,& + & proc_to_comm, p2ptag, p2pstat(mpi_status_size),& + & iret, err_act, i, idx_pt, totsnd_, totrcv_,& + & snd_pt, rcv_pt, pnti, data_, n + integer, allocatable, dimension(:) :: bsdidx, brvidx,& + & sdsz, rvsz, prcid, rvhd, sdhd + integer :: int_err(5) + logical :: swap_mpi, swap_sync, swap_send, swap_recv,& + & albf,do_send,do_recv + logical, parameter :: usersend=.false. + + integer, pointer, dimension(:) :: sndbuf, rcvbuf + character(len=20) :: name + + info = 0 + name='psi_swap_tran' + call psb_erractionsave(err_act) + + call psb_info(ictxt,me,np) + if (np == -1) then + info = 2010 + 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),& @@ -642,9 +722,9 @@ subroutine psi_iswaptranv(flag,beta,y,desc_a,work,info,data) snd_pt = 1 rcv_pt = 1 do i=1, totxch - proc_to_comm = d_idx(pnti+psb_proc_id_) - nerv = d_idx(pnti+psb_n_elem_recv_) - nesd = d_idx(pnti+nerv+psb_n_elem_send_) + proc_to_comm = idx(pnti+psb_proc_id_) + nerv = idx(pnti+psb_n_elem_recv_) + nesd = idx(pnti+nerv+psb_n_elem_send_) call psb_get_rank(prcid(proc_to_comm),ictxt,proc_to_comm) brvidx(proc_to_comm) = rcv_pt @@ -668,15 +748,14 @@ subroutine psi_iswaptranv(flag,beta,y,desc_a,work,info,data) end if - - idxr = max(idxr,1) - idxs = max(idxs,1) - if((idxr+idxs) < size(work)) then - sndbuf => work(1:idxs) - rcvbuf => work(idxs+1:idxs+idxr) + 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(idxs),rcvbuf(idxr), stat=info) + allocate(sndbuf(totsnd_),rcvbuf(totrcv_), stat=info) if(info /= 0) then call psb_errpush(4000,name) goto 9999 @@ -692,12 +771,12 @@ subroutine psi_iswaptranv(flag,beta,y,desc_a,work,info,data) snd_pt = 1 rcv_pt = 1 do i=1, totxch - proc_to_comm = d_idx(pnti+psb_proc_id_) - nerv = d_idx(pnti+psb_n_elem_recv_) - nesd = d_idx(pnti+nerv+psb_n_elem_send_) + 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_gth(nerv,d_idx(idx_pt:idx_pt+nerv-1),& + call psi_gth(nerv,idx(idx_pt:idx_pt+nerv-1),& & y,rcvbuf(rcv_pt:rcv_pt+nerv-1)) rcv_pt = rcv_pt + nerv @@ -727,9 +806,9 @@ subroutine psi_iswaptranv(flag,beta,y,desc_a,work,info,data) snd_pt = 1 rcv_pt = 1 do i=1, totxch - proc_to_comm = d_idx(pnti+psb_proc_id_) - nerv = d_idx(pnti+psb_n_elem_recv_) - nesd = d_idx(pnti+nerv+psb_n_elem_send_) + 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 (nerv>0) call psb_snd(ictxt,& @@ -741,6 +820,11 @@ subroutine psi_iswaptranv(flag,beta,y,desc_a,work,info,data) & sndbuf(snd_pt:snd_pt+nesd-1), proc_to_comm) if (nerv>0) call psb_snd(ictxt,& & rcvbuf(rcv_pt:rcv_pt+nerv-1), proc_to_comm) + else if (proc_to_comm == me) then + if (nesd /= nerv) then + write(0,*) 'Fatal error in swaptran: mismatch on self sendf',nerv,nesd + end if + sndbuf(snd_pt:snd_pt+nesd-1) = rcvbuf(rcv_pt:rcv_pt+nerv-1) end if rcv_pt = rcv_pt + nerv snd_pt = snd_pt + nesd @@ -756,11 +840,11 @@ subroutine psi_iswaptranv(flag,beta,y,desc_a,work,info,data) snd_pt = 1 rcv_pt = 1 do i=1, totxch - proc_to_comm = d_idx(pnti+psb_proc_id_) - nerv = d_idx(pnti+psb_n_elem_recv_) - nesd = d_idx(pnti+nerv+psb_n_elem_send_) + proc_to_comm = idx(pnti+psb_proc_id_) + nerv = idx(pnti+psb_n_elem_recv_) + nesd = idx(pnti+nerv+psb_n_elem_send_) call psb_get_rank(prcid(i),ictxt,proc_to_comm) - if (nesd>0) then + if ((nesd>0).and.(proc_to_comm/=me)) then p2ptag = krecvid(ictxt,proc_to_comm,me) call mpi_irecv(sndbuf(snd_pt),nesd,& & mpi_integer,prcid(i),& @@ -779,11 +863,11 @@ subroutine psi_iswaptranv(flag,beta,y,desc_a,work,info,data) snd_pt = 1 rcv_pt = 1 do i=1, totxch - proc_to_comm = d_idx(pnti+psb_proc_id_) - nerv = d_idx(pnti+psb_n_elem_recv_) - nesd = d_idx(pnti+nerv+psb_n_elem_send_) + 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) then + if ((nerv>0).and.(proc_to_comm/=me)) then p2ptag=ksendid(ictxt,proc_to_comm,me) if (usersend) then call mpi_rsend(rcvbuf(rcv_pt),nerv,& @@ -811,11 +895,11 @@ subroutine psi_iswaptranv(flag,beta,y,desc_a,work,info,data) pnti = 1 do i=1, totxch - proc_to_comm = d_idx(pnti+psb_proc_id_) - nerv = d_idx(pnti+psb_n_elem_recv_) - nesd = d_idx(pnti+nerv+psb_n_elem_send_) + proc_to_comm = idx(pnti+psb_proc_id_) + nerv = idx(pnti+psb_n_elem_recv_) + nesd = idx(pnti+nerv+psb_n_elem_send_) p2ptag = krecvid(ictxt,proc_to_comm,me) - + if ((proc_to_comm /= me).and.(nesd>0)) then call mpi_wait(rvhd(i),p2pstat,iret) if(iret /= mpi_success) then @@ -824,6 +908,11 @@ subroutine psi_iswaptranv(flag,beta,y,desc_a,work,info,data) call psb_errpush(info,name,i_err=int_err) goto 9999 end if + else if (proc_to_comm == me) then + if (nesd /= nerv) then + write(0,*) 'Fatal error in swaptran: mismatch on self sendf',nerv,nesd + end if + sndbuf(snd_pt:snd_pt+nesd-1) = rcvbuf(rcv_pt:rcv_pt+nerv-1) end if pnti = pnti + nerv + nesd + 3 end do @@ -835,9 +924,9 @@ subroutine psi_iswaptranv(flag,beta,y,desc_a,work,info,data) snd_pt = 1 rcv_pt = 1 do i=1, totxch - proc_to_comm = d_idx(pnti+psb_proc_id_) - nerv = d_idx(pnti+psb_n_elem_recv_) - nesd = d_idx(pnti+nerv+psb_n_elem_send_) + 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_snd(ictxt,& & rcvbuf(rcv_pt:rcv_pt+nerv-1), proc_to_comm) rcv_pt = rcv_pt + nerv @@ -852,9 +941,9 @@ subroutine psi_iswaptranv(flag,beta,y,desc_a,work,info,data) snd_pt = 1 rcv_pt = 1 do i=1, totxch - proc_to_comm = d_idx(pnti+psb_proc_id_) - nerv = d_idx(pnti+psb_n_elem_recv_) - nesd = d_idx(pnti+nerv+psb_n_elem_send_) + 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_rcv(ictxt,& & sndbuf(snd_pt:snd_pt+nesd-1), proc_to_comm) rcv_pt = rcv_pt + nerv @@ -872,11 +961,11 @@ subroutine psi_iswaptranv(flag,beta,y,desc_a,work,info,data) snd_pt = 1 rcv_pt = 1 do i=1, totxch - proc_to_comm = d_idx(pnti+psb_proc_id_) - nerv = d_idx(pnti+psb_n_elem_recv_) - nesd = d_idx(pnti+nerv+psb_n_elem_send_) + 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+nerv+psb_n_elem_send_ - call psi_sct(nesd,d_idx(idx_pt:idx_pt+nesd-1),& + call psi_sct(nesd,idx(idx_pt:idx_pt+nesd-1),& & sndbuf(snd_pt:snd_pt+nesd-1),beta,y) rcv_pt = rcv_pt + nerv snd_pt = snd_pt + nesd @@ -912,4 +1001,4 @@ subroutine psi_iswaptranv(flag,beta,y,desc_a,work,info,data) return end if return -end subroutine psi_iswaptranv +end subroutine psi_itranidxv diff --git a/base/internals/psi_zswapdata.F90 b/base/internals/psi_zswapdata.F90 index d442ab18..2341af0f 100644 --- a/base/internals/psi_zswapdata.F90 +++ b/base/internals/psi_zswapdata.F90 @@ -40,7 +40,7 @@ ! D real(kind(1.d0)) ! Z complex(kind(1.d0)) ! 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))); +! 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 @@ -86,7 +86,6 @@ subroutine psi_zswapdatam(flag,n,beta,y,desc_a,work,info,data) use psb_error_mod use psb_descriptor_type use psb_penv_mod -!!$ use psi_gthsct_mod #ifdef MPI_MOD use mpi #endif @@ -103,19 +102,8 @@ subroutine psi_zswapdatam(flag,n,beta,y,desc_a,work,info,data) integer, optional :: data ! locals - integer :: ictxt, np, me, nesd, nerv,& - & proc_to_comm, p2ptag, icomm, p2pstat(mpi_status_size),& - & idxs, idxr, iret, err_act, totxch, i, idx_pt,& - & snd_pt, rcv_pt, pnti, data_ - integer, allocatable, dimension(:) :: bsdidx, brvidx,& - & sdsz, rvsz, prcid, rvhd, sdhd + integer :: ictxt, np, me, icomm, idxs, idxr, totxch, data_, err_act integer, pointer :: d_idx(:) - integer :: int_err(5) - logical :: swap_mpi, swap_sync, swap_send, swap_recv,& - & albf,do_send,do_recv - logical, parameter :: usersend=.false. - - complex(kind(1.d0)), pointer, dimension(:) :: sndbuf, rcvbuf character(len=20) :: name info = 0 @@ -124,7 +112,6 @@ subroutine psi_zswapdatam(flag,n,beta,y,desc_a,work,info,data) ictxt = psb_cd_get_context(desc_a) icomm = psb_cd_get_mpic(desc_a) - call psb_info(ictxt,me,np) if (np == -1) then info = 2010 @@ -138,13 +125,6 @@ subroutine psi_zswapdatam(flag,n,beta,y,desc_a,work,info,data) 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 - if(present(data)) then data_ = data else @@ -157,8 +137,78 @@ subroutine psi_zswapdatam(flag,n,beta,y,desc_a,work,info,data) goto 9999 end if - idxr = idxr * n - idxs = idxs * n + + call psi_swapdata(ictxt,icomm,flag,n,beta,y,d_idx,totxch,idxs,idxr,work,info) + if (info /= 0) goto 9999 + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + if (err_act == psb_act_abort_) then + call psb_error(ictxt) + return + end if + return +end subroutine psi_zswapdatam + +subroutine psi_zswapidxm(ictxt,icomm,flag,n,beta,y,idx,totxch,totsnd,totrcv,work,info) + + use psi_mod, psb_protect_name => psi_zswapidxm + use psb_error_mod + use psb_descriptor_type + use psb_penv_mod +#ifdef MPI_MOD + use mpi +#endif + implicit none +#ifdef MPI_H + include 'mpif.h' +#endif + + integer, intent(in) :: ictxt,icomm,flag,n + integer, intent(out) :: info + complex(kind(1.d0)) :: y(:,:), beta + complex(kind(1.d0)), target :: work(:) + integer, intent(in) :: idx(:),totxch,totsnd, totrcv + + ! locals + integer :: np, me, nesd, nerv,& + & proc_to_comm, p2ptag, p2pstat(mpi_status_size),& + & iret, err_act, i, idx_pt, totsnd_, totrcv_,& + & snd_pt, rcv_pt, pnti, data_ + integer, allocatable, dimension(:) :: bsdidx, brvidx,& + & sdsz, rvsz, prcid, rvhd, sdhd + integer :: int_err(5) + logical :: swap_mpi, swap_sync, swap_send, swap_recv,& + & albf,do_send,do_recv + logical, parameter :: usersend=.false. + + complex(kind(1.d0)), pointer, dimension(:) :: sndbuf, rcvbuf + character(len=20) :: name + + info = 0 + name='psi_swap_data' + call psb_erractionsave(err_act) + + call psb_info(ictxt,me,np) + if (np == -1) then + info = 2010 + 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),& @@ -180,9 +230,9 @@ subroutine psi_zswapdatam(flag,n,beta,y,desc_a,work,info,data) snd_pt = 1 rcv_pt = 1 do i=1, totxch - proc_to_comm = d_idx(pnti+psb_proc_id_) - nerv = d_idx(pnti+psb_n_elem_recv_) - nesd = d_idx(pnti+nerv+psb_n_elem_send_) + proc_to_comm = idx(pnti+psb_proc_id_) + nerv = idx(pnti+psb_n_elem_recv_) + nesd = idx(pnti+nerv+psb_n_elem_send_) call psb_get_rank(prcid(proc_to_comm),ictxt,proc_to_comm) brvidx(proc_to_comm) = rcv_pt @@ -190,7 +240,7 @@ subroutine psi_zswapdatam(flag,n,beta,y,desc_a,work,info,data) 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 @@ -205,14 +255,14 @@ subroutine psi_zswapdatam(flag,n,beta,y,desc_a,work,info,data) end if end if - idxr = max(idxr,1) - idxs = max(idxs,1) - if((idxr+idxs) < size(work)) then - sndbuf => work(1:idxs) - rcvbuf => work(idxs+1:idxs+idxr) + 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(idxs),rcvbuf(idxr), stat=info) + allocate(sndbuf(totsnd_),rcvbuf(totrcv_), stat=info) if(info /= 0) then call psb_errpush(4000,name) goto 9999 @@ -226,10 +276,10 @@ subroutine psi_zswapdatam(flag,n,beta,y,desc_a,work,info,data) pnti = 1 snd_pt = 1 do i=1, totxch - nerv = d_idx(pnti+psb_n_elem_recv_) - nesd = d_idx(pnti+nerv+psb_n_elem_send_) + 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,d_idx(idx_pt:idx_pt+nesd-1),& + 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 @@ -260,10 +310,10 @@ subroutine psi_zswapdatam(flag,n,beta,y,desc_a,work,info,data) snd_pt = 1 rcv_pt = 1 do i=1, totxch - proc_to_comm = d_idx(pnti+psb_proc_id_) - nerv = d_idx(pnti+psb_n_elem_recv_) - nesd = d_idx(pnti+nerv+psb_n_elem_send_) - + 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) @@ -274,8 +324,13 @@ subroutine psi_zswapdatam(flag,n,beta,y,desc_a,work,info,data) & 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(0,*) 'Fatal error in swapdata: mismatch on self sendf',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 @@ -286,21 +341,21 @@ subroutine psi_zswapdatam(flag,n,beta,y,desc_a,work,info,data) 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 = d_idx(pnti+psb_proc_id_) - nerv = d_idx(pnti+psb_n_elem_recv_) - nesd = d_idx(pnti+nerv+psb_n_elem_send_) + proc_to_comm = idx(pnti+psb_proc_id_) + nerv = idx(pnti+psb_n_elem_recv_) + nesd = idx(pnti+nerv+psb_n_elem_send_) call psb_get_rank(prcid(i),ictxt,proc_to_comm) - if (nerv > 0) then + if ((nerv>0).and.(proc_to_comm/=me)) then p2ptag = krecvid(ictxt,proc_to_comm,me) call mpi_irecv(rcvbuf(rcv_pt),n*nerv,& & mpi_double_complex,prcid(i),& & p2ptag, icomm,rvhd(i),iret) - endif + end if rcv_pt = rcv_pt + n*nerv snd_pt = snd_pt + n*nesd pnti = pnti + nerv + nesd + 3 @@ -315,12 +370,12 @@ subroutine psi_zswapdatam(flag,n,beta,y,desc_a,work,info,data) snd_pt = 1 rcv_pt = 1 do i=1, totxch - proc_to_comm = d_idx(pnti+psb_proc_id_) - nerv = d_idx(pnti+psb_n_elem_recv_) - nesd = d_idx(pnti+nerv+psb_n_elem_send_) + proc_to_comm = idx(pnti+psb_proc_id_) + nerv = idx(pnti+psb_n_elem_recv_) + nesd = idx(pnti+nerv+psb_n_elem_send_) p2ptag=ksendid(ictxt,proc_to_comm,me) - if (nesd>0) then + if ((nesd>0).and.(proc_to_comm/=me)) then if (usersend) then call mpi_rsend(sndbuf(snd_pt),n*nesd,& & mpi_double_complex,prcid(i),& @@ -330,7 +385,7 @@ subroutine psi_zswapdatam(flag,n,beta,y,desc_a,work,info,data) & mpi_double_complex,prcid(i),& & p2ptag,icomm,iret) end if - + if(iret /= mpi_success) then int_err(1) = iret info=400 @@ -348,9 +403,9 @@ subroutine psi_zswapdatam(flag,n,beta,y,desc_a,work,info,data) pnti = 1 do i=1, totxch - proc_to_comm = d_idx(pnti+psb_proc_id_) - nerv = d_idx(pnti+psb_n_elem_recv_) - nesd = d_idx(pnti+nerv+psb_n_elem_send_) + proc_to_comm = idx(pnti+psb_proc_id_) + nerv = idx(pnti+psb_n_elem_recv_) + nesd = idx(pnti+nerv+psb_n_elem_send_) p2ptag = krecvid(ictxt,proc_to_comm,me) @@ -362,10 +417,15 @@ subroutine psi_zswapdatam(flag,n,beta,y,desc_a,work,info,data) call psb_errpush(info,name,i_err=int_err) goto 9999 end if + else if (proc_to_comm == me) then + if (nesd /= nerv) then + write(0,*) 'Fatal error in swapdata: mismatch on self sendf',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 @@ -374,9 +434,9 @@ subroutine psi_zswapdatam(flag,n,beta,y,desc_a,work,info,data) snd_pt = 1 rcv_pt = 1 do i=1, totxch - proc_to_comm = d_idx(pnti+psb_proc_id_) - nerv = d_idx(pnti+psb_n_elem_recv_) - nesd = d_idx(pnti+nerv+psb_n_elem_send_) + 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) @@ -393,9 +453,9 @@ subroutine psi_zswapdatam(flag,n,beta,y,desc_a,work,info,data) snd_pt = 1 rcv_pt = 1 do i=1, totxch - proc_to_comm = d_idx(pnti+psb_proc_id_) - nerv = d_idx(pnti+psb_n_elem_recv_) - nesd = d_idx(pnti+nerv+psb_n_elem_send_) + 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 @@ -414,11 +474,11 @@ subroutine psi_zswapdatam(flag,n,beta,y,desc_a,work,info,data) snd_pt = 1 rcv_pt = 1 do i=1, totxch - proc_to_comm = d_idx(pnti+psb_proc_id_) - nerv = d_idx(pnti+psb_n_elem_recv_) - nesd = d_idx(pnti+nerv+psb_n_elem_send_) + 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,d_idx(idx_pt:idx_pt+nerv-1),& + 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 @@ -454,39 +514,8 @@ subroutine psi_zswapdatam(flag,n,beta,y,desc_a,work,info,data) return end if return -end subroutine psi_zswapdatam - +end subroutine psi_zswapidxm -!!$ -!!$ Parallel Sparse BLAS v2.0 -!!$ (C) Copyright 2006 Salvatore Filippone University of Rome Tor Vergata -!!$ Alfredo Buttari University of Rome Tor Vergata -!!$ -!!$ 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_zswapdatav @@ -544,7 +573,6 @@ subroutine psi_zswapdatav(flag,beta,y,desc_a,work,info,data) use psb_error_mod use psb_descriptor_type use psb_penv_mod -!!$ use psi_gthsct_mod #ifdef MPI_MOD use mpi #endif @@ -552,7 +580,7 @@ subroutine psi_zswapdatav(flag,beta,y,desc_a,work,info,data) #ifdef MPI_H include 'mpif.h' #endif - + integer, intent(in) :: flag integer, intent(out) :: info complex(kind(1.d0)) :: y(:), beta @@ -561,20 +589,8 @@ subroutine psi_zswapdatav(flag,beta,y,desc_a,work,info,data) integer, optional :: data ! locals - integer :: ictxt, np, me, nesd, nerv,& - & proc_to_comm, p2ptag, icomm, p2pstat(mpi_status_size),& - & idxs, idxr, iret, err_act, totxch, i, & - & idx_pt, snd_pt, rcv_pt, n, pnti, data_ - - integer, allocatable, dimension(:) :: bsdidx, brvidx,& - & sdsz, rvsz, prcid, rvhd, sdhd + integer :: ictxt, np, me, icomm, idxs, idxr, totxch, data_, err_act integer, pointer :: d_idx(:) - integer :: int_err(5) - logical :: swap_mpi, swap_sync, swap_send, swap_recv,& - & albf,do_send,do_recv - logical, parameter :: usersend=.false. - - complex(kind(1.d0)), pointer, dimension(:) :: sndbuf, rcvbuf character(len=20) :: name info = 0 @@ -595,16 +611,7 @@ subroutine psi_zswapdatav(flag,beta,y,desc_a,work,info,data) goto 9999 endif - icomm = desc_a%matrix_data(psb_mpi_c_) - 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 + icomm = psb_cd_get_mpic(desc_a) if(present(data)) then data_ = data @@ -618,8 +625,80 @@ subroutine psi_zswapdatav(flag,beta,y,desc_a,work,info,data) goto 9999 end if - idxr = idxr * n - idxs = idxs * n + call psi_swapdata(ictxt,icomm,flag,beta,y,d_idx,totxch,idxs,idxr,work,info) + if (info /= 0) goto 9999 + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + if (err_act == psb_act_abort_) then + call psb_error(ictxt) + return + end if + return +end subroutine psi_zswapdatav + + +subroutine psi_zswapidxv(ictxt,icomm,flag,beta,y,idx,totxch,totsnd,totrcv,work,info) + + use psi_mod, psb_protect_name => psi_zswapidxv + use psb_error_mod + use psb_descriptor_type + use psb_penv_mod +#ifdef MPI_MOD + use mpi +#endif + implicit none +#ifdef MPI_H + include 'mpif.h' +#endif + + integer, intent(in) :: ictxt,icomm,flag + integer, intent(out) :: info + complex(kind(1.d0)) :: y(:), beta + complex(kind(1.d0)), target :: work(:) + integer, intent(in) :: idx(:),totxch,totsnd, totrcv + + ! locals + integer :: np, me, nesd, nerv,& + & proc_to_comm, p2ptag, p2pstat(mpi_status_size),& + & iret, err_act, i, totsnd_, totrcv_,& + & idx_pt, snd_pt, rcv_pt, n, pnti, data_ + + integer, allocatable, dimension(:) :: bsdidx, brvidx,& + & sdsz, rvsz, prcid, rvhd, sdhd + integer :: int_err(5) + logical :: swap_mpi, swap_sync, swap_send, swap_recv,& + & albf,do_send,do_recv + logical, parameter :: usersend=.false. + + complex(kind(1.d0)), pointer, dimension(:) :: sndbuf, rcvbuf + character(len=20) :: name + + info = 0 + name='psi_swap_datav' + call psb_erractionsave(err_act) + + call psb_info(ictxt,me,np) + if (np == -1) then + info = 2010 + 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),& @@ -640,9 +719,9 @@ subroutine psi_zswapdatav(flag,beta,y,desc_a,work,info,data) snd_pt = 1 rcv_pt = 1 do i=1, totxch - proc_to_comm = d_idx(pnti+psb_proc_id_) - nerv = d_idx(pnti+psb_n_elem_recv_) - nesd = d_idx(pnti+nerv+psb_n_elem_send_) + proc_to_comm = idx(pnti+psb_proc_id_) + nerv = idx(pnti+psb_n_elem_recv_) + nesd = idx(pnti+nerv+psb_n_elem_send_) call psb_get_rank(prcid(proc_to_comm),ictxt,proc_to_comm) brvidx(proc_to_comm) = rcv_pt @@ -666,14 +745,14 @@ subroutine psi_zswapdatav(flag,beta,y,desc_a,work,info,data) end if - idxr = max(idxr,1) - idxs = max(idxs,1) - if((idxr+idxs) < size(work)) then - sndbuf => work(1:idxs) - rcvbuf => work(idxs+1:idxs+idxr) + 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(idxs),rcvbuf(idxr), stat=info) + allocate(sndbuf(totsnd_),rcvbuf(totrcv_), stat=info) if(info /= 0) then call psb_errpush(4000,name) goto 9999 @@ -688,10 +767,10 @@ subroutine psi_zswapdatav(flag,beta,y,desc_a,work,info,data) pnti = 1 snd_pt = 1 do i=1, totxch - nerv = d_idx(pnti+psb_n_elem_recv_) - nesd = d_idx(pnti+nerv+psb_n_elem_send_) + 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,d_idx(idx_pt:idx_pt+nesd-1),& + 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 @@ -719,9 +798,9 @@ subroutine psi_zswapdatav(flag,beta,y,desc_a,work,info,data) snd_pt = 1 rcv_pt = 1 do i=1, totxch - proc_to_comm = d_idx(pnti+psb_proc_id_) - nerv = d_idx(pnti+psb_n_elem_recv_) - nesd = d_idx(pnti+nerv+psb_n_elem_send_) + 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,& @@ -733,6 +812,11 @@ subroutine psi_zswapdatav(flag,beta,y,desc_a,work,info,data) & 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(0,*) 'Fatal error in swapdata: mismatch on self sendf',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 @@ -747,16 +831,17 @@ subroutine psi_zswapdatav(flag,beta,y,desc_a,work,info,data) snd_pt = 1 rcv_pt = 1 do i=1, totxch - proc_to_comm = d_idx(pnti+psb_proc_id_) - nerv = d_idx(pnti+psb_n_elem_recv_) - nesd = d_idx(pnti+nerv+psb_n_elem_send_) + proc_to_comm = idx(pnti+psb_proc_id_) + nerv = idx(pnti+psb_n_elem_recv_) + nesd = idx(pnti+nerv+psb_n_elem_send_) + call psb_get_rank(prcid(i),ictxt,proc_to_comm) - if (nerv>0) then + if ((nerv>0).and.(proc_to_comm /= me)) then p2ptag = krecvid(ictxt,proc_to_comm,me) call mpi_irecv(rcvbuf(rcv_pt),nerv,& & mpi_double_complex,prcid(i),& & p2ptag, icomm,rvhd(i),iret) - endif + end if rcv_pt = rcv_pt + nerv snd_pt = snd_pt + nesd pnti = pnti + nerv + nesd + 3 @@ -770,13 +855,13 @@ subroutine psi_zswapdatav(flag,beta,y,desc_a,work,info,data) snd_pt = 1 rcv_pt = 1 do i=1, totxch - proc_to_comm = d_idx(pnti+psb_proc_id_) - nerv = d_idx(pnti+psb_n_elem_recv_) - nesd = d_idx(pnti+nerv+psb_n_elem_send_) + proc_to_comm = idx(pnti+psb_proc_id_) + nerv = idx(pnti+psb_n_elem_recv_) + nesd = idx(pnti+nerv+psb_n_elem_send_) p2ptag=ksendid(ictxt,proc_to_comm,me) - if (nesd>0) then + if ((nesd>0).and.(proc_to_comm /= me)) then if (usersend) then call mpi_rsend(sndbuf(snd_pt),nesd,& & mpi_double_complex,prcid(i),& @@ -802,9 +887,9 @@ subroutine psi_zswapdatav(flag,beta,y,desc_a,work,info,data) pnti = 1 do i=1, totxch - proc_to_comm = d_idx(pnti+psb_proc_id_) - nerv = d_idx(pnti+psb_n_elem_recv_) - nesd = d_idx(pnti+nerv+psb_n_elem_send_) + proc_to_comm = idx(pnti+psb_proc_id_) + nerv = idx(pnti+psb_n_elem_recv_) + nesd = idx(pnti+nerv+psb_n_elem_send_) p2ptag = krecvid(ictxt,proc_to_comm,me) @@ -816,6 +901,11 @@ subroutine psi_zswapdatav(flag,beta,y,desc_a,work,info,data) call psb_errpush(info,name,i_err=int_err) goto 9999 end if + else if (proc_to_comm == me) then + if (nesd /= nerv) then + write(0,*) 'Fatal error in swapdata: mismatch on self sendf',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 @@ -827,9 +917,9 @@ subroutine psi_zswapdatav(flag,beta,y,desc_a,work,info,data) snd_pt = 1 rcv_pt = 1 do i=1, totxch - proc_to_comm = d_idx(pnti+psb_proc_id_) - nerv = d_idx(pnti+psb_n_elem_recv_) - nesd = d_idx(pnti+nerv+psb_n_elem_send_) + 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 @@ -843,9 +933,9 @@ subroutine psi_zswapdatav(flag,beta,y,desc_a,work,info,data) snd_pt = 1 rcv_pt = 1 do i=1, totxch - proc_to_comm = d_idx(pnti+psb_proc_id_) - nerv = d_idx(pnti+psb_n_elem_recv_) - nesd = d_idx(pnti+nerv+psb_n_elem_send_) + 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 @@ -855,18 +945,17 @@ subroutine psi_zswapdatav(flag,beta,y,desc_a,work,info,data) end if - if (do_recv) then pnti = 1 snd_pt = 1 rcv_pt = 1 do i=1, totxch - proc_to_comm = d_idx(pnti+psb_proc_id_) - nerv = d_idx(pnti+psb_n_elem_recv_) - nesd = d_idx(pnti+nerv+psb_n_elem_send_) + 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,d_idx(idx_pt:idx_pt+nerv-1),& + 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 @@ -901,4 +990,4 @@ subroutine psi_zswapdatav(flag,beta,y,desc_a,work,info,data) return end if return -end subroutine psi_zswapdatav +end subroutine psi_zswapidxv diff --git a/base/internals/psi_zswaptran.F90 b/base/internals/psi_zswaptran.F90 index 7a892608..121461be 100644 --- a/base/internals/psi_zswaptran.F90 +++ b/base/internals/psi_zswaptran.F90 @@ -90,7 +90,6 @@ subroutine psi_zswaptranm(flag,n,beta,y,desc_a,work,info,data) use psb_error_mod use psb_descriptor_type use psb_penv_mod -!!$ use psi_gthsct_mod #ifdef MPI_MOD use mpi #endif @@ -107,19 +106,9 @@ subroutine psi_zswaptranm(flag,n,beta,y,desc_a,work,info,data) integer, optional :: data ! locals - integer :: ictxt, np, me, nesd, nerv,& - & proc_to_comm, p2ptag, icomm, p2pstat(mpi_status_size),& - & idxs, idxr, iret, err_act, totxch, i, idx_pt,& - & snd_pt, rcv_pt, pnti, data_ - integer, allocatable, dimension(:) :: bsdidx, brvidx,& - & sdsz, rvsz, prcid, rvhd, sdhd + integer :: ictxt, np, me, icomm, idxs, idxr, err_act, totxch, data_ integer, pointer :: d_idx(:) integer :: int_err(5) - logical :: swap_mpi, swap_sync, swap_send, swap_recv,& - & albf,do_send,do_recv - logical, parameter :: usersend=.false. - - complex(kind(1.d0)), pointer, dimension(:) :: sndbuf, rcvbuf character(len=20) :: name info = 0 @@ -142,13 +131,6 @@ subroutine psi_zswaptranm(flag,n,beta,y,desc_a,work,info,data) 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 - if(present(data)) then data_ = data else @@ -161,8 +143,77 @@ subroutine psi_zswaptranm(flag,n,beta,y,desc_a,work,info,data) goto 9999 end if - idxr = idxr * n - idxs = idxs * n + call psi_swaptran(ictxt,icomm,flag,n,beta,y,d_idx,totxch,idxs,idxr,work,info) + if (info /= 0) goto 9999 + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + if (err_act == psb_act_abort_) then + call psb_error(ictxt) + return + end if + return +end subroutine psi_zswaptranm + +subroutine psi_ztranidxm(ictxt,icomm,flag,n,beta,y,idx,totxch,totsnd,totrcv,work,info) + + use psi_mod, psb_protect_name => psi_ztranidxm + use psb_error_mod + use psb_descriptor_type + use psb_penv_mod +#ifdef MPI_MOD + use mpi +#endif + implicit none +#ifdef MPI_H + include 'mpif.h' +#endif + + integer, intent(in) :: ictxt,icomm,flag,n + integer, intent(out) :: info + complex(kind(1.d0)) :: y(:,:), beta + complex(kind(1.d0)), target :: work(:) + integer, intent(in) :: idx(:),totxch,totsnd, totrcv + + ! locals + integer :: np, me, nesd, nerv,& + & proc_to_comm, p2ptag, p2pstat(mpi_status_size),& + & iret, err_act, i, idx_pt, totsnd_, totrcv_,& + & snd_pt, rcv_pt, pnti, data_ + integer, allocatable, dimension(:) :: bsdidx, brvidx,& + & sdsz, rvsz, prcid, rvhd, sdhd + integer :: int_err(5) + logical :: swap_mpi, swap_sync, swap_send, swap_recv,& + & albf,do_send,do_recv + logical, parameter :: usersend=.false. + + complex(kind(1.d0)), pointer, dimension(:) :: sndbuf, rcvbuf + character(len=20) :: name + + info = 0 + name='psi_swap_tran' + call psb_erractionsave(err_act) + + call psb_info(ictxt,me,np) + if (np == -1) then + info = 2010 + 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),& @@ -183,9 +234,9 @@ subroutine psi_zswaptranm(flag,n,beta,y,desc_a,work,info,data) snd_pt = 1 rcv_pt = 1 do i=1, totxch - proc_to_comm = d_idx(pnti+psb_proc_id_) - nerv = d_idx(pnti+psb_n_elem_recv_) - nesd = d_idx(pnti+nerv+psb_n_elem_send_) + proc_to_comm = idx(pnti+psb_proc_id_) + nerv = idx(pnti+psb_n_elem_recv_) + nesd = idx(pnti+nerv+psb_n_elem_send_) call psb_get_rank(prcid(proc_to_comm),ictxt,proc_to_comm) @@ -210,14 +261,14 @@ subroutine psi_zswaptranm(flag,n,beta,y,desc_a,work,info,data) end if - idxr = max(idxr,1) - idxs = max(idxs,1) - if((idxr+idxs) < size(work)) then - sndbuf => work(1:idxs) - rcvbuf => work(idxs+1:idxs+idxr) + 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(idxs),rcvbuf(idxr), stat=info) + allocate(sndbuf(totsnd_),rcvbuf(totrcv_), stat=info) if(info /= 0) then call psb_errpush(4000,name) goto 9999 @@ -232,12 +283,12 @@ subroutine psi_zswaptranm(flag,n,beta,y,desc_a,work,info,data) snd_pt = 1 rcv_pt = 1 do i=1, totxch - proc_to_comm = d_idx(pnti+psb_proc_id_) - nerv = d_idx(pnti+psb_n_elem_recv_) - nesd = d_idx(pnti+nerv+psb_n_elem_send_) + 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_gth(nerv,n,d_idx(idx_pt:idx_pt+nerv-1),& + call psi_gth(nerv,n,idx(idx_pt:idx_pt+nerv-1),& & y,rcvbuf(rcv_pt:rcv_pt+n*nerv-1)) rcv_pt = rcv_pt + n*nerv @@ -269,9 +320,9 @@ subroutine psi_zswaptranm(flag,n,beta,y,desc_a,work,info,data) snd_pt = 1 rcv_pt = 1 do i=1, totxch - proc_to_comm = d_idx(pnti+psb_proc_id_) - nerv = d_idx(pnti+psb_n_elem_recv_) - nesd = d_idx(pnti+nerv+psb_n_elem_send_) + 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 (nerv>0) call psb_snd(ictxt,& @@ -283,6 +334,11 @@ subroutine psi_zswaptranm(flag,n,beta,y,desc_a,work,info,data) & sndbuf(snd_pt:snd_pt+n*nesd-1), proc_to_comm) if (nerv>0) call psb_snd(ictxt,& & rcvbuf(rcv_pt:rcv_pt+n*nerv-1), proc_to_comm) + else if (proc_to_comm == me) then + if (nesd /= nerv) then + write(0,*) 'Fatal error in swaptran: mismatch on self sendf',nerv,nesd + end if + sndbuf(snd_pt:snd_pt+n*nesd-1) = rcvbuf(rcv_pt:rcv_pt+n*nerv-1) end if rcv_pt = rcv_pt + n*nerv snd_pt = snd_pt + n*nesd @@ -298,11 +354,11 @@ subroutine psi_zswaptranm(flag,n,beta,y,desc_a,work,info,data) snd_pt = 1 rcv_pt = 1 do i=1, totxch - proc_to_comm = d_idx(pnti+psb_proc_id_) - nerv = d_idx(pnti+psb_n_elem_recv_) - nesd = d_idx(pnti+nerv+psb_n_elem_send_) + proc_to_comm = idx(pnti+psb_proc_id_) + nerv = idx(pnti+psb_n_elem_recv_) + nesd = idx(pnti+nerv+psb_n_elem_send_) call psb_get_rank(prcid(i),ictxt,proc_to_comm) - if (nesd>0) then + if ((nesd>0).and.(proc_to_comm/=me)) then p2ptag = krecvid(ictxt,proc_to_comm,me) call mpi_irecv(sndbuf(snd_pt),n*nesd,& & mpi_double_complex,prcid(i),& @@ -321,11 +377,11 @@ subroutine psi_zswaptranm(flag,n,beta,y,desc_a,work,info,data) snd_pt = 1 rcv_pt = 1 do i=1, totxch - proc_to_comm = d_idx(pnti+psb_proc_id_) - nerv = d_idx(pnti+psb_n_elem_recv_) - nesd = d_idx(pnti+nerv+psb_n_elem_send_) + 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) then + if ((nerv>0).and.(proc_to_comm/=me)) then p2ptag=ksendid(ictxt,proc_to_comm,me) if (usersend) then call mpi_rsend(rcvbuf(rcv_pt),n*nerv,& @@ -353,9 +409,9 @@ subroutine psi_zswaptranm(flag,n,beta,y,desc_a,work,info,data) pnti = 1 do i=1, totxch - proc_to_comm = d_idx(pnti+psb_proc_id_) - nerv = d_idx(pnti+psb_n_elem_recv_) - nesd = d_idx(pnti+nerv+psb_n_elem_send_) + proc_to_comm = idx(pnti+psb_proc_id_) + nerv = idx(pnti+psb_n_elem_recv_) + nesd = idx(pnti+nerv+psb_n_elem_send_) p2ptag = krecvid(ictxt,proc_to_comm,me) @@ -367,6 +423,11 @@ subroutine psi_zswaptranm(flag,n,beta,y,desc_a,work,info,data) call psb_errpush(info,name,i_err=int_err) goto 9999 end if + else if (proc_to_comm == me) then + if (nesd /= nerv) then + write(0,*) 'Fatal error in swaptran: mismatch on self sendf',nerv,nesd + end if + sndbuf(snd_pt:snd_pt+n*nesd-1) = rcvbuf(rcv_pt:rcv_pt+n*nerv-1) end if pnti = pnti + nerv + nesd + 3 end do @@ -378,9 +439,9 @@ subroutine psi_zswaptranm(flag,n,beta,y,desc_a,work,info,data) snd_pt = 1 rcv_pt = 1 do i=1, totxch - proc_to_comm = d_idx(pnti+psb_proc_id_) - nerv = d_idx(pnti+psb_n_elem_recv_) - nesd = d_idx(pnti+nerv+psb_n_elem_send_) + 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_snd(ictxt,& & rcvbuf(rcv_pt:rcv_pt+n*nerv-1), proc_to_comm) rcv_pt = rcv_pt + n*nerv @@ -395,9 +456,9 @@ subroutine psi_zswaptranm(flag,n,beta,y,desc_a,work,info,data) snd_pt = 1 rcv_pt = 1 do i=1, totxch - proc_to_comm = d_idx(pnti+psb_proc_id_) - nerv = d_idx(pnti+psb_n_elem_recv_) - nesd = d_idx(pnti+nerv+psb_n_elem_send_) + 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_rcv(ictxt,& & sndbuf(snd_pt:snd_pt+n*nesd-1), proc_to_comm) rcv_pt = rcv_pt + n*nerv @@ -413,11 +474,11 @@ subroutine psi_zswaptranm(flag,n,beta,y,desc_a,work,info,data) snd_pt = 1 rcv_pt = 1 do i=1, totxch - proc_to_comm = d_idx(pnti+psb_proc_id_) - nerv = d_idx(pnti+psb_n_elem_recv_) - nesd = d_idx(pnti+nerv+psb_n_elem_send_) + 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+nerv+psb_n_elem_send_ - call psi_sct(nesd,n,d_idx(idx_pt:idx_pt+nesd-1),& + call psi_sct(nesd,n,idx(idx_pt:idx_pt+nesd-1),& & sndbuf(snd_pt:snd_pt+n*nesd-1),beta,y) rcv_pt = rcv_pt + n*nerv snd_pt = snd_pt + n*nesd @@ -453,39 +514,7 @@ subroutine psi_zswaptranm(flag,n,beta,y,desc_a,work,info,data) return end if return -end subroutine psi_zswaptranm - - -!!$ -!!$ Parallel Sparse BLAS v2.0 -!!$ (C) Copyright 2006 Salvatore Filippone University of Rome Tor Vergata -!!$ Alfredo Buttari University of Rome Tor Vergata -!!$ -!!$ 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. -!!$ -!!$ +end subroutine psi_ztranidxm ! ! ! Subroutine: psi_zswaptranv @@ -547,7 +576,6 @@ subroutine psi_zswaptranv(flag,beta,y,desc_a,work,info,data) use psb_error_mod use psb_descriptor_type use psb_penv_mod -!!$ use psi_gthsct_mod #ifdef MPI_MOD use mpi #endif @@ -564,28 +592,17 @@ subroutine psi_zswaptranv(flag,beta,y,desc_a,work,info,data) integer, optional :: data ! locals - integer :: ictxt, np, me, nesd, nerv,& - & proc_to_comm, p2ptag, icomm, p2pstat(mpi_status_size),& - & idxs, idxr, iret, err_act, totxch, i, & - & idx_pt, snd_pt, rcv_pt, n, pnti, data_ - - integer, allocatable, dimension(:) :: bsdidx, brvidx,& - & sdsz, rvsz, prcid, rvhd, sdhd + integer :: ictxt, np, me, icomm, idxs, idxr, totxch, err_act, data_ integer, pointer :: d_idx(:) integer :: int_err(5) - logical :: swap_mpi, swap_sync, swap_send, swap_recv,& - & albf,do_send,do_recv - logical, parameter :: usersend=.false. - - complex(kind(1.d0)), pointer, dimension(:) :: sndbuf, rcvbuf character(len=20) :: name info = 0 name='psi_swap_tranv' call psb_erractionsave(err_act) - ictxt=psb_cd_get_context(desc_a) - + ictxt = psb_cd_get_context(desc_a) + icomm = psb_cd_get_mpic(desc_a) call psb_info(ictxt,me,np) if (np == -1) then info = 2010 @@ -599,30 +616,91 @@ subroutine psi_zswaptranv(flag,beta,y,desc_a,work,info,data) goto 9999 endif - icomm = desc_a%matrix_data(psb_mpi_c_) - 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 - - if(present(data)) then + if (present(data)) then data_ = data else data_ = psb_comm_halo_ end if - + call psb_cd_get_list(data_,desc_a,d_idx,totxch,idxr,idxs,info) if (info /= 0) then call psb_errpush(4001,name,a_err='psb_cd_get_list') goto 9999 end if - idxr = idxr * n - idxs = idxs * n + call psi_swaptran(ictxt,icomm,flag,beta,y,d_idx,totxch,idxs,idxr,work,info) + if (info /= 0) goto 9999 + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + if (err_act == psb_act_abort_) then + call psb_error(ictxt) + return + end if + return +end subroutine psi_zswaptranv + + + +subroutine psi_ztranidxv(ictxt,icomm,flag,beta,y,idx,totxch,totsnd,totrcv,work,info) + + use psi_mod, psb_protect_name => psi_ztranidxv + use psb_error_mod + use psb_descriptor_type + use psb_penv_mod +#ifdef MPI_MOD + use mpi +#endif + implicit none +#ifdef MPI_H + include 'mpif.h' +#endif + + integer, intent(in) :: ictxt,icomm,flag + integer, intent(out) :: info + complex(kind(1.d0)) :: y(:), beta + complex(kind(1.d0)), target :: work(:) + integer, intent(in) :: idx(:),totxch,totsnd, totrcv + + ! locals + integer :: np, me, nesd, nerv,& + & proc_to_comm, p2ptag, p2pstat(mpi_status_size),& + & iret, err_act, i, idx_pt, totsnd_, totrcv_,& + & snd_pt, rcv_pt, pnti, data_, n + integer, allocatable, dimension(:) :: bsdidx, brvidx,& + & sdsz, rvsz, prcid, rvhd, sdhd + integer :: int_err(5) + logical :: swap_mpi, swap_sync, swap_send, swap_recv,& + & albf,do_send,do_recv + logical, parameter :: usersend=.false. + + complex(kind(1.d0)), pointer, dimension(:) :: sndbuf, rcvbuf + character(len=20) :: name + + info = 0 + name='psi_swap_tran' + call psb_erractionsave(err_act) + + call psb_info(ictxt,me,np) + if (np == -1) then + info = 2010 + 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),& @@ -644,9 +722,9 @@ subroutine psi_zswaptranv(flag,beta,y,desc_a,work,info,data) snd_pt = 1 rcv_pt = 1 do i=1, totxch - proc_to_comm = d_idx(pnti+psb_proc_id_) - nerv = d_idx(pnti+psb_n_elem_recv_) - nesd = d_idx(pnti+nerv+psb_n_elem_send_) + proc_to_comm = idx(pnti+psb_proc_id_) + nerv = idx(pnti+psb_n_elem_recv_) + nesd = idx(pnti+nerv+psb_n_elem_send_) call psb_get_rank(prcid(proc_to_comm),ictxt,proc_to_comm) brvidx(proc_to_comm) = rcv_pt @@ -670,15 +748,14 @@ subroutine psi_zswaptranv(flag,beta,y,desc_a,work,info,data) end if - - idxr = max(idxr,1) - idxs = max(idxs,1) - if((idxr+idxs) < size(work)) then - sndbuf => work(1:idxs) - rcvbuf => work(idxs+1:idxs+idxr) + 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(idxs),rcvbuf(idxr), stat=info) + allocate(sndbuf(totsnd_),rcvbuf(totrcv_), stat=info) if(info /= 0) then call psb_errpush(4000,name) goto 9999 @@ -694,12 +771,12 @@ subroutine psi_zswaptranv(flag,beta,y,desc_a,work,info,data) snd_pt = 1 rcv_pt = 1 do i=1, totxch - proc_to_comm = d_idx(pnti+psb_proc_id_) - nerv = d_idx(pnti+psb_n_elem_recv_) - nesd = d_idx(pnti+nerv+psb_n_elem_send_) + 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_gth(nerv,d_idx(idx_pt:idx_pt+nerv-1),& + call psi_gth(nerv,idx(idx_pt:idx_pt+nerv-1),& & y,rcvbuf(rcv_pt:rcv_pt+nerv-1)) rcv_pt = rcv_pt + nerv @@ -729,9 +806,9 @@ subroutine psi_zswaptranv(flag,beta,y,desc_a,work,info,data) snd_pt = 1 rcv_pt = 1 do i=1, totxch - proc_to_comm = d_idx(pnti+psb_proc_id_) - nerv = d_idx(pnti+psb_n_elem_recv_) - nesd = d_idx(pnti+nerv+psb_n_elem_send_) + 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 (nerv>0) call psb_snd(ictxt,& @@ -743,6 +820,11 @@ subroutine psi_zswaptranv(flag,beta,y,desc_a,work,info,data) & sndbuf(snd_pt:snd_pt+nesd-1), proc_to_comm) if (nerv>0) call psb_snd(ictxt,& & rcvbuf(rcv_pt:rcv_pt+nerv-1), proc_to_comm) + else if (proc_to_comm == me) then + if (nesd /= nerv) then + write(0,*) 'Fatal error in swaptran: mismatch on self sendf',nerv,nesd + end if + sndbuf(snd_pt:snd_pt+nesd-1) = rcvbuf(rcv_pt:rcv_pt+nerv-1) end if rcv_pt = rcv_pt + nerv snd_pt = snd_pt + nesd @@ -758,11 +840,11 @@ subroutine psi_zswaptranv(flag,beta,y,desc_a,work,info,data) snd_pt = 1 rcv_pt = 1 do i=1, totxch - proc_to_comm = d_idx(pnti+psb_proc_id_) - nerv = d_idx(pnti+psb_n_elem_recv_) - nesd = d_idx(pnti+nerv+psb_n_elem_send_) + proc_to_comm = idx(pnti+psb_proc_id_) + nerv = idx(pnti+psb_n_elem_recv_) + nesd = idx(pnti+nerv+psb_n_elem_send_) call psb_get_rank(prcid(i),ictxt,proc_to_comm) - if (nesd>0) then + if ((nesd>0).and.(proc_to_comm/=me)) then p2ptag = krecvid(ictxt,proc_to_comm,me) call mpi_irecv(sndbuf(snd_pt),nesd,& & mpi_double_complex,prcid(i),& @@ -781,11 +863,11 @@ subroutine psi_zswaptranv(flag,beta,y,desc_a,work,info,data) snd_pt = 1 rcv_pt = 1 do i=1, totxch - proc_to_comm = d_idx(pnti+psb_proc_id_) - nerv = d_idx(pnti+psb_n_elem_recv_) - nesd = d_idx(pnti+nerv+psb_n_elem_send_) + 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) then + if ((nerv>0).and.(proc_to_comm/=me)) then p2ptag=ksendid(ictxt,proc_to_comm,me) if (usersend) then call mpi_rsend(rcvbuf(rcv_pt),nerv,& @@ -813,11 +895,11 @@ subroutine psi_zswaptranv(flag,beta,y,desc_a,work,info,data) pnti = 1 do i=1, totxch - proc_to_comm = d_idx(pnti+psb_proc_id_) - nerv = d_idx(pnti+psb_n_elem_recv_) - nesd = d_idx(pnti+nerv+psb_n_elem_send_) + proc_to_comm = idx(pnti+psb_proc_id_) + nerv = idx(pnti+psb_n_elem_recv_) + nesd = idx(pnti+nerv+psb_n_elem_send_) p2ptag = krecvid(ictxt,proc_to_comm,me) - + if ((proc_to_comm /= me).and.(nesd>0)) then call mpi_wait(rvhd(i),p2pstat,iret) if(iret /= mpi_success) then @@ -826,10 +908,15 @@ subroutine psi_zswaptranv(flag,beta,y,desc_a,work,info,data) call psb_errpush(info,name,i_err=int_err) goto 9999 end if + else if (proc_to_comm == me) then + if (nesd /= nerv) then + write(0,*) 'Fatal error in swaptran: mismatch on self sendf',nerv,nesd + end if + sndbuf(snd_pt:snd_pt+nesd-1) = rcvbuf(rcv_pt:rcv_pt+nerv-1) end if pnti = pnti + nerv + nesd + 3 end do - + else if (swap_send) then @@ -837,9 +924,9 @@ subroutine psi_zswaptranv(flag,beta,y,desc_a,work,info,data) snd_pt = 1 rcv_pt = 1 do i=1, totxch - proc_to_comm = d_idx(pnti+psb_proc_id_) - nerv = d_idx(pnti+psb_n_elem_recv_) - nesd = d_idx(pnti+nerv+psb_n_elem_send_) + 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_snd(ictxt,& & rcvbuf(rcv_pt:rcv_pt+nerv-1), proc_to_comm) rcv_pt = rcv_pt + nerv @@ -854,9 +941,9 @@ subroutine psi_zswaptranv(flag,beta,y,desc_a,work,info,data) snd_pt = 1 rcv_pt = 1 do i=1, totxch - proc_to_comm = d_idx(pnti+psb_proc_id_) - nerv = d_idx(pnti+psb_n_elem_recv_) - nesd = d_idx(pnti+nerv+psb_n_elem_send_) + 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_rcv(ictxt,& & sndbuf(snd_pt:snd_pt+nesd-1), proc_to_comm) rcv_pt = rcv_pt + nerv @@ -874,11 +961,11 @@ subroutine psi_zswaptranv(flag,beta,y,desc_a,work,info,data) snd_pt = 1 rcv_pt = 1 do i=1, totxch - proc_to_comm = d_idx(pnti+psb_proc_id_) - nerv = d_idx(pnti+psb_n_elem_recv_) - nesd = d_idx(pnti+nerv+psb_n_elem_send_) + 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+nerv+psb_n_elem_send_ - call psi_sct(nesd,d_idx(idx_pt:idx_pt+nesd-1),& + call psi_sct(nesd,idx(idx_pt:idx_pt+nesd-1),& & sndbuf(snd_pt:snd_pt+nesd-1),beta,y) rcv_pt = rcv_pt + nerv snd_pt = snd_pt + nesd @@ -914,4 +1001,4 @@ subroutine psi_zswaptranv(flag,beta,y,desc_a,work,info,data) return end if return -end subroutine psi_zswaptranv +end subroutine psi_ztranidxv diff --git a/base/modules/Makefile b/base/modules/Makefile index 94492927..33d75436 100644 --- a/base/modules/Makefile +++ b/base/modules/Makefile @@ -3,9 +3,8 @@ include ../../Make.inc MODULES = psb_realloc_mod.o psb_string_mod.o psb_spmat_type.o \ psb_desc_type.o psb_spsb_mod.o psb_sort_mod.o\ psb_serial_mod.o psb_tools_mod.o \ - psb_error_mod.o \ - psb_const_mod.o \ - psb_comm_mod.o psb_psblas_mod.o psi_mod.o \ + psb_error_mod.o psb_const_mod.o psb_inter_desc_type.o \ + psb_comm_mod.o psb_psblas_mod.o psi_serial_mod.o psi_mod.o \ psb_check_mod.o psb_gps_mod.o # psb_methd_mod.o psb_prec_type.o psb_prec_mod.o blacs_env.o @@ -28,16 +27,18 @@ psb_realloc_mod.o : psb_error_mod.o psb_spmat_type.o : psb_realloc_mod.o psb_error_mod.o psb_const_mod.o psb_string_mod.o psb_error_mod.o: psb_const_mod.o psb_penv_mod.o: psb_const_mod.o psb_error_mod.o psb_realloc_mod.o -psi_mod.o: psb_penv_mod.o psb_error_mod.o psb_desc_type.o psb_const_mod.o +psi_serial_mod.o: psb_const_mod.o +psi_mod.o: psb_penv_mod.o psb_error_mod.o psb_desc_type.o psb_const_mod.o psi_serial_mod.o psb_desc_type.o: psb_const_mod.o psb_error_mod.o psb_penv_mod.o +psb_inter_desc_type.o: psb_desc_type.o psb_spmat_type.o psb_error_mod.o psb_serial_mod.o psb_comm_mod.o psb_check_mod.o: psb_desc_type.o -psb_serial_mod: psb_spmat_type.o psb_string_mod.o psb_sort_mod.o +psb_serial_mod.o: psb_spmat_type.o psb_string_mod.o psb_sort_mod.o psi_serial_mod.o psb_sort_mod.o: psb_error_mod.o psb_realloc_mod.o psb_methd_mod.o: psb_serial_mod.o psb_desc_type.o psb_prec_type.o -psb_tools_mod.o: psb_spmat_type.o psb_desc_type.o psi_mod.o psb_gps_mod.o +psb_tools_mod.o: psb_spmat_type.o psb_desc_type.o psi_mod.o psb_gps_mod.o psb_inter_desc_type.o psb_gps_mod.o: psb_realloc_mod.o -$(LIBMOD): $(MODULES:.o:$(.mod)) +psb_base_mod.o: $(MODULES) mpfobjs: (make $(MPFOBJS) F90="$(MPF90)" FC="$(MPF90)" FCOPT="$(F90COPT)") diff --git a/base/modules/psb_base_mod.f90 b/base/modules/psb_base_mod.f90 index 51c88cb5..a06facf8 100644 --- a/base/modules/psb_base_mod.f90 +++ b/base/modules/psb_base_mod.f90 @@ -34,9 +34,11 @@ module psb_base_mod use psb_penv_mod use psb_check_mod use psb_descriptor_type + use psb_inter_descriptor_type use psb_serial_mod use psb_comm_mod use psb_psblas_mod + use psb_gps_mod use psb_tools_mod end module psb_base_mod diff --git a/base/modules/psb_desc_type.f90 b/base/modules/psb_desc_type.f90 index f4e3fdfd..b8309285 100644 --- a/base/modules/psb_desc_type.f90 +++ b/base/modules/psb_desc_type.f90 @@ -55,8 +55,13 @@ module psb_descriptor_type integer, parameter :: psb_no_comm_=-1 integer, parameter :: psb_comm_halo_=1, psb_comm_ovr_=2 integer, parameter :: psb_comm_ext_=3, psb_comm_mov_=4 - integer, parameter :: psb_ovt_xhal_ = 123, psb_ovt_asov_=psb_ovt_xhal_+1 - + ! Types of mapping between descriptors. + integer, parameter :: psb_map_xhal_ = 123 + integer, parameter :: psb_map_asov_ = psb_map_xhal_+1 + integer, parameter :: psb_map_aggr_ = psb_map_asov_+1 + integer, parameter :: psb_map_gen_linear_ = psb_map_aggr_+1 + + integer, parameter :: psb_ovt_xhal_ = psb_map_xhal_, psb_ovt_asov_=psb_map_asov_ ! ! Entries and values in desc%matrix_data ! @@ -81,10 +86,12 @@ module psb_descriptor_type integer, parameter :: psb_desc_bld_=psb_desc_asb_+1 integer, parameter :: psb_desc_repl_=3199 integer, parameter :: psb_desc_upd_=psb_desc_bld_+1 - integer, parameter :: psb_desc_normal_=3299 - integer, parameter :: psb_desc_large_=psb_desc_normal_+1 + ! these two are reserved for descriptors which are + ! "overlap-extensions" of base descriptors. integer, parameter :: psb_cd_ovl_bld_=3399 integer, parameter :: psb_cd_ovl_asb_=psb_cd_ovl_bld_+1 + integer, parameter :: psb_desc_normal_=3299 + integer, parameter :: psb_desc_large_=psb_desc_normal_+1 ! ! Constants for hashing into desc%hashv(:) and desc%glb_lc(:,:) ! @@ -292,18 +299,44 @@ module psb_descriptor_type integer, allocatable :: hashv(:), glb_lc(:,:), ptree(:) integer, allocatable :: lprm(:) integer, allocatable :: idx_space(:) + type(psb_desc_type), pointer :: base_desc => null() end type psb_desc_type interface psb_sizeof module procedure psb_cd_sizeof end interface + + interface psb_is_ok_desc + module procedure psb_is_ok_desc + end interface + + interface psb_is_asb_desc + module procedure psb_is_asb_desc + end interface + + interface psb_is_upd_desc + module procedure psb_is_upd_desc + end interface + + interface psb_is_ovl_desc + module procedure psb_is_ovl_desc + end interface + + interface psb_is_bld_desc + module procedure psb_is_bld_desc + end interface + + interface psb_is_large_desc + module procedure psb_is_large_desc + end interface + integer, private, save :: cd_large_threshold=psb_default_large_threshold contains - + function psb_cd_sizeof(desc) implicit none !....Parameters... @@ -401,6 +434,13 @@ contains end function psb_is_upd_desc + logical function psb_is_ovl_desc(desc) + type(psb_desc_type), intent(in) :: desc + + psb_is_ovl_desc = psb_is_ovl_dec(psb_cd_get_dectype(desc)) + + end function psb_is_ovl_desc + logical function psb_is_asb_desc(desc) type(psb_desc_type), intent(in) :: desc @@ -414,6 +454,7 @@ contains integer :: dectype psb_is_ok_dec = ((dectype == psb_desc_asb_).or.(dectype == psb_desc_bld_).or.& + &(dectype == psb_cd_ovl_asb_).or.(dectype == psb_cd_ovl_bld_).or.& &(dectype == psb_desc_upd_).or.& &(dectype== psb_desc_repl_)) end function psb_is_ok_dec @@ -421,7 +462,7 @@ contains logical function psb_is_bld_dec(dectype) integer :: dectype - psb_is_bld_dec = (dectype == psb_desc_bld_) + psb_is_bld_dec = (dectype == psb_desc_bld_).or.(dectype == psb_cd_ovl_bld_) end function psb_is_bld_dec logical function psb_is_upd_dec(dectype) @@ -436,10 +477,18 @@ contains integer :: dectype psb_is_asb_dec = (dectype == psb_desc_asb_).or.& - & (dectype== psb_desc_repl_) + & (dectype== psb_desc_repl_).or.(dectype == psb_cd_ovl_asb_) end function psb_is_asb_dec + logical function psb_is_ovl_dec(dectype) + integer :: dectype + + psb_is_ovl_dec = (dectype == psb_cd_ovl_bld_).or.& + & (dectype == psb_cd_ovl_asb_) + + end function psb_is_ovl_dec + integer function psb_cd_get_local_rows(desc) type(psb_desc_type), intent(in) :: desc @@ -573,6 +622,61 @@ contains return end subroutine psb_cd_set_bld + subroutine psb_cd_set_ovl_asb(desc,info) + ! + ! Change state of a descriptor into ovl_build. + implicit none + type(psb_desc_type), intent(inout) :: desc + integer :: info + + + if (psb_is_asb_desc(desc)) desc%matrix_data(psb_dec_type_) = psb_cd_ovl_asb_ + + end subroutine psb_cd_set_ovl_asb + + subroutine psb_cd_set_ovl_bld(desc,info) + ! + ! Change state of a descriptor into ovl_build. + implicit none + type(psb_desc_type), intent(inout) :: desc + integer :: info + + call psb_cd_set_bld(desc,info) + if (info == 0) desc%matrix_data(psb_dec_type_) = psb_cd_ovl_bld_ + + end subroutine psb_cd_set_ovl_bld + + + subroutine psb_get_xch_idx(idx,totxch,totsnd,totrcv) + implicit none + integer, intent(in) :: idx(:) + integer, intent(out) :: totxch,totsnd,totrcv + + integer :: ip, nerv, nesd + character(len=20), parameter :: name='psb_get_xch_idx' + + totxch = 0 + totsnd = 0 + totrcv = 0 + ip = 1 + + do + if (ip > size(idx)) then + write(0,*) trim(name),': Warning: out of size of input vector ' + exit + end if + if (idx(ip) == -1) exit + totxch = totxch+1 + nerv = idx(ip+psb_n_elem_recv_) + nesd = idx(ip+nerv+psb_n_elem_send_) + totsnd = totsnd + nesd + totrcv = totrcv + nerv + ip = ip+nerv+nesd+3 + end do + + end subroutine psb_get_xch_idx + + subroutine psb_cd_get_list(data,desc,ipnt,totxch,idxr,idxs,info) use psb_const_mod @@ -600,33 +704,30 @@ contains select case(data) case(psb_comm_halo_) ipnt => desc%halo_index - totxch = desc%matrix_data(psb_thal_xch_) - idxr = desc%matrix_data(psb_thal_rcv_) - idxs = desc%matrix_data(psb_thal_snd_) - case(psb_comm_ovr_) ipnt => desc%ovrlap_index - totxch = desc%matrix_data(psb_tovr_xch_) - idxr = desc%matrix_data(psb_tovr_rcv_) - idxs = desc%matrix_data(psb_tovr_snd_) - case(psb_comm_ext_) ipnt => desc%ext_index - totxch = desc%matrix_data(psb_text_xch_) - idxr = desc%matrix_data(psb_text_rcv_) - idxs = desc%matrix_data(psb_text_snd_) - + if (.not.associated(desc%base_desc)) then + write(0,*) trim(name),& + & ': Warning: trying to get ext_index on a descriptor ',& + & 'which does not have a base_desc!' + end if + if (.not.psb_is_ovl_desc(desc)) then + write(0,*) trim(name),& + & ': Warning: trying to get ext_index on a descriptor ',& + & 'which is not overlap-extended!' + end if case(psb_comm_mov_) ipnt => desc%ovr_mst_idx - totxch = desc%matrix_data(psb_tmov_xch_) - idxr = desc%matrix_data(psb_tmov_rcv_) - idxs = desc%matrix_data(psb_tmov_snd_) - case default info=4010 call psb_errpush(info,name,a_err='wrong Data selector') goto 9999 end select + call psb_get_xch_idx(ipnt,totxch,idxs,idxr) + + call psb_erractionrestore(err_act) return diff --git a/base/modules/psb_inter_desc_type.f90 b/base/modules/psb_inter_desc_type.f90 new file mode 100644 index 00000000..a01598a1 --- /dev/null +++ b/base/modules/psb_inter_desc_type.f90 @@ -0,0 +1,655 @@ +module psb_inter_descriptor_type + use psb_spmat_type + use psb_descriptor_type + + + + ! Inter-descriptor mapping data structures. + integer, parameter :: psb_map_kind_ = 1 + integer, parameter :: psb_map_data_ = 2 + integer, parameter :: psb_map_integer_ = 1 + integer, parameter :: psb_map_double_ = 2 + integer, parameter :: psb_map_complex_ = 3 + + integer, parameter :: psb_fw_tmp_kind_ = 5 + integer, parameter :: psb_fw_tmp_sz_ = 6 + integer, parameter :: psb_bk_tmp_kind_ = 7 + integer, parameter :: psb_bk_tmp_sz_ = 8 + integer, parameter :: psb_itd_data_size_=20 + + + type psb_d_map_type + type(psb_dspmat_type) :: map_fw, map_bk + end type psb_d_map_type + type psb_z_map_type + type(psb_zspmat_type) :: map_fw, map_bk + end type psb_z_map_type + + type psb_inter_desc_type + integer, allocatable :: itd_data(:) + type(psb_desc_type), pointer :: desc_1=>null(), desc_2=>null() + integer, allocatable :: exch_fw_idx(:), exch_bk_idx(:) + type(psb_d_map_type) :: dmap + type(psb_z_map_type) :: zmap + end type psb_inter_desc_type + + interface psb_forward_map + module procedure psb_d_forward_map, psb_z_forward_map + end interface + + interface psb_backward_map + module procedure psb_d_backward_map, psb_z_backward_map + end interface + + interface psb_is_ok_desc + module procedure psb_is_ok_inter_desc + end interface + + interface psb_is_asb_desc + module procedure psb_is_asb_inter_desc + end interface + + interface psb_inter_desc + module procedure psb_d_inter_desc, psb_d_inter_desc_noidx,& + & psb_z_inter_desc, psb_z_inter_desc_noidx + end interface + + interface psb_sizeof + module procedure psb_itd_sizeof,& + & psb_d_map_sizeof, psb_z_map_sizeof + end interface + + +contains + + function psb_cd_get_map_kind(desc) + implicit none + type(psb_inter_desc_type), intent(in) :: desc + Integer :: psb_cd_get_map_kind + if (psb_is_ok_desc(desc)) then + psb_cd_get_map_kind = desc%itd_data(psb_map_kind_) + else + psb_cd_get_map_kind = -1 + end if + end function psb_cd_get_map_kind + + subroutine psb_cd_set_map_kind(map_kind,desc) + implicit none + integer, intent(in) :: map_kind + type(psb_inter_desc_type), intent(inout) :: desc + + desc%itd_data(psb_map_kind_) = map_kind + + end subroutine psb_cd_set_map_kind + + function psb_cd_get_map_data(desc) + implicit none + type(psb_inter_desc_type), intent(in) :: desc + Integer :: psb_cd_get_map_data + if (psb_is_ok_desc(desc)) then + psb_cd_get_map_data = desc%itd_data(psb_map_data_) + else + psb_cd_get_map_data = -1 + end if + end function psb_cd_get_map_data + + subroutine psb_cd_set_map_data(map_data,desc) + implicit none + integer, intent(in) :: map_data + type(psb_inter_desc_type), intent(inout) :: desc + + + desc%itd_data(psb_map_data_) = map_data + + end subroutine psb_cd_set_map_data + + + function psb_cd_get_fw_tmp_sz(desc) + implicit none + type(psb_inter_desc_type), intent(in) :: desc + Integer :: psb_cd_get_fw_tmp_sz + + psb_cd_get_fw_tmp_sz = desc%itd_data(psb_fw_tmp_sz_) + end function psb_cd_get_fw_tmp_sz + + function psb_cd_get_bk_tmp_sz(desc) + implicit none + type(psb_inter_desc_type), intent(in) :: desc + Integer :: psb_cd_get_bk_tmp_sz + + psb_cd_get_bk_tmp_sz = desc%itd_data(psb_bk_tmp_sz_) + end function psb_cd_get_bk_tmp_sz + + subroutine psb_cd_set_fw_tmp_sz(isz,desc) + implicit none + type(psb_inter_desc_type), intent(inout) :: desc + integer, intent(in) :: isz + + desc%itd_data(psb_fw_tmp_sz_) =isz + end subroutine psb_cd_set_fw_tmp_sz + + subroutine psb_cd_set_bk_tmp_sz(isz,desc) + implicit none + type(psb_inter_desc_type), intent(inout) :: desc + integer, intent(in) :: isz + + desc%itd_data(psb_bk_tmp_sz_) =isz + + end subroutine psb_cd_set_bk_tmp_sz + + + logical function psb_is_asb_inter_desc(desc) + type(psb_inter_desc_type), intent(in) :: desc + + psb_is_asb_inter_desc = .false. + if (.not.allocated(desc%itd_data)) return + if (.not.associated(desc%desc_1)) return + if (.not.associated(desc%desc_2)) return + psb_is_asb_inter_desc = & + & psb_is_asb_desc(desc%desc_1).and.psb_is_asb_desc(desc%desc_2) + + end function psb_is_asb_inter_desc + + logical function psb_is_ok_inter_desc(desc) + type(psb_inter_desc_type), intent(in) :: desc + + psb_is_ok_inter_desc = .false. + if (.not.allocated(desc%itd_data)) return + select case(desc%itd_data(psb_map_data_)) + case(psb_map_integer_, psb_map_double_, psb_map_complex_) + ! Ok go ahead + case default + ! Since it's false so far, simply return + return + end select + if (.not.associated(desc%desc_1)) return + if (.not.associated(desc%desc_2)) return + psb_is_ok_inter_desc = & + & psb_is_ok_desc(desc%desc_1).and.psb_is_ok_desc(desc%desc_2) + + end function psb_is_ok_inter_desc + + + function psb_d_map_sizeof(map) + implicit none + type(psb_d_map_type), intent(in) :: map + Integer :: psb_d_map_sizeof + integer :: val + + val = 0 + + val = val + psb_sizeof(map%map_fw) + val = val + psb_sizeof(map%map_bk) + psb_d_map_sizeof = val + end function psb_d_map_sizeof + + function psb_z_map_sizeof(map) + implicit none + type(psb_z_map_type), intent(in) :: map + Integer :: psb_z_map_sizeof + integer :: val + + val = 0 + + val = val + psb_sizeof(map%map_fw) + val = val + psb_sizeof(map%map_bk) + psb_z_map_sizeof = val + end function psb_z_map_sizeof + + function psb_itd_sizeof(desc) + + implicit none + type(psb_inter_desc_type), intent(in) :: desc + Integer :: psb_itd_sizeof + integer :: val + + val = 0 + + if (allocated(desc%itd_data)) val = val + 4*size(desc%itd_data) + if (allocated(desc%exch_fw_idx)) val = val + 4*size(desc%exch_fw_idx) + if (allocated(desc%exch_bk_idx)) val = val + 4*size(desc%exch_bk_idx) + val = val + psb_sizeof(desc%dmap) + val = val + psb_sizeof(desc%zmap) + psb_itd_sizeof = val + end function psb_itd_sizeof + + function psb_d_inter_desc(map_kind,desc1, desc2, map_fw, map_bk, idx_fw, idx_bk) + use psb_serial_mod + use psi_mod + implicit none + type(psb_inter_desc_type) :: psb_d_inter_desc + type(psb_desc_type), target :: desc1, desc2 + type(psb_dspmat_type), intent(in) :: map_fw, map_bk + integer, intent(in) :: map_kind,idx_fw(:), idx_bk(:) + ! + type(psb_inter_desc_type) :: this + integer :: info + character(len=20), parameter :: name='psb_inter_desc' + + info = 0 + if (psb_is_ok_desc(desc1)) then + this%desc_1=>desc1 + else + info = 2 + endif + if (psb_is_ok_desc(desc2)) then + this%desc_2=>desc2 + else + info = 3 + endif + + if (info == 0) call psb_sp_clone(map_fw,this%dmap%map_fw,info) + if (info == 0) call psb_sp_clone(map_bk,this%dmap%map_bk,info) + if (info == 0) call psb_safe_cpy(idx_fw,this%exch_fw_idx,info) + if (info == 0) call psb_safe_cpy(idx_bk,this%exch_bk_idx,info) + if (info == 0) call psb_realloc(psb_itd_data_size_,this%itd_data,info) + if (info == 0) then + call psb_cd_set_map_kind(map_kind, this) + call psb_cd_set_map_data(psb_map_double_, this) +!!$ call psb_cd_set_fw_tmp_sz(map_fw%k, this) +!!$ call psb_cd_set_bk_tmp_sz(map_bk%k, this) + end if + if (info /= 0) then + write(0,*) trim(name),' Invalid descriptor input' + return + end if + + psb_d_inter_desc = this + + end function psb_d_inter_desc + + function psb_d_inter_desc_noidx(map_kind,desc1, desc2, map_fw, map_bk) + use psb_serial_mod + use psi_mod + implicit none + type(psb_inter_desc_type) :: psb_d_inter_desc_noidx + type(psb_desc_type), target :: desc1, desc2 + type(psb_dspmat_type), intent(in) :: map_fw, map_bk + integer, intent(in) :: map_kind + ! + type(psb_inter_desc_type) :: this + integer :: info + character(len=20), parameter :: name='psb_inter_desc' + + info = 0 + select case(map_kind) + case (psb_map_aggr_) + ! OK + case default + write(0,*) 'Bad mapp kind into psb_inter_desc ',map_kind + info = 1 + end select + + if (psb_is_ok_desc(desc1)) then + this%desc_1=>desc1 + else + info = 2 + endif + if (psb_is_ok_desc(desc2)) then + this%desc_2=>desc2 + else + info = 3 + endif + + if (info == 0) call psb_sp_clone(map_fw,this%dmap%map_fw,info) + if (info == 0) call psb_sp_clone(map_bk,this%dmap%map_bk,info) + if (info == 0) call psb_realloc(psb_itd_data_size_,this%itd_data,info) + if (info == 0) then + call psb_cd_set_map_kind(map_kind, this) + call psb_cd_set_map_data(psb_map_double_, this) +!!$ call psb_cd_set_fw_tmp_sz(map_fw%k, this) +!!$ call psb_cd_set_bk_tmp_sz(map_bk%k, this) + end if + if (info /= 0) then + write(0,*) trim(name),' Invalid descriptor input' + return + end if + + psb_d_inter_desc_noidx = this + + end function psb_d_inter_desc_noidx + + function psb_z_inter_desc(map_kind,desc1, desc2, map_fw, map_bk, idx_fw, idx_bk) + use psb_serial_mod + use psi_mod + implicit none + type(psb_inter_desc_type) :: psb_z_inter_desc + type(psb_desc_type), target :: desc1, desc2 + type(psb_zspmat_type), intent(in) :: map_fw, map_bk + integer, intent(in) :: map_kind,idx_fw(:), idx_bk(:) + ! + type(psb_inter_desc_type) :: this + integer :: info + character(len=20), parameter :: name='psb_inter_desc' + + info = 0 + if (psb_is_ok_desc(desc1)) then + this%desc_1=>desc1 + else + info = 2 + endif + if (psb_is_ok_desc(desc2)) then + this%desc_2=>desc2 + else + info = 3 + endif + + if (info == 0) call psb_sp_clone(map_fw,this%zmap%map_fw,info) + if (info == 0) call psb_sp_clone(map_bk,this%zmap%map_bk,info) + if (info == 0) call psb_safe_cpy(idx_fw,this%exch_fw_idx,info) + if (info == 0) call psb_safe_cpy(idx_bk,this%exch_bk_idx,info) + if (info == 0) call psb_realloc(psb_itd_data_size_,this%itd_data,info) + if (info == 0) then + call psb_cd_set_map_kind(map_kind, this) + call psb_cd_set_map_data(psb_map_complex_, this) +!!$ call psb_cd_set_fw_tmp_sz(map_fw%k, this) +!!$ call psb_cd_set_bk_tmp_sz(map_bk%k, this) + end if + if (info /= 0) then + write(0,*) trim(name),' Invalid descriptor input' + return + end if + + psb_z_inter_desc = this + + end function psb_z_inter_desc + + function psb_z_inter_desc_noidx(map_kind,desc1, desc2, map_fw, map_bk) + use psb_serial_mod + use psi_mod + implicit none + type(psb_inter_desc_type) :: psb_z_inter_desc_noidx + type(psb_desc_type), target :: desc1, desc2 + type(psb_zspmat_type), intent(in) :: map_fw, map_bk + integer, intent(in) :: map_kind + ! + type(psb_inter_desc_type) :: this + integer :: info + character(len=20), parameter :: name='psb_inter_desc' + + info = 0 + select case(map_kind) + case (psb_map_aggr_) + ! OK + case default + write(0,*) 'Bad mapp kind into psb_inter_desc ',map_kind + info = 1 + end select + + if (psb_is_ok_desc(desc1)) then + this%desc_1=>desc1 + else + info = 2 + endif + if (psb_is_ok_desc(desc2)) then + this%desc_2=>desc2 + else + info = 3 + endif + + if (info == 0) call psb_sp_clone(map_fw,this%zmap%map_fw,info) + if (info == 0) call psb_sp_clone(map_bk,this%zmap%map_bk,info) + if (info == 0) call psb_realloc(psb_itd_data_size_,this%itd_data,info) + if (info == 0) then + call psb_cd_set_map_kind(map_kind, this) + call psb_cd_set_map_data(psb_map_complex_, this) +!!$ call psb_cd_set_fw_tmp_sz(map_fw%k, this) +!!$ call psb_cd_set_bk_tmp_sz(map_bk%k, this) + end if + if (info /= 0) then + write(0,*) trim(name),' Invalid descriptor input' + return + end if + + psb_z_inter_desc_noidx = this + + end function psb_z_inter_desc_noidx + + + + + ! + ! Takes a vector X from space desc%desc_1 and maps it onto + ! desc%desc_2 under desc%map_fw possibly with communication + ! due to exch_fw_idx + ! + subroutine psb_d_forward_map(alpha,x,beta,y,desc,info,work) + use psb_comm_mod + use psb_serial_mod + use psi_mod + implicit none + type(psb_inter_desc_type), intent(in) :: desc + real(kind(1.d0)), intent(in) :: alpha,beta + real(kind(1.d0)), intent(inout) :: x(:) + real(kind(1.d0)), intent(out) :: y(:) + integer, intent(out) :: info + real(kind(1.d0)), optional :: work(:) + + ! + real(kind(1.d0)), allocatable :: xt(:) + integer :: itsz, i, j,totxch,totsnd,totrcv,& + & map_kind, map_data + character(len=20), parameter :: name='psb_forward_map' + + info = 0 + if (.not.psb_is_asb_desc(desc)) then + write(0,*) trim(name),' Invalid descriptor inupt' + info = 1 + return + end if + + itsz = psb_cd_get_fw_tmp_sz(desc) + map_kind = psb_cd_get_map_kind(desc) + map_data = psb_cd_get_map_data(desc) + if (map_data /= psb_map_double_) then + write(0,*) trim(name),' Invalid descriptor inupt: map_data', & + & map_data,psb_map_double_ + info = 1 + return + endif + + select case(map_kind) + case(psb_map_aggr_) + ! Ok, we just need to call a halo update and a matrix-vector product. + call psb_halo(x,desc%desc_1,info,work=work) + if (info == 0) call psb_csmm(alpha,desc%dmap%map_fw,x,beta,y,info) + + if (info /= 0) then + write(0,*) trim(name),' Error from inner routines',info + info = -1 + end if + + + case default + write(0,*) trim(name),' Invalid descriptor inupt' + info = 1 + return + end select + + end subroutine psb_d_forward_map + + + ! + ! Takes a vector X from space desc%desc_2 and maps it onto + ! desc%desc_1 under desc%map_bk possibly with communication + ! due to exch_bk_idx + ! + subroutine psb_d_backward_map(alpha,x,beta,y,desc,info,work) + use psb_comm_mod + use psb_serial_mod + use psi_mod + implicit none + type(psb_inter_desc_type), intent(in) :: desc + real(kind(1.d0)), intent(in) :: alpha,beta + real(kind(1.d0)), intent(inout) :: x(:) + real(kind(1.d0)), intent(out) :: y(:) + integer, intent(out) :: info + real(kind(1.d0)), optional :: work(:) + + ! + real(kind(1.d0)), allocatable :: xt(:) + integer :: itsz, i, j,totxch,totsnd,totrcv,& + & map_kind, map_data + character(len=20), parameter :: name='psb_backward_map' + + info = 0 + if (.not.psb_is_asb_desc(desc)) then + write(0,*) trim(name),' Invalid descriptor inupt' + info = 1 + return + end if + + itsz = psb_cd_get_bk_tmp_sz(desc) + map_kind = psb_cd_get_map_kind(desc) + map_data = psb_cd_get_map_data(desc) + if (map_data /= psb_map_double_) then + write(0,*) trim(name),' Invalid descriptor inupt: map_data',& + & map_data,psb_map_double_ + info = 1 + return + endif + + select case(map_kind) + case(psb_map_aggr_) + ! Ok, we just need to call a halo update and a matrix-vector product. + call psb_halo(x,desc%desc_2,info,work=work) + if (info == 0) call psb_csmm(alpha,desc%dmap%map_bk,x,beta,y,info) + + if (info /= 0) then + write(0,*) trim(name),' Error from inner routines',info + info = -1 + end if + + + case default + write(0,*) trim(name),' Invalid descriptor inupt' + info = 1 + return + end select + + end subroutine psb_d_backward_map + + + ! + ! Takes a vector X from space desc%desc_1 and maps it onto + ! desc%desc_2 under desc%map_fw possibly with communication + ! due to exch_fw_idx + ! + subroutine psb_z_forward_map(alpha,x,beta,y,desc,info,work) + use psb_comm_mod + use psb_serial_mod + use psi_mod + implicit none + type(psb_inter_desc_type), intent(in) :: desc + complex(kind(1.d0)), intent(in) :: alpha,beta + complex(kind(1.d0)), intent(inout) :: x(:) + complex(kind(1.d0)), intent(out) :: y(:) + integer, intent(out) :: info + complex(kind(1.d0)), optional :: work(:) + + ! + complex(kind(1.d0)), allocatable :: xt(:) + integer :: itsz, i, j,totxch,totsnd,totrcv,& + & map_kind, map_data + character(len=20), parameter :: name='psb_forward_map' + + info = 0 + if (.not.psb_is_asb_desc(desc)) then + write(0,*) trim(name),' Invalid descriptor inupt' + info = 1 + return + end if + + itsz = psb_cd_get_fw_tmp_sz(desc) + map_kind = psb_cd_get_map_kind(desc) + map_data = psb_cd_get_map_data(desc) + if (map_data /= psb_map_complex_) then + write(0,*) trim(name),' Invalid descriptor inupt: map_data',& + & map_data,psb_map_complex_ + info = 1 + return + endif + + select case(map_kind) + case(psb_map_aggr_) + ! Ok, we just need to call a halo update and a matrix-vector product. + call psb_halo(x,desc%desc_1,info,work=work) + if (info == 0) call psb_csmm(alpha,desc%zmap%map_fw,x,beta,y,info) + + if (info /= 0) then + write(0,*) trim(name),' Error from inner routines',info + info = -1 + end if + + + case default + write(0,*) trim(name),' Invalid descriptor inupt' + info = 1 + return + end select + + end subroutine psb_z_forward_map + + + ! + ! Takes a vector X from space desc%desc_2 and maps it onto + ! desc%desc_1 under desc%map_bk possibly with communication + ! due to exch_bk_idx + ! + subroutine psb_z_backward_map(alpha,x,beta,y,desc,info,work) + use psb_comm_mod + use psb_serial_mod + use psi_mod + implicit none + type(psb_inter_desc_type), intent(in) :: desc + complex(kind(1.d0)), intent(in) :: alpha,beta + complex(kind(1.d0)), intent(inout) :: x(:) + complex(kind(1.d0)), intent(out) :: y(:) + integer, intent(out) :: info + complex(kind(1.d0)), optional :: work(:) + + ! + complex(kind(1.d0)), allocatable :: xt(:) + integer :: itsz, i, j,totxch,totsnd,totrcv,& + & map_kind, map_data + character(len=20), parameter :: name='psb_backward_map' + + info = 0 + if (.not.psb_is_asb_desc(desc)) then + write(0,*) trim(name),' Invalid descriptor inupt' + info = 1 + return + end if + + itsz = psb_cd_get_bk_tmp_sz(desc) + map_kind = psb_cd_get_map_kind(desc) + map_data = psb_cd_get_map_data(desc) + if (map_data /= psb_map_complex_) then + write(0,*) trim(name),' Invalid descriptor inupt: map_data',& + & map_data,psb_map_complex_ + info = 1 + return + endif + + select case(map_kind) + case(psb_map_aggr_) + ! Ok, we just need to call a halo update and a matrix-vector product. + call psb_halo(x,desc%desc_2,info,work=work) + if (info == 0) call psb_csmm(alpha,desc%zmap%map_bk,x,beta,y,info) + + if (info /= 0) then + write(0,*) trim(name),' Error from inner routines',info + info = -1 + end if + + + case default + write(0,*) trim(name),' Invalid descriptor inupt' + info = 1 + return + end select + + end subroutine psb_z_backward_map + + +end module psb_inter_descriptor_type diff --git a/base/modules/psb_realloc_mod.F90 b/base/modules/psb_realloc_mod.F90 index 3c19a6f2..73449d79 100644 --- a/base/modules/psb_realloc_mod.F90 +++ b/base/modules/psb_realloc_mod.F90 @@ -57,9 +57,14 @@ module psb_realloc_mod module procedure psb_ztransfer2d end interface + Interface psb_safe_ab_cpy + module procedure psb_i_ab_cpy1d,psb_i_ab_cpy2d, & + & psb_d_ab_cpy1d, psb_d_ab_cpy2d, psb_z_ab_cpy1d, psb_z_ab_cpy2d + end Interface + Interface psb_safe_cpy - module procedure psb_icpy1d,psb_icpy2d, & - & psb_dcpy1d, psb_dcpy2d, psb_zcpy1d, psb_zcpy2d + module procedure psb_i_cpy1d,psb_i_cpy2d, & + & psb_d_cpy1d, psb_d_cpy2d, psb_z_cpy1d, psb_z_cpy2d end Interface ! @@ -80,7 +85,7 @@ module psb_realloc_mod Contains - subroutine psb_icpy1d(vin,vout,info) + subroutine psb_i_ab_cpy1d(vin,vout,info) use psb_error_mod ! ...Subroutine Arguments @@ -93,7 +98,7 @@ Contains character(len=20) :: name, char_err logical, parameter :: debug=.false. - name='psb_cpy1d' + name='psb_safe_ab_cpy' call psb_erractionsave(err_act) if(psb_get_errstatus() /= 0) return @@ -125,9 +130,9 @@ Contains end if return - end subroutine psb_icpy1d + end subroutine psb_i_ab_cpy1d - subroutine psb_icpy2d(vin,vout,info) + subroutine psb_i_ab_cpy2d(vin,vout,info) use psb_error_mod ! ...Subroutine Arguments @@ -140,7 +145,7 @@ Contains character(len=20) :: name, char_err logical, parameter :: debug=.false. - name='psb_cpy1d' + name='psb_safe_ab_cpy' call psb_erractionsave(err_act) if(psb_get_errstatus() /= 0) return @@ -174,9 +179,9 @@ Contains end if return - end subroutine psb_icpy2d + end subroutine psb_i_ab_cpy2d - subroutine psb_dcpy1d(vin,vout,info) + subroutine psb_d_ab_cpy1d(vin,vout,info) use psb_error_mod ! ...Subroutine Arguments @@ -189,7 +194,7 @@ Contains character(len=20) :: name, char_err logical, parameter :: debug=.false. - name='psb_cpy1d' + name='psb_safe_ab_cpy' call psb_erractionsave(err_act) if(psb_get_errstatus() /= 0) return @@ -221,9 +226,9 @@ Contains end if return - end subroutine psb_dcpy1d + end subroutine psb_d_ab_cpy1d - subroutine psb_dcpy2d(vin,vout,info) + subroutine psb_d_ab_cpy2d(vin,vout,info) use psb_error_mod ! ...Subroutine Arguments @@ -236,7 +241,7 @@ Contains character(len=20) :: name, char_err logical, parameter :: debug=.false. - name='psb_cpy1d' + name='psb_safe_ab_cpy' call psb_erractionsave(err_act) if(psb_get_errstatus() /= 0) return @@ -270,9 +275,9 @@ Contains end if return - end subroutine psb_dcpy2d + end subroutine psb_d_ab_cpy2d - subroutine psb_zcpy1d(vin,vout,info) + subroutine psb_z_ab_cpy1d(vin,vout,info) use psb_error_mod ! ...Subroutine Arguments @@ -285,7 +290,7 @@ Contains character(len=20) :: name, char_err logical, parameter :: debug=.false. - name='psb_cpy1d' + name='psb_safe_ab_cpy' call psb_erractionsave(err_act) if(psb_get_errstatus() /= 0) return @@ -317,9 +322,9 @@ Contains end if return - end subroutine psb_zcpy1d + end subroutine psb_z_ab_cpy1d - subroutine psb_zcpy2d(vin,vout,info) + subroutine psb_z_ab_cpy2d(vin,vout,info) use psb_error_mod ! ...Subroutine Arguments @@ -332,7 +337,7 @@ Contains character(len=20) :: name, char_err logical, parameter :: debug=.false. - name='psb_cpy1d' + name='psb_safe_ab_cpy' call psb_erractionsave(err_act) if(psb_get_errstatus() /= 0) return @@ -366,7 +371,285 @@ Contains end if return - end subroutine psb_zcpy2d + end subroutine psb_z_ab_cpy2d + + + subroutine psb_i_cpy1d(vin,vout,info) + use psb_error_mod + + ! ...Subroutine Arguments + Integer, intent(in) :: vin(:) + Integer, allocatable, intent(out) :: vout(:) + integer :: info + ! ...Local Variables + + Integer :: isz,err_act,lb + character(len=20) :: name, char_err + logical, parameter :: debug=.false. + + name='psb_safe_cpy' + call psb_erractionsave(err_act) + + if(psb_get_errstatus() /= 0) return + info = 0 + isz = size(vin) + lb = lbound(vin,1) + call psb_realloc(isz,vout,info,lb=lb) + if (info /= 0) then + info=4010 + char_err='psb_realloc' + call psb_errpush(info,name,a_err=char_err) + goto 9999 + else + vout(:) = vin(:) + endif + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + + if (err_act == psb_act_ret_) then + return + else + call psb_error() + end if + return + + end subroutine psb_i_cpy1d + + subroutine psb_i_cpy2d(vin,vout,info) + use psb_error_mod + + ! ...Subroutine Arguments + Integer, intent(in) :: vin(:,:) + Integer, allocatable, intent(out) :: vout(:,:) + integer :: info + ! ...Local Variables + + Integer :: isz1, isz2,err_act, lb1, lb2 + character(len=20) :: name, char_err + logical, parameter :: debug=.false. + + name='psb_safe_cpy' + call psb_erractionsave(err_act) + + if(psb_get_errstatus() /= 0) return + info = 0 + isz1 = size(vin,1) + isz2 = size(vin,2) + lb1 = lbound(vin,1) + lb2 = lbound(vin,2) + call psb_realloc(isz1,isz2,vout,info,lb1=lb1,lb2=lb2) + if (info /= 0) then + info=4010 + char_err='psb_realloc' + call psb_errpush(info,name,a_err=char_err) + goto 9999 + else + vout(:,:) = vin(:,:) + endif + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + + if (err_act == psb_act_ret_) then + return + else + call psb_error() + end if + return + + end subroutine psb_i_cpy2d + + subroutine psb_d_cpy1d(vin,vout,info) + use psb_error_mod + + ! ...Subroutine Arguments + real(kind(1.d0)), intent(in) :: vin(:) + real(kind(1.d0)), allocatable, intent(out) :: vout(:) + integer :: info + ! ...Local Variables + + Integer :: isz,err_act,lb + character(len=20) :: name, char_err + logical, parameter :: debug=.false. + + name='psb_safe_cpy' + call psb_erractionsave(err_act) + + if(psb_get_errstatus() /= 0) return + info = 0 + isz = size(vin) + lb = lbound(vin,1) + call psb_realloc(isz,vout,info,lb=lb) + if (info /= 0) then + info=4010 + char_err='psb_realloc' + call psb_errpush(info,name,a_err=char_err) + goto 9999 + else + vout(:) = vin(:) + endif + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + + if (err_act == psb_act_ret_) then + return + else + call psb_error() + end if + return + + end subroutine psb_d_cpy1d + + subroutine psb_d_cpy2d(vin,vout,info) + use psb_error_mod + + ! ...Subroutine Arguments + real(kind(1.d0)), intent(in) :: vin(:,:) + real(kind(1.d0)), allocatable, intent(out) :: vout(:,:) + integer :: info + ! ...Local Variables + + Integer :: isz1, isz2,err_act, lb1, lb2 + character(len=20) :: name, char_err + logical, parameter :: debug=.false. + + name='psb_safe_cpy' + call psb_erractionsave(err_act) + + if(psb_get_errstatus() /= 0) return + info = 0 + isz1 = size(vin,1) + isz2 = size(vin,2) + lb1 = lbound(vin,1) + lb2 = lbound(vin,2) + call psb_realloc(isz1,isz2,vout,info,lb1=lb1,lb2=lb2) + if (info /= 0) then + info=4010 + char_err='psb_realloc' + call psb_errpush(info,name,a_err=char_err) + goto 9999 + else + vout(:,:) = vin(:,:) + endif + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + + if (err_act == psb_act_ret_) then + return + else + call psb_error() + end if + return + + end subroutine psb_d_cpy2d + + subroutine psb_z_cpy1d(vin,vout,info) + use psb_error_mod + + ! ...Subroutine Arguments + complex(kind(1.d0)), intent(in) :: vin(:) + complex(kind(1.d0)), allocatable, intent(out) :: vout(:) + integer :: info + ! ...Local Variables + + Integer :: isz,err_act,lb + character(len=20) :: name, char_err + logical, parameter :: debug=.false. + + name='psb_safe_cpy' + call psb_erractionsave(err_act) + + if(psb_get_errstatus() /= 0) return + info = 0 + isz = size(vin) + lb = lbound(vin,1) + call psb_realloc(isz,vout,info,lb=lb) + if (info /= 0) then + info=4010 + char_err='psb_realloc' + call psb_errpush(info,name,a_err=char_err) + goto 9999 + else + vout(:) = vin(:) + endif + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + + if (err_act == psb_act_ret_) then + return + else + call psb_error() + end if + return + + end subroutine psb_z_cpy1d + + subroutine psb_z_cpy2d(vin,vout,info) + use psb_error_mod + + ! ...Subroutine Arguments + complex(kind(1.d0)), intent(in) :: vin(:,:) + complex(kind(1.d0)), allocatable, intent(out) :: vout(:,:) + integer :: info + ! ...Local Variables + + Integer :: isz1, isz2,err_act, lb1, lb2 + character(len=20) :: name, char_err + logical, parameter :: debug=.false. + + name='psb_safe_cpy' + call psb_erractionsave(err_act) + + if(psb_get_errstatus() /= 0) return + info = 0 + isz1 = size(vin,1) + isz2 = size(vin,2) + lb1 = lbound(vin,1) + lb2 = lbound(vin,2) + call psb_realloc(isz1,isz2,vout,info,lb1=lb1,lb2=lb2) + if (info /= 0) then + info=4010 + char_err='psb_realloc' + call psb_errpush(info,name,a_err=char_err) + goto 9999 + else + vout(:,:) = vin(:,:) + endif + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + + if (err_act == psb_act_ret_) then + return + else + call psb_error() + end if + return + + end subroutine psb_z_cpy2d + function psb_isize1d(vin) integer :: psb_isize1d diff --git a/base/modules/psb_serial_mod.f90 b/base/modules/psb_serial_mod.f90 index 2a2e10c7..637f249d 100644 --- a/base/modules/psb_serial_mod.f90 +++ b/base/modules/psb_serial_mod.f90 @@ -33,6 +33,11 @@ module psb_serial_mod use psb_string_mod use psb_sort_mod + use psi_serial_mod, & + & psb_gth => psi_gth,& + & psb_sct => psi_sct + + interface psb_csdp subroutine psb_dcsdp(a, b,info,ifc,check,trans,unitd,upd,dupl) use psb_spmat_type diff --git a/base/modules/psb_spmat_type.f90 b/base/modules/psb_spmat_type.f90 index 44b07f21..8a95a67b 100644 --- a/base/modules/psb_spmat_type.f90 +++ b/base/modules/psb_spmat_type.f90 @@ -610,11 +610,11 @@ contains INFO = 0 call psb_nullify_sp(b) - call psb_safe_cpy(a%aspk,b%aspk,info) - if (info == 0) call psb_safe_cpy(a%ia1,b%ia1,info) - if (info == 0) call psb_safe_cpy(a%ia2,b%ia2,info) - if (info == 0) call psb_safe_cpy(a%pl,b%pl,info) - if (info == 0) call psb_safe_cpy(a%pr,b%pr,info) + call psb_safe_ab_cpy(a%aspk,b%aspk,info) + if (info == 0) call psb_safe_ab_cpy(a%ia1,b%ia1,info) + if (info == 0) call psb_safe_ab_cpy(a%ia2,b%ia2,info) + if (info == 0) call psb_safe_ab_cpy(a%pl,b%pl,info) + if (info == 0) call psb_safe_ab_cpy(a%pr,b%pr,info) if (info /= 0) then info=2023 return @@ -1141,11 +1141,11 @@ contains INFO = 0 call psb_nullify_sp(b) - call psb_safe_cpy(a%aspk,b%aspk,info) - if (info == 0) call psb_safe_cpy(a%ia1,b%ia1,info) - if (info == 0) call psb_safe_cpy(a%ia2,b%ia2,info) - if (info == 0) call psb_safe_cpy(a%pl,b%pl,info) - if (info == 0) call psb_safe_cpy(a%pr,b%pr,info) + call psb_safe_ab_cpy(a%aspk,b%aspk,info) + if (info == 0) call psb_safe_ab_cpy(a%ia1,b%ia1,info) + if (info == 0) call psb_safe_ab_cpy(a%ia2,b%ia2,info) + if (info == 0) call psb_safe_ab_cpy(a%pl,b%pl,info) + if (info == 0) call psb_safe_ab_cpy(a%pr,b%pr,info) if (info /= 0) then info=2023 return diff --git a/base/modules/psb_tools_mod.f90 b/base/modules/psb_tools_mod.f90 index 108f7ed1..6f20a0a6 100644 --- a/base/modules/psb_tools_mod.f90 +++ b/base/modules/psb_tools_mod.f90 @@ -30,7 +30,6 @@ !!$ Module psb_tools_mod use psb_const_mod - use psb_gps_mod use psb_spmat_type interface psb_geall @@ -337,22 +336,22 @@ Module psb_tools_mod Subroutine psb_dcdovr(a,desc_a,novr,desc_ov,info,extype) use psb_descriptor_type Use psb_spmat_type - integer, intent(in) :: novr - Type(psb_dspmat_type), Intent(in) :: a - Type(psb_desc_type), Intent(in) :: desc_a - Type(psb_desc_type), Intent(out) :: desc_ov - integer, intent(out) :: info - integer, intent(in),optional :: extype + integer, intent(in) :: novr + Type(psb_dspmat_type), Intent(in) :: a + Type(psb_desc_type), Intent(in), target :: desc_a + Type(psb_desc_type), Intent(out) :: desc_ov + integer, intent(out) :: info + integer, intent(in),optional :: extype end Subroutine psb_dcdovr Subroutine psb_zcdovr(a,desc_a,novr,desc_ov,info,extype) use psb_descriptor_type Use psb_spmat_type - integer, intent(in) :: novr - Type(psb_zspmat_type), Intent(in) :: a - Type(psb_desc_type), Intent(in) :: desc_a - Type(psb_desc_type), Intent(out) :: desc_ov - integer, intent(out) :: info - integer, intent(in),optional :: extype + integer, intent(in) :: novr + Type(psb_zspmat_type), Intent(in) :: a + Type(psb_desc_type), Intent(in), target :: desc_a + Type(psb_desc_type), Intent(out) :: desc_ov + integer, intent(out) :: info + integer, intent(in),optional :: extype end Subroutine psb_zcdovr end interface @@ -490,6 +489,27 @@ Module psb_tools_mod logical, intent(in), optional :: owned character, intent(in), optional :: iact end subroutine psb_glob_to_loc + subroutine psb_glob_to_loc2s(x,y,desc_a,info,iact,owned) + use psb_descriptor_type + implicit none + type(psb_desc_type), intent(in) :: desc_a + integer,intent(in) :: x + integer,intent(out) :: y + integer, intent(out) :: info + character, intent(in), optional :: iact + logical, intent(in), optional :: owned + + end subroutine psb_glob_to_loc2s + + subroutine psb_glob_to_locs(x,desc_a,info,iact,owned) + use psb_descriptor_type + implicit none + type(psb_desc_type), intent(in) :: desc_a + integer,intent(inout) :: x + integer, intent(out) :: info + character, intent(in), optional :: iact + logical, intent(in), optional :: owned + end subroutine psb_glob_to_locs end interface interface psb_loc_to_glob @@ -508,8 +528,25 @@ Module psb_tools_mod integer, intent(out) :: info character, intent(in), optional :: iact end subroutine psb_loc_to_glob - end interface + subroutine psb_loc_to_glob2s(x,y,desc_a,info,iact) + use psb_descriptor_type + implicit none + type(psb_desc_type), intent(in) :: desc_a + integer,intent(in) :: x + integer,intent(out) :: y + integer, intent(out) :: info + character, intent(in), optional :: iact + + end subroutine psb_loc_to_glob2s + subroutine psb_loc_to_globs(x,desc_a,info,iact) + use psb_descriptor_type + type(psb_desc_type), intent(in) :: desc_a + integer,intent(inout) :: x + integer, intent(out) :: info + character, intent(in), optional :: iact + end subroutine psb_loc_to_globs + end interface interface psb_get_boundary module procedure psb_get_boundary @@ -606,6 +643,8 @@ contains goto 999 endif + desc_a%base_desc => null() + if (present(parts)) then if (.not.present(mg)) then info=581 @@ -678,14 +717,14 @@ contains subroutine psb_cdasb(desc_a,info) use psb_descriptor_type - interface - subroutine psb_icdasb(desc_a,info,ext_hv) - use psb_descriptor_type - Type(psb_desc_type), intent(inout) :: desc_a - integer, intent(out) :: info - logical, intent(in),optional :: ext_hv - end subroutine psb_icdasb - end interface + interface + subroutine psb_icdasb(desc_a,info,ext_hv) + use psb_descriptor_type + Type(psb_desc_type), intent(inout) :: desc_a + integer, intent(out) :: info + logical, intent(in),optional :: ext_hv + end subroutine psb_icdasb + end interface Type(psb_desc_type), intent(inout) :: desc_a integer, intent(out) :: info @@ -694,4 +733,5 @@ contains end subroutine psb_cdasb + end module psb_tools_mod diff --git a/base/modules/psi_mod.f90 b/base/modules/psi_mod.f90 index 9ac4a5c9..bfd9089c 100644 --- a/base/modules/psi_mod.f90 +++ b/base/modules/psi_mod.f90 @@ -28,9 +28,11 @@ !!$ POSSIBILITY OF SUCH DAMAGE. !!$ !!$ - module psi_mod + use psi_serial_mod + + interface subroutine psi_compute_size(desc_data,& & index_in, dl_lda, info) @@ -111,6 +113,20 @@ module psi_mod type(psb_desc_type), target :: desc_a integer, optional :: data end subroutine psi_dswapdatav + subroutine psi_dswapidxm(ictxt,icomm,flag,n,beta,y,idx,totxch,totsnd,totrcv,work,info) + integer, intent(in) :: ictxt,icomm,flag, n + integer, intent(out) :: info + real(kind(1.d0)) :: y(:,:), beta + real(kind(1.d0)),target :: work(:) + integer, intent(in) :: idx(:),totxch,totsnd,totrcv + end subroutine psi_dswapidxm + subroutine psi_dswapidxv(ictxt,icomm,flag,beta,y,idx,totxch,totsnd,totrcv,work,info) + integer, intent(in) :: ictxt,icomm,flag + integer, intent(out) :: info + real(kind(1.d0)) :: y(:), beta + real(kind(1.d0)),target :: work(:) + integer, intent(in) :: idx(:),totxch,totsnd,totrcv + end subroutine psi_dswapidxv subroutine psi_iswapdatam(flag,n,beta,y,desc_a,work,info,data) use psb_descriptor_type integer, intent(in) :: flag, n @@ -129,6 +145,20 @@ module psi_mod type(psb_desc_type), target :: desc_a integer, optional :: data end subroutine psi_iswapdatav + subroutine psi_iswapidxm(ictxt,icomm,flag,n,beta,y,idx,totxch,totsnd,totrcv,work,info) + integer, intent(in) :: ictxt,icomm,flag, n + integer, intent(out) :: info + integer :: y(:,:), beta + integer,target :: work(:) + integer, intent(in) :: idx(:),totxch,totsnd,totrcv + end subroutine psi_iswapidxm + subroutine psi_iswapidxv(ictxt,icomm,flag,beta,y,idx,totxch,totsnd,totrcv,work,info) + integer, intent(in) :: ictxt,icomm,flag + integer, intent(out) :: info + integer :: y(:), beta + integer,target :: work(:) + integer, intent(in) :: idx(:),totxch,totsnd,totrcv + end subroutine psi_iswapidxv subroutine psi_zswapdatam(flag,n,beta,y,desc_a,work,info,data) use psb_descriptor_type integer, intent(in) :: flag, n @@ -147,6 +177,20 @@ module psi_mod type(psb_desc_type), target :: desc_a integer, optional :: data end subroutine psi_zswapdatav + subroutine psi_zswapidxm(ictxt,icomm,flag,n,beta,y,idx,totxch,totsnd,totrcv,work,info) + integer, intent(in) :: ictxt,icomm,flag, n + integer, intent(out) :: info + complex(kind(1.d0)) :: y(:,:), beta + complex(kind(1.d0)),target :: work(:) + integer, intent(in) :: idx(:),totxch,totsnd,totrcv + end subroutine psi_zswapidxm + subroutine psi_zswapidxv(ictxt,icomm,flag,beta,y,idx,totxch,totsnd,totrcv,work,info) + integer, intent(in) :: ictxt,icomm,flag + integer, intent(out) :: info + complex(kind(1.d0)) :: y(:), beta + complex(kind(1.d0)),target :: work(:) + integer, intent(in) :: idx(:),totxch,totsnd,totrcv + end subroutine psi_zswapidxv end interface @@ -169,6 +213,20 @@ module psi_mod type(psb_desc_type), target :: desc_a integer, optional :: data end subroutine psi_dswaptranv + subroutine psi_dtranidxm(ictxt,icomm,flag,n,beta,y,idx,totxch,totsnd,totrcv,work,info) + integer, intent(in) :: ictxt,icomm,flag, n + integer, intent(out) :: info + real(kind(1.d0)) :: y(:,:), beta + real(kind(1.d0)),target :: work(:) + integer, intent(in) :: idx(:),totxch,totsnd,totrcv + end subroutine psi_dtranidxm + subroutine psi_dtranidxv(ictxt,icomm,flag,beta,y,idx,totxch,totsnd,totrcv,work,info) + integer, intent(in) :: ictxt,icomm,flag + integer, intent(out) :: info + real(kind(1.d0)) :: y(:), beta + real(kind(1.d0)),target :: work(:) + integer, intent(in) :: idx(:),totxch,totsnd,totrcv + end subroutine psi_dtranidxv subroutine psi_iswaptranm(flag,n,beta,y,desc_a,work,info,data) use psb_descriptor_type integer, intent(in) :: flag, n @@ -187,6 +245,20 @@ module psi_mod type(psb_desc_type), target :: desc_a integer, optional :: data end subroutine psi_iswaptranv + subroutine psi_itranidxm(ictxt,icomm,flag,n,beta,y,idx,totxch,totsnd,totrcv,work,info) + integer, intent(in) :: ictxt,icomm,flag, n + integer, intent(out) :: info + integer :: y(:,:), beta + integer, target :: work(:) + integer, intent(in) :: idx(:),totxch,totsnd,totrcv + end subroutine psi_itranidxm + subroutine psi_itranidxv(ictxt,icomm,flag,beta,y,idx,totxch,totsnd,totrcv,work,info) + integer, intent(in) :: ictxt,icomm,flag + integer, intent(out) :: info + integer :: y(:), beta + integer, target :: work(:) + integer, intent(in) :: idx(:),totxch,totsnd,totrcv + end subroutine psi_itranidxv subroutine psi_zswaptranm(flag,n,beta,y,desc_a,work,info,data) use psb_descriptor_type integer, intent(in) :: flag, n @@ -205,6 +277,20 @@ module psi_mod type(psb_desc_type), target :: desc_a integer, optional :: data end subroutine psi_zswaptranv + subroutine psi_ztranidxm(ictxt,icomm,flag,n,beta,y,idx,totxch,totsnd,totrcv,work,info) + integer, intent(in) :: ictxt,icomm,flag, n + integer, intent(out) :: info + complex(kind(1.d0)) :: y(:,:), beta + complex(kind(1.d0)), target :: work(:) + integer, intent(in) :: idx(:),totxch,totsnd,totrcv + end subroutine psi_ztranidxm + subroutine psi_ztranidxv(ictxt,icomm,flag,beta,y,idx,totxch,totsnd,totrcv,work,info) + integer, intent(in) :: ictxt,icomm,flag + integer, intent(out) :: info + complex(kind(1.d0)) :: y(:), beta + complex(kind(1.d0)), target :: work(:) + integer, intent(in) :: idx(:),totxch,totsnd,totrcv + end subroutine psi_ztranidxv end interface interface @@ -344,18 +430,6 @@ module psi_mod & psi_zovrl_restrr1, psi_zovrl_restrr2 end interface - interface psi_gth - module procedure psi_igthm, psi_igthv,& - & psi_dgthm, psi_dgthv,& - & psi_zgthm, psi_zgthv - end interface - - interface psi_sct - module procedure psi_isctm, psi_isctv,& - & psi_dsctm, psi_dsctv,& - & psi_zsctm, psi_zsctv - end interface - contains @@ -1677,319 +1751,319 @@ contains end subroutine psi_iovrl_restrr2 - subroutine psi_dgthm(n,k,idx,x,y) - - use psb_const_mod - implicit none - - integer :: n, k, idx(:) - real(kind(1.d0)) :: x(:,:), y(:) - - ! Locals - integer :: i, j, pt - - pt=0 - do j=1,k - do i=1,n - pt=pt+1 - y(pt)=x(idx(i),j) - end do - end do - - end subroutine psi_dgthm - - subroutine psi_dgthv(n,idx,x,y) - - use psb_const_mod - implicit none - - integer :: n, idx(:) - real(kind(1.d0)) :: x(:), y(:) - - ! Locals - integer :: i - - do i=1,n - y(i)=x(idx(i)) - end do - - end subroutine psi_dgthv - - - subroutine psi_dsctm(n,k,idx,x,beta,y) - - use psb_const_mod - implicit none - - integer :: n, k, idx(:) - real(kind(1.d0)) :: beta, x(:), y(:,:) - - ! Locals - integer :: i, j, pt - - if (beta == dzero) then - pt=0 - do j=1,k - do i=1,n - pt=pt+1 - y(idx(i),j) = x(pt) - end do - end do - else if (beta == done) then - pt=0 - do j=1,k - do i=1,n - pt=pt+1 - y(idx(i),j) = y(idx(i),j)+x(pt) - end do - end do - else - pt=0 - do j=1,k - do i=1,n - pt=pt+1 - y(idx(i),j) = beta*y(idx(i),j)+x(pt) - end do - end do - end if - end subroutine psi_dsctm - - subroutine psi_dsctv(n,idx,x,beta,y) - - use psb_const_mod - implicit none - - integer :: n, idx(:) - real(kind(1.d0)) :: beta, x(:), y(:) - - ! Locals - integer :: i - - if (beta == dzero) then - do i=1,n - y(idx(i)) = x(i) - end do - else if (beta == done) then - do i=1,n - y(idx(i)) = y(idx(i))+x(i) - end do - else - do i=1,n - y(idx(i)) = beta*y(idx(i))+x(i) - end do - end if - end subroutine psi_dsctv - - - subroutine psi_igthm(n,k,idx,x,y) - - use psb_const_mod - implicit none - - integer :: n, k, idx(:) - integer :: x(:,:), y(:) - - ! Locals - integer :: i, j, pt - - pt=0 - do j=1,k - do i=1,n - pt=pt+1 - y(pt)=x(idx(i),j) - end do - end do - - end subroutine psi_igthm - - - subroutine psi_igthv(n,idx,x,y) - - use psb_const_mod - implicit none - - integer :: n, idx(:) - integer :: x(:), y(:) - - ! Locals - integer :: i - - do i=1,n - y(i)=x(idx(i)) - end do - - end subroutine psi_igthv - - - - subroutine psi_isctm(n,k,idx,x,beta,y) - - use psb_const_mod - implicit none - - integer :: n, k, idx(:) - integer :: beta, x(:), y(:,:) - - ! Locals - integer :: i, j, pt - - if (beta == izero) then - pt=0 - do j=1,k - do i=1,n - pt=pt+1 - y(idx(i),j) = x(pt) - end do - end do - else if (beta == ione) then - pt=0 - do j=1,k - do i=1,n - pt=pt+1 - y(idx(i),j) = y(idx(i),j)+x(pt) - end do - end do - else - pt=0 - do j=1,k - do i=1,n - pt=pt+1 - y(idx(i),j) = beta*y(idx(i),j)+x(pt) - end do - end do - end if - end subroutine psi_isctm - - subroutine psi_isctv(n,idx,x,beta,y) - - use psb_const_mod - implicit none - - integer :: n, idx(:) - integer :: beta, x(:), y(:) - - ! Locals - integer :: i - - if (beta == izero) then - do i=1,n - y(idx(i)) = x(i) - end do - else if (beta == ione) then - do i=1,n - y(idx(i)) = y(idx(i))+x(i) - end do - else - do i=1,n - y(idx(i)) = beta*y(idx(i))+x(i) - end do - end if - end subroutine psi_isctv - - - subroutine psi_zgthm(n,k,idx,x,y) - - use psb_const_mod - implicit none - - integer :: n, k, idx(:) - complex(kind(1.d0)) :: x(:,:), y(:) - - ! Locals - integer :: i, j, pt - - pt=0 - do j=1,k - do i=1,n - pt=pt+1 - y(pt)=x(idx(i),j) - end do - end do - - end subroutine psi_zgthm - - - subroutine psi_zgthv(n,idx,x,y) - - use psb_const_mod - implicit none - - integer :: n, idx(:) - complex(kind(1.d0)) :: x(:), y(:) - - ! Locals - integer :: i - - do i=1,n - y(i)=x(idx(i)) - end do - - end subroutine psi_zgthv - - subroutine psi_zsctm(n,k,idx,x,beta,y) - - use psb_const_mod - implicit none - - integer :: n, k, idx(:) - complex(kind(1.d0)) :: beta, x(:), y(:,:) - - ! Locals - integer :: i, j, pt - - if (beta == zzero) then - pt=0 - do j=1,k - do i=1,n - pt=pt+1 - y(idx(i),j) = x(pt) - end do - end do - else if (beta == zone) then - pt=0 - do j=1,k - do i=1,n - pt=pt+1 - y(idx(i),j) = y(idx(i),j)+x(pt) - end do - end do - else - pt=0 - do j=1,k - do i=1,n - pt=pt+1 - y(idx(i),j) = beta*y(idx(i),j)+x(pt) - end do - end do - end if - end subroutine psi_zsctm - - - subroutine psi_zsctv(n,idx,x,beta,y) - - use psb_const_mod - implicit none - - integer :: n, idx(:) - complex(kind(1.d0)) :: beta, x(:), y(:) - - ! Locals - integer :: i - - if (beta == zzero) then - do i=1,n - y(idx(i)) = x(i) - end do - else if (beta == zone) then - do i=1,n - y(idx(i)) = y(idx(i))+x(i) - end do - else - do i=1,n - y(idx(i)) = beta*y(idx(i))+x(i) - end do - end if - end subroutine psi_zsctv +!!$ subroutine psi_dgthzm(n,k,idx,x,y) +!!$ +!!$ use psb_const_mod +!!$ implicit none +!!$ +!!$ integer :: n, k, idx(:) +!!$ real(kind(1.d0)) :: x(:,:), y(:) +!!$ +!!$ ! Locals +!!$ integer :: i, j, pt +!!$ +!!$ pt=0 +!!$ do j=1,k +!!$ do i=1,n +!!$ pt=pt+1 +!!$ y(pt)=x(idx(i),j) +!!$ end do +!!$ end do +!!$ +!!$ end subroutine psi_dgthzm +!!$ +!!$ subroutine psi_dgthzv(n,idx,x,y) +!!$ +!!$ use psb_const_mod +!!$ implicit none +!!$ +!!$ integer :: n, idx(:) +!!$ real(kind(1.d0)) :: x(:), y(:) +!!$ +!!$ ! Locals +!!$ integer :: i +!!$ +!!$ do i=1,n +!!$ y(i)=x(idx(i)) +!!$ end do +!!$ +!!$ end subroutine psi_dgthzv +!!$ +!!$ +!!$ subroutine psi_dsctm(n,k,idx,x,beta,y) +!!$ +!!$ use psb_const_mod +!!$ implicit none +!!$ +!!$ integer :: n, k, idx(:) +!!$ real(kind(1.d0)) :: beta, x(:), y(:,:) +!!$ +!!$ ! Locals +!!$ integer :: i, j, pt +!!$ +!!$ if (beta == dzero) then +!!$ pt=0 +!!$ do j=1,k +!!$ do i=1,n +!!$ pt=pt+1 +!!$ y(idx(i),j) = x(pt) +!!$ end do +!!$ end do +!!$ else if (beta == done) then +!!$ pt=0 +!!$ do j=1,k +!!$ do i=1,n +!!$ pt=pt+1 +!!$ y(idx(i),j) = y(idx(i),j)+x(pt) +!!$ end do +!!$ end do +!!$ else +!!$ pt=0 +!!$ do j=1,k +!!$ do i=1,n +!!$ pt=pt+1 +!!$ y(idx(i),j) = beta*y(idx(i),j)+x(pt) +!!$ end do +!!$ end do +!!$ end if +!!$ end subroutine psi_dsctm +!!$ +!!$ subroutine psi_dsctv(n,idx,x,beta,y) +!!$ +!!$ use psb_const_mod +!!$ implicit none +!!$ +!!$ integer :: n, idx(:) +!!$ real(kind(1.d0)) :: beta, x(:), y(:) +!!$ +!!$ ! Locals +!!$ integer :: i +!!$ +!!$ if (beta == dzero) then +!!$ do i=1,n +!!$ y(idx(i)) = x(i) +!!$ end do +!!$ else if (beta == done) then +!!$ do i=1,n +!!$ y(idx(i)) = y(idx(i))+x(i) +!!$ end do +!!$ else +!!$ do i=1,n +!!$ y(idx(i)) = beta*y(idx(i))+x(i) +!!$ end do +!!$ end if +!!$ end subroutine psi_dsctv +!!$ +!!$ +!!$ subroutine psi_igthzm(n,k,idx,x,y) +!!$ +!!$ use psb_const_mod +!!$ implicit none +!!$ +!!$ integer :: n, k, idx(:) +!!$ integer :: x(:,:), y(:) +!!$ +!!$ ! Locals +!!$ integer :: i, j, pt +!!$ +!!$ pt=0 +!!$ do j=1,k +!!$ do i=1,n +!!$ pt=pt+1 +!!$ y(pt)=x(idx(i),j) +!!$ end do +!!$ end do +!!$ +!!$ end subroutine psi_igthzm +!!$ +!!$ +!!$ subroutine psi_igthzv(n,idx,x,y) +!!$ +!!$ use psb_const_mod +!!$ implicit none +!!$ +!!$ integer :: n, idx(:) +!!$ integer :: x(:), y(:) +!!$ +!!$ ! Locals +!!$ integer :: i +!!$ +!!$ do i=1,n +!!$ y(i)=x(idx(i)) +!!$ end do +!!$ +!!$ end subroutine psi_igthzv +!!$ +!!$ +!!$ +!!$ subroutine psi_isctm(n,k,idx,x,beta,y) +!!$ +!!$ use psb_const_mod +!!$ implicit none +!!$ +!!$ integer :: n, k, idx(:) +!!$ integer :: beta, x(:), y(:,:) +!!$ +!!$ ! Locals +!!$ integer :: i, j, pt +!!$ +!!$ if (beta == izero) then +!!$ pt=0 +!!$ do j=1,k +!!$ do i=1,n +!!$ pt=pt+1 +!!$ y(idx(i),j) = x(pt) +!!$ end do +!!$ end do +!!$ else if (beta == ione) then +!!$ pt=0 +!!$ do j=1,k +!!$ do i=1,n +!!$ pt=pt+1 +!!$ y(idx(i),j) = y(idx(i),j)+x(pt) +!!$ end do +!!$ end do +!!$ else +!!$ pt=0 +!!$ do j=1,k +!!$ do i=1,n +!!$ pt=pt+1 +!!$ y(idx(i),j) = beta*y(idx(i),j)+x(pt) +!!$ end do +!!$ end do +!!$ end if +!!$ end subroutine psi_isctm +!!$ +!!$ subroutine psi_isctv(n,idx,x,beta,y) +!!$ +!!$ use psb_const_mod +!!$ implicit none +!!$ +!!$ integer :: n, idx(:) +!!$ integer :: beta, x(:), y(:) +!!$ +!!$ ! Locals +!!$ integer :: i +!!$ +!!$ if (beta == izero) then +!!$ do i=1,n +!!$ y(idx(i)) = x(i) +!!$ end do +!!$ else if (beta == ione) then +!!$ do i=1,n +!!$ y(idx(i)) = y(idx(i))+x(i) +!!$ end do +!!$ else +!!$ do i=1,n +!!$ y(idx(i)) = beta*y(idx(i))+x(i) +!!$ end do +!!$ end if +!!$ end subroutine psi_isctv +!!$ +!!$ +!!$ subroutine psi_zgthzm(n,k,idx,x,y) +!!$ +!!$ use psb_const_mod +!!$ implicit none +!!$ +!!$ integer :: n, k, idx(:) +!!$ complex(kind(1.d0)) :: x(:,:), y(:) +!!$ +!!$ ! Locals +!!$ integer :: i, j, pt +!!$ +!!$ pt=0 +!!$ do j=1,k +!!$ do i=1,n +!!$ pt=pt+1 +!!$ y(pt)=x(idx(i),j) +!!$ end do +!!$ end do +!!$ +!!$ end subroutine psi_zgthzm +!!$ +!!$ +!!$ subroutine psi_zgthzv(n,idx,x,y) +!!$ +!!$ use psb_const_mod +!!$ implicit none +!!$ +!!$ integer :: n, idx(:) +!!$ complex(kind(1.d0)) :: x(:), y(:) +!!$ +!!$ ! Locals +!!$ integer :: i +!!$ +!!$ do i=1,n +!!$ y(i)=x(idx(i)) +!!$ end do +!!$ +!!$ end subroutine psi_zgthzv +!!$ +!!$ subroutine psi_zsctm(n,k,idx,x,beta,y) +!!$ +!!$ use psb_const_mod +!!$ implicit none +!!$ +!!$ integer :: n, k, idx(:) +!!$ complex(kind(1.d0)) :: beta, x(:), y(:,:) +!!$ +!!$ ! Locals +!!$ integer :: i, j, pt +!!$ +!!$ if (beta == zzero) then +!!$ pt=0 +!!$ do j=1,k +!!$ do i=1,n +!!$ pt=pt+1 +!!$ y(idx(i),j) = x(pt) +!!$ end do +!!$ end do +!!$ else if (beta == zone) then +!!$ pt=0 +!!$ do j=1,k +!!$ do i=1,n +!!$ pt=pt+1 +!!$ y(idx(i),j) = y(idx(i),j)+x(pt) +!!$ end do +!!$ end do +!!$ else +!!$ pt=0 +!!$ do j=1,k +!!$ do i=1,n +!!$ pt=pt+1 +!!$ y(idx(i),j) = beta*y(idx(i),j)+x(pt) +!!$ end do +!!$ end do +!!$ end if +!!$ end subroutine psi_zsctm +!!$ +!!$ +!!$ subroutine psi_zsctv(n,idx,x,beta,y) +!!$ +!!$ use psb_const_mod +!!$ implicit none +!!$ +!!$ integer :: n, idx(:) +!!$ complex(kind(1.d0)) :: beta, x(:), y(:) +!!$ +!!$ ! Locals +!!$ integer :: i +!!$ +!!$ if (beta == zzero) then +!!$ do i=1,n +!!$ y(idx(i)) = x(i) +!!$ end do +!!$ else if (beta == zone) then +!!$ do i=1,n +!!$ y(idx(i)) = y(idx(i))+x(i) +!!$ end do +!!$ else +!!$ do i=1,n +!!$ y(idx(i)) = beta*y(idx(i))+x(i) +!!$ end do +!!$ end if +!!$ end subroutine psi_zsctv subroutine psi_bld_ovr_mst(me,ovrlap_elem,mst_idx,info) use psb_const_mod diff --git a/base/modules/psi_serial_mod.f90 b/base/modules/psi_serial_mod.f90 new file mode 100644 index 00000000..6d7b4dd5 --- /dev/null +++ b/base/modules/psi_serial_mod.f90 @@ -0,0 +1,518 @@ +!!$ +!!$ Parallel Sparse BLAS v2.0 +!!$ (C) Copyright 2006 Salvatore Filippone University of Rome Tor Vergata +!!$ Alfredo Buttari University of Rome Tor Vergata +!!$ +!!$ 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. +!!$ +!!$ +module psi_serial_mod + + + interface psi_gth + module procedure & + & psi_igthv, psi_dgthv, psi_zgthv,& + & psi_igthzv, psi_dgthzv, psi_zgthzv,& + & psi_igthzmv, psi_dgthzmv, psi_zgthzmv + end interface + + interface psi_sct + module procedure psi_isctmv, psi_isctv,& + & psi_dsctmv, psi_dsctv,& + & psi_zsctmv, psi_zsctv + end interface + + +contains + + + subroutine psi_igthv(n,idx,alpha,x,beta,y) + + use psb_const_mod + implicit none + + integer :: n, idx(:) + integer :: x(:), y(:), alpha, beta + + ! Locals + integer :: i + if (beta == izero) then + if (alpha == ione) then + do i=1,n + y(i) = x(idx(i)) + end do + else if (alpha == -ione) then + do i=1,n + y(i) = -x(idx(i)) + end do + else + do i=1,n + y(i) = alpha*x(idx(i)) + end do + end if + else + if (beta == ione) then + ! Do nothing + else if (beta == -ione) then + y(1:n) = -y(1:n) + else + y(1:n) = beta*y(1:n) + end if + + if (alpha == ione) then + do i=1,n + y(i) = y(i) + x(idx(i)) + end do + else if (alpha == -ione) then + do i=1,n + y(i) = y(i) - x(idx(i)) + end do + else + do i=1,n + y(i) = y(i) + alpha*x(idx(i)) + end do + end if + end if + + end subroutine psi_igthv + + subroutine psi_dgthv(n,idx,alpha,x,beta,y) + + use psb_const_mod + implicit none + + integer :: n, idx(:) + real(kind(1.d0)) :: x(:), y(:), alpha, beta + + ! Locals + integer :: i + if (beta == dzero) then + if (alpha == done) then + do i=1,n + y(i) = x(idx(i)) + end do + else if (alpha == -done) then + do i=1,n + y(i) = -x(idx(i)) + end do + else + do i=1,n + y(i) = alpha*x(idx(i)) + end do + end if + else + if (beta == done) then + ! Do nothing + else if (beta == -done) then + y(1:n) = -y(1:n) + else + y(1:n) = beta*y(1:n) + end if + + if (alpha == done) then + do i=1,n + y(i) = y(i) + x(idx(i)) + end do + else if (alpha == -done) then + do i=1,n + y(i) = y(i) - x(idx(i)) + end do + else + do i=1,n + y(i) = y(i) + alpha*x(idx(i)) + end do + end if + end if + + end subroutine psi_dgthv + + subroutine psi_zgthv(n,idx,alpha,x,beta,y) + + use psb_const_mod + implicit none + + integer :: n, idx(:) + complex(kind(1.d0)) :: x(:), y(:),alpha,beta + + ! Locals + integer :: i + if (beta == zzero) then + if (alpha == zone) then + do i=1,n + y(i) = x(idx(i)) + end do + else if (alpha == -zone) then + do i=1,n + y(i) = -x(idx(i)) + end do + else + do i=1,n + y(i) = alpha*x(idx(i)) + end do + end if + else + if (beta == zone) then + ! Do nothing + else if (beta == -zone) then + y(1:n) = -y(1:n) + else + y(1:n) = beta*y(1:n) + end if + + if (alpha == zone) then + do i=1,n + y(i) = y(i) + x(idx(i)) + end do + else if (alpha == -zone) then + do i=1,n + y(i) = y(i) - x(idx(i)) + end do + else + do i=1,n + y(i) = y(i) + alpha*x(idx(i)) + end do + end if + end if + + end subroutine psi_zgthv + + + + + subroutine psi_dgthzmv(n,k,idx,x,y) + + use psb_const_mod + implicit none + + integer :: n, k, idx(:) + real(kind(1.d0)) :: x(:,:), y(:) + + ! Locals + integer :: i, j, pt + + pt=0 + do j=1,k + do i=1,n + pt=pt+1 + y(pt)=x(idx(i),j) + end do + end do + + end subroutine psi_dgthzmv + + + subroutine psi_igthzmv(n,k,idx,x,y) + + use psb_const_mod + implicit none + + integer :: n, k, idx(:) + integer :: x(:,:), y(:) + + ! Locals + integer :: i, j, pt + + pt=0 + do j=1,k + do i=1,n + pt=pt+1 + y(pt)=x(idx(i),j) + end do + end do + + end subroutine psi_igthzmv + + + subroutine psi_zgthzmv(n,k,idx,x,y) + + use psb_const_mod + implicit none + + integer :: n, k, idx(:) + complex(kind(1.d0)) :: x(:,:), y(:) + + ! Locals + integer :: i, j, pt + + pt=0 + do j=1,k + do i=1,n + pt=pt+1 + y(pt)=x(idx(i),j) + end do + end do + + end subroutine psi_zgthzmv + + subroutine psi_dgthzv(n,idx,x,y) + + use psb_const_mod + implicit none + + integer :: n, idx(:) + real(kind(1.d0)) :: x(:), y(:) + + ! Locals + integer :: i + + do i=1,n + y(i)=x(idx(i)) + end do + + end subroutine psi_dgthzv + + subroutine psi_igthzv(n,idx,x,y) + + use psb_const_mod + implicit none + + integer :: n, idx(:) + integer :: x(:), y(:) + + ! Locals + integer :: i + + do i=1,n + y(i)=x(idx(i)) + end do + + end subroutine psi_igthzv + + subroutine psi_zgthzv(n,idx,x,y) + + use psb_const_mod + implicit none + + integer :: n, idx(:) + complex(kind(1.d0)) :: x(:), y(:) + + ! Locals + integer :: i + + do i=1,n + y(i)=x(idx(i)) + end do + + end subroutine psi_zgthzv + + + subroutine psi_dsctmv(n,k,idx,x,beta,y) + + use psb_const_mod + implicit none + + integer :: n, k, idx(:) + real(kind(1.d0)) :: beta, x(:), y(:,:) + + ! Locals + integer :: i, j, pt + + if (beta == dzero) then + pt=0 + do j=1,k + do i=1,n + pt=pt+1 + y(idx(i),j) = x(pt) + end do + end do + else if (beta == done) then + pt=0 + do j=1,k + do i=1,n + pt=pt+1 + y(idx(i),j) = y(idx(i),j)+x(pt) + end do + end do + else + pt=0 + do j=1,k + do i=1,n + pt=pt+1 + y(idx(i),j) = beta*y(idx(i),j)+x(pt) + end do + end do + end if + end subroutine psi_dsctmv + + subroutine psi_dsctv(n,idx,x,beta,y) + + use psb_const_mod + implicit none + + integer :: n, idx(:) + real(kind(1.d0)) :: beta, x(:), y(:) + + ! Locals + integer :: i + + if (beta == dzero) then + do i=1,n + y(idx(i)) = x(i) + end do + else if (beta == done) then + do i=1,n + y(idx(i)) = y(idx(i))+x(i) + end do + else + do i=1,n + y(idx(i)) = beta*y(idx(i)) + end do + do i=1,n + y(idx(i)) = y(idx(i))+x(i) + end do + end if + end subroutine psi_dsctv + + subroutine psi_isctmv(n,k,idx,x,beta,y) + + use psb_const_mod + implicit none + + integer :: n, k, idx(:) + integer :: beta, x(:), y(:,:) + + ! Locals + integer :: i, j, pt + + if (beta == izero) then + pt=0 + do j=1,k + do i=1,n + pt=pt+1 + y(idx(i),j) = x(pt) + end do + end do + else if (beta == ione) then + pt=0 + do j=1,k + do i=1,n + pt=pt+1 + y(idx(i),j) = y(idx(i),j)+x(pt) + end do + end do + else + pt=0 + do j=1,k + do i=1,n + pt=pt+1 + y(idx(i),j) = beta*y(idx(i),j)+x(pt) + end do + end do + end if + end subroutine psi_isctmv + + subroutine psi_isctv(n,idx,x,beta,y) + + use psb_const_mod + implicit none + + integer :: n, idx(:) + integer :: beta, x(:), y(:) + + ! Locals + integer :: i + + if (beta == izero) then + do i=1,n + y(idx(i)) = x(i) + end do + else if (beta == ione) then + do i=1,n + y(idx(i)) = y(idx(i))+x(i) + end do + else + do i=1,n + y(idx(i)) = beta*y(idx(i))+x(i) + end do + end if + end subroutine psi_isctv + + subroutine psi_zsctmv(n,k,idx,x,beta,y) + + use psb_const_mod + implicit none + + integer :: n, k, idx(:) + complex(kind(1.d0)) :: beta, x(:), y(:,:) + + ! Locals + integer :: i, j, pt + + if (beta == zzero) then + pt=0 + do j=1,k + do i=1,n + pt=pt+1 + y(idx(i),j) = x(pt) + end do + end do + else if (beta == zone) then + pt=0 + do j=1,k + do i=1,n + pt=pt+1 + y(idx(i),j) = y(idx(i),j)+x(pt) + end do + end do + else + pt=0 + do j=1,k + do i=1,n + pt=pt+1 + y(idx(i),j) = beta*y(idx(i),j)+x(pt) + end do + end do + end if + end subroutine psi_zsctmv + + + subroutine psi_zsctv(n,idx,x,beta,y) + + use psb_const_mod + implicit none + + integer :: n, idx(:) + complex(kind(1.d0)) :: beta, x(:), y(:) + + ! Locals + integer :: i + + if (beta == zzero) then + do i=1,n + y(idx(i)) = x(i) + end do + else if (beta == zone) then + do i=1,n + y(idx(i)) = y(idx(i))+x(i) + end do + else + do i=1,n + y(idx(i)) = beta*y(idx(i))+x(i) + end do + end if + end subroutine psi_zsctv + + +end module psi_serial_mod diff --git a/base/tools/psb_cdcpy.f90 b/base/tools/psb_cdcpy.f90 index 5fff9a6a..0090f437 100644 --- a/base/tools/psb_cdcpy.f90 +++ b/base/tools/psb_cdcpy.f90 @@ -78,19 +78,19 @@ subroutine psb_cdcpy(desc_in, desc_out, info) goto 9999 endif - call psb_safe_cpy(desc_in%matrix_data,desc_out%matrix_data,info) - if (info == 0) call psb_safe_cpy(desc_in%halo_index,desc_out%halo_index,info) - if (info == 0) call psb_safe_cpy(desc_in%ext_index,desc_out%ext_index,info) - if (info == 0) call psb_safe_cpy(desc_in%ovrlap_index,desc_out%ovrlap_index,info) - if (info == 0) call psb_safe_cpy(desc_in%bnd_elem,desc_out%bnd_elem,info) - if (info == 0) call psb_safe_cpy(desc_in%ovrlap_elem,desc_out%ovrlap_elem,info) - if (info == 0) call psb_safe_cpy(desc_in%ovr_mst_idx,desc_out%ovr_mst_idx,info) - if (info == 0) call psb_safe_cpy(desc_in%loc_to_glob,desc_out%loc_to_glob,info) - if (info == 0) call psb_safe_cpy(desc_in%glob_to_loc,desc_out%glob_to_loc,info) - if (info == 0) call psb_safe_cpy(desc_in%lprm,desc_out%lprm,info) - if (info == 0) call psb_safe_cpy(desc_in%idx_space,desc_out%idx_space,info) - if (info == 0) call psb_safe_cpy(desc_in%hashv,desc_out%hashv,info) - if (info == 0) call psb_safe_cpy(desc_in%glb_lc,desc_out%glb_lc,info) + call psb_safe_ab_cpy(desc_in%matrix_data,desc_out%matrix_data,info) + if (info == 0) call psb_safe_ab_cpy(desc_in%halo_index,desc_out%halo_index,info) + if (info == 0) call psb_safe_ab_cpy(desc_in%ext_index,desc_out%ext_index,info) + if (info == 0) call psb_safe_ab_cpy(desc_in%ovrlap_index,desc_out%ovrlap_index,info) + if (info == 0) call psb_safe_ab_cpy(desc_in%bnd_elem,desc_out%bnd_elem,info) + if (info == 0) call psb_safe_ab_cpy(desc_in%ovrlap_elem,desc_out%ovrlap_elem,info) + if (info == 0) call psb_safe_ab_cpy(desc_in%ovr_mst_idx,desc_out%ovr_mst_idx,info) + if (info == 0) call psb_safe_ab_cpy(desc_in%loc_to_glob,desc_out%loc_to_glob,info) + if (info == 0) call psb_safe_ab_cpy(desc_in%glob_to_loc,desc_out%glob_to_loc,info) + if (info == 0) call psb_safe_ab_cpy(desc_in%lprm,desc_out%lprm,info) + if (info == 0) call psb_safe_ab_cpy(desc_in%idx_space,desc_out%idx_space,info) + if (info == 0) call psb_safe_ab_cpy(desc_in%hashv,desc_out%hashv,info) + if (info == 0) call psb_safe_ab_cpy(desc_in%glb_lc,desc_out%glb_lc,info) if (info == 0) then if (allocated(desc_in%ptree)) then diff --git a/base/tools/psb_dcdovr.F90 b/base/tools/psb_dcdovr.F90 index 9f909a20..0e52f22f 100644 --- a/base/tools/psb_dcdovr.F90 +++ b/base/tools/psb_dcdovr.F90 @@ -78,11 +78,11 @@ Subroutine psb_dcdovr(a,desc_a,novr,desc_ov,info, extype) ! .. Array Arguments .. integer, intent(in) :: novr - Type(psb_dspmat_type), Intent(in) :: a - Type(psb_desc_type), Intent(in) :: desc_a - Type(psb_desc_type), Intent(out) :: desc_ov - integer, intent(out) :: info - integer, intent(in),optional :: extype + Type(psb_dspmat_type), Intent(in) :: a + Type(psb_desc_type), Intent(in), target :: desc_a + Type(psb_desc_type), Intent(out) :: desc_ov + integer, intent(out) :: info + integer, intent(in),optional :: extype interface subroutine psb_icdasb(desc_a,info,ext_hv) @@ -192,7 +192,8 @@ Subroutine psb_dcdovr(a,desc_a,novr,desc_ov,info, extype) l_tmp_ovr_idx = novr*(3*Max(2*index_dim,1)+1) l_tmp_halo = novr*(3*Size(desc_a%halo_index)) - call psb_cd_set_bld(desc_ov,info) + call psb_cd_set_ovl_bld(desc_ov,info) + desc_ov%base_desc => desc_a If (debug_level >= psb_debug_outer_) then Write(debug_unit,*) me,' ',trim(name),':Start',lworks,lworkr @@ -229,16 +230,16 @@ Subroutine psb_dcdovr(a,desc_a,novr,desc_ov,info, extype) call psb_errpush(4010,name,a_err='Allocate') goto 9999 end if - halo(:) = desc_a%halo_index(:) - tmp_ovr_idx(:) = -1 - orig_ovr(:) = -1 - tmp_halo(:) = -1 - counter_e = 1 - tot_recv = 0 - counter_t = 1 - counter_h = 1 - counter_o = 1 - cntov_o = 1 + halo(:) = desc_a%halo_index(:) + tmp_ovr_idx(:) = -1 + orig_ovr(:) = -1 + tmp_halo(:) = -1 + counter_e = 1 + tot_recv = 0 + counter_t = 1 + counter_h = 1 + counter_o = 1 + cntov_o = 1 ! Init overlap with desc_a%ovrlap (if any) counter = 1 Do While (desc_a%ovrlap_index(counter) /= -1) @@ -739,6 +740,7 @@ Subroutine psb_dcdovr(a,desc_a,novr,desc_ov,info, extype) call psb_icdasb(desc_ov,info,ext_hv=.true.) + call psb_cd_set_ovl_asb(desc_ov,info) if (info == 0) call psb_sp_free(blk,info) if (info /= 0) then diff --git a/base/tools/psb_glob_to_loc.f90 b/base/tools/psb_glob_to_loc.f90 index b45ad262..27a6808c 100644 --- a/base/tools/psb_glob_to_loc.f90 +++ b/base/tools/psb_glob_to_loc.f90 @@ -249,3 +249,38 @@ subroutine psb_glob_to_loc(x,desc_a,info,iact,owned) end subroutine psb_glob_to_loc +subroutine psb_glob_to_loc2s(x,y,desc_a,info,iact,owned) + use psb_descriptor_type + use psb_tools_mod, psb_protect_name => psb_glob_to_loc2s + implicit none + type(psb_desc_type), intent(in) :: desc_a + integer,intent(in) :: x + integer,intent(out) :: y + integer, intent(out) :: info + character, intent(in), optional :: iact + logical, intent(in), optional :: owned + + integer :: iv1(1), iv2(1) + + iv1(1) = x + call psb_glob_to_loc(iv1,iv2,desc_a,info,iact,owned) + y = iv2(1) +end subroutine psb_glob_to_loc2s + +subroutine psb_glob_to_locs(x,desc_a,info,iact,owned) + use psb_descriptor_type + use psb_tools_mod, psb_protect_name => psb_glob_to_locs + implicit none + type(psb_desc_type), intent(in) :: desc_a + integer,intent(inout) :: x + integer, intent(out) :: info + character, intent(in), optional :: iact + logical, intent(in), optional :: owned + integer :: iv1(1) + + iv1(1) = x + call psb_glob_to_loc(iv1,desc_a,info,iact,owned) + x = iv1(1) + +end subroutine psb_glob_to_locs + diff --git a/base/tools/psb_loc_to_glob.f90 b/base/tools/psb_loc_to_glob.f90 index 25d01faf..1fe493aa 100644 --- a/base/tools/psb_loc_to_glob.f90 +++ b/base/tools/psb_loc_to_glob.f90 @@ -249,3 +249,36 @@ subroutine psb_loc_to_glob(x,desc_a,info,iact) end subroutine psb_loc_to_glob +subroutine psb_loc_to_glob2s(x,y,desc_a,info,iact) + use psb_descriptor_type + use psb_tools_mod, psb_protect_name => psb_loc_to_glob2s + implicit none + type(psb_desc_type), intent(in) :: desc_a + integer,intent(in) :: x + integer,intent(out) :: y + integer, intent(out) :: info + character, intent(in), optional :: iact + + integer :: iv1(1), iv2(1) + + iv1(1) = x + call psb_loc_to_glob(iv1,iv2,desc_a,info,iact) + y = iv2(1) +end subroutine psb_loc_to_glob2s + +subroutine psb_loc_to_globs(x,desc_a,info,iact) + use psb_descriptor_type + use psb_tools_mod, psb_protect_name => psb_loc_to_globs + implicit none + type(psb_desc_type), intent(in) :: desc_a + integer,intent(inout) :: x + integer, intent(out) :: info + character, intent(in), optional :: iact + integer :: iv1(1) + + iv1(1) = x + call psb_loc_to_glob(iv1,desc_a,info,iact) + x = iv1(1) + +end subroutine psb_loc_to_globs + diff --git a/base/tools/psb_zcdovr.F90 b/base/tools/psb_zcdovr.F90 index 0093f75b..7116e6b8 100644 --- a/base/tools/psb_zcdovr.F90 +++ b/base/tools/psb_zcdovr.F90 @@ -77,11 +77,11 @@ Subroutine psb_zcdovr(a,desc_a,novr,desc_ov,info, extype) ! .. Array Arguments .. integer, intent(in) :: novr - Type(psb_zspmat_type), Intent(in) :: a - Type(psb_desc_type), Intent(in) :: desc_a - Type(psb_desc_type), Intent(out) :: desc_ov - integer, intent(out) :: info - integer, intent(in),optional :: extype + Type(psb_zspmat_type), Intent(in) :: a + Type(psb_desc_type), Intent(in), target :: desc_a + Type(psb_desc_type), Intent(out) :: desc_ov + integer, intent(out) :: info + integer, intent(in),optional :: extype interface subroutine psb_icdasb(desc_a,info,ext_hv) @@ -191,7 +191,8 @@ Subroutine psb_zcdovr(a,desc_a,novr,desc_ov,info, extype) l_tmp_ovr_idx = novr*(3*Max(2*index_dim,1)+1) l_tmp_halo = novr*(3*Size(desc_a%halo_index)) - call psb_cd_set_bld(desc_ov,info) + call psb_cd_set_ovl_bld(desc_ov,info) + desc_ov%base_desc => desc_a If (debug_level >= psb_debug_outer_) then Write(debug_unit,*) me,' ',trim(name),':Start',lworks,lworkr @@ -386,7 +387,8 @@ Subroutine psb_zcdovr(a,desc_a,novr,desc_ov,info, extype) Enddo if (debug_level >= psb_debug_outer_) & - & write(debug_unit,*) me,' ',trim(name),':Checktmp_o_i Loop Mid1',tmp_ovr_idx(1:10) + & write(debug_unit,*) me,' ',trim(name),& + & ':Checktmp_o_i Loop Mid1',tmp_ovr_idx(1:10) counter = counter+n_elem_recv ! @@ -460,16 +462,19 @@ Subroutine psb_zcdovr(a,desc_a,novr,desc_ov,info, extype) tot_elem=i endif if (debug_level >= psb_debug_outer_) & - & write(debug_unit,*) me,' ',trim(name),':Checktmp_o_i Loop Mid2',tmp_ovr_idx(1:10) + & write(debug_unit,*) me,' ',trim(name),& + & ':Checktmp_o_i Loop Mid2',tmp_ovr_idx(1:10) sdsz(proc+1) = tot_elem idxs = idxs + tot_elem end if counter = counter+n_elem_send+3 if (debug_level >= psb_debug_outer_) & - & write(debug_unit,*) me,' ',trim(name),':Checktmp_o_i Loop End',tmp_ovr_idx(1:10) + & write(debug_unit,*) me,' ',trim(name),& + & ':Checktmp_o_i Loop End',tmp_ovr_idx(1:10) Enddo if (debug_level >= psb_debug_outer_) & - & write(debug_unit,*) me,' ',trim(name),':End phase 1', m, n_col, tot_recv + & write(debug_unit,*) me,' ',trim(name),& + & ':End phase 1', m, n_col, tot_recv if (i_ovr <= novr) then ! @@ -734,6 +739,7 @@ Subroutine psb_zcdovr(a,desc_a,novr,desc_ov,info, extype) call psb_icdasb(desc_ov,info,ext_hv=.true.) + call psb_cd_set_ovl_asb(desc_ov,info) if (info == 0) call psb_sp_free(blk,info) if (info /= 0) then