!!$ !!$ Parallel Sparse BLAS version 2.2 !!$ (C) Copyright 2006/2007/2008 !!$ 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. !!$ !!$ ! ! File: psi_fnd_owner.f90 ! ! Subroutine: psi_fnd_owner ! Figure out who owns global indices. ! ! Arguments: ! nv - integer Number of indices required on the calling ! process ! idx(:) - integer Required indices on the calling process. ! Note: the indices should be unique! ! iprc(:) - integer, allocatable Output: process identifiers for the corresponding ! indices ! desc_a - type(psb_desc_type). The communication descriptor. ! info - integer. return code. ! subroutine psi_fnd_owner(nv,idx,iprc,desc,info) use psb_descriptor_type use psb_serial_mod use psb_const_mod use psb_error_mod use psb_penv_mod use psb_realloc_mod use psi_mod, psb_protect_name => psi_fnd_owner #ifdef MPI_MOD use mpi #endif implicit none #ifdef MPI_H include 'mpif.h' #endif integer, intent(in) :: nv integer, intent(in) :: idx(:) integer, allocatable, intent(out) :: iprc(:) type(psb_desc_type), intent(in) :: desc integer, intent(out) :: info integer, allocatable :: hsz(:),hidx(:),helem(:),hproc(:),& & sdsz(:),sdidx(:), rvsz(:), rvidx(:),answers(:,:),idxsrch(:,:) integer :: i,n_row,n_col,err_act,ih,icomm,hsize,ip,isz,k,j,& & last_ih, last_j integer :: ictxt,np,me logical, parameter :: gettime=.false. real(psb_dpk_) :: t0, t1, t2, t3, t4, tamx, tidx character(len=20) :: name info = psb_success_ name = 'psi_fnd_owner' call psb_erractionsave(err_act) ictxt = psb_cd_get_context(desc) icomm = psb_cd_get_mpic(desc) n_row = psb_cd_get_local_rows(desc) n_col = psb_cd_get_local_cols(desc) ! check on blacs grid call psb_info(ictxt, me, np) if (np == -1) then info = psb_err_context_error_ call psb_errpush(info,name) goto 9999 endif if (nv < 0 ) then info = psb_err_context_error_ call psb_errpush(info,name) goto 9999 endif if (.not.(psb_is_ok_desc(desc))) then call psb_errpush(psb_err_from_subroutine_,name,a_err='invalid desc') goto 9999 end if if (gettime) then t0 = psb_wtime() end if ! ! The basic idea is very simple. ! First we collect (to all) all the requests. Allocate(hidx(np+1),hsz(np),& & sdsz(0:np-1),sdidx(0:np-1),& & rvsz(0:np-1),rvidx(0:np-1),& & stat=info) if (info /= psb_success_) then call psb_errpush(psb_err_from_subroutine_,name,a_err='Allocate') goto 9999 end if hsz = 0 hsz(me+1) = nv call psb_amx(ictxt,hsz) hidx(1) = 0 do i=1, np hidx(i+1) = hidx(i) + hsz(i) end do hsize = hidx(np+1) Allocate(helem(hsize),hproc(hsize),stat=info) if (info /= psb_success_) then call psb_errpush(psb_err_from_subroutine_,name,a_err='Allocate') goto 9999 end if if (gettime) then t3 = psb_wtime() end if call mpi_allgatherv(idx,hsz(me+1),mpi_integer,& & hproc,hsz,hidx,mpi_integer,& & icomm,info) if (gettime) then tamx = psb_wtime() - t3 end if ! Second, we figure out locally whether we own the indices (whoever is ! asking for them). if (gettime) then t3 = psb_wtime() end if call psi_idx_cnv(hsize,hproc,helem,desc,info,owned=.true.) if (gettime) then tidx = psb_wtime()-t3 end if if (info == psb_err_iarray_outside_bounds_) info = psb_success_ if (info /= psb_success_) then call psb_errpush(psb_err_from_subroutine_,name,a_err='psi_idx_cnv') goto 9999 end if ! Third: we build the answers for those indices we own, ! with a section for each process asking. hidx = hidx +1 j = 0 do ip = 0, np-1 sdidx(ip) = j sdsz(ip) = 0 do i=hidx(ip+1), hidx(ip+1+1)-1 if ((0 < helem(i)).and. (helem(i) <= n_row)) then j = j + 1 hproc(j) = hproc(i) sdsz(ip) = sdsz(ip) + 1 end if end do end do if (gettime) then t3 = psb_wtime() end if ! Collect all the answers with alltoallv (need sizes) call mpi_alltoall(sdsz,1,mpi_integer,rvsz,1,mpi_integer,icomm,info) isz = sum(rvsz) allocate(answers(isz,2),idxsrch(nv,2),stat=info) if (info /= psb_success_) then call psb_errpush(psb_err_from_subroutine_,name,a_err='Allocate') goto 9999 end if j = 0 do ip=0, np-1 rvidx(ip) = j j = j + rvsz(ip) end do call mpi_alltoallv(hproc,sdsz,sdidx,mpi_integer,& & answers(:,1),rvsz,rvidx,mpi_integer,& & icomm,info) if (gettime) then tamx = psb_wtime() - t3 + tamx end if j = 1 do ip = 0,np-1 do k=1,rvsz(ip) answers(j,2) = ip j = j + 1 end do end do ! Sort the answers and the requests, so we can ! match them efficiently call psb_msort(answers(:,1),ix=answers(:,2),& & flag=psb_sort_keep_idx_) idxsrch(1:nv,1) = idx(1:nv) call psb_msort(idxsrch(1:nv,1),ix=idxsrch(1:nv,2)) ! Now extract the answers for our local query call psb_realloc(nv,iprc,info) if (info /= psb_success_) then call psb_errpush(psb_err_from_subroutine_,name,a_err='psb_realloc') goto 9999 end if last_ih = -1 last_j = -1 j = 1 do i=1, nv ih = idxsrch(i,1) if (ih == last_ih) then iprc(idxsrch(i,2)) = answers(last_j,2) else do if (j > size(answers,1)) then ! Last resort attempt. j = psb_ibsrch(ih,size(answers,1),answers(:,1)) if (j == -1) then write(0,*) me,'psi_fnd_owner: searching for ',ih, & & 'not found : ',size(answers,1),':',answers(:,1) info = psb_err_internal_error_ call psb_errpush(psb_err_internal_error_,name,a_err='out bounds srch ih') goto 9999 end if end if if (answers(j,1) == ih) exit if (answers(j,1) > ih) then k = j j = psb_ibsrch(ih,k,answers(1:k,1)) if (j == -1) then write(0,*) me,'psi_fnd_owner: searching for ',ih, & & 'not found : ',size(answers,1),':',answers(:,1) info = psb_err_internal_error_ call psb_errpush(psb_err_internal_error_,name,a_err='out bounds srch ih') goto 9999 end if end if j = j + 1 end do ! Note that the answers here are given in order ! of sending process, so we are implicitly getting ! the max process index in case of overlap. last_ih = ih do last_j = j iprc(idxsrch(i,2)) = answers(j,2) j = j + 1 if (j > size(answers,1)) exit if (answers(j,1) /= ih) exit end do end if end do if (gettime) then call psb_barrier(ictxt) t1 = psb_wtime() t1 = t1 -t0 - tamx - tidx call psb_amx(ictxt,tamx) call psb_amx(ictxt,tidx) call psb_amx(ictxt,t1) if (me == psb_root_) then write(*,'(" fnd_owner idx time : ",es10.4)') tidx write(*,'(" fnd_owner amx time : ",es10.4)') tamx write(*,'(" fnd_owner remainedr : ",es10.4)') t1 endif end if 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(ictxt) end if return end subroutine psi_fnd_owner