! debug out-of-bound
! #define DBG_OOB
! debug personal
! #define DBG_PER

!
!                Parallel Sparse BLAS  version 3.5
!      (C) Copyright 2006-2018
!        Salvatore Filippone
!        Alfredo Buttari
!
!    Redistribution and use in source and binary forms, with or without
!    modification, are permitted provided that the following conditions
!    are met:
!      1. Redistributions of source code must retain the above copyright
!         notice, this list of conditions and the following disclaimer.
!      2. Redistributions in binary form must reproduce the above copyright
!         notice, this list of conditions, and the following disclaimer in the
!         documentation and/or other materials provided with the distribution.
!      3. The name of the PSBLAS group or the names of its contributors may
!         not be used to endorse or promote products derived from this
!         software without specific written permission.
!
!    THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
!    ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
!    TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
!    PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE PSBLAS GROUP OR ITS CONTRIBUTORS
!    BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
!    CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
!    SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
!    INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
!    CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
!    ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
!    POSSIBILITY OF SUCH DAMAGE.
!
!
!
! File: psi_dswapdata.F90
!
! Subroutine: psi_dswapdatam
!   Does the data exchange among processes. Essentially this is doing
!   a variable all-to-all data exchange (ALLTOALLV in MPI parlance), but
!   it is capable of pruning empty exchanges, which are very likely in out
!   application environment. All the variants have the same structure
!   In all these subroutines X may be:    I    Integer
!                                         S    real(psb_spk_)
!                                         D    real(psb_dpk_)
!                                         C    complex(psb_spk_)
!                                         Z    complex(psb_dpk_)
!   Basically the operation is as follows: on each process, we identify
!   sections SND(Y) and RCV(Y); then we do a send on (PACK(SND(Y)));
!   then we receive, and we do an update with Y = UNPACK(RCV(Y)) + BETA * Y
!   but only on the elements involved in the UNPACK operation.
!   Thus: for halo data exchange, the receive section is confined in the
!   halo indices, and BETA=0, whereas for overlap exchange the receive section
!   is scattered in the owned indices, and BETA=1.
!
! Arguments:
!    flag     - integer                 Choose the algorithm for data exchange:
!                                       this is chosen through bit fields.
!                                        swap_mpi  = iand(flag,psb_swap_mpi_)  /= 0
!                                        swap_sync = iand(flag,psb_swap_sync_) /= 0
!                                        swap_send = iand(flag,psb_swap_send_) /= 0
!                                        swap_recv = iand(flag,psb_swap_recv_) /= 0
!                                       if (swap_mpi):  use underlying MPI_ALLTOALLV.
!                                       if (swap_sync): use PSB_SND and PSB_RCV in
!                                                       synchronized pairs
!                                       if (swap_send .and. swap_recv): use mpi_irecv
!                                                       and mpi_send
!                                       if (swap_send): use psb_snd (but need another
!                                                       call with swap_recv to complete)
!                                       if (swap_recv): use psb_rcv (completing a
!                                                       previous call with swap_send)
!
!
!    n        - integer                 Number of columns in Y
!    beta     - X                       Choose overwrite or sum.
!    y(:,:)   - X                       The data area
!    desc_a   - type(psb_desc_type).  The communication descriptor.
!    work(:)  - X                       Buffer space. If not sufficient, will do
!                                       our own internal allocation.
!    info     - integer.                return code.
!    data     - integer                 which list is to be used to exchange data
!                                       default psb_comm_halo_
!                                       psb_comm_halo_    use halo_index
!                                       psb_comm_ext_     use ext_index
!                                       psb_comm_ovrl_    use ovrl_index
!                                       psb_comm_mov_     use ovr_mst_idx
!
!
!
!
! Subroutine: psi_dswapdata_vect
!   Data exchange among processes.
!
!   Takes care of Y an exanspulated vector.
!
!
!
subroutine psi_dswapdata_vect(flag,beta,y,desc_a,work,info,data)

  use psi_mod, psb_protect_name => psi_dswapdata_vect
  use psb_d_base_vect_mod
  use psb_error_mod
  use psb_desc_mod
  use psb_penv_mod
#ifdef MPI_MOD
  use mpi
#endif
  implicit none
#ifdef MPI_H
  include 'mpif.h'
#endif
#include "scorep/SCOREP_User.inc"

  integer(psb_ipk_), intent(in)  :: flag
  integer(psb_ipk_), intent(out) :: info
  class(psb_d_base_vect_type)    :: y
  real(psb_dpk_)                 :: beta
  real(psb_dpk_), target         :: work(:)
  type(psb_desc_type), target    :: desc_a
  integer(psb_ipk_), optional    :: data

  ! locals
  integer(psb_ipk_) :: ictxt, np, me, icomm, idxs, idxr, totxch, data_, err_act
  class(psb_i_base_vect_type), pointer :: d_vidx
  character(len=20)  :: name
  logical :: swap_persistent, swap_nonpersistent, do_alltoallv

  info=psb_success_
  name='psi_swap_datav'
  call psb_erractionsave(err_act)

  ictxt = desc_a%get_context()

  icomm = desc_a%get_mpic()
  ! TODO: get_mpic should be used to get dist_graph_comm, but for now this works
  swap_persistent = iand(flag,psb_swap_persistent_) /= 0
  swap_nonpersistent = iand(flag,psb_swap_nonpersistent_) /= 0
  do_alltoallv = swap_persistent .or. swap_nonpersistent
  if (do_alltoallv) then
    if (allocated(desc_a%dist_graph_comm)) then
      ! print *, "desc_a%dist_graph_comm", desc_a%dist_graph_comm
      icomm = desc_a%dist_graph_comm
    else
      print *, "ERROR: desc_a%dist_graph_comm not allocated but swap persistent flag passed"
      goto 9999
    end if
  end if

  call psb_info(ictxt,me,np)


  if (np == -1) then
    info=psb_err_context_error_
    call psb_errpush(info,name)
    goto 9999
  endif

  if (.not.psb_is_asb_desc(desc_a)) then
    info=psb_err_invalid_cd_state_
    call psb_errpush(info,name)
    goto 9999
  endif

  if(present(data)) then
    data_ = data
  else
    data_ = psb_comm_halo_
  end if

  call desc_a%get_list(data_,d_vidx,totxch,idxr,idxs,info)
  if (info /= psb_success_) then
    call psb_errpush(psb_err_internal_error_,name,a_err='psb_cd_get_list')
    goto 9999
  end if

  call psi_swapdata(ictxt,icomm,flag,beta,y,d_vidx,totxch,idxs,idxr,work,info)
  if (info /= psb_success_) goto 9999

  call psb_erractionrestore(err_act)
  return

9999 call psb_error_handler(ictxt,err_act)

    return
end subroutine psi_dswapdata_vect


!
!
! Subroutine: psi_dswap_vidx_vect
!   Data exchange among processes.
!
!   Takes care of Y an exanspulated vector. Relies on the gather/scatter methods
!   of vectors.
!
!   The real workhorse: the outer routine will only choose the index list
!   this one takes the index list and does the actual exchange.
!
!
!
subroutine psi_dswap_vidx_vect(iictxt,iicomm,flag,beta,y,idx, &
     & totxch,totsnd,totrcv,work,info)

  use psi_mod, psb_protect_name => psi_dswap_vidx_vect
  use psb_error_mod
  use psb_realloc_mod
  use psb_desc_mod
  use psb_penv_mod
  use psb_d_base_vect_mod
#ifdef MPI_MOD
  use mpi
#endif
  implicit none
#ifdef MPI_H
  include 'mpif.h'
#endif

  integer(psb_ipk_), intent(in)         :: iictxt,iicomm,flag
  integer(psb_ipk_), intent(out)        :: info
  class(psb_d_base_vect_type) :: y
  real(psb_dpk_)           :: beta
  real(psb_dpk_), target   :: work(:)
  class(psb_i_base_vect_type), intent(inout) :: idx
  integer(psb_ipk_), intent(in)              :: totxch,totsnd, totrcv

  ! locals
  integer(psb_mpk_) :: ictxt, icomm, np, me,&
       & proc_to_comm, p2ptag, p2pstat(mpi_status_size), iret
  integer(psb_mpk_), allocatable :: prcid(:)
  integer(psb_ipk_) :: nesd, nerv,&
       & err_act, i, j, ii, ri, si, idx_pt, totsnd_, totrcv_,&
       & snd_pt, rcv_pt, pnti, n, ierr
  logical :: swap_mpi, swap_sync, swap_send, swap_recv,&
      & albf,do_send,do_recv, &
      do_alltoallv, swap_persistent, swap_nonpersistent
  logical, parameter :: usersend=.false., debug=.false.
  character(len=20)  :: name

  !remove
  integer :: status(MPI_STATUS_SIZE)
  real(psb_dpk_) :: start_req_t, finish_req_t, start_comm_t, end_comm_t, &
      start_t, finish_t
  logical :: weight
  real(psb_dpk_)           :: tmp
  ! score-p declaration
  ! SCOREP_USER_REGION_DEFINE( MPIX_Neighbor_alltoallv_init_region )
  ! SCOREP_USER_REGION_DEFINE( artless_test )


  info=psb_success_
  name='psi_swap_datav'
  call psb_erractionsave(err_act)
  ictxt = iictxt
  icomm = iicomm

  call psb_info(ictxt,me,np)

  if (np == -1) then
    info=psb_err_context_error_
    call psb_errpush(info,name)
    goto 9999
  endif

  n=1
  swap_mpi        = iand(flag,psb_swap_mpi_) /= 0
  swap_sync       = iand(flag,psb_swap_sync_) /= 0
  swap_send       = iand(flag,psb_swap_send_) /= 0
  swap_recv       = iand(flag,psb_swap_recv_) /= 0
  swap_persistent     = iand(flag,psb_swap_persistent_) /= 0
  swap_nonpersistent  = iand(flag,psb_swap_nonpersistent_) /= 0
  do_send = swap_mpi .or. swap_sync .or. swap_send
  do_recv = swap_mpi .or. swap_sync .or. swap_recv
  do_alltoallv = swap_persistent .or. swap_nonpersistent

  totrcv_ = totrcv * n
  totsnd_ = totsnd * n
  call idx%sync()


  ! TODO check do_persistent twice, here and calling psi_dswapdata_vect,
  !      redudent but needed for now
  if (do_alltoallv) then
    if (.not. allocated(y%p)) then
      allocate(y%p)
    end if
    start_t = MPI_Wtime()
    if (.not. allocated(y%p%sndbuf)) then
      ! allocate time counters
      allocate(y%p%alltoall_comm_time, y%p%total_time)
      y%p%alltoall_comm_time = 0
      y%p%total_time = 0
      ! allocate buffers
      allocate(y%p%sndbuf(totsnd), y%p%rcvbuf(totrcv))
      allocate(y%p%rcv_count, y%p%snd_count)
      ! get number of neighbors and graph so we know who to communicate with
      call MPI_Dist_graph_neighbors_count(icomm, y%p%rcv_count, &
          y%p%snd_count, weight, ierr) ! should weight go into psb_d_persis_vect_type??
      allocate(y%p%rcv_from(y%p%rcv_count), y%p%rcv_ws(y%p%rcv_count))
      allocate(y%p%snd_to(y%p%snd_count), y%p%snd_ws(y%p%snd_count))
      call MPI_Dist_graph_neighbors(icomm, y%p%rcv_count, y%p%rcv_from, &
          y%p%rcv_ws, y%p%snd_count, y%p%snd_to, y%p%snd_ws, ierr)
      allocate(y%p%snd_counts(y%p%snd_count), y%p%rcv_counts(y%p%rcv_count), &
          y%p%snd_displs(y%p%snd_count), y%p%rcv_displs(y%p%rcv_count))

      y%p%snd_counts=0
      y%p%rcv_counts=0
      y%p%snd_displs=0
      y%p%rcv_displs=0

      si = 1 ! sndbuf index
      ri = 1
      pnti   = 1
      snd_pt = 1
      rcv_pt = 1
      do i=1, totxch
        proc_to_comm = idx%v(pnti+psb_proc_id_)
        nerv = idx%v(pnti+psb_n_elem_recv_)
        nesd = idx%v(pnti+nerv+psb_n_elem_send_)
        idx_pt = 1+pnti+psb_n_elem_recv_

        snd_pt = 1+pnti+nerv+psb_n_elem_send_
        rcv_pt = 1+pnti+psb_n_elem_recv_

        ! array of snd/rcv counts and displacements needed for request creation
        do ii=1, y%p%snd_count
          if (y%p%rcv_from(ii) == proc_to_comm) then
            y%p%rcv_counts(ii) = nerv
            y%p%rcv_displs(ii) = ri - 1
          end if
          if (y%p%snd_to(ii) == proc_to_comm) then
            y%p%snd_counts(ii) = nesd
            y%p%snd_displs(ii) = si - 1
          end if
        end do

        ! pack sndbuf
        do ii=0,nesd-1
          y%p%sndbuf(si) = y%v(idx%v(ii+snd_pt))
          si = si + 1
        end do
        ri = ri + nerv
        pnti = pnti + nerv + nesd + 3
      end do

! #ifdef DBG_PER
      ! print *, me, ": totxch", totxch,"total to send", totsnd, "total to recv", totrcv
      ! print *, me, ": rcv_count = ", y%p%rcv_count, "snd_count = ", y%p%snd_count
      ! print *, me, ":y%p%snd_displs = ",y%p%snd_displs, &
      !     "y%p%rcv_displs",y%p%rcv_displs
! #endif

      ! ----------------------------
      ! print *, me, ": totsnd=",totsnd ,"y%p%snd_to =", y%p%snd_to, &
      !     "y%p%snd_counts = ", y%p%snd_counts
      ! print *, me, ": totrcv=",totrcv,"y%p%rcv_from =", y%p%rcv_from, &
      !     "y%p%rcv_counts = ", y%p%rcv_counts
      ! ----------------------------


      ! ----------------------------
      ! Persistant Request Creation
      ! ----------------------------
      if (swap_persistent) then
        allocate(y%p%init_request)
        y%p%init_request = -1
        start_req_t = MPI_Wtime()
        call MPIX_Neighbor_alltoallv_init(y%p%sndbuf, y%p%snd_counts, y%p%snd_displs, &
            MPI_DOUBLE_PRECISION, y%p%rcvbuf, y%p%rcv_counts, y%p%rcv_displs, &
            MPI_DOUBLE_PRECISION, icomm, MPI_INFO_NULL, &
            y%p%init_request, ierr)
        finish_req_t = MPI_Wtime()
        allocate(y%p%request_create_time)
        y%p%request_create_time = finish_req_t - start_req_t
      end if
      ! ----------------------------
      ! if (ierr .ne. 0) then
      !   print *, "ERROR: MPIX_Neighbor_alltoallvinit ierr = ", ierr
      !   goto 9999
      ! end if
    else ! send and recv buffers exist, need to pack send buffer
      si = 1 ! sndbuf index
      pnti   = 1
      snd_pt = 1
      rcv_pt = 1
      do i=1, totxch
        proc_to_comm = idx%v(pnti+psb_proc_id_)
        nerv = idx%v(pnti+psb_n_elem_recv_)
        nesd = idx%v(pnti+nerv+psb_n_elem_send_)
        idx_pt = 1+pnti+psb_n_elem_recv_

        snd_pt = 1+pnti+nerv+psb_n_elem_send_
        rcv_pt = 1+pnti+psb_n_elem_recv_

        ! pack sndbuf
        do ii=0,nesd-1
          y%p%sndbuf(si) = y%v(idx%v(ii+snd_pt))
          si = si + 1
        end do
        pnti   = pnti + nerv + nesd + 3
      end do
    end if ! end packing phase

    ! ---- start communication ----
    start_comm_t = MPI_Wtime()
    ! ---- if persistant ----
    if (swap_persistent) then
      call MPI_Start(y%p%init_request, ierr)
      call MPI_Wait(y%p%init_request, status, ierr)
    end if

    ! ---- if non-persistent ----
    if (swap_nonpersistent) then
      call MPI_Neighbor_alltoallv(y%p%sndbuf, y%p%snd_counts, y%p%snd_displs, &
          MPI_DOUBLE_PRECISION, y%p%rcvbuf, y%p%rcv_counts, y%p%rcv_displs, &
          MPI_DOUBLE_PRECISION, icomm, ierr)
    end if
    end_comm_t = MPI_Wtime() - start_comm_t
    y%p%alltoall_comm_time = y%p%alltoall_comm_time + end_comm_t

    ! ----scatter/unpack buffer----
    pnti   = 1
    rcv_pt = 1
    do i=1, totxch
      ri     = 1
      proc_to_comm = idx%v(pnti+psb_proc_id_)
      nerv = idx%v(pnti+psb_n_elem_recv_)
      nesd = idx%v(pnti+nerv+psb_n_elem_send_)
      idx_pt = 1+pnti+psb_n_elem_recv_
      rcv_pt = 1+pnti+psb_n_elem_recv_

      do j=1, y%p%rcv_count
        ! make sure proc_to_comm is the rcv_from(j) so we get the right number
        ! of receive elements
        if (proc_to_comm .eq. y%p%rcv_from(j)) then
          if (nerv .ne. y%p%rcv_counts(j)) then
            print *, me, ":Error: psi_dswapdata.F90:", &
                "number of persistent elements received differs from nerv"
            goto 9999
          end if

          do ii=0,nerv-1
            y%v(idx%v(ii+rcv_pt)) = y%p%rcvbuf(ri + y%p%rcv_displs(j))
            ri = ri + 1
          end do
        end if
      end do
      pnti = pnti + nerv + nesd + 3
    end do

#ifdef DBG_PER
    print*, me, ": end scatter/unpack buffer"
#endif
    finish_t = MPI_Wtime()
    y%p%total_time = y%p%total_time + finish_t - start_t
    ! print *, me,":ARTLESS total_time = ", y%p%total_time
  end if ! ----- persistent communication is complete -----


  if (debug) write(*,*) me,'Internal buffer'
  if (do_send) then
    if (allocated(y%comid)) then
      if (any(y%comid /= mpi_request_null)) then
        !
        ! Unfinished communication? Something is wrong....
        !
        info=psb_err_mpi_error_
        call psb_errpush(info,name,m_err=(/-2/))
        goto 9999
      end if
    end if
    if (debug) write(*,*) me,'do_send start'
    call y%new_buffer(ione*size(idx%v),info)
    call y%new_comid(totxch,info)
    y%comid = mpi_request_null
    call psb_realloc(totxch,prcid,info)
    ! First I post all the non blocking receives
    pnti   = 1
    do i=1, totxch
      proc_to_comm = idx%v(pnti+psb_proc_id_)
      nerv = idx%v(pnti+psb_n_elem_recv_)
      nesd = idx%v(pnti+nerv+psb_n_elem_send_)

      rcv_pt = 1+pnti+psb_n_elem_recv_
      call psb_get_rank(prcid(i),ictxt,proc_to_comm)
      if ((nerv>0).and.(proc_to_comm /= me)) then
        if (debug) write(*,*) me,'Posting receive from',prcid(i),rcv_pt
        p2ptag = psb_double_swap_tag
        call mpi_irecv(y%combuf(rcv_pt),nerv,&
             & psb_mpi_r_dpk_,prcid(i),&
             & p2ptag, icomm,y%comid(i,2),iret)
      end if
      pnti   = pnti + nerv + nesd + 3
    end do
    if (debug) write(*,*) me,' Gather '
    !
    ! Then gather for sending.
    !
    pnti   = 1
    do i=1, totxch
      nerv = idx%v(pnti+psb_n_elem_recv_)
      nesd = idx%v(pnti+nerv+psb_n_elem_send_)
      snd_pt = 1+pnti+nerv+psb_n_elem_send_
      rcv_pt = 1+pnti+psb_n_elem_recv_
      idx_pt = snd_pt
      call y%gth(idx_pt,nesd,idx)
      pnti   = pnti + nerv + nesd + 3
    end do

    !
    ! Then wait
    !
    call y%device_wait()

    if (debug) write(*,*) me,' isend'
    !
    ! Then send
    !

    pnti   = 1
    snd_pt = 1
    rcv_pt = 1
    p2ptag = psb_double_swap_tag
    do i=1, totxch
      proc_to_comm = idx%v(pnti+psb_proc_id_)
      nerv = idx%v(pnti+psb_n_elem_recv_)
      nesd = idx%v(pnti+nerv+psb_n_elem_send_)
      snd_pt = 1+pnti+nerv+psb_n_elem_send_
      rcv_pt = 1+pnti+psb_n_elem_recv_

      if ((nesd>0).and.(proc_to_comm /= me)) then
        call mpi_isend(y%combuf(snd_pt),nesd,&
             & psb_mpi_r_dpk_,prcid(i),&
             & p2ptag,icomm,y%comid(i,1),iret)
      end if

      if(iret /= mpi_success) then
        info=psb_err_mpi_error_
        call psb_errpush(info,name,m_err=(/iret/))
        goto 9999
      end if

      pnti   = pnti + nerv + nesd + 3
    end do
  end if

  if (do_recv) then
    if (debug) write(*,*) me,' do_Recv'
    if (.not.allocated(y%comid)) then
      !
      ! No matching send? Something is wrong....
      !
      info=psb_err_mpi_error_
      call psb_errpush(info,name,m_err=(/-2/))
      goto 9999
    end if
    call psb_realloc(totxch,prcid,info)

    if (debug) write(*,*) me,' wait'
    pnti   = 1
    p2ptag = psb_double_swap_tag
    do i=1, totxch
      proc_to_comm = idx%v(pnti+psb_proc_id_)
      nerv = idx%v(pnti+psb_n_elem_recv_)
      nesd = idx%v(pnti+nerv+psb_n_elem_send_)
      snd_pt = 1+pnti+nerv+psb_n_elem_send_
      rcv_pt = 1+pnti+psb_n_elem_recv_

      if (proc_to_comm /= me)then
        if (nesd>0) then
          call mpi_wait(y%comid(i,1),p2pstat,iret)
          if(iret /= mpi_success) then
            info=psb_err_mpi_error_
            call psb_errpush(info,name,m_err=(/iret/))
            goto 9999
          end if
        end if
        if (nerv>0) then
          call mpi_wait(y%comid(i,2),p2pstat,iret)
          if(iret /= mpi_success) then
            info=psb_err_mpi_error_
            call psb_errpush(info,name,m_err=(/iret/))
            goto 9999
          end if
        end if
      else if (proc_to_comm == me) then
        if (nesd /= nerv) then
          write(psb_err_unit,*) &
               & 'Fatal error in swapdata: mismatch on self send',&
               & nerv,nesd
        end if
        y%combuf(rcv_pt:rcv_pt+nerv-1) = y%combuf(snd_pt:snd_pt+nesd-1)
      end if
      pnti   = pnti + nerv + nesd + 3
    end do

    if (debug) write(*,*) me,' scatter'
    pnti   = 1
    snd_pt = 1
    rcv_pt = 1
    do i=1, totxch
      proc_to_comm = idx%v(pnti+psb_proc_id_)
      nerv = idx%v(pnti+psb_n_elem_recv_)
      nesd = idx%v(pnti+nerv+psb_n_elem_send_)
      idx_pt = 1+pnti+psb_n_elem_recv_
      snd_pt = 1+pnti+nerv+psb_n_elem_send_
      rcv_pt = 1+pnti+psb_n_elem_recv_

      if (debug) write(0,*)me,' Received from: ',prcid(i),&
           & y%combuf(rcv_pt:rcv_pt+nerv-1)
      call y%sct(rcv_pt,nerv,idx,beta)
      pnti   = pnti + nerv + nesd + 3
    end do
    !
    ! Waited for everybody, clean up
    !
    y%comid = mpi_request_null

    !
    ! Then wait for device
    !
    if (debug) write(*,*) me,' wait'
    call y%device_wait()
    if (debug) write(*,*) me,' free buffer'
    call y%maybe_free_buffer(info)
    if (info == 0) call y%free_comid(info)
    if (info /= 0) then
      call psb_errpush(psb_err_alloc_dealloc_,name)
      goto 9999
    end if
    if (debug) write(*,*) me,' done'
  end if


  call psb_erractionrestore(err_act)
  return

9999 call psb_error_handler(iictxt,err_act)

  return
end subroutine psi_dswap_vidx_vect

!
!
! Subroutine: psi_dswapdata_multivect
!   Data exchange among processes.
!
!   Takes care of Y an encaspulated vector.
!
!
subroutine psi_dswapdata_multivect(flag,beta,y,desc_a,work,info,data)

  use psi_mod, psb_protect_name => psi_dswapdata_multivect
  use psb_d_base_multivect_mod
  use psb_error_mod
  use psb_desc_mod
  use psb_penv_mod
#ifdef MPI_MOD
  use mpi
#endif
  implicit none
#ifdef MPI_H
  include 'mpif.h'
#endif

  integer(psb_ipk_), intent(in)         :: flag
  integer(psb_ipk_), intent(out)        :: info
  class(psb_d_base_multivect_type) :: y
  real(psb_dpk_)           :: beta
  real(psb_dpk_), target   :: work(:)
  type(psb_desc_type), target  :: desc_a
  integer(psb_ipk_), optional           :: data

  ! locals
  integer(psb_ipk_) :: ictxt, np, me, icomm, idxs, idxr, totxch, data_, err_act
  class(psb_i_base_vect_type), pointer :: d_vidx
  character(len=20)  :: name

  info=psb_success_
  name='psi_swap_datav'
  call psb_erractionsave(err_act)

  ictxt = desc_a%get_context()
  icomm = desc_a%get_mpic()
  call psb_info(ictxt,me,np)
  print *, me, ": psi_dswapdata.F90:psi_dswapdata_multivect"
  if (np == -1) then
    info=psb_err_context_error_
    call psb_errpush(info,name)
    goto 9999
  endif

  if (.not.psb_is_asb_desc(desc_a)) then
    info=psb_err_invalid_cd_state_
    call psb_errpush(info,name)
    goto 9999
  endif

  if(present(data)) then
    data_ = data
  else
    data_ = psb_comm_halo_
  end if

  call desc_a%get_list(data_,d_vidx,totxch,idxr,idxs,info)
  if (info /= psb_success_) then
    call psb_errpush(psb_err_internal_error_,name,a_err='psb_cd_get_list')
    goto 9999
  end if

  call psi_swapdata(ictxt,icomm,flag,beta,y,d_vidx,totxch,idxs,idxr,work,info)
  if (info /= psb_success_) goto 9999

  call psb_erractionrestore(err_act)
  return

9999 call psb_error_handler(ictxt,err_act)

    return
end subroutine psi_dswapdata_multivect


!
!
! Subroutine: psi_dswap_vidx_multivect
!   Data exchange among processes.
!
!   Takes care of Y an encapsulated multivector. Relies on the gather/scatter methods
!   of multivectors.
!
!   The real workhorse: the outer routine will only choose the index list
!   this one takes the index list and does the actual exchange.
!
!
!
subroutine psi_dswap_vidx_multivect(iictxt,iicomm,flag,beta,y,idx, &
     & totxch,totsnd,totrcv,work,info)

  use psi_mod, psb_protect_name => psi_dswap_vidx_multivect
  use psb_error_mod
  use psb_realloc_mod
  use psb_desc_mod
  use psb_penv_mod
  use psb_d_base_multivect_mod
#ifdef MPI_MOD
  use mpi
#endif
  implicit none
#ifdef MPI_H
  include 'mpif.h'
#endif

  integer(psb_ipk_), intent(in)         :: iictxt,iicomm,flag
  integer(psb_ipk_), intent(out)        :: info
  class(psb_d_base_multivect_type) :: y
  real(psb_dpk_)         :: beta
  real(psb_dpk_), target :: work(:)
  class(psb_i_base_vect_type), intent(inout) :: idx
  integer(psb_ipk_), intent(in)              :: totxch,totsnd, totrcv

  ! locals
  integer(psb_mpk_) :: ictxt, icomm, np, me,&
       & proc_to_comm, p2ptag, p2pstat(mpi_status_size), iret
  integer(psb_mpk_), allocatable :: prcid(:)
  integer(psb_ipk_) :: nesd, nerv,&
       & err_act, i, idx_pt, totsnd_, totrcv_,&
       & snd_pt, rcv_pt, pnti, n
  logical :: swap_mpi, swap_sync, swap_send, swap_recv,&
       & albf,do_send,do_recv
  logical, parameter :: usersend=.false., debug=.false.
  character(len=20)  :: name

  info=psb_success_
  name='psi_swap_datav'
  call psb_erractionsave(err_act)
  ictxt = iictxt
  icomm = iicomm

  call psb_info(ictxt,me,np)
  print *, me, ": psi_dswapdata.F90:psi_dswap_vidx_multivect"
  if (np == -1) then
    info=psb_err_context_error_
    call psb_errpush(info,name)
    goto 9999
  endif

  n = y%get_ncols()

  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

  call idx%sync()

  if (debug) write(*,*) me,'Internal buffer'
  if (do_send) then
    if (allocated(y%comid)) then
      if (any(y%comid /= mpi_request_null)) then
        !
        ! Unfinished communication? Something is wrong....
        !
        info=psb_err_mpi_error_
        call psb_errpush(info,name,m_err=(/-2/))
        goto 9999
      end if
    end if
    if (debug) write(*,*) me,'do_send start'
    call y%new_buffer(ione*size(idx%v),info)
    call y%new_comid(totxch,info)
    y%comid = mpi_request_null
    call psb_realloc(totxch,prcid,info)
    ! First I post all the non blocking receives
    pnti   = 1
    snd_pt = totrcv_+1
    rcv_pt = 1
    do i=1, totxch
      proc_to_comm = idx%v(pnti+psb_proc_id_)
      nerv = idx%v(pnti+psb_n_elem_recv_)
      nesd = idx%v(pnti+nerv+psb_n_elem_send_)
      call psb_get_rank(prcid(i),ictxt,proc_to_comm)
      if ((nerv>0).and.(proc_to_comm /= me)) then
        if (debug) write(*,*) me,'Posting receive from',prcid(i),rcv_pt
        p2ptag = psb_double_swap_tag
        call mpi_irecv(y%combuf(rcv_pt),n*nerv,&
             & psb_mpi_r_dpk_,prcid(i),&
             & p2ptag, icomm,y%comid(i,2),iret)
      end if
      rcv_pt = rcv_pt + n*nerv
      snd_pt = snd_pt + n*nesd
      pnti   = pnti + nerv + nesd + 3
    end do
    if (debug) write(*,*) me,' Gather '
    !
    ! Then gather for sending.
    !
    pnti   = 1
    snd_pt = totrcv_+1
    rcv_pt = 1
    do i=1, totxch
      nerv = idx%v(pnti+psb_n_elem_recv_)
      nesd = idx%v(pnti+nerv+psb_n_elem_send_)
      idx_pt = 1+pnti+nerv+psb_n_elem_send_
      call y%gth(idx_pt,snd_pt,nesd,idx)
      rcv_pt = rcv_pt + n*nerv
      snd_pt = snd_pt + n*nesd
      pnti   = pnti + nerv + nesd + 3
    end do

    !
    ! Then wait for device
    !
    call y%device_wait()

    if (debug) write(*,*) me,' isend'
    !
    ! Then send
    !

    pnti   = 1
    snd_pt = totrcv_+1
    rcv_pt = 1
    p2ptag = psb_double_swap_tag
    do i=1, totxch
      proc_to_comm = idx%v(pnti+psb_proc_id_)
      nerv = idx%v(pnti+psb_n_elem_recv_)
      nesd = idx%v(pnti+nerv+psb_n_elem_send_)

      if ((nesd>0).and.(proc_to_comm /= me)) then
        call mpi_isend(y%combuf(snd_pt),n*nesd,&
             & psb_mpi_r_dpk_,prcid(i),&
             & p2ptag,icomm,y%comid(i,1),iret)
      end if

      if(iret /= mpi_success) then
        info=psb_err_mpi_error_
        call psb_errpush(info,name,m_err=(/iret/))
        goto 9999
      end if
      rcv_pt = rcv_pt + n*nerv
      snd_pt = snd_pt + n*nesd
      pnti   = pnti + nerv + nesd + 3
    end do
  end if

  if (do_recv) then
    if (debug) write(*,*) me,' do_Recv'
    if (.not.allocated(y%comid)) then
      !
      ! No matching send? Something is wrong....
      !
      info=psb_err_mpi_error_
      call psb_errpush(info,name,m_err=(/-2/))
      goto 9999
    end if
    call psb_realloc(totxch,prcid,info)

    if (debug) write(*,*) me,' wait'
    pnti   = 1
    snd_pt = totrcv_+1
    rcv_pt = 1
    p2ptag = psb_double_swap_tag
    do i=1, totxch
      proc_to_comm = idx%v(pnti+psb_proc_id_)
      nerv = idx%v(pnti+psb_n_elem_recv_)
      nesd = idx%v(pnti+nerv+psb_n_elem_send_)
      if (proc_to_comm /= me)then
        if (nesd>0) then
          call mpi_wait(y%comid(i,1),p2pstat,iret)
          if(iret /= mpi_success) then
            info=psb_err_mpi_error_
            call psb_errpush(info,name,m_err=(/iret/))
            goto 9999
          end if
        end if
        if (nerv>0) then
          call mpi_wait(y%comid(i,2),p2pstat,iret)
          if(iret /= mpi_success) then
            info=psb_err_mpi_error_
            call psb_errpush(info,name,m_err=(/iret/))
            goto 9999
          end if
        end if
      else if (proc_to_comm == me) then
        if (nesd /= nerv) then
          write(psb_err_unit,*) &
               & 'Fatal error in swapdata: mismatch on self send',&
               & nerv,nesd
        end if
        y%combuf(rcv_pt:rcv_pt+n*nerv-1) = y%combuf(snd_pt:snd_pt+n*nesd-1)
      end if
      rcv_pt = rcv_pt + n*nerv
      snd_pt = snd_pt + n*nesd
      pnti   = pnti + nerv + nesd + 3
    end do

    if (debug) write(*,*) me,' scatter'
    pnti   = 1
    snd_pt = totrcv_+1
    rcv_pt = 1
    do i=1, totxch
      proc_to_comm = idx%v(pnti+psb_proc_id_)
      nerv = idx%v(pnti+psb_n_elem_recv_)
      nesd = idx%v(pnti+nerv+psb_n_elem_send_)
      idx_pt = 1+pnti+psb_n_elem_recv_

      if (debug) write(0,*)me,' Received from: ',prcid(i),&
           & y%combuf(rcv_pt:rcv_pt+n*nerv-1)
      call y%sct(idx_pt,rcv_pt,nerv,idx,beta)
      rcv_pt = rcv_pt + n*nerv
      snd_pt = snd_pt + n*nesd
      pnti   = pnti + nerv + nesd + 3
    end do
    !
    ! Waited for com, cleanup comid
    !
    y%comid = mpi_request_null

    !
    ! Then wait for device
    !
    if (debug) write(*,*) me,' wait'
    call y%device_wait()
    if (debug) write(*,*) me,' free buffer'
    call y%free_buffer(info)
    if (info == 0) call y%free_comid(info)
    if (info /= 0) then
      call psb_errpush(psb_err_alloc_dealloc_,name)
      goto 9999
    end if
    if (debug) write(*,*) me,' done'
  end if


  call psb_erractionrestore(err_act)
  return

9999 call psb_error_handler(iictxt,err_act)

  return
end subroutine psi_dswap_vidx_multivect