psblas3/base/internals/psi_desc_index.F90

!   
!                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_desc_index.f90
!
! Subroutine: psi_desc_index
!    Converts a list of data exchanges from build format to assembled format. 
!    See below for a description of the formats.
!
! Arguments:
! desc_a       - type(psb_desc_type) The descriptor; in this context only the index 
!                                     mapping parts are used.
! index_in(:)  - integer             The index list, build format  
! index_out(:) - integer(psb_ipk_), allocatable  The index list, assembled format
! nxch         - integer             The number of data exchanges on the calling process
! nsnd         - integer             Total send buffer size       on the calling process
! nrcv         - integer             Total receive buffer size    on the calling process
!
! The format of the index lists. Copied from base/modules/psb_desc_type
!
!  7. The data exchange is based on lists of local indices to be exchanged; all the 
!     lists have the same format, as follows:
!     the list is  composed of variable dimension blocks, one for each process to 
!     communicate with; each block contains indices of local elements to be 
!     exchanged. We do choose the order of communications:  do not change 
!     the sequence of blocks unless you know what you're doing, or you'll 
!     risk a deadlock. NOTE: This is the format when the state is PSB_ASB_.
!     See below for BLD. The end-of-list is marked with a -1. 
!
!  notation        stored in		          explanation
!  --------------- --------------------------- -----------------------------------
!  process_id      index_v(p+proc_id_)      identifier of process with which 
!                                                data is  exchanged.
!  n_elements_recv index_v(p+n_elem_recv_)  number of elements to receive.
!  elements_recv   index_v(p+elem_recv_+i)  indexes of local elements to
!					          receive. these are stored in the
!					          array from location p+elem_recv_ to
!					          location p+elem_recv_+
!						  index_v(p+n_elem_recv_)-1.
!  n_elements_send index_v(p+n_elem_send_)  number of elements to send.
!  elements_send   index_v(p+elem_send_+i)  indexes of local elements to
!					          send. these are stored in the
!					          array from location p+elem_send_ to
!					          location p+elem_send_+
!						  index_v(p+n_elem_send_)-1.
!
!     This organization is valid for both halo and overlap indices; overlap entries
!     need to be updated to ensure that a variable at a given global index 
!     (assigned to multiple processes) has the same value. The way to resolve the 
!     issue is to exchange the data and then sum (or average) the values. See
!     psb_ovrl subroutine. 
!  
!  8. When the descriptor is in the BLD state the INDEX vectors contains only 
!     the indices to be received, organized as  a sequence 
!     of entries of the form (proc,N,(lx1,lx2,...,lxn)) with owning process,
!     number of indices (most often N=1), list of local indices. 
!     This is because we only know the list of halo indices to be received 
!     as we go about building the sparse matrix pattern, and we want the build 
!     phase to be loosely synchronized. Thus we record the indices we have to ask 
!     for, and at the time we call PSB_CDASB we match all the requests to figure 
!     out who should be sending what to whom.
!     However this implies that we know who owns the indices; if we are in the 
!     LARGE case (as described above) this is actually only true for the OVERLAP list 
!     that is filled in at CDALL time, and not for the HALO; thus the HALO list 
!     is rebuilt during the CDASB process (in the psi_ldsc_pre_halo subroutine). 
!
!
subroutine psi_i_desc_index(desc,index_in,dep_list,&
     & length_dl,nsnd,nrcv,desc_index,info)
  use psb_desc_mod
  use psb_realloc_mod
  use psb_error_mod
  use psb_const_mod
#ifdef MPI_MOD
  use mpi
#endif
  use psb_penv_mod
  use psb_timers_mod
  use psi_mod, psb_protect_name => psi_i_desc_index
  implicit none
#ifdef MPI_H
  include 'mpif.h'
#endif

  !    ...array parameters.....
  type(psb_desc_type) :: desc
  integer(psb_ipk_) :: index_in(:)
  integer(psb_ipk_) :: dep_list(:)
  integer(psb_ipk_),allocatable  :: desc_index(:)
  integer(psb_ipk_) :: length_dl,nsnd,nrcv,info
  !    ....local scalars...        
  integer(psb_ipk_) :: j,me,np,i,proc
  !    ...parameters...
  type(psb_ctxt_type) :: ctxt
  integer(psb_ipk_), parameter  :: no_comm=-1
  !     ...local arrays..
  integer(psb_lpk_),allocatable  :: sndbuf(:), rcvbuf(:)

  integer(psb_mpk_),allocatable  :: brvindx(:),rvsz(:),&
       & bsdindx(:),sdsz(:)

  integer(psb_mpk_) :: proc_to_comm, p2ptag, p2pstat(mpi_status_size),&
       & iret, sz
  integer(psb_mpk_), allocatable :: prcid(:), rvhd(:)

  integer(psb_ipk_) :: ihinsz,ntot,k,err_act,nidx,&
       & idxr, idxs, iszs, iszr, nesd, nerv, ixp, idx
  integer(psb_mpk_) :: icomm, minfo

  logical, parameter :: do_timings=.false., oldstyle=.false., debug=.false. 
  integer(psb_ipk_), save :: idx_phase1=-1, idx_phase2=-1, idx_phase3=-1, idx_phase4=-1
  logical, parameter :: usempi=.false.
  integer(psb_ipk_) :: debug_level, debug_unit
  character(len=20) :: name

  info = psb_success_
  name='psi_desc_index'
  call psb_erractionsave(err_act)
  debug_unit  = psb_get_debug_unit()
  debug_level = psb_get_debug_level() 

  ctxt = desc%get_context()
  icomm = desc%get_mpic()
  call psb_info(ctxt,me,np) 
  if (np == -1) then
    info = psb_err_context_error_
    call psb_errpush(info,name)
    goto 9999
  endif

  if (debug_level >= psb_debug_inner_) then 
    write(debug_unit,*) me,' ',trim(name),': start'
    call psb_barrier(ctxt)
  endif
  if ((do_timings).and.(idx_phase1==-1)) &
       & idx_phase1 = psb_get_timer_idx("I_DSC_IDX: phase1   ")
  if ((do_timings).and.(idx_phase2==-1)) &
       & idx_phase2 = psb_get_timer_idx("I_DSC_IDX: phase2   ")
  if ((do_timings).and.(idx_phase3==-1)) &
       & idx_phase3 = psb_get_timer_idx("I_DSC_IDX: phase3   ")
  if ((do_timings).and.(idx_phase4==-1)) &
       & idx_phase4 = psb_get_timer_idx("I_DSC_IDX: phase4   ")

  if (do_timings) call psb_tic(idx_phase1)
  ! 
  !     first, find out the sizes to be exchanged.
  !     note: things marked here as sndbuf/rcvbuf (for mpi) corresponds to things  
  !     to be received/sent (in the final psblas descriptor).
  !     be careful of the inversion
  !
  allocate(sdsz(np),rvsz(np),bsdindx(np),brvindx(np),stat=info)
  if(info /= psb_success_) then
    info=psb_err_alloc_dealloc_
    call psb_errpush(info,name)
    goto 9999
  end if

  sdsz(:)    = 0
  rvsz(:)    = 0
  bsdindx(:) = 0
  brvindx(:) = 0
  i = 1
  do 
    if (index_in(i) == -1) exit
    proc = index_in(i)
    i = i + 1 
    nerv = index_in(i)
    sdsz(proc+1) = sdsz(proc+1) + nerv
    i = i + nerv + 1 
  end do
  ihinsz=i
  call mpi_alltoall(sdsz,1,psb_mpi_mpk_,rvsz,1,psb_mpi_mpk_,icomm,minfo)
  if (minfo /= psb_success_) then
    call psb_errpush(psb_err_from_subroutine_,name,a_err='mpi_alltoall')
    goto 9999
  end if

  i    = 1
  idxs = 0
  idxr = 0
  do i=1, length_dl
    proc = dep_list(i)
    bsdindx(proc+1) = idxs
    idxs = idxs + sdsz(proc+1)
    brvindx(proc+1) = idxr
    idxr = idxr + rvsz(proc+1)
  end do
  iszs = sum(sdsz)
  iszr = sum(rvsz)
  nsnd = iszr
  nrcv = iszs  

  if ((iszs /= idxs).or.(iszr /= idxr)) then 
    write(psb_err_unit,*) me, trim(name),': Warning: strange results?', &
         & iszs,idxs,iszr,idxr
  end if
  if (debug_level >= psb_debug_inner_) then 
    write(debug_unit,*) me,' ',trim(name),': computed sizes ',iszr,iszs
    call psb_barrier(ctxt)
  endif

  ntot = (3*(count((sdsz>0).or.(rvsz>0)))+ iszs + iszr) + 1

  call psb_ensure_size(ntot,desc_index,info)

  if (info /= psb_success_) then 
    call psb_errpush(psb_err_from_subroutine_,name,a_err='psb_realloc')
    goto 9999
  end if

  if (debug_level >= psb_debug_inner_) then 
    write(debug_unit,*) me,' ',trim(name),': computed allocated workspace ',iszr,iszs
    call psb_barrier(ctxt)
  endif
  allocate(sndbuf(iszs),rcvbuf(iszr),stat=info)
  if(info /= psb_success_) then
    info=psb_err_alloc_dealloc_
    call psb_errpush(info,name)
    goto 9999
  end if
  if (do_timings) call psb_toc(idx_phase1)
  if (do_timings) call psb_tic(idx_phase2)

  ! 
  !     Second build the lists of requests
  !   
  i = 1
  do 
    if (i > ihinsz) then 
!!$      write(psb_err_unit,*) me,' did not find index_in end??? ',i,ihinsz
      exit
    end if
    if (index_in(i) == -1) exit
    proc = index_in(i)
    i = i + 1 
    nerv = index_in(i)
    !  
    !   note that here bsdinx is zero-based, hence the following loop
    !          
    call desc%indxmap%l2g(index_in(i+1:i+nerv),&
         & sndbuf(bsdindx(proc+1)+1:bsdindx(proc+1)+nerv),&
         & info) 

    if (info /= psb_success_) then
      call psb_errpush(psb_err_from_subroutine_,name,a_err='l2g')
      goto 9999
    end if

    bsdindx(proc+1) = bsdindx(proc+1) + nerv
    i = i + nerv + 1 
  end do

  if (debug_level >= psb_debug_inner_) then 
    write(debug_unit,*) me,' ',trim(name),': prepared send buffer '
    call psb_barrier(ctxt)
  endif
  !
  !   now have to regenerate bsdindx
  !  
  idxs = 0
  idxr = 0
  do i=1, length_dl
    proc = dep_list(i)
    bsdindx(proc+1) = idxs
    idxs = idxs + sdsz(proc+1)
    brvindx(proc+1) = idxr
    idxr = idxr + rvsz(proc+1)
  end do
  if (do_timings) call psb_toc(idx_phase2)
  if (do_timings) call psb_tic(idx_phase3)

  if (usempi) then   
    call mpi_alltoallv(sndbuf,sdsz,bsdindx,psb_mpi_lpk_,&
         & rcvbuf,rvsz,brvindx,psb_mpi_lpk_,icomm,minfo)
    if (minfo /= psb_success_) then
      call psb_errpush(psb_err_from_subroutine_,name,a_err='mpi_alltoallv')
      goto 9999
    end if
  else
    if (.true.) then 
      allocate(prcid(length_dl),rvhd(length_dl))
      prcid = -1
      ixp   = 1 
      do i=1, length_dl
        proc = dep_list(i)
        prcid(ixp) = psb_get_mpi_rank(ctxt,proc)
        sz         = rvsz(proc+1)
        if (sz > 0) then
          p2ptag   = psb_long_tag
          idx      = brvindx(proc+1)
          call mpi_irecv(rcvbuf(idx+1:idx+sz),sz,&
               & psb_mpi_lpk_, prcid(ixp), p2ptag, icomm,&
               & rvhd(ixp),iret)
        end if
        ixp = ixp + 1
      end do
      ixp   = 1 
      do i=1, length_dl
        proc       = dep_list(i)
        prcid(ixp) = psb_get_mpi_rank(ctxt,proc)
        sz         = sdsz(proc+1)
        if (sz > 0) then
          p2ptag   = psb_long_tag
          idx      = bsdindx(proc+1)
          call mpi_send(sndbuf(idx+1:idx+sz),sz,&
               & psb_mpi_lpk_, prcid(ixp), p2ptag, &
               & icomm,iret)
        end if
        ixp = ixp + 1
      end do
      ixp   = 1 
      do i=1, length_dl
        proc = dep_list(i)
        prcid(ixp) = psb_get_mpi_rank(ctxt,proc)
        sz         = rvsz(proc+1)
        if (sz > 0) then
          call mpi_wait(rvhd(ixp),p2pstat,iret)
        end if
        ixp = ixp + 1
      end do
    else

      do i=1, length_dl
        proc       = dep_list(i)
        sz         = sdsz(proc+1)
        idx      = bsdindx(proc+1)
        if (sz > 0) then
          call psb_snd(ctxt,sndbuf(idx+1:idx+sz), proc)
        end if
      end do

      do i=1, length_dl
        proc = dep_list(i)
        sz         = rvsz(proc+1)
        idx      = brvindx(proc+1)
        if (sz > 0) then
          call psb_rcv(ctxt,rcvbuf(idx+1:idx+sz),proc)
        end if
      end do

    end if
  end if
  if (do_timings) call psb_toc(idx_phase3)
  if (do_timings) call psb_tic(idx_phase4)

  !
  !  at this point we can finally build the output desc_index. beware
  !  of snd/rcv inversion. 
  !      
  i = 1
  do k = 1, length_dl
    proc = dep_list(k)
    desc_index(i) = proc
    i = i + 1 
    nerv = sdsz(proc+1) 
    desc_index(i) = nerv
    call desc%indxmap%g2l(sndbuf(bsdindx(proc+1)+1:bsdindx(proc+1)+nerv),&
         &  desc_index(i+1:i+nerv),info)
    i = i + nerv + 1 
    nesd = rvsz(proc+1) 
    desc_index(i) = nesd
    call desc%indxmap%g2l(rcvbuf(brvindx(proc+1)+1:brvindx(proc+1)+nesd),&
         &  desc_index(i+1:i+nesd),info)
    i = i + nesd + 1 
  end do
  desc_index(i) = - 1 

  deallocate(sdsz,rvsz,bsdindx,brvindx,sndbuf,rcvbuf,stat=info)
  if (info /= psb_success_) then 
    info=psb_err_alloc_dealloc_
    call psb_errpush(info,name)
    goto 9999
  end if
  if (do_timings) call psb_toc(idx_phase4)
  if (debug_level >= psb_debug_inner_) then 
    write(debug_unit,*) me,' ',trim(name),': done'
    call psb_barrier(ctxt)
  endif

  call psb_erractionrestore(err_act)
  return

9999 call psb_error_handler(ctxt,err_act)

  return

end subroutine psi_i_desc_index