diff --git a/src/tools/Makefile b/src/tools/Makefile index 9fdba620..ea54dad7 100644 --- a/src/tools/Makefile +++ b/src/tools/Makefile @@ -2,8 +2,8 @@ include ../../Make.inc FOBJS = psb_dallc.o psb_dasb.o psb_dcsrp.o psb_cdprt.o \ psb_dfree.o psb_dgelp.o psb_dins.o \ - psb_cdall.o psb_cdalv.o psb_cdasb.o psb_cdcpy.o \ - psb_cddec.o psb_cdfree.o psb_cdins.o psb_dcdovr.o \ + psb_cdall.o psb_cdalv.o psb_cd_inloc.o psb_cdcpy.o \ + psb_cddec.o psb_cdfree.o psb_cdins.o \ psb_cdren.o psb_cdrep.o psb_cdtransfer.o psb_get_overlap.o\ psb_dspalloc.o psb_dspasb.o \ psb_dspcnv.o psb_dspfree.o psb_dspins.o psb_dsprn.o \ @@ -13,7 +13,8 @@ FOBJS = psb_dallc.o psb_dasb.o psb_dcsrp.o psb_cdprt.o \ psb_zspalloc.o psb_zspasb.o psb_zspcnv.o psb_zspfree.o\ psb_zspins.o psb_zsprn.o psb_zcdovr.o psb_zgelp.o -MPFOBJS = psb_dcdovrbld.o psb_dsphalo.o psb_zcdovrbld.o psb_zsphalo.o +MPFOBJS = psb_dcdovrbld.o psb_dsphalo.o psb_zcdovrbld.o psb_zsphalo.o psb_cdasb.o \ + psb_dcdovr.o psb_zcdovr.o INCDIRS = -I ../../lib -I . LIBDIR = ../../lib @@ -26,7 +27,8 @@ lib: mpfobjs $(FOBJS) mpfobjs: (make $(MPFOBJS) F90="$(MPF90)" FC="$(MPF90)" FCOPT="$(F90COPT)") - +psb_glob_to_loc.o: psb_glob_to_loc.f90 + $(F90) $(F90COPT) $(INCDIRS) -c $< $(F90INLINE) clean: /bin/rm -f $(MPFOBJS) $(FOBJS) diff --git a/src/tools/psb_cdall.f90 b/src/tools/psb_cdall.f90 index 5e7ac988..deeb6d8e 100644 --- a/src/tools/psb_cdall.f90 +++ b/src/tools/psb_cdall.f90 @@ -269,6 +269,7 @@ subroutine psb_cdall(m, n, parts, ictxt, desc_a, info) endif ! estimate local cols number loc_col=min(2*loc_row,m) + allocate(desc_a%loc_to_glob(loc_col),& &desc_a%lprm(1),stat=info) if (info /= 0) then diff --git a/src/tools/psb_cdalv.f90 b/src/tools/psb_cdalv.f90 index f41113b6..67a4e338 100644 --- a/src/tools/psb_cdalv.f90 +++ b/src/tools/psb_cdalv.f90 @@ -56,7 +56,7 @@ subroutine psb_cdalv(m, v, ictxt, desc_a, info, flag) !locals Integer :: counter,i,j,np,me,loc_row,err,& - & loc_col,nprocs,n,itmpov, k,& + & loc_col,nprocs,n,itmpov, k,glx,gidx,gle,& & l_ov_ix,l_ov_el,idx, flag_, err_act integer :: int_err(5),exch(2) Integer, allocatable :: temp_ovrlap(:), ov_idx(:),ov_el(:) @@ -88,8 +88,8 @@ subroutine psb_cdalv(m, v, ictxt, desc_a, info, flag) endif if (info /= 0) then - call psb_errpush(info,name,i_err=int_err) - goto 9999 + call psb_errpush(info,name,i_err=int_err) + goto 9999 end if if (debug) write(*,*) 'psb_cdall: doing global checks' @@ -122,7 +122,7 @@ subroutine psb_cdalv(m, v, ictxt, desc_a, info, flag) else flag_=0 endif - + if ((flag_<0).or.(flag_>1)) then info = 6 err=info @@ -132,8 +132,13 @@ subroutine psb_cdalv(m, v, ictxt, desc_a, info, flag) !count local rows number ! allocate work vector - allocate(desc_a%glob_to_loc(m),desc_a%matrix_data(psb_mdata_size_),& - &temp_ovrlap(m),stat=info) + if (m > psb_cd_get_large_threshold()) then + allocate(desc_a%matrix_data(psb_mdata_size_),& + &temp_ovrlap(m),stat=info) + else + allocate(desc_a%glob_to_loc(m),desc_a%matrix_data(psb_mdata_size_),& + &temp_ovrlap(m),stat=info) + end if if (info /= 0) then info=2025 int_err(1)=m @@ -146,34 +151,127 @@ subroutine psb_cdalv(m, v, ictxt, desc_a, info, flag) counter = 0 itmpov = 0 temp_ovrlap(:) = -1 - do i=1,m - - if (((v(i)-flag_) > np-1).or.((v(i)-flag_) < 0)) then - info=580 - int_err(1)=3 - int_err(2)=v(i) - flag_ - int_err(3)=i - exit + + if ( m >psb_cd_get_large_threshold()) then + desc_a%matrix_data(psb_dec_type_) = psb_desc_large_bld_ + + do i=1,m + + if (((v(i)-flag_) > np-1).or.((v(i)-flag_) < 0)) then + info=580 + int_err(1)=3 + int_err(2)=v(i) - flag_ + int_err(3)=i + exit + end if + + if ((v(i)-flag_) == me) then + ! this point belongs to me + counter=counter+1 + end if + enddo + + + loc_row=counter + ! check on parts function + if (debug) write(*,*) 'PSB_CDALL: End main loop:' ,loc_row,itmpov,info + + if (info /= 0) then + call psb_errpush(info,name,i_err=int_err) + goto 9999 end if + if (debug) write(*,*) 'PSB_CDALL: error check:' ,err - if ((v(i)-flag_) == me) then - ! this point belongs to me - counter=counter+1 - desc_a%glob_to_loc(i) = counter - else - desc_a%glob_to_loc(i) = -(np+(v(i)-flag_)+1) + ! estimate local cols number + loc_col = min(2*loc_row,m) + + allocate(desc_a%loc_to_glob(loc_col), desc_a%lprm(1),& + & desc_a%ptree(2),stat=info) + if (info == 0) call InitPairSearchTree(desc_a%ptree,info) + if (info /= 0) then + info=2025 + int_err(1)=loc_col + call psb_errpush(info,name,i_err=int_err) + goto 9999 end if - enddo - loc_row=counter - ! check on parts function - if (debug) write(*,*) 'PSB_CDALL: End main loop:' ,loc_row,itmpov,info + ! set LOC_TO_GLOB array to all "-1" values + desc_a%lprm(1) = 0 + desc_a%loc_to_glob(:) = -1 + k = 0 + do i=1,m + if ((v(i)-flag_) == me) then + k = k + 1 + desc_a%loc_to_glob(k) = i + call SearchInsKeyVal(desc_a%ptree,i,k,glx,info) + endif + enddo + if (k /= loc_row) then + write(0,*) 'Large cd init: ',k,loc_row + end if + + if (info /= 0) then + info=4000 + call psb_errpush(info,name) + goto 9999 + endif + + else + + desc_a%matrix_data(psb_dec_type_) = psb_desc_bld_ + do i=1,m + + if (((v(i)-flag_) > np-1).or.((v(i)-flag_) < 0)) then + info=580 + int_err(1)=3 + int_err(2)=v(i) - flag_ + int_err(3)=i + exit + end if + + if ((v(i)-flag_) == me) then + ! this point belongs to me + counter=counter+1 + desc_a%glob_to_loc(i) = counter + else + desc_a%glob_to_loc(i) = -(np+(v(i)-flag_)+1) + end if + enddo + + + loc_row=counter + ! check on parts function + if (debug) write(*,*) 'PSB_CDALL: End main loop:' ,loc_row,itmpov,info + + if (info /= 0) then + call psb_errpush(info,name,i_err=int_err) + goto 9999 + end if + if (debug) write(*,*) 'PSB_CDALL: error check:' ,err + + ! estimate local cols number + loc_col = min(2*loc_row,m) + + allocate(desc_a%loc_to_glob(loc_col),& + &desc_a%lprm(1),stat=info) + if (info /= 0) then + info=2025 + int_err(1)=loc_col + call psb_errpush(info,name,i_err=int_err) + goto 9999 + end if + + ! set LOC_TO_GLOB array to all "-1" values + desc_a%lprm(1) = 0 + desc_a%loc_to_glob(:) = -1 + do i=1,m + k = desc_a%glob_to_loc(i) + if (k.gt.0) then + desc_a%loc_to_glob(k) = i + endif + enddo - if (info /= 0) then - call psb_errpush(info,name,i_err=int_err) - goto 9999 end if - if (debug) write(*,*) 'PSB_CDALL: error check:' ,err l_ov_ix=0 l_ov_el=0 @@ -234,27 +332,6 @@ subroutine psb_cdalv(m, v, ictxt, desc_a, info, flag) goto 9999 endif - ! estimate local cols number - loc_col = min(2*loc_row,m) - allocate(desc_a%loc_to_glob(loc_col),& - &desc_a%lprm(1),stat=info) - if (info /= 0) then - info=2025 - int_err(1)=loc_col - call psb_errpush(info,name,i_err=int_err) - goto 9999 - end if - - ! set LOC_TO_GLOB array to all "-1" values - desc_a%lprm(1) = 0 - desc_a%loc_to_glob(:) = -1 - do i=1,m - k = desc_a%glob_to_loc(i) - if (k.gt.0) then - desc_a%loc_to_glob(k) = i - endif - enddo - ! set fields in desc_a%MATRIX_DATA.... desc_a%matrix_data(psb_n_row_) = loc_row desc_a%matrix_data(psb_n_col_) = loc_row @@ -268,22 +345,20 @@ subroutine psb_cdalv(m, v, ictxt, desc_a, info, flag) desc_a%halo_index(:) = -1 - - desc_a%matrix_data(psb_m_) = m desc_a%matrix_data(psb_n_) = n - desc_a%matrix_data(psb_dec_type_) = psb_desc_bld_ desc_a%matrix_data(psb_ctxt_) = ictxt call psb_get_mpicomm(ictxt,desc_a%matrix_data(psb_mpi_c_)) + call psb_erractionrestore(err_act) return - + 9999 continue call psb_erractionrestore(err_act) if (err_act == act_abort) then - call psb_error(ictxt) - return + call psb_error(ictxt) + return end if return diff --git a/src/tools/psb_cdasb.f90 b/src/tools/psb_cdasb.f90 index f93b03c7..762e3c42 100644 --- a/src/tools/psb_cdasb.f90 +++ b/src/tools/psb_cdasb.f90 @@ -1,4 +1,4 @@ -!!$ +!!$ !!$ Parallel Sparse BLAS v2.0 !!$ (C) Copyright 2006 Salvatore Filippone University of Rome Tor Vergata !!$ Alfredo Buttari University of Rome Tor Vergata @@ -45,6 +45,7 @@ subroutine psb_cdasb(desc_a,info) use psb_penv_mod implicit none + include 'mpif.h' !...Parameters.... type(psb_desc_type), intent(inout) :: desc_a integer, intent(out) :: info @@ -52,10 +53,11 @@ subroutine psb_cdasb(desc_a,info) !....Locals.... integer :: int_err(5), itemp(2) - integer,allocatable :: ovrlap_index(:),halo_index(:) - integer :: i,err,np,me,& - & lovrlap,lhalo,nhalo,novrlap,max_size,max_halo,n_col,ldesc_halo,& - & ldesc_ovrlap, dectype, err_act + integer,allocatable :: ovrlap_index(:),halo_index(:) + + integer :: i,j,err,np,me,lovrlap,lhalo,nhalo,novrlap,max_size,& + & max_halo,n_col,ldesc_halo, ldesc_ovrlap, dectype, err_act, & + & key, ih, nh, idx, nk,icomm,hsize integer :: ictxt,n_row logical, parameter :: debug=.false., debugwrt=.false. character(len=20) :: name,ch_err @@ -70,6 +72,7 @@ subroutine psb_cdasb(desc_a,info) dectype = psb_cd_get_dectype(desc_a) n_row = psb_cd_get_local_rows(desc_a) n_col = psb_cd_get_local_cols(desc_a) + call psb_get_mpicomm(ictxt,icomm ) ! check on blacs grid call psb_info(ictxt, me, np) @@ -102,8 +105,13 @@ subroutine psb_cdasb(desc_a,info) call psb_errpush(info,name,i_err=int_err) goto 9999 endif + call psb_realloc(psb_cd_get_local_cols(desc_a),desc_a%loc_to_glob,info) + if (psb_is_large_desc(desc_a)) then + call psi_ldsc_pre_halo(desc_a,info) + end if + call psb_transfer(desc_a%ovrlap_index,ovrlap_index,info) call psb_transfer(desc_a%halo_index,halo_index,info) @@ -113,8 +121,6 @@ subroutine psb_cdasb(desc_a,info) goto 9999 end if - ! Ok, register into MATRIX_DATA & free temporary work areas - desc_a%matrix_data(psb_dec_type_) = psb_desc_asb_ deallocate(ovrlap_index, halo_index, stat=info) if (info /= 0) then @@ -122,7 +128,15 @@ subroutine psb_cdasb(desc_a,info) call psb_errpush(info,name) goto 9999 end if - + + if (psb_is_large_dec(dectype) ) then + desc_a%matrix_data(psb_dec_type_) = psb_desc_large_asb_ +!!$ write(0,*) 'Done large dec asmbly',desc_a%matrix_data(psb_dec_type_),& +!!$ & psb_desc_large_asb_,psb_is_asb_dec(desc_a%matrix_data(psb_dec_type_)) + else + ! Ok, register into MATRIX_DATA & free temporary work areas + desc_a%matrix_data(psb_dec_type_) = psb_desc_asb_ + endif else info = 600 call psb_errpush(info,name) diff --git a/src/tools/psb_cdcpy.f90 b/src/tools/psb_cdcpy.f90 index 2ac5980b..ae27c94d 100644 --- a/src/tools/psb_cdcpy.f90 +++ b/src/tools/psb_cdcpy.f90 @@ -54,7 +54,7 @@ subroutine psb_cdcpy(desc_in, desc_out, info) integer, intent(out) :: info !locals - integer :: np,me,ictxt, isz, err_act + integer :: np,me,ictxt, isz, err_act,idx,gidx,lidx logical, parameter :: debug=.false.,debugprt=.false. character(len=20) :: name, char_err if (debug) write(0,*) me,'Entered CDCPY' @@ -63,37 +63,60 @@ subroutine psb_cdcpy(desc_in, desc_out, info) call psb_erractionsave(err_act) name = 'psb_cdcpy' - - ictxt=desc_in%matrix_data(psb_ctxt_) + ictxt = psb_cd_get_context(desc_in) ! check on blacs grid call psb_info(ictxt, me, np) if (debug) write(0,*) me,'Entered CDCPY' if (np == -1) then - info = 2010 - call psb_errpush(info,name) - goto 9999 + info = 2010 + call psb_errpush(info,name) + goto 9999 endif +!!$ call psb_cdfree(desc_out,info) +!!$ call psb_nullify_desc(desc_out) call psb_safe_cpy(desc_in%matrix_data,desc_out%matrix_data,info) if (info == 0) call psb_safe_cpy(desc_in%halo_index,desc_out%halo_index,info) -!!$ if (info == 0) call psb_safe_cpy(desc_in%halo_pt,desc_out%halo_pt,info) if (info == 0) call psb_safe_cpy(desc_in%ovrlap_index,desc_out%ovrlap_index,info) -!!$ if (info == 0) call psb_safe_cpy(desc_in%ovrlap_pt,desc_out%ovrlap_pt,info) if (info == 0) call psb_safe_cpy(desc_in%bnd_elem,desc_out%bnd_elem,info) if (info == 0) call psb_safe_cpy(desc_in%ovrlap_elem,desc_out%ovrlap_elem,info) if (info == 0) call psb_safe_cpy(desc_in%loc_to_glob,desc_out%loc_to_glob,info) if (info == 0) call psb_safe_cpy(desc_in%glob_to_loc,desc_out%glob_to_loc,info) if (info == 0) call psb_safe_cpy(desc_in%lprm,desc_out%lprm,info) if (info == 0) call psb_safe_cpy(desc_in%idx_space,desc_out%idx_space,info) - + if (info == 0) call psb_safe_cpy(desc_in%hashv,desc_out%hashv,info) + if (info == 0) call psb_safe_cpy(desc_in%glb_lc,desc_out%glb_lc,info) + + if (info == 0) then + if (allocated(desc_in%ptree)) then + allocate(desc_out%ptree(2),stat=info) + if (info /= 0) then + info=4000 + goto 9999 + endif + if (.true.) then + call ClonePairSearchTree(desc_in%ptree,desc_out%ptree) + else + call InitPairSearchTree(desc_out%ptree,info) + do idx=1, psb_cd_get_local_cols(desc_out) + gidx = desc_out%loc_to_glob(idx) + call SearchInsKeyVal(desc_out%ptree,gidx,idx,lidx,info) + if (lidx /= idx) then + write(0,*) 'Warning from cdcpy: mismatch in PTREE ',idx,lidx + endif + enddo + end if + end if + end if + if (info /= 0) then info = 4010 call psb_errpush(info,name) goto 9999 endif - + call psb_erractionrestore(err_act) return diff --git a/src/tools/psb_cddec.f90 b/src/tools/psb_cddec.f90 index 27ec0622..ddfcd950 100644 --- a/src/tools/psb_cddec.f90 +++ b/src/tools/psb_cddec.f90 @@ -203,7 +203,6 @@ subroutine psb_cddec(nloc, ictxt, desc_a, info) endif enddo - tovr = -1 thalo = -1 @@ -215,7 +214,6 @@ subroutine psb_cddec(nloc, ictxt, desc_a, info) goto 9999 end if - desc_a%matrix_data(psb_dec_type_) = psb_desc_asb_ call psb_erractionrestore(err_act) diff --git a/src/tools/psb_cdfree.f90 b/src/tools/psb_cdfree.f90 index 113d6383..dbb13bbd 100644 --- a/src/tools/psb_cdfree.f90 +++ b/src/tools/psb_cdfree.f90 @@ -63,7 +63,7 @@ subroutine psb_cdfree(desc_a,info) end if ictxt=psb_cd_get_context(desc_a) - deallocate(desc_a%matrix_data) + call psb_info(ictxt, me, np) ! ....verify blacs grid correctness.. if (np == -1) then @@ -74,7 +74,7 @@ subroutine psb_cdfree(desc_a,info) !...deallocate desc_a.... if(.not.allocated(desc_a%loc_to_glob)) then - info=295 + info=296 call psb_errpush(info,name) goto 9999 end if @@ -87,22 +87,24 @@ subroutine psb_cdfree(desc_a,info) goto 9999 end if - if (.not.allocated(desc_a%glob_to_loc)) then - info=295 - call psb_errpush(info,name) - goto 9999 - end if - - !deallocate glob_to_loc field - deallocate(desc_a%glob_to_loc,stat=info) - if (info /= 0) then - info=2052 - call psb_errpush(info,name) - goto 9999 - end if + if (.not.psb_is_large_desc(desc_a)) then + if (.not.allocated(desc_a%glob_to_loc)) then + info=297 + call psb_errpush(info,name) + goto 9999 + end if + + !deallocate glob_to_loc field + deallocate(desc_a%glob_to_loc,stat=info) + if (info /= 0) then + info=2052 + call psb_errpush(info,name) + goto 9999 + end if + endif if (.not.allocated(desc_a%halo_index)) then - info=295 + info=298 call psb_errpush(info,name) goto 9999 end if @@ -131,7 +133,7 @@ subroutine psb_cdfree(desc_a,info) end if if (.not.allocated(desc_a%ovrlap_index)) then - info=295 + info=299 call psb_errpush(info,name) goto 9999 end if @@ -160,6 +162,36 @@ subroutine psb_cdfree(desc_a,info) goto 9999 end if + if (allocated(desc_a%hashv)) then + deallocate(desc_a%hashv,stat=info) + if (info /= 0) then + info=2058 + call psb_errpush(info,name) + goto 9999 + end if + end if + + if (allocated(desc_a%glb_lc)) then + deallocate(desc_a%glb_lc,stat=info) + if (info /= 0) then + info=2059 + call psb_errpush(info,name) + goto 9999 + end if + end if + + if (allocated(desc_a%ptree)) then + call FreePairSearchTree(desc_a%ptree) + deallocate(desc_a%ptree,stat=info) + if (info /= 0) then + info=2059 + call psb_errpush(info,name) + goto 9999 + end if + end if + + + if (allocated(desc_a%idx_space)) then deallocate(desc_a%idx_space,stat=info) if (info /= 0) then @@ -169,6 +201,8 @@ subroutine psb_cdfree(desc_a,info) end if end if + deallocate(desc_a%matrix_data) + call psb_nullify_desc(desc_a) call psb_erractionrestore(err_act) diff --git a/src/tools/psb_cdins.f90 b/src/tools/psb_cdins.f90 index 5585c8cd..398087e9 100644 --- a/src/tools/psb_cdins.f90 +++ b/src/tools/psb_cdins.f90 @@ -39,34 +39,37 @@ ! ja - integer,dimension(:). The column indices of the points. ! desc_a - type(). The communication descriptor to be freed. ! info - integer. Eventually returns an error code. -subroutine psb_cdins(nz,ia,ja,desc_a,info) +subroutine psb_cdins(nz,ia,ja,desc_a,info,ila,jla) use psb_descriptor_type use psb_serial_mod use psb_const_mod use psb_error_mod use psb_penv_mod + use psi_mod implicit none !....PARAMETERS... Type(psb_desc_type), intent(inout) :: desc_a Integer, intent(in) :: nz,ia(:),ja(:) integer, intent(out) :: info + integer, optional, intent(out) :: ila(:), jla(:) !LOCALS..... integer :: i,ictxt,row,k,dectype,mglob, nglob,err integer :: np, me, isize - integer :: pnt_halo,nrow,ncol, nh, ip,jp, err_act + integer :: pnt_halo,nrow,ncol, nh, ip,jp, err_act,lip,ljp,nxt logical, parameter :: debug=.false. integer, parameter :: relocsz=200 + integer, allocatable :: ila_(:), jla_(:) character(len=20) :: name,ch_err info = 0 name = 'psb_cdins' call psb_erractionsave(err_act) - ictxt = psb_cd_get_context(desc_a) + ictxt = psb_cd_get_context(desc_a) dectype = psb_cd_get_dectype(desc_a) mglob = psb_cd_get_global_rows(desc_a) nglob = psb_cd_get_global_cols(desc_a) @@ -97,71 +100,39 @@ subroutine psb_cdins(nz,ia,ja,desc_a,info) call psb_errpush(info,name) goto 9999 end if - - if (.not.allocated(desc_a%halo_index)) then - allocate(desc_a%halo_index(relocsz)) - desc_a%halo_index(:) = -1 - endif - pnt_halo=1 - do while (desc_a%halo_index(pnt_halo) /= -1 ) - pnt_halo = pnt_halo + 1 - end do - isize = size(desc_a%halo_index) - - do i = 1, nz - ip = ia(i) - jp = ja(i) - if ((ip < 1 ).or.(ip>mglob).or.(jp<1).or.(jp>mglob)) then - ! write(0,*) 'wrong input ',i,ip,jp - info = 1133 + if (present(ila)) then + if (size(ila) < nz) then + info = 1111 call psb_errpush(info,name) goto 9999 + end if + end if + if (present(jla)) then + if (size(jla) < nz) then + info = 1111 + call psb_errpush(info,name) + goto 9999 + end if + end if + + if (present(ila).and.present(jla)) then + call psi_idx_cnv(nz,ia,ila,desc_a,info,owned=.true.) + call psi_idx_ins_cnv(nz,ja,jla,desc_a,info,mask=(ila(1:nz)>0)) + + else + if (present(ila).or.present(jla)) then + write(0,*) 'Inconsistent call : ',present(ila),present(jla) endif - if ((1<=desc_a%glob_to_loc(ip)).and.(desc_a%glob_to_loc(ip))<=nrow) then - k = desc_a%glob_to_loc(jp) - if (k.lt.-np) then - k = k + np - k = - k - 1 - ncol = ncol + 1 - desc_a%glob_to_loc(jp) = ncol - isize = size(desc_a%loc_to_glob) - if (ncol > isize) then - nh = ncol + max(nz,relocsz) - call psb_realloc(nh,desc_a%loc_to_glob,info,pad=-1) - if (me==0) then - if (debug) write(0,*) 'done realloc ',nh - end if - if (info /= 0) then - info=4010 - ch_err='psb_realloc' - call psb_errpush(info,name) - goto 9999 - end if - isize = nh - endif - desc_a%loc_to_glob(ncol) = jp - isize = size(desc_a%halo_index) - if ((pnt_halo+3).gt.isize) then - nh = isize + max(nz,relocsz) - call psb_realloc(nh,desc_a%halo_index,info,pad=-1) - if (info /= 0) then - info=4010 - ch_err='psb_realloc' - call psb_errpush(info,name) - goto 9999 - end if - isize = nh - endif - desc_a%halo_index(pnt_halo) = k - desc_a%halo_index(pnt_halo+1) = 1 - desc_a%halo_index(pnt_halo+2) = ncol - pnt_halo = pnt_halo + 3 - endif - else - ! currently we ignore items not belonging to us. - endif - enddo - desc_a%matrix_data(psb_n_col_) = ncol + allocate(ila_(nz),jla_(nz),stat=info) + if (info /= 0) then + info = 4000 + call psb_errpush(info,name) + goto 9999 + end if + call psi_idx_cnv(nz,ia,ila_,desc_a,info,owned=.true.) + call psi_idx_ins_cnv(nz,ja,jla_,desc_a,info,mask=(ila_(1:nz)>0)) + deallocate(ila_,jla_) + end if call psb_erractionrestore(err_act) return diff --git a/src/tools/psb_cdprt.f90 b/src/tools/psb_cdprt.f90 index 2d14eef3..8422f2cb 100644 --- a/src/tools/psb_cdprt.f90 +++ b/src/tools/psb_cdprt.f90 @@ -70,7 +70,7 @@ subroutine psb_cdprt(iout,desc_p,glob,short) & write(iout,*) 'Loc_to_glob ',desc_p%loc_to_glob(1:n_row), ': ',& & desc_p%loc_to_glob(n_row+1:n_col) - if (.not.lshort) write(iout,*) 'glob_to_loc ',desc_p%glob_to_loc(1:m) +!!$ if (.not.lshort) write(iout,*) 'glob_to_loc ',desc_p%glob_to_loc(1:m) write(iout,*) 'Halo_index' counter = 1 Do diff --git a/src/tools/psb_cdren.f90 b/src/tools/psb_cdren.f90 index 6a5af445..cfb0f7a7 100644 --- a/src/tools/psb_cdren.f90 +++ b/src/tools/psb_cdren.f90 @@ -75,10 +75,10 @@ subroutine psb_cdren(trans,iperm,desc_a,info) time(1) = mpi_wtime() - ictxt=psb_cd_get_context(desc_a) - dectype=psb_cd_get_dectype(desc_a) - n_row = psb_cd_get_local_rows(desc_a) - n_col = psb_cd_get_local_cols(desc_a) + ictxt = psb_cd_get_context(desc_a) + dectype = psb_cd_get_dectype(desc_a) + n_row = psb_cd_get_local_rows(desc_a) + n_col = psb_cd_get_local_cols(desc_a) ! check on blacs grid call psb_info(ictxt, me, np) diff --git a/src/tools/psb_cdrep.f90 b/src/tools/psb_cdrep.f90 index 90e2ec06..7ed846b6 100644 --- a/src/tools/psb_cdrep.f90 +++ b/src/tools/psb_cdrep.f90 @@ -186,7 +186,6 @@ subroutine psb_cdrep(m, ictxt, desc_a, info) endif - desc_a%matrix_data(psb_m_) = m desc_a%matrix_data(psb_n_) = n desc_a%matrix_data(psb_n_row_) = m @@ -194,20 +193,19 @@ subroutine psb_cdrep(m, ictxt, desc_a, info) desc_a%matrix_data(psb_dec_type_) = psb_desc_bld_ desc_a%matrix_data(psb_ctxt_) = ictxt call psb_get_mpicomm(ictxt,desc_a%matrix_data(psb_mpi_c_)) + do i=1,m desc_a%glob_to_loc(i) = i desc_a%loc_to_glob(i) = i enddo - - tovr = -1 thalo = -1 desc_a%lprm(:) = 0 call psi_cnv_dsc(thalo,tovr,desc_a,info) if (info /= 0) then - call psb_errpush(4010,name,a_err='psi_bld_cdesc') + call psb_errpush(4010,name,a_err='psi_cvn_dsc') goto 9999 end if diff --git a/src/tools/psb_cdtransfer.f90 b/src/tools/psb_cdtransfer.f90 index 4298abe3..1a4e8d0e 100644 --- a/src/tools/psb_cdtransfer.f90 +++ b/src/tools/psb_cdtransfer.f90 @@ -50,7 +50,7 @@ subroutine psb_cdtransfer(desc_in, desc_out, info) !....parameters... type(psb_desc_type), intent(inout) :: desc_in - type(psb_desc_type), intent(out) :: desc_out + type(psb_desc_type), intent(inout) :: desc_out integer, intent(out) :: info !locals @@ -63,41 +63,59 @@ subroutine psb_cdtransfer(desc_in, desc_out, info) info = 0 call psb_erractionsave(err_act) name = 'psb_cdtransfer' - - ictxt=desc_in%matrix_data(psb_ctxt_) + + ictxt=psb_cd_get_context(desc_in) call psb_info(ictxt, me, np) if (debug) write(0,*) me,'Entered CDTRANSFER' if (np == -1) then - info = 2010 - call psb_errpush(info,name) - goto 9999 + info = 2010 + call psb_errpush(info,name) + goto 9999 endif -!!$ call psb_nullify_desc(desc_out) -!!$ - call psb_transfer( desc_in%matrix_data , desc_out%matrix_data , info) - call psb_transfer( desc_in%halo_index , desc_out%halo_index , info) - call psb_transfer( desc_in%bnd_elem , desc_out%bnd_elem , info) - call psb_transfer( desc_in%ovrlap_elem , desc_out%ovrlap_elem , info) - call psb_transfer( desc_in%ovrlap_index, desc_out%ovrlap_index , info) - call psb_transfer( desc_in%loc_to_glob , desc_out%loc_to_glob , info) - call psb_transfer( desc_in%glob_to_loc , desc_out%glob_to_loc , info) - call psb_transfer( desc_in%lprm , desc_out%lprm , info) - call psb_transfer( desc_in%idx_space , desc_out%idx_space , info) - + if (info == 0) call psb_transfer( desc_in%halo_index , desc_out%halo_index , info) + if (info == 0) call psb_transfer( desc_in%bnd_elem , desc_out%bnd_elem , info) + if (info == 0) call psb_transfer( desc_in%ovrlap_elem , desc_out%ovrlap_elem , info) + if (info == 0) call psb_transfer( desc_in%ovrlap_index, desc_out%ovrlap_index , info) + if (info == 0) call psb_transfer( desc_in%loc_to_glob , desc_out%loc_to_glob , info) + if (info == 0) call psb_transfer( desc_in%glob_to_loc , desc_out%glob_to_loc , info) + if (info == 0) call psb_transfer( desc_in%lprm , desc_out%lprm , info) + if (info == 0) call psb_transfer( desc_in%idx_space , desc_out%idx_space , info) + if (info == 0) call psb_transfer( desc_in%hashv , desc_out%hashv , info) + if (info == 0) call psb_transfer( desc_in%glb_lc , desc_out%glb_lc , info) + ! Why doesn't transfer work below? Dunno..... + if (info == 0) call psb_transfer( desc_in%ptree , desc_out%ptree , info) + if (.false.) then + if (info == 0) then + if (allocated(desc_in%ptree)) then + allocate(desc_out%ptree(2),stat=info) + if (info /= 0) then + info=4000 + goto 9999 + endif + desc_out%ptree(1:2) = desc_in%ptree(1:2) + deallocate(desc_in%ptree,stat=info) + end if + end if + endif + if (info /= 0) then + info = 4010 + call psb_errpush(info,name) + goto 9999 + endif call psb_erractionrestore(err_act) return - + 9999 continue call psb_erractionrestore(err_act) - + if (err_act == act_ret) then - return + return else - call psb_error(ictxt) + call psb_error(ictxt) end if return diff --git a/src/tools/psb_dallc.f90 b/src/tools/psb_dallc.f90 index 77f9519c..7a929ad8 100644 --- a/src/tools/psb_dallc.f90 +++ b/src/tools/psb_dallc.f90 @@ -51,13 +51,13 @@ subroutine psb_dalloc(x, desc_a, info, n) !....parameters... real(kind(1.d0)), allocatable, intent(out) :: x(:,:) type(psb_desc_type), intent(in) :: desc_a - integer :: info + integer,intent(out) :: info integer, optional, intent(in) :: n !locals integer :: np,me,err,n_col,n_row,i,j,err_act - integer :: ictxt,dectype,n_ + integer :: ictxt,n_ integer :: int_err(5), exch(3) character(len=20) :: name, ch_err @@ -77,7 +77,6 @@ subroutine psb_dalloc(x, desc_a, info, n) goto 9999 endif - dectype=psb_cd_get_dectype(desc_a) !... check m and n parameters.... if (.not.psb_is_ok_desc(desc_a)) then info = 3110 @@ -199,11 +198,11 @@ subroutine psb_dallocv(x, desc_a,info,n) !....parameters... real(kind(1.d0)), allocatable, intent(out) :: x(:) type(psb_desc_type), intent(in) :: desc_a - integer :: info + integer,intent(out) :: info integer, optional, intent(in) :: n !locals - integer :: np,me,n_col,n_row,dectype,i,err_act + integer :: np,me,n_col,n_row,i,err_act integer :: ictxt logical, parameter :: debug=.false. character(len=20) :: name, ch_err @@ -223,9 +222,6 @@ subroutine psb_dallocv(x, desc_a,info,n) goto 9999 endif - dectype=psb_cd_get_dectype(desc_a) - if (debug) write(0,*) 'dall: dectype',dectype - if (debug) write(0,*) 'dall: is_ok? dectype',psb_is_ok_desc(desc_a) !... check m and n parameters.... if (.not.psb_is_ok_desc(desc_a)) then info = 3110 @@ -261,6 +257,8 @@ subroutine psb_dallocv(x, desc_a,info,n) do i=1,n_row x(i) = 0.0d0 end do + else + write(0,*) 'Did not allocate anything because of dectype',psb_cd_get_dectype(desc_a) endif call psb_erractionrestore(err_act) diff --git a/src/tools/psb_dasb.f90 b/src/tools/psb_dasb.f90 index 54728287..b1c54bd8 100644 --- a/src/tools/psb_dasb.f90 +++ b/src/tools/psb_dasb.f90 @@ -53,7 +53,7 @@ subroutine psb_dasb(x, desc_a, info) ! local variables integer :: ictxt,np,me,nrow,ncol, err_act - integer :: int_err(5), i1sz, i2sz, dectype + integer :: int_err(5), i1sz, i2sz real(kind(1.d0)),parameter :: one=1 logical, parameter :: debug=.false. character(len=20) :: name, ch_err @@ -69,14 +69,13 @@ subroutine psb_dasb(x, desc_a, info) goto 9999 endif - ictxt=psb_cd_get_context(desc_a) - dectype=psb_cd_get_dectype(desc_a) + ictxt = psb_cd_get_context(desc_a) call psb_info(ictxt, me, np) if (debug) write(*,*) 'asb start: ',np,me,& - &psb_cd_get_dectype(desc_a) + & psb_cd_get_dectype(desc_a) ! ....verify blacs grid correctness.. if (np == -1) then info = 2010 @@ -84,16 +83,16 @@ subroutine psb_dasb(x, desc_a, info) goto 9999 else if (.not.psb_is_asb_desc(desc_a)) then if (debug) write(*,*) 'asb error ',& - &dectype + & psb_cd_get_dectype(desc_a) info = 3110 call psb_errpush(info,name) goto 9999 endif ! check size - ictxt=psb_cd_get_context(desc_a) - nrow=psb_cd_get_local_rows(desc_a) - ncol=psb_cd_get_local_cols(desc_a) + ictxt = psb_cd_get_context(desc_a) + nrow = psb_cd_get_local_rows(desc_a) + ncol = psb_cd_get_local_cols(desc_a) i1sz = size(x,dim=1) i2sz = size(x,dim=2) if (debug) write(*,*) 'asb: ',i1sz,i2sz,nrow,ncol @@ -184,7 +183,7 @@ subroutine psb_dasbv(x, desc_a, info) ! local variables integer :: ictxt,np,me - integer :: int_err(5), i1sz,nrow,ncol, dectype, err_act + integer :: int_err(5), i1sz,nrow,ncol, err_act real(kind(1.d0)),parameter :: one=1 logical, parameter :: debug=.false. character(len=20) :: name,ch_err @@ -193,8 +192,7 @@ subroutine psb_dasbv(x, desc_a, info) int_err(1) = 0 name = 'psb_dasbv' - ictxt=psb_cd_get_context(desc_a) - dectype=psb_cd_get_dectype(desc_a) + ictxt = psb_cd_get_context(desc_a) call psb_info(ictxt, me, np) @@ -209,8 +207,8 @@ subroutine psb_dasbv(x, desc_a, info) goto 9999 endif - nrow=psb_cd_get_local_rows(desc_a) - ncol=psb_cd_get_local_cols(desc_a) + nrow = psb_cd_get_local_rows(desc_a) + ncol = psb_cd_get_local_cols(desc_a) if (debug) write(*,*) name,' sizes: ',nrow,ncol i1sz = size(x) if (debug) write(*,*) 'dasb: sizes ',i1sz,ncol diff --git a/src/tools/psb_dcdovr.f90 b/src/tools/psb_dcdovr.f90 index cef4a852..5196e2d8 100644 --- a/src/tools/psb_dcdovr.f90 +++ b/src/tools/psb_dcdovr.f90 @@ -40,7 +40,8 @@ ! a - type(). The input sparse matrix. ! desc_a - type(). The input communication descriptor. ! norv - integer. The number of overlap levels. -! desc_ov - type(). The auxiliary output communication descriptor. +! desc_ov - type(). The auxiliary output communication +! descriptor. ! info - integer. Eventually returns an error code. ! Subroutine psb_dcdovr(a,desc_a,novr,desc_ov,info) @@ -51,6 +52,10 @@ Subroutine psb_dcdovr(a,desc_a,novr,desc_ov,info) Use psb_prec_mod use psb_error_mod use psb_penv_mod + use psb_tools_mod + use psb_realloc_mod + use psi_mod + use mpi Implicit None ! .. Array Arguments .. @@ -60,40 +65,23 @@ Subroutine psb_dcdovr(a,desc_a,novr,desc_ov,info) Type(psb_desc_type), Intent(inout) :: desc_ov integer, intent(out) :: info - real(kind(1.d0)) :: t1,t2,t3,mpi_wtime - external mpi_wtime integer icomm, err_act - interface psb_cdcpy - subroutine psb_cdcpy(desc_in,desc_out,info) - use psb_descriptor_type - type(psb_desc_type), intent(in) :: desc_in - type(psb_desc_type), intent(out) :: desc_out - integer, intent(out) :: info - end subroutine psb_cdcpy - end interface - - interface psb_cdovrbld - subroutine psb_dcdovrbld(n_ovr,desc_p,desc_a,a,l_tmp_halo,& - & l_tmp_ovr_idx,lworks,lworkr,info) - use psb_prec_type - use psb_spmat_type - type(psb_dspmat_type),intent(in) :: a - type(psb_desc_type),intent(in) :: desc_a - type(psb_desc_type),intent(inout) :: desc_p - integer,intent(in) :: n_ovr - integer, intent(in) :: l_tmp_halo,l_tmp_ovr_idx - integer, intent(inout) :: lworks, lworkr - integer, intent(out) :: info - end subroutine psb_dcdovrbld - end interface - - - ! .. Local Scalars .. Integer :: i, j, k, np, me,m,nnzero,& - & ictxt, lovr, lworks,lworkr, n_col, int_err(5),& + & ictxt, lovr, lworks,lworkr, n_row,n_col, int_err(5),& & index_dim,elem_dim, l_tmp_ovr_idx,l_tmp_halo, nztot,nhalo + Integer :: counter,counter_h, counter_o, counter_e,idx,gidx,proc,n_elem_recv,& + & n_elem_send,tot_recv,tot_elem,& + & counter_t,n_elem,i_ovr,jj,proc_id,isz, mglob, glx, & + & idxr, idxs, lx, iszr, iszs, nxch, nsnd, nrcv,lidx + + type(psb_dspmat_type) :: blk + Integer, allocatable :: tmp_halo(:),tmp_ovr_idx(:) + Integer,allocatable :: halo(:),works(:),workr(:),t_halo_in(:),& + & t_halo_out(:),temp(:),maskr(:) + Integer,allocatable :: brvindx(:),rvsz(:), bsdindx(:),sdsz(:) + Logical,Parameter :: debug=.false. character(len=20) :: name, ch_err @@ -107,10 +95,12 @@ Subroutine psb_dcdovr(a,desc_a,novr,desc_ov,info) If(debug) Write(0,*)'in psb_cdovr',novr - m=psb_cd_get_local_rows(desc_a) - nnzero=Size(a%aspk) - n_col=psb_cd_get_local_cols(desc_a) - nhalo = n_col-m + m = psb_cd_get_local_rows(desc_a) + nnzero = Size(a%aspk) + n_row = psb_cd_get_local_rows(desc_a) + n_col = psb_cd_get_local_cols(desc_a) + nhalo = n_col-m + If(debug) Write(0,*)'IN CDOVR1',novr ,m,nnzero,n_col if (novr<0) then info=10 @@ -120,34 +110,32 @@ Subroutine psb_dcdovr(a,desc_a,novr,desc_ov,info) goto 9999 endif + if (debug) write(0,*) 'Calling desccpy' + call psb_cdcpy(desc_a,desc_ov,info) + if (info /= 0) then + info=4010 + ch_err='psb_cdcpy' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + if (debug) write(0,*) 'From desccpy' if (novr==0) then ! ! Just copy the input. ! - if (debug) write(0,*) 'Calling desccpy' - call psb_cdcpy(desc_a,desc_ov,info) - if (info /= 0) then - info=4010 - ch_err='psb_cdcpy' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - if (debug) write(0,*) 'From desccpy' return endif call psb_get_mpicomm(ictxt,icomm ) - If(debug) then + If(debug)then Write(0,*)'BEGIN cdovr',me,nhalo call psb_barrier(ictxt) endif - t1 = mpi_wtime() - ! ! Ok, since we are only estimating, do it as follows: - ! LOVR= (NNZ/NROW)*N_HALO*N_OVR This assumes that the local average + ! LOVR= (NNZ/NROW)*N_HALO*NOVR This assumes that the local average ! nonzeros per row is the same as the global. ! nztot = psb_sp_get_nnzeros(a) @@ -161,54 +149,552 @@ Subroutine psb_dcdovr(a,desc_a,novr,desc_ov,info) goto 9999 endif If(debug)Write(0,*)'ovr_est done',me,novr,lovr + index_dim = size(desc_a%halo_index) elem_dim = size(desc_a%halo_index) - call psb_realloc(psb_mdata_size_,desc_ov%matrix_data,info) - if (info==0) call psb_realloc(novr*(Max(elem_dim,1)+3),desc_ov%ovrlap_elem,info) - if (info /= 0) then - info=4000 - call psb_errpush(info,name) - goto 9999 - end if - l_tmp_ovr_idx=novr*(3*Max(2*index_dim,1)+1) - l_tmp_halo=novr*(3*Size(desc_a%halo_index)) + l_tmp_ovr_idx = novr*(3*Max(2*index_dim,1)+1) + l_tmp_halo = novr*(3*Size(desc_a%halo_index)) - desc_ov%matrix_data(:) = desc_a%matrix_data(:) - desc_ov%matrix_data(psb_dec_type_) = psb_desc_bld_ - - Allocate(desc_ov%loc_to_glob(Size(desc_a%loc_to_glob)),& - & desc_ov%glob_to_loc(Size(desc_a%glob_to_loc)),stat=info) - if (info /= 0) then - info=4000 - call psb_errpush(info,name) - goto 9999 + if (psb_is_large_desc(desc_a)) then + desc_ov%matrix_data(psb_dec_type_) = psb_desc_large_bld_ + else + desc_ov%matrix_data(psb_dec_type_) = psb_desc_bld_ end if - - desc_ov%loc_to_glob(:) = desc_a%loc_to_glob(:) - desc_ov%glob_to_loc(:) = desc_a%glob_to_loc(:) If(debug) then Write(0,*)'Start cdovrbld',me,lworks,lworkr call psb_barrier(ictxt) endif - ! - ! The real work goes on in here.... - ! - Call psb_cdovrbld(novr,desc_ov,desc_a,a,& - & l_tmp_halo,l_tmp_ovr_idx,lworks,lworkr,info) - if (info /= 0) then - info=4010 - ch_err='psb_cdovrbld' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 + if (.false.) then + ! + ! The real work goes on in here.... + ! + Call psb_cdovrbld(novr,desc_ov,desc_a,a,& + & l_tmp_halo,l_tmp_ovr_idx,lworks,lworkr,info) + + if (info /= 0) then + info=4010 + ch_err='psb_cdovrbld' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + + If(debug) then + Write(0,*)'Done cdovrbld',me,lworks,lworkr + call psb_barrier(ictxt) + endif + + else + + + Allocate(brvindx(np+1),rvsz(np),sdsz(np),bsdindx(np+1),stat=info) + if (info /= 0) then + call psb_errpush(4010,name,a_err='Allocate') + goto 9999 + end if + + Allocate(works(lworks),workr(lworkr),t_halo_in(l_tmp_halo),& + & t_halo_out(l_tmp_halo), temp(lworkr),stat=info) + if (info /= 0) then + call psb_errpush(4010,name,a_err='Allocate') + goto 9999 + end if + + + call psb_sp_all(blk,max(lworks,lworkr),info) + if (info /= 0) then + info=4010 + ch_err='psb_sp_all' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + + blk%fida='COO' + + Allocate(tmp_ovr_idx(l_tmp_ovr_idx),tmp_halo(l_tmp_halo),& + & halo(size(desc_a%halo_index)),stat=info) + if (info /= 0) then + call psb_errpush(4010,name,a_err='Allocate') + goto 9999 + end if + halo(:) = desc_a%halo_index(:) + desc_ov%ovrlap_elem(:) = -1 + tmp_ovr_idx(:) = -1 + tmp_halo(:) = -1 + counter_e = 1 + tot_recv = 0 + counter_h = 1 + counter_o = 1 + + ! Init overlap with desc_a%ovrlap (if any) + counter = 1 + Do While (desc_a%ovrlap_index(counter) /= -1) + proc = desc_a%ovrlap_index(counter+psb_proc_id_) + n_elem_recv = desc_a%ovrlap_index(counter+psb_n_elem_recv_) + n_elem_send = desc_a%ovrlap_index(counter+n_elem_recv+psb_n_elem_send_) + + Do j=0,n_elem_recv-1 + + idx = desc_a%ovrlap_index(counter+psb_elem_recv_+j) + If(idx > Size(desc_ov%loc_to_glob)) then + info=-3 + call psb_errpush(info,name) + goto 9999 + endif + + gidx = desc_ov%loc_to_glob(idx) + + call psb_check_size((counter_o+3),tmp_ovr_idx,info,pad=-1) + if (info /= 0) then + info=4010 + call psb_errpush(info,name,a_err='psb_check_size') + goto 9999 + end if + + tmp_ovr_idx(counter_o)=proc + tmp_ovr_idx(counter_o+1)=1 + tmp_ovr_idx(counter_o+2)=gidx + tmp_ovr_idx(counter_o+3)=-1 + counter_o=counter_o+3 + end Do + counter=counter+n_elem_recv+n_elem_send+2 + end Do + + + + + ! + ! A picture is in order to understand what goes on here. + ! I is the internal part; H is halo, R row, C column. The final + ! matrix with N levels of overlap looks like this + ! + ! I | Hc1 | 0 | 0 | + ! -------|-----|-----|-----| + ! Hr1 | Hd1 | Hc2 | 0 | + ! -------|-----|-----|-----| + ! 0 | Hr2 | Hd2 | Hc2 | + ! -------|-----|-----|-----| + ! 0 | 0 | Hr3 | Hd3 | Hc3 + ! + ! At the start we already have I and Hc1, so we know the row + ! indices that will make up Hr1, and also who owns them. As we + ! actually get those rows, we receive the column indices in Hc2; + ! these define the row indices for Hr2, and so on. When we have + ! reached the desired level HrN, we may ignore HcN. + ! + ! + Do i_ovr = 1, novr + + if (debug) write(0,*) me,'Running on overlap level ',i_ovr,' of ',novr + + ! + ! At this point, halo contains a valid halo corresponding to the + ! matrix enlarged with the elements in the frontier for I_OVR-1. + ! At the start, this is just the halo for A; the rows for indices in + ! the first halo will contain column indices defining the second halo + ! level and so on. + ! + bsdindx(:) = 0 + sdsz(:) = 0 + brvindx(:) = 0 + rvsz(:) = 0 + idxr = 0 + idxs = 0 + counter = 1 + counter_t = 1 + + + Do While (halo(counter) /= -1) + tot_elem=0 + proc=halo(counter+psb_proc_id_) + n_elem_recv=halo(counter+psb_n_elem_recv_) + n_elem_send=halo(counter+n_elem_recv+psb_n_elem_send_) + If ((counter+n_elem_recv+n_elem_send) > Size(halo)) then + info = -1 + call psb_errpush(info,name) + goto 9999 + end If + tot_recv=tot_recv+n_elem_recv + if (debug) write(0,*) me,' CDOVRBLD tot_recv:',proc,n_elem_recv,tot_recv + ! + ! + ! The format of the halo vector exists in two forms: 1. Temporary + ! 2. Assembled. In this loop we are using the (assembled) halo_in and + ! copying it into (temporary) halo_out; this is because tmp_halo will + ! be enlarged with the new column indices received, and will reassemble + ! everything for the next iteration. + ! + + ! + ! add recv elements in halo_index into ovrlap_index + ! + Do j=0,n_elem_recv-1 + If((counter+psb_elem_recv_+j)>Size(halo)) then + info=-2 + call psb_errpush(info,name) + goto 9999 + end If + + idx = halo(counter+psb_elem_recv_+j) + If(idx > Size(desc_ov%loc_to_glob)) then + info=-3 + call psb_errpush(info,name) + goto 9999 + endif + + gidx = desc_ov%loc_to_glob(idx) + + call psb_check_size((counter_o+3),tmp_ovr_idx,info,pad=-1) + if (info /= 0) then + info=4010 + call psb_errpush(info,name,a_err='psb_check_size') + goto 9999 + end if + + tmp_ovr_idx(counter_o)=proc + tmp_ovr_idx(counter_o+1)=1 + tmp_ovr_idx(counter_o+2)=gidx + tmp_ovr_idx(counter_o+3)=-1 + counter_o=counter_o+3 + if (.not.psb_is_large_desc(desc_ov)) then + call psb_check_size((counter_h+3),tmp_halo,info,pad=-1) + if (info /= 0) then + info=4010 + call psb_errpush(info,name,a_err='psb_check_size') + goto 9999 + end if + + tmp_halo(counter_h)=proc + tmp_halo(counter_h+1)=1 + tmp_halo(counter_h+2)=idx + tmp_halo(counter_h+3)=-1 + + counter_h=counter_h+3 + end if + + Enddo + if (debug) write(0,*) me,'Checktmp_o_i Loop Mid1',tmp_ovr_idx(1:10) + counter = counter+n_elem_recv + + ! + ! add send elements in halo_index into ovrlap_index + ! + Do j=0,n_elem_send-1 + + idx = halo(counter+psb_elem_send_+j) + gidx = desc_ov%loc_to_glob(idx) + if (idx > psb_cd_get_local_rows(Desc_a)) & + & write(0,*) me,i_ovr,'Out of local rows ',& + & idx,psb_cd_get_local_rows(Desc_a) + + call psb_check_size((counter_o+3),tmp_ovr_idx,info,pad=-1) + if (info /= 0) then + info=4010 + call psb_errpush(info,name,a_err='psb_check_size') + goto 9999 + end if + + tmp_ovr_idx(counter_o)=proc + tmp_ovr_idx(counter_o+1)=1 + tmp_ovr_idx(counter_o+2)=gidx + tmp_ovr_idx(counter_o+3)=-1 + counter_o=counter_o+3 + + ! + ! Prepare to exchange the halo rows with the other proc. + ! + If (i_ovr < (novr)) Then + n_elem = psb_sp_get_nnz_row(idx,a) + + call psb_check_size((idxs+tot_elem+n_elem),works,info) + if (info /= 0) then + info=4010 + call psb_errpush(info,name,a_err='psb_check_size') + goto 9999 + end if + + If((n_elem) > size(blk%ia2)) Then + isz = max((3*size(blk%ia2))/2,(n_elem)) + if (debug) write(0,*) me,'Realloc blk',isz + call psb_sp_reall(blk,isz,info) + if (info /= 0) then + info=4010 + ch_err='psb_sp_reall' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + End If + + call psb_sp_getblk(idx,a,blk,info) + if (info /= 0) then + info=4010 + ch_err='psb_sp_getblk' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if +!!$ write(0,*) me,'Iteration: ',j,i_ovr + Do jj=1,n_elem + works(idxs+tot_elem+jj)=desc_ov%loc_to_glob(blk%ia2(jj)) + End Do + tot_elem=tot_elem+n_elem + End If + + Enddo + + + if (i_ovr < novr) then + if (tot_elem > 1) then + call imsr(tot_elem,works(idxs+1)) + lx = works(idxs+1) + i = 1 + j = 1 + do + j = j + 1 + if (j > tot_elem) exit + if (works(idxs+j) /= lx) then + i = i + 1 + works(idxs+i) = works(idxs+j) + lx = works(idxs+i) + end if + end do + tot_elem = i + endif + if (debug) write(0,*) me,'Checktmp_o_i Loop Mid2',tmp_ovr_idx(1:10) + sdsz(proc+1) = tot_elem + idxs = idxs + tot_elem + end if + counter = counter+n_elem_send+3 + if (debug) write(0,*) me,'Checktmp_o_i Loop End',tmp_ovr_idx(1:10) + Enddo + if (debug) write(0,*)me,'End phase 1 CDOVRBLD', m, n_col, tot_recv + + if (i_ovr < novr) then + ! + ! Exchange data requests with everybody else: so far we have + ! accumulated RECV requests, we have an all-to-all to build + ! matchings SENDs. + ! + call mpi_alltoall(sdsz,1,mpi_integer,rvsz,1,mpi_integer,icomm,info) + if (info /= 0) then + info=4010 + ch_err='mpi_alltoall' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + idxs = 0 + idxr = 0 + counter = 1 + Do + proc=halo(counter) + if (proc == -1) exit + n_elem_recv = halo(counter+psb_n_elem_recv_) + counter = counter+n_elem_recv + n_elem_send = halo(counter+psb_n_elem_send_) + + bsdindx(proc+1) = idxs + idxs = idxs + sdsz(proc+1) + brvindx(proc+1) = idxr + idxr = idxr + rvsz(proc+1) + counter = counter+n_elem_send+3 + Enddo + + iszr=sum(rvsz) + if (max(iszr,1) > lworkr) then + call psb_realloc(max(iszr,1),workr,info) + if (info /= 0) then + info=4010 + ch_err='psb_realloc' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + lworkr=max(iszr,1) + end if + + call mpi_alltoallv(works,sdsz,bsdindx,mpi_integer,& + & workr,rvsz,brvindx,mpi_integer,icomm,info) + if (info /= 0) then + info=4010 + ch_err='mpi_alltoallv' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + + if (debug) write(0,*) 'ISZR :',iszr + + if (psb_is_large_desc(desc_a)) then + call psb_check_size(iszr,maskr,info) + if (info /= 0) then + info=4010 + call psb_errpush(info,name,a_err='psb_check_size') + goto 9999 + end if + call psi_idx_cnv(iszr,workr,maskr,desc_ov,info) + iszs = count(maskr(1:iszr)<=0) + if (iszs > size(works)) call psb_realloc(iszs,works,info) + j = 0 + do i=1,iszr + if (maskr(i) < 0) then + j=j+1 + works(j) = workr(i) + end if + end do + ! + ! fnd_owner on desc_a because we want the procs who + ! owned the rows from the beginning! + ! + call psi_fnd_owner(iszs,works,temp,desc_a,info) + n_col=psb_cd_get_local_cols(desc_ov) + + do i=1,iszs + idx = works(i) + n_col = psb_cd_get_local_cols(desc_ov) + call psi_idx_ins_cnv(idx,lidx,desc_ov,info) + if (psb_cd_get_local_cols(desc_ov) > n_col ) then + ! + ! This is a new index. Assigning a local index as + ! we receive them guarantees that all indices for HALO(I) + ! will be less than those for HALO(J) whenever I 0 thanks to the caller routine. ! ! Parameters: -! n_ovr - integer. The number of overlap levels +! n_ovr - integer. The number of overlap levels ! desc_p - type(). The communication descriptor ! for the preconditioner. ! desc_a - type(). The communication descriptor. @@ -54,34 +54,30 @@ Subroutine psb_dcdovrbld(n_ovr,desc_p,desc_a,a,& use psb_serial_mod Use psi_mod use psb_realloc_mod - use psb_tools_mod, only : psb_cdprt use psb_error_mod use psb_const_mod use psb_penv_mod + use psb_tools_mod + use mpi Implicit None - include 'mpif.h' type(psb_dspmat_type),intent(in) :: a type(psb_desc_type),intent(in) :: desc_a type(psb_desc_type),intent(inout) :: desc_p - integer,intent(in) :: n_ovr + integer,intent(in) :: n_ovr ! Input estimates for allocation sizes. Could we do without these two??? - Integer, Intent(in) :: l_tmp_halo,l_tmp_ovr_idx + Integer, Intent(in) :: l_tmp_halo,l_tmp_ovr_idx Integer, Intent(inout) :: lworks, lworkr integer, intent(out) :: info - type(psb_dspmat_type) :: blk - - - Integer, allocatable :: tmp_halo(:),tmp_ovr_idx(:) - + type(psb_dspmat_type) :: blk + Integer, allocatable :: tmp_halo(:),tmp_ovr_idx(:) Integer :: counter,counter_h, counter_o, counter_e,j,idx,gidx,proc,n_elem_recv,& & n_elem_send,tot_recv,tot_elem,n_col,m,ictxt,np,me,dl_lda,lwork,& & counter_t,n_elem,i_ovr,jj,i,proc_id,isz, mglob, glx,n_row, & - & idxr, idxs, lx, iszr, err_act, icomm, nxch, nsnd, nrcv - - Integer,allocatable :: halo(:),length_dl(:),works(:),workr(:),t_halo_in(:),& - & t_halo_out(:),work(:),dep_list(:),temp(:) + & idxr, idxs, lx, iszr, iszs, err_act, icomm, nxch, nsnd, nrcv,lidx + Integer,allocatable :: halo(:),works(:),workr(:),t_halo_in(:),& + & t_halo_out(:),temp(:),maskr(:) Integer,allocatable :: brvindx(:),rvsz(:), bsdindx(:),sdsz(:) Logical,Parameter :: debug=.false. @@ -114,11 +110,8 @@ Subroutine psb_dcdovrbld(n_ovr,desc_p,desc_a,a,& n_col = psb_cd_get_local_cols(desc_a) if (debug) write(0,*) me,' On entry to CDOVRBLD n_col:',n_col - dl_lda=np*5 - lwork=5*(5*np+2)*np+10 Allocate(works(lworks),workr(lworkr),t_halo_in(l_tmp_halo),& - & t_halo_out(l_tmp_halo), work(lwork),& - & length_dl(np+1),dep_list(dl_lda*np),temp(lworkr),stat=info) + & t_halo_out(l_tmp_halo), temp(lworkr),stat=info) if (info /= 0) then call psb_errpush(4010,name,a_err='Allocate') goto 9999 @@ -141,7 +134,7 @@ Subroutine psb_dcdovrbld(n_ovr,desc_p,desc_a,a,& call psb_errpush(4010,name,a_err='Allocate') goto 9999 end if - halo(:) = desc_a%halo_index(:) + halo(:) = desc_a%halo_index(:) desc_p%ovrlap_elem(:) = -1 tmp_ovr_idx(:) = -1 tmp_halo(:) = -1 @@ -150,7 +143,40 @@ Subroutine psb_dcdovrbld(n_ovr,desc_p,desc_a,a,& counter_h = 1 counter_o = 1 - ! See comment in main loop below. + ! Init overlap with desc_a%ovrlap (if any) + counter = 1 + Do While (desc_a%ovrlap_index(counter) /= -1) + proc = desc_a%ovrlap_index(counter+psb_proc_id_) + n_elem_recv = desc_a%ovrlap_index(counter+psb_n_elem_recv_) + n_elem_send = desc_a%ovrlap_index(counter+n_elem_recv+psb_n_elem_send_) + + Do j=0,n_elem_recv-1 + + idx = desc_a%ovrlap_index(counter+psb_elem_recv_+j) + If(idx > Size(desc_p%loc_to_glob)) then + info=-3 + call psb_errpush(info,name) + goto 9999 + endif + + gidx = desc_p%loc_to_glob(idx) + + call psb_check_size((counter_o+3),tmp_ovr_idx,info,pad=-1) + if (info /= 0) then + info=4010 + call psb_errpush(info,name,a_err='psb_check_size') + goto 9999 + end if + + tmp_ovr_idx(counter_o)=proc + tmp_ovr_idx(counter_o+1)=1 + tmp_ovr_idx(counter_o+2)=gidx + tmp_ovr_idx(counter_o+3)=-1 + counter_o=counter_o+3 + end Do + counter=counter+n_elem_recv+n_elem_send+2 + end Do + ! ! A picture is in order to understand what goes on here. @@ -175,7 +201,6 @@ Subroutine psb_dcdovrbld(n_ovr,desc_p,desc_a,a,& Do i_ovr=1,n_ovr if (debug) write(0,*) me,'Running on overlap level ',i_ovr,' of ',n_ovr -!!$ t_halo_in(:) = -1 ! ! At this point, halo contains a valid halo corresponding to the @@ -192,7 +217,8 @@ Subroutine psb_dcdovrbld(n_ovr,desc_p,desc_a,a,& idxs = 0 counter = 1 counter_t = 1 - + + t1 = mpi_wtime() Do While (halo(counter) /= -1) tot_elem=0 @@ -200,7 +226,7 @@ Subroutine psb_dcdovrbld(n_ovr,desc_p,desc_a,a,& n_elem_recv=halo(counter+psb_n_elem_recv_) n_elem_send=halo(counter+n_elem_recv+psb_n_elem_send_) If ((counter+n_elem_recv+n_elem_send) > Size(halo)) then - info=-1 + info = -1 call psb_errpush(info,name) goto 9999 end If @@ -210,7 +236,7 @@ Subroutine psb_dcdovrbld(n_ovr,desc_p,desc_a,a,& ! ! The format of the halo vector exists in two forms: 1. Temporary ! 2. Assembled. In this loop we are using the (assembled) halo_in and - ! copying it into (temporary) tmp_halo; this is because tmp_halo will + ! copying it into (temporary) halo_out; this is because tmp_halo will ! be enlarged with the new column indices received, and will reassemble ! everything for the next iteration. ! @@ -225,7 +251,7 @@ Subroutine psb_dcdovrbld(n_ovr,desc_p,desc_a,a,& call psb_errpush(info,name) goto 9999 end If - idx = halo(counter+psb_elem_recv_+j) + idx = halo(counter+psb_elem_recv_+j) If(idx > Size(desc_p%loc_to_glob)) then info=-3 @@ -235,42 +261,33 @@ Subroutine psb_dcdovrbld(n_ovr,desc_p,desc_a,a,& gidx = desc_p%loc_to_glob(idx) - If((counter_o+2) > Size(tmp_ovr_idx)) Then - isz = max((3*Size(tmp_ovr_idx))/2,(counter_o+3)) - if (debug) write(0,*) me,'Realloc tmp_ovr',isz - call psb_realloc(isz,tmp_ovr_idx,info,pad=-1) - if (info /= 0) then - info=4010 - ch_err='psb_realloc' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - End If + call psb_check_size((counter_o+3),tmp_ovr_idx,info,pad=-1) + if (info /= 0) then + info=4010 + call psb_errpush(info,name,a_err='psb_check_size') + goto 9999 + end if tmp_ovr_idx(counter_o)=proc tmp_ovr_idx(counter_o+1)=1 tmp_ovr_idx(counter_o+2)=gidx tmp_ovr_idx(counter_o+3)=-1 counter_o=counter_o+3 - - If((counter_h+2) > Size(tmp_halo)) Then - isz = max((3*Size(tmp_halo))/2,(counter_h+3)) - if (debug) write(0,*) me,'Realloc tmp_halo',isz - call psb_realloc(isz,tmp_halo,info) + if (.not.psb_is_large_desc(desc_p)) then + call psb_check_size((counter_h+3),tmp_halo,info,pad=-1) if (info /= 0) then info=4010 - ch_err='psb_realloc' - call psb_errpush(info,name,a_err=ch_err) + call psb_errpush(info,name,a_err='psb_check_size') goto 9999 end if - End If - - tmp_halo(counter_h)=proc - tmp_halo(counter_h+1)=1 - tmp_halo(counter_h+2)=idx - tmp_halo(counter_h+3)=-1 - - counter_h=counter_h+3 + + tmp_halo(counter_h)=proc + tmp_halo(counter_h+1)=1 + tmp_halo(counter_h+2)=idx + tmp_halo(counter_h+3)=-1 + + counter_h=counter_h+3 + end if Enddo if (debug) write(0,*) me,'Checktmp_o_i Loop Mid1',tmp_ovr_idx(1:10) @@ -283,18 +300,15 @@ Subroutine psb_dcdovrbld(n_ovr,desc_p,desc_a,a,& idx = halo(counter+psb_elem_send_+j) gidx = desc_p%loc_to_glob(idx) - - If((counter_o+2) > Size(tmp_ovr_idx)) Then - isz = max((3*Size(tmp_ovr_idx))/2,(counter_o+3)) - if (debug) write(0,*) me,'Realloc tmp_ovr',isz - call psb_realloc(isz,tmp_ovr_idx,info) - if (info /= 0) then - info=4010 - ch_err='psrealloc' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - End If + if (idx > psb_cd_get_local_rows(Desc_a)) & + & write(0,*) me,i_ovr,'Out of local rows ',idx,psb_cd_get_local_rows(Desc_a) + + call psb_check_size((counter_o+3),tmp_ovr_idx,info,pad=-1) + if (info /= 0) then + info=4010 + call psb_errpush(info,name,a_err='psb_check_size') + goto 9999 + end if tmp_ovr_idx(counter_o)=proc tmp_ovr_idx(counter_o+1)=1 @@ -308,18 +322,12 @@ Subroutine psb_dcdovrbld(n_ovr,desc_p,desc_a,a,& If (i_ovr < (n_ovr)) Then n_elem = psb_sp_get_nnz_row(idx,a) - If((idxs+tot_elem+n_elem) > lworks) Then - isz = max((3*lworks)/2,(idxs+tot_elem+n_elem)) - if (debug) write(0,*) me,'Realloc works',isz - call psb_realloc(isz,works,info) - if (info /= 0) then - info=4010 - ch_err='psb_realloc' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - lworks = isz - End If + call psb_check_size((idxs+tot_elem+n_elem),works,info) + if (info /= 0) then + info=4010 + call psb_errpush(info,name,a_err='psb_check_size') + goto 9999 + end if If((n_elem) > size(blk%ia2)) Then isz = max((3*size(blk%ia2))/2,(n_elem)) @@ -428,61 +436,109 @@ Subroutine psb_dcdovrbld(n_ovr,desc_p,desc_a,a,& goto 9999 end if - Do i=1,iszr - idx=workr(i) - if (idx <1) then - write(0,*) me,'Error in CDOVRBLD ',idx,i,iszr -!!$ write(0,*) me, ' WORKR :',workr(1:iszr) - else If (desc_p%glob_to_loc(idx) < -np) Then - ! - ! This is a new index. Assigning a local index as - ! we receive the guarantees that all indices for HALO(I) - ! will be less than those for HALO(J) whenever I Size(desc_p%loc_to_glob)) Then - isz = 3*n_col/2 - if (debug) write(0,*) me,'Realloc loc_to_glob' - call psb_realloc(isz,desc_p%loc_to_glob,info) + if (debug) write(0,*) 'ISZR :',iszr + + if (psb_is_large_desc(desc_a)) then + call psb_check_size(iszr,maskr,info) + if (info /= 0) then + info=4010 + call psb_errpush(info,name,a_err='psb_check_size') + goto 9999 + end if + call psi_idx_cnv(iszr,workr,maskr,desc_p,info) + iszs = count(maskr(1:iszr)<=0) + if (iszs > size(works)) call psb_realloc(iszs,works,info) + j = 0 + do i=1,iszr + if (maskr(i) < 0) then + j=j+1 + works(j) = workr(i) + end if + end do + ! + ! fnd_owner on desc_a because we want the procs who + ! owned the rows from the beginning! + ! + call psi_fnd_owner(iszs,works,temp,desc_a,info) + n_col=psb_cd_get_local_cols(desc_p) + + do i=1,iszs + idx = works(i) + n_col = psb_cd_get_local_cols(desc_p) + call psi_idx_ins_cnv(idx,lidx,desc_p,info) + if (psb_cd_get_local_cols(desc_p) > n_col ) then + ! + ! This is a new index. Assigning a local index as + ! we receive them guarantees that all indices for HALO(I) + ! will be less than those for HALO(J) whenever I Size(t_halo_in)) then - isz = max((3*Size(t_halo_in))/2,(counter_t+3+1000)) - if (debug) write(0,*) me,'Realloc ovr_el',isz - call psb_realloc(isz,t_halo_in,info) + + t_halo_in(counter_t)=proc_id + t_halo_in(counter_t+1)=1 + t_halo_in(counter_t+2)=lidx + counter_t=counter_t+3 + if (.false.) write(0,*) me,' CDOVRBLD: Added t_halo_in ',& + &proc_id,lidx,idx + endif + end Do + + else + + Do i=1,iszr + idx=workr(i) + if (idx <1) then + write(0,*) me,'Error in CDOVRBLD level',i_ovr,' : ',idx,i,iszr + else If (desc_p%glob_to_loc(idx) < -np) Then + ! + ! This is a new index. Assigning a local index as + ! we receive them guarantees that all indices for HALO(I) + ! will be less than those for HALO(J) whenever I). The communication descriptor. ! info - integer. Eventually returns an error code subroutine psb_dinsvi(m, irw, val, x, desc_a, info, dupl) @@ -45,6 +45,7 @@ subroutine psb_dinsvi(m, irw, val, x, desc_a, info, dupl) use psb_const_mod use psb_error_mod use psb_penv_mod + use psi_mod implicit none ! m rows number of submatrix belonging to val to be inserted @@ -64,25 +65,26 @@ subroutine psb_dinsvi(m, irw, val, x, desc_a, info, dupl) integer :: ictxt,i,loc_row,glob_row,& & loc_rows,loc_cols,mglob,err_act, int_err(5) integer :: np, me, dupl_ - character(len=20) :: name + integer, allocatable :: irl(:) + character(len=20) :: name if(psb_get_errstatus() /= 0) return info=0 call psb_erractionsave(err_act) name = 'psb_dinsvi' - if (.not.allocated(desc_a%glob_to_loc)) then - info=3110 - call psb_errpush(info,name) - return - end if +!!$ if (.not.allocated(desc_a%glob_to_loc)) then +!!$ info=3110 +!!$ call psb_errpush(info,name) +!!$ return +!!$ end if if ((.not.allocated(desc_a%matrix_data))) then int_err(1)=3110 call psb_errpush(info,name) return end if - ictxt=psb_cd_get_context(desc_a) + ictxt = psb_cd_get_context(desc_a) call psb_info(ictxt, me, np) if (np == -1) then @@ -111,15 +113,26 @@ subroutine psb_dinsvi(m, irw, val, x, desc_a, info, dupl) goto 9999 endif - loc_rows=psb_cd_get_local_rows(desc_a) - loc_cols=psb_cd_get_local_cols(desc_a) + if (m==0) return + loc_rows = psb_cd_get_local_rows(desc_a) + loc_cols = psb_cd_get_local_cols(desc_a) mglob = psb_cd_get_global_rows(desc_a) + + allocate(irl(m),stat=info) + if (info /= 0) then + info = 4000 + call psb_errpush(info,name) + goto 9999 + endif + if (present(dupl)) then dupl_ = dupl else dupl_ = psb_dupl_ovwrt_ endif + + call psi_idx_cnv(m,irw,irl,desc_a,info,owned=.true.) select case(dupl_) case(psb_dupl_ovwrt_) @@ -127,15 +140,10 @@ subroutine psb_dinsvi(m, irw, val, x, desc_a, info, dupl) !loop over all val's rows ! row actual block row - glob_row=irw(i) - if ((glob_row>0).and.(glob_row <= mglob)) then - - loc_row=desc_a%glob_to_loc(glob_row) - if (loc_row.ge.1) then - ! this row belongs to me - ! copy i-th row of block val in x - x(loc_row) = val(i) - end if + if (irl(i) > 0) then + ! this row belongs to me + ! copy i-th row of block val in x + x(irl(i)) = val(i) end if enddo @@ -144,16 +152,10 @@ subroutine psb_dinsvi(m, irw, val, x, desc_a, info, dupl) do i = 1, m !loop over all val's rows - ! row actual block row - glob_row=irw(i) - if ((glob_row>0).and.(glob_row <= mglob)) then - - loc_row=desc_a%glob_to_loc(glob_row) - if (loc_row.ge.1) then + if (irl(i) > 0) then ! this row belongs to me ! copy i-th row of block val in x - x(loc_row) = x(loc_row) + val(i) - end if + x(irl(i)) = x(irl(i)) + val(i) end if enddo @@ -162,6 +164,7 @@ subroutine psb_dinsvi(m, irw, val, x, desc_a, info, dupl) call psb_errpush(info,name) goto 9999 end select + deallocate(irl) call psb_erractionrestore(err_act) return @@ -215,8 +218,8 @@ end subroutine psb_dinsvi ! m - integer. Number of rows of submatrix belonging to ! val to be inserted. ! irw - integer(:) Row indices of rows of val (global numbering) -! val - real, pointer, dimension(:,:). The source dense submatrix. -! x - real, pointer, dimension(:,:). The destination dense matrix. +! val - real, dimension(:,:). The source dense submatrix. +! x - real, dimension(:,:). The destination dense matrix. ! desc_a - type(). The communication descriptor. ! info - integer. Eventually returns an error code subroutine psb_dinsi(m, irw, val, x, desc_a, info, dupl) @@ -225,12 +228,11 @@ subroutine psb_dinsi(m, irw, val, x, desc_a, info, dupl) use psb_const_mod use psb_error_mod use psb_penv_mod + use psi_mod implicit none ! m rows number of submatrix belonging to val to be inserted - ! ival first row of submatrix belonging to val to be inserted - ! ix x global-row corresponding to position at which val submatrix ! must be inserted @@ -238,7 +240,7 @@ subroutine psb_dinsi(m, irw, val, x, desc_a, info, dupl) integer, intent(in) :: m integer, intent(in) :: irw(:) real(kind(1.d0)), intent(in) :: val(:,:) - real(kind(1.d0)),pointer :: x(:,:) + real(kind(1.d0)), intent(inout) :: x(:,:) type(psb_desc_type), intent(in) :: desc_a integer, intent(out) :: info integer, optional, intent(in) :: dupl @@ -247,6 +249,7 @@ subroutine psb_dinsi(m, irw, val, x, desc_a, info, dupl) integer :: ictxt,i,loc_row,glob_row,j,n,& & loc_rows,loc_cols,mglob,err_act, int_err(5) integer :: np,me,dupl_ + integer, allocatable :: irl(:) character(len=20) :: name if(psb_get_errstatus() /= 0) return @@ -254,18 +257,18 @@ subroutine psb_dinsi(m, irw, val, x, desc_a, info, dupl) call psb_erractionsave(err_act) name = 'psb_dinsi' - if (.not.allocated(desc_a%glob_to_loc)) then - info=3110 - call psb_errpush(info,name) - return - end if +!!$ if (.not.allocated(desc_a%glob_to_loc)) then +!!$ info=3110 +!!$ call psb_errpush(info,name) +!!$ return +!!$ end if if ((.not.allocated(desc_a%matrix_data))) then int_err(1)=3110 call psb_errpush(info,name) return end if - ictxt=psb_cd_get_context(desc_a) + ictxt = psb_cd_get_context(desc_a) call psb_info(ictxt, me, np) if (np == -1) then @@ -293,9 +296,10 @@ subroutine psb_dinsi(m, irw, val, x, desc_a, info, dupl) call psb_errpush(info,name,int_err) goto 9999 endif + if (m==0) return - loc_rows=psb_cd_get_local_rows(desc_a) - loc_cols=psb_cd_get_local_cols(desc_a) + loc_rows = psb_cd_get_local_rows(desc_a) + loc_cols = psb_cd_get_local_cols(desc_a) mglob = psb_cd_get_global_rows(desc_a) n = min(size(val,2),size(x,2)) @@ -306,23 +310,28 @@ subroutine psb_dinsi(m, irw, val, x, desc_a, info, dupl) dupl_ = psb_dupl_ovwrt_ endif + allocate(irl(m),stat=info) + if (info /= 0) then + info = 4000 + call psb_errpush(info,name) + goto 9999 + endif + + call psi_idx_cnv(m,irw,irl,desc_a,info,owned=.true.) + select case(dupl_) case(psb_dupl_ovwrt_) do i = 1, m !loop over all val's rows ! row actual block row - glob_row=irw(i) - if ((glob_row>0).and.(glob_row <= mglob)) then - - loc_row=desc_a%glob_to_loc(glob_row) - if (loc_row.ge.1) then - ! this row belongs to me - ! copy i-th row of block val in x - do j=1,n - x(loc_row,j) = val(i,j) - end do - end if + loc_row = irl(i) + if (loc_row > 0) then + ! this row belongs to me + ! copy i-th row of block val in x + do j=1,n + x(loc_row,j) = val(i,j) + end do end if enddo @@ -332,17 +341,13 @@ subroutine psb_dinsi(m, irw, val, x, desc_a, info, dupl) !loop over all val's rows ! row actual block row - glob_row=irw(i) - if ((glob_row>0).and.(glob_row <= mglob)) then - - loc_row=desc_a%glob_to_loc(glob_row) - if (loc_row.ge.1) then - ! this row belongs to me - ! copy i-th row of block val in x - do j=1,n - x(loc_row,j) = x(loc_row,j) + val(i,j) - end do - end if + loc_row = irl(i) + if (loc_row > 0) then + ! this row belongs to me + ! copy i-th row of block val in x + do j=1,n + x(loc_row,j) = x(loc_row,j) + val(i,j) + end do end if enddo @@ -351,6 +356,7 @@ subroutine psb_dinsi(m, irw, val, x, desc_a, info, dupl) call psb_errpush(info,name) goto 9999 end select + deallocate(irl) call psb_erractionrestore(err_act) return diff --git a/src/tools/psb_dspalloc.f90 b/src/tools/psb_dspalloc.f90 index 7d2ab7ba..741972c8 100644 --- a/src/tools/psb_dspalloc.f90 +++ b/src/tools/psb_dspalloc.f90 @@ -68,8 +68,9 @@ subroutine psb_dspalloc(a, desc_a, info, nnz) call psb_erractionsave(err_act) name = 'psb_dspalloc' - ictxt = psb_cd_get_context(desc_a) + ictxt = psb_cd_get_context(desc_a) dectype = psb_cd_get_dectype(desc_a) + call psb_info(ictxt, me, np) ! ....verify blacs grid correctness.. if (np == -1) then @@ -127,7 +128,7 @@ subroutine psb_dspalloc(a, desc_a, info, nnz) a%infoa(psb_state_) = psb_spmat_bld_ if (debug) write(0,*) 'spall: ', & - &psb_cd_get_dectype(desc_a),psb_desc_bld_ + & psb_cd_get_dectype(desc_a),psb_desc_bld_ return diff --git a/src/tools/psb_dspasb.f90 b/src/tools/psb_dspasb.f90 index 17a9fd80..826d246f 100644 --- a/src/tools/psb_dspasb.f90 +++ b/src/tools/psb_dspasb.f90 @@ -76,7 +76,7 @@ subroutine psb_dspasb(a,desc_a, info, afmt, upd, dupl) name = 'psb_spasb' call psb_erractionsave(err_act) - ictxt = psb_cd_get_context(desc_a) + ictxt = psb_cd_get_context(desc_a) dscstate = psb_cd_get_dectype(desc_a) n_row = psb_cd_get_local_rows(desc_a) n_col = psb_cd_get_local_cols(desc_a) @@ -105,10 +105,10 @@ subroutine psb_dspasb(a,desc_a, info, afmt, upd, dupl) ! ! First case: we come from a fresh build. ! - n_row = psb_cd_get_local_rows(desc_a) n_col = psb_cd_get_local_cols(desc_a) + ! ! Second step: handle the local matrix part. ! diff --git a/src/tools/psb_dspcnv.f90 b/src/tools/psb_dspcnv.f90 index 7d578f90..a8056f88 100644 --- a/src/tools/psb_dspcnv.f90 +++ b/src/tools/psb_dspcnv.f90 @@ -121,8 +121,7 @@ subroutine psb_dspcnv(a,b,desc_a,info) time(1) = mpi_wtime() - - ictxt = psb_cd_get_context(desc_a) + ictxt = psb_cd_get_context(desc_a) dectype = psb_cd_get_dectype(desc_a) n_row = psb_cd_get_local_rows(desc_a) n_col = psb_cd_get_local_cols(desc_a) diff --git a/src/tools/psb_dspfree.f90 b/src/tools/psb_dspfree.f90 index 978a025c..62751831 100644 --- a/src/tools/psb_dspfree.f90 +++ b/src/tools/psb_dspfree.f90 @@ -65,7 +65,7 @@ subroutine psb_dspfree(a, desc_a,info) call psb_errpush(info,name) return else - ictxt=psb_cd_get_context(desc_a) + ictxt = psb_cd_get_context(desc_a) end if !...deallocate a.... diff --git a/src/tools/psb_dsphalo.f90 b/src/tools/psb_dsphalo.f90 index 0fd34575..eb77ff0f 100644 --- a/src/tools/psb_dsphalo.f90 +++ b/src/tools/psb_dsphalo.f90 @@ -73,7 +73,7 @@ Subroutine psb_dsphalo(a,desc_a,blk,info,rwcnv,clcnv,outfmt) & rvsz(:), bsdindx(:),sdsz(:) logical :: rwcnv_,clcnv_ character(len=5) :: outfmt_ - Logical,Parameter :: debug=.false., usea2av=.true. + Logical,Parameter :: debug=.false., debugprt=.false. real(kind(1.d0)) :: t1,t2,t3,t4,t5,t6,t7,t8,t9 character(len=20) :: name, ch_err @@ -100,7 +100,8 @@ Subroutine psb_dsphalo(a,desc_a,blk,info,rwcnv,clcnv,outfmt) outfmt_ = 'CSR' endif - ictxt=psb_cd_get_context(desc_a) + ictxt = psb_cd_get_context(desc_a) + Call psb_info(ictxt, me, np) t1 = mpi_wtime() @@ -225,8 +226,6 @@ Subroutine psb_dsphalo(a,desc_a,blk,info,rwcnv,clcnv,outfmt) counter = counter+n_el_send+3 Enddo nz = tmp%infoa(psb_nnz_) -!!$ call csprt(20+me,tmp,head='% SPHALO border SEND .') -!!$ close(20+me) if (rwcnv_) call psb_loc_to_glob(tmp%ia1(1:nz),desc_a,info,iact='I') if (clcnv_) call psb_loc_to_glob(tmp%ia2(1:nz),desc_a,info,iact='I') @@ -236,8 +235,12 @@ Subroutine psb_dsphalo(a,desc_a,blk,info,rwcnv,clcnv,outfmt) call psb_errpush(info,name,a_err=ch_err) goto 9999 end if -!!$ call csprt(30+me,tmp,head='% SPHALO border SEND .') -!!$ close(30+me) + if (debugprt) then + open(30+me) + call psb_csprt(30+me,tmp,head='% SPHALO border SEND .') + call flush(30+me) + close(30+me) + end if call mpi_alltoallv(tmp%aspk,sdsz,bsdindx,mpi_double_precision,& @@ -261,6 +264,7 @@ Subroutine psb_dsphalo(a,desc_a,blk,info,rwcnv,clcnv,outfmt) ! if (rwcnv_) call psb_glob_to_loc(blk%ia1(1:iszr),desc_a,info,iact='I') if (clcnv_) call psb_glob_to_loc(blk%ia2(1:iszr),desc_a,info,iact='I') + if (info /= 0) then info=4010 ch_err='psbglob_to_loc' @@ -268,6 +272,14 @@ Subroutine psb_dsphalo(a,desc_a,blk,info,rwcnv,clcnv,outfmt) goto 9999 end if + if (debugprt) then + blk%fida='COO' + blk%infoa(psb_nnz_)=iszr + open(40+me) + call psb_csprt(40+me,blk,head='% SPHALO border .') + call flush(40+me) + close(40+me) + end if l1 = 0 Do i=1,iszr !!$ write(0,*) work5(i),work6(i) @@ -283,9 +295,13 @@ Subroutine psb_dsphalo(a,desc_a,blk,info,rwcnv,clcnv,outfmt) Enddo blk%fida='COO' blk%infoa(psb_nnz_)=l1 -!!$ open(50+me) -!!$ call csprt(50+me,blk,head='% SPHALO border .') -!!$ close(50+me) + if (debugprt) then + open(50+me) + call psb_csprt(50+me,blk,head='% SPHALO border .') + call flush(50+me) + close(50+me) + call psb_barrier(ictxt) + end if t4 = mpi_wtime() if(debug) Write(0,*)me,'End first loop',counter,l1,blk%m diff --git a/src/tools/psb_dspins.f90 b/src/tools/psb_dspins.f90 index 03619411..418e1262 100644 --- a/src/tools/psb_dspins.f90 +++ b/src/tools/psb_dspins.f90 @@ -63,22 +63,34 @@ subroutine psb_dspins(nz,ia,ja,val,a,desc_a,info,rebuild) logical, intent(in), optional :: rebuild !locals..... - integer :: nrow, err_act, dectype,mglob,ncol, spstate + integer :: nrow, err_act,mglob,ncol, spstate integer :: ictxt,np,me logical, parameter :: debug=.false. integer, parameter :: relocsz=200 logical :: rebuild_ + integer, allocatable :: ila(:),jla(:) interface psb_cdins - subroutine psb_cdins(nz,ia,ja,desc_a,info) + subroutine psb_cdins(nz,ia,ja,desc_a,info,ila,jla) use psb_descriptor_type implicit none type(psb_desc_type), intent(inout) :: desc_a - integer, intent(in) :: nz,ia(:),ja(:) - integer, intent(out) :: info + integer, intent(in) :: nz,ia(:),ja(:) + integer, intent(out) :: info + integer, optional, intent(out) :: ila(:), jla(:) end subroutine psb_cdins end interface + interface psb_glob_to_loc + subroutine psb_glob_to_loc(x,desc_a,info,iact) + use psb_descriptor_type + implicit none + type(psb_desc_type), intent(in) :: desc_a + integer, intent(inout) :: x(:) + integer, intent(out) :: info + character, intent(in), optional :: iact + end subroutine psb_glob_to_loc + end interface character(len=20) :: name, ch_err info = 0 @@ -87,8 +99,7 @@ subroutine psb_dspins(nz,ia,ja,val,a,desc_a,info,rebuild) ictxt = psb_cd_get_context(desc_a) - dectype = psb_cd_get_dectype(desc_a) - mglob = psb_cd_get_global_rows(desc_a) + mglob = psb_cd_get_global_rows(desc_a) call psb_info(ictxt, me, np) @@ -98,7 +109,7 @@ subroutine psb_dspins(nz,ia,ja,val,a,desc_a,info,rebuild) goto 9999 endif - if (nz <= 0) then + if (nz < 0) then info = 1111 call psb_errpush(info,name) goto 9999 @@ -119,6 +130,7 @@ subroutine psb_dspins(nz,ia,ja,val,a,desc_a,info,rebuild) call psb_errpush(info,name) goto 9999 end if + if (nz==0) return if (present(rebuild)) then rebuild_ = rebuild @@ -128,39 +140,104 @@ subroutine psb_dspins(nz,ia,ja,val,a,desc_a,info,rebuild) spstate = a%infoa(psb_state_) if (psb_is_bld_desc(desc_a)) then - call psb_cdins(nz,ia,ja,desc_a,info) - if (info /= 0) then - info=4010 - ch_err='psb_cdins' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 + if (psb_is_large_desc(desc_a)) then + + allocate(ila(nz),jla(nz),stat=info) + if (info /= 0) then + ch_err='allocate' + call psb_errpush(4013,name,a_err=ch_err,i_err=(/info,0,0,0,0/)) + goto 9999 + end if + call psb_cdins(nz,ia,ja,desc_a,info,ila=ila,jla=jla) + if (info /= 0) then + ch_err='psb_cdins' + call psb_errpush(4013,name,a_err=ch_err,i_err=(/info,0,0,0,0/)) + goto 9999 + end if + nrow = psb_cd_get_local_rows(desc_a) + ncol = psb_cd_get_local_cols(desc_a) + + if (spstate == psb_spmat_bld_) then + call psb_coins(nz,ila,jla,val,a,1,nrow,1,ncol,info) + if (info /= 0) then + info=4010 + ch_err='psb_coins' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + else + info = 1123 + call psb_errpush(info,name) + goto 9999 + end if + + else + + call psb_cdins(nz,ia,ja,desc_a,info) + if (info /= 0) then + ch_err='psb_cdins' + call psb_errpush(4013,name,a_err=ch_err,i_err=(/info,0,0,0,0/)) + goto 9999 + end if + nrow = psb_cd_get_local_rows(desc_a) + ncol = psb_cd_get_local_cols(desc_a) + + if (spstate == psb_spmat_bld_) then + call psb_coins(nz,ia,ja,val,a,1,nrow,1,ncol,info,gtl=desc_a%glob_to_loc) + if (info /= 0) then + info=4010 + ch_err='psb_coins' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + else + info = 1123 + call psb_errpush(info,name) + goto 9999 + end if + end if nrow = psb_cd_get_local_rows(desc_a) ncol = psb_cd_get_local_cols(desc_a) - if (spstate == psb_spmat_bld_) then - call psb_coins(nz,ia,ja,val,a,1,nrow,1,ncol,info,gtl=desc_a%glob_to_loc) + else if (psb_is_asb_desc(desc_a)) then + + if (psb_is_large_desc(desc_a)) then + + allocate(ila(nz),jla(nz),stat=info) + if (info /= 0) then + ch_err='allocate' + call psb_errpush(4013,name,a_err=ch_err,i_err=(/info,0,0,0,0/)) + goto 9999 + end if + + ila(1:nz) = ia(1:nz) + jla(1:nz) = ja(1:nz) + call psb_glob_to_loc(ila(1:nz),desc_a,info,iact='I') + call psb_glob_to_loc(jla(1:nz),desc_a,info,iact='I') + nrow = psb_cd_get_local_rows(desc_a) + ncol = psb_cd_get_local_cols(desc_a) + + call psb_coins(nz,ia,ja,val,a,1,nrow,1,ncol,& + & info,rebuild=rebuild_) if (info /= 0) then info=4010 ch_err='psb_coins' call psb_errpush(info,name,a_err=ch_err) goto 9999 end if + else - info = 1123 - call psb_errpush(info,name) - goto 9999 - end if - else if (psb_is_asb_desc(desc_a)) then - nrow = psb_cd_get_local_rows(desc_a) - ncol = psb_cd_get_local_cols(desc_a) - call psb_coins(nz,ia,ja,val,a,1,nrow,1,ncol,& - & info,gtl=desc_a%glob_to_loc,rebuild=rebuild_) - if (info /= 0) then - info=4010 - ch_err='psb_coins' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 + nrow = psb_cd_get_local_rows(desc_a) + ncol = psb_cd_get_local_cols(desc_a) + call psb_coins(nz,ia,ja,val,a,1,nrow,1,ncol,& + & info,gtl=desc_a%glob_to_loc,rebuild=rebuild_) + if (info /= 0) then + info=4010 + ch_err='psb_coins' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if end if else info = 1122 diff --git a/src/tools/psb_dsprn.f90 b/src/tools/psb_dsprn.f90 index 53dc0a66..e905b157 100644 --- a/src/tools/psb_dsprn.f90 +++ b/src/tools/psb_dsprn.f90 @@ -70,6 +70,7 @@ Subroutine psb_dsprn(a, desc_a,info,clear) call psb_erractionsave(err_act) ictxt = psb_cd_get_context(desc_a) + call psb_info(ictxt, me, np) if (debug) & &write(*,*) 'starting spalloc ',ictxt,np,me diff --git a/src/tools/psb_glob_to_loc.f90 b/src/tools/psb_glob_to_loc.f90 index 4fa6c273..b35d5015 100644 --- a/src/tools/psb_glob_to_loc.f90 +++ b/src/tools/psb_glob_to_loc.f90 @@ -46,6 +46,8 @@ subroutine psb_glob_to_loc2(x,y,desc_a,info,iact) use psb_const_mod use psb_error_mod use psb_string_mod + use psb_penv_mod + use psi_mod implicit none !...parameters.... @@ -77,47 +79,23 @@ subroutine psb_glob_to_loc2(x,y,desc_a,info,iact) int_err=0 real_val = 0.d0 - n=size(x) - do i=1,n - if ((x(i).gt.psb_cd_get_global_rows(desc_a)).or.& - & (x(i).le.zero)) then - if (act == 'I') then - y(i)=-3*psb_cd_get_global_rows(desc_a) - else - info=140 - int_err(1)=x(i) - int_err(2)=psb_cd_get_global_rows(desc_a) - exit - end if - else - tmp=desc_a%glob_to_loc(x(i)) - if((tmp.gt.zero).or.(tmp.le.psb_cd_get_local_cols(desc_a))) then - y(i)=tmp - else if (tmp.le.zero) then - info = 150 - int_err(1)=tmp - exit - else if (tmp.gt.psb_cd_get_local_cols(desc_a)) then - info = 140 - int_err(1)=tmp - int_err(2)=psb_cd_get_local_cols(desc_a) - exit - end if - end if - enddo - - if (info /= 0) then - select case(act) - case('E','I') - call psb_erractionrestore(err_act) - return - case('W') + n = size(x) + call psi_idx_cnv(n,x,y,desc_a,info) + + select case(act) + case('E','I') + call psb_erractionrestore(err_act) + return + case('W') + if ((info /= 0).or.(count(y(1:n)<0) >0)) then write(0,'("Error ",i5," in subroutine glob_to_loc")') info - case('A') + end if + case('A') + if ((info /= 0).or.(count(y(1:n)<0) >0)) then call psb_errpush(info,name) goto 9999 - end select - endif + end if + end select call psb_erractionrestore(err_act) return @@ -177,31 +155,37 @@ end subroutine psb_glob_to_loc2 ! subroutine psb_glob_to_loc(x,desc_a,info,iact) + use psb_penv_mod use psb_descriptor_type use psb_const_mod use psb_error_mod use psb_string_mod + use psi_mod implicit none !...parameters.... - type(psb_desc_type), intent(in) :: desc_a - integer, intent(inout) :: x(:) - integer, intent(out) :: info - character, intent(in), optional :: iact + type(psb_desc_type), intent(in) :: desc_a + integer, intent(inout) :: x(:) + integer, intent(out) :: info + character, intent(in), optional :: iact !....locals.... - integer :: n, i, tmp - character :: act - integer :: int_err(5), err_act - real(kind(1.d0)) :: real_val - integer, parameter :: zero=0 + integer :: n, i, tmp, nk, key, idx, ih, nh, lb, ub, lm + character :: act + integer :: int_err(5), err_act, dectype + real(kind(1.d0)) :: real_val, t0, t1,t2 + integer, parameter :: zero=0 character(len=20) :: name + integer :: ictxt, iam, np if(psb_get_errstatus() /= 0) return info=0 name = 'glob_to_loc' + ictxt = desc_a%matrix_data(psb_ctxt_) + call psb_info(ictxt,iam,np) call psb_erractionsave(err_act) + dectype = desc_a%matrix_data(psb_dec_type_) if (present(iact)) then act=iact else @@ -210,48 +194,24 @@ subroutine psb_glob_to_loc(x,desc_a,info,iact) act = toupper(act) - real_val = 0.d0 - n=size(x) - do i=1,n - if ((x(i).gt.psb_cd_get_global_rows(desc_a)).or.& - & (x(i).le.zero)) then - if(act == 'I') then - x(i)=-3*psb_cd_get_global_rows(desc_a) - else - info=140 - int_err(1)=x(i) - int_err(2)=psb_cd_get_global_rows(desc_a) - exit - end if - else - tmp=desc_a%glob_to_loc(x(i)) - if((tmp.gt.zero).or.(tmp.le.psb_cd_get_local_cols(desc_a))) then - x(i)=tmp - else if (tmp.le.zero) then - info = 150 - int_err(1)=tmp - exit - else if (tmp.ge.psb_cd_get_local_cols(desc_a)) then - info = 140 - int_err(1)=tmp - int_err(2)=psb_cd_get_local_cols(desc_a) - exit - end if - end if - enddo - - if (info /= 0) then - select case(act) - case('E','I') - call psb_erractionrestore(err_act) - return - case('W') + n = size(x) + call psi_idx_cnv(n,x,desc_a,info) + + select case(act) + case('E','I') + call psb_erractionrestore(err_act) + return + case('W') + if ((info /= 0).or.(count(x(1:n)<0) >0)) then write(0,'("Error ",i5," in subroutine glob_to_loc")') info - case('A') + end if + case('A') + if ((info /= 0).or.(count(x(1:n)<0) >0)) then + write(0,*) count(x(1:n)<0) call psb_errpush(info,name) goto 9999 - end select - endif + end if + end select call psb_erractionrestore(err_act) return @@ -266,5 +226,69 @@ subroutine psb_glob_to_loc(x,desc_a,info,iact) end if return +contains + + subroutine inlbsrch(ipos,key,n,v) + implicit none + integer ipos, key, n + integer v(n) + + integer lb, ub, m + + + lb = 1 + ub = n + ipos = -1 + + do + if (lb > ub) return + m = (lb+ub)/2 + if (key.eq.v(m)) then + ipos = m + return + else if (key.lt.v(m)) then + ub = m-1 + else + lb = m + 1 + end if + enddo + return + end subroutine inlbsrch + + subroutine inner_cnv(n,x,hashsize,hashmask,hashv,glb_lc) + integer :: n, hashsize,hashmask,x(:), hashv(0:),glb_lc(:,:) + integer :: i, ih, key, idx,nh,tmp,lb,ub,lm + do i=1, n + key = x(i) + ih = iand(key,hashmask) + idx = hashv(ih) + nh = hashv(ih+1) - hashv(ih) + if (nh > 0) then + tmp = -1 + lb = idx + ub = idx+nh-1 + do + if (lb>ub) exit + lm = (lb+ub)/2 + if (key==glb_lc(lm,1)) then + tmp = lm + exit + else if (key 0) then + x(i) = glb_lc(tmp,2) + else + x(i) = tmp + end if + end do + end subroutine inner_cnv + end subroutine psb_glob_to_loc diff --git a/src/tools/psb_ialloc.f90 b/src/tools/psb_ialloc.f90 index 669117bc..def8b647 100644 --- a/src/tools/psb_ialloc.f90 +++ b/src/tools/psb_ialloc.f90 @@ -55,9 +55,9 @@ subroutine psb_ialloc(x, desc_a, info, n) !locals - integer :: np,me,n_col,n_row,i,j,err_act - integer :: ictxt,dectype,n_ - integer :: int_err(5),exch(3) + integer :: np,me,err,n_col,n_row,i,j,err_act + integer :: ictxt,n_ + integer :: int_err(5), exch(3) character(len=20) :: name, ch_err if(psb_get_errstatus() /= 0) return @@ -74,7 +74,6 @@ subroutine psb_ialloc(x, desc_a, info, n) goto 9999 endif - dectype=psb_cd_get_dectype(desc_a) !... check m and n parameters.... if (.not.psb_is_ok_desc(desc_a)) then info = 3110 @@ -194,7 +193,7 @@ subroutine psb_iallocv(x, desc_a, info,n) use psb_penv_mod implicit none - + !....parameters... integer, allocatable, intent(out) :: x(:) type(psb_desc_type), intent(in) :: desc_a @@ -202,7 +201,7 @@ subroutine psb_iallocv(x, desc_a, info,n) integer, optional, intent(in) :: n !locals - integer :: np,me,n_col,n_row,dectype,err_act + integer :: np,me,n_col,n_row,err_act integer :: ictxt, n_ integer :: int_err(5) logical, parameter :: debug=.false. @@ -212,7 +211,7 @@ subroutine psb_iallocv(x, desc_a, info,n) info=0 name='psb_iallocv' call psb_erractionsave(err_act) - + ictxt=psb_cd_get_context(desc_a) call psb_info(ictxt, me, np) @@ -223,9 +222,6 @@ subroutine psb_iallocv(x, desc_a, info,n) goto 9999 endif - dectype=psb_cd_get_dectype(desc_a) - if (debug) write(0,*) 'dall: dectype',dectype - if (debug) write(0,*) 'dall: is_ok? dectype',psb_is_ok_desc(desc_a) !... check m and n parameters.... if (.not.psb_is_ok_desc(desc_a)) then info = 3110 @@ -237,25 +233,25 @@ subroutine psb_iallocv(x, desc_a, info,n) !....allocate x ..... if (psb_is_asb_desc(desc_a).or.psb_is_upd_desc(desc_a)) then - n_col = max(1,psb_cd_get_local_cols(desc_a)) - allocate(x(n_col),stat=info) - if (info.ne.0) then - info=2025 - int_err(1)=n_col - call psb_errpush(info,name,int_err) - goto 9999 - endif + n_col = max(1,psb_cd_get_local_cols(desc_a)) + allocate(x(n_col),stat=info) + if (info.ne.0) then + info=2025 + int_err(1)=n_col + call psb_errpush(info,name,int_err) + goto 9999 + endif else if (psb_is_bld_desc(desc_a)) then - n_row = max(1,psb_cd_get_local_rows(desc_a)) - allocate(x(n_row),stat=info) - if (info.ne.0) then - info=2025 - int_err(1)=n_row - call psb_errpush(info,name,int_err) - goto 9999 - endif + n_row = max(1,psb_cd_get_local_rows(desc_a)) + allocate(x(n_row),stat=info) + if (info.ne.0) then + info=2025 + int_err(1)=n_row + call psb_errpush(info,name,int_err) + goto 9999 + endif endif - + x = 0 call psb_erractionrestore(err_act) diff --git a/src/tools/psb_iasb.f90 b/src/tools/psb_iasb.f90 index 018fdba4..2ab47475 100644 --- a/src/tools/psb_iasb.f90 +++ b/src/tools/psb_iasb.f90 @@ -53,7 +53,7 @@ subroutine psb_iasb(x, desc_a, info) ! local variables integer :: ictxt,np,me,nrow,ncol,err_act - integer :: int_err(5), i1sz, i2sz, dectype + integer :: int_err(5), i1sz, i2sz logical, parameter :: debug=.false. character(len=20) :: name @@ -68,14 +68,13 @@ subroutine psb_iasb(x, desc_a, info) return endif - ictxt=psb_cd_get_context(desc_a) - dectype=psb_cd_get_dectype(desc_a) + ictxt = psb_cd_get_context(desc_a) call psb_info(ictxt, me, np) if (debug) write(*,*) 'asb start: ',np,me,& - &psb_cd_get_dectype(desc_a) + & psb_cd_get_dectype(desc_a) ! ....verify blacs grid correctness.. if (np == -1) then info = 2010 @@ -83,16 +82,16 @@ subroutine psb_iasb(x, desc_a, info) goto 9999 else if (.not.psb_is_asb_desc(desc_a)) then if (debug) write(*,*) 'asb error ',& - &dectype + & psb_cd_get_dectype(desc_a) info = 3110 call psb_errpush(info,name) goto 9999 endif ! check size - ictxt=psb_cd_get_context(desc_a) - nrow=psb_cd_get_local_rows(desc_a) - ncol=psb_cd_get_local_cols(desc_a) + ictxt = psb_cd_get_context(desc_a) + nrow = psb_cd_get_local_rows(desc_a) + ncol = psb_cd_get_local_cols(desc_a) i1sz = size(x,dim=1) i2sz = size(x,dim=2) if (debug) write(*,*) 'asb: ',i1sz,i2sz,nrow,ncol @@ -177,7 +176,7 @@ subroutine psb_iasbv(x, desc_a, info) ! local variables integer :: ictxt,np,me, err_act - integer :: int_err(5), i1sz,nrow,ncol, dectype + integer :: int_err(5), i1sz,nrow,ncol logical, parameter :: debug=.false. character(len=20) :: name @@ -187,8 +186,7 @@ subroutine psb_iasbv(x, desc_a, info) name = 'psb_iasbv' - ictxt=psb_cd_get_context(desc_a) - dectype=psb_cd_get_dectype(desc_a) + ictxt = psb_cd_get_context(desc_a) call psb_info(ictxt, me, np) @@ -203,8 +201,9 @@ subroutine psb_iasbv(x, desc_a, info) goto 9999 endif - nrow=psb_cd_get_local_rows(desc_a) - ncol=psb_cd_get_local_cols(desc_a) + nrow = psb_cd_get_local_rows(desc_a) + ncol = psb_cd_get_local_cols(desc_a) + if (debug) write(*,*) name,' sizes: ',nrow,ncol i1sz = size(x) if (debug) write(*,*) 'dasb: sizes ',i1sz,ncol diff --git a/src/tools/psb_iins.f90 b/src/tools/psb_iins.f90 index 55cfc063..fb6e3a9c 100644 --- a/src/tools/psb_iins.f90 +++ b/src/tools/psb_iins.f90 @@ -36,7 +36,7 @@ ! val to be inserted. ! irw - integer(:) Row indices of rows of val (global numbering) ! val - integer, dimension(:). The source dense submatrix. -! x - integer, pointer, dimension(:). The destination dense matrix. +! x - integer, dimension(:). The destination dense matrix. ! desc_a - type(). The communication descriptor. ! info - integer. Eventually returns an error code subroutine psb_iinsvi(m, irw, val, x, desc_a, info, dupl) @@ -45,6 +45,7 @@ subroutine psb_iinsvi(m, irw, val, x, desc_a, info, dupl) use psb_const_mod use psb_error_mod use psb_penv_mod + use psi_mod implicit none ! m rows number of submatrix belonging to val to be inserted @@ -62,20 +63,21 @@ subroutine psb_iinsvi(m, irw, val, x, desc_a, info, dupl) !locals..... integer :: ictxt,i,loc_row,glob_row,& - & loc_rows,loc_cols,mglob,err_act, int_err(5), err - integer :: np, me,dupl_ - character(len=20) :: name + & loc_rows,loc_cols,mglob,err_act, int_err(5) + integer :: np, me, dupl_ + integer, allocatable :: irl(:) + character(len=20) :: name if(psb_get_errstatus() /= 0) return info=0 call psb_erractionsave(err_act) - name = 'psb_dinsvv' + name = 'psb_insvi' - if (.not.allocated(desc_a%glob_to_loc)) then - info=3110 - call psb_errpush(info,name) - return - end if +!!$ if (.not.allocated(desc_a%glob_to_loc)) then +!!$ info=3110 +!!$ call psb_errpush(info,name) +!!$ return +!!$ end if if ((.not.allocated(desc_a%matrix_data))) then int_err(1)=3110 call psb_errpush(info,name) @@ -111,15 +113,25 @@ subroutine psb_iinsvi(m, irw, val, x, desc_a, info, dupl) goto 9999 endif - loc_rows=psb_cd_get_local_rows(desc_a) - loc_cols=psb_cd_get_local_cols(desc_a) + if (m==0) return + loc_rows = psb_cd_get_local_rows(desc_a) + loc_cols = psb_cd_get_local_cols(desc_a) mglob = psb_cd_get_global_rows(desc_a) + allocate(irl(m),stat=info) + if (info /= 0) then + info = 4000 + call psb_errpush(info,name) + goto 9999 + endif + if (present(dupl)) then dupl_ = dupl else dupl_ = psb_dupl_ovwrt_ endif + + call psi_idx_cnv(m,irw,irl,desc_a,info,owned=.true.) select case(dupl_) case(psb_dupl_ovwrt_) @@ -127,15 +139,10 @@ subroutine psb_iinsvi(m, irw, val, x, desc_a, info, dupl) !loop over all val's rows ! row actual block row - glob_row=irw(i) - if ((glob_row>0).and.(glob_row <= mglob)) then - - loc_row=desc_a%glob_to_loc(glob_row) - if (loc_row.ge.1) then - ! this row belongs to me - ! copy i-th row of block val in x - x(loc_row) = val(i) - end if + if (irl(i) > 0) then + ! this row belongs to me + ! copy i-th row of block val in x + x(irl(i)) = val(i) end if enddo @@ -144,16 +151,10 @@ subroutine psb_iinsvi(m, irw, val, x, desc_a, info, dupl) do i = 1, m !loop over all val's rows - ! row actual block row - glob_row=irw(i) - if ((glob_row>0).and.(glob_row <= mglob)) then - - loc_row=desc_a%glob_to_loc(glob_row) - if (loc_row.ge.1) then + if (irl(i) > 0) then ! this row belongs to me ! copy i-th row of block val in x - x(loc_row) = x(loc_row) + val(i) - end if + x(irl(i)) = x(irl(i)) + val(i) end if enddo @@ -162,6 +163,7 @@ subroutine psb_iinsvi(m, irw, val, x, desc_a, info, dupl) call psb_errpush(info,name) goto 9999 end select + deallocate(irl) call psb_erractionrestore(err_act) return @@ -216,7 +218,7 @@ end subroutine psb_iinsvi ! val to be inserted. ! irw - integer(:) Row indices of rows of val (global numbering) ! val - integer, dimension(:,:). The source dense submatrix. -! x - integer, pointer, dimension(:,:). The destination dense matrix. +! x - integer, dimension(:,:). The destination dense matrix. ! desc_a - type(). The communication descriptor. ! info - integer. Eventually returns an error code subroutine psb_iinsi(m,irw, val, x, desc_a, info, dupl) @@ -225,6 +227,7 @@ subroutine psb_iinsi(m,irw, val, x, desc_a, info, dupl) use psb_const_mod use psb_error_mod use psb_penv_mod + use psi_mod implicit none ! m rows number of submatrix belonging to val to be inserted @@ -244,19 +247,20 @@ subroutine psb_iinsi(m,irw, val, x, desc_a, info, dupl) !locals..... integer :: ictxt,i,loc_row,glob_row,j,n,& & loc_rows,loc_cols,mglob,err_act, int_err(5) - integer :: np, me,dupl_ + integer :: np,me,dupl_ + integer, allocatable :: irl(:) character(len=20) :: name if(psb_get_errstatus() /= 0) return info=0 call psb_erractionsave(err_act) - name = 'psb_dinsvv' + name = 'psb_iinsi' - if (.not.allocated(desc_a%glob_to_loc)) then - info=3110 - call psb_errpush(info,name) - return - end if +!!$ if (.not.allocated(desc_a%glob_to_loc)) then +!!$ info=3110 +!!$ call psb_errpush(info,name) +!!$ return +!!$ end if if ((.not.allocated(desc_a%matrix_data))) then int_err(1)=3110 call psb_errpush(info,name) @@ -291,9 +295,10 @@ subroutine psb_iinsi(m,irw, val, x, desc_a, info, dupl) call psb_errpush(info,name,int_err) goto 9999 endif + if (m==0) return - loc_rows=psb_cd_get_local_rows(desc_a) - loc_cols=psb_cd_get_local_cols(desc_a) + loc_rows = psb_cd_get_local_rows(desc_a) + loc_cols = psb_cd_get_local_cols(desc_a) mglob = psb_cd_get_global_rows(desc_a) n = min(size(val,2),size(x,2)) @@ -304,23 +309,28 @@ subroutine psb_iinsi(m,irw, val, x, desc_a, info, dupl) dupl_ = psb_dupl_ovwrt_ endif + allocate(irl(m),stat=info) + if (info /= 0) then + info = 4000 + call psb_errpush(info,name) + goto 9999 + endif + + call psi_idx_cnv(m,irw,irl,desc_a,info,owned=.true.) + select case(dupl_) case(psb_dupl_ovwrt_) do i = 1, m !loop over all val's rows ! row actual block row - glob_row=irw(i) - if ((glob_row>0).and.(glob_row <= mglob)) then - - loc_row=desc_a%glob_to_loc(glob_row) - if (loc_row.ge.1) then - ! this row belongs to me - ! copy i-th row of block val in x - do j=1,n - x(loc_row,j) = val(i,j) - end do - end if + loc_row = irl(i) + if (loc_row > 0) then + ! this row belongs to me + ! copy i-th row of block val in x + do j=1,n + x(loc_row,j) = val(i,j) + end do end if enddo @@ -330,17 +340,13 @@ subroutine psb_iinsi(m,irw, val, x, desc_a, info, dupl) !loop over all val's rows ! row actual block row - glob_row=irw(i) - if ((glob_row>0).and.(glob_row <= mglob)) then - - loc_row=desc_a%glob_to_loc(glob_row) - if (loc_row.ge.1) then - ! this row belongs to me - ! copy i-th row of block val in x - do j=1,n - x(loc_row,j) = x(loc_row,j) + val(i,j) - end do - end if + loc_row = irl(i) + if (loc_row > 0) then + ! this row belongs to me + ! copy i-th row of block val in x + do j=1,n + x(loc_row,j) = x(loc_row,j) + val(i,j) + end do end if enddo @@ -349,6 +355,7 @@ subroutine psb_iinsi(m,irw, val, x, desc_a, info, dupl) call psb_errpush(info,name) goto 9999 end select + deallocate(irl) call psb_erractionrestore(err_act) return diff --git a/src/tools/psb_zallc.f90 b/src/tools/psb_zallc.f90 index 96faa1e3..36ac2911 100644 --- a/src/tools/psb_zallc.f90 +++ b/src/tools/psb_zallc.f90 @@ -56,7 +56,7 @@ subroutine psb_zalloc(x, desc_a, info, n) !locals integer :: np,me,err,n_col,n_row,i,j,err_act - integer :: ictxt,dectype,n_ + integer :: ictxt,n_ integer :: int_err(5),exch(3) character(len=20) :: name, ch_err @@ -76,7 +76,6 @@ subroutine psb_zalloc(x, desc_a, info, n) goto 9999 endif - dectype=psb_cd_get_dectype(desc_a) !... check m and n parameters.... if (.not.psb_is_ok_desc(desc_a)) then info = 3110 @@ -202,7 +201,7 @@ subroutine psb_zallocv(x, desc_a,info,n) integer, optional, intent(in) :: n !locals - integer :: np,me,n_col,n_row,dectype,i,err_act + integer :: np,me,n_col,n_row,i,err_act integer :: ictxt, n_ logical, parameter :: debug=.false. character(len=20) :: name, ch_err @@ -222,9 +221,7 @@ subroutine psb_zallocv(x, desc_a,info,n) goto 9999 endif - dectype=psb_cd_get_dectype(desc_a) - if (debug) write(0,*) 'dall: dectype',dectype - if (debug) write(0,*) 'dall: is_ok? dectype',psb_is_ok_desc(desc_a) + if (debug) write(0,*) 'dall: is_ok?',psb_is_ok_desc(desc_a) !... check m and n parameters.... if (.not.psb_is_ok_desc(desc_a)) then info = 3110 diff --git a/src/tools/psb_zasb.f90 b/src/tools/psb_zasb.f90 index c347c78e..4c73369e 100644 --- a/src/tools/psb_zasb.f90 +++ b/src/tools/psb_zasb.f90 @@ -53,7 +53,7 @@ subroutine psb_zasb(x, desc_a, info) ! local variables integer :: ictxt,np,me,nrow,ncol, err_act - integer :: int_err(5), i1sz, i2sz, dectype + integer :: int_err(5), i1sz, i2sz logical, parameter :: debug=.false. character(len=20) :: name, ch_err @@ -69,13 +69,11 @@ subroutine psb_zasb(x, desc_a, info) endif ictxt=psb_cd_get_context(desc_a) - dectype=psb_cd_get_dectype(desc_a) call psb_info(ictxt, me, np) - if (debug) write(*,*) 'asb start: ',np,me,& - &psb_cd_get_dectype(desc_a) + & psb_cd_get_dectype(desc_a) ! ....verify blacs grid correctness.. if (np == -1) then info = 2010 @@ -83,16 +81,16 @@ subroutine psb_zasb(x, desc_a, info) goto 9999 else if (.not.psb_is_asb_desc(desc_a)) then if (debug) write(*,*) 'asb error ',& - &dectype + & psb_cd_get_dectype(desc_a) info = 3110 call psb_errpush(info,name,i_err=int_err) goto 9999 endif ! check size - ictxt=psb_cd_get_context(desc_a) - nrow=psb_cd_get_local_rows(desc_a) - ncol=psb_cd_get_local_cols(desc_a) + ictxt = psb_cd_get_context(desc_a) + nrow = psb_cd_get_local_rows(desc_a) + ncol = psb_cd_get_local_cols(desc_a) i1sz = size(x,dim=1) i2sz = size(x,dim=2) if (debug) write(*,*) 'asb: ',i1sz,i2sz,nrow,ncol @@ -182,7 +180,7 @@ subroutine psb_zasbv(x, desc_a, info) ! local variables integer :: ictxt,np,me - integer :: int_err(5), i1sz,nrow,ncol, dectype, err_act + integer :: int_err(5), i1sz,nrow,ncol, err_act logical, parameter :: debug=.false. character(len=20) :: name,ch_err @@ -192,7 +190,6 @@ subroutine psb_zasbv(x, desc_a, info) name = 'psb_zasbv' ictxt=psb_cd_get_context(desc_a) - dectype=psb_cd_get_dectype(desc_a) call psb_info(ictxt, me, np) diff --git a/src/tools/psb_zcdovr.f90 b/src/tools/psb_zcdovr.f90 index 35f2c028..8eeb1628 100644 --- a/src/tools/psb_zcdovr.f90 +++ b/src/tools/psb_zcdovr.f90 @@ -51,6 +51,10 @@ Subroutine psb_zcdovr(a,desc_a,novr,desc_ov,info) Use psb_prec_mod use psb_error_mod use psb_penv_mod + use psb_tools_mod + use psb_realloc_mod + use psi_mod + use mpi Implicit None ! .. Array Arguments .. @@ -60,42 +64,25 @@ Subroutine psb_zcdovr(a,desc_a,novr,desc_ov,info) Type(psb_desc_type), Intent(inout) :: desc_ov integer, intent(out) :: info - real(kind(1.d0)) :: t1,t2,t3,mpi_wtime - external mpi_wtime integer icomm, err_act - interface psb_cdcpy - subroutine psb_cdcpy(desc_in,desc_out,info) - use psb_descriptor_type - type(psb_desc_type), intent(in) :: desc_in - type(psb_desc_type), intent(out) :: desc_out - integer, intent(out) :: info - end subroutine psb_cdcpy - end interface - - interface psb_cdovrbld - subroutine psb_zcdovrbld(n_ovr,desc_p,desc_a,a,l_tmp_halo,& - & l_tmp_ovr_idx,lworks,lworkr,info) - use psb_prec_type - use psb_spmat_type - type(psb_zspmat_type),intent(in) :: a - type(psb_desc_type),intent(in) :: desc_a - type(psb_desc_type),intent(inout) :: desc_p - integer,intent(in) :: n_ovr - integer, intent(in) :: l_tmp_halo,l_tmp_ovr_idx - integer, intent(inout) :: lworks, lworkr - integer, intent(out) :: info - end subroutine psb_zcdovrbld - end interface - - - ! .. Local Scalars .. Integer :: i, j, k, np, me,m,nnzero,& - & ictxt, lovr, lworks,lworkr, n_col, int_err(5),& + & ictxt, lovr, lworks,lworkr, n_row,n_col, int_err(5),& & index_dim,elem_dim, l_tmp_ovr_idx,l_tmp_halo, nztot,nhalo - Logical, parameter :: debug=.false. - character(len=20) :: name, ch_err + Integer :: counter,counter_h, counter_o, counter_e,idx,gidx,proc,n_elem_recv,& + & n_elem_send,tot_recv,tot_elem,& + & counter_t,n_elem,i_ovr,jj,proc_id,isz, mglob, glx, & + & idxr, idxs, lx, iszr, iszs, nxch, nsnd, nrcv,lidx + + type(psb_zspmat_type) :: blk + Integer, allocatable :: tmp_halo(:),tmp_ovr_idx(:) + Integer,allocatable :: halo(:),works(:),workr(:),t_halo_in(:),& + & t_halo_out(:),temp(:),maskr(:) + Integer,allocatable :: brvindx(:),rvsz(:), bsdindx(:),sdsz(:) + + Logical,Parameter :: debug=.false. + character(len=20) :: name, ch_err name='psb_cdovr' info = 0 @@ -107,10 +94,12 @@ Subroutine psb_zcdovr(a,desc_a,novr,desc_ov,info) If(debug) Write(0,*)'in psb_cdovr',novr - m=psb_cd_get_local_rows(desc_a) - nnzero=Size(a%aspk) - n_col=psb_cd_get_local_cols(desc_a) - nhalo = n_col-m + m = psb_cd_get_local_rows(desc_a) + nnzero = Size(a%aspk) + n_row = psb_cd_get_local_rows(desc_a) + n_col = psb_cd_get_local_cols(desc_a) + nhalo = n_col-m + If(debug) Write(0,*)'IN CDOVR1',novr ,m,nnzero,n_col if (novr<0) then info=10 @@ -120,19 +109,19 @@ Subroutine psb_zcdovr(a,desc_a,novr,desc_ov,info) goto 9999 endif + if (debug) write(0,*) 'Calling desccpy' + call psb_cdcpy(desc_a,desc_ov,info) + if (info /= 0) then + info=4010 + ch_err='psb_cdcpy' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + if (debug) write(0,*) 'From desccpy' if (novr==0) then ! ! Just copy the input. ! - if (debug) write(0,*) 'Calling desccpy' - call psb_cdcpy(desc_a,desc_ov,info) - if (info /= 0) then - info=4010 - ch_err='psb_cdcpy' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - if (debug) write(0,*) 'From desccpy' return endif @@ -142,12 +131,10 @@ Subroutine psb_zcdovr(a,desc_a,novr,desc_ov,info) Write(0,*)'BEGIN cdovr',me,nhalo call psb_barrier(ictxt) endif - t1 = mpi_wtime() - ! ! Ok, since we are only estimating, do it as follows: - ! LOVR= (NNZ/NROW)*N_HALO*N_OVR This assumes that the local average + ! LOVR= (NNZ/NROW)*N_HALO*NOVR This assumes that the local average ! nonzeros per row is the same as the global. ! nztot = psb_sp_get_nnzeros(a) @@ -156,59 +143,557 @@ Subroutine psb_zcdovr(a,desc_a,novr,desc_ov,info) lworks = ((nztot+m-1)/m)*nhalo lworkr = ((nztot+m-1)/m)*nhalo else - info=-1 - call psb_errpush(info,name) - goto 9999 + info=-1 + call psb_errpush(info,name) + goto 9999 endif If(debug)Write(0,*)'ovr_est done',me,novr,lovr + index_dim = size(desc_a%halo_index) elem_dim = size(desc_a%halo_index) - call psb_realloc(psb_mdata_size_,desc_ov%matrix_data,info) - if (info==0) call psb_realloc(novr*(Max(elem_dim,1)+3),desc_ov%ovrlap_elem,info) - if (info /= 0) then - info=4000 - call psb_errpush(info,name) - goto 9999 - end if - - l_tmp_ovr_idx=novr*(3*Max(2*index_dim,1)+1) - l_tmp_halo=novr*(3*Size(desc_a%halo_index)) - desc_ov%matrix_data(:) = desc_a%matrix_data(:) - desc_ov%matrix_data(psb_dec_type_) = psb_desc_bld_ + l_tmp_ovr_idx = novr*(3*Max(2*index_dim,1)+1) + l_tmp_halo = novr*(3*Size(desc_a%halo_index)) - Allocate(desc_ov%loc_to_glob(Size(desc_a%loc_to_glob)),& - & desc_ov%glob_to_loc(Size(desc_a%glob_to_loc)),stat=info) - if (info /= 0) then - info=4000 - call psb_errpush(info,name) - goto 9999 + if (psb_is_large_desc(desc_a)) then + desc_ov%matrix_data(psb_dec_type_) = psb_desc_large_bld_ + else + desc_ov%matrix_data(psb_dec_type_) = psb_desc_bld_ end if - - desc_ov%loc_to_glob(:) = desc_a%loc_to_glob(:) - desc_ov%glob_to_loc(:) = desc_a%glob_to_loc(:) If(debug) then Write(0,*)'Start cdovrbld',me,lworks,lworkr call psb_barrier(ictxt) endif - ! - ! The real work goes on in here.... - ! - Call psb_cdovrbld(novr,desc_ov,desc_a,a,& - & l_tmp_halo,l_tmp_ovr_idx,lworks,lworkr,info) - if (info /= 0) then - info=4010 - ch_err='psb_cdovrbld' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 + if (.false.) then + ! + ! The real work goes on in here.... + ! + Call psb_cdovrbld(novr,desc_ov,desc_a,a,& + & l_tmp_halo,l_tmp_ovr_idx,lworks,lworkr,info) + + if (info /= 0) then + info=4010 + ch_err='psb_cdovrbld' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + + If(debug) then + Write(0,*)'Done cdovrbld',me,lworks,lworkr + call psb_barrier(ictxt) + endif + + else + + + Allocate(brvindx(np+1),rvsz(np),sdsz(np),bsdindx(np+1),stat=info) + if (info /= 0) then + call psb_errpush(4010,name,a_err='Allocate') + goto 9999 + end if + + Allocate(works(lworks),workr(lworkr),t_halo_in(l_tmp_halo),& + & t_halo_out(l_tmp_halo), temp(lworkr),stat=info) + if (info /= 0) then + call psb_errpush(4010,name,a_err='Allocate') + goto 9999 + end if + + + call psb_sp_all(blk,max(lworks,lworkr),info) + if (info /= 0) then + info=4010 + ch_err='psb_sp_all' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + + blk%fida='COO' + + Allocate(tmp_ovr_idx(l_tmp_ovr_idx),tmp_halo(l_tmp_halo),& + & halo(size(desc_a%halo_index)),stat=info) + if (info /= 0) then + call psb_errpush(4010,name,a_err='Allocate') + goto 9999 + end if + halo(:) = desc_a%halo_index(:) + desc_ov%ovrlap_elem(:) = -1 + tmp_ovr_idx(:) = -1 + tmp_halo(:) = -1 + counter_e = 1 + tot_recv = 0 + counter_h = 1 + counter_o = 1 + + ! Init overlap with desc_a%ovrlap (if any) + counter = 1 + Do While (desc_a%ovrlap_index(counter) /= -1) + proc = desc_a%ovrlap_index(counter+psb_proc_id_) + n_elem_recv = desc_a%ovrlap_index(counter+psb_n_elem_recv_) + n_elem_send = desc_a%ovrlap_index(counter+n_elem_recv+psb_n_elem_send_) + + Do j=0,n_elem_recv-1 + + idx = desc_a%ovrlap_index(counter+psb_elem_recv_+j) + If(idx > Size(desc_ov%loc_to_glob)) then + info=-3 + call psb_errpush(info,name) + goto 9999 + endif + + gidx = desc_ov%loc_to_glob(idx) + + call psb_check_size((counter_o+3),tmp_ovr_idx,info,pad=-1) + if (info /= 0) then + info=4010 + call psb_errpush(info,name,a_err='psb_check_size') + goto 9999 + end if + + tmp_ovr_idx(counter_o)=proc + tmp_ovr_idx(counter_o+1)=1 + tmp_ovr_idx(counter_o+2)=gidx + tmp_ovr_idx(counter_o+3)=-1 + counter_o=counter_o+3 + end Do + counter=counter+n_elem_recv+n_elem_send+2 + end Do + + + + + ! + ! A picture is in order to understand what goes on here. + ! I is the internal part; H is halo, R row, C column. The final + ! matrix with N levels of overlap looks like this + ! + ! I | Hc1 | 0 | 0 | + ! -------|-----|-----|-----| + ! Hr1 | Hd1 | Hc2 | 0 | + ! -------|-----|-----|-----| + ! 0 | Hr2 | Hd2 | Hc2 | + ! -------|-----|-----|-----| + ! 0 | 0 | Hr3 | Hd3 | Hc3 + ! + ! At the start we already have I and Hc1, so we know the row + ! indices that will make up Hr1, and also who owns them. As we + ! actually get those rows, we receive the column indices in Hc2; + ! these define the row indices for Hr2, and so on. When we have + ! reached the desired level HrN, we may ignore HcN. + ! + ! + Do i_ovr = 1, novr + + if (debug) write(0,*) me,'Running on overlap level ',i_ovr,' of ',novr + + ! + ! At this point, halo contains a valid halo corresponding to the + ! matrix enlarged with the elements in the frontier for I_OVR-1. + ! At the start, this is just the halo for A; the rows for indices in + ! the first halo will contain column indices defining the second halo + ! level and so on. + ! + bsdindx(:) = 0 + sdsz(:) = 0 + brvindx(:) = 0 + rvsz(:) = 0 + idxr = 0 + idxs = 0 + counter = 1 + counter_t = 1 + + + Do While (halo(counter) /= -1) + tot_elem=0 + proc=halo(counter+psb_proc_id_) + n_elem_recv=halo(counter+psb_n_elem_recv_) + n_elem_send=halo(counter+n_elem_recv+psb_n_elem_send_) + If ((counter+n_elem_recv+n_elem_send) > Size(halo)) then + info = -1 + call psb_errpush(info,name) + goto 9999 + end If + tot_recv=tot_recv+n_elem_recv + if (debug) write(0,*) me,' CDOVRBLD tot_recv:',proc,n_elem_recv,tot_recv + ! + ! + ! The format of the halo vector exists in two forms: 1. Temporary + ! 2. Assembled. In this loop we are using the (assembled) halo_in and + ! copying it into (temporary) halo_out; this is because tmp_halo will + ! be enlarged with the new column indices received, and will reassemble + ! everything for the next iteration. + ! + + ! + ! add recv elements in halo_index into ovrlap_index + ! + Do j=0,n_elem_recv-1 + If((counter+psb_elem_recv_+j)>Size(halo)) then + info=-2 + call psb_errpush(info,name) + goto 9999 + end If + + idx = halo(counter+psb_elem_recv_+j) + If(idx > Size(desc_ov%loc_to_glob)) then + info=-3 + call psb_errpush(info,name) + goto 9999 + endif + + gidx = desc_ov%loc_to_glob(idx) + + call psb_check_size((counter_o+3),tmp_ovr_idx,info,pad=-1) + if (info /= 0) then + info=4010 + call psb_errpush(info,name,a_err='psb_check_size') + goto 9999 + end if + + tmp_ovr_idx(counter_o)=proc + tmp_ovr_idx(counter_o+1)=1 + tmp_ovr_idx(counter_o+2)=gidx + tmp_ovr_idx(counter_o+3)=-1 + counter_o=counter_o+3 + if (.not.psb_is_large_desc(desc_ov)) then + call psb_check_size((counter_h+3),tmp_halo,info,pad=-1) + if (info /= 0) then + info=4010 + call psb_errpush(info,name,a_err='psb_check_size') + goto 9999 + end if + + tmp_halo(counter_h)=proc + tmp_halo(counter_h+1)=1 + tmp_halo(counter_h+2)=idx + tmp_halo(counter_h+3)=-1 + + counter_h=counter_h+3 + end if + + Enddo + if (debug) write(0,*) me,'Checktmp_o_i Loop Mid1',tmp_ovr_idx(1:10) + counter = counter+n_elem_recv + + ! + ! add send elements in halo_index into ovrlap_index + ! + Do j=0,n_elem_send-1 + + idx = halo(counter+psb_elem_send_+j) + gidx = desc_ov%loc_to_glob(idx) + if (idx > psb_cd_get_local_rows(Desc_a)) & + & write(0,*) me,i_ovr,'Out of local rows ',& + & idx,psb_cd_get_local_rows(Desc_a) + + call psb_check_size((counter_o+3),tmp_ovr_idx,info,pad=-1) + if (info /= 0) then + info=4010 + call psb_errpush(info,name,a_err='psb_check_size') + goto 9999 + end if + + tmp_ovr_idx(counter_o)=proc + tmp_ovr_idx(counter_o+1)=1 + tmp_ovr_idx(counter_o+2)=gidx + tmp_ovr_idx(counter_o+3)=-1 + counter_o=counter_o+3 + + ! + ! Prepare to exchange the halo rows with the other proc. + ! + If (i_ovr < (novr)) Then + n_elem = psb_sp_get_nnz_row(idx,a) + + call psb_check_size((idxs+tot_elem+n_elem),works,info) + if (info /= 0) then + info=4010 + call psb_errpush(info,name,a_err='psb_check_size') + goto 9999 + end if + + If((n_elem) > size(blk%ia2)) Then + isz = max((3*size(blk%ia2))/2,(n_elem)) + if (debug) write(0,*) me,'Realloc blk',isz + call psb_sp_reall(blk,isz,info) + if (info /= 0) then + info=4010 + ch_err='psb_sp_reall' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + End If + + call psb_sp_getblk(idx,a,blk,info) + if (info /= 0) then + info=4010 + ch_err='psb_sp_getblk' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if +!!$ write(0,*) me,'Iteration: ',j,i_ovr + Do jj=1,n_elem + works(idxs+tot_elem+jj)=desc_ov%loc_to_glob(blk%ia2(jj)) + End Do + tot_elem=tot_elem+n_elem + End If + + Enddo + + + if (i_ovr < novr) then + if (tot_elem > 1) then + call imsr(tot_elem,works(idxs+1)) + lx = works(idxs+1) + i = 1 + j = 1 + do + j = j + 1 + if (j > tot_elem) exit + if (works(idxs+j) /= lx) then + i = i + 1 + works(idxs+i) = works(idxs+j) + lx = works(idxs+i) + end if + end do + tot_elem = i + endif + if (debug) write(0,*) me,'Checktmp_o_i Loop Mid2',tmp_ovr_idx(1:10) + sdsz(proc+1) = tot_elem + idxs = idxs + tot_elem + end if + counter = counter+n_elem_send+3 + if (debug) write(0,*) me,'Checktmp_o_i Loop End',tmp_ovr_idx(1:10) + Enddo + if (debug) write(0,*)me,'End phase 1 CDOVRBLD', m, n_col, tot_recv + + if (i_ovr < novr) then + ! + ! Exchange data requests with everybody else: so far we have + ! accumulated RECV requests, we have an all-to-all to build + ! matchings SENDs. + ! + call mpi_alltoall(sdsz,1,mpi_integer,rvsz,1,mpi_integer,icomm,info) + if (info /= 0) then + info=4010 + ch_err='mpi_alltoall' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + idxs = 0 + idxr = 0 + counter = 1 + Do + proc=halo(counter) + if (proc == -1) exit + n_elem_recv = halo(counter+psb_n_elem_recv_) + counter = counter+n_elem_recv + n_elem_send = halo(counter+psb_n_elem_send_) + + bsdindx(proc+1) = idxs + idxs = idxs + sdsz(proc+1) + brvindx(proc+1) = idxr + idxr = idxr + rvsz(proc+1) + counter = counter+n_elem_send+3 + Enddo + + iszr=sum(rvsz) + if (max(iszr,1) > lworkr) then + call psb_realloc(max(iszr,1),workr,info) + if (info /= 0) then + info=4010 + ch_err='psb_realloc' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + lworkr=max(iszr,1) + end if + + call mpi_alltoallv(works,sdsz,bsdindx,mpi_integer,& + & workr,rvsz,brvindx,mpi_integer,icomm,info) + if (info /= 0) then + info=4010 + ch_err='mpi_alltoallv' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + + if (debug) write(0,*) 'ISZR :',iszr + + if (psb_is_large_desc(desc_a)) then + call psb_check_size(iszr,maskr,info) + if (info /= 0) then + info=4010 + call psb_errpush(info,name,a_err='psb_check_size') + goto 9999 + end if + call psi_idx_cnv(iszr,workr,maskr,desc_ov,info) + iszs = count(maskr(1:iszr)<=0) + if (iszs > size(works)) call psb_realloc(iszs,works,info) + j = 0 + do i=1,iszr + if (maskr(i) < 0) then + j=j+1 + works(j) = workr(i) + end if + end do + ! + ! fnd_owner on desc_a because we want the procs who + ! owned the rows from the beginning! + ! + call psi_fnd_owner(iszs,works,temp,desc_a,info) + n_col=psb_cd_get_local_cols(desc_ov) + + do i=1,iszs + idx = works(i) + n_col = psb_cd_get_local_cols(desc_ov) + call psi_idx_ins_cnv(idx,lidx,desc_ov,info) + if (psb_cd_get_local_cols(desc_ov) > n_col ) then + ! + ! This is a new index. Assigning a local index as + ! we receive them guarantees that all indices for HALO(I) + ! will be less than those for HALO(J) whenever I 0 thanks to the caller routine. ! ! Parameters: -! n_ovr - integer. The number of overlap levels +! n_ovr - integer. The number of overlap levels ! desc_p - type(). The communication descriptor ! for the preconditioner. ! desc_a - type(). The communication descriptor. @@ -54,34 +54,30 @@ Subroutine psb_zcdovrbld(n_ovr,desc_p,desc_a,a,& use psb_serial_mod Use psi_mod use psb_realloc_mod - use psb_tools_mod, only : psb_cdprt use psb_error_mod use psb_const_mod use psb_penv_mod + use psb_tools_mod + use mpi Implicit None - include 'mpif.h' type(psb_zspmat_type),intent(in) :: a type(psb_desc_type),intent(in) :: desc_a type(psb_desc_type),intent(inout) :: desc_p - integer,intent(in) :: n_ovr + integer,intent(in) :: n_ovr ! Input estimates for allocation sizes. Could we do without these two??? - Integer, Intent(in) :: l_tmp_halo,l_tmp_ovr_idx + Integer, Intent(in) :: l_tmp_halo,l_tmp_ovr_idx Integer, Intent(inout) :: lworks, lworkr integer, intent(out) :: info - type(psb_zspmat_type) :: blk - - + type(psb_zspmat_type) :: blk Integer, allocatable :: tmp_halo(:),tmp_ovr_idx(:) - Integer :: counter,counter_h, counter_o, counter_e,j,idx,gidx,proc,n_elem_recv,& & n_elem_send,tot_recv,tot_elem,n_col,m,ictxt,np,me,dl_lda,lwork,& & counter_t,n_elem,i_ovr,jj,i,proc_id,isz, mglob, glx,n_row, & - & idxr, idxs, lx, iszr, err_act, icomm, nxch, nsnd, nrcv - - Integer,allocatable :: halo(:),length_dl(:),works(:),workr(:),t_halo_in(:),& - & t_halo_out(:),work(:),dep_list(:),temp(:) + & idxr, idxs, lx, iszr, iszs, err_act, icomm, nxch, nsnd, nrcv,lidx + Integer,allocatable :: halo(:),works(:),workr(:),t_halo_in(:),& + & t_halo_out(:),temp(:),maskr(:) Integer,allocatable :: brvindx(:),rvsz(:), bsdindx(:),sdsz(:) Logical,Parameter :: debug=.false. @@ -114,11 +110,8 @@ Subroutine psb_zcdovrbld(n_ovr,desc_p,desc_a,a,& n_col = psb_cd_get_local_cols(desc_a) if (debug) write(0,*) me,' On entry to CDOVRBLD n_col:',n_col - dl_lda=np*5 - lwork=5*(5*np+2)*np+10 Allocate(works(lworks),workr(lworkr),t_halo_in(l_tmp_halo),& - & t_halo_out(l_tmp_halo), work(lwork),& - & length_dl(np+1),dep_list(dl_lda*np),temp(lworkr),stat=info) + & t_halo_out(l_tmp_halo), temp(lworkr),stat=info) if (info /= 0) then call psb_errpush(4010,name,a_err='Allocate') goto 9999 @@ -141,7 +134,7 @@ Subroutine psb_zcdovrbld(n_ovr,desc_p,desc_a,a,& call psb_errpush(4010,name,a_err='Allocate') goto 9999 end if - halo(:) = desc_a%halo_index(:) + halo(:) = desc_a%halo_index(:) desc_p%ovrlap_elem(:) = -1 tmp_ovr_idx(:) = -1 tmp_halo(:) = -1 @@ -150,7 +143,40 @@ Subroutine psb_zcdovrbld(n_ovr,desc_p,desc_a,a,& counter_h = 1 counter_o = 1 - ! See comment in main loop below. + ! Init overlap with desc_a%ovrlap (if any) + counter = 1 + Do While (desc_a%ovrlap_index(counter) /= -1) + proc = desc_a%ovrlap_index(counter+psb_proc_id_) + n_elem_recv = desc_a%ovrlap_index(counter+psb_n_elem_recv_) + n_elem_send = desc_a%ovrlap_index(counter+n_elem_recv+psb_n_elem_send_) + + Do j=0,n_elem_recv-1 + + idx = desc_a%ovrlap_index(counter+psb_elem_recv_+j) + If(idx > Size(desc_p%loc_to_glob)) then + info=-3 + call psb_errpush(info,name) + goto 9999 + endif + + gidx = desc_p%loc_to_glob(idx) + + call psb_check_size((counter_o+3),tmp_ovr_idx,info,pad=-1) + if (info /= 0) then + info=4010 + call psb_errpush(info,name,a_err='psb_check_size') + goto 9999 + end if + + tmp_ovr_idx(counter_o)=proc + tmp_ovr_idx(counter_o+1)=1 + tmp_ovr_idx(counter_o+2)=gidx + tmp_ovr_idx(counter_o+3)=-1 + counter_o=counter_o+3 + end Do + counter=counter+n_elem_recv+n_elem_send+2 + end Do + ! ! A picture is in order to understand what goes on here. @@ -175,7 +201,6 @@ Subroutine psb_zcdovrbld(n_ovr,desc_p,desc_a,a,& Do i_ovr=1,n_ovr if (debug) write(0,*) me,'Running on overlap level ',i_ovr,' of ',n_ovr -!!$ t_halo_in(:) = -1 ! ! At this point, halo contains a valid halo corresponding to the @@ -192,7 +217,8 @@ Subroutine psb_zcdovrbld(n_ovr,desc_p,desc_a,a,& idxs = 0 counter = 1 counter_t = 1 - + + t1 = mpi_wtime() Do While (halo(counter) /= -1) tot_elem=0 @@ -200,7 +226,7 @@ Subroutine psb_zcdovrbld(n_ovr,desc_p,desc_a,a,& n_elem_recv=halo(counter+psb_n_elem_recv_) n_elem_send=halo(counter+n_elem_recv+psb_n_elem_send_) If ((counter+n_elem_recv+n_elem_send) > Size(halo)) then - info=-1 + info = -1 call psb_errpush(info,name) goto 9999 end If @@ -210,7 +236,7 @@ Subroutine psb_zcdovrbld(n_ovr,desc_p,desc_a,a,& ! ! The format of the halo vector exists in two forms: 1. Temporary ! 2. Assembled. In this loop we are using the (assembled) halo_in and - ! copying it into (temporary) tmp_halo; this is because tmp_halo will + ! copying it into (temporary) halo_out; this is because tmp_halo will ! be enlarged with the new column indices received, and will reassemble ! everything for the next iteration. ! @@ -225,7 +251,7 @@ Subroutine psb_zcdovrbld(n_ovr,desc_p,desc_a,a,& call psb_errpush(info,name) goto 9999 end If - idx = halo(counter+psb_elem_recv_+j) + idx = halo(counter+psb_elem_recv_+j) If(idx > Size(desc_p%loc_to_glob)) then info=-3 @@ -235,42 +261,33 @@ Subroutine psb_zcdovrbld(n_ovr,desc_p,desc_a,a,& gidx = desc_p%loc_to_glob(idx) - If((counter_o+2) > Size(tmp_ovr_idx)) Then - isz = max((3*Size(tmp_ovr_idx))/2,(counter_o+3)) - if (debug) write(0,*) me,'Realloc tmp_ovr',isz - call psb_realloc(isz,tmp_ovr_idx,info,pad=-1) - if (info /= 0) then - info=4010 - ch_err='psb_realloc' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - End If + call psb_check_size((counter_o+3),tmp_ovr_idx,info,pad=-1) + if (info /= 0) then + info=4010 + call psb_errpush(info,name,a_err='psb_check_size') + goto 9999 + end if tmp_ovr_idx(counter_o)=proc tmp_ovr_idx(counter_o+1)=1 tmp_ovr_idx(counter_o+2)=gidx tmp_ovr_idx(counter_o+3)=-1 counter_o=counter_o+3 - - If((counter_h+2) > Size(tmp_halo)) Then - isz = max((3*Size(tmp_halo))/2,(counter_h+3)) - if (debug) write(0,*) me,'Realloc tmp_halo',isz - call psb_realloc(isz,tmp_halo,info) + if (.not.psb_is_large_desc(desc_p)) then + call psb_check_size((counter_h+3),tmp_halo,info,pad=-1) if (info /= 0) then info=4010 - ch_err='psb_realloc' - call psb_errpush(info,name,a_err=ch_err) + call psb_errpush(info,name,a_err='psb_check_size') goto 9999 end if - End If - - tmp_halo(counter_h)=proc - tmp_halo(counter_h+1)=1 - tmp_halo(counter_h+2)=idx - tmp_halo(counter_h+3)=-1 - - counter_h=counter_h+3 + + tmp_halo(counter_h)=proc + tmp_halo(counter_h+1)=1 + tmp_halo(counter_h+2)=idx + tmp_halo(counter_h+3)=-1 + + counter_h=counter_h+3 + end if Enddo if (debug) write(0,*) me,'Checktmp_o_i Loop Mid1',tmp_ovr_idx(1:10) @@ -283,18 +300,15 @@ Subroutine psb_zcdovrbld(n_ovr,desc_p,desc_a,a,& idx = halo(counter+psb_elem_send_+j) gidx = desc_p%loc_to_glob(idx) - - If((counter_o+2) > Size(tmp_ovr_idx)) Then - isz = max((3*Size(tmp_ovr_idx))/2,(counter_o+3)) - if (debug) write(0,*) me,'Realloc tmp_ovr',isz - call psb_realloc(isz,tmp_ovr_idx,info) - if (info /= 0) then - info=4010 - ch_err='psrealloc' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - End If + if (idx > psb_cd_get_local_rows(Desc_a)) & + & write(0,*) me,i_ovr,'Out of local rows ',idx,psb_cd_get_local_rows(Desc_a) + + call psb_check_size((counter_o+3),tmp_ovr_idx,info,pad=-1) + if (info /= 0) then + info=4010 + call psb_errpush(info,name,a_err='psb_check_size') + goto 9999 + end if tmp_ovr_idx(counter_o)=proc tmp_ovr_idx(counter_o+1)=1 @@ -308,18 +322,12 @@ Subroutine psb_zcdovrbld(n_ovr,desc_p,desc_a,a,& If (i_ovr < (n_ovr)) Then n_elem = psb_sp_get_nnz_row(idx,a) - If((idxs+tot_elem+n_elem) > lworks) Then - isz = max((3*lworks)/2,(idxs+tot_elem+n_elem)) - if (debug) write(0,*) me,'Realloc works',isz - call psb_realloc(isz,works,info) - if (info /= 0) then - info=4010 - ch_err='psb_realloc' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - lworks = isz - End If + call psb_check_size((idxs+tot_elem+n_elem),works,info) + if (info /= 0) then + info=4010 + call psb_errpush(info,name,a_err='psb_check_size') + goto 9999 + end if If((n_elem) > size(blk%ia2)) Then isz = max((3*size(blk%ia2))/2,(n_elem)) @@ -428,61 +436,109 @@ Subroutine psb_zcdovrbld(n_ovr,desc_p,desc_a,a,& goto 9999 end if - Do i=1,iszr - idx=workr(i) - if (idx <1) then - write(0,*) me,'Error in CDOVRBLD ',idx,i,iszr -!!$ write(0,*) me, ' WORKR :',workr(1:iszr) - else If (desc_p%glob_to_loc(idx) < -np) Then - ! - ! This is a new index. Assigning a local index as - ! we receive the guarantees that all indices for HALO(I) - ! will be less than those for HALO(J) whenever I Size(desc_p%loc_to_glob)) Then - isz = 3*n_col/2 - if (debug) write(0,*) me,'Realloc loc_to_glob' - call psb_realloc(isz,desc_p%loc_to_glob,info) + if (debug) write(0,*) 'ISZR :',iszr + + if (psb_is_large_desc(desc_a)) then + call psb_check_size(iszr,maskr,info) + if (info /= 0) then + info=4010 + call psb_errpush(info,name,a_err='psb_check_size') + goto 9999 + end if + call psi_idx_cnv(iszr,workr,maskr,desc_p,info) + iszs = count(maskr(1:iszr)<=0) + if (iszs > size(works)) call psb_realloc(iszs,works,info) + j = 0 + do i=1,iszr + if (maskr(i) < 0) then + j=j+1 + works(j) = workr(i) + end if + end do + ! + ! fnd_owner on desc_a because we want the procs who + ! owned the rows from the beginning! + ! + call psi_fnd_owner(iszs,works,temp,desc_a,info) + n_col=psb_cd_get_local_cols(desc_p) + + do i=1,iszs + idx = works(i) + n_col = psb_cd_get_local_cols(desc_p) + call psi_idx_ins_cnv(idx,lidx,desc_p,info) + if (psb_cd_get_local_cols(desc_p) > n_col ) then + ! + ! This is a new index. Assigning a local index as + ! we receive them guarantees that all indices for HALO(I) + ! will be less than those for HALO(J) whenever I Size(t_halo_in)) then - isz = max((3*Size(t_halo_in))/2,(counter_t+3+1000)) - if (debug) write(0,*) me,'Realloc ovr_el',isz - call psb_realloc(isz,t_halo_in,info) + + t_halo_in(counter_t)=proc_id + t_halo_in(counter_t+1)=1 + t_halo_in(counter_t+2)=lidx + counter_t=counter_t+3 + if (.false.) write(0,*) me,' CDOVRBLD: Added t_halo_in ',& + &proc_id,lidx,idx + endif + end Do + + else + + Do i=1,iszr + idx=workr(i) + if (idx <1) then + write(0,*) me,'Error in CDOVRBLD level',i_ovr,' : ',idx,i,iszr + else If (desc_p%glob_to_loc(idx) < -np) Then + ! + ! This is a new index. Assigning a local index as + ! we receive them guarantees that all indices for HALO(I) + ! will be less than those for HALO(J) whenever I). The communication descriptor. ! info - integer. Eventually returns an error code subroutine psb_zinsvi(m, irw, val, x, desc_a, info, dupl) @@ -45,6 +45,7 @@ subroutine psb_zinsvi(m, irw, val, x, desc_a, info, dupl) use psb_const_mod use psb_error_mod use psb_penv_mod + use psi_mod implicit none ! m rows number of submatrix belonging to val to be inserted @@ -62,21 +63,22 @@ subroutine psb_zinsvi(m, irw, val, x, desc_a, info, dupl) integer, optional, intent(in) :: dupl !locals..... - integer :: ictxt,i,loc_row,glob_row,row,& + integer :: ictxt,i,loc_row,glob_row,& & loc_rows,loc_cols,mglob,err_act, int_err(5) - integer :: np,me,dupl_ - character(len=20) :: name + integer :: np, me, dupl_ + integer, allocatable :: irl(:) + character(len=20) :: name if(psb_get_errstatus() /= 0) return info=0 call psb_erractionsave(err_act) - name = 'psb_dinsvv' + name = 'psb_zinsvi' - if (.not.allocated(desc_a%glob_to_loc)) then - info=3110 - call psb_errpush(info,name) - return - end if +!!$ if (.not.allocated(desc_a%glob_to_loc)) then +!!$ info=3110 +!!$ call psb_errpush(info,name) +!!$ return +!!$ end if if ((.not.allocated(desc_a%matrix_data))) then int_err(1)=3110 call psb_errpush(info,name) @@ -112,15 +114,25 @@ subroutine psb_zinsvi(m, irw, val, x, desc_a, info, dupl) goto 9999 endif - loc_rows=psb_cd_get_local_rows(desc_a) - loc_cols=psb_cd_get_local_cols(desc_a) + if (m==0) return + loc_rows = psb_cd_get_local_rows(desc_a) + loc_cols = psb_cd_get_local_cols(desc_a) mglob = psb_cd_get_global_rows(desc_a) + allocate(irl(m),stat=info) + if (info /= 0) then + info = 4000 + call psb_errpush(info,name) + goto 9999 + endif + if (present(dupl)) then dupl_ = dupl else dupl_ = psb_dupl_ovwrt_ endif + + call psi_idx_cnv(m,irw,irl,desc_a,info,owned=.true.) select case(dupl_) case(psb_dupl_ovwrt_) @@ -128,15 +140,10 @@ subroutine psb_zinsvi(m, irw, val, x, desc_a, info, dupl) !loop over all val's rows ! row actual block row - glob_row=irw(i) - if ((glob_row>0).and.(glob_row <= mglob)) then - - loc_row=desc_a%glob_to_loc(glob_row) - if (loc_row.ge.1) then - ! this row belongs to me - ! copy i-th row of block val in x - x(loc_row) = val(i) - end if + if (irl(i) > 0) then + ! this row belongs to me + ! copy i-th row of block val in x + x(irl(i)) = val(i) end if enddo @@ -145,16 +152,10 @@ subroutine psb_zinsvi(m, irw, val, x, desc_a, info, dupl) do i = 1, m !loop over all val's rows - ! row actual block row - glob_row=irw(i) - if ((glob_row>0).and.(glob_row <= mglob)) then - - loc_row=desc_a%glob_to_loc(glob_row) - if (loc_row.ge.1) then + if (irl(i) > 0) then ! this row belongs to me ! copy i-th row of block val in x - x(loc_row) = x(loc_row) + val(i) - end if + x(irl(i)) = x(irl(i)) + val(i) end if enddo @@ -163,6 +164,7 @@ subroutine psb_zinsvi(m, irw, val, x, desc_a, info, dupl) call psb_errpush(info,name) goto 9999 end select + deallocate(irl) call psb_erractionrestore(err_act) return @@ -217,7 +219,7 @@ end subroutine psb_zinsvi ! val to be inserted. ! irw - integer(:) Row indices of rows of val (global numbering) ! val - complex, dimension(:,:). The source dense submatrix. -! x - complex, pointer, dimension(:,:). The destination dense matrix. +! x - complex, dimension(:,:). The destination dense matrix. ! desc_a - type(). The communication descriptor. ! info - integer. Eventually returns an error code subroutine psb_zinsi(m,irw, val, x, desc_a, info, dupl) @@ -226,6 +228,7 @@ subroutine psb_zinsi(m,irw, val, x, desc_a, info, dupl) use psb_const_mod use psb_error_mod use psb_penv_mod + use psi_mod implicit none ! m rows number of submatrix belonging to val to be inserted @@ -246,18 +249,19 @@ subroutine psb_zinsi(m,irw, val, x, desc_a, info, dupl) integer :: ictxt,i,loc_row,glob_row,j,n,& & loc_rows,loc_cols,mglob,err_act, int_err(5) integer :: np,me,dupl_ + integer, allocatable :: irl(:) character(len=20) :: name if(psb_get_errstatus() /= 0) return info=0 call psb_erractionsave(err_act) - name = 'psb_dinsvv' + name = 'psb_zinsi' - if (.not.allocated(desc_a%glob_to_loc)) then - info=3110 - call psb_errpush(info,name) - return - end if +!!$ if (.not.allocated(desc_a%glob_to_loc)) then +!!$ info=3110 +!!$ call psb_errpush(info,name) +!!$ return +!!$ end if if ((.not.allocated(desc_a%matrix_data))) then int_err(1)=3110 call psb_errpush(info,name) @@ -274,7 +278,7 @@ subroutine psb_zinsi(m,irw, val, x, desc_a, info, dupl) endif !... check parameters.... - if (m.lt.0) then + if (m < 0) then info = 10 int_err(1) = 1 int_err(2) = m @@ -292,9 +296,10 @@ subroutine psb_zinsi(m,irw, val, x, desc_a, info, dupl) call psb_errpush(info,name,int_err) goto 9999 endif + if (m==0) return - loc_rows=psb_cd_get_local_rows(desc_a) - loc_cols=psb_cd_get_local_cols(desc_a) + loc_rows = psb_cd_get_local_rows(desc_a) + loc_cols = psb_cd_get_local_cols(desc_a) mglob = psb_cd_get_global_rows(desc_a) n = min(size(val,2),size(x,2)) @@ -305,23 +310,28 @@ subroutine psb_zinsi(m,irw, val, x, desc_a, info, dupl) dupl_ = psb_dupl_ovwrt_ endif + allocate(irl(m),stat=info) + if (info /= 0) then + info = 4000 + call psb_errpush(info,name) + goto 9999 + endif + + call psi_idx_cnv(m,irw,irl,desc_a,info,owned=.true.) + select case(dupl_) case(psb_dupl_ovwrt_) do i = 1, m !loop over all val's rows ! row actual block row - glob_row=irw(i) - if ((glob_row>0).and.(glob_row <= mglob)) then - - loc_row=desc_a%glob_to_loc(glob_row) - if (loc_row.ge.1) then - ! this row belongs to me - ! copy i-th row of block val in x - do j=1,n - x(loc_row,j) = val(i,j) - end do - end if + loc_row = irl(i) + if (loc_row > 0) then + ! this row belongs to me + ! copy i-th row of block val in x + do j=1,n + x(loc_row,j) = val(i,j) + end do end if enddo @@ -331,17 +341,13 @@ subroutine psb_zinsi(m,irw, val, x, desc_a, info, dupl) !loop over all val's rows ! row actual block row - glob_row=irw(i) - if ((glob_row>0).and.(glob_row <= mglob)) then - - loc_row=desc_a%glob_to_loc(glob_row) - if (loc_row.ge.1) then - ! this row belongs to me - ! copy i-th row of block val in x - do j=1,n - x(loc_row,j) = x(loc_row,j) + val(i,j) - end do - end if + loc_row = irl(i) + if (loc_row > 0) then + ! this row belongs to me + ! copy i-th row of block val in x + do j=1,n + x(loc_row,j) = x(loc_row,j) + val(i,j) + end do end if enddo @@ -350,6 +356,7 @@ subroutine psb_zinsi(m,irw, val, x, desc_a, info, dupl) call psb_errpush(info,name) goto 9999 end select + deallocate(irl) call psb_erractionrestore(err_act) return diff --git a/src/tools/psb_zspalloc.f90 b/src/tools/psb_zspalloc.f90 index 65e0bf15..cfcce5a6 100644 --- a/src/tools/psb_zspalloc.f90 +++ b/src/tools/psb_zspalloc.f90 @@ -56,7 +56,7 @@ subroutine psb_zspalloc(a, desc_a, info, nnz) integer, optional, intent(in) :: nnz !locals - integer :: ictxt, dectype + integer :: ictxt integer :: np,me,loc_row,& & length_ia1,length_ia2, err_act,m,n integer :: int_err(5) @@ -69,7 +69,6 @@ subroutine psb_zspalloc(a, desc_a, info, nnz) name = 'psb_zspalloc' ictxt = psb_cd_get_context(desc_a) - dectype = psb_cd_get_dectype(desc_a) call psb_info(ictxt, me, np) ! ....verify blacs grid correctness.. if (np == -1) then @@ -127,7 +126,7 @@ subroutine psb_zspalloc(a, desc_a, info, nnz) a%infoa(psb_state_) = psb_spmat_bld_ if (debug) write(0,*) 'spall: ', & - &psb_cd_get_dectype(desc_a),psb_desc_bld_ + & psb_cd_get_dectype(desc_a),psb_desc_bld_ return diff --git a/src/tools/psb_zspasb.f90 b/src/tools/psb_zspasb.f90 index 33e7324f..c4d00a5a 100644 --- a/src/tools/psb_zspasb.f90 +++ b/src/tools/psb_zspasb.f90 @@ -65,7 +65,7 @@ subroutine psb_zspasb(a,desc_a, info, afmt, upd, dupl) integer :: int_err(5) type(psb_zspmat_type) :: atemp integer :: np,me,n_col,iout, err_act - integer :: dscstate, spstate + integer :: spstate integer :: upd_, dupl_ integer :: ictxt,n_row logical, parameter :: debug=.false., debugwrt=.false. @@ -76,8 +76,7 @@ subroutine psb_zspasb(a,desc_a, info, afmt, upd, dupl) name = 'psb_spasb' call psb_erractionsave(err_act) - ictxt = psb_cd_get_context(desc_a) - dscstate = psb_cd_get_dectype(desc_a) + ictxt = psb_cd_get_context(desc_a) n_row = psb_cd_get_local_rows(desc_a) n_col = psb_cd_get_local_cols(desc_a) @@ -89,9 +88,9 @@ subroutine psb_zspasb(a,desc_a, info, afmt, upd, dupl) goto 9999 endif - if (.not.psb_is_asb_dec(dscstate)) then + if (.not.psb_is_asb_desc(desc_a)) then info = 600 - int_err(1) = dscstate + int_err(1) = psb_cd_get_dectype(desc_a) call psb_errpush(info,name) goto 9999 endif @@ -100,7 +99,7 @@ subroutine psb_zspasb(a,desc_a, info, afmt, upd, dupl) !check on errors encountered in psdspins - spstate = a%infoa(psb_state_) + spstate = a%infoa(psb_state_) if (spstate == psb_spmat_bld_) then ! ! First case: we come from a fresh build. diff --git a/src/tools/psb_zspcnv.f90 b/src/tools/psb_zspcnv.f90 index 4ca55a95..c3d24568 100644 --- a/src/tools/psb_zspcnv.f90 +++ b/src/tools/psb_zspcnv.f90 @@ -122,8 +122,7 @@ subroutine psb_zspcnv(a,b,desc_a,info) time(1) = mpi_wtime() - - ictxt = psb_cd_get_context(desc_a) + ictxt = psb_cd_get_context(desc_a) dectype = psb_cd_get_dectype(desc_a) n_row = psb_cd_get_local_rows(desc_a) n_col = psb_cd_get_local_cols(desc_a) diff --git a/src/tools/psb_zspfree.f90 b/src/tools/psb_zspfree.f90 index dd9d0aff..eb7a0f9b 100644 --- a/src/tools/psb_zspfree.f90 +++ b/src/tools/psb_zspfree.f90 @@ -61,11 +61,11 @@ subroutine psb_zspfree(a, desc_a,info) call psb_erractionsave(err_act) if (.not.allocated(desc_a%matrix_data)) then - info=295 + info = 295 call psb_errpush(info,name) return else - ictxt=psb_cd_get_context(desc_a) + ictxt = psb_cd_get_context(desc_a) end if !...deallocate a.... diff --git a/src/tools/psb_zsphalo.f90 b/src/tools/psb_zsphalo.f90 index aa7aa483..7e72c78e 100644 --- a/src/tools/psb_zsphalo.f90 +++ b/src/tools/psb_zsphalo.f90 @@ -100,7 +100,8 @@ Subroutine psb_zsphalo(a,desc_a,blk,info,rwcnv,clcnv,outfmt) outfmt_ = 'CSR' endif - ictxt=psb_cd_get_context(desc_a) + ictxt = psb_cd_get_context(desc_a) + Call psb_info(ictxt, me, np) t1 = mpi_wtime() diff --git a/src/tools/psb_zspins.f90 b/src/tools/psb_zspins.f90 index 262bd414..712bb241 100644 --- a/src/tools/psb_zspins.f90 +++ b/src/tools/psb_zspins.f90 @@ -52,6 +52,7 @@ subroutine psb_zspins(nz,ia,ja,val,a,desc_a,info,rebuild) use psb_const_mod use psb_error_mod use psb_penv_mod + use psb_tools_mod implicit none !....parameters... @@ -63,22 +64,34 @@ subroutine psb_zspins(nz,ia,ja,val,a,desc_a,info,rebuild) logical, intent(in), optional :: rebuild !locals..... - integer :: nrow,err_act, dectype,mglob,ncol, spstate - integer :: ictxt,np, me + integer :: nrow, err_act,mglob,ncol, spstate + integer :: ictxt,np,me logical, parameter :: debug=.false. integer, parameter :: relocsz=200 logical :: rebuild_ + integer, allocatable :: ila(:),jla(:) - interface psb_cdins - subroutine psb_cdins(nz,ia,ja,desc_a,info) - use psb_descriptor_type - implicit none - type(psb_desc_type), intent(inout) :: desc_a - integer, intent(in) :: nz,ia(:),ja(:) - integer, intent(out) :: info - end subroutine psb_cdins - end interface - +!!$ interface psb_cdins +!!$ subroutine psb_cdins(nz,ia,ja,desc_a,info,ila,jla) +!!$ use psb_descriptor_type +!!$ implicit none +!!$ type(psb_desc_type), intent(inout) :: desc_a +!!$ integer, intent(in) :: nz,ia(:),ja(:) +!!$ integer, intent(out) :: info +!!$ integer, optional, intent(out) :: ila(:), jla(:) +!!$ end subroutine psb_cdins +!!$ end interface +!!$ +!!$ interface psb_glob_to_loc +!!$ subroutine psb_glob_to_loc(x,desc_a,info,iact) +!!$ use psb_descriptor_type +!!$ implicit none +!!$ type(psb_desc_type), intent(in) :: desc_a +!!$ integer, intent(inout) :: x(:) +!!$ integer, intent(out) :: info +!!$ character, intent(in), optional :: iact +!!$ end subroutine psb_glob_to_loc +!!$ end interface character(len=20) :: name, ch_err info = 0 @@ -87,8 +100,7 @@ subroutine psb_zspins(nz,ia,ja,val,a,desc_a,info,rebuild) ictxt = psb_cd_get_context(desc_a) - dectype = psb_cd_get_dectype(desc_a) - mglob = psb_cd_get_global_rows(desc_a) + mglob = psb_cd_get_global_rows(desc_a) call psb_info(ictxt, me, np) @@ -98,7 +110,7 @@ subroutine psb_zspins(nz,ia,ja,val,a,desc_a,info,rebuild) goto 9999 endif - if (nz <= 0) then + if (nz < 0) then info = 1111 call psb_errpush(info,name) goto 9999 @@ -119,6 +131,7 @@ subroutine psb_zspins(nz,ia,ja,val,a,desc_a,info,rebuild) call psb_errpush(info,name) goto 9999 end if + if (nz==0) return if (present(rebuild)) then rebuild_ = rebuild @@ -128,39 +141,102 @@ subroutine psb_zspins(nz,ia,ja,val,a,desc_a,info,rebuild) spstate = a%infoa(psb_state_) if (psb_is_bld_desc(desc_a)) then - call psb_cdins(nz,ia,ja,desc_a,info) - if (info /= 0) then - info=4010 - ch_err='psb_cdins' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 + if (psb_is_large_desc(desc_a)) then + + allocate(ila(nz),jla(nz),stat=info) + if (info /= 0) then + ch_err='allocate' + call psb_errpush(4013,name,a_err=ch_err,i_err=(/info,0,0,0,0/)) + goto 9999 + end if + call psb_cdins(nz,ia,ja,desc_a,info,ila=ila,jla=jla) + if (info /= 0) then + ch_err='psb_cdins' + call psb_errpush(4013,name,a_err=ch_err,i_err=(/info,0,0,0,0/)) + goto 9999 + end if + nrow = psb_cd_get_local_rows(desc_a) + ncol = psb_cd_get_local_cols(desc_a) + + if (spstate == psb_spmat_bld_) then + call psb_coins(nz,ila,jla,val,a,1,nrow,1,ncol,info) + if (info /= 0) then + info=4010 + ch_err='psb_coins' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + else + info = 1123 + call psb_errpush(info,name) + goto 9999 + end if + + else + + call psb_cdins(nz,ia,ja,desc_a,info) + if (info /= 0) then + ch_err='psb_cdins' + call psb_errpush(4013,name,a_err=ch_err,i_err=(/info,0,0,0,0/)) + goto 9999 + end if + nrow = psb_cd_get_local_rows(desc_a) + ncol = psb_cd_get_local_cols(desc_a) + + if (spstate == psb_spmat_bld_) then + call psb_coins(nz,ia,ja,val,a,1,nrow,1,ncol,info,gtl=desc_a%glob_to_loc) + if (info /= 0) then + info=4010 + ch_err='psb_coins' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + else + info = 1123 + call psb_errpush(info,name) + goto 9999 + end if + end if - nrow = psb_cd_get_local_rows(desc_a) - ncol = psb_cd_get_local_cols(desc_a) - if (spstate == psb_spmat_bld_) then - call psb_coins(nz,ia,ja,val,a,1,nrow,1,ncol,info,gtl=desc_a%glob_to_loc) + else if (psb_is_asb_desc(desc_a)) then + + if (psb_is_large_desc(desc_a)) then + + allocate(ila(nz),jla(nz),stat=info) + if (info /= 0) then + ch_err='allocate' + call psb_errpush(4013,name,a_err=ch_err,i_err=(/info,0,0,0,0/)) + goto 9999 + end if + + ila(1:nz) = ia(1:nz) + jla(1:nz) = ja(1:nz) + call psb_glob_to_loc(ila(1:nz),desc_a,info,iact='I') + call psb_glob_to_loc(jla(1:nz),desc_a,info,iact='I') + nrow = psb_cd_get_local_rows(desc_a) + ncol = psb_cd_get_local_cols(desc_a) + + call psb_coins(nz,ia,ja,val,a,1,nrow,1,ncol,& + & info,rebuild=rebuild_) if (info /= 0) then info=4010 ch_err='psb_coins' call psb_errpush(info,name,a_err=ch_err) goto 9999 end if + else - info = 1123 - call psb_errpush(info,name) - goto 9999 - end if - else if (psb_is_asb_desc(desc_a)) then - nrow = psb_cd_get_local_rows(desc_a) - ncol = psb_cd_get_local_cols(desc_a) - call psb_coins(nz,ia,ja,val,a,1,nrow,1,ncol,& - & info,gtl=desc_a%glob_to_loc,rebuild=rebuild_) - if (info /= 0) then - info=4010 - ch_err='psb_coins' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 + nrow = psb_cd_get_local_rows(desc_a) + ncol = psb_cd_get_local_cols(desc_a) + call psb_coins(nz,ia,ja,val,a,1,nrow,1,ncol,& + & info,gtl=desc_a%glob_to_loc,rebuild=rebuild_) + if (info /= 0) then + info=4010 + ch_err='psb_coins' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if end if else info = 1122