diff --git a/Makefile b/Makefile index 829d4fc1..41caa3cd 100644 --- a/Makefile +++ b/Makefile @@ -15,15 +15,15 @@ libd: (if test ! -d include ; then mkdir include; fi; $(INSTALL_DATA) Make.inc include/Make.inc.psblas) (if test ! -d modules ; then mkdir modules; fi;) based: - cd base && $(MAKE) lib + $(MAKE) -C base lib precd: - cd prec && $(MAKE) lib + $(MAKE) -C prec lib kryld: - cd krylov && $(MAKE) lib + $(MAKE) -C krylov lib utild: - cd util&& $(MAKE) lib + $(MAKE) -C util lib cbindd: - cd cbind&& $(MAKE) lib + $(MAKE) -C cbind lib install: all mkdir -p $(INSTALL_INCLUDEDIR) &&\ @@ -42,11 +42,11 @@ install: all /bin/cp -fr test/pargen test/fileread test/kernel $(INSTALL_SAMPLESDIR) && \ mkdir -p $(INSTALL_SAMPLESDIR)/cbind && /bin/cp -fr cbind/test/pargen/* $(INSTALL_SAMPLESDIR)/cbind clean: - cd base && $(MAKE) clean - cd prec && $(MAKE) clean - cd krylov && $(MAKE) clean - cd util && $(MAKE) clean - cd cbind && $(MAKE) clean + $(MAKE) -C base clean + $(MAKE) -C prec clean + $(MAKE) -C krylov clean + $(MAKE) -C util clean + $(MAKE) -C cbind clean check: all make check -C test/serial diff --git a/base/Makefile b/base/Makefile index 3304176c..039a5963 100644 --- a/base/Makefile +++ b/base/Makefile @@ -13,25 +13,25 @@ lib: mods sr cm in pb tl sr cm in pb tl: mods mods: - cd modules && $(MAKE) lib LIBNAME=$(BASELIBNAME) F90="$(MPF90)" F90COPT="$(F90COPT) $(MPI_OPT)" + $(MAKE) -C modules lib LIBNAME=$(BASELIBNAME) F90="$(MPF90)" F90COPT="$(F90COPT) $(MPI_OPT)" sr: - cd serial && $(MAKE) lib LIBNAME=$(BASELIBNAME) + $(MAKE) -C serial lib LIBNAME=$(BASELIBNAME) cm: - cd comm && $(MAKE) lib LIBNAME=$(BASELIBNAME) + $(MAKE) -C comm lib LIBNAME=$(BASELIBNAME) in: - cd internals && $(MAKE) lib LIBNAME=$(BASELIBNAME) + $(MAKE) -C internals lib LIBNAME=$(BASELIBNAME) pb: - cd psblas && $(MAKE) lib LIBNAME=$(BASELIBNAME) + $(MAKE) -C psblas lib LIBNAME=$(BASELIBNAME) tl: - cd tools && $(MAKE) lib LIBNAME=$(BASELIBNAME) + $(MAKE) -C tools lib LIBNAME=$(BASELIBNAME) clean: - (cd modules; $(MAKE) clean) - (cd comm; $(MAKE) clean) - (cd internals; $(MAKE) clean) - (cd tools; $(MAKE) clean) - (cd serial; $(MAKE) clean) - (cd psblas; $(MAKE) clean) + ($(MAKE) -C modules clean) + ($(MAKE) -C comm clean) + ($(MAKE) -C internals clean) + ($(MAKE) -C tools clean) + ($(MAKE) -C serial clean) + ($(MAKE) -C psblas clean) veryclean: clean /bin/rm -f $(HERE)/$(LIBNAME) $(LIBMOD) *$(.mod) diff --git a/base/comm/Makefile b/base/comm/Makefile index f27d38cf..950a95a0 100644 --- a/base/comm/Makefile +++ b/base/comm/Makefile @@ -32,7 +32,7 @@ lib: interns mpfobjs $(OBJS) $(RANLIB) $(LIBDIR)/$(LIBNAME) interns: - cd internals && $(MAKE) lib + $(MAKE) -C internals lib mpfobjs: $(MAKE) $(MPFOBJS) FC="$(MPFC)" diff --git a/base/internals/psi_graph_fnd_owner.F90 b/base/internals/psi_graph_fnd_owner.F90 index 2ecdb24c..b951945e 100644 --- a/base/internals/psi_graph_fnd_owner.F90 +++ b/base/internals/psi_graph_fnd_owner.F90 @@ -162,7 +162,7 @@ subroutine psi_graph_fnd_owner(idx,iprc,idxmap,info) ! call psb_safe_ab_cpy(idxmap%p_adjcncy,ladj,info) nadj = psb_size(ladj) - ! This makes ladj allocated with size 0 just in case + ! This makes ladj allocated with size 0 if needed, as opposed to unallocated call psb_realloc(nadj,ladj,info) ! ! Throughout the subroutine, nreqst is the number of local inquiries diff --git a/base/internals/psi_indx_map_fnd_owner.F90 b/base/internals/psi_indx_map_fnd_owner.F90 index 0aee0806..25f900f0 100644 --- a/base/internals/psi_indx_map_fnd_owner.F90 +++ b/base/internals/psi_indx_map_fnd_owner.F90 @@ -74,7 +74,8 @@ subroutine psi_indx_map_fnd_owner(idx,iprc,idxmap,info) integer(psb_ipk_), allocatable :: hhidx(:) integer(psb_mpk_) :: icomm, minfo, iictxt - integer(psb_ipk_) :: i, err_act, hsize, nv + integer(psb_ipk_) :: i, err_act, hsize + integer(psb_lpk_) :: nv integer(psb_lpk_) :: mglob integer(psb_ipk_) :: ictxt,np,me, nresp logical, parameter :: gettime=.false. @@ -140,7 +141,66 @@ subroutine psi_indx_map_fnd_owner(idx,iprc,idxmap,info) else - call psi_graph_fnd_owner(idx,iprc,idxmap,info) + if (allocated(idxmap%halo_owner)) then + ! + ! Maybe we are coming here after a REINIT event. + ! In this case, reuse the existing information as much as possible. + ! + block + integer(psb_ipk_), allocatable :: tprc(:), lidx(:) + integer(psb_lpk_), allocatable :: tidx(:) + integer(psb_lpk_) :: k1, k2, nh + allocate(lidx(nv),stat=info) + if (info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,a_err='Allocate') + goto 9999 + end if + ! + ! Get local answers, if any + ! + call idxmap%g2l(idx,lidx,info,owned=.false.) + call idxmap%qry_halo_owner(lidx,iprc,info) + + nh = count(iprc<0) + !write(0,*) me,'Going through new impl from ',nv,' to ',nh + allocate(tidx(nh),tprc(nh),stat=info) + if (info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,a_err='Allocate') + goto 9999 + end if + ! + ! Prepare remote queries + ! + k2 = 0 + do k1 = 1, nv + if (iprc(k1) < 0) then + k2 = k2 + 1 + if (k2 > nh) then + info = psb_err_internal_error_ + call psb_errpush(info,name,a_err='Wrong auxiliary count') + goto 9999 + end if + tidx(k2) = idx(k1) + end if + end do + call psi_graph_fnd_owner(tidx,tprc,idxmap,info) + k2 = 0 + do k1 = 1, nv + if (iprc(k1) < 0) then + k2 = k2 + 1 + if (k2 > nh) then + info = psb_err_internal_error_ + call psb_errpush(info,name,a_err='Wrong auxiliary count') + goto 9999 + end if + iprc(k1) = tprc(k2) + end if + end do + end block + else + call psi_graph_fnd_owner(idx,iprc,idxmap,info) + end if + end if diff --git a/base/modules/desc/psb_hash_map_mod.f90 b/base/modules/desc/psb_hash_map_mod.f90 index 2f4e8653..ae109d5a 100644 --- a/base/modules/desc/psb_hash_map_mod.f90 +++ b/base/modules/desc/psb_hash_map_mod.f90 @@ -94,8 +94,8 @@ module psb_hash_map_mod procedure, pass(idxmap) :: lg2lv1_ins => hash_g2lv1_ins procedure, pass(idxmap) :: lg2lv2_ins => hash_g2lv2_ins - procedure, pass(idxmap) :: hash_cpy - generic, public :: assignment(=) => hash_cpy +!!$ procedure, pass(idxmap) :: hash_cpy +!!$ generic, public :: assignment(=) => hash_cpy procedure, pass(idxmap) :: bld_g2l_map => hash_bld_g2l_map end type psb_hash_map @@ -1464,7 +1464,7 @@ contains end if select type (outmap) - type is (psb_hash_map) + type is (psb_hash_map) call idxmap%psb_indx_map%cpy(outmap%psb_indx_map,info) if (info == psb_success_) then outmap%hashvsize = idxmap%hashvsize @@ -1478,7 +1478,7 @@ contains & call psb_safe_ab_cpy(idxmap%glb_lc,outmap%glb_lc,info) if (info == psb_success_)& & call psb_hash_copy(idxmap%hash,outmap%hash,info) -!!$ outmap = idxmap + class default ! This should be impossible info = -1 @@ -1499,28 +1499,30 @@ contains end subroutine hash_clone - subroutine hash_cpy(outmap,idxmap) - use psb_penv_mod - use psb_error_mod - use psb_realloc_mod - implicit none - class(psb_hash_map), intent(in) :: idxmap - type(psb_hash_map), intent(out) :: outmap - integer(psb_ipk_) :: info - - info = psb_success_ - outmap%psb_indx_map = idxmap%psb_indx_map - outmap%hashvsize = idxmap%hashvsize - outmap%hashvmask = idxmap%hashvmask - if (info == psb_success_)& - & call psb_safe_ab_cpy(idxmap%loc_to_glob,outmap%loc_to_glob,info) - if (info == psb_success_)& - & call psb_safe_ab_cpy(idxmap%hashv,outmap%hashv,info) - if (info == psb_success_)& - & call psb_safe_ab_cpy(idxmap%glb_lc,outmap%glb_lc,info) - if (info == psb_success_)& - & call psb_hash_copy(idxmap%hash,outmap%hash,info) - end subroutine hash_cpy +!!$ subroutine hash_cpy(outmap,idxmap) +!!$ use psb_penv_mod +!!$ use psb_error_mod +!!$ use psb_realloc_mod +!!$ implicit none +!!$ class(psb_hash_map), intent(in) :: idxmap +!!$ type(psb_hash_map), intent(out) :: outmap +!!$ integer(psb_ipk_) :: info +!!$ +!!$ info = psb_success_ +!!$ call idxmap%psb_indx_map%cpy(outmap%psb_indx_map,info) +!!$ if (info == psb_success_) then +!!$ outmap%hashvsize = idxmap%hashvsize +!!$ outmap%hashvmask = idxmap%hashvmask +!!$ end if +!!$ if (info == psb_success_)& +!!$ & call psb_safe_ab_cpy(idxmap%loc_to_glob,outmap%loc_to_glob,info) +!!$ if (info == psb_success_)& +!!$ & call psb_safe_ab_cpy(idxmap%hashv,outmap%hashv,info) +!!$ if (info == psb_success_)& +!!$ & call psb_safe_ab_cpy(idxmap%glb_lc,outmap%glb_lc,info) +!!$ if (info == psb_success_)& +!!$ & call psb_hash_copy(idxmap%hash,outmap%hash,info) +!!$ end subroutine hash_cpy subroutine hash_reinit(idxmap,info) use psb_penv_mod @@ -1533,7 +1535,7 @@ contains integer(psb_lpk_) :: lk integer(psb_lpk_) :: ntot integer(psb_ipk_) :: ictxt, me, np - integer(psb_ipk_), allocatable :: lidx(:) + integer(psb_ipk_), allocatable :: lidx(:), tadj(:), th_own(:) integer(psb_lpk_), allocatable :: gidx(:) character(len=20) :: name='hash_reinit' logical, parameter :: debug=.false. @@ -1549,13 +1551,16 @@ contains lidx = (/(k,k=1,nc)/) gidx = (/(lk,lk=1,nc)/) call idxmap%l2gip(gidx,info) - + tadj = idxmap%get_p_adjcncy() + call idxmap%get_halo_owner(th_own,info) + call idxmap%free() call hash_init_vlu(idxmap,ictxt,ntot,nr,gidx(1:nr),info) if (nc>nr) then call idxmap%g2lip_ins(gidx(nr+1:nc),info,lidx=lidx(nr+1:nc)) end if - + call idxmap%set_p_adjcncy(tadj) + call idxmap%set_halo_owner(th_own,info) if (info /= psb_success_) then info = psb_err_from_subroutine_ diff --git a/base/modules/desc/psb_indx_map_mod.f90 b/base/modules/desc/psb_indx_map_mod.f90 index f5aeaa4d..dd72280b 100644 --- a/base/modules/desc/psb_indx_map_mod.f90 +++ b/base/modules/desc/psb_indx_map_mod.f90 @@ -219,9 +219,9 @@ module psb_indx_map_mod procedure, pass(idxmap) :: set_halo_owner => base_set_halo_owner procedure, pass(idxmap) :: get_halo_owner => base_get_halo_owner - procedure, pass(idxmap) :: fnd_halo_owner_s => base_fnd_halo_owner_s - procedure, pass(idxmap) :: fnd_halo_owner_v => base_fnd_halo_owner_v - generic, public :: fnd_halo_owner => fnd_halo_owner_s, fnd_halo_owner_v + procedure, pass(idxmap) :: qry_halo_owner_s => base_qry_halo_owner_s + procedure, pass(idxmap) :: qry_halo_owner_v => base_qry_halo_owner_v + generic, public :: qry_halo_owner => qry_halo_owner_s, qry_halo_owner_v procedure, pass(idxmap) :: fnd_owner => psi_indx_map_fnd_owner procedure, pass(idxmap) :: init_vl => base_init_vl @@ -245,7 +245,7 @@ module psb_indx_map_mod & base_lg2lv2_ins, base_init_vl, base_is_null,& & base_row_extendable, base_clone, base_cpy, base_reinit, & & base_set_halo_owner, base_get_halo_owner, & - & base_fnd_halo_owner_s, base_fnd_halo_owner_v,& + & base_qry_halo_owner_s, base_qry_halo_owner_v,& & base_get_p_adjcncy, base_set_p_adjcncy, base_xtnd_p_adjcncy !> Function: psi_indx_map_fnd_owner @@ -1486,6 +1486,7 @@ contains end subroutine base_set_halo_owner subroutine base_get_halo_owner(idxmap,v,info) + use psb_realloc_mod use psb_penv_mod use psb_error_mod implicit none @@ -1494,13 +1495,15 @@ contains integer(psb_ipk_), intent(out) :: info integer(psb_ipk_) :: nh - nh = min(size(v),size(idxmap%halo_owner)) - v(1:nh) = idxmap%halo_owner(1:nh) + nh = size(idxmap%halo_owner) + !v = idxmap%halo_owner(1:nh) + call psb_safe_ab_cpy(idxmap%halo_owner,v,info) end subroutine base_get_halo_owner - subroutine base_fnd_halo_owner_s(idxmap,xin,xout,info) + subroutine base_qry_halo_owner_s(idxmap,xin,xout,info) use psb_penv_mod use psb_error_mod + use psb_realloc_mod implicit none class(psb_indx_map), intent(inout) :: idxmap integer(psb_ipk_), intent(in) :: xin @@ -1510,24 +1513,26 @@ contains integer(psb_ipk_) :: i, j, nr, nc, nh nr = idxmap%local_rows nc = idxmap%local_cols + nc = min(idxmap%local_cols, (nr+psb_size(idxmap%halo_owner))) xout = -1 if (.not.allocated(idxmap%halo_owner)) then !write(0,*) 'Halo_owner not allocated!', nr, nc, xin return end if if ((nr hlstart) then - call p_desc_c%indxmap%fnd_halo_owner(j,proc,info) + call p_desc_c%indxmap%qry_halo_owner(j,proc,info) sdsz(proc+1) = sdsz(proc+1) +1 end if end do @@ -264,19 +263,19 @@ subroutine psb_lc_coo_glob_transpose(ain,desc_r,info,atrans,desc_c,desc_rx) acoo%ja(nzd) = acoo%ja(k) acoo%val(nzd) = acoo%val(k) else - call p_desc_c%indxmap%fnd_halo_owner(j,proc,info) + call p_desc_c%indxmap%qry_halo_owner(j,proc,info) tsdx(proc+1) = tsdx(proc+1) +1 iasnd(tsdx(proc+1)) = acoo%ia(k) jasnd(tsdx(proc+1)) = acoo%ja(k) valsnd(tsdx(proc+1)) = acoo%val(k) end if end do + call acoo%set_nzeros(nzd) ! ! Put halo entries in global numbering ! call desc_r%indxmap%l2gip(jasnd(1:iszs),info) call p_desc_c%indxmap%l2gip(iasnd(1:iszs),info) - call acoo%set_nzeros(nzd) ! And exchange data. ! Normally we'll use our SIMPLE A2AV and not MPI, because ! the communication pattern is sparse, so ours is more @@ -407,19 +406,293 @@ subroutine psb_c_coo_glob_transpose(ain,desc_r,info,atrans,desc_c,desc_rx) type(psb_desc_type), intent(out), optional :: desc_rx integer(psb_ipk_), intent(out) :: info - type(psb_lc_coo_sparse_mat) :: atcoo + integer(psb_ipk_) :: ictxt, np,me + integer(psb_ipk_) :: counter,proc, err_act, j + integer(psb_ipk_) :: i, k, idx, r, ipx,mat_recv, iszs, iszr,idxs,idxr,nz,& + & l1, nsnds, nrcvs, nr,nc,nzl, hlstart, nzd + integer(psb_lpk_) :: nzt, lszr + integer(psb_mpk_) :: icomm, minfo + integer(psb_mpk_), allocatable :: brvindx(:), & + & rvsz(:), bsdindx(:), sdsz(:), tsdx(:), trvx(:) + integer(psb_ipk_), allocatable :: halo_owner(:) + integer(psb_lpk_), allocatable :: iasnd(:), jasnd(:),iarcv(:),jarcv(:) + complex(psb_spk_), allocatable :: valsnd(:) + type(psb_c_coo_sparse_mat), allocatable :: acoo + logical :: rowcnv_,colcnv_,rowscale_,colscale_,outcol_glob_ + type(psb_desc_type), pointer :: p_desc_c + character(len=5) :: outfmt_ + integer(psb_ipk_) :: debug_level, debug_unit + character(len=20) :: name, ch_err - if (present(atrans)) then - call ain%cp_to_lcoo(atcoo,info) + if(psb_get_errstatus() /= 0) return + info=psb_success_ + name='mld_glob_transpose' + call psb_erractionsave(err_act) + if (psb_errstatus_fatal()) then + info = psb_err_internal_error_ ; goto 9999 + end if + debug_unit = psb_get_debug_unit() + debug_level = psb_get_debug_level() + + ictxt = desc_r%get_context() + icomm = desc_r%get_mpic() + + Call psb_info(ictxt, me, np) + + if (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),': Start' + + if (present(desc_c)) then + p_desc_c => desc_c + else + p_desc_c => desc_r + end if + + Allocate(brvindx(np+1),& + & rvsz(np),sdsz(np),bsdindx(np+1), acoo,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 + l1 = 0 + brvindx(1) = 0 + bsdindx(1) = 0 + counter=1 + idx = 0 + idxs = 0 + idxr = 0 + + if (present(atrans)) then + call ain%cp_to_coo(acoo,info) + else + call ain%mv_to_coo(acoo,info) + end if + + + ! + ! Compute number of entries in the + ! halo part, sorted by destination process + ! + nr = desc_r%get_local_rows() + nc = p_desc_c%get_local_cols() + nzl = acoo%get_nzeros() + hlstart = p_desc_c%get_local_rows() + do k=1, nzl + j = acoo%ja(k) + if (j > hlstart) then + call p_desc_c%indxmap%qry_halo_owner(j,proc,info) + sdsz(proc+1) = sdsz(proc+1) +1 + end if + end do + + ! + ! Exchange sizes, so as to know sends/receives. + ! This is different from the halo exchange because the + ! number of entries was not precomputed in the descriptor, + ! which was vector-oriented and not matrix-entry-oriented + ! + call mpi_alltoall(sdsz,1,psb_mpi_mpk_,& + & rvsz,1,psb_mpi_mpk_,icomm,minfo) + + if (info /= psb_success_) then + info=psb_err_from_subroutine_ + call psb_errpush(info,name,a_err='mpi_alltoall') + goto 9999 + end if + nsnds = count(sdsz /= 0) + nrcvs = count(rvsz /= 0) + if (debug_level >= psb_debug_outer_)& + & write(debug_unit,*) me,' ',trim(name),': Done initial alltoall',nsnds,nrcvs + + idxs = 0 + idxr = 0 + counter = 1 + Do proc = 0, np-1 + bsdindx(proc+1) = idxs + idxs = idxs + sdsz(proc+1) + brvindx(proc+1) = idxr + idxr = idxr + rvsz(proc+1) + Enddo + + tsdx = bsdindx + trvx = brvindx + + iszr = sum(rvsz) + iszs = sum(sdsz) + + if (debug_level >= psb_debug_outer_)& + & write(debug_unit,*) me,' ',trim(name),': Sizes:',& + & ' Send:',sdsz(:),' Receive:',rvsz(:) + + call psb_ensure_size(max(iszs,1),iasnd,info) + if (info == psb_success_) call psb_ensure_size(max(iszs,1),jasnd,info) + if (info == psb_success_) call psb_ensure_size(max(iszs,1),valsnd,info) + if (info == psb_success_) call psb_ensure_size(max(iszr,1),iarcv,info) + if (info == psb_success_) call psb_ensure_size(max(iszr,1),jarcv,info) + if (info /= psb_success_) then + info=psb_err_from_subroutine_ + call psb_errpush(info,name,a_err='ensure_size') + goto 9999 + end if + + ! + ! Now, transpose the matrix, then split between itself + ! and the send buffers + ! + call acoo%transp() + if (acoo%get_nzeros()/= nzl) then + write(0,*) me,'Something strange upon transpose: ',nzl,acoo%get_nzeros() + end if + nzl = acoo%get_nzeros() + hlstart = p_desc_c%get_local_rows() + + + nzd = 0 + do k = 1, nzl + j = acoo%ia(k) + if (j<=hlstart) then + nzd = nzd + 1 + acoo%ia(nzd) = acoo%ia(k) + acoo%ja(nzd) = acoo%ja(k) + acoo%val(nzd) = acoo%val(k) + else + call p_desc_c%indxmap%qry_halo_owner(j,proc,info) + tsdx(proc+1) = tsdx(proc+1) +1 + iasnd(tsdx(proc+1)) = acoo%ia(k) + jasnd(tsdx(proc+1)) = acoo%ja(k) + valsnd(tsdx(proc+1)) = acoo%val(k) + end if + end do + call acoo%set_nzeros(nzd) + ! + ! Put halo entries in global numbering + ! + call desc_r%indxmap%l2gip(jasnd(1:iszs),info) + call p_desc_c%indxmap%l2gip(iasnd(1:iszs),info) + ! And exchange data. + ! Normally we'll use our SIMPLE A2AV and not MPI, because + ! the communication pattern is sparse, so ours is more + ! efficient. Using ACOO for the receive buffers. + nzl = acoo%get_nzeros() + call acoo%ensure_size(nzl+iszr) + + select case(psb_get_sp_a2av_alg()) + case(psb_sp_a2av_smpl_triad_) + call psb_simple_triad_a2av(valsnd,iasnd,jasnd,sdsz,bsdindx,& + & acoo%val(nzl+1:nzl+iszr),iarcv(1:iszr),& + & jarcv(1:iszr),rvsz,brvindx,ictxt,info) + case(psb_sp_a2av_smpl_v_) + call psb_simple_a2av(valsnd,sdsz,bsdindx,& + & acoo%val(nzl+1:nzl+iszr),rvsz,brvindx,ictxt,info) + if (info == psb_success_) call psb_simple_a2av(iasnd,sdsz,bsdindx,& + & iarcv(1:iszr),rvsz,brvindx,ictxt,info) + if (info == psb_success_) call psb_simple_a2av(jasnd,sdsz,bsdindx,& + & jarcv(1:iszr),rvsz,brvindx,ictxt,info) + case(psb_sp_a2av_mpi_) + + call mpi_alltoallv(valsnd,sdsz,bsdindx,psb_mpi_c_spk_,& + & acoo%val(nzl+1:nzl+iszr),rvsz,brvindx,psb_mpi_c_spk_,icomm,minfo) + if (minfo == mpi_success) & + & call mpi_alltoallv(iasnd,sdsz,bsdindx,psb_mpi_lpk_,& + & iarcv(1:iszr),rvsz,brvindx,psb_mpi_lpk_,icomm,minfo) + if (minfo == mpi_success) & + & call mpi_alltoallv(jasnd,sdsz,bsdindx,psb_mpi_lpk_,& + & jarcv(1:iszr),rvsz,brvindx,psb_mpi_lpk_,icomm,minfo) + if (minfo /= mpi_success) info = minfo + case default + info = psb_err_internal_error_ + call psb_errpush(info,name,a_err='wrong A2AV alg selector') + goto 9999 + end select + + if (info /= psb_success_) then + info=psb_err_from_subroutine_ + call psb_errpush(info,name,a_err='alltoallv') + goto 9999 + end if + + if (debug_level >= psb_debug_outer_)& + & write(debug_unit,*) me,' ',trim(name),': Done alltoallv' + + if (present(desc_rx)) then + ! + ! Extend the appropriate descriptor; started as R but on + ! transpose it now describes C + ! + call desc_r%clone(desc_rx,info) + call psb_cd_reinit(desc_rx,info) + ! + ! Take out non-local rows + ! + call psb_glob_to_loc(iarcv(1:iszr),p_desc_c,info,iact='I',owned=.true.) + lszr = iszr + call psb_coo_clean_negidx_inner(lszr,iarcv(1:iszr),jarcv(1:iszr),& + & acoo%val(nzl+1:nzl+iszr),nzt,info) + call desc_rx%g2lip_ins(jarcv(1:nzt),info) + call psb_cdasb(desc_rx,info) + acoo%ia(nzl+1:nzl+nzt) = iarcv(1:nzt) + acoo%ja(nzl+1:nzl+nzt) = jarcv(1:nzt) + nzl = nzl + nzt + call acoo%set_nzeros(nzl) + nzl = acoo%get_nzeros() + call acoo%set_sorted(.false.) + ! + ! Insert to extend descriptor + ! + call acoo%set_nrows(p_desc_c%get_local_rows()) + call acoo%set_ncols(desc_rx%get_local_cols()) + !write(0,*) me,' Trans RX ',acoo%get_nrows(),acoo%get_ncols(),acoo%get_nzeros() else - call ain%mv_to_lcoo(atcoo,info) + ! + ! + ! Take out non-local rows + ! + call psb_glob_to_loc(iarcv(1:iszr),p_desc_c,info,iact='I',owned=.true.) + call psb_glob_to_loc(jarcv(1:iszr),desc_r,info,iact='I') + lszr = iszr + call psb_coo_clean_negidx_inner(lszr,iarcv(1:iszr),jarcv(1:iszr),& + & acoo%val(nzl+1:nzl+iszr),nzt,info) + acoo%ia(nzl+1:nzl+nzt) = iarcv(1:nzt) + acoo%ja(nzl+1:nzl+nzt) = jarcv(1:nzt) + nzl = nzl + nzt + call acoo%set_nzeros(nzl) + nzl = acoo%get_nzeros() + call acoo%set_sorted(.false.) + + call acoo%set_nrows(p_desc_c%get_local_rows()) + call acoo%set_ncols(desc_r%get_local_cols()) + !write(0,*) me,' Trans R- ',acoo%get_nrows(),acoo%get_ncols(),acoo%get_nzeros() end if - if (info == 0) call psb_glob_transpose(atcoo,desc_r,info,desc_c=desc_c,desc_rx=desc_rx) + +!!$ write(0,*) me,' Sanity check after rx%g2l :',count(acoo%ja(1:nzl)<0) + + + call acoo%fix(info) + nzl = acoo%get_nzeros() + if (present(atrans)) then - call atrans%mv_from_lcoo(atcoo,info) + call atrans%mv_from_coo(acoo,info) else - call ain%mv_from_lcoo(atcoo,info) + call ain%mv_from_coo(acoo,info) end if + Deallocate(brvindx,bsdindx,rvsz,sdsz,& + & iasnd,jasnd,valsnd,& + & iarcv,jarcv,stat=info) + if (debug_level >= psb_debug_outer_)& + & write(debug_unit,*) me,' ',trim(name),': Done' + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(ictxt,err_act) + + return end subroutine psb_c_coo_glob_transpose subroutine psb_c_simple_glob_transpose_ip(ain,desc_a,info) @@ -434,7 +707,7 @@ subroutine psb_c_simple_glob_transpose_ip(ain,desc_a,info) ! that the same DESC_A works for both A and A^T, which ! essentially means that A has a symmetric pattern. ! - type(psb_lc_coo_sparse_mat) :: tmpc1, tmpc2 + type(psb_c_coo_sparse_mat) :: tmpc1, tmpc2 integer(psb_ipk_) :: nz1, nz2, nzh, nz integer(psb_ipk_) :: ictxt, me, np integer(psb_lpk_) :: i, j, k, nrow, ncol, nlz @@ -485,7 +758,7 @@ subroutine psb_c_simple_glob_transpose(ain,aout,desc_a,info) ! that the same DESC_A works for both A and A^T, which ! essentially means that A has a symmetric pattern. ! - type(psb_lc_coo_sparse_mat) :: tmpc1, tmpc2 + type(psb_c_coo_sparse_mat) :: tmpc1, tmpc2 integer(psb_ipk_) :: nz1, nz2, nzh, nz integer(psb_ipk_) :: ictxt, me, np integer(psb_lpk_) :: i, j, k, nrow, ncol, nlz diff --git a/base/tools/psb_cd_reinit.f90 b/base/tools/psb_cd_reinit.f90 index d579ba95..15dbbc59 100644 --- a/base/tools/psb_cd_reinit.f90 +++ b/base/tools/psb_cd_reinit.f90 @@ -69,6 +69,9 @@ Subroutine psb_cd_reinit(desc,info) call psb_move_alloc(tmp_halo,desc%halo_index,info) call psb_move_alloc(tmp_ext,desc%ext_index,info) call desc%indxmap%reinit(info) +!!$ if (me == 0) write(0,*) 'On cdreinit status :',& +!!$ & allocated(desc%indxmap%p_adjcncy),allocated(desc%indxmap%halo_owner), & +!!$ & desc%get_fmt() ! call psb_cd_set_bld(desc,info) end if diff --git a/base/tools/psb_cspasb.f90 b/base/tools/psb_cspasb.f90 index ae5c0af2..073fcbbd 100644 --- a/base/tools/psb_cspasb.f90 +++ b/base/tools/psb_cspasb.f90 @@ -122,11 +122,11 @@ subroutine psb_cspasb(a,desc_a, info, afmt, upd, dupl, mold) end if - IF (debug_level >= psb_debug_ext_) then + if (debug_level >= psb_debug_ext_) then ch_err=a%get_fmt() write(debug_unit, *) me,' ',trim(name),': From SPCNV',& & info,' ',ch_err - end IF + end if if (psb_errstatus_fatal()) then info=psb_err_from_subroutine_ diff --git a/base/tools/psb_d_glob_transpose.F90 b/base/tools/psb_d_glob_transpose.F90 index 697e7ad2..638ba174 100644 --- a/base/tools/psb_d_glob_transpose.F90 +++ b/base/tools/psb_d_glob_transpose.F90 @@ -113,7 +113,6 @@ subroutine psb_ld_coo_glob_transpose(ain,desc_r,info,atrans,desc_c,desc_rx) integer(psb_ipk_) :: ictxt, np,me integer(psb_ipk_) :: counter,proc, err_act, j integer(psb_lpk_) :: i, k, idx, r, ipx,mat_recv, iszs, iszr,idxs,idxr,nz,& - & irmin,icmin,irmax,icmax,& & l1, nsnds, nrcvs, nr,nc,nzl, hlstart, nzt, nzd integer(psb_mpk_) :: icomm, minfo integer(psb_mpk_), allocatable :: brvindx(:), & @@ -190,7 +189,7 @@ subroutine psb_ld_coo_glob_transpose(ain,desc_r,info,atrans,desc_c,desc_rx) do k=1, nzl j = acoo%ja(k) if (j > hlstart) then - call p_desc_c%indxmap%fnd_halo_owner(j,proc,info) + call p_desc_c%indxmap%qry_halo_owner(j,proc,info) sdsz(proc+1) = sdsz(proc+1) +1 end if end do @@ -264,19 +263,19 @@ subroutine psb_ld_coo_glob_transpose(ain,desc_r,info,atrans,desc_c,desc_rx) acoo%ja(nzd) = acoo%ja(k) acoo%val(nzd) = acoo%val(k) else - call p_desc_c%indxmap%fnd_halo_owner(j,proc,info) + call p_desc_c%indxmap%qry_halo_owner(j,proc,info) tsdx(proc+1) = tsdx(proc+1) +1 iasnd(tsdx(proc+1)) = acoo%ia(k) jasnd(tsdx(proc+1)) = acoo%ja(k) valsnd(tsdx(proc+1)) = acoo%val(k) end if end do + call acoo%set_nzeros(nzd) ! ! Put halo entries in global numbering ! call desc_r%indxmap%l2gip(jasnd(1:iszs),info) call p_desc_c%indxmap%l2gip(iasnd(1:iszs),info) - call acoo%set_nzeros(nzd) ! And exchange data. ! Normally we'll use our SIMPLE A2AV and not MPI, because ! the communication pattern is sparse, so ours is more @@ -407,19 +406,293 @@ subroutine psb_d_coo_glob_transpose(ain,desc_r,info,atrans,desc_c,desc_rx) type(psb_desc_type), intent(out), optional :: desc_rx integer(psb_ipk_), intent(out) :: info - type(psb_ld_coo_sparse_mat) :: atcoo + integer(psb_ipk_) :: ictxt, np,me + integer(psb_ipk_) :: counter,proc, err_act, j + integer(psb_ipk_) :: i, k, idx, r, ipx,mat_recv, iszs, iszr,idxs,idxr,nz,& + & l1, nsnds, nrcvs, nr,nc,nzl, hlstart, nzd + integer(psb_lpk_) :: nzt, lszr + integer(psb_mpk_) :: icomm, minfo + integer(psb_mpk_), allocatable :: brvindx(:), & + & rvsz(:), bsdindx(:), sdsz(:), tsdx(:), trvx(:) + integer(psb_ipk_), allocatable :: halo_owner(:) + integer(psb_lpk_), allocatable :: iasnd(:), jasnd(:),iarcv(:),jarcv(:) + real(psb_dpk_), allocatable :: valsnd(:) + type(psb_d_coo_sparse_mat), allocatable :: acoo + logical :: rowcnv_,colcnv_,rowscale_,colscale_,outcol_glob_ + type(psb_desc_type), pointer :: p_desc_c + character(len=5) :: outfmt_ + integer(psb_ipk_) :: debug_level, debug_unit + character(len=20) :: name, ch_err - if (present(atrans)) then - call ain%cp_to_lcoo(atcoo,info) + if(psb_get_errstatus() /= 0) return + info=psb_success_ + name='mld_glob_transpose' + call psb_erractionsave(err_act) + if (psb_errstatus_fatal()) then + info = psb_err_internal_error_ ; goto 9999 + end if + debug_unit = psb_get_debug_unit() + debug_level = psb_get_debug_level() + + ictxt = desc_r%get_context() + icomm = desc_r%get_mpic() + + Call psb_info(ictxt, me, np) + + if (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),': Start' + + if (present(desc_c)) then + p_desc_c => desc_c + else + p_desc_c => desc_r + end if + + Allocate(brvindx(np+1),& + & rvsz(np),sdsz(np),bsdindx(np+1), acoo,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 + l1 = 0 + brvindx(1) = 0 + bsdindx(1) = 0 + counter=1 + idx = 0 + idxs = 0 + idxr = 0 + + if (present(atrans)) then + call ain%cp_to_coo(acoo,info) + else + call ain%mv_to_coo(acoo,info) + end if + + + ! + ! Compute number of entries in the + ! halo part, sorted by destination process + ! + nr = desc_r%get_local_rows() + nc = p_desc_c%get_local_cols() + nzl = acoo%get_nzeros() + hlstart = p_desc_c%get_local_rows() + do k=1, nzl + j = acoo%ja(k) + if (j > hlstart) then + call p_desc_c%indxmap%qry_halo_owner(j,proc,info) + sdsz(proc+1) = sdsz(proc+1) +1 + end if + end do + + ! + ! Exchange sizes, so as to know sends/receives. + ! This is different from the halo exchange because the + ! number of entries was not precomputed in the descriptor, + ! which was vector-oriented and not matrix-entry-oriented + ! + call mpi_alltoall(sdsz,1,psb_mpi_mpk_,& + & rvsz,1,psb_mpi_mpk_,icomm,minfo) + + if (info /= psb_success_) then + info=psb_err_from_subroutine_ + call psb_errpush(info,name,a_err='mpi_alltoall') + goto 9999 + end if + nsnds = count(sdsz /= 0) + nrcvs = count(rvsz /= 0) + if (debug_level >= psb_debug_outer_)& + & write(debug_unit,*) me,' ',trim(name),': Done initial alltoall',nsnds,nrcvs + + idxs = 0 + idxr = 0 + counter = 1 + Do proc = 0, np-1 + bsdindx(proc+1) = idxs + idxs = idxs + sdsz(proc+1) + brvindx(proc+1) = idxr + idxr = idxr + rvsz(proc+1) + Enddo + + tsdx = bsdindx + trvx = brvindx + + iszr = sum(rvsz) + iszs = sum(sdsz) + + if (debug_level >= psb_debug_outer_)& + & write(debug_unit,*) me,' ',trim(name),': Sizes:',& + & ' Send:',sdsz(:),' Receive:',rvsz(:) + + call psb_ensure_size(max(iszs,1),iasnd,info) + if (info == psb_success_) call psb_ensure_size(max(iszs,1),jasnd,info) + if (info == psb_success_) call psb_ensure_size(max(iszs,1),valsnd,info) + if (info == psb_success_) call psb_ensure_size(max(iszr,1),iarcv,info) + if (info == psb_success_) call psb_ensure_size(max(iszr,1),jarcv,info) + if (info /= psb_success_) then + info=psb_err_from_subroutine_ + call psb_errpush(info,name,a_err='ensure_size') + goto 9999 + end if + + ! + ! Now, transpose the matrix, then split between itself + ! and the send buffers + ! + call acoo%transp() + if (acoo%get_nzeros()/= nzl) then + write(0,*) me,'Something strange upon transpose: ',nzl,acoo%get_nzeros() + end if + nzl = acoo%get_nzeros() + hlstart = p_desc_c%get_local_rows() + + + nzd = 0 + do k = 1, nzl + j = acoo%ia(k) + if (j<=hlstart) then + nzd = nzd + 1 + acoo%ia(nzd) = acoo%ia(k) + acoo%ja(nzd) = acoo%ja(k) + acoo%val(nzd) = acoo%val(k) + else + call p_desc_c%indxmap%qry_halo_owner(j,proc,info) + tsdx(proc+1) = tsdx(proc+1) +1 + iasnd(tsdx(proc+1)) = acoo%ia(k) + jasnd(tsdx(proc+1)) = acoo%ja(k) + valsnd(tsdx(proc+1)) = acoo%val(k) + end if + end do + call acoo%set_nzeros(nzd) + ! + ! Put halo entries in global numbering + ! + call desc_r%indxmap%l2gip(jasnd(1:iszs),info) + call p_desc_c%indxmap%l2gip(iasnd(1:iszs),info) + ! And exchange data. + ! Normally we'll use our SIMPLE A2AV and not MPI, because + ! the communication pattern is sparse, so ours is more + ! efficient. Using ACOO for the receive buffers. + nzl = acoo%get_nzeros() + call acoo%ensure_size(nzl+iszr) + + select case(psb_get_sp_a2av_alg()) + case(psb_sp_a2av_smpl_triad_) + call psb_simple_triad_a2av(valsnd,iasnd,jasnd,sdsz,bsdindx,& + & acoo%val(nzl+1:nzl+iszr),iarcv(1:iszr),& + & jarcv(1:iszr),rvsz,brvindx,ictxt,info) + case(psb_sp_a2av_smpl_v_) + call psb_simple_a2av(valsnd,sdsz,bsdindx,& + & acoo%val(nzl+1:nzl+iszr),rvsz,brvindx,ictxt,info) + if (info == psb_success_) call psb_simple_a2av(iasnd,sdsz,bsdindx,& + & iarcv(1:iszr),rvsz,brvindx,ictxt,info) + if (info == psb_success_) call psb_simple_a2av(jasnd,sdsz,bsdindx,& + & jarcv(1:iszr),rvsz,brvindx,ictxt,info) + case(psb_sp_a2av_mpi_) + + call mpi_alltoallv(valsnd,sdsz,bsdindx,psb_mpi_r_dpk_,& + & acoo%val(nzl+1:nzl+iszr),rvsz,brvindx,psb_mpi_r_dpk_,icomm,minfo) + if (minfo == mpi_success) & + & call mpi_alltoallv(iasnd,sdsz,bsdindx,psb_mpi_lpk_,& + & iarcv(1:iszr),rvsz,brvindx,psb_mpi_lpk_,icomm,minfo) + if (minfo == mpi_success) & + & call mpi_alltoallv(jasnd,sdsz,bsdindx,psb_mpi_lpk_,& + & jarcv(1:iszr),rvsz,brvindx,psb_mpi_lpk_,icomm,minfo) + if (minfo /= mpi_success) info = minfo + case default + info = psb_err_internal_error_ + call psb_errpush(info,name,a_err='wrong A2AV alg selector') + goto 9999 + end select + + if (info /= psb_success_) then + info=psb_err_from_subroutine_ + call psb_errpush(info,name,a_err='alltoallv') + goto 9999 + end if + + if (debug_level >= psb_debug_outer_)& + & write(debug_unit,*) me,' ',trim(name),': Done alltoallv' + + if (present(desc_rx)) then + ! + ! Extend the appropriate descriptor; started as R but on + ! transpose it now describes C + ! + call desc_r%clone(desc_rx,info) + call psb_cd_reinit(desc_rx,info) + ! + ! Take out non-local rows + ! + call psb_glob_to_loc(iarcv(1:iszr),p_desc_c,info,iact='I',owned=.true.) + lszr = iszr + call psb_coo_clean_negidx_inner(lszr,iarcv(1:iszr),jarcv(1:iszr),& + & acoo%val(nzl+1:nzl+iszr),nzt,info) + call desc_rx%g2lip_ins(jarcv(1:nzt),info) + call psb_cdasb(desc_rx,info) + acoo%ia(nzl+1:nzl+nzt) = iarcv(1:nzt) + acoo%ja(nzl+1:nzl+nzt) = jarcv(1:nzt) + nzl = nzl + nzt + call acoo%set_nzeros(nzl) + nzl = acoo%get_nzeros() + call acoo%set_sorted(.false.) + ! + ! Insert to extend descriptor + ! + call acoo%set_nrows(p_desc_c%get_local_rows()) + call acoo%set_ncols(desc_rx%get_local_cols()) + !write(0,*) me,' Trans RX ',acoo%get_nrows(),acoo%get_ncols(),acoo%get_nzeros() else - call ain%mv_to_lcoo(atcoo,info) + ! + ! + ! Take out non-local rows + ! + call psb_glob_to_loc(iarcv(1:iszr),p_desc_c,info,iact='I',owned=.true.) + call psb_glob_to_loc(jarcv(1:iszr),desc_r,info,iact='I') + lszr = iszr + call psb_coo_clean_negidx_inner(lszr,iarcv(1:iszr),jarcv(1:iszr),& + & acoo%val(nzl+1:nzl+iszr),nzt,info) + acoo%ia(nzl+1:nzl+nzt) = iarcv(1:nzt) + acoo%ja(nzl+1:nzl+nzt) = jarcv(1:nzt) + nzl = nzl + nzt + call acoo%set_nzeros(nzl) + nzl = acoo%get_nzeros() + call acoo%set_sorted(.false.) + + call acoo%set_nrows(p_desc_c%get_local_rows()) + call acoo%set_ncols(desc_r%get_local_cols()) + !write(0,*) me,' Trans R- ',acoo%get_nrows(),acoo%get_ncols(),acoo%get_nzeros() end if - if (info == 0) call psb_glob_transpose(atcoo,desc_r,info,desc_c=desc_c,desc_rx=desc_rx) + +!!$ write(0,*) me,' Sanity check after rx%g2l :',count(acoo%ja(1:nzl)<0) + + + call acoo%fix(info) + nzl = acoo%get_nzeros() + if (present(atrans)) then - call atrans%mv_from_lcoo(atcoo,info) + call atrans%mv_from_coo(acoo,info) else - call ain%mv_from_lcoo(atcoo,info) + call ain%mv_from_coo(acoo,info) end if + Deallocate(brvindx,bsdindx,rvsz,sdsz,& + & iasnd,jasnd,valsnd,& + & iarcv,jarcv,stat=info) + if (debug_level >= psb_debug_outer_)& + & write(debug_unit,*) me,' ',trim(name),': Done' + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(ictxt,err_act) + + return end subroutine psb_d_coo_glob_transpose subroutine psb_d_simple_glob_transpose_ip(ain,desc_a,info) @@ -434,7 +707,7 @@ subroutine psb_d_simple_glob_transpose_ip(ain,desc_a,info) ! that the same DESC_A works for both A and A^T, which ! essentially means that A has a symmetric pattern. ! - type(psb_ld_coo_sparse_mat) :: tmpc1, tmpc2 + type(psb_d_coo_sparse_mat) :: tmpc1, tmpc2 integer(psb_ipk_) :: nz1, nz2, nzh, nz integer(psb_ipk_) :: ictxt, me, np integer(psb_lpk_) :: i, j, k, nrow, ncol, nlz @@ -485,7 +758,7 @@ subroutine psb_d_simple_glob_transpose(ain,aout,desc_a,info) ! that the same DESC_A works for both A and A^T, which ! essentially means that A has a symmetric pattern. ! - type(psb_ld_coo_sparse_mat) :: tmpc1, tmpc2 + type(psb_d_coo_sparse_mat) :: tmpc1, tmpc2 integer(psb_ipk_) :: nz1, nz2, nzh, nz integer(psb_ipk_) :: ictxt, me, np integer(psb_lpk_) :: i, j, k, nrow, ncol, nlz diff --git a/base/tools/psb_dspasb.f90 b/base/tools/psb_dspasb.f90 index c2434c09..542e6901 100644 --- a/base/tools/psb_dspasb.f90 +++ b/base/tools/psb_dspasb.f90 @@ -122,11 +122,11 @@ subroutine psb_dspasb(a,desc_a, info, afmt, upd, dupl, mold) end if - IF (debug_level >= psb_debug_ext_) then + if (debug_level >= psb_debug_ext_) then ch_err=a%get_fmt() write(debug_unit, *) me,' ',trim(name),': From SPCNV',& & info,' ',ch_err - end IF + end if if (psb_errstatus_fatal()) then info=psb_err_from_subroutine_ diff --git a/base/tools/psb_s_glob_transpose.F90 b/base/tools/psb_s_glob_transpose.F90 index 9b6fd434..875f92db 100644 --- a/base/tools/psb_s_glob_transpose.F90 +++ b/base/tools/psb_s_glob_transpose.F90 @@ -113,7 +113,6 @@ subroutine psb_ls_coo_glob_transpose(ain,desc_r,info,atrans,desc_c,desc_rx) integer(psb_ipk_) :: ictxt, np,me integer(psb_ipk_) :: counter,proc, err_act, j integer(psb_lpk_) :: i, k, idx, r, ipx,mat_recv, iszs, iszr,idxs,idxr,nz,& - & irmin,icmin,irmax,icmax,& & l1, nsnds, nrcvs, nr,nc,nzl, hlstart, nzt, nzd integer(psb_mpk_) :: icomm, minfo integer(psb_mpk_), allocatable :: brvindx(:), & @@ -190,7 +189,7 @@ subroutine psb_ls_coo_glob_transpose(ain,desc_r,info,atrans,desc_c,desc_rx) do k=1, nzl j = acoo%ja(k) if (j > hlstart) then - call p_desc_c%indxmap%fnd_halo_owner(j,proc,info) + call p_desc_c%indxmap%qry_halo_owner(j,proc,info) sdsz(proc+1) = sdsz(proc+1) +1 end if end do @@ -264,19 +263,19 @@ subroutine psb_ls_coo_glob_transpose(ain,desc_r,info,atrans,desc_c,desc_rx) acoo%ja(nzd) = acoo%ja(k) acoo%val(nzd) = acoo%val(k) else - call p_desc_c%indxmap%fnd_halo_owner(j,proc,info) + call p_desc_c%indxmap%qry_halo_owner(j,proc,info) tsdx(proc+1) = tsdx(proc+1) +1 iasnd(tsdx(proc+1)) = acoo%ia(k) jasnd(tsdx(proc+1)) = acoo%ja(k) valsnd(tsdx(proc+1)) = acoo%val(k) end if end do + call acoo%set_nzeros(nzd) ! ! Put halo entries in global numbering ! call desc_r%indxmap%l2gip(jasnd(1:iszs),info) call p_desc_c%indxmap%l2gip(iasnd(1:iszs),info) - call acoo%set_nzeros(nzd) ! And exchange data. ! Normally we'll use our SIMPLE A2AV and not MPI, because ! the communication pattern is sparse, so ours is more @@ -407,19 +406,293 @@ subroutine psb_s_coo_glob_transpose(ain,desc_r,info,atrans,desc_c,desc_rx) type(psb_desc_type), intent(out), optional :: desc_rx integer(psb_ipk_), intent(out) :: info - type(psb_ls_coo_sparse_mat) :: atcoo + integer(psb_ipk_) :: ictxt, np,me + integer(psb_ipk_) :: counter,proc, err_act, j + integer(psb_ipk_) :: i, k, idx, r, ipx,mat_recv, iszs, iszr,idxs,idxr,nz,& + & l1, nsnds, nrcvs, nr,nc,nzl, hlstart, nzd + integer(psb_lpk_) :: nzt, lszr + integer(psb_mpk_) :: icomm, minfo + integer(psb_mpk_), allocatable :: brvindx(:), & + & rvsz(:), bsdindx(:), sdsz(:), tsdx(:), trvx(:) + integer(psb_ipk_), allocatable :: halo_owner(:) + integer(psb_lpk_), allocatable :: iasnd(:), jasnd(:),iarcv(:),jarcv(:) + real(psb_spk_), allocatable :: valsnd(:) + type(psb_s_coo_sparse_mat), allocatable :: acoo + logical :: rowcnv_,colcnv_,rowscale_,colscale_,outcol_glob_ + type(psb_desc_type), pointer :: p_desc_c + character(len=5) :: outfmt_ + integer(psb_ipk_) :: debug_level, debug_unit + character(len=20) :: name, ch_err - if (present(atrans)) then - call ain%cp_to_lcoo(atcoo,info) + if(psb_get_errstatus() /= 0) return + info=psb_success_ + name='mld_glob_transpose' + call psb_erractionsave(err_act) + if (psb_errstatus_fatal()) then + info = psb_err_internal_error_ ; goto 9999 + end if + debug_unit = psb_get_debug_unit() + debug_level = psb_get_debug_level() + + ictxt = desc_r%get_context() + icomm = desc_r%get_mpic() + + Call psb_info(ictxt, me, np) + + if (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),': Start' + + if (present(desc_c)) then + p_desc_c => desc_c + else + p_desc_c => desc_r + end if + + Allocate(brvindx(np+1),& + & rvsz(np),sdsz(np),bsdindx(np+1), acoo,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 + l1 = 0 + brvindx(1) = 0 + bsdindx(1) = 0 + counter=1 + idx = 0 + idxs = 0 + idxr = 0 + + if (present(atrans)) then + call ain%cp_to_coo(acoo,info) + else + call ain%mv_to_coo(acoo,info) + end if + + + ! + ! Compute number of entries in the + ! halo part, sorted by destination process + ! + nr = desc_r%get_local_rows() + nc = p_desc_c%get_local_cols() + nzl = acoo%get_nzeros() + hlstart = p_desc_c%get_local_rows() + do k=1, nzl + j = acoo%ja(k) + if (j > hlstart) then + call p_desc_c%indxmap%qry_halo_owner(j,proc,info) + sdsz(proc+1) = sdsz(proc+1) +1 + end if + end do + + ! + ! Exchange sizes, so as to know sends/receives. + ! This is different from the halo exchange because the + ! number of entries was not precomputed in the descriptor, + ! which was vector-oriented and not matrix-entry-oriented + ! + call mpi_alltoall(sdsz,1,psb_mpi_mpk_,& + & rvsz,1,psb_mpi_mpk_,icomm,minfo) + + if (info /= psb_success_) then + info=psb_err_from_subroutine_ + call psb_errpush(info,name,a_err='mpi_alltoall') + goto 9999 + end if + nsnds = count(sdsz /= 0) + nrcvs = count(rvsz /= 0) + if (debug_level >= psb_debug_outer_)& + & write(debug_unit,*) me,' ',trim(name),': Done initial alltoall',nsnds,nrcvs + + idxs = 0 + idxr = 0 + counter = 1 + Do proc = 0, np-1 + bsdindx(proc+1) = idxs + idxs = idxs + sdsz(proc+1) + brvindx(proc+1) = idxr + idxr = idxr + rvsz(proc+1) + Enddo + + tsdx = bsdindx + trvx = brvindx + + iszr = sum(rvsz) + iszs = sum(sdsz) + + if (debug_level >= psb_debug_outer_)& + & write(debug_unit,*) me,' ',trim(name),': Sizes:',& + & ' Send:',sdsz(:),' Receive:',rvsz(:) + + call psb_ensure_size(max(iszs,1),iasnd,info) + if (info == psb_success_) call psb_ensure_size(max(iszs,1),jasnd,info) + if (info == psb_success_) call psb_ensure_size(max(iszs,1),valsnd,info) + if (info == psb_success_) call psb_ensure_size(max(iszr,1),iarcv,info) + if (info == psb_success_) call psb_ensure_size(max(iszr,1),jarcv,info) + if (info /= psb_success_) then + info=psb_err_from_subroutine_ + call psb_errpush(info,name,a_err='ensure_size') + goto 9999 + end if + + ! + ! Now, transpose the matrix, then split between itself + ! and the send buffers + ! + call acoo%transp() + if (acoo%get_nzeros()/= nzl) then + write(0,*) me,'Something strange upon transpose: ',nzl,acoo%get_nzeros() + end if + nzl = acoo%get_nzeros() + hlstart = p_desc_c%get_local_rows() + + + nzd = 0 + do k = 1, nzl + j = acoo%ia(k) + if (j<=hlstart) then + nzd = nzd + 1 + acoo%ia(nzd) = acoo%ia(k) + acoo%ja(nzd) = acoo%ja(k) + acoo%val(nzd) = acoo%val(k) + else + call p_desc_c%indxmap%qry_halo_owner(j,proc,info) + tsdx(proc+1) = tsdx(proc+1) +1 + iasnd(tsdx(proc+1)) = acoo%ia(k) + jasnd(tsdx(proc+1)) = acoo%ja(k) + valsnd(tsdx(proc+1)) = acoo%val(k) + end if + end do + call acoo%set_nzeros(nzd) + ! + ! Put halo entries in global numbering + ! + call desc_r%indxmap%l2gip(jasnd(1:iszs),info) + call p_desc_c%indxmap%l2gip(iasnd(1:iszs),info) + ! And exchange data. + ! Normally we'll use our SIMPLE A2AV and not MPI, because + ! the communication pattern is sparse, so ours is more + ! efficient. Using ACOO for the receive buffers. + nzl = acoo%get_nzeros() + call acoo%ensure_size(nzl+iszr) + + select case(psb_get_sp_a2av_alg()) + case(psb_sp_a2av_smpl_triad_) + call psb_simple_triad_a2av(valsnd,iasnd,jasnd,sdsz,bsdindx,& + & acoo%val(nzl+1:nzl+iszr),iarcv(1:iszr),& + & jarcv(1:iszr),rvsz,brvindx,ictxt,info) + case(psb_sp_a2av_smpl_v_) + call psb_simple_a2av(valsnd,sdsz,bsdindx,& + & acoo%val(nzl+1:nzl+iszr),rvsz,brvindx,ictxt,info) + if (info == psb_success_) call psb_simple_a2av(iasnd,sdsz,bsdindx,& + & iarcv(1:iszr),rvsz,brvindx,ictxt,info) + if (info == psb_success_) call psb_simple_a2av(jasnd,sdsz,bsdindx,& + & jarcv(1:iszr),rvsz,brvindx,ictxt,info) + case(psb_sp_a2av_mpi_) + + call mpi_alltoallv(valsnd,sdsz,bsdindx,psb_mpi_r_spk_,& + & acoo%val(nzl+1:nzl+iszr),rvsz,brvindx,psb_mpi_r_spk_,icomm,minfo) + if (minfo == mpi_success) & + & call mpi_alltoallv(iasnd,sdsz,bsdindx,psb_mpi_lpk_,& + & iarcv(1:iszr),rvsz,brvindx,psb_mpi_lpk_,icomm,minfo) + if (minfo == mpi_success) & + & call mpi_alltoallv(jasnd,sdsz,bsdindx,psb_mpi_lpk_,& + & jarcv(1:iszr),rvsz,brvindx,psb_mpi_lpk_,icomm,minfo) + if (minfo /= mpi_success) info = minfo + case default + info = psb_err_internal_error_ + call psb_errpush(info,name,a_err='wrong A2AV alg selector') + goto 9999 + end select + + if (info /= psb_success_) then + info=psb_err_from_subroutine_ + call psb_errpush(info,name,a_err='alltoallv') + goto 9999 + end if + + if (debug_level >= psb_debug_outer_)& + & write(debug_unit,*) me,' ',trim(name),': Done alltoallv' + + if (present(desc_rx)) then + ! + ! Extend the appropriate descriptor; started as R but on + ! transpose it now describes C + ! + call desc_r%clone(desc_rx,info) + call psb_cd_reinit(desc_rx,info) + ! + ! Take out non-local rows + ! + call psb_glob_to_loc(iarcv(1:iszr),p_desc_c,info,iact='I',owned=.true.) + lszr = iszr + call psb_coo_clean_negidx_inner(lszr,iarcv(1:iszr),jarcv(1:iszr),& + & acoo%val(nzl+1:nzl+iszr),nzt,info) + call desc_rx%g2lip_ins(jarcv(1:nzt),info) + call psb_cdasb(desc_rx,info) + acoo%ia(nzl+1:nzl+nzt) = iarcv(1:nzt) + acoo%ja(nzl+1:nzl+nzt) = jarcv(1:nzt) + nzl = nzl + nzt + call acoo%set_nzeros(nzl) + nzl = acoo%get_nzeros() + call acoo%set_sorted(.false.) + ! + ! Insert to extend descriptor + ! + call acoo%set_nrows(p_desc_c%get_local_rows()) + call acoo%set_ncols(desc_rx%get_local_cols()) + !write(0,*) me,' Trans RX ',acoo%get_nrows(),acoo%get_ncols(),acoo%get_nzeros() else - call ain%mv_to_lcoo(atcoo,info) + ! + ! + ! Take out non-local rows + ! + call psb_glob_to_loc(iarcv(1:iszr),p_desc_c,info,iact='I',owned=.true.) + call psb_glob_to_loc(jarcv(1:iszr),desc_r,info,iact='I') + lszr = iszr + call psb_coo_clean_negidx_inner(lszr,iarcv(1:iszr),jarcv(1:iszr),& + & acoo%val(nzl+1:nzl+iszr),nzt,info) + acoo%ia(nzl+1:nzl+nzt) = iarcv(1:nzt) + acoo%ja(nzl+1:nzl+nzt) = jarcv(1:nzt) + nzl = nzl + nzt + call acoo%set_nzeros(nzl) + nzl = acoo%get_nzeros() + call acoo%set_sorted(.false.) + + call acoo%set_nrows(p_desc_c%get_local_rows()) + call acoo%set_ncols(desc_r%get_local_cols()) + !write(0,*) me,' Trans R- ',acoo%get_nrows(),acoo%get_ncols(),acoo%get_nzeros() end if - if (info == 0) call psb_glob_transpose(atcoo,desc_r,info,desc_c=desc_c,desc_rx=desc_rx) + +!!$ write(0,*) me,' Sanity check after rx%g2l :',count(acoo%ja(1:nzl)<0) + + + call acoo%fix(info) + nzl = acoo%get_nzeros() + if (present(atrans)) then - call atrans%mv_from_lcoo(atcoo,info) + call atrans%mv_from_coo(acoo,info) else - call ain%mv_from_lcoo(atcoo,info) + call ain%mv_from_coo(acoo,info) end if + Deallocate(brvindx,bsdindx,rvsz,sdsz,& + & iasnd,jasnd,valsnd,& + & iarcv,jarcv,stat=info) + if (debug_level >= psb_debug_outer_)& + & write(debug_unit,*) me,' ',trim(name),': Done' + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(ictxt,err_act) + + return end subroutine psb_s_coo_glob_transpose subroutine psb_s_simple_glob_transpose_ip(ain,desc_a,info) @@ -434,7 +707,7 @@ subroutine psb_s_simple_glob_transpose_ip(ain,desc_a,info) ! that the same DESC_A works for both A and A^T, which ! essentially means that A has a symmetric pattern. ! - type(psb_ls_coo_sparse_mat) :: tmpc1, tmpc2 + type(psb_s_coo_sparse_mat) :: tmpc1, tmpc2 integer(psb_ipk_) :: nz1, nz2, nzh, nz integer(psb_ipk_) :: ictxt, me, np integer(psb_lpk_) :: i, j, k, nrow, ncol, nlz @@ -485,7 +758,7 @@ subroutine psb_s_simple_glob_transpose(ain,aout,desc_a,info) ! that the same DESC_A works for both A and A^T, which ! essentially means that A has a symmetric pattern. ! - type(psb_ls_coo_sparse_mat) :: tmpc1, tmpc2 + type(psb_s_coo_sparse_mat) :: tmpc1, tmpc2 integer(psb_ipk_) :: nz1, nz2, nzh, nz integer(psb_ipk_) :: ictxt, me, np integer(psb_lpk_) :: i, j, k, nrow, ncol, nlz diff --git a/base/tools/psb_sspasb.f90 b/base/tools/psb_sspasb.f90 index 58241a92..01187497 100644 --- a/base/tools/psb_sspasb.f90 +++ b/base/tools/psb_sspasb.f90 @@ -122,11 +122,11 @@ subroutine psb_sspasb(a,desc_a, info, afmt, upd, dupl, mold) end if - IF (debug_level >= psb_debug_ext_) then + if (debug_level >= psb_debug_ext_) then ch_err=a%get_fmt() write(debug_unit, *) me,' ',trim(name),': From SPCNV',& & info,' ',ch_err - end IF + end if if (psb_errstatus_fatal()) then info=psb_err_from_subroutine_ diff --git a/base/tools/psb_z_glob_transpose.F90 b/base/tools/psb_z_glob_transpose.F90 index 1d087bd3..8f7cadd5 100644 --- a/base/tools/psb_z_glob_transpose.F90 +++ b/base/tools/psb_z_glob_transpose.F90 @@ -113,7 +113,6 @@ subroutine psb_lz_coo_glob_transpose(ain,desc_r,info,atrans,desc_c,desc_rx) integer(psb_ipk_) :: ictxt, np,me integer(psb_ipk_) :: counter,proc, err_act, j integer(psb_lpk_) :: i, k, idx, r, ipx,mat_recv, iszs, iszr,idxs,idxr,nz,& - & irmin,icmin,irmax,icmax,& & l1, nsnds, nrcvs, nr,nc,nzl, hlstart, nzt, nzd integer(psb_mpk_) :: icomm, minfo integer(psb_mpk_), allocatable :: brvindx(:), & @@ -190,7 +189,7 @@ subroutine psb_lz_coo_glob_transpose(ain,desc_r,info,atrans,desc_c,desc_rx) do k=1, nzl j = acoo%ja(k) if (j > hlstart) then - call p_desc_c%indxmap%fnd_halo_owner(j,proc,info) + call p_desc_c%indxmap%qry_halo_owner(j,proc,info) sdsz(proc+1) = sdsz(proc+1) +1 end if end do @@ -264,19 +263,19 @@ subroutine psb_lz_coo_glob_transpose(ain,desc_r,info,atrans,desc_c,desc_rx) acoo%ja(nzd) = acoo%ja(k) acoo%val(nzd) = acoo%val(k) else - call p_desc_c%indxmap%fnd_halo_owner(j,proc,info) + call p_desc_c%indxmap%qry_halo_owner(j,proc,info) tsdx(proc+1) = tsdx(proc+1) +1 iasnd(tsdx(proc+1)) = acoo%ia(k) jasnd(tsdx(proc+1)) = acoo%ja(k) valsnd(tsdx(proc+1)) = acoo%val(k) end if end do + call acoo%set_nzeros(nzd) ! ! Put halo entries in global numbering ! call desc_r%indxmap%l2gip(jasnd(1:iszs),info) call p_desc_c%indxmap%l2gip(iasnd(1:iszs),info) - call acoo%set_nzeros(nzd) ! And exchange data. ! Normally we'll use our SIMPLE A2AV and not MPI, because ! the communication pattern is sparse, so ours is more @@ -407,19 +406,293 @@ subroutine psb_z_coo_glob_transpose(ain,desc_r,info,atrans,desc_c,desc_rx) type(psb_desc_type), intent(out), optional :: desc_rx integer(psb_ipk_), intent(out) :: info - type(psb_lz_coo_sparse_mat) :: atcoo + integer(psb_ipk_) :: ictxt, np,me + integer(psb_ipk_) :: counter,proc, err_act, j + integer(psb_ipk_) :: i, k, idx, r, ipx,mat_recv, iszs, iszr,idxs,idxr,nz,& + & l1, nsnds, nrcvs, nr,nc,nzl, hlstart, nzd + integer(psb_lpk_) :: nzt, lszr + integer(psb_mpk_) :: icomm, minfo + integer(psb_mpk_), allocatable :: brvindx(:), & + & rvsz(:), bsdindx(:), sdsz(:), tsdx(:), trvx(:) + integer(psb_ipk_), allocatable :: halo_owner(:) + integer(psb_lpk_), allocatable :: iasnd(:), jasnd(:),iarcv(:),jarcv(:) + complex(psb_dpk_), allocatable :: valsnd(:) + type(psb_z_coo_sparse_mat), allocatable :: acoo + logical :: rowcnv_,colcnv_,rowscale_,colscale_,outcol_glob_ + type(psb_desc_type), pointer :: p_desc_c + character(len=5) :: outfmt_ + integer(psb_ipk_) :: debug_level, debug_unit + character(len=20) :: name, ch_err - if (present(atrans)) then - call ain%cp_to_lcoo(atcoo,info) + if(psb_get_errstatus() /= 0) return + info=psb_success_ + name='mld_glob_transpose' + call psb_erractionsave(err_act) + if (psb_errstatus_fatal()) then + info = psb_err_internal_error_ ; goto 9999 + end if + debug_unit = psb_get_debug_unit() + debug_level = psb_get_debug_level() + + ictxt = desc_r%get_context() + icomm = desc_r%get_mpic() + + Call psb_info(ictxt, me, np) + + if (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),': Start' + + if (present(desc_c)) then + p_desc_c => desc_c + else + p_desc_c => desc_r + end if + + Allocate(brvindx(np+1),& + & rvsz(np),sdsz(np),bsdindx(np+1), acoo,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 + l1 = 0 + brvindx(1) = 0 + bsdindx(1) = 0 + counter=1 + idx = 0 + idxs = 0 + idxr = 0 + + if (present(atrans)) then + call ain%cp_to_coo(acoo,info) + else + call ain%mv_to_coo(acoo,info) + end if + + + ! + ! Compute number of entries in the + ! halo part, sorted by destination process + ! + nr = desc_r%get_local_rows() + nc = p_desc_c%get_local_cols() + nzl = acoo%get_nzeros() + hlstart = p_desc_c%get_local_rows() + do k=1, nzl + j = acoo%ja(k) + if (j > hlstart) then + call p_desc_c%indxmap%qry_halo_owner(j,proc,info) + sdsz(proc+1) = sdsz(proc+1) +1 + end if + end do + + ! + ! Exchange sizes, so as to know sends/receives. + ! This is different from the halo exchange because the + ! number of entries was not precomputed in the descriptor, + ! which was vector-oriented and not matrix-entry-oriented + ! + call mpi_alltoall(sdsz,1,psb_mpi_mpk_,& + & rvsz,1,psb_mpi_mpk_,icomm,minfo) + + if (info /= psb_success_) then + info=psb_err_from_subroutine_ + call psb_errpush(info,name,a_err='mpi_alltoall') + goto 9999 + end if + nsnds = count(sdsz /= 0) + nrcvs = count(rvsz /= 0) + if (debug_level >= psb_debug_outer_)& + & write(debug_unit,*) me,' ',trim(name),': Done initial alltoall',nsnds,nrcvs + + idxs = 0 + idxr = 0 + counter = 1 + Do proc = 0, np-1 + bsdindx(proc+1) = idxs + idxs = idxs + sdsz(proc+1) + brvindx(proc+1) = idxr + idxr = idxr + rvsz(proc+1) + Enddo + + tsdx = bsdindx + trvx = brvindx + + iszr = sum(rvsz) + iszs = sum(sdsz) + + if (debug_level >= psb_debug_outer_)& + & write(debug_unit,*) me,' ',trim(name),': Sizes:',& + & ' Send:',sdsz(:),' Receive:',rvsz(:) + + call psb_ensure_size(max(iszs,1),iasnd,info) + if (info == psb_success_) call psb_ensure_size(max(iszs,1),jasnd,info) + if (info == psb_success_) call psb_ensure_size(max(iszs,1),valsnd,info) + if (info == psb_success_) call psb_ensure_size(max(iszr,1),iarcv,info) + if (info == psb_success_) call psb_ensure_size(max(iszr,1),jarcv,info) + if (info /= psb_success_) then + info=psb_err_from_subroutine_ + call psb_errpush(info,name,a_err='ensure_size') + goto 9999 + end if + + ! + ! Now, transpose the matrix, then split between itself + ! and the send buffers + ! + call acoo%transp() + if (acoo%get_nzeros()/= nzl) then + write(0,*) me,'Something strange upon transpose: ',nzl,acoo%get_nzeros() + end if + nzl = acoo%get_nzeros() + hlstart = p_desc_c%get_local_rows() + + + nzd = 0 + do k = 1, nzl + j = acoo%ia(k) + if (j<=hlstart) then + nzd = nzd + 1 + acoo%ia(nzd) = acoo%ia(k) + acoo%ja(nzd) = acoo%ja(k) + acoo%val(nzd) = acoo%val(k) + else + call p_desc_c%indxmap%qry_halo_owner(j,proc,info) + tsdx(proc+1) = tsdx(proc+1) +1 + iasnd(tsdx(proc+1)) = acoo%ia(k) + jasnd(tsdx(proc+1)) = acoo%ja(k) + valsnd(tsdx(proc+1)) = acoo%val(k) + end if + end do + call acoo%set_nzeros(nzd) + ! + ! Put halo entries in global numbering + ! + call desc_r%indxmap%l2gip(jasnd(1:iszs),info) + call p_desc_c%indxmap%l2gip(iasnd(1:iszs),info) + ! And exchange data. + ! Normally we'll use our SIMPLE A2AV and not MPI, because + ! the communication pattern is sparse, so ours is more + ! efficient. Using ACOO for the receive buffers. + nzl = acoo%get_nzeros() + call acoo%ensure_size(nzl+iszr) + + select case(psb_get_sp_a2av_alg()) + case(psb_sp_a2av_smpl_triad_) + call psb_simple_triad_a2av(valsnd,iasnd,jasnd,sdsz,bsdindx,& + & acoo%val(nzl+1:nzl+iszr),iarcv(1:iszr),& + & jarcv(1:iszr),rvsz,brvindx,ictxt,info) + case(psb_sp_a2av_smpl_v_) + call psb_simple_a2av(valsnd,sdsz,bsdindx,& + & acoo%val(nzl+1:nzl+iszr),rvsz,brvindx,ictxt,info) + if (info == psb_success_) call psb_simple_a2av(iasnd,sdsz,bsdindx,& + & iarcv(1:iszr),rvsz,brvindx,ictxt,info) + if (info == psb_success_) call psb_simple_a2av(jasnd,sdsz,bsdindx,& + & jarcv(1:iszr),rvsz,brvindx,ictxt,info) + case(psb_sp_a2av_mpi_) + + call mpi_alltoallv(valsnd,sdsz,bsdindx,psb_mpi_c_dpk_,& + & acoo%val(nzl+1:nzl+iszr),rvsz,brvindx,psb_mpi_c_dpk_,icomm,minfo) + if (minfo == mpi_success) & + & call mpi_alltoallv(iasnd,sdsz,bsdindx,psb_mpi_lpk_,& + & iarcv(1:iszr),rvsz,brvindx,psb_mpi_lpk_,icomm,minfo) + if (minfo == mpi_success) & + & call mpi_alltoallv(jasnd,sdsz,bsdindx,psb_mpi_lpk_,& + & jarcv(1:iszr),rvsz,brvindx,psb_mpi_lpk_,icomm,minfo) + if (minfo /= mpi_success) info = minfo + case default + info = psb_err_internal_error_ + call psb_errpush(info,name,a_err='wrong A2AV alg selector') + goto 9999 + end select + + if (info /= psb_success_) then + info=psb_err_from_subroutine_ + call psb_errpush(info,name,a_err='alltoallv') + goto 9999 + end if + + if (debug_level >= psb_debug_outer_)& + & write(debug_unit,*) me,' ',trim(name),': Done alltoallv' + + if (present(desc_rx)) then + ! + ! Extend the appropriate descriptor; started as R but on + ! transpose it now describes C + ! + call desc_r%clone(desc_rx,info) + call psb_cd_reinit(desc_rx,info) + ! + ! Take out non-local rows + ! + call psb_glob_to_loc(iarcv(1:iszr),p_desc_c,info,iact='I',owned=.true.) + lszr = iszr + call psb_coo_clean_negidx_inner(lszr,iarcv(1:iszr),jarcv(1:iszr),& + & acoo%val(nzl+1:nzl+iszr),nzt,info) + call desc_rx%g2lip_ins(jarcv(1:nzt),info) + call psb_cdasb(desc_rx,info) + acoo%ia(nzl+1:nzl+nzt) = iarcv(1:nzt) + acoo%ja(nzl+1:nzl+nzt) = jarcv(1:nzt) + nzl = nzl + nzt + call acoo%set_nzeros(nzl) + nzl = acoo%get_nzeros() + call acoo%set_sorted(.false.) + ! + ! Insert to extend descriptor + ! + call acoo%set_nrows(p_desc_c%get_local_rows()) + call acoo%set_ncols(desc_rx%get_local_cols()) + !write(0,*) me,' Trans RX ',acoo%get_nrows(),acoo%get_ncols(),acoo%get_nzeros() else - call ain%mv_to_lcoo(atcoo,info) + ! + ! + ! Take out non-local rows + ! + call psb_glob_to_loc(iarcv(1:iszr),p_desc_c,info,iact='I',owned=.true.) + call psb_glob_to_loc(jarcv(1:iszr),desc_r,info,iact='I') + lszr = iszr + call psb_coo_clean_negidx_inner(lszr,iarcv(1:iszr),jarcv(1:iszr),& + & acoo%val(nzl+1:nzl+iszr),nzt,info) + acoo%ia(nzl+1:nzl+nzt) = iarcv(1:nzt) + acoo%ja(nzl+1:nzl+nzt) = jarcv(1:nzt) + nzl = nzl + nzt + call acoo%set_nzeros(nzl) + nzl = acoo%get_nzeros() + call acoo%set_sorted(.false.) + + call acoo%set_nrows(p_desc_c%get_local_rows()) + call acoo%set_ncols(desc_r%get_local_cols()) + !write(0,*) me,' Trans R- ',acoo%get_nrows(),acoo%get_ncols(),acoo%get_nzeros() end if - if (info == 0) call psb_glob_transpose(atcoo,desc_r,info,desc_c=desc_c,desc_rx=desc_rx) + +!!$ write(0,*) me,' Sanity check after rx%g2l :',count(acoo%ja(1:nzl)<0) + + + call acoo%fix(info) + nzl = acoo%get_nzeros() + if (present(atrans)) then - call atrans%mv_from_lcoo(atcoo,info) + call atrans%mv_from_coo(acoo,info) else - call ain%mv_from_lcoo(atcoo,info) + call ain%mv_from_coo(acoo,info) end if + Deallocate(brvindx,bsdindx,rvsz,sdsz,& + & iasnd,jasnd,valsnd,& + & iarcv,jarcv,stat=info) + if (debug_level >= psb_debug_outer_)& + & write(debug_unit,*) me,' ',trim(name),': Done' + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(ictxt,err_act) + + return end subroutine psb_z_coo_glob_transpose subroutine psb_z_simple_glob_transpose_ip(ain,desc_a,info) @@ -434,7 +707,7 @@ subroutine psb_z_simple_glob_transpose_ip(ain,desc_a,info) ! that the same DESC_A works for both A and A^T, which ! essentially means that A has a symmetric pattern. ! - type(psb_lz_coo_sparse_mat) :: tmpc1, tmpc2 + type(psb_z_coo_sparse_mat) :: tmpc1, tmpc2 integer(psb_ipk_) :: nz1, nz2, nzh, nz integer(psb_ipk_) :: ictxt, me, np integer(psb_lpk_) :: i, j, k, nrow, ncol, nlz @@ -485,7 +758,7 @@ subroutine psb_z_simple_glob_transpose(ain,aout,desc_a,info) ! that the same DESC_A works for both A and A^T, which ! essentially means that A has a symmetric pattern. ! - type(psb_lz_coo_sparse_mat) :: tmpc1, tmpc2 + type(psb_z_coo_sparse_mat) :: tmpc1, tmpc2 integer(psb_ipk_) :: nz1, nz2, nzh, nz integer(psb_ipk_) :: ictxt, me, np integer(psb_lpk_) :: i, j, k, nrow, ncol, nlz diff --git a/base/tools/psb_zspasb.f90 b/base/tools/psb_zspasb.f90 index 9db28550..6cdfc61f 100644 --- a/base/tools/psb_zspasb.f90 +++ b/base/tools/psb_zspasb.f90 @@ -122,11 +122,11 @@ subroutine psb_zspasb(a,desc_a, info, afmt, upd, dupl, mold) end if - IF (debug_level >= psb_debug_ext_) then + if (debug_level >= psb_debug_ext_) then ch_err=a%get_fmt() write(debug_unit, *) me,' ',trim(name),': From SPCNV',& & info,' ',ch_err - end IF + end if if (psb_errstatus_fatal()) then info=psb_err_from_subroutine_ diff --git a/docs/html/userhtmlsu2.html b/docs/html/userhtmlsu2.html index 9c452948..1a266fc2 100644 --- a/docs/html/userhtmlsu2.html +++ b/docs/html/userhtmlsu2.html @@ -114,7 +114,7 @@ class="description">Each process has its own value(s) independently. src="userhtml0x.png" alt="psb_version_string_ " class="math-display" >

whose current value is 3.4.0 +class="cmtt-10">3.7.0 diff --git a/docs/html/userhtmlsu3.html b/docs/html/userhtmlsu3.html index 15721cf1..e01fd793 100644 --- a/docs/html/userhtmlsu3.html +++ b/docs/html/userhtmlsu3.html @@ -128,11 +128,13 @@ href="userhtml10.html#fn3x0">3 .

  • Call the iterative method of choice, e.g. psb_bicgstab
    • -

    This is the structure of the sample programs in the directory Call the iterative driver psb_krylov with the method of choice, e.g. + bicgstab. +

    This is the structure of the sample programs in the directory test/pargen/. -

    For a simulation in which the same discretization mesh is used over multiple time +

    For a simulation in which the same discretization mesh is used over multiple time steps, the following structure may be more appropriate:

    1. prec%build class="enumerate" id="x9-6043x5">Call the iterative method of choice, e.g. psb_bicgstab
    -

    The insertion routines will be called as many times as needed; they only need to be +

    The insertion routines will be called as many times as needed; they only need to be called on the data that is actually allocated to the current process, i.e. each process generates its own data. -

    In principle there is no specific order in the calls to

    In principle there is no specific order in the calls to psb_spins, nor is there a requirement to build a matrix row in its entirety before calling the routine; this allows the application programmer to walk through the discretization mesh element by element, generating the main part of a given matrix row but also contributions to the rows corresponding to neighbouring elements. -

    From a functional point of view it is even possible to execute one call for each +

    From a functional point of view it is even possible to execute one call for each nonzero coefficient; however this would have a substantial computational overhead. It is therefore advisable to pack a certain amount of data into each call to the insertion routine, say touching on a few tens of rows; the best @@ -209,10 +211,10 @@ process and pass it in a single call to psb_spins; this, however, would entail a doubling of memory occupation, and thus would be almost always far from optimal. -

    +

    2.3.1 User-defined index mappings
    -

    PSBLAS supports user-defined global to local index mappings, subject to the +

    PSBLAS supports user-defined global to local index mappings, subject to the constraints outlined in sec. 2.3:

      @@ -230,7 +232,7 @@ class="cmmi-10">…ncol i;
    -

    but otherwise the mapping is arbitrary. The user application is responsible to ensure +

    but otherwise the mapping is arbitrary. The user application is responsible to ensure consistency of this mapping; some errors may be caught by the library, but this is not guaranteed. The application structure to support this usage is as follows: @@ -271,12 +273,12 @@ class="cmtt-10">irw, respectively, are already local indice -

    [

    [next] [prev] [prev-tail] [front] [up]

    -

    diff --git a/docs/html/userhtmlsu4.html b/docs/html/userhtmlsu4.html index ffc95516..62c023bf 100644 --- a/docs/html/userhtmlsu4.html +++ b/docs/html/userhtmlsu4.html @@ -10,23 +10,23 @@ -

    [

    [prev] [prev-tail] [tail] [up]

    2.4 Programming model

    -

    The PSBLAS librarary is based on the Single Program Multiple Data (SPMD) +

    The PSBLAS librarary is based on the Single Program Multiple Data (SPMD) programming model: each process participating in the computation performs the same actions on a chunk of data. Parallelism is thus data-driven. -

    Because of this structure, many subroutines coordinate their action across the +

    Because of this structure, many subroutines coordinate their action across the various processes, thus providing an implicit synchronization point, and therefore must be called simultaneously by all processes participating in the computation. This is certainly true for the data allocation and assembly routines, for all the computational routines and for some of the tools routines. -

    However there are many cases where no synchronization, and indeed no +

    However there are many cases where no synchronization, and indeed no communication among processes, is implied; for instance, all the routines in sec. 3 are only acting on the local data structures, and thus may be called independently. @@ -34,7 +34,7 @@ The most important case is that of the coefficient insertion routines: since the number of coefficients in the sparse and dense matrices varies among the processors, and since the user is free to choose an arbitrary order in builiding the matrix entries, these routines cannot imply a synchronization. -

    Throughout this user’s guide each subroutine will be clearly indicated +

    Throughout this user’s guide each subroutine will be clearly indicated as:

    > stream 0 g 0 G @@ -4801,7 +4801,7 @@ BT 0 g 0 G -69.503 -22.397 Td [(7.)]TJ 0 g 0 G - [-500(Call)-333(the)-334(iterativ)28(e)-333(metho)-28(d)-333(of)-334(c)28(hoice,)-333(e.g.)]TJ/F30 9.9626 Tf 189.595 0 Td [(psb_bicgstab)]TJ/F8 9.9626 Tf -201.772 -21.778 Td [(This)-333(is)-334(the)-333(structure)-333(of)-334(the)-333(sample)-333(programs)-334(in)-333(the)-333(directory)]TJ/F30 9.9626 Tf 269.435 0 Td [(test/pargen/)]TJ/F8 9.9626 Tf 62.764 0 Td [(.)]TJ -317.255 -12.573 Td [(F)83(or)-291(a)-292(sim)28(ulation)-292(in)-291(whic)27(h)-291(the)-292(same)-292(discretization)-291(mes)-1(h)-291(is)-292(used)-291(o)27(v)28(er)-292(m)28(ultiple)]TJ -14.944 -11.955 Td [(time)-333(ste)-1(p)1(s)-1(,)-333(the)-333(follo)28(wing)-334(structure)-333(ma)28(y)-333(b)-28(e)-334(more)-333(appropriate:)]TJ + [-500(Call)-222(the)-222(iterativ)27(e)-222(driv)28(er)]TJ/F30 9.9626 Tf 114.654 0 Td [(psb_krylov)]TJ/F8 9.9626 Tf 54.517 0 Td [(with)-222(the)-222(metho)-28(d)-222(of)-223(c)28(hoice,)-244(e.g.)]TJ/F30 9.9626 Tf 133.518 0 Td [(bicgstab)]TJ/F8 9.9626 Tf 41.843 0 Td [(.)]TJ -356.709 -21.778 Td [(This)-333(is)-334(the)-333(structure)-333(of)-334(the)-333(sample)-333(programs)-334(in)-333(the)-333(directory)]TJ/F30 9.9626 Tf 269.435 0 Td [(test/pargen/)]TJ/F8 9.9626 Tf 62.764 0 Td [(.)]TJ -317.255 -12.573 Td [(F)83(or)-291(a)-292(sim)28(ulation)-292(in)-291(whic)27(h)-291(the)-292(same)-292(discretization)-291(mes)-1(h)-291(is)-292(used)-291(o)27(v)28(er)-292(m)28(ultiple)]TJ -14.944 -11.955 Td [(time)-333(ste)-1(p)1(s)-1(,)-333(the)-333(follo)28(wing)-334(structure)-333(ma)28(y)-333(b)-28(e)-334(more)-333(appropriate:)]TJ 0 g 0 G 12.177 -21.779 Td [(1.)]TJ 0 g 0 G @@ -27548,10 +27548,10 @@ endstream endobj 1997 0 obj << -/Length1 2460 -/Length2 17290 +/Length1 2478 +/Length2 17406 /Length3 0 -/Length 19750 +/Length 19884 >> stream %!PS-AdobeFont-1.0: CMTT10 003.002 @@ -27571,7 +27571,7 @@ FontDirectory/CMTT10 known{/CMTT10 findfont dup/UniqueID known{dup 11 dict begin /FontType 1 def /FontMatrix [0.001 0 0 0.001 0 0 ]readonly def -/FontName /TJSMYH+CMTT10 def +/FontName /LTGWDM+CMTT10 def /FontBBox {-4 -233 537 696 }readonly def /PaintType 0 def /FontInfo 9 dict dup begin @@ -27642,6 +27642,7 @@ dup 43 /plus put dup 113 /q put dup 114 /r put dup 115 /s put +dup 55 /seven put dup 54 /six put dup 47 /slash put dup 116 /t put @@ -27692,46 +27693,39 @@ KV "h/KU%K~0[ P`CU/ݘM3ʽuǂ&gk7XpMǟA{C%u%۩|V#GLߊb{%OY3]Za_W!l?vyf*8*l4]e> zf(Y`ۘ8?HPىJ'.dڰ?X𤄳*y!v48Iv 1Lnဝ?hF,CO 2fn_Z,HcO냹nҧcH; n7ZBM]%ǖ6v &<ˋ9ݖdžx}-`{yq.*Za}Bub}ڧaƓNh \WK#aB۳;yROF˩^2X1 <)4mTUdը4 ȵc3P1O\'(L?]J{||-UOjk\|<Ju!I59٩_;Mz>T>K1W^F:ÖR6&beˇM$/YETR䙕ļxsp iP/YLa5{ K -yo,F'=Cnxl;`%mwoh`o(A\=#TЇZްT6ԞPHvMcwVES"}G$6D>?gx6O+=|d2˟;l >ƫ9΅ʥ1v`Tb]^ެj-YOC=镹jQIOnu&$iVJsй^8!? -6.D>LSh $;zQ -}A,R:Tc7ǁq 3*.<VU5ܮ 8 zЩbtpt$2{DT2*KJwnlR8vAZo`DLok}f.wi а+hI)uxs٦@w&ڿئжIܭ}h -Ak(C1<BO$ݿ&8 9_LZXz.F ~!TMzd -ƻx^wSzwA)K9`Q^4[,.#M&,3"_;-SHHpF4kCwTJZ&~>nu -i1$mHL&X6|rH :C.O56 "ஸRyu(in;27KQ(9ZH kUi -f"fM `3ޘʚKd|Y&ƖPhg?LZ ?f|=Y r.rB)e$K"<2?8Aѫccܘ!Jl  r~tSTLu'M\Aj+Ƅi 5OG8go^'.l aG7{59عsK$L?2{0Tl述Yń~1k1{SkZXq,RmQ:gPTS{*8W2D]H}S?)iBV5ķ4):9Wfcs;*VzJ™A>dvrhqX cb&:jm+[16dtW1Nī}gNӴ?'RՄ۽dXlNn')[)<_$LoyT2x w Oo$‹]Ⱦ;F3Yo3 - " -?[mzlX- ޿R΅ QMo=Cm`!1 -YU#AyUZv8~ |3s&ˤJLn՟ǝ]?Zʺ}jF `i>{G A-& K }(x]d1 :dm~ dߊċMo׃v%ۿYsD 5xA8+I귩IIpYat!&XmE<.D`ьZ(gVB]-㺁tjc!k(匀JK'gc14cu?;oTo''E֣"xL_ =$Icĉb WU"HꢟĘsGo]k/{,]G+78&vXi:B,!ѐG"*qDŽQʨ!24_!/t"k?h0>u| zՖclŽ!uÍqsR5H$puYqXo#R_\5_?!IgF)/eau'AVdr}Y'kfVu9 !p죍e S`*kBkO}<|OA)C )ʝBങ(t* -l !2N/E>`S,RQ#bM,TYb}+`6kGe9!6$?]hJ.?Y @w԰dlptgURɡ-(r _,R‚nuX)"V@.2*I/`>/ {V97 >eŀMN\N#/ -L|$!d JKؐ}-T;0_, k>!:1OSt0J@uGXSBgc/ԃ N<^9 Ք\I4jՓ9;k a:EM -Ple+pkTl^,Xe_UJe@y}+blKWZ:,cyf0\0ZSE>˧=zpG\y+[[yZt3ߍ.xci^;'ArM;]=8?V6yʿUGF-m4@wc"VaU-#cHd#]EJ:I/ï]S0_Cٍm0 W,Di$yw*˕ck#:G hj-RXvtd?XKakZ|^2biBMd0F(glӆIcԙKYV:_g7. >i@3&үI~SMeBa򢳕=܇ O8Gq A;8alU`>v_#*A" (N?Yz{ZU'4."Qkoswa :7F{ *hw,btzdth.}_p:\lGo]v%J4 7Y=}`%WgE -GuEDԥ;d!/1kIfYYZ-Mi]$!r -LTV=o -۽P߸Bx=$-6){wk(@øOHUB ?rl6c>uҔIhկU}z$ߓ$j WB_= D#0 ExZiܤI}կjW|ƜڈMy-[؅٠ޥU$bxʄ&PznXAz^ |Y:˚pxSe9k -umޠDv -Ibף~)bX|O&g I'u!77X*() ݥ}iTl,Kͯl -ܹP/o=K5)k9 ,!P_"#A4얹l1CU_fEiNat u*H6VV F+D̙C~Zk+M - JT<5JŵE]+D>B8]8Th'/!7?TknG Ro7y>")Ɉ@aN J?$O@_,a4ƀy@+ ;`FZvǺֺeeG X~&eu`Oףehy@W)%u734T^vhά$ sv,$ 뷓xR| ,W]ץ#0*4>xf'FIc؋et7;e;˯y *R/r3}Ks9 6{W?i#UgeV}yݬLq#nʏƧWd1vۻEg]Èeq<RZ=~Wy4x*-%VVնQ -ЇMhaQdөy?fx>> QN[sv .Xpru -ݴ+I^DO ;6wۇQFoA8Oɻ)?A=+(*ajL2)㪼PR*C*(}syXq%+8x.Mڽ4cUb{3]*suSsr{AQ*k.V_a ;\rߒab2 :GO/>qF'R$5*:8~K[̟/Dh_%1}Vp}#\/9`S/ |Va؉s026,GWd=pHo6<2n^ibkʗu!sԳrb{vb)`7c<GN7%6$PB[N/f U@]k~'/_H7`"oNF^Jف}-c\žsY(A7,z*\֟_~ |O -hNN9vG'Ug4R1L}ݷ7Ar"s]D̜*f8>5#!]J@oxqrh߬ڹ jB  Jže)5薷71ӬQ^uܐ\#:<+1382ɜiB+WT^"a~J7Xi0`'_&Eg̦n3n eŃdTF<>p .2ĩ ?j٦} jXXcT7J B)Nc6+'d -Fm.]mD!]Yplk`J sTBVe졽. k90TʶņrKj,fh32zz9=d=ҡ9|4ophg\TI6Òҋ3ݪRT=$?ìAd`UƏ`ߌ${sѥ * .VlyƍK jfk yؒ:wi{zNكZ]d_"HJ?e`:0H-Q1vCt-T"w)׵{pδmF6gM_鴓0(K(m`Wkmy⯁+0ugIf>W -ƅϛ7-DiJG_)utU -",%ULڀbM!qxE~Als.4fH"\¸Ա _wgB8LMLUyasjҞӽ DȚ -7nOWHg~7I' -ȷ[MVw 7SMR4; J?jѿJ\덈7Kp (WcB3 @_zdxdUR .̺5msHIxIrt0Y"t/Y n.YUwAYʞa QI4sͷ`Zq{x73KS6FpJş`m,݊!iD'@).cAFkaG34Lkubu`wmS0-.)Z/ԥwhC -M-.Fx*Iw jO i㟧:eȶIſ:܊Hݧ-rh\'̒R_TgGP ZquSAA2,ᇉ^ϒO?靹Mw6BIxpDy,4ՓAG.WA'R(ZǏh8Ͽ?3']rmRKbw`wBSvgMqL&"d4E({qtf[A\ʻt1ǞRRy# dEtju Y*2?.T2>Gvwtj3^!E3廢XUZ _|lmE_r/-ϳٽ'ILRQd[a˴TC -,Ojʊb?כ:m~U]†^E_J ciG:OKá]QHqz2n_y#uB4oXL7:nT5Fˆ)L81C [bPN.u{2e$ ^%8=%I%u-^"m.۸~:k@jucd"kDmZTMc"Gtd$. $  */a+wz(SNN4jU@7Ku7$-椟xPU?v0&T t8'9==!E'рl>bd"}r hCAi\88;B8RpgFaT9%:$+Jhn ߢ&[K45?ۑR<~cLsoEpdអػTݰj\#2 ؚF)^ 5.6 -k4ݗي3J)bo߮k @I_wz -eHxw.=lO:Ї(=V{/k?Ib"9Rdxq+'ͺh'{1etb|ў!ݱ -s7M9-@ЮRB#:މףlV̌eUwJ"1 CW|M`/n$wΛ_7S `kssމo1,ؽII䝊s>9y U.&ˁ:.mQG-+NM91 +yo,F'=Cnxl;`%mwoh`o(A\=#TЇZްT6ԞPHvMcwVES"}G$6D>?gx6O+=|d2˟;l >ƫ9΅ʥ1v`Tb]^ެj-YOC=镹jQIOnu&$iVK)9!xz^lŔ޲|Ehp?zQiK3WZ >SLVDMj O<>qc·sB@?<8I;mT}u_.>)uƋz k:\"וvSA". !䳩ҦMn mc&lw HZ.sU)>JV؉+%# pbx=x%|)M= JpY઒PcXB)UW+ϭ؟[ͱ^h8ˢI+5qfÛikMuf&!9.S&H&m8b6oq:MHHH"'KU JsW/yAOXxsA͙Ȕjb0p9;@cN2V$[5$(i_֍7y-Jbi$\uFu)7Xnem'nLSB$Ii) 䁨nV'w(MؚJΣ +ʭ ~{nN;BD󃚼F:f~bk?ج԰Ӻ8w~,eS=k n#)57ҊvSHl"qj0+2͕Oj*p!? G!5wޑ#NU+uDeX*gFOCaoc[e?xhY4aZK,~,x8b!ׁ~pK5IIciE#i h Kw hy<1в_XRٛ1_:r^>zD/)h@܁0ѥ {C6'5 +'DC((!,ۖ3+à ]F8Kܞ!5dmΎx=lmnhK̉AhNY{- *I+A,O 2-ko̶U€rW٢":V2a-qWY(1ߤh1I + ϒB~dž\e.TDp6YФҕ0 "fh*fo%H#~=w[1)-.UI902AdeɄWP +>I] dZFg3"|9 3يRU-}à&T ! Ú17 bcb2Ѕ&l ަWaXmj {QM/T2D( G-.c_)#6- )W$qzEG}^뽴%K_D$/S'y|Ol~ߓ Zp#o(HX/D'nX{6{cZ/A M&]RGwBf`wv$49K2\K#2#!owj_AQ B!ZIfl!h<8N `”Ϗ#i2*9B|{IWZO.NNI٥Rg1_5=cO)`ZShYYAΓJZHgwҧ~*=}<O-Ʋ!1ϒsхd ҉i* + ?YO8P`+Y>E6ea6觐wڈ-.QYli 5c~U<ñ9, +͎{v*ax}N})ޒJVrwZ<b_zi]#1eš82VmuK俲&|cΛIK@QO)lb +W8w7pm=Ib- (1Lmx,Y3c$s!"e5sŕ."8(a +]bWvX?GM$tpZElA}gdzb|?ysgD@F_ +5ZR,$-h _éʱJQ eE̮_3Ϣ 0QﱟEG=c=5]6&n6gU͜8ԪV)~k{'VT :dƖLLm8NŐ9l(;EiQl0{,XA$jxꀡJ`?=(Ybc. 7lJ Ύybe,[Vs)3r3őLuNWګ~㻛sLusԘ|ٙg[DRE--0]5ubR/#N%pA'0;2XY%d(&8nQb(Q]֖W_" s6/$\NiVڄp%ugëKѶ6g_8u?RBc l7LZ!)IdgiщE9(~َ O3X K𨝳h_:h6`,$ ǭإ#j x+vbuKYEX|1ߥ!Gev^J[4= ]gi|H5΄ДHxa>Smi W}'^ɍsXA֭SAEi}F ŕ- pwQZErb4>-YƼx@@*+ g좮 +B !YR/_9#U%SU)pvtha._b z-©O +fy&0^Ƃ"|Wi$**OY N}QsTTI]96,3[[-W0(L;&}i)=M -U0> Qk 3"zN*m4= +#-Xf(Py eGzib7Fc(%%'?z"EBdbܒnn0`E\gIuwD*{tiK'hP 2 ѮS9#+sфu<=M~,{L_$pziμ^w8Gb`Dpm~ׄFk αwtC?^?܊yTA{Cu@+q[0XªfV/vAmo+n;. +ŘKmȼo"݁B1Vl&E#^2!ܠ5(K%V.u9"+ɘHl UG;$ Т sl~bUp:0g/+͟Ph C{)K`][bY "{9ԲF1&R:)WQt}A즕Līk]s၌Oh?8xY i"U7%A}SnTrǁ1=gj>0a`y-W%{@Ӻ6w.,[q f JfK,y|Ʃf!,]hN|~b|Di(M!|0 77ړ5t-RV&BݟnJrOUOnesE}OZҘl")e1*J6]#ĉިNyìS$gz fgl$_QUi1zc qҍwVgw쩵Joz @$ Ɠ g jV9Ѥ+o•gf,Vv!u/H["Z I,ziéHfLMfɋ )E>,3l=em)=R$Uu3G•~gIfSEqhһ!=iz6wNT꘍ċ;@Vwjx+y\~t%%6;Վ4Q=PD?8mI&lEدIr>?^a/(63en T1,dc jP}:k:E2KHN~L/}8N;~YeHHzDuv*֎'ZHzZ`Bt(":l+4jRF)60!+bO5q@7kh훑T_:XH Uu)h7b\\$Ɠ 0-D}Q5$>?*\˓9dtX%7떾&4gF}Ϫ@#_~۽gQ!w"MG sPŠ[æ1Ws7SڽbVZ,Sg 5w_%q;u~"ilNf S7:-J"ɛΜ kbb|y* İ/J$5?,JQ1!7xt`O)PPm޾v"}߰Vp$J !Q?bh#PvE ]-G\{L:_ALk0g<tWY/Mm^~Z_ࢣ'aaew F}I3}ɷYB<5, +X\kO4۷p +N6эH䝲f 8c5S?n"mb$d!oʶ*/$R[XUDAGaG3vwXEts*GZehbqĥ,V&zQ}t]<\I|:}vK=x)0>:v[aMHV[&,(E?mK-\yט<1NK(|[\z]R#t  z~Xzwǵ֜gOKPbt寓7$36Mz|{}5)y(>qP* o;;JkVxɞz.2Ð؇Fm6my|s6ӺSבq^yf͡p&8~6IzT!:AΔX.2K+X ʌ7 F9MoNW~+uקҔO] +\P/׬-kt|e ɶG vy+*٘9UwF2<Ȟ4GuAЖQ(jxQb O6K _"AK+e EpC2sg0B\gTRBe] ,B[T$ ~rU|Br8_ϲ׫bGZ8\;eKw8׷ +ȣ9*Vhj?m.Fp$kCEI/rNȖ BxM&JF +Q?anY3˦|*ZM퇣egG:r}f(t َ*љ;׽~YA KB+sxhafW!txXN_[\2{TD(IN\|:zg +M|?cs &r:,|/Pk^%Zz&;M/)<(=\ms?Sa4(qW%0jA໲H={D7ָ͡Ci@SLY\L*! $wZ̕'<0grDj$M~!U!.tzgP1s}Ĵd{q~ vF;D䈜e\_sJb$r*k_y~|5abXRq.RrDosf ά23Gp? hA3BrM94 +(( oʒތos=j(:@J5 v+-F7~Mղv'ׁ71=tG_HNEzNiK?BxNp-ƐPEJ9l kWF K2 ?͌[8/mN%!ȄL4$GtA-R"[sl 50}c@P8zyuU D w8xEyWli۬ڙz -%`*;b0UU?ohUxNWߚ@^b1@7>Y`eQ c?n(t%c8V\yEOMlǠIM*;S\Fw*-Ƽy+iy݉`t8TȡjJ.=hS4@U:q0qy 1NzDSmu0~,F( ' rBs~mb[5,f$4]ĔRH,*@'?n(+ 3\em \:!j[uvԛI[B]ײ f􇻨S2lYugnp@y5{gFNrYϥ^p`7Vr843jHz)3Ƀrê +!:nêS\w#UD&zNU>r8EZcRͮ<W)KƲ֨eVඨ&J#GpfoݶxD\ zE npA3tRoLšE(IVz>1v~|C#iE^d~ywGқl^4ScB2gLb}VQ|^\U!!j1>%CFGp +N_ܿ$D;dy#Y$Y$9:Qa9sMҸL~NcH=G &_l}@3ž9S`lv':etwR UMwh_2qTx y#X)a|M"!ipv?J p=ЀUit]LjJxSzp֕nݗ9?3 +NKty97]5QN7dׄP=LHjz@L{ކcϲ,YT+_ᶂJ)Qo=` +-цwӱ (V!Y.fvȈ'}پ't3 ~T(`].ǎiN½a! W K1h%dٓcjP&, +GJ- ^J|L'NIzv4FY?h-oavF^egtokN`*zjksPUb /=d5OmLMgUFe9,*} G;ha;tx0E-^f0eP>8ߕHCv ]\KZ\g"ף=#P#ۖCO[ӦKwqJܶ*?\]t+(5lH7Hxs(x@#m["`aɳ99ꘋ';a^v ~"t ;uF) W؟U{ dմcLZ@jJ5MKpH=Jb)єFH:4fs,j/MqR s+vҝm!(; »WdPa3U۩A%Zt.d7yo8!G%J'Of{{ +B ,0 I/g iSR\1+(ԭȴGF~NQiΓ^A6<N4Fޤ _?a 88-# o'9 YItH0D"PQC*q-e!~*t&{+N>2|I/=5z],FtMcٗf)s`Sk{o6 HR"|u̜ >1@{}@0 [)$UŬ4f&{X߽#|cCfkP+q#S+a @ߎu$ +-2[ :ÖȡB~tCM7(pBg;"mAgM:"k]:q>B 4nؒ݉QUr\Ӭijd=ɑRPo]4Rf !X@ bL,r<-nm 5׊L8(~^lꦟUxewo|[n[.كUid\SZvNm SyXd8P$~7,bgB f9Zs|",P?75mt^41ptNӟ.a Ăe%dk[`lx߬N%dΫa4vUZ_K+ endstream endobj 1999 0 obj @@ -27913,19 +27907,19 @@ endobj /Type /ObjStm /N 100 /First 958 -/Length 13084 +/Length 13090 >> stream -1991 0 1994 231 1996 604 1998 867 2000 1346 2002 1578 584 1862 582 2003 1753 2144 1726 2285 -797 2427 845 2568 813 2709 1911 2849 585 2989 815 3129 812 3267 817 3405 1301 3544 814 3684 -1171 3824 777 3963 583 4104 811 4245 873 4386 1013 4526 586 4667 778 4780 874 4893 929 5006 -964 5119 993 5232 1036 5349 1091 5469 1139 5589 1197 5709 1250 5829 1302 5949 1358 6069 1399 6189 -1437 6309 1486 6429 1528 6549 1562 6669 1602 6789 1643 6909 1676 7029 1707 7149 1744 7269 1806 7389 -1847 7509 1884 7629 1924 7749 2003 7851 2004 7966 2005 8087 2006 8208 2007 8329 2008 8423 2009 8519 -574 8588 570 8648 566 8759 562 8833 558 8921 554 9009 550 9097 546 9185 542 9259 538 9384 -534 9458 530 9546 526 9634 522 9722 518 9810 514 9884 510 10009 506 10083 502 10171 498 10259 -494 10333 490 10458 486 10532 482 10620 478 10708 474 10796 470 10884 466 10972 462 11060 458 11148 -454 11236 450 11324 446 11412 442 11500 438 11588 434 11676 430 11764 426 11852 422 11926 418 12052 +1991 0 1994 231 1996 604 1998 867 2000 1352 2002 1584 584 1868 582 2009 1753 2150 1726 2291 +797 2433 845 2574 813 2715 1911 2855 585 2995 815 3135 812 3273 817 3411 1301 3550 814 3690 +1171 3830 777 3969 583 4110 811 4251 873 4392 1013 4532 586 4673 778 4786 874 4899 929 5012 +964 5125 993 5238 1036 5355 1091 5475 1139 5595 1197 5715 1250 5835 1302 5955 1358 6075 1399 6195 +1437 6315 1486 6435 1528 6555 1562 6675 1602 6795 1643 6915 1676 7035 1707 7155 1744 7275 1806 7395 +1847 7515 1884 7635 1924 7755 2003 7857 2004 7972 2005 8093 2006 8214 2007 8335 2008 8429 2009 8525 +574 8594 570 8654 566 8765 562 8839 558 8927 554 9015 550 9103 546 9191 542 9265 538 9390 +534 9464 530 9552 526 9640 522 9728 518 9816 514 9890 510 10015 506 10089 502 10177 498 10265 +494 10339 490 10464 486 10538 482 10626 478 10714 474 10802 470 10890 466 10978 462 11066 458 11154 +454 11242 450 11330 446 11418 442 11506 438 11594 434 11682 430 11770 426 11858 422 11932 418 12058 % 1991 0 obj << /Type /FontDescriptor @@ -27974,7 +27968,7 @@ stream % 1998 0 obj << /Type /FontDescriptor -/FontName /TJSMYH+CMTT10 +/FontName /LTGWDM+CMTT10 /Flags 4 /FontBBox [-4 -233 537 696] /Ascent 611 @@ -27983,7 +27977,7 @@ stream /ItalicAngle 0 /StemV 69 /XHeight 431 -/CharSet (/A/B/C/D/E/F/I/K/L/M/N/O/P/R/S/T/U/W/Y/a/ampersand/asciitilde/asterisk/b/backslash/bracketleft/bracketright/c/colon/comma/d/e/equal/f/four/g/h/hyphen/i/j/k/l/m/n/o/one/p/parenleft/parenright/percent/period/plus/q/r/s/six/slash/t/three/two/u/underscore/v/w/x/y/z/zero) +/CharSet (/A/B/C/D/E/F/I/K/L/M/N/O/P/R/S/T/U/W/Y/a/ampersand/asciitilde/asterisk/b/backslash/bracketleft/bracketright/c/colon/comma/d/e/equal/f/four/g/h/hyphen/i/j/k/l/m/n/o/one/p/parenleft/parenright/percent/period/plus/q/r/s/seven/six/slash/t/three/two/u/underscore/v/w/x/y/z/zero) /FontFile 1997 0 R >> % 2000 0 obj @@ -28190,7 +28184,7 @@ stream << /Type /Font /Subtype /Type1 -/BaseFont /TJSMYH+CMTT10 +/BaseFont /LTGWDM+CMTT10 /FontDescriptor 1998 0 R /FirstChar 37 /LastChar 126 @@ -30143,8 +30137,8 @@ endobj 2131 0 obj << /Title (Parallel Sparse BLAS V. 3.7.0) /Subject (Parallel Sparse Basic Linear Algebra Subroutines) /Keywords (Computer Science Linear Algebra Fluid Dynamics Parallel Linux MPI PSBLAS Iterative Solvers Preconditioners) /Creator (pdfLaTeX) /Producer ($Id$) /Author()/Title()/Subject()/Creator(LaTeX with hyperref package)/Producer(pdfTeX-1.40.19)/Keywords() -/CreationDate (D:20191218141557Z) -/ModDate (D:20191218141557Z) +/CreationDate (D:20200425083929+02'00') +/ModDate (D:20200425083929+02'00') /Trapped /False /PTEX.Fullbanner (This is pdfTeX, Version 3.14159265-2.6-1.40.19 (TeX Live 2018) kpathsea version 6.3.0) >> @@ -30279,7 +30273,7 @@ endobj /W [1 3 1] /Root 2130 0 R /Info 2131 0 R -/ID [<6B1E8680CFEF4C06C14E8290AA58B896> <6B1E8680CFEF4C06C14E8290AA58B896>] +/ID [ ] /Length 10665 >> stream @@ -30294,34 +30288,35 @@ stream     !"#$%&'().,_5*+/0123456789:;<=>?@ABCDEFGHIJKLMNOPQRSTUVWXYZ[\`^-]abc     !#$&(-+_* "%')/13456789:;=,W<.02?@ABCDEFGHJ>h]>]5]+]?]=]2]3]<]/]0PQR WUK -TLMNO]1XY^VcZ \]];c@S@_5@@`ab].A]-@@@ @*@@@@ -@z@  @ @@@@@@@@)@'@b@@@@@@@@@ @!@"@#@$@%@&@*@4@(@+@-@.@/@0@1@2@6@;@5Y@7@:@=@>@?@F@<@@@A@B@C@D@E@H@I@T@Gڌ@J@K@L@M@N@O@P@Q@R@S@V@W@]@U@X@c@^@_9<b L d.]:9$ "*%&/+,401:5V6>;o<@AH?$BGKOI#JQRUP2MS W_]VAX\ac^S`G  -  L (̦!"#$%&'*+,-./3)0264589:>7;=G?@ABCDEF!IJKOH1LNQRSWPJTV`XgYZ[\]^_bcffavjffffff f -f ffmf f"fffffffffffff fWfff"f#f'f!f$f&f)f*f.f(]f+f-f0f1f2f3f<f:f/3Of4f6f7f8f9f=f>f?f@fAfBfCfEf;P2fD#fGfMfKfFkfHfJfNfOfPfQfRfSfUfL|fTfWfXfZfV'fYf\f]f^f_f`fafbf[fc  -  $%+ #Is!%&'6.$L;(*+,-u7:/h89012345`=;<%?@AH>BDEFGJ9KLQIMNOP]TR)S^U8\VWXYZ[H/{`abI_6cIIII -:4III&I -III  vI I IIIIII $8III 4III I NIII"I%I! hI#I$I-I& FI'I(I)I*I+I,'I/I2I. I0I4I9I3 1I5I6I7I8I;I<I@I: I=I?IBICIFIA NIDIHIKIG bIIIMINIQIL lIO(IVIR - ISITIUIXIYI_IW -EIZI\I]I^IcI` -#Ia +TLMNO]1XY^VcZ \]];c@S@_5@@`ab].]-@@@ @*@@@@ +@z@  @ @@@@@@@@)@'@b@@@@@@@@@ @!@"@#@$@%@&@*@4@(e@+@-@.@/@0@1@2@6@;@5@7@:@=@>@?@F@<@@@A@B@C@D@E@H@I@T@G@J@K@L@M@N@O@P@Q@R@S@V@W@]@U{@X@c@^+@_9b MM d]:$9 "*%L&/+-,40ρ1:56>;<@AH?BGKOI$JQRUP2S W_]VBX\ac^Tm`k  +E %  (!"#$%&'*+,-./3)a0264589:>75;=G? M@ABCDEF!IJKOH2OLNQRSWPK<TV`XgyYZ[\]^_bcffavfftff8ff f +f fff f"ffpfffffffffff ffff"f#f'f! f$f&f)f*f.f(f+f-f0f1f2f3f<f:f/3f4f6f7f8f9f=f>f?f@fAfBfCfEf;PfD#fGfMfKfFl9fHfJfNfOfPfQfRfSfUfLfTfWfXfZfVfYf\f]f^f_f`fafbf[#fcΌg  +   $&Z,d #I!%&'6.$L(*+,-v7:/iV89012345j=;q<%?@AH>BDEFGJ: +KLQIMNOP]TR*(S^U8\VWXYZ[H|`abI_cIIII +:II I&I +III  I I IIIIII $III 5CIII I OIII"I%I! iAI#I$I-I& I'I(I)I*I+I,'I/I2I. `I0I4I9I3 I5I6I7I8I;I<I@I: ;I=I?IBICIFIA IDIHIKIG IIIMINIQIL IO(IVIR + +NISITIUIXIYI_IW +IZI\I]I^IcI` +$Ia  -c  s  -   -)" - !$%(# -#&*+/) -y,.1250 -G378;6 -9@< -=>?*BCDGA EIJMH KQN %8OSVR 2T[W LXYZ]`\ O^+b!a fqc! 1!!! !!! -!!  >! ! !!! :!!!!! !!!!!"! )!! !!,!$!'!# !%!-!( !)!*!+!,!/!0!1!4!. l!2!8!5 %!6!D!9 <!:!;!<!=!>!?!@!A!B!C!H!E X*!F-!N!I Y!J!L!M!T!O o*!P!R!S![!U !V!X!Y!Z!_!\ !]!c!` f!a } . a  ܅      z!  /'" 1s#%&-( FX)+,4. [ /12 395 s68>: ;=@AD? 5B0`abE _FGHIJKLMNOPQRSTUVWXYZ[\]^_c  -w   4 Z1! #("/$&'*/)Dj+-.50R1347<6a8:;B=r]>@A2FCDHIMGJLOPQRUNSXVlWZ[\]`Y^b]ac3q]]]]]] ] -] ]] ] ]]]]];]]]]] ]]]] ]8]]"]$]!U ]#4]']%g]&]7](h+])]*],]4]6]A]8]9]@]B]C]D]E]F]G]H]I]J]K]L]M]N]O]P]Q]R]S]T]U?]V]W7]X^]Y|B]Z][]\C]]!]^]_]`f]a -]bI]cSש'lVn56789:;5  -    !"#$%&'()*+,-./0123456789:;<=>?@ABCDEFGHIJKLMNOPQRSTUVWXYZ[\]^_`abc>v>>>>>>>>> > -> > > >>>>>>>A" +d*  t  +9   +X)" + !$%(# +&*+/) +,.1250 +378;6 +9@< +=>?*BCDGA EIJMH 7KQN %OSVR 3xT[W MMXYZ]`\ PK^+b!a fc! !!! !!! +!!  ! ! !!! !!!!! _!!!!!"! !! !!,!$!'!# 1!%!-!( 2!)!*!+!,!/!0!1!4!. !2!8!5 &L!6!D!9 =!:!;!<!=!>!?!@!A!B!C!H!E X!F-!N!I Y!J!L!M!T!O o!P!R!S![!U !V!X!Y!Z!_!\ e!]!c!` !a  .         ! | /'" 1#%&-( F)+,4. [/12 395 tH68>: r;=@AD? B0`abE FGHIJKLMNOPQRSTUVWXYZ[\]^_c <x +    W1!) #("0$&'*/)D+-.50Sr1347<6b48:;B=r>@A2FCTDHIMGPJLOPQRUNnSXVWZ[\]`Yh^b]aӋc3q]]]]]] ] *] ]] ] ]]]]]]]]]] ]]]] ]9u]]"]$]!U]#4]']%h]&]7](h])]*],]4]6]A]8G]9]@]B]C]D]E]F]G]H]I]J]K]L]M]N]O]P]Q]R]S]T]U]V]W7}]X_]Y|]ZK][I]\]]"!]^]_]`]a +]bJN]cR"'W\m56789:;6  +    !"#$%&'()*+,-./0123456789:;<=>?@ABCDEFGHIJKLMNOPQRSTUVWXYZ[\]^_`abc>>>>>>>>>> > +> > > >>>>>>>F3 endstream endobj startxref -1348130 +1348403 %%EOF diff --git a/docs/src/Makefile b/docs/src/Makefile index 78c2bc85..4e917009 100644 --- a/docs/src/Makefile +++ b/docs/src/Makefile @@ -112,7 +112,7 @@ DATE = $(shell date +%Y-%m-%d) ## Programs # LATEX = latex -LTX2HTML = latex2html -local_icons +#LTX2HTML = latex2html -local_icons LTX2HTML = htlatex PDFLATEX = pdflatex ACRO = evince diff --git a/docs/src/intro.tex b/docs/src/intro.tex index 0ee391a5..92ce6040 100644 --- a/docs/src/intro.tex +++ b/docs/src/intro.tex @@ -293,7 +293,7 @@ are classified as: To finish our general description, we define a version string with the constant \[ \verb|psb_version_string_|\] -whose current value is \verb|3.4.0| +whose current value is \verb|3.7.0| \subsection{Application structure} \label{sec:appstruct} @@ -360,7 +360,8 @@ follows: build it with \verb|prec%build|\footnote{The subroutine style {\tt psb\_precinit} and {\tt psb\_precbl} are still supported for backward compatibility}. -\item Call the iterative method of choice, e.g. \verb|psb_bicgstab| +\item Call the iterative driver \verb|psb_krylov| with the method of + choice, e.g. \verb|bicgstab|. \end{enumerate} This is the structure of the sample programs in the directory \verb|test/pargen/|. diff --git a/prec/Makefile b/prec/Makefile index 0b9079c5..5b9551f2 100644 --- a/prec/Makefile +++ b/prec/Makefile @@ -27,7 +27,7 @@ lib: $(OBJS) impld /bin/cp -p $(CPUPDFLAG) *$(.mod) $(MODDIR) impld: $(OBJS) - cd impl && $(MAKE) + $(MAKE) -C impl $(OBJS): $(MODDIR)/$(BASEMODNAME)$(.mod)