Major restructuring.

1. Changed descriptor allocators to distinguish between small/large
   descriptors: large descriptors don't have GLOB_TO_LOC, and use an
   AVL tree to keep track of indices while building, and a set of
   sorted lists with a hash based on low order bits for assembled
   descriptors.
2. Changed CDINS to use inner routines, hiding usage of GLOB_TO_LOC
   vs. AVL tree. Added an option to return converted indices.
3. Changed SPINS to use the new version of CDINS with the converted
   indices.
4. Changed xxINS to use new internals for selecting local indices.
5. Folded CDOVRBLD into CDOVR. Perhaps we can delede OVRBLD and rename
   OVR into OVRBLD. While doing so, changed the implementation to
   distinguish large vs small descriptors. Also changed to call cdasb
   at the end, to minimize code rewriting.
psblas3-type-indexed
Salvatore Filippone 20 years ago
parent b220064d37
commit 822eb9f59f

@ -2,8 +2,8 @@ include ../../Make.inc
FOBJS = psb_dallc.o psb_dasb.o psb_dcsrp.o psb_cdprt.o \ FOBJS = psb_dallc.o psb_dasb.o psb_dcsrp.o psb_cdprt.o \
psb_dfree.o psb_dgelp.o psb_dins.o \ psb_dfree.o psb_dgelp.o psb_dins.o \
psb_cdall.o psb_cdalv.o psb_cdasb.o psb_cdcpy.o \ psb_cdall.o psb_cdalv.o psb_cd_inloc.o psb_cdcpy.o \
psb_cddec.o psb_cdfree.o psb_cdins.o psb_dcdovr.o \ psb_cddec.o psb_cdfree.o psb_cdins.o \
psb_cdren.o psb_cdrep.o psb_cdtransfer.o psb_get_overlap.o\ psb_cdren.o psb_cdrep.o psb_cdtransfer.o psb_get_overlap.o\
psb_dspalloc.o psb_dspasb.o \ psb_dspalloc.o psb_dspasb.o \
psb_dspcnv.o psb_dspfree.o psb_dspins.o psb_dsprn.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_zspalloc.o psb_zspasb.o psb_zspcnv.o psb_zspfree.o\
psb_zspins.o psb_zsprn.o psb_zcdovr.o psb_zgelp.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 . INCDIRS = -I ../../lib -I .
LIBDIR = ../../lib LIBDIR = ../../lib
@ -26,7 +27,8 @@ lib: mpfobjs $(FOBJS)
mpfobjs: mpfobjs:
(make $(MPFOBJS) F90="$(MPF90)" FC="$(MPF90)" FCOPT="$(F90COPT)") (make $(MPFOBJS) F90="$(MPF90)" FC="$(MPF90)" FCOPT="$(F90COPT)")
psb_glob_to_loc.o: psb_glob_to_loc.f90
$(F90) $(F90COPT) $(INCDIRS) -c $< $(F90INLINE)
clean: clean:
/bin/rm -f $(MPFOBJS) $(FOBJS) /bin/rm -f $(MPFOBJS) $(FOBJS)

@ -269,6 +269,7 @@ subroutine psb_cdall(m, n, parts, ictxt, desc_a, info)
endif endif
! estimate local cols number ! estimate local cols number
loc_col=min(2*loc_row,m) loc_col=min(2*loc_row,m)
allocate(desc_a%loc_to_glob(loc_col),& allocate(desc_a%loc_to_glob(loc_col),&
&desc_a%lprm(1),stat=info) &desc_a%lprm(1),stat=info)
if (info /= 0) then if (info /= 0) then

@ -56,7 +56,7 @@ subroutine psb_cdalv(m, v, ictxt, desc_a, info, flag)
!locals !locals
Integer :: counter,i,j,np,me,loc_row,err,& 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 & l_ov_ix,l_ov_el,idx, flag_, err_act
integer :: int_err(5),exch(2) integer :: int_err(5),exch(2)
Integer, allocatable :: temp_ovrlap(:), ov_idx(:),ov_el(:) Integer, allocatable :: temp_ovrlap(:), ov_idx(:),ov_el(:)
@ -88,8 +88,8 @@ subroutine psb_cdalv(m, v, ictxt, desc_a, info, flag)
endif endif
if (info /= 0) then if (info /= 0) then
call psb_errpush(info,name,i_err=int_err) call psb_errpush(info,name,i_err=int_err)
goto 9999 goto 9999
end if end if
if (debug) write(*,*) 'psb_cdall: doing global checks' if (debug) write(*,*) 'psb_cdall: doing global checks'
@ -132,8 +132,13 @@ subroutine psb_cdalv(m, v, ictxt, desc_a, info, flag)
!count local rows number !count local rows number
! allocate work vector ! allocate work vector
allocate(desc_a%glob_to_loc(m),desc_a%matrix_data(psb_mdata_size_),& if (m > psb_cd_get_large_threshold()) then
&temp_ovrlap(m),stat=info) 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 if (info /= 0) then
info=2025 info=2025
int_err(1)=m int_err(1)=m
@ -146,34 +151,127 @@ subroutine psb_cdalv(m, v, ictxt, desc_a, info, flag)
counter = 0 counter = 0
itmpov = 0 itmpov = 0
temp_ovrlap(:) = -1 temp_ovrlap(:) = -1
do i=1,m
if ( m >psb_cd_get_large_threshold()) then
if (((v(i)-flag_) > np-1).or.((v(i)-flag_) < 0)) then desc_a%matrix_data(psb_dec_type_) = psb_desc_large_bld_
info=580
int_err(1)=3 do i=1,m
int_err(2)=v(i) - flag_
int_err(3)=i if (((v(i)-flag_) > np-1).or.((v(i)-flag_) < 0)) then
exit 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
! 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 end if
if ((v(i)-flag_) == me) then ! set LOC_TO_GLOB array to all "-1" values
! this point belongs to me desc_a%lprm(1) = 0
counter=counter+1 desc_a%loc_to_glob(:) = -1
desc_a%glob_to_loc(i) = counter k = 0
else do i=1,m
desc_a%glob_to_loc(i) = -(np+(v(i)-flag_)+1) 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 end if
enddo
loc_row=counter if (info /= 0) then
! check on parts function info=4000
if (debug) write(*,*) 'PSB_CDALL: End main loop:' ,loc_row,itmpov,info 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 end if
if (debug) write(*,*) 'PSB_CDALL: error check:' ,err
l_ov_ix=0 l_ov_ix=0
l_ov_el=0 l_ov_el=0
@ -234,27 +332,6 @@ subroutine psb_cdalv(m, v, ictxt, desc_a, info, flag)
goto 9999 goto 9999
endif 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.... ! set fields in desc_a%MATRIX_DATA....
desc_a%matrix_data(psb_n_row_) = loc_row desc_a%matrix_data(psb_n_row_) = loc_row
desc_a%matrix_data(psb_n_col_) = 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%halo_index(:) = -1
desc_a%matrix_data(psb_m_) = m desc_a%matrix_data(psb_m_) = m
desc_a%matrix_data(psb_n_) = n desc_a%matrix_data(psb_n_) = n
desc_a%matrix_data(psb_dec_type_) = psb_desc_bld_
desc_a%matrix_data(psb_ctxt_) = ictxt desc_a%matrix_data(psb_ctxt_) = ictxt
call psb_get_mpicomm(ictxt,desc_a%matrix_data(psb_mpi_c_)) call psb_get_mpicomm(ictxt,desc_a%matrix_data(psb_mpi_c_))
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
return return
9999 continue 9999 continue
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
if (err_act == act_abort) then if (err_act == act_abort) then
call psb_error(ictxt) call psb_error(ictxt)
return return
end if end if
return return

@ -45,6 +45,7 @@ subroutine psb_cdasb(desc_a,info)
use psb_penv_mod use psb_penv_mod
implicit none implicit none
include 'mpif.h'
!...Parameters.... !...Parameters....
type(psb_desc_type), intent(inout) :: desc_a type(psb_desc_type), intent(inout) :: desc_a
integer, intent(out) :: info integer, intent(out) :: info
@ -52,10 +53,11 @@ subroutine psb_cdasb(desc_a,info)
!....Locals.... !....Locals....
integer :: int_err(5), itemp(2) integer :: int_err(5), itemp(2)
integer,allocatable :: ovrlap_index(:),halo_index(:) integer,allocatable :: ovrlap_index(:),halo_index(:)
integer :: i,err,np,me,&
& lovrlap,lhalo,nhalo,novrlap,max_size,max_halo,n_col,ldesc_halo,& integer :: i,j,err,np,me,lovrlap,lhalo,nhalo,novrlap,max_size,&
& ldesc_ovrlap, dectype, err_act & max_halo,n_col,ldesc_halo, ldesc_ovrlap, dectype, err_act, &
& key, ih, nh, idx, nk,icomm,hsize
integer :: ictxt,n_row integer :: ictxt,n_row
logical, parameter :: debug=.false., debugwrt=.false. logical, parameter :: debug=.false., debugwrt=.false.
character(len=20) :: name,ch_err character(len=20) :: name,ch_err
@ -70,6 +72,7 @@ subroutine psb_cdasb(desc_a,info)
dectype = psb_cd_get_dectype(desc_a) dectype = psb_cd_get_dectype(desc_a)
n_row = psb_cd_get_local_rows(desc_a) n_row = psb_cd_get_local_rows(desc_a)
n_col = psb_cd_get_local_cols(desc_a) n_col = psb_cd_get_local_cols(desc_a)
call psb_get_mpicomm(ictxt,icomm )
! check on blacs grid ! check on blacs grid
call psb_info(ictxt, me, np) 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) call psb_errpush(info,name,i_err=int_err)
goto 9999 goto 9999
endif endif
call psb_realloc(psb_cd_get_local_cols(desc_a),desc_a%loc_to_glob,info) 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%ovrlap_index,ovrlap_index,info)
call psb_transfer(desc_a%halo_index,halo_index,info) call psb_transfer(desc_a%halo_index,halo_index,info)
@ -113,8 +121,6 @@ subroutine psb_cdasb(desc_a,info)
goto 9999 goto 9999
end if 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) deallocate(ovrlap_index, halo_index, stat=info)
if (info /= 0) then if (info /= 0) then
@ -123,6 +129,14 @@ subroutine psb_cdasb(desc_a,info)
goto 9999 goto 9999
end if 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 else
info = 600 info = 600
call psb_errpush(info,name) call psb_errpush(info,name)

@ -54,7 +54,7 @@ subroutine psb_cdcpy(desc_in, desc_out, info)
integer, intent(out) :: info integer, intent(out) :: info
!locals !locals
integer :: np,me,ictxt, isz, err_act integer :: np,me,ictxt, isz, err_act,idx,gidx,lidx
logical, parameter :: debug=.false.,debugprt=.false. logical, parameter :: debug=.false.,debugprt=.false.
character(len=20) :: name, char_err character(len=20) :: name, char_err
if (debug) write(0,*) me,'Entered CDCPY' if (debug) write(0,*) me,'Entered CDCPY'
@ -63,30 +63,53 @@ subroutine psb_cdcpy(desc_in, desc_out, info)
call psb_erractionsave(err_act) call psb_erractionsave(err_act)
name = 'psb_cdcpy' name = 'psb_cdcpy'
ictxt = psb_cd_get_context(desc_in)
ictxt=desc_in%matrix_data(psb_ctxt_)
! check on blacs grid ! check on blacs grid
call psb_info(ictxt, me, np) call psb_info(ictxt, me, np)
if (debug) write(0,*) me,'Entered CDCPY' if (debug) write(0,*) me,'Entered CDCPY'
if (np == -1) then if (np == -1) then
info = 2010 info = 2010
call psb_errpush(info,name) call psb_errpush(info,name)
goto 9999 goto 9999
endif 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) 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_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_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%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%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%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%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%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%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 if (info /= 0) then
info = 4010 info = 4010

@ -203,7 +203,6 @@ subroutine psb_cddec(nloc, ictxt, desc_a, info)
endif endif
enddo enddo
tovr = -1 tovr = -1
thalo = -1 thalo = -1
@ -215,7 +214,6 @@ subroutine psb_cddec(nloc, ictxt, desc_a, info)
goto 9999 goto 9999
end if end if
desc_a%matrix_data(psb_dec_type_) = psb_desc_asb_ desc_a%matrix_data(psb_dec_type_) = psb_desc_asb_
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)

@ -63,7 +63,7 @@ subroutine psb_cdfree(desc_a,info)
end if end if
ictxt=psb_cd_get_context(desc_a) ictxt=psb_cd_get_context(desc_a)
deallocate(desc_a%matrix_data)
call psb_info(ictxt, me, np) call psb_info(ictxt, me, np)
! ....verify blacs grid correctness.. ! ....verify blacs grid correctness..
if (np == -1) then if (np == -1) then
@ -74,7 +74,7 @@ subroutine psb_cdfree(desc_a,info)
!...deallocate desc_a.... !...deallocate desc_a....
if(.not.allocated(desc_a%loc_to_glob)) then if(.not.allocated(desc_a%loc_to_glob)) then
info=295 info=296
call psb_errpush(info,name) call psb_errpush(info,name)
goto 9999 goto 9999
end if end if
@ -87,22 +87,24 @@ subroutine psb_cdfree(desc_a,info)
goto 9999 goto 9999
end if end if
if (.not.allocated(desc_a%glob_to_loc)) then if (.not.psb_is_large_desc(desc_a)) then
info=295 if (.not.allocated(desc_a%glob_to_loc)) then
call psb_errpush(info,name) info=297
goto 9999 call psb_errpush(info,name)
end if goto 9999
end if
!deallocate glob_to_loc field !deallocate glob_to_loc field
deallocate(desc_a%glob_to_loc,stat=info) deallocate(desc_a%glob_to_loc,stat=info)
if (info /= 0) then if (info /= 0) then
info=2052 info=2052
call psb_errpush(info,name) call psb_errpush(info,name)
goto 9999 goto 9999
end if end if
endif
if (.not.allocated(desc_a%halo_index)) then if (.not.allocated(desc_a%halo_index)) then
info=295 info=298
call psb_errpush(info,name) call psb_errpush(info,name)
goto 9999 goto 9999
end if end if
@ -131,7 +133,7 @@ subroutine psb_cdfree(desc_a,info)
end if end if
if (.not.allocated(desc_a%ovrlap_index)) then if (.not.allocated(desc_a%ovrlap_index)) then
info=295 info=299
call psb_errpush(info,name) call psb_errpush(info,name)
goto 9999 goto 9999
end if end if
@ -160,6 +162,36 @@ subroutine psb_cdfree(desc_a,info)
goto 9999 goto 9999
end if 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 if (allocated(desc_a%idx_space)) then
deallocate(desc_a%idx_space,stat=info) deallocate(desc_a%idx_space,stat=info)
if (info /= 0) then if (info /= 0) then
@ -169,6 +201,8 @@ subroutine psb_cdfree(desc_a,info)
end if end if
end if end if
deallocate(desc_a%matrix_data)
call psb_nullify_desc(desc_a) call psb_nullify_desc(desc_a)
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)

@ -39,34 +39,37 @@
! ja - integer,dimension(:). The column indices of the points. ! ja - integer,dimension(:). The column indices of the points.
! desc_a - type(<psb_desc_type>). The communication descriptor to be freed. ! desc_a - type(<psb_desc_type>). The communication descriptor to be freed.
! info - integer. Eventually returns an error code. ! 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_descriptor_type
use psb_serial_mod use psb_serial_mod
use psb_const_mod use psb_const_mod
use psb_error_mod use psb_error_mod
use psb_penv_mod use psb_penv_mod
use psi_mod
implicit none implicit none
!....PARAMETERS... !....PARAMETERS...
Type(psb_desc_type), intent(inout) :: desc_a Type(psb_desc_type), intent(inout) :: desc_a
Integer, intent(in) :: nz,ia(:),ja(:) Integer, intent(in) :: nz,ia(:),ja(:)
integer, intent(out) :: info integer, intent(out) :: info
integer, optional, intent(out) :: ila(:), jla(:)
!LOCALS..... !LOCALS.....
integer :: i,ictxt,row,k,dectype,mglob, nglob,err integer :: i,ictxt,row,k,dectype,mglob, nglob,err
integer :: np, me, isize 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. logical, parameter :: debug=.false.
integer, parameter :: relocsz=200 integer, parameter :: relocsz=200
integer, allocatable :: ila_(:), jla_(:)
character(len=20) :: name,ch_err character(len=20) :: name,ch_err
info = 0 info = 0
name = 'psb_cdins' name = 'psb_cdins'
call psb_erractionsave(err_act) 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) dectype = psb_cd_get_dectype(desc_a)
mglob = psb_cd_get_global_rows(desc_a) mglob = psb_cd_get_global_rows(desc_a)
nglob = psb_cd_get_global_cols(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) call psb_errpush(info,name)
goto 9999 goto 9999
end if end if
if (present(ila)) then
if (.not.allocated(desc_a%halo_index)) then if (size(ila) < nz) then
allocate(desc_a%halo_index(relocsz)) info = 1111
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
call psb_errpush(info,name) call psb_errpush(info,name)
goto 9999 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 endif
if ((1<=desc_a%glob_to_loc(ip)).and.(desc_a%glob_to_loc(ip))<=nrow) then allocate(ila_(nz),jla_(nz),stat=info)
k = desc_a%glob_to_loc(jp) if (info /= 0) then
if (k.lt.-np) then info = 4000
k = k + np call psb_errpush(info,name)
k = - k - 1 goto 9999
ncol = ncol + 1 end if
desc_a%glob_to_loc(jp) = ncol call psi_idx_cnv(nz,ia,ila_,desc_a,info,owned=.true.)
isize = size(desc_a%loc_to_glob) call psi_idx_ins_cnv(nz,ja,jla_,desc_a,info,mask=(ila_(1:nz)>0))
if (ncol > isize) then deallocate(ila_,jla_)
nh = ncol + max(nz,relocsz) end if
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
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
return return

@ -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), ': ',& & write(iout,*) 'Loc_to_glob ',desc_p%loc_to_glob(1:n_row), ': ',&
& desc_p%loc_to_glob(n_row+1:n_col) & 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' write(iout,*) 'Halo_index'
counter = 1 counter = 1
Do Do

@ -75,10 +75,10 @@ subroutine psb_cdren(trans,iperm,desc_a,info)
time(1) = mpi_wtime() 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) dectype = psb_cd_get_dectype(desc_a)
n_row = psb_cd_get_local_rows(desc_a) n_row = psb_cd_get_local_rows(desc_a)
n_col = psb_cd_get_local_cols(desc_a) n_col = psb_cd_get_local_cols(desc_a)
! check on blacs grid ! check on blacs grid
call psb_info(ictxt, me, np) call psb_info(ictxt, me, np)

@ -186,7 +186,6 @@ subroutine psb_cdrep(m, ictxt, desc_a, info)
endif endif
desc_a%matrix_data(psb_m_) = m desc_a%matrix_data(psb_m_) = m
desc_a%matrix_data(psb_n_) = n desc_a%matrix_data(psb_n_) = n
desc_a%matrix_data(psb_n_row_) = m 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_dec_type_) = psb_desc_bld_
desc_a%matrix_data(psb_ctxt_) = ictxt desc_a%matrix_data(psb_ctxt_) = ictxt
call psb_get_mpicomm(ictxt,desc_a%matrix_data(psb_mpi_c_)) call psb_get_mpicomm(ictxt,desc_a%matrix_data(psb_mpi_c_))
do i=1,m do i=1,m
desc_a%glob_to_loc(i) = i desc_a%glob_to_loc(i) = i
desc_a%loc_to_glob(i) = i desc_a%loc_to_glob(i) = i
enddo enddo
tovr = -1 tovr = -1
thalo = -1 thalo = -1
desc_a%lprm(:) = 0 desc_a%lprm(:) = 0
call psi_cnv_dsc(thalo,tovr,desc_a,info) call psi_cnv_dsc(thalo,tovr,desc_a,info)
if (info /= 0) then 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 goto 9999
end if end if

@ -50,7 +50,7 @@ subroutine psb_cdtransfer(desc_in, desc_out, info)
!....parameters... !....parameters...
type(psb_desc_type), intent(inout) :: desc_in 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 integer, intent(out) :: info
!locals !locals
@ -64,29 +64,47 @@ subroutine psb_cdtransfer(desc_in, desc_out, info)
call psb_erractionsave(err_act) call psb_erractionsave(err_act)
name = 'psb_cdtransfer' name = 'psb_cdtransfer'
ictxt=desc_in%matrix_data(psb_ctxt_) ictxt=psb_cd_get_context(desc_in)
call psb_info(ictxt, me, np) call psb_info(ictxt, me, np)
if (debug) write(0,*) me,'Entered CDTRANSFER' if (debug) write(0,*) me,'Entered CDTRANSFER'
if (np == -1) then if (np == -1) then
info = 2010 info = 2010
call psb_errpush(info,name) call psb_errpush(info,name)
goto 9999 goto 9999
endif endif
!!$ call psb_nullify_desc(desc_out)
!!$
call psb_transfer( desc_in%matrix_data , desc_out%matrix_data , info) call psb_transfer( desc_in%matrix_data , desc_out%matrix_data , info)
call psb_transfer( desc_in%halo_index , desc_out%halo_index , info) if (info == 0) call psb_transfer( desc_in%halo_index , desc_out%halo_index , info)
call psb_transfer( desc_in%bnd_elem , desc_out%bnd_elem , info) if (info == 0) call psb_transfer( desc_in%bnd_elem , desc_out%bnd_elem , info)
call psb_transfer( desc_in%ovrlap_elem , desc_out%ovrlap_elem , info) if (info == 0) call psb_transfer( desc_in%ovrlap_elem , desc_out%ovrlap_elem , info)
call psb_transfer( desc_in%ovrlap_index, desc_out%ovrlap_index , info) if (info == 0) 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) if (info == 0) 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) if (info == 0) call psb_transfer( desc_in%glob_to_loc , desc_out%glob_to_loc , info)
call psb_transfer( desc_in%lprm , desc_out%lprm , info) if (info == 0) 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%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) call psb_erractionrestore(err_act)
return return
@ -95,9 +113,9 @@ subroutine psb_cdtransfer(desc_in, desc_out, info)
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
if (err_act == act_ret) then if (err_act == act_ret) then
return return
else else
call psb_error(ictxt) call psb_error(ictxt)
end if end if
return return

@ -51,13 +51,13 @@ subroutine psb_dalloc(x, desc_a, info, n)
!....parameters... !....parameters...
real(kind(1.d0)), allocatable, intent(out) :: x(:,:) real(kind(1.d0)), allocatable, intent(out) :: x(:,:)
type(psb_desc_type), intent(in) :: desc_a type(psb_desc_type), intent(in) :: desc_a
integer :: info integer,intent(out) :: info
integer, optional, intent(in) :: n integer, optional, intent(in) :: n
!locals !locals
integer :: np,me,err,n_col,n_row,i,j,err_act 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) integer :: int_err(5), exch(3)
character(len=20) :: name, ch_err character(len=20) :: name, ch_err
@ -77,7 +77,6 @@ subroutine psb_dalloc(x, desc_a, info, n)
goto 9999 goto 9999
endif endif
dectype=psb_cd_get_dectype(desc_a)
!... check m and n parameters.... !... check m and n parameters....
if (.not.psb_is_ok_desc(desc_a)) then if (.not.psb_is_ok_desc(desc_a)) then
info = 3110 info = 3110
@ -199,11 +198,11 @@ subroutine psb_dallocv(x, desc_a,info,n)
!....parameters... !....parameters...
real(kind(1.d0)), allocatable, intent(out) :: x(:) real(kind(1.d0)), allocatable, intent(out) :: x(:)
type(psb_desc_type), intent(in) :: desc_a type(psb_desc_type), intent(in) :: desc_a
integer :: info integer,intent(out) :: info
integer, optional, intent(in) :: n integer, optional, intent(in) :: n
!locals !locals
integer :: np,me,n_col,n_row,dectype,i,err_act integer :: np,me,n_col,n_row,i,err_act
integer :: ictxt integer :: ictxt
logical, parameter :: debug=.false. logical, parameter :: debug=.false.
character(len=20) :: name, ch_err character(len=20) :: name, ch_err
@ -223,9 +222,6 @@ subroutine psb_dallocv(x, desc_a,info,n)
goto 9999 goto 9999
endif 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.... !... check m and n parameters....
if (.not.psb_is_ok_desc(desc_a)) then if (.not.psb_is_ok_desc(desc_a)) then
info = 3110 info = 3110
@ -261,6 +257,8 @@ subroutine psb_dallocv(x, desc_a,info,n)
do i=1,n_row do i=1,n_row
x(i) = 0.0d0 x(i) = 0.0d0
end do end do
else
write(0,*) 'Did not allocate anything because of dectype',psb_cd_get_dectype(desc_a)
endif endif
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)

@ -53,7 +53,7 @@ subroutine psb_dasb(x, desc_a, info)
! local variables ! local variables
integer :: ictxt,np,me,nrow,ncol, err_act 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 real(kind(1.d0)),parameter :: one=1
logical, parameter :: debug=.false. logical, parameter :: debug=.false.
character(len=20) :: name, ch_err character(len=20) :: name, ch_err
@ -69,14 +69,13 @@ subroutine psb_dasb(x, desc_a, info)
goto 9999 goto 9999
endif endif
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) call psb_info(ictxt, me, np)
if (debug) write(*,*) 'asb start: ',np,me,& if (debug) write(*,*) 'asb start: ',np,me,&
&psb_cd_get_dectype(desc_a) & psb_cd_get_dectype(desc_a)
! ....verify blacs grid correctness.. ! ....verify blacs grid correctness..
if (np == -1) then if (np == -1) then
info = 2010 info = 2010
@ -84,16 +83,16 @@ subroutine psb_dasb(x, desc_a, info)
goto 9999 goto 9999
else if (.not.psb_is_asb_desc(desc_a)) then else if (.not.psb_is_asb_desc(desc_a)) then
if (debug) write(*,*) 'asb error ',& if (debug) write(*,*) 'asb error ',&
&dectype & psb_cd_get_dectype(desc_a)
info = 3110 info = 3110
call psb_errpush(info,name) call psb_errpush(info,name)
goto 9999 goto 9999
endif endif
! check size ! check size
ictxt=psb_cd_get_context(desc_a) ictxt = psb_cd_get_context(desc_a)
nrow=psb_cd_get_local_rows(desc_a) nrow = psb_cd_get_local_rows(desc_a)
ncol=psb_cd_get_local_cols(desc_a) ncol = psb_cd_get_local_cols(desc_a)
i1sz = size(x,dim=1) i1sz = size(x,dim=1)
i2sz = size(x,dim=2) i2sz = size(x,dim=2)
if (debug) write(*,*) 'asb: ',i1sz,i2sz,nrow,ncol if (debug) write(*,*) 'asb: ',i1sz,i2sz,nrow,ncol
@ -184,7 +183,7 @@ subroutine psb_dasbv(x, desc_a, info)
! local variables ! local variables
integer :: ictxt,np,me 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 real(kind(1.d0)),parameter :: one=1
logical, parameter :: debug=.false. logical, parameter :: debug=.false.
character(len=20) :: name,ch_err character(len=20) :: name,ch_err
@ -193,8 +192,7 @@ subroutine psb_dasbv(x, desc_a, info)
int_err(1) = 0 int_err(1) = 0
name = 'psb_dasbv' name = 'psb_dasbv'
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) call psb_info(ictxt, me, np)
@ -209,8 +207,8 @@ subroutine psb_dasbv(x, desc_a, info)
goto 9999 goto 9999
endif endif
nrow=psb_cd_get_local_rows(desc_a) nrow = psb_cd_get_local_rows(desc_a)
ncol=psb_cd_get_local_cols(desc_a) ncol = psb_cd_get_local_cols(desc_a)
if (debug) write(*,*) name,' sizes: ',nrow,ncol if (debug) write(*,*) name,' sizes: ',nrow,ncol
i1sz = size(x) i1sz = size(x)
if (debug) write(*,*) 'dasb: sizes ',i1sz,ncol if (debug) write(*,*) 'dasb: sizes ',i1sz,ncol

@ -40,7 +40,8 @@
! a - type(<psb_dspmat_type>). The input sparse matrix. ! a - type(<psb_dspmat_type>). The input sparse matrix.
! desc_a - type(<psb_desc_type>). The input communication descriptor. ! desc_a - type(<psb_desc_type>). The input communication descriptor.
! norv - integer. The number of overlap levels. ! norv - integer. The number of overlap levels.
! desc_ov - type(<psb_desc_type>). The auxiliary output communication descriptor. ! desc_ov - type(<psb_desc_type>). The auxiliary output communication
! descriptor.
! info - integer. Eventually returns an error code. ! info - integer. Eventually returns an error code.
! !
Subroutine psb_dcdovr(a,desc_a,novr,desc_ov,info) 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_prec_mod
use psb_error_mod use psb_error_mod
use psb_penv_mod use psb_penv_mod
use psb_tools_mod
use psb_realloc_mod
use psi_mod
use mpi
Implicit None Implicit None
! .. Array Arguments .. ! .. Array Arguments ..
@ -60,40 +65,23 @@ Subroutine psb_dcdovr(a,desc_a,novr,desc_ov,info)
Type(psb_desc_type), Intent(inout) :: desc_ov Type(psb_desc_type), Intent(inout) :: desc_ov
integer, intent(out) :: info integer, intent(out) :: info
real(kind(1.d0)) :: t1,t2,t3,mpi_wtime
external mpi_wtime
integer icomm, err_act 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 .. ! .. Local Scalars ..
Integer :: i, j, k, np, me,m,nnzero,& 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 & 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. Logical,Parameter :: debug=.false.
character(len=20) :: name, ch_err 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 If(debug) Write(0,*)'in psb_cdovr',novr
m=psb_cd_get_local_rows(desc_a) m = psb_cd_get_local_rows(desc_a)
nnzero=Size(a%aspk) nnzero = Size(a%aspk)
n_col=psb_cd_get_local_cols(desc_a) n_row = psb_cd_get_local_rows(desc_a)
nhalo = n_col-m 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(debug) Write(0,*)'IN CDOVR1',novr ,m,nnzero,n_col
if (novr<0) then if (novr<0) then
info=10 info=10
@ -120,34 +110,32 @@ Subroutine psb_dcdovr(a,desc_a,novr,desc_ov,info)
goto 9999 goto 9999
endif 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 if (novr==0) then
! !
! Just copy the input. ! 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 return
endif endif
call psb_get_mpicomm(ictxt,icomm ) call psb_get_mpicomm(ictxt,icomm )
If(debug) then If(debug)then
Write(0,*)'BEGIN cdovr',me,nhalo Write(0,*)'BEGIN cdovr',me,nhalo
call psb_barrier(ictxt) call psb_barrier(ictxt)
endif endif
t1 = mpi_wtime()
! !
! Ok, since we are only estimating, do it as follows: ! 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. ! nonzeros per row is the same as the global.
! !
nztot = psb_sp_get_nnzeros(a) nztot = psb_sp_get_nnzeros(a)
@ -161,54 +149,552 @@ Subroutine psb_dcdovr(a,desc_a,novr,desc_ov,info)
goto 9999 goto 9999
endif endif
If(debug)Write(0,*)'ovr_est done',me,novr,lovr If(debug)Write(0,*)'ovr_est done',me,novr,lovr
index_dim = size(desc_a%halo_index) index_dim = size(desc_a%halo_index)
elem_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(:) l_tmp_ovr_idx = novr*(3*Max(2*index_dim,1)+1)
desc_ov%matrix_data(psb_dec_type_) = psb_desc_bld_ l_tmp_halo = novr*(3*Size(desc_a%halo_index))
Allocate(desc_ov%loc_to_glob(Size(desc_a%loc_to_glob)),& if (psb_is_large_desc(desc_a)) then
& desc_ov%glob_to_loc(Size(desc_a%glob_to_loc)),stat=info) desc_ov%matrix_data(psb_dec_type_) = psb_desc_large_bld_
if (info /= 0) then else
info=4000 desc_ov%matrix_data(psb_dec_type_) = psb_desc_bld_
call psb_errpush(info,name)
goto 9999
end if 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 If(debug) then
Write(0,*)'Start cdovrbld',me,lworks,lworkr Write(0,*)'Start cdovrbld',me,lworks,lworkr
call psb_barrier(ictxt) call psb_barrier(ictxt)
endif 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 if (.false.) then
info=4010 !
ch_err='psb_cdovrbld' ! The real work goes on in here....
call psb_errpush(info,name,a_err=ch_err) !
goto 9999 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<J
!
proc_id = temp(i)
call psb_check_size((counter_t+3),t_halo_in,info,pad=-1)
if (info /= 0) then
info=4010
call psb_errpush(info,name,a_err='psb_check_size')
goto 9999
end if
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_ov%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<J
!
n_col=n_col+1
proc_id=-desc_ov%glob_to_loc(idx)-np-1
call psb_check_size(n_col,desc_ov%loc_to_glob,info,pad=-1)
if (info /= 0) then
info=4010
call psb_errpush(info,name,a_err='psb_check_size')
goto 9999
end if
desc_ov%glob_to_loc(idx)=n_col
desc_ov%loc_to_glob(n_col)=idx
call psb_check_size((counter_t+3),t_halo_in,info,pad=-1)
if (info /= 0) then
info=4010
call psb_errpush(info,name,a_err='psb_check_size')
goto 9999
end if
t_halo_in(counter_t)=proc_id
t_halo_in(counter_t+1)=1
t_halo_in(counter_t+2)=n_col
counter_t=counter_t+3
if (debug) write(0,*) me,' CDOVRBLD: Added into t_halo_in from recv',&
&proc_id,n_col,idx
else if (desc_ov%glob_to_loc(idx) < 0) Then
if (debug) write(0,*) me,'Wrong input to cdovrbld??',&
&idx,desc_ov%glob_to_loc(idx)
End If
End Do
desc_ov%matrix_data(psb_n_col_)=n_col
end if
end if
!!$ desc_ov%matrix_data(psb_n_row_)=desc_ov%matrix_data(psb_n_col_)
!
! Ok, now we have a temporary halo with all the info for the
! next round. If we need to keep going, convert the halo format
! from temporary to final, so that we can work out the next iteration.
! This uses one of the convert_comm internals, i.e. we are doing
! the equivalent of a partial call to convert_comm
!
If (i_ovr < (novr)) Then
t_halo_in(counter_t)=-1
if (debug) write(0,*) me,'Checktmp_o_i 1',tmp_ovr_idx(1:10)
if (debug) write(0,*) me,'Calling Crea_Halo'
call psi_crea_index(desc_ov,t_halo_in,t_halo_out,.false.,&
& nxch,nsnd,nrcv,info)
if (debug) then
write(0,*) me,'Done Crea_Index'
call psb_barrier(ictxt)
end if
if (debug) write(0,*) me,'Checktmp_o_i 2',tmp_ovr_idx(1:10)
if (debug) write(0,*) me,'Done Crea_Halo'
call psb_transfer(t_halo_out,halo,info)
!
! At this point we have built the halo necessary for I_OVR+1.
!
End If
if (debug) write(0,*) me,'Checktmp_o_i ',tmp_ovr_idx(1:10)
End Do
desc_ov%matrix_data(psb_m_)=psb_cd_get_global_rows(desc_a)
desc_ov%matrix_data(psb_n_)=psb_cd_get_global_cols(desc_a)
tmp_halo(counter_h:)=-1
tmp_ovr_idx(counter_o:)=-1
!
! At this point we have gathered all the indices in the halo at
! N levels of overlap. Just call cnv_dsc. This is
! the same routine as gets called inside CDASB.
!
if (debug) then
write(0,*) 'psb_cdovrbld: converting indexes'
call psb_barrier(ictxt)
end if
!.... convert comunication stuctures....
! Note that we have to keep local_rows until the very end,
! because otherwise the halo build mechanism of cdasb
! will not work.
call psb_transfer(tmp_ovr_idx,desc_ov%ovrlap_index,info)
call psb_transfer(tmp_halo,desc_ov%halo_index,info)
call psb_cdasb(desc_ov,info)
desc_ov%matrix_data(psb_n_row_)=desc_ov%matrix_data(psb_n_col_)
if (debug) then
write(0,*) me,'Done CDASB'
call psb_barrier(ictxt)
end if
if (info == 0) call psb_sp_free(blk,info)
if (info /= 0) then
ch_err='sp_free'
call psb_errpush(4013,name,a_err=ch_err,i_err=(/info,0,0,0,0/))
goto 9999
end if
end if end if
desc_ov%matrix_data(psb_dec_type_) = psb_desc_asb_
If(debug) then
Write(0,*)'Done cdovrbld',me,lworks,lworkr
call psb_barrier(ictxt)
endif
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
return return

@ -54,34 +54,30 @@ Subroutine psb_dcdovrbld(n_ovr,desc_p,desc_a,a,&
use psb_serial_mod use psb_serial_mod
Use psi_mod Use psi_mod
use psb_realloc_mod use psb_realloc_mod
use psb_tools_mod, only : psb_cdprt
use psb_error_mod use psb_error_mod
use psb_const_mod use psb_const_mod
use psb_penv_mod use psb_penv_mod
use psb_tools_mod
use mpi
Implicit None Implicit None
include 'mpif.h'
type(psb_dspmat_type),intent(in) :: a type(psb_dspmat_type),intent(in) :: a
type(psb_desc_type),intent(in) :: desc_a type(psb_desc_type),intent(in) :: desc_a
type(psb_desc_type),intent(inout) :: desc_p 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??? ! 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(inout) :: lworks, lworkr
integer, intent(out) :: info integer, intent(out) :: info
type(psb_dspmat_type) :: blk type(psb_dspmat_type) :: blk
Integer, allocatable :: tmp_halo(:),tmp_ovr_idx(:)
Integer, allocatable :: tmp_halo(:),tmp_ovr_idx(:)
Integer :: counter,counter_h, counter_o, counter_e,j,idx,gidx,proc,n_elem_recv,& 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,& & 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, & & 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 & idxr, idxs, lx, iszr, iszs, err_act, icomm, nxch, nsnd, nrcv,lidx
Integer,allocatable :: halo(:),works(:),workr(:),t_halo_in(:),&
Integer,allocatable :: halo(:),length_dl(:),works(:),workr(:),t_halo_in(:),& & t_halo_out(:),temp(:),maskr(:)
& t_halo_out(:),work(:),dep_list(:),temp(:)
Integer,allocatable :: brvindx(:),rvsz(:), bsdindx(:),sdsz(:) Integer,allocatable :: brvindx(:),rvsz(:), bsdindx(:),sdsz(:)
Logical,Parameter :: debug=.false. 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) n_col = psb_cd_get_local_cols(desc_a)
if (debug) write(0,*) me,' On entry to CDOVRBLD n_col:',n_col 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),& Allocate(works(lworks),workr(lworkr),t_halo_in(l_tmp_halo),&
& t_halo_out(l_tmp_halo), work(lwork),& & t_halo_out(l_tmp_halo), temp(lworkr),stat=info)
& length_dl(np+1),dep_list(dl_lda*np),temp(lworkr),stat=info)
if (info /= 0) then if (info /= 0) then
call psb_errpush(4010,name,a_err='Allocate') call psb_errpush(4010,name,a_err='Allocate')
goto 9999 goto 9999
@ -141,7 +134,7 @@ Subroutine psb_dcdovrbld(n_ovr,desc_p,desc_a,a,&
call psb_errpush(4010,name,a_err='Allocate') call psb_errpush(4010,name,a_err='Allocate')
goto 9999 goto 9999
end if end if
halo(:) = desc_a%halo_index(:) halo(:) = desc_a%halo_index(:)
desc_p%ovrlap_elem(:) = -1 desc_p%ovrlap_elem(:) = -1
tmp_ovr_idx(:) = -1 tmp_ovr_idx(:) = -1
tmp_halo(:) = -1 tmp_halo(:) = -1
@ -150,7 +143,40 @@ Subroutine psb_dcdovrbld(n_ovr,desc_p,desc_a,a,&
counter_h = 1 counter_h = 1
counter_o = 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. ! 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 Do i_ovr=1,n_ovr
if (debug) write(0,*) me,'Running on overlap level ',i_ovr,' of ',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 ! At this point, halo contains a valid halo corresponding to the
@ -193,6 +218,7 @@ Subroutine psb_dcdovrbld(n_ovr,desc_p,desc_a,a,&
counter = 1 counter = 1
counter_t = 1 counter_t = 1
t1 = mpi_wtime() t1 = mpi_wtime()
Do While (halo(counter) /= -1) Do While (halo(counter) /= -1)
tot_elem=0 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_recv=halo(counter+psb_n_elem_recv_)
n_elem_send=halo(counter+n_elem_recv+psb_n_elem_send_) n_elem_send=halo(counter+n_elem_recv+psb_n_elem_send_)
If ((counter+n_elem_recv+n_elem_send) > Size(halo)) then If ((counter+n_elem_recv+n_elem_send) > Size(halo)) then
info=-1 info = -1
call psb_errpush(info,name) call psb_errpush(info,name)
goto 9999 goto 9999
end If 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 ! 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 ! 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 ! be enlarged with the new column indices received, and will reassemble
! everything for the next iteration. ! everything for the next iteration.
! !
@ -225,7 +251,7 @@ Subroutine psb_dcdovrbld(n_ovr,desc_p,desc_a,a,&
call psb_errpush(info,name) call psb_errpush(info,name)
goto 9999 goto 9999
end If end If
idx = halo(counter+psb_elem_recv_+j)
idx = halo(counter+psb_elem_recv_+j) idx = halo(counter+psb_elem_recv_+j)
If(idx > Size(desc_p%loc_to_glob)) then If(idx > Size(desc_p%loc_to_glob)) then
info=-3 info=-3
@ -235,42 +261,33 @@ Subroutine psb_dcdovrbld(n_ovr,desc_p,desc_a,a,&
gidx = desc_p%loc_to_glob(idx) gidx = desc_p%loc_to_glob(idx)
If((counter_o+2) > Size(tmp_ovr_idx)) Then call psb_check_size((counter_o+3),tmp_ovr_idx,info,pad=-1)
isz = max((3*Size(tmp_ovr_idx))/2,(counter_o+3)) if (info /= 0) then
if (debug) write(0,*) me,'Realloc tmp_ovr',isz info=4010
call psb_realloc(isz,tmp_ovr_idx,info,pad=-1) call psb_errpush(info,name,a_err='psb_check_size')
if (info /= 0) then goto 9999
info=4010 end if
ch_err='psb_realloc'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
End If
tmp_ovr_idx(counter_o)=proc tmp_ovr_idx(counter_o)=proc
tmp_ovr_idx(counter_o+1)=1 tmp_ovr_idx(counter_o+1)=1
tmp_ovr_idx(counter_o+2)=gidx tmp_ovr_idx(counter_o+2)=gidx
tmp_ovr_idx(counter_o+3)=-1 tmp_ovr_idx(counter_o+3)=-1
counter_o=counter_o+3 counter_o=counter_o+3
if (.not.psb_is_large_desc(desc_p)) then
If((counter_h+2) > Size(tmp_halo)) Then call psb_check_size((counter_h+3),tmp_halo,info,pad=-1)
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 (info /= 0) then if (info /= 0) then
info=4010 info=4010
ch_err='psb_realloc' call psb_errpush(info,name,a_err='psb_check_size')
call psb_errpush(info,name,a_err=ch_err)
goto 9999 goto 9999
end if end if
End If
tmp_halo(counter_h)=proc tmp_halo(counter_h)=proc
tmp_halo(counter_h+1)=1 tmp_halo(counter_h+1)=1
tmp_halo(counter_h+2)=idx tmp_halo(counter_h+2)=idx
tmp_halo(counter_h+3)=-1 tmp_halo(counter_h+3)=-1
counter_h=counter_h+3 counter_h=counter_h+3
end if
Enddo Enddo
if (debug) write(0,*) me,'Checktmp_o_i Loop Mid1',tmp_ovr_idx(1:10) 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) idx = halo(counter+psb_elem_send_+j)
gidx = desc_p%loc_to_glob(idx) gidx = desc_p%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)
If((counter_o+2) > Size(tmp_ovr_idx)) Then call psb_check_size((counter_o+3),tmp_ovr_idx,info,pad=-1)
isz = max((3*Size(tmp_ovr_idx))/2,(counter_o+3)) if (info /= 0) then
if (debug) write(0,*) me,'Realloc tmp_ovr',isz info=4010
call psb_realloc(isz,tmp_ovr_idx,info) call psb_errpush(info,name,a_err='psb_check_size')
if (info /= 0) then goto 9999
info=4010 end if
ch_err='psrealloc'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
End If
tmp_ovr_idx(counter_o)=proc tmp_ovr_idx(counter_o)=proc
tmp_ovr_idx(counter_o+1)=1 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 If (i_ovr < (n_ovr)) Then
n_elem = psb_sp_get_nnz_row(idx,a) n_elem = psb_sp_get_nnz_row(idx,a)
If((idxs+tot_elem+n_elem) > lworks) Then call psb_check_size((idxs+tot_elem+n_elem),works,info)
isz = max((3*lworks)/2,(idxs+tot_elem+n_elem)) if (info /= 0) then
if (debug) write(0,*) me,'Realloc works',isz info=4010
call psb_realloc(isz,works,info) call psb_errpush(info,name,a_err='psb_check_size')
if (info /= 0) then goto 9999
info=4010 end if
ch_err='psb_realloc'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
lworks = isz
End If
If((n_elem) > size(blk%ia2)) Then If((n_elem) > size(blk%ia2)) Then
isz = max((3*size(blk%ia2))/2,(n_elem)) 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 goto 9999
end if end if
Do i=1,iszr if (debug) write(0,*) 'ISZR :',iszr
idx=workr(i)
if (idx <1) then if (psb_is_large_desc(desc_a)) then
write(0,*) me,'Error in CDOVRBLD ',idx,i,iszr call psb_check_size(iszr,maskr,info)
!!$ write(0,*) me, ' WORKR :',workr(1:iszr) if (info /= 0) then
else If (desc_p%glob_to_loc(idx) < -np) Then info=4010
! call psb_errpush(info,name,a_err='psb_check_size')
! This is a new index. Assigning a local index as goto 9999
! we receive the guarantees that all indices for HALO(I) end if
! will be less than those for HALO(J) whenever I<J call psi_idx_cnv(iszr,workr,maskr,desc_p,info)
! iszs = count(maskr(1:iszr)<=0)
n_col=n_col+1 if (iszs > size(works)) call psb_realloc(iszs,works,info)
proc_id=-desc_p%glob_to_loc(idx)-np-1 j = 0
If(n_col > Size(desc_p%loc_to_glob)) Then do i=1,iszr
isz = 3*n_col/2 if (maskr(i) < 0) then
if (debug) write(0,*) me,'Realloc loc_to_glob' j=j+1
call psb_realloc(isz,desc_p%loc_to_glob,info) 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<J
!
proc_id = temp(i)
call psb_check_size((counter_t+3),t_halo_in,info,pad=-1)
if (info /= 0) then if (info /= 0) then
info=4010 info=4010
ch_err='psrealloc' call psb_errpush(info,name,a_err='psb_check_size')
call psb_errpush(info,name,a_err=ch_err)
goto 9999 goto 9999
end if end if
End If
desc_p%glob_to_loc(idx)=n_col t_halo_in(counter_t)=proc_id
desc_p%loc_to_glob(n_col)=idx t_halo_in(counter_t+1)=1
If ((counter_t+3) > Size(t_halo_in)) then t_halo_in(counter_t+2)=lidx
isz = max((3*Size(t_halo_in))/2,(counter_t+3+1000)) counter_t=counter_t+3
if (debug) write(0,*) me,'Realloc ovr_el',isz if (.false.) write(0,*) me,' CDOVRBLD: Added t_halo_in ',&
call psb_realloc(isz,t_halo_in,info) &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<J
!
n_col=n_col+1
proc_id=-desc_p%glob_to_loc(idx)-np-1
call psb_check_size(n_col,desc_p%loc_to_glob,info,pad=-1)
if (info /= 0) then if (info /= 0) then
info=4010 info=4010
ch_err='psrealloc' call psb_errpush(info,name,a_err='psb_check_size')
call psb_errpush(info,name,a_err=ch_err)
goto 9999 goto 9999
end if end if
end If
desc_p%glob_to_loc(idx)=n_col
t_halo_in(counter_t)=proc_id desc_p%loc_to_glob(n_col)=idx
t_halo_in(counter_t+1)=1
t_halo_in(counter_t+2)=n_col call psb_check_size((counter_t+3),t_halo_in,info,pad=-1)
counter_t=counter_t+3 if (info /= 0) then
if (debug) write(0,*) me,' CDOVRBLD: Added into t_halo_in from recv',& info=4010
&proc_id,n_col,idx call psb_errpush(info,name,a_err='psb_check_size')
else if (desc_p%glob_to_loc(idx) < 0) Then goto 9999
if (debug) write(0,*) me,'Wrong input to cdovrbld??',& end if
&idx,desc_p%glob_to_loc(idx)
End If t_halo_in(counter_t)=proc_id
End Do t_halo_in(counter_t+1)=1
t_halo_in(counter_t+2)=n_col
counter_t=counter_t+3
if (debug) write(0,*) me,' CDOVRBLD: Added into t_halo_in from recv',&
&proc_id,n_col,idx
else if (desc_p%glob_to_loc(idx) < 0) Then
if (debug) write(0,*) me,'Wrong input to cdovrbld??',&
&idx,desc_p%glob_to_loc(idx)
End If
End Do
desc_p%matrix_data(psb_n_col_)=n_col
end if
end if end if
t2 = mpi_wtime() t2 = mpi_wtime()
n_row=m+tot_recv !!$ desc_p%matrix_data(psb_n_row_)=desc_p%matrix_data(psb_n_col_)
desc_p%matrix_data(psb_n_row_)=n_row
desc_p%matrix_data(psb_n_col_)=n_col
! !
! Ok, now we have a temporary halo with all the info for the ! Ok, now we have a temporary halo with all the info for the
! next round. If we need to keep going, convert the halo format ! next round. If we need to keep going, convert the halo format
@ -493,18 +549,6 @@ Subroutine psb_dcdovrbld(n_ovr,desc_p,desc_a,a,&
If (i_ovr < (n_ovr)) Then If (i_ovr < (n_ovr)) Then
If(lwork < (counter_t/3+np*3)) Then
isz = max((3*lwork)/2,(counter_t/3+np*3))
if (debug) write(0,*) me,'Realloc work',isz
deallocate(work)
allocate(work(isz),stat=info)
if (info /= 0) then
call psb_errpush(4010,name,a_err='Allocate')
goto 9999
end if
lwork=size(work)
Endif
t_halo_in(counter_t)=-1 t_halo_in(counter_t)=-1
if (debug) write(0,*) me,'Checktmp_o_i 1',tmp_ovr_idx(1:10) if (debug) write(0,*) me,'Checktmp_o_i 1',tmp_ovr_idx(1:10)
@ -518,7 +562,7 @@ Subroutine psb_dcdovrbld(n_ovr,desc_p,desc_a,a,&
end if end if
if (debug) write(0,*) me,'Checktmp_o_i 2',tmp_ovr_idx(1:10) if (debug) write(0,*) me,'Checktmp_o_i 2',tmp_ovr_idx(1:10)
if (debug) write(0,*) me,'Done Crea_Halo' if (debug) write(0,*) me,'Done Crea_Halo'
call psb_transfer( t_halo_out,halo,info) call psb_transfer(t_halo_out,halo,info)
! !
! At this point we have built the halo necessary for I_OVR+1. ! At this point we have built the halo necessary for I_OVR+1.
! !
@ -534,14 +578,14 @@ Subroutine psb_dcdovrbld(n_ovr,desc_p,desc_a,a,&
desc_p%matrix_data(psb_m_)=psb_cd_get_global_rows(desc_a) desc_p%matrix_data(psb_m_)=psb_cd_get_global_rows(desc_a)
desc_p%matrix_data(psb_n_)=psb_cd_get_global_cols(desc_a) desc_p%matrix_data(psb_n_)=psb_cd_get_global_cols(desc_a)
tmp_halo(counter_h)=-1 tmp_halo(counter_h:)=-1
tmp_ovr_idx(counter_o)=-1 tmp_ovr_idx(counter_o:)=-1
! !
! At this point we have gathered all the indices in the halo at ! At this point we have gathered all the indices in the halo at
! N levels of overlap. Just call convert_comm. This is ! N levels of overlap. Just call cnv_dsc. This is
! the same routine as gets called inside SPASB. ! the same routine as gets called inside CDASB.
! !
if (debug) then if (debug) then
@ -549,45 +593,19 @@ Subroutine psb_dcdovrbld(n_ovr,desc_p,desc_a,a,&
call psb_barrier(ictxt) call psb_barrier(ictxt)
end if end if
!.... convert comunication stuctures.... !.... convert comunication stuctures....
! Note that we have to keep local_rows until the very end,
call psi_cnv_dsc(tmp_halo,tmp_ovr_idx,desc_p,info) ! because otherwise the halo build mechanism of cdasb
! will not work.
! Ok, register into MATRIX_DATA & free temporary work areas call psb_transfer(tmp_ovr_idx,desc_p%ovrlap_index,info)
desc_p%matrix_data(psb_dec_type_) = psb_desc_asb_ call psb_transfer(tmp_halo,desc_p%halo_index,info)
call psb_cdasb(desc_p,info)
allocate(desc_p%lprm(1),stat=info) desc_p%matrix_data(psb_n_row_)=desc_p%matrix_data(psb_n_col_)
if (info /= 0) then
call psb_errpush(4010,name,a_err='Allocate')
goto 9999
end if
desc_p%lprm(1) = 0
if (debug) then if (debug) then
write(0,*) me,'Done Convert_comm' write(0,*) me,'Done CDASB'
call psb_barrier(ictxt) call psb_barrier(ictxt)
end if end if
if (.false.) then
call psb_cdprt(70+me,desc_p,.false.)
end if
if (debug) write(0,*) me,'Done ConvertComm'
!!$ write(0,*) 't_halo_out ',allocated(t_halo_out)
!!$ Deallocate(works,workr,t_halo_in,work,&
!!$ & length_dl,dep_list,stat=info)
!!$ if (info /= 0) then
!!$ ch_err='Deallocate 1'
!!$ call psb_errpush(4013,name,a_err=ch_err,i_err=(/info,0,0,0,0/))
!!$ goto 9999
!!$ end if
!!$ Deallocate(tmp_ovr_idx,tmp_halo,&
!!$ & brvindx,rvsz,sdsz,bsdindx,temp,halo,stat=info)
!!$ if (info /= 0) then
!!$ ch_err='Deallocate 2'
!!$ call psb_errpush(4013,name,a_err=ch_err,i_err=(/info,0,0,0,0/))
!!$ goto 9999
!!$ end if
if (info == 0) call psb_sp_free(blk,info) if (info == 0) call psb_sp_free(blk,info)
if (info /= 0) then if (info /= 0) then
ch_err='sp_free' ch_err='sp_free'

@ -89,10 +89,10 @@ subroutine psb_dcsrp(trans,iperm,a, desc_a, info)
time(1) = mpi_wtime() 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) dectype = psb_cd_get_dectype(desc_a)
n_row = psb_cd_get_local_rows(desc_a) n_row = psb_cd_get_local_rows(desc_a)
n_col = psb_cd_get_local_cols(desc_a) n_col = psb_cd_get_local_cols(desc_a)
if(psb_get_errstatus() /= 0) return if(psb_get_errstatus() /= 0) return
info=0 info=0

@ -65,7 +65,7 @@ subroutine psb_dfree(x, desc_a, info)
goto 9999 goto 9999
end if end if
ictxt=psb_cd_get_context(desc_a) ictxt = psb_cd_get_context(desc_a)
call psb_info(ictxt, me, np) call psb_info(ictxt, me, np)
! ....verify blacs grid correctness.. ! ....verify blacs grid correctness..
@ -139,7 +139,7 @@ subroutine psb_dfreev(x, desc_a, info)
call psb_errpush(info,name) call psb_errpush(info,name)
return return
end if end if
ictxt=psb_cd_get_context(desc_a) ictxt = psb_cd_get_context(desc_a)
call psb_info(ictxt, me, np) call psb_info(ictxt, me, np)
if (np == -1) then if (np == -1) then

@ -85,8 +85,8 @@ subroutine psb_dgelp(trans,iperm,x,desc_a,info)
info=0 info=0
call psb_erractionsave(err_act) 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) dectype = psb_cd_get_dectype(desc_a)
nrow = psb_cd_get_local_rows(desc_a) nrow = psb_cd_get_local_rows(desc_a)
ncol = psb_cd_get_local_cols(desc_a) ncol = psb_cd_get_local_cols(desc_a)
i1sz = size(x,dim=1) i1sz = size(x,dim=1)
@ -95,7 +95,7 @@ subroutine psb_dgelp(trans,iperm,x,desc_a,info)
call psb_info(ictxt, me, np) call psb_info(ictxt, me, np)
if (debug) write(*,*) 'asb start: ',np,me,& if (debug) write(*,*) 'asb start: ',np,me,&
&psb_cd_get_dectype(desc_a) & psb_cd_get_dectype(desc_a)
! ....verify blacs grid correctness.. ! ....verify blacs grid correctness..
if (np == -1) then if (np == -1) then
info = 2010 info = 2010
@ -231,10 +231,10 @@ subroutine psb_dgelpv(trans,iperm,x,desc_a,info)
i1sz = size(x) i1sz = size(x)
ictxt=psb_cd_get_context(desc_a) ictxt = psb_cd_get_context(desc_a)
dectype=psb_cd_get_dectype(desc_a) dectype = psb_cd_get_dectype(desc_a)
nrow=psb_cd_get_local_rows(desc_a) nrow = psb_cd_get_local_rows(desc_a)
ncol=psb_cd_get_local_cols(desc_a) ncol = psb_cd_get_local_cols(desc_a)
call psb_info(ictxt, me, np) call psb_info(ictxt, me, np)

@ -35,8 +35,8 @@
! m - integer. Number of rows of submatrix belonging to ! m - integer. Number of rows of submatrix belonging to
! val to be inserted. ! val to be inserted.
! irw - integer(:) Row indices of rows of val (global numbering) ! irw - integer(:) Row indices of rows of val (global numbering)
! val - real, pointer, dimension(:). The source dense submatrix. ! val - real, dimension(:). The source dense submatrix.
! x - real, pointer, dimension(:). The destination dense matrix. ! x - real, dimension(:). The destination dense matrix.
! desc_a - type(<psb_desc_type>). The communication descriptor. ! desc_a - type(<psb_desc_type>). The communication descriptor.
! info - integer. Eventually returns an error code ! info - integer. Eventually returns an error code
subroutine psb_dinsvi(m, irw, val, x, desc_a, info, dupl) 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_const_mod
use psb_error_mod use psb_error_mod
use psb_penv_mod use psb_penv_mod
use psi_mod
implicit none implicit none
! m rows number of submatrix belonging to val to be inserted ! 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,& integer :: ictxt,i,loc_row,glob_row,&
& loc_rows,loc_cols,mglob,err_act, int_err(5) & loc_rows,loc_cols,mglob,err_act, int_err(5)
integer :: np, me, dupl_ integer :: np, me, dupl_
character(len=20) :: name integer, allocatable :: irl(:)
character(len=20) :: name
if(psb_get_errstatus() /= 0) return if(psb_get_errstatus() /= 0) return
info=0 info=0
call psb_erractionsave(err_act) call psb_erractionsave(err_act)
name = 'psb_dinsvi' name = 'psb_dinsvi'
if (.not.allocated(desc_a%glob_to_loc)) then !!$ if (.not.allocated(desc_a%glob_to_loc)) then
info=3110 !!$ info=3110
call psb_errpush(info,name) !!$ call psb_errpush(info,name)
return !!$ return
end if !!$ end if
if ((.not.allocated(desc_a%matrix_data))) then if ((.not.allocated(desc_a%matrix_data))) then
int_err(1)=3110 int_err(1)=3110
call psb_errpush(info,name) call psb_errpush(info,name)
return return
end if end if
ictxt=psb_cd_get_context(desc_a) ictxt = psb_cd_get_context(desc_a)
call psb_info(ictxt, me, np) call psb_info(ictxt, me, np)
if (np == -1) then if (np == -1) then
@ -111,31 +113,37 @@ subroutine psb_dinsvi(m, irw, val, x, desc_a, info, dupl)
goto 9999 goto 9999
endif endif
loc_rows=psb_cd_get_local_rows(desc_a) if (m==0) return
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) 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 if (present(dupl)) then
dupl_ = dupl dupl_ = dupl
else else
dupl_ = psb_dupl_ovwrt_ dupl_ = psb_dupl_ovwrt_
endif endif
call psi_idx_cnv(m,irw,irl,desc_a,info,owned=.true.)
select case(dupl_) select case(dupl_)
case(psb_dupl_ovwrt_) case(psb_dupl_ovwrt_)
do i = 1, m do i = 1, m
!loop over all val's rows !loop over all val's rows
! row actual block row ! row actual block row
glob_row=irw(i) if (irl(i) > 0) then
if ((glob_row>0).and.(glob_row <= mglob)) then ! this row belongs to me
! copy i-th row of block val in x
loc_row=desc_a%glob_to_loc(glob_row) x(irl(i)) = val(i)
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
end if end if
enddo enddo
@ -144,16 +152,10 @@ subroutine psb_dinsvi(m, irw, val, x, desc_a, info, dupl)
do i = 1, m do i = 1, m
!loop over all val's rows !loop over all val's rows
! row actual block row if (irl(i) > 0) then
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 ! this row belongs to me
! copy i-th row of block val in x ! copy i-th row of block val in x
x(loc_row) = x(loc_row) + val(i) x(irl(i)) = x(irl(i)) + val(i)
end if
end if end if
enddo enddo
@ -162,6 +164,7 @@ subroutine psb_dinsvi(m, irw, val, x, desc_a, info, dupl)
call psb_errpush(info,name) call psb_errpush(info,name)
goto 9999 goto 9999
end select end select
deallocate(irl)
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
return return
@ -215,8 +218,8 @@ end subroutine psb_dinsvi
! m - integer. Number of rows of submatrix belonging to ! m - integer. Number of rows of submatrix belonging to
! val to be inserted. ! val to be inserted.
! irw - integer(:) Row indices of rows of val (global numbering) ! irw - integer(:) Row indices of rows of val (global numbering)
! val - real, pointer, dimension(:,:). The source dense submatrix. ! val - real, dimension(:,:). The source dense submatrix.
! x - real, pointer, dimension(:,:). The destination dense matrix. ! x - real, dimension(:,:). The destination dense matrix.
! desc_a - type(<psb_desc_type>). The communication descriptor. ! desc_a - type(<psb_desc_type>). The communication descriptor.
! info - integer. Eventually returns an error code ! info - integer. Eventually returns an error code
subroutine psb_dinsi(m, irw, val, x, desc_a, info, dupl) 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_const_mod
use psb_error_mod use psb_error_mod
use psb_penv_mod use psb_penv_mod
use psi_mod
implicit none implicit none
! m rows number of submatrix belonging to val to be inserted ! 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 ! ix x global-row corresponding to position at which val submatrix
! must be inserted ! 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) :: m
integer, intent(in) :: irw(:) integer, intent(in) :: irw(:)
real(kind(1.d0)), intent(in) :: val(:,:) 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 type(psb_desc_type), intent(in) :: desc_a
integer, intent(out) :: info integer, intent(out) :: info
integer, optional, intent(in) :: dupl 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,& integer :: ictxt,i,loc_row,glob_row,j,n,&
& loc_rows,loc_cols,mglob,err_act, int_err(5) & 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 character(len=20) :: name
if(psb_get_errstatus() /= 0) return 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) call psb_erractionsave(err_act)
name = 'psb_dinsi' name = 'psb_dinsi'
if (.not.allocated(desc_a%glob_to_loc)) then !!$ if (.not.allocated(desc_a%glob_to_loc)) then
info=3110 !!$ info=3110
call psb_errpush(info,name) !!$ call psb_errpush(info,name)
return !!$ return
end if !!$ end if
if ((.not.allocated(desc_a%matrix_data))) then if ((.not.allocated(desc_a%matrix_data))) then
int_err(1)=3110 int_err(1)=3110
call psb_errpush(info,name) call psb_errpush(info,name)
return return
end if end if
ictxt=psb_cd_get_context(desc_a) ictxt = psb_cd_get_context(desc_a)
call psb_info(ictxt, me, np) call psb_info(ictxt, me, np)
if (np == -1) then 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) call psb_errpush(info,name,int_err)
goto 9999 goto 9999
endif endif
if (m==0) return
loc_rows=psb_cd_get_local_rows(desc_a) loc_rows = psb_cd_get_local_rows(desc_a)
loc_cols=psb_cd_get_local_cols(desc_a) loc_cols = psb_cd_get_local_cols(desc_a)
mglob = psb_cd_get_global_rows(desc_a) mglob = psb_cd_get_global_rows(desc_a)
n = min(size(val,2),size(x,2)) 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_ dupl_ = psb_dupl_ovwrt_
endif 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_) select case(dupl_)
case(psb_dupl_ovwrt_) case(psb_dupl_ovwrt_)
do i = 1, m do i = 1, m
!loop over all val's rows !loop over all val's rows
! row actual block row ! row actual block row
glob_row=irw(i) loc_row = irl(i)
if ((glob_row>0).and.(glob_row <= mglob)) then if (loc_row > 0) then
! this row belongs to me
loc_row=desc_a%glob_to_loc(glob_row) ! copy i-th row of block val in x
if (loc_row.ge.1) then do j=1,n
! this row belongs to me x(loc_row,j) = val(i,j)
! copy i-th row of block val in x end do
do j=1,n
x(loc_row,j) = val(i,j)
end do
end if
end if end if
enddo enddo
@ -332,17 +341,13 @@ subroutine psb_dinsi(m, irw, val, x, desc_a, info, dupl)
!loop over all val's rows !loop over all val's rows
! row actual block row ! row actual block row
glob_row=irw(i) loc_row = irl(i)
if ((glob_row>0).and.(glob_row <= mglob)) then if (loc_row > 0) then
! this row belongs to me
loc_row=desc_a%glob_to_loc(glob_row) ! copy i-th row of block val in x
if (loc_row.ge.1) then do j=1,n
! this row belongs to me x(loc_row,j) = x(loc_row,j) + val(i,j)
! copy i-th row of block val in x end do
do j=1,n
x(loc_row,j) = x(loc_row,j) + val(i,j)
end do
end if
end if end if
enddo enddo
@ -351,6 +356,7 @@ subroutine psb_dinsi(m, irw, val, x, desc_a, info, dupl)
call psb_errpush(info,name) call psb_errpush(info,name)
goto 9999 goto 9999
end select end select
deallocate(irl)
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
return return

@ -68,8 +68,9 @@ subroutine psb_dspalloc(a, desc_a, info, nnz)
call psb_erractionsave(err_act) call psb_erractionsave(err_act)
name = 'psb_dspalloc' name = 'psb_dspalloc'
ictxt = psb_cd_get_context(desc_a) ictxt = psb_cd_get_context(desc_a)
dectype = psb_cd_get_dectype(desc_a) dectype = psb_cd_get_dectype(desc_a)
call psb_info(ictxt, me, np) call psb_info(ictxt, me, np)
! ....verify blacs grid correctness.. ! ....verify blacs grid correctness..
if (np == -1) then if (np == -1) then
@ -127,7 +128,7 @@ subroutine psb_dspalloc(a, desc_a, info, nnz)
a%infoa(psb_state_) = psb_spmat_bld_ a%infoa(psb_state_) = psb_spmat_bld_
if (debug) write(0,*) 'spall: ', & if (debug) write(0,*) 'spall: ', &
&psb_cd_get_dectype(desc_a),psb_desc_bld_ & psb_cd_get_dectype(desc_a),psb_desc_bld_
return return

@ -76,7 +76,7 @@ subroutine psb_dspasb(a,desc_a, info, afmt, upd, dupl)
name = 'psb_spasb' name = 'psb_spasb'
call psb_erractionsave(err_act) 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) dscstate = psb_cd_get_dectype(desc_a)
n_row = psb_cd_get_local_rows(desc_a) n_row = psb_cd_get_local_rows(desc_a)
n_col = psb_cd_get_local_cols(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. ! First case: we come from a fresh build.
! !
n_row = psb_cd_get_local_rows(desc_a) n_row = psb_cd_get_local_rows(desc_a)
n_col = psb_cd_get_local_cols(desc_a) n_col = psb_cd_get_local_cols(desc_a)
! !
! Second step: handle the local matrix part. ! Second step: handle the local matrix part.
! !

@ -121,8 +121,7 @@ subroutine psb_dspcnv(a,b,desc_a,info)
time(1) = mpi_wtime() 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) dectype = psb_cd_get_dectype(desc_a)
n_row = psb_cd_get_local_rows(desc_a) n_row = psb_cd_get_local_rows(desc_a)
n_col = psb_cd_get_local_cols(desc_a) n_col = psb_cd_get_local_cols(desc_a)

@ -65,7 +65,7 @@ subroutine psb_dspfree(a, desc_a,info)
call psb_errpush(info,name) call psb_errpush(info,name)
return return
else else
ictxt=psb_cd_get_context(desc_a) ictxt = psb_cd_get_context(desc_a)
end if end if
!...deallocate a.... !...deallocate a....

@ -73,7 +73,7 @@ Subroutine psb_dsphalo(a,desc_a,blk,info,rwcnv,clcnv,outfmt)
& rvsz(:), bsdindx(:),sdsz(:) & rvsz(:), bsdindx(:),sdsz(:)
logical :: rwcnv_,clcnv_ logical :: rwcnv_,clcnv_
character(len=5) :: outfmt_ 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 real(kind(1.d0)) :: t1,t2,t3,t4,t5,t6,t7,t8,t9
character(len=20) :: name, ch_err character(len=20) :: name, ch_err
@ -100,7 +100,8 @@ Subroutine psb_dsphalo(a,desc_a,blk,info,rwcnv,clcnv,outfmt)
outfmt_ = 'CSR' outfmt_ = 'CSR'
endif endif
ictxt=psb_cd_get_context(desc_a) ictxt = psb_cd_get_context(desc_a)
Call psb_info(ictxt, me, np) Call psb_info(ictxt, me, np)
t1 = mpi_wtime() t1 = mpi_wtime()
@ -225,8 +226,6 @@ Subroutine psb_dsphalo(a,desc_a,blk,info,rwcnv,clcnv,outfmt)
counter = counter+n_el_send+3 counter = counter+n_el_send+3
Enddo Enddo
nz = tmp%infoa(psb_nnz_) 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 (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') 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) call psb_errpush(info,name,a_err=ch_err)
goto 9999 goto 9999
end if end if
!!$ call csprt(30+me,tmp,head='% SPHALO border SEND .') if (debugprt) then
!!$ close(30+me) 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,& 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 (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 (clcnv_) call psb_glob_to_loc(blk%ia2(1:iszr),desc_a,info,iact='I')
if (info /= 0) then if (info /= 0) then
info=4010 info=4010
ch_err='psbglob_to_loc' ch_err='psbglob_to_loc'
@ -268,6 +272,14 @@ Subroutine psb_dsphalo(a,desc_a,blk,info,rwcnv,clcnv,outfmt)
goto 9999 goto 9999
end if 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 l1 = 0
Do i=1,iszr Do i=1,iszr
!!$ write(0,*) work5(i),work6(i) !!$ write(0,*) work5(i),work6(i)
@ -283,9 +295,13 @@ Subroutine psb_dsphalo(a,desc_a,blk,info,rwcnv,clcnv,outfmt)
Enddo Enddo
blk%fida='COO' blk%fida='COO'
blk%infoa(psb_nnz_)=l1 blk%infoa(psb_nnz_)=l1
!!$ open(50+me) if (debugprt) then
!!$ call csprt(50+me,blk,head='% SPHALO border .') open(50+me)
!!$ close(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() t4 = mpi_wtime()
if(debug) Write(0,*)me,'End first loop',counter,l1,blk%m if(debug) Write(0,*)me,'End first loop',counter,l1,blk%m

@ -63,22 +63,34 @@ subroutine psb_dspins(nz,ia,ja,val,a,desc_a,info,rebuild)
logical, intent(in), optional :: rebuild logical, intent(in), optional :: rebuild
!locals..... !locals.....
integer :: nrow, err_act, dectype,mglob,ncol, spstate integer :: nrow, err_act,mglob,ncol, spstate
integer :: ictxt,np,me integer :: ictxt,np,me
logical, parameter :: debug=.false. logical, parameter :: debug=.false.
integer, parameter :: relocsz=200 integer, parameter :: relocsz=200
logical :: rebuild_ logical :: rebuild_
integer, allocatable :: ila(:),jla(:)
interface psb_cdins 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 use psb_descriptor_type
implicit none implicit none
type(psb_desc_type), intent(inout) :: desc_a type(psb_desc_type), intent(inout) :: desc_a
integer, intent(in) :: nz,ia(:),ja(:) integer, intent(in) :: nz,ia(:),ja(:)
integer, intent(out) :: info integer, intent(out) :: info
integer, optional, intent(out) :: ila(:), jla(:)
end subroutine psb_cdins end subroutine psb_cdins
end interface 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 character(len=20) :: name, ch_err
info = 0 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) 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) call psb_info(ictxt, me, np)
@ -98,7 +109,7 @@ subroutine psb_dspins(nz,ia,ja,val,a,desc_a,info,rebuild)
goto 9999 goto 9999
endif endif
if (nz <= 0) then if (nz < 0) then
info = 1111 info = 1111
call psb_errpush(info,name) call psb_errpush(info,name)
goto 9999 goto 9999
@ -119,6 +130,7 @@ subroutine psb_dspins(nz,ia,ja,val,a,desc_a,info,rebuild)
call psb_errpush(info,name) call psb_errpush(info,name)
goto 9999 goto 9999
end if end if
if (nz==0) return
if (present(rebuild)) then if (present(rebuild)) then
rebuild_ = rebuild rebuild_ = rebuild
@ -128,39 +140,104 @@ subroutine psb_dspins(nz,ia,ja,val,a,desc_a,info,rebuild)
spstate = a%infoa(psb_state_) spstate = a%infoa(psb_state_)
if (psb_is_bld_desc(desc_a)) then if (psb_is_bld_desc(desc_a)) then
call psb_cdins(nz,ia,ja,desc_a,info) if (psb_is_large_desc(desc_a)) then
if (info /= 0) then
info=4010 allocate(ila(nz),jla(nz),stat=info)
ch_err='psb_cdins' if (info /= 0) then
call psb_errpush(info,name,a_err=ch_err) ch_err='allocate'
goto 9999 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 end if
nrow = psb_cd_get_local_rows(desc_a) nrow = psb_cd_get_local_rows(desc_a)
ncol = psb_cd_get_local_cols(desc_a) ncol = psb_cd_get_local_cols(desc_a)
if (spstate == psb_spmat_bld_) then else if (psb_is_asb_desc(desc_a)) then
call psb_coins(nz,ia,ja,val,a,1,nrow,1,ncol,info,gtl=desc_a%glob_to_loc)
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 if (info /= 0) then
info=4010 info=4010
ch_err='psb_coins' ch_err='psb_coins'
call psb_errpush(info,name,a_err=ch_err) call psb_errpush(info,name,a_err=ch_err)
goto 9999 goto 9999
end if end if
else else
info = 1123 nrow = psb_cd_get_local_rows(desc_a)
call psb_errpush(info,name) ncol = psb_cd_get_local_cols(desc_a)
goto 9999 call psb_coins(nz,ia,ja,val,a,1,nrow,1,ncol,&
end if & info,gtl=desc_a%glob_to_loc,rebuild=rebuild_)
else if (psb_is_asb_desc(desc_a)) then if (info /= 0) then
nrow = psb_cd_get_local_rows(desc_a) info=4010
ncol = psb_cd_get_local_cols(desc_a) ch_err='psb_coins'
call psb_coins(nz,ia,ja,val,a,1,nrow,1,ncol,& call psb_errpush(info,name,a_err=ch_err)
& info,gtl=desc_a%glob_to_loc,rebuild=rebuild_) goto 9999
if (info /= 0) then end if
info=4010
ch_err='psb_coins'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if end if
else else
info = 1122 info = 1122

@ -70,6 +70,7 @@ Subroutine psb_dsprn(a, desc_a,info,clear)
call psb_erractionsave(err_act) call psb_erractionsave(err_act)
ictxt = psb_cd_get_context(desc_a) ictxt = psb_cd_get_context(desc_a)
call psb_info(ictxt, me, np) call psb_info(ictxt, me, np)
if (debug) & if (debug) &
&write(*,*) 'starting spalloc ',ictxt,np,me &write(*,*) 'starting spalloc ',ictxt,np,me

@ -46,6 +46,8 @@ subroutine psb_glob_to_loc2(x,y,desc_a,info,iact)
use psb_const_mod use psb_const_mod
use psb_error_mod use psb_error_mod
use psb_string_mod use psb_string_mod
use psb_penv_mod
use psi_mod
implicit none implicit none
!...parameters.... !...parameters....
@ -77,47 +79,23 @@ subroutine psb_glob_to_loc2(x,y,desc_a,info,iact)
int_err=0 int_err=0
real_val = 0.d0 real_val = 0.d0
n=size(x) n = size(x)
do i=1,n call psi_idx_cnv(n,x,y,desc_a,info)
if ((x(i).gt.psb_cd_get_global_rows(desc_a)).or.&
& (x(i).le.zero)) then select case(act)
if (act == 'I') then case('E','I')
y(i)=-3*psb_cd_get_global_rows(desc_a) call psb_erractionrestore(err_act)
else return
info=140 case('W')
int_err(1)=x(i) if ((info /= 0).or.(count(y(1:n)<0) >0)) then
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')
write(0,'("Error ",i5," in subroutine glob_to_loc")') info 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) call psb_errpush(info,name)
goto 9999 goto 9999
end select end if
endif end select
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
return return
@ -177,31 +155,37 @@ end subroutine psb_glob_to_loc2
! !
subroutine psb_glob_to_loc(x,desc_a,info,iact) subroutine psb_glob_to_loc(x,desc_a,info,iact)
use psb_penv_mod
use psb_descriptor_type use psb_descriptor_type
use psb_const_mod use psb_const_mod
use psb_error_mod use psb_error_mod
use psb_string_mod use psb_string_mod
use psi_mod
implicit none implicit none
!...parameters.... !...parameters....
type(psb_desc_type), intent(in) :: desc_a type(psb_desc_type), intent(in) :: desc_a
integer, intent(inout) :: x(:) integer, intent(inout) :: x(:)
integer, intent(out) :: info integer, intent(out) :: info
character, intent(in), optional :: iact character, intent(in), optional :: iact
!....locals.... !....locals....
integer :: n, i, tmp integer :: n, i, tmp, nk, key, idx, ih, nh, lb, ub, lm
character :: act character :: act
integer :: int_err(5), err_act integer :: int_err(5), err_act, dectype
real(kind(1.d0)) :: real_val real(kind(1.d0)) :: real_val, t0, t1,t2
integer, parameter :: zero=0 integer, parameter :: zero=0
character(len=20) :: name character(len=20) :: name
integer :: ictxt, iam, np
if(psb_get_errstatus() /= 0) return if(psb_get_errstatus() /= 0) return
info=0 info=0
name = 'glob_to_loc' name = 'glob_to_loc'
ictxt = desc_a%matrix_data(psb_ctxt_)
call psb_info(ictxt,iam,np)
call psb_erractionsave(err_act) call psb_erractionsave(err_act)
dectype = desc_a%matrix_data(psb_dec_type_)
if (present(iact)) then if (present(iact)) then
act=iact act=iact
else else
@ -210,48 +194,24 @@ subroutine psb_glob_to_loc(x,desc_a,info,iact)
act = toupper(act) act = toupper(act)
real_val = 0.d0 n = size(x)
n=size(x) call psi_idx_cnv(n,x,desc_a,info)
do i=1,n
if ((x(i).gt.psb_cd_get_global_rows(desc_a)).or.& select case(act)
& (x(i).le.zero)) then case('E','I')
if(act == 'I') then call psb_erractionrestore(err_act)
x(i)=-3*psb_cd_get_global_rows(desc_a) return
else case('W')
info=140 if ((info /= 0).or.(count(x(1:n)<0) >0)) then
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')
write(0,'("Error ",i5," in subroutine glob_to_loc")') info 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) call psb_errpush(info,name)
goto 9999 goto 9999
end select end if
endif end select
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
return return
@ -266,5 +226,69 @@ subroutine psb_glob_to_loc(x,desc_a,info,iact)
end if end if
return 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<glb_lc(lm,1)) then
ub = lm - 1
else
lb = lm + 1
end if
end do
else
tmp = -1
end if
if (tmp > 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 end subroutine psb_glob_to_loc

@ -55,9 +55,9 @@ subroutine psb_ialloc(x, desc_a, info, n)
!locals !locals
integer :: np,me,n_col,n_row,i,j,err_act 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) integer :: int_err(5), exch(3)
character(len=20) :: name, ch_err character(len=20) :: name, ch_err
if(psb_get_errstatus() /= 0) return if(psb_get_errstatus() /= 0) return
@ -74,7 +74,6 @@ subroutine psb_ialloc(x, desc_a, info, n)
goto 9999 goto 9999
endif endif
dectype=psb_cd_get_dectype(desc_a)
!... check m and n parameters.... !... check m and n parameters....
if (.not.psb_is_ok_desc(desc_a)) then if (.not.psb_is_ok_desc(desc_a)) then
info = 3110 info = 3110
@ -202,7 +201,7 @@ subroutine psb_iallocv(x, desc_a, info,n)
integer, optional, intent(in) :: n integer, optional, intent(in) :: n
!locals !locals
integer :: np,me,n_col,n_row,dectype,err_act integer :: np,me,n_col,n_row,err_act
integer :: ictxt, n_ integer :: ictxt, n_
integer :: int_err(5) integer :: int_err(5)
logical, parameter :: debug=.false. logical, parameter :: debug=.false.
@ -223,9 +222,6 @@ subroutine psb_iallocv(x, desc_a, info,n)
goto 9999 goto 9999
endif 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.... !... check m and n parameters....
if (.not.psb_is_ok_desc(desc_a)) then if (.not.psb_is_ok_desc(desc_a)) then
info = 3110 info = 3110
@ -237,23 +233,23 @@ subroutine psb_iallocv(x, desc_a, info,n)
!....allocate x ..... !....allocate x .....
if (psb_is_asb_desc(desc_a).or.psb_is_upd_desc(desc_a)) then 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)) n_col = max(1,psb_cd_get_local_cols(desc_a))
allocate(x(n_col),stat=info) allocate(x(n_col),stat=info)
if (info.ne.0) then if (info.ne.0) then
info=2025 info=2025
int_err(1)=n_col int_err(1)=n_col
call psb_errpush(info,name,int_err) call psb_errpush(info,name,int_err)
goto 9999 goto 9999
endif endif
else if (psb_is_bld_desc(desc_a)) then else if (psb_is_bld_desc(desc_a)) then
n_row = max(1,psb_cd_get_local_rows(desc_a)) n_row = max(1,psb_cd_get_local_rows(desc_a))
allocate(x(n_row),stat=info) allocate(x(n_row),stat=info)
if (info.ne.0) then if (info.ne.0) then
info=2025 info=2025
int_err(1)=n_row int_err(1)=n_row
call psb_errpush(info,name,int_err) call psb_errpush(info,name,int_err)
goto 9999 goto 9999
endif endif
endif endif
x = 0 x = 0

@ -53,7 +53,7 @@ subroutine psb_iasb(x, desc_a, info)
! local variables ! local variables
integer :: ictxt,np,me,nrow,ncol,err_act 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. logical, parameter :: debug=.false.
character(len=20) :: name character(len=20) :: name
@ -68,14 +68,13 @@ subroutine psb_iasb(x, desc_a, info)
return return
endif endif
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) call psb_info(ictxt, me, np)
if (debug) write(*,*) 'asb start: ',np,me,& if (debug) write(*,*) 'asb start: ',np,me,&
&psb_cd_get_dectype(desc_a) & psb_cd_get_dectype(desc_a)
! ....verify blacs grid correctness.. ! ....verify blacs grid correctness..
if (np == -1) then if (np == -1) then
info = 2010 info = 2010
@ -83,16 +82,16 @@ subroutine psb_iasb(x, desc_a, info)
goto 9999 goto 9999
else if (.not.psb_is_asb_desc(desc_a)) then else if (.not.psb_is_asb_desc(desc_a)) then
if (debug) write(*,*) 'asb error ',& if (debug) write(*,*) 'asb error ',&
&dectype & psb_cd_get_dectype(desc_a)
info = 3110 info = 3110
call psb_errpush(info,name) call psb_errpush(info,name)
goto 9999 goto 9999
endif endif
! check size ! check size
ictxt=psb_cd_get_context(desc_a) ictxt = psb_cd_get_context(desc_a)
nrow=psb_cd_get_local_rows(desc_a) nrow = psb_cd_get_local_rows(desc_a)
ncol=psb_cd_get_local_cols(desc_a) ncol = psb_cd_get_local_cols(desc_a)
i1sz = size(x,dim=1) i1sz = size(x,dim=1)
i2sz = size(x,dim=2) i2sz = size(x,dim=2)
if (debug) write(*,*) 'asb: ',i1sz,i2sz,nrow,ncol if (debug) write(*,*) 'asb: ',i1sz,i2sz,nrow,ncol
@ -177,7 +176,7 @@ subroutine psb_iasbv(x, desc_a, info)
! local variables ! local variables
integer :: ictxt,np,me, err_act 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. logical, parameter :: debug=.false.
character(len=20) :: name character(len=20) :: name
@ -187,8 +186,7 @@ subroutine psb_iasbv(x, desc_a, info)
name = 'psb_iasbv' name = 'psb_iasbv'
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) call psb_info(ictxt, me, np)
@ -203,8 +201,9 @@ subroutine psb_iasbv(x, desc_a, info)
goto 9999 goto 9999
endif endif
nrow=psb_cd_get_local_rows(desc_a) nrow = psb_cd_get_local_rows(desc_a)
ncol=psb_cd_get_local_cols(desc_a) ncol = psb_cd_get_local_cols(desc_a)
if (debug) write(*,*) name,' sizes: ',nrow,ncol if (debug) write(*,*) name,' sizes: ',nrow,ncol
i1sz = size(x) i1sz = size(x)
if (debug) write(*,*) 'dasb: sizes ',i1sz,ncol if (debug) write(*,*) 'dasb: sizes ',i1sz,ncol

@ -36,7 +36,7 @@
! val to be inserted. ! val to be inserted.
! irw - integer(:) Row indices of rows of val (global numbering) ! irw - integer(:) Row indices of rows of val (global numbering)
! val - integer, dimension(:). The source dense submatrix. ! 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(<psb_desc_type>). The communication descriptor. ! desc_a - type(<psb_desc_type>). The communication descriptor.
! info - integer. Eventually returns an error code ! info - integer. Eventually returns an error code
subroutine psb_iinsvi(m, irw, val, x, desc_a, info, dupl) 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_const_mod
use psb_error_mod use psb_error_mod
use psb_penv_mod use psb_penv_mod
use psi_mod
implicit none implicit none
! m rows number of submatrix belonging to val to be inserted ! 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..... !locals.....
integer :: ictxt,i,loc_row,glob_row,& integer :: ictxt,i,loc_row,glob_row,&
& loc_rows,loc_cols,mglob,err_act, int_err(5), err & loc_rows,loc_cols,mglob,err_act, int_err(5)
integer :: np, me,dupl_ integer :: np, me, dupl_
character(len=20) :: name integer, allocatable :: irl(:)
character(len=20) :: name
if(psb_get_errstatus() /= 0) return if(psb_get_errstatus() /= 0) return
info=0 info=0
call psb_erractionsave(err_act) call psb_erractionsave(err_act)
name = 'psb_dinsvv' name = 'psb_insvi'
if (.not.allocated(desc_a%glob_to_loc)) then !!$ if (.not.allocated(desc_a%glob_to_loc)) then
info=3110 !!$ info=3110
call psb_errpush(info,name) !!$ call psb_errpush(info,name)
return !!$ return
end if !!$ end if
if ((.not.allocated(desc_a%matrix_data))) then if ((.not.allocated(desc_a%matrix_data))) then
int_err(1)=3110 int_err(1)=3110
call psb_errpush(info,name) call psb_errpush(info,name)
@ -111,31 +113,36 @@ subroutine psb_iinsvi(m, irw, val, x, desc_a, info, dupl)
goto 9999 goto 9999
endif endif
loc_rows=psb_cd_get_local_rows(desc_a) if (m==0) return
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) 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 if (present(dupl)) then
dupl_ = dupl dupl_ = dupl
else else
dupl_ = psb_dupl_ovwrt_ dupl_ = psb_dupl_ovwrt_
endif endif
call psi_idx_cnv(m,irw,irl,desc_a,info,owned=.true.)
select case(dupl_) select case(dupl_)
case(psb_dupl_ovwrt_) case(psb_dupl_ovwrt_)
do i = 1, m do i = 1, m
!loop over all val's rows !loop over all val's rows
! row actual block row ! row actual block row
glob_row=irw(i) if (irl(i) > 0) then
if ((glob_row>0).and.(glob_row <= mglob)) then ! this row belongs to me
! copy i-th row of block val in x
loc_row=desc_a%glob_to_loc(glob_row) x(irl(i)) = val(i)
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
end if end if
enddo enddo
@ -144,16 +151,10 @@ subroutine psb_iinsvi(m, irw, val, x, desc_a, info, dupl)
do i = 1, m do i = 1, m
!loop over all val's rows !loop over all val's rows
! row actual block row if (irl(i) > 0) then
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 ! this row belongs to me
! copy i-th row of block val in x ! copy i-th row of block val in x
x(loc_row) = x(loc_row) + val(i) x(irl(i)) = x(irl(i)) + val(i)
end if
end if end if
enddo enddo
@ -162,6 +163,7 @@ subroutine psb_iinsvi(m, irw, val, x, desc_a, info, dupl)
call psb_errpush(info,name) call psb_errpush(info,name)
goto 9999 goto 9999
end select end select
deallocate(irl)
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
return return
@ -216,7 +218,7 @@ end subroutine psb_iinsvi
! val to be inserted. ! val to be inserted.
! irw - integer(:) Row indices of rows of val (global numbering) ! irw - integer(:) Row indices of rows of val (global numbering)
! val - integer, dimension(:,:). The source dense submatrix. ! 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(<psb_desc_type>). The communication descriptor. ! desc_a - type(<psb_desc_type>). The communication descriptor.
! info - integer. Eventually returns an error code ! info - integer. Eventually returns an error code
subroutine psb_iinsi(m,irw, val, x, desc_a, info, dupl) 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_const_mod
use psb_error_mod use psb_error_mod
use psb_penv_mod use psb_penv_mod
use psi_mod
implicit none implicit none
! m rows number of submatrix belonging to val to be inserted ! 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..... !locals.....
integer :: ictxt,i,loc_row,glob_row,j,n,& integer :: ictxt,i,loc_row,glob_row,j,n,&
& loc_rows,loc_cols,mglob,err_act, int_err(5) & 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 character(len=20) :: name
if(psb_get_errstatus() /= 0) return if(psb_get_errstatus() /= 0) return
info=0 info=0
call psb_erractionsave(err_act) call psb_erractionsave(err_act)
name = 'psb_dinsvv' name = 'psb_iinsi'
if (.not.allocated(desc_a%glob_to_loc)) then !!$ if (.not.allocated(desc_a%glob_to_loc)) then
info=3110 !!$ info=3110
call psb_errpush(info,name) !!$ call psb_errpush(info,name)
return !!$ return
end if !!$ end if
if ((.not.allocated(desc_a%matrix_data))) then if ((.not.allocated(desc_a%matrix_data))) then
int_err(1)=3110 int_err(1)=3110
call psb_errpush(info,name) 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) call psb_errpush(info,name,int_err)
goto 9999 goto 9999
endif endif
if (m==0) return
loc_rows=psb_cd_get_local_rows(desc_a) loc_rows = psb_cd_get_local_rows(desc_a)
loc_cols=psb_cd_get_local_cols(desc_a) loc_cols = psb_cd_get_local_cols(desc_a)
mglob = psb_cd_get_global_rows(desc_a) mglob = psb_cd_get_global_rows(desc_a)
n = min(size(val,2),size(x,2)) 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_ dupl_ = psb_dupl_ovwrt_
endif 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_) select case(dupl_)
case(psb_dupl_ovwrt_) case(psb_dupl_ovwrt_)
do i = 1, m do i = 1, m
!loop over all val's rows !loop over all val's rows
! row actual block row ! row actual block row
glob_row=irw(i) loc_row = irl(i)
if ((glob_row>0).and.(glob_row <= mglob)) then if (loc_row > 0) then
! this row belongs to me
loc_row=desc_a%glob_to_loc(glob_row) ! copy i-th row of block val in x
if (loc_row.ge.1) then do j=1,n
! this row belongs to me x(loc_row,j) = val(i,j)
! copy i-th row of block val in x end do
do j=1,n
x(loc_row,j) = val(i,j)
end do
end if
end if end if
enddo enddo
@ -330,17 +340,13 @@ subroutine psb_iinsi(m,irw, val, x, desc_a, info, dupl)
!loop over all val's rows !loop over all val's rows
! row actual block row ! row actual block row
glob_row=irw(i) loc_row = irl(i)
if ((glob_row>0).and.(glob_row <= mglob)) then if (loc_row > 0) then
! this row belongs to me
loc_row=desc_a%glob_to_loc(glob_row) ! copy i-th row of block val in x
if (loc_row.ge.1) then do j=1,n
! this row belongs to me x(loc_row,j) = x(loc_row,j) + val(i,j)
! copy i-th row of block val in x end do
do j=1,n
x(loc_row,j) = x(loc_row,j) + val(i,j)
end do
end if
end if end if
enddo enddo
@ -349,6 +355,7 @@ subroutine psb_iinsi(m,irw, val, x, desc_a, info, dupl)
call psb_errpush(info,name) call psb_errpush(info,name)
goto 9999 goto 9999
end select end select
deallocate(irl)
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
return return

@ -56,7 +56,7 @@ subroutine psb_zalloc(x, desc_a, info, n)
!locals !locals
integer :: np,me,err,n_col,n_row,i,j,err_act 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) integer :: int_err(5),exch(3)
character(len=20) :: name, ch_err character(len=20) :: name, ch_err
@ -76,7 +76,6 @@ subroutine psb_zalloc(x, desc_a, info, n)
goto 9999 goto 9999
endif endif
dectype=psb_cd_get_dectype(desc_a)
!... check m and n parameters.... !... check m and n parameters....
if (.not.psb_is_ok_desc(desc_a)) then if (.not.psb_is_ok_desc(desc_a)) then
info = 3110 info = 3110
@ -202,7 +201,7 @@ subroutine psb_zallocv(x, desc_a,info,n)
integer, optional, intent(in) :: n integer, optional, intent(in) :: n
!locals !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_ integer :: ictxt, n_
logical, parameter :: debug=.false. logical, parameter :: debug=.false.
character(len=20) :: name, ch_err character(len=20) :: name, ch_err
@ -222,9 +221,7 @@ subroutine psb_zallocv(x, desc_a,info,n)
goto 9999 goto 9999
endif endif
dectype=psb_cd_get_dectype(desc_a) if (debug) write(0,*) 'dall: is_ok?',psb_is_ok_desc(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.... !... check m and n parameters....
if (.not.psb_is_ok_desc(desc_a)) then if (.not.psb_is_ok_desc(desc_a)) then
info = 3110 info = 3110

@ -53,7 +53,7 @@ subroutine psb_zasb(x, desc_a, info)
! local variables ! local variables
integer :: ictxt,np,me,nrow,ncol, err_act 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. logical, parameter :: debug=.false.
character(len=20) :: name, ch_err character(len=20) :: name, ch_err
@ -69,13 +69,11 @@ subroutine psb_zasb(x, desc_a, info)
endif endif
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) call psb_info(ictxt, me, np)
if (debug) write(*,*) 'asb start: ',np,me,& if (debug) write(*,*) 'asb start: ',np,me,&
&psb_cd_get_dectype(desc_a) & psb_cd_get_dectype(desc_a)
! ....verify blacs grid correctness.. ! ....verify blacs grid correctness..
if (np == -1) then if (np == -1) then
info = 2010 info = 2010
@ -83,16 +81,16 @@ subroutine psb_zasb(x, desc_a, info)
goto 9999 goto 9999
else if (.not.psb_is_asb_desc(desc_a)) then else if (.not.psb_is_asb_desc(desc_a)) then
if (debug) write(*,*) 'asb error ',& if (debug) write(*,*) 'asb error ',&
&dectype & psb_cd_get_dectype(desc_a)
info = 3110 info = 3110
call psb_errpush(info,name,i_err=int_err) call psb_errpush(info,name,i_err=int_err)
goto 9999 goto 9999
endif endif
! check size ! check size
ictxt=psb_cd_get_context(desc_a) ictxt = psb_cd_get_context(desc_a)
nrow=psb_cd_get_local_rows(desc_a) nrow = psb_cd_get_local_rows(desc_a)
ncol=psb_cd_get_local_cols(desc_a) ncol = psb_cd_get_local_cols(desc_a)
i1sz = size(x,dim=1) i1sz = size(x,dim=1)
i2sz = size(x,dim=2) i2sz = size(x,dim=2)
if (debug) write(*,*) 'asb: ',i1sz,i2sz,nrow,ncol if (debug) write(*,*) 'asb: ',i1sz,i2sz,nrow,ncol
@ -182,7 +180,7 @@ subroutine psb_zasbv(x, desc_a, info)
! local variables ! local variables
integer :: ictxt,np,me 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. logical, parameter :: debug=.false.
character(len=20) :: name,ch_err character(len=20) :: name,ch_err
@ -192,7 +190,6 @@ subroutine psb_zasbv(x, desc_a, info)
name = 'psb_zasbv' name = 'psb_zasbv'
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) call psb_info(ictxt, me, np)

@ -51,6 +51,10 @@ Subroutine psb_zcdovr(a,desc_a,novr,desc_ov,info)
Use psb_prec_mod Use psb_prec_mod
use psb_error_mod use psb_error_mod
use psb_penv_mod use psb_penv_mod
use psb_tools_mod
use psb_realloc_mod
use psi_mod
use mpi
Implicit None Implicit None
! .. Array Arguments .. ! .. Array Arguments ..
@ -60,42 +64,25 @@ Subroutine psb_zcdovr(a,desc_a,novr,desc_ov,info)
Type(psb_desc_type), Intent(inout) :: desc_ov Type(psb_desc_type), Intent(inout) :: desc_ov
integer, intent(out) :: info integer, intent(out) :: info
real(kind(1.d0)) :: t1,t2,t3,mpi_wtime
external mpi_wtime
integer icomm, err_act 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 .. ! .. Local Scalars ..
Integer :: i, j, k, np, me,m,nnzero,& 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 & index_dim,elem_dim, l_tmp_ovr_idx,l_tmp_halo, nztot,nhalo
Logical, parameter :: debug=.false. Integer :: counter,counter_h, counter_o, counter_e,idx,gidx,proc,n_elem_recv,&
character(len=20) :: name, ch_err & 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' name='psb_cdovr'
info = 0 info = 0
@ -107,10 +94,12 @@ Subroutine psb_zcdovr(a,desc_a,novr,desc_ov,info)
If(debug) Write(0,*)'in psb_cdovr',novr If(debug) Write(0,*)'in psb_cdovr',novr
m=psb_cd_get_local_rows(desc_a) m = psb_cd_get_local_rows(desc_a)
nnzero=Size(a%aspk) nnzero = Size(a%aspk)
n_col=psb_cd_get_local_cols(desc_a) n_row = psb_cd_get_local_rows(desc_a)
nhalo = n_col-m 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(debug) Write(0,*)'IN CDOVR1',novr ,m,nnzero,n_col
if (novr<0) then if (novr<0) then
info=10 info=10
@ -120,19 +109,19 @@ Subroutine psb_zcdovr(a,desc_a,novr,desc_ov,info)
goto 9999 goto 9999
endif 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 if (novr==0) then
! !
! Just copy the input. ! 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 return
endif endif
@ -142,12 +131,10 @@ Subroutine psb_zcdovr(a,desc_a,novr,desc_ov,info)
Write(0,*)'BEGIN cdovr',me,nhalo Write(0,*)'BEGIN cdovr',me,nhalo
call psb_barrier(ictxt) call psb_barrier(ictxt)
endif endif
t1 = mpi_wtime()
! !
! Ok, since we are only estimating, do it as follows: ! 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. ! nonzeros per row is the same as the global.
! !
nztot = psb_sp_get_nnzeros(a) 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 lworks = ((nztot+m-1)/m)*nhalo
lworkr = ((nztot+m-1)/m)*nhalo lworkr = ((nztot+m-1)/m)*nhalo
else else
info=-1 info=-1
call psb_errpush(info,name) call psb_errpush(info,name)
goto 9999 goto 9999
endif endif
If(debug)Write(0,*)'ovr_est done',me,novr,lovr If(debug)Write(0,*)'ovr_est done',me,novr,lovr
index_dim = size(desc_a%halo_index) index_dim = size(desc_a%halo_index)
elem_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_ovr_idx = novr*(3*Max(2*index_dim,1)+1)
l_tmp_halo=novr*(3*Size(desc_a%halo_index)) l_tmp_halo = novr*(3*Size(desc_a%halo_index))
desc_ov%matrix_data(:) = desc_a%matrix_data(:) if (psb_is_large_desc(desc_a)) then
desc_ov%matrix_data(psb_dec_type_) = psb_desc_bld_ desc_ov%matrix_data(psb_dec_type_) = psb_desc_large_bld_
else
Allocate(desc_ov%loc_to_glob(Size(desc_a%loc_to_glob)),& desc_ov%matrix_data(psb_dec_type_) = psb_desc_bld_
& 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
end if 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 If(debug) then
Write(0,*)'Start cdovrbld',me,lworks,lworkr Write(0,*)'Start cdovrbld',me,lworks,lworkr
call psb_barrier(ictxt) call psb_barrier(ictxt)
endif 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 if (.false.) then
info=4010 !
ch_err='psb_cdovrbld' ! The real work goes on in here....
call psb_errpush(info,name,a_err=ch_err) !
goto 9999 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<J
!
proc_id = temp(i)
call psb_check_size((counter_t+3),t_halo_in,info,pad=-1)
if (info /= 0) then
info=4010
call psb_errpush(info,name,a_err='psb_check_size')
goto 9999
end if
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_ov%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<J
!
n_col=n_col+1
proc_id=-desc_ov%glob_to_loc(idx)-np-1
call psb_check_size(n_col,desc_ov%loc_to_glob,info,pad=-1)
if (info /= 0) then
info=4010
call psb_errpush(info,name,a_err='psb_check_size')
goto 9999
end if
desc_ov%glob_to_loc(idx)=n_col
desc_ov%loc_to_glob(n_col)=idx
call psb_check_size((counter_t+3),t_halo_in,info,pad=-1)
if (info /= 0) then
info=4010
call psb_errpush(info,name,a_err='psb_check_size')
goto 9999
end if
t_halo_in(counter_t)=proc_id
t_halo_in(counter_t+1)=1
t_halo_in(counter_t+2)=n_col
counter_t=counter_t+3
if (debug) write(0,*) me,' CDOVRBLD: Added into t_halo_in from recv',&
&proc_id,n_col,idx
else if (desc_ov%glob_to_loc(idx) < 0) Then
if (debug) write(0,*) me,'Wrong input to cdovrbld??',&
&idx,desc_ov%glob_to_loc(idx)
End If
End Do
desc_ov%matrix_data(psb_n_col_)=n_col
end if
end if
!!$ desc_ov%matrix_data(psb_n_row_)=desc_ov%matrix_data(psb_n_col_)
!
! Ok, now we have a temporary halo with all the info for the
! next round. If we need to keep going, convert the halo format
! from temporary to final, so that we can work out the next iteration.
! This uses one of the convert_comm internals, i.e. we are doing
! the equivalent of a partial call to convert_comm
!
If (i_ovr < (novr)) Then
t_halo_in(counter_t)=-1
if (debug) write(0,*) me,'Checktmp_o_i 1',tmp_ovr_idx(1:10)
if (debug) write(0,*) me,'Calling Crea_Halo'
call psi_crea_index(desc_ov,t_halo_in,t_halo_out,.false.,&
& nxch,nsnd,nrcv,info)
if (debug) then
write(0,*) me,'Done Crea_Index'
call psb_barrier(ictxt)
end if
if (debug) write(0,*) me,'Checktmp_o_i 2',tmp_ovr_idx(1:10)
if (debug) write(0,*) me,'Done Crea_Halo'
call psb_transfer(t_halo_out,halo,info)
!
! At this point we have built the halo necessary for I_OVR+1.
!
End If
if (debug) write(0,*) me,'Checktmp_o_i ',tmp_ovr_idx(1:10)
End Do
desc_ov%matrix_data(psb_m_)=psb_cd_get_global_rows(desc_a)
desc_ov%matrix_data(psb_n_)=psb_cd_get_global_cols(desc_a)
tmp_halo(counter_h:)=-1
tmp_ovr_idx(counter_o:)=-1
!
! At this point we have gathered all the indices in the halo at
! N levels of overlap. Just call cnv_dsc. This is
! the same routine as gets called inside CDASB.
!
if (debug) then
write(0,*) 'psb_cdovrbld: converting indexes'
call psb_barrier(ictxt)
end if
!.... convert comunication stuctures....
! Note that we have to keep local_rows until the very end,
! because otherwise the halo build mechanism of cdasb
! will not work.
call psb_transfer(tmp_ovr_idx,desc_ov%ovrlap_index,info)
call psb_transfer(tmp_halo,desc_ov%halo_index,info)
call psb_cdasb(desc_ov,info)
desc_ov%matrix_data(psb_n_row_)=desc_ov%matrix_data(psb_n_col_)
if (debug) then
write(0,*) me,'Done CDASB'
call psb_barrier(ictxt)
end if
if (info == 0) call psb_sp_free(blk,info)
if (info /= 0) then
ch_err='sp_free'
call psb_errpush(4013,name,a_err=ch_err,i_err=(/info,0,0,0,0/))
goto 9999
end if
end if end if
desc_ov%matrix_data(psb_dec_type_) = psb_desc_asb_
If(debug) then
Write(0,*)'Done cdovrbld',me,lworks,lworkr
call psb_barrier(ictxt)
endif
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
return return
@ -216,8 +701,8 @@ Subroutine psb_zcdovr(a,desc_a,novr,desc_ov,info)
9999 continue 9999 continue
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
if (err_act == act_abort) then if (err_act == act_abort) then
call psb_error(ictxt) call psb_error(ictxt)
return return
end if end if
Return Return

@ -54,34 +54,30 @@ Subroutine psb_zcdovrbld(n_ovr,desc_p,desc_a,a,&
use psb_serial_mod use psb_serial_mod
Use psi_mod Use psi_mod
use psb_realloc_mod use psb_realloc_mod
use psb_tools_mod, only : psb_cdprt
use psb_error_mod use psb_error_mod
use psb_const_mod use psb_const_mod
use psb_penv_mod use psb_penv_mod
use psb_tools_mod
use mpi
Implicit None Implicit None
include 'mpif.h'
type(psb_zspmat_type),intent(in) :: a type(psb_zspmat_type),intent(in) :: a
type(psb_desc_type),intent(in) :: desc_a type(psb_desc_type),intent(in) :: desc_a
type(psb_desc_type),intent(inout) :: desc_p 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??? ! 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(inout) :: lworks, lworkr
integer, intent(out) :: info integer, intent(out) :: info
type(psb_zspmat_type) :: blk type(psb_zspmat_type) :: blk
Integer, allocatable :: tmp_halo(:),tmp_ovr_idx(:) Integer, allocatable :: tmp_halo(:),tmp_ovr_idx(:)
Integer :: counter,counter_h, counter_o, counter_e,j,idx,gidx,proc,n_elem_recv,& 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,& & 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, & & 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 & idxr, idxs, lx, iszr, iszs, err_act, icomm, nxch, nsnd, nrcv,lidx
Integer,allocatable :: halo(:),works(:),workr(:),t_halo_in(:),&
Integer,allocatable :: halo(:),length_dl(:),works(:),workr(:),t_halo_in(:),& & t_halo_out(:),temp(:),maskr(:)
& t_halo_out(:),work(:),dep_list(:),temp(:)
Integer,allocatable :: brvindx(:),rvsz(:), bsdindx(:),sdsz(:) Integer,allocatable :: brvindx(:),rvsz(:), bsdindx(:),sdsz(:)
Logical,Parameter :: debug=.false. 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) n_col = psb_cd_get_local_cols(desc_a)
if (debug) write(0,*) me,' On entry to CDOVRBLD n_col:',n_col 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),& Allocate(works(lworks),workr(lworkr),t_halo_in(l_tmp_halo),&
& t_halo_out(l_tmp_halo), work(lwork),& & t_halo_out(l_tmp_halo), temp(lworkr),stat=info)
& length_dl(np+1),dep_list(dl_lda*np),temp(lworkr),stat=info)
if (info /= 0) then if (info /= 0) then
call psb_errpush(4010,name,a_err='Allocate') call psb_errpush(4010,name,a_err='Allocate')
goto 9999 goto 9999
@ -141,7 +134,7 @@ Subroutine psb_zcdovrbld(n_ovr,desc_p,desc_a,a,&
call psb_errpush(4010,name,a_err='Allocate') call psb_errpush(4010,name,a_err='Allocate')
goto 9999 goto 9999
end if end if
halo(:) = desc_a%halo_index(:) halo(:) = desc_a%halo_index(:)
desc_p%ovrlap_elem(:) = -1 desc_p%ovrlap_elem(:) = -1
tmp_ovr_idx(:) = -1 tmp_ovr_idx(:) = -1
tmp_halo(:) = -1 tmp_halo(:) = -1
@ -150,7 +143,40 @@ Subroutine psb_zcdovrbld(n_ovr,desc_p,desc_a,a,&
counter_h = 1 counter_h = 1
counter_o = 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. ! 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 Do i_ovr=1,n_ovr
if (debug) write(0,*) me,'Running on overlap level ',i_ovr,' of ',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 ! At this point, halo contains a valid halo corresponding to the
@ -193,6 +218,7 @@ Subroutine psb_zcdovrbld(n_ovr,desc_p,desc_a,a,&
counter = 1 counter = 1
counter_t = 1 counter_t = 1
t1 = mpi_wtime() t1 = mpi_wtime()
Do While (halo(counter) /= -1) Do While (halo(counter) /= -1)
tot_elem=0 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_recv=halo(counter+psb_n_elem_recv_)
n_elem_send=halo(counter+n_elem_recv+psb_n_elem_send_) n_elem_send=halo(counter+n_elem_recv+psb_n_elem_send_)
If ((counter+n_elem_recv+n_elem_send) > Size(halo)) then If ((counter+n_elem_recv+n_elem_send) > Size(halo)) then
info=-1 info = -1
call psb_errpush(info,name) call psb_errpush(info,name)
goto 9999 goto 9999
end If 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 ! 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 ! 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 ! be enlarged with the new column indices received, and will reassemble
! everything for the next iteration. ! everything for the next iteration.
! !
@ -225,7 +251,7 @@ Subroutine psb_zcdovrbld(n_ovr,desc_p,desc_a,a,&
call psb_errpush(info,name) call psb_errpush(info,name)
goto 9999 goto 9999
end If end If
idx = halo(counter+psb_elem_recv_+j)
idx = halo(counter+psb_elem_recv_+j) idx = halo(counter+psb_elem_recv_+j)
If(idx > Size(desc_p%loc_to_glob)) then If(idx > Size(desc_p%loc_to_glob)) then
info=-3 info=-3
@ -235,42 +261,33 @@ Subroutine psb_zcdovrbld(n_ovr,desc_p,desc_a,a,&
gidx = desc_p%loc_to_glob(idx) gidx = desc_p%loc_to_glob(idx)
If((counter_o+2) > Size(tmp_ovr_idx)) Then call psb_check_size((counter_o+3),tmp_ovr_idx,info,pad=-1)
isz = max((3*Size(tmp_ovr_idx))/2,(counter_o+3)) if (info /= 0) then
if (debug) write(0,*) me,'Realloc tmp_ovr',isz info=4010
call psb_realloc(isz,tmp_ovr_idx,info,pad=-1) call psb_errpush(info,name,a_err='psb_check_size')
if (info /= 0) then goto 9999
info=4010 end if
ch_err='psb_realloc'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
End If
tmp_ovr_idx(counter_o)=proc tmp_ovr_idx(counter_o)=proc
tmp_ovr_idx(counter_o+1)=1 tmp_ovr_idx(counter_o+1)=1
tmp_ovr_idx(counter_o+2)=gidx tmp_ovr_idx(counter_o+2)=gidx
tmp_ovr_idx(counter_o+3)=-1 tmp_ovr_idx(counter_o+3)=-1
counter_o=counter_o+3 counter_o=counter_o+3
if (.not.psb_is_large_desc(desc_p)) then
If((counter_h+2) > Size(tmp_halo)) Then call psb_check_size((counter_h+3),tmp_halo,info,pad=-1)
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 (info /= 0) then if (info /= 0) then
info=4010 info=4010
ch_err='psb_realloc' call psb_errpush(info,name,a_err='psb_check_size')
call psb_errpush(info,name,a_err=ch_err)
goto 9999 goto 9999
end if end if
End If
tmp_halo(counter_h)=proc tmp_halo(counter_h)=proc
tmp_halo(counter_h+1)=1 tmp_halo(counter_h+1)=1
tmp_halo(counter_h+2)=idx tmp_halo(counter_h+2)=idx
tmp_halo(counter_h+3)=-1 tmp_halo(counter_h+3)=-1
counter_h=counter_h+3 counter_h=counter_h+3
end if
Enddo Enddo
if (debug) write(0,*) me,'Checktmp_o_i Loop Mid1',tmp_ovr_idx(1:10) 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) idx = halo(counter+psb_elem_send_+j)
gidx = desc_p%loc_to_glob(idx) gidx = desc_p%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)
If((counter_o+2) > Size(tmp_ovr_idx)) Then call psb_check_size((counter_o+3),tmp_ovr_idx,info,pad=-1)
isz = max((3*Size(tmp_ovr_idx))/2,(counter_o+3)) if (info /= 0) then
if (debug) write(0,*) me,'Realloc tmp_ovr',isz info=4010
call psb_realloc(isz,tmp_ovr_idx,info) call psb_errpush(info,name,a_err='psb_check_size')
if (info /= 0) then goto 9999
info=4010 end if
ch_err='psrealloc'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
End If
tmp_ovr_idx(counter_o)=proc tmp_ovr_idx(counter_o)=proc
tmp_ovr_idx(counter_o+1)=1 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 If (i_ovr < (n_ovr)) Then
n_elem = psb_sp_get_nnz_row(idx,a) n_elem = psb_sp_get_nnz_row(idx,a)
If((idxs+tot_elem+n_elem) > lworks) Then call psb_check_size((idxs+tot_elem+n_elem),works,info)
isz = max((3*lworks)/2,(idxs+tot_elem+n_elem)) if (info /= 0) then
if (debug) write(0,*) me,'Realloc works',isz info=4010
call psb_realloc(isz,works,info) call psb_errpush(info,name,a_err='psb_check_size')
if (info /= 0) then goto 9999
info=4010 end if
ch_err='psb_realloc'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
lworks = isz
End If
If((n_elem) > size(blk%ia2)) Then If((n_elem) > size(blk%ia2)) Then
isz = max((3*size(blk%ia2))/2,(n_elem)) 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 goto 9999
end if end if
Do i=1,iszr if (debug) write(0,*) 'ISZR :',iszr
idx=workr(i)
if (idx <1) then if (psb_is_large_desc(desc_a)) then
write(0,*) me,'Error in CDOVRBLD ',idx,i,iszr call psb_check_size(iszr,maskr,info)
!!$ write(0,*) me, ' WORKR :',workr(1:iszr) if (info /= 0) then
else If (desc_p%glob_to_loc(idx) < -np) Then info=4010
! call psb_errpush(info,name,a_err='psb_check_size')
! This is a new index. Assigning a local index as goto 9999
! we receive the guarantees that all indices for HALO(I) end if
! will be less than those for HALO(J) whenever I<J call psi_idx_cnv(iszr,workr,maskr,desc_p,info)
! iszs = count(maskr(1:iszr)<=0)
n_col=n_col+1 if (iszs > size(works)) call psb_realloc(iszs,works,info)
proc_id=-desc_p%glob_to_loc(idx)-np-1 j = 0
If(n_col > Size(desc_p%loc_to_glob)) Then do i=1,iszr
isz = 3*n_col/2 if (maskr(i) < 0) then
if (debug) write(0,*) me,'Realloc loc_to_glob' j=j+1
call psb_realloc(isz,desc_p%loc_to_glob,info) 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<J
!
proc_id = temp(i)
call psb_check_size((counter_t+3),t_halo_in,info,pad=-1)
if (info /= 0) then if (info /= 0) then
info=4010 info=4010
ch_err='psrealloc' call psb_errpush(info,name,a_err='psb_check_size')
call psb_errpush(info,name,a_err=ch_err)
goto 9999 goto 9999
end if end if
End If
desc_p%glob_to_loc(idx)=n_col t_halo_in(counter_t)=proc_id
desc_p%loc_to_glob(n_col)=idx t_halo_in(counter_t+1)=1
If ((counter_t+3) > Size(t_halo_in)) then t_halo_in(counter_t+2)=lidx
isz = max((3*Size(t_halo_in))/2,(counter_t+3+1000)) counter_t=counter_t+3
if (debug) write(0,*) me,'Realloc ovr_el',isz if (.false.) write(0,*) me,' CDOVRBLD: Added t_halo_in ',&
call psb_realloc(isz,t_halo_in,info) &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<J
!
n_col=n_col+1
proc_id=-desc_p%glob_to_loc(idx)-np-1
call psb_check_size(n_col,desc_p%loc_to_glob,info,pad=-1)
if (info /= 0) then if (info /= 0) then
info=4010 info=4010
ch_err='psrealloc' call psb_errpush(info,name,a_err='psb_check_size')
call psb_errpush(info,name,a_err=ch_err)
goto 9999 goto 9999
end if end if
end If
desc_p%glob_to_loc(idx)=n_col
t_halo_in(counter_t)=proc_id desc_p%loc_to_glob(n_col)=idx
t_halo_in(counter_t+1)=1
t_halo_in(counter_t+2)=n_col call psb_check_size((counter_t+3),t_halo_in,info,pad=-1)
counter_t=counter_t+3 if (info /= 0) then
if (debug) write(0,*) me,' CDOVRBLD: Added into t_halo_in from recv',& info=4010
&proc_id,n_col,idx call psb_errpush(info,name,a_err='psb_check_size')
else if (desc_p%glob_to_loc(idx) < 0) Then goto 9999
if (debug) write(0,*) me,'Wrong input to cdovrbld??',& end if
&idx,desc_p%glob_to_loc(idx)
End If t_halo_in(counter_t)=proc_id
End Do t_halo_in(counter_t+1)=1
t_halo_in(counter_t+2)=n_col
counter_t=counter_t+3
if (debug) write(0,*) me,' CDOVRBLD: Added into t_halo_in from recv',&
&proc_id,n_col,idx
else if (desc_p%glob_to_loc(idx) < 0) Then
if (debug) write(0,*) me,'Wrong input to cdovrbld??',&
&idx,desc_p%glob_to_loc(idx)
End If
End Do
desc_p%matrix_data(psb_n_col_)=n_col
end if
end if end if
t2 = mpi_wtime() t2 = mpi_wtime()
n_row=m+tot_recv !!$ desc_p%matrix_data(psb_n_row_)=desc_p%matrix_data(psb_n_col_)
desc_p%matrix_data(psb_n_row_)=n_row
desc_p%matrix_data(psb_n_col_)=n_col
! !
! Ok, now we have a temporary halo with all the info for the ! Ok, now we have a temporary halo with all the info for the
! next round. If we need to keep going, convert the halo format ! next round. If we need to keep going, convert the halo format
@ -493,18 +549,6 @@ Subroutine psb_zcdovrbld(n_ovr,desc_p,desc_a,a,&
If (i_ovr < (n_ovr)) Then If (i_ovr < (n_ovr)) Then
If(lwork < (counter_t/3+np*3)) Then
isz = max((3*lwork)/2,(counter_t/3+np*3))
if (debug) write(0,*) me,'Realloc work',isz
deallocate(work)
allocate(work(isz),stat=info)
if (info /= 0) then
call psb_errpush(4010,name,a_err='Allocate')
goto 9999
end if
lwork=size(work)
Endif
t_halo_in(counter_t)=-1 t_halo_in(counter_t)=-1
if (debug) write(0,*) me,'Checktmp_o_i 1',tmp_ovr_idx(1:10) if (debug) write(0,*) me,'Checktmp_o_i 1',tmp_ovr_idx(1:10)
@ -518,7 +562,7 @@ Subroutine psb_zcdovrbld(n_ovr,desc_p,desc_a,a,&
end if end if
if (debug) write(0,*) me,'Checktmp_o_i 2',tmp_ovr_idx(1:10) if (debug) write(0,*) me,'Checktmp_o_i 2',tmp_ovr_idx(1:10)
if (debug) write(0,*) me,'Done Crea_Halo' if (debug) write(0,*) me,'Done Crea_Halo'
call psb_transfer( t_halo_out,halo,info) call psb_transfer(t_halo_out,halo,info)
! !
! At this point we have built the halo necessary for I_OVR+1. ! At this point we have built the halo necessary for I_OVR+1.
! !
@ -534,14 +578,14 @@ Subroutine psb_zcdovrbld(n_ovr,desc_p,desc_a,a,&
desc_p%matrix_data(psb_m_)=psb_cd_get_global_rows(desc_a) desc_p%matrix_data(psb_m_)=psb_cd_get_global_rows(desc_a)
desc_p%matrix_data(psb_n_)=psb_cd_get_global_cols(desc_a) desc_p%matrix_data(psb_n_)=psb_cd_get_global_cols(desc_a)
tmp_halo(counter_h)=-1 tmp_halo(counter_h:)=-1
tmp_ovr_idx(counter_o)=-1 tmp_ovr_idx(counter_o:)=-1
! !
! At this point we have gathered all the indices in the halo at ! At this point we have gathered all the indices in the halo at
! N levels of overlap. Just call convert_comm. This is ! N levels of overlap. Just call cnv_dsc. This is
! the same routine as gets called inside SPASB. ! the same routine as gets called inside CDASB.
! !
if (debug) then if (debug) then
@ -549,45 +593,19 @@ Subroutine psb_zcdovrbld(n_ovr,desc_p,desc_a,a,&
call psb_barrier(ictxt) call psb_barrier(ictxt)
end if end if
!.... convert comunication stuctures.... !.... convert comunication stuctures....
! Note that we have to keep local_rows until the very end,
call psi_cnv_dsc(tmp_halo,tmp_ovr_idx,desc_p,info) ! because otherwise the halo build mechanism of cdasb
! will not work.
! Ok, register into MATRIX_DATA & free temporary work areas call psb_transfer(tmp_ovr_idx,desc_p%ovrlap_index,info)
desc_p%matrix_data(psb_dec_type_) = psb_desc_asb_ call psb_transfer(tmp_halo,desc_p%halo_index,info)
call psb_cdasb(desc_p,info)
allocate(desc_p%lprm(1),stat=info) desc_p%matrix_data(psb_n_row_)=desc_p%matrix_data(psb_n_col_)
if (info /= 0) then
call psb_errpush(4010,name,a_err='Allocate')
goto 9999
end if
desc_p%lprm(1) = 0
if (debug) then if (debug) then
write(0,*) me,'Done Convert_comm' write(0,*) me,'Done CDASB'
call psb_barrier(ictxt) call psb_barrier(ictxt)
end if end if
if (.false.) then
call psb_cdprt(70+me,desc_p,.false.)
end if
if (debug) write(0,*) me,'Done ConvertComm'
!!$ write(0,*) 't_halo_out ',allocated(t_halo_out)
!!$ Deallocate(works,workr,t_halo_in,work,&
!!$ & length_dl,dep_list,stat=info)
!!$ if (info /= 0) then
!!$ ch_err='Deallocate 1'
!!$ call psb_errpush(4013,name,a_err=ch_err,i_err=(/info,0,0,0,0/))
!!$ goto 9999
!!$ end if
!!$ Deallocate(tmp_ovr_idx,tmp_halo,&
!!$ & brvindx,rvsz,sdsz,bsdindx,temp,halo,stat=info)
!!$ if (info /= 0) then
!!$ ch_err='Deallocate 2'
!!$ call psb_errpush(4013,name,a_err=ch_err,i_err=(/info,0,0,0,0/))
!!$ goto 9999
!!$ end if
if (info == 0) call psb_sp_free(blk,info) if (info == 0) call psb_sp_free(blk,info)
if (info /= 0) then if (info /= 0) then
ch_err='sp_free' ch_err='sp_free'

@ -88,10 +88,10 @@ subroutine psb_zcsrp(trans,iperm,a, desc_a, info)
time(1) = mpi_wtime() 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) dectype = psb_cd_get_dectype(desc_a)
n_row = psb_cd_get_local_rows(desc_a) n_row = psb_cd_get_local_rows(desc_a)
n_col = psb_cd_get_local_cols(desc_a) n_col = psb_cd_get_local_cols(desc_a)
if(psb_get_errstatus() /= 0) return if(psb_get_errstatus() /= 0) return
info=0 info=0

@ -57,7 +57,7 @@ subroutine psb_zgelp(trans,iperm,x,desc_a,info)
! local variables ! local variables
integer :: ictxt,np,me,nrow,ncol integer :: ictxt,np,me,nrow,ncol
complex(kind(1.d0)),pointer :: dtemp(:) complex(kind(1.d0)),pointer :: dtemp(:)
integer :: int_err(5), i1sz, i2sz, dectype, i, err_act integer :: int_err(5), i1sz, i2sz, i, err_act
real(kind(1.d0)),parameter :: one=1 real(kind(1.d0)),parameter :: one=1
logical, parameter :: debug=.false. logical, parameter :: debug=.false.
@ -86,8 +86,7 @@ subroutine psb_zgelp(trans,iperm,x,desc_a,info)
info=0 info=0
call psb_erractionsave(err_act) 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)
nrow = psb_cd_get_local_rows(desc_a) nrow = psb_cd_get_local_rows(desc_a)
ncol = psb_cd_get_local_cols(desc_a) ncol = psb_cd_get_local_cols(desc_a)
i1sz = size(x,dim=1) i1sz = size(x,dim=1)
@ -96,7 +95,7 @@ subroutine psb_zgelp(trans,iperm,x,desc_a,info)
call psb_info(ictxt, me, np) call psb_info(ictxt, me, np)
if (debug) write(*,*) 'asb start: ',np,me,& if (debug) write(*,*) 'asb start: ',np,me,&
&psb_cd_get_dectype(desc_a) & psb_cd_get_dectype(desc_a)
! ....verify blacs grid correctness.. ! ....verify blacs grid correctness..
if (np == -1) then if (np == -1) then
info = 2010 info = 2010
@ -200,7 +199,7 @@ subroutine psb_zgelpv(trans,iperm,x,desc_a,info)
! local variables ! local variables
integer :: ictxt,np,me integer :: ictxt,np,me
integer :: int_err(5), i1sz,nrow,ncol,dectype, i, err_act integer :: int_err(5), i1sz,nrow,ncol, i, err_act
complex(kind(1.d0)),pointer :: dtemp(:) complex(kind(1.d0)),pointer :: dtemp(:)
real(kind(1.d0)),parameter :: one=1 real(kind(1.d0)),parameter :: one=1
logical, parameter :: debug=.false. logical, parameter :: debug=.false.
@ -232,10 +231,9 @@ subroutine psb_zgelpv(trans,iperm,x,desc_a,info)
i1sz = size(x) i1sz = size(x)
ictxt=psb_cd_get_context(desc_a) ictxt = psb_cd_get_context(desc_a)
dectype=psb_cd_get_dectype(desc_a) nrow = psb_cd_get_local_rows(desc_a)
nrow=psb_cd_get_local_rows(desc_a) ncol = psb_cd_get_local_cols(desc_a)
ncol=psb_cd_get_local_cols(desc_a)
call psb_info(ictxt, me, np) call psb_info(ictxt, me, np)

@ -36,7 +36,7 @@
! val to be inserted. ! val to be inserted.
! irw - integer(:) Row indices of rows of val (global numbering) ! irw - integer(:) Row indices of rows of val (global numbering)
! val - complex, dimension(:). The source dense submatrix. ! 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(<psb_desc_type>). The communication descriptor. ! desc_a - type(<psb_desc_type>). The communication descriptor.
! info - integer. Eventually returns an error code ! info - integer. Eventually returns an error code
subroutine psb_zinsvi(m, irw, val, x, desc_a, info, dupl) 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_const_mod
use psb_error_mod use psb_error_mod
use psb_penv_mod use psb_penv_mod
use psi_mod
implicit none implicit none
! m rows number of submatrix belonging to val to be inserted ! 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 integer, optional, intent(in) :: dupl
!locals..... !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) & loc_rows,loc_cols,mglob,err_act, int_err(5)
integer :: np,me,dupl_ integer :: np, me, dupl_
character(len=20) :: name integer, allocatable :: irl(:)
character(len=20) :: name
if(psb_get_errstatus() /= 0) return if(psb_get_errstatus() /= 0) return
info=0 info=0
call psb_erractionsave(err_act) call psb_erractionsave(err_act)
name = 'psb_dinsvv' name = 'psb_zinsvi'
if (.not.allocated(desc_a%glob_to_loc)) then !!$ if (.not.allocated(desc_a%glob_to_loc)) then
info=3110 !!$ info=3110
call psb_errpush(info,name) !!$ call psb_errpush(info,name)
return !!$ return
end if !!$ end if
if ((.not.allocated(desc_a%matrix_data))) then if ((.not.allocated(desc_a%matrix_data))) then
int_err(1)=3110 int_err(1)=3110
call psb_errpush(info,name) call psb_errpush(info,name)
@ -112,31 +114,36 @@ subroutine psb_zinsvi(m, irw, val, x, desc_a, info, dupl)
goto 9999 goto 9999
endif endif
loc_rows=psb_cd_get_local_rows(desc_a) if (m==0) return
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) 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 if (present(dupl)) then
dupl_ = dupl dupl_ = dupl
else else
dupl_ = psb_dupl_ovwrt_ dupl_ = psb_dupl_ovwrt_
endif endif
call psi_idx_cnv(m,irw,irl,desc_a,info,owned=.true.)
select case(dupl_) select case(dupl_)
case(psb_dupl_ovwrt_) case(psb_dupl_ovwrt_)
do i = 1, m do i = 1, m
!loop over all val's rows !loop over all val's rows
! row actual block row ! row actual block row
glob_row=irw(i) if (irl(i) > 0) then
if ((glob_row>0).and.(glob_row <= mglob)) then ! this row belongs to me
! copy i-th row of block val in x
loc_row=desc_a%glob_to_loc(glob_row) x(irl(i)) = val(i)
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
end if end if
enddo enddo
@ -145,16 +152,10 @@ subroutine psb_zinsvi(m, irw, val, x, desc_a, info, dupl)
do i = 1, m do i = 1, m
!loop over all val's rows !loop over all val's rows
! row actual block row if (irl(i) > 0) then
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 ! this row belongs to me
! copy i-th row of block val in x ! copy i-th row of block val in x
x(loc_row) = x(loc_row) + val(i) x(irl(i)) = x(irl(i)) + val(i)
end if
end if end if
enddo enddo
@ -163,6 +164,7 @@ subroutine psb_zinsvi(m, irw, val, x, desc_a, info, dupl)
call psb_errpush(info,name) call psb_errpush(info,name)
goto 9999 goto 9999
end select end select
deallocate(irl)
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
return return
@ -217,7 +219,7 @@ end subroutine psb_zinsvi
! val to be inserted. ! val to be inserted.
! irw - integer(:) Row indices of rows of val (global numbering) ! irw - integer(:) Row indices of rows of val (global numbering)
! val - complex, dimension(:,:). The source dense submatrix. ! 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(<psb_desc_type>). The communication descriptor. ! desc_a - type(<psb_desc_type>). The communication descriptor.
! info - integer. Eventually returns an error code ! info - integer. Eventually returns an error code
subroutine psb_zinsi(m,irw, val, x, desc_a, info, dupl) 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_const_mod
use psb_error_mod use psb_error_mod
use psb_penv_mod use psb_penv_mod
use psi_mod
implicit none implicit none
! m rows number of submatrix belonging to val to be inserted ! 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,& integer :: ictxt,i,loc_row,glob_row,j,n,&
& loc_rows,loc_cols,mglob,err_act, int_err(5) & 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 character(len=20) :: name
if(psb_get_errstatus() /= 0) return if(psb_get_errstatus() /= 0) return
info=0 info=0
call psb_erractionsave(err_act) call psb_erractionsave(err_act)
name = 'psb_dinsvv' name = 'psb_zinsi'
if (.not.allocated(desc_a%glob_to_loc)) then !!$ if (.not.allocated(desc_a%glob_to_loc)) then
info=3110 !!$ info=3110
call psb_errpush(info,name) !!$ call psb_errpush(info,name)
return !!$ return
end if !!$ end if
if ((.not.allocated(desc_a%matrix_data))) then if ((.not.allocated(desc_a%matrix_data))) then
int_err(1)=3110 int_err(1)=3110
call psb_errpush(info,name) call psb_errpush(info,name)
@ -274,7 +278,7 @@ subroutine psb_zinsi(m,irw, val, x, desc_a, info, dupl)
endif endif
!... check parameters.... !... check parameters....
if (m.lt.0) then if (m < 0) then
info = 10 info = 10
int_err(1) = 1 int_err(1) = 1
int_err(2) = m 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) call psb_errpush(info,name,int_err)
goto 9999 goto 9999
endif endif
if (m==0) return
loc_rows=psb_cd_get_local_rows(desc_a) loc_rows = psb_cd_get_local_rows(desc_a)
loc_cols=psb_cd_get_local_cols(desc_a) loc_cols = psb_cd_get_local_cols(desc_a)
mglob = psb_cd_get_global_rows(desc_a) mglob = psb_cd_get_global_rows(desc_a)
n = min(size(val,2),size(x,2)) 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_ dupl_ = psb_dupl_ovwrt_
endif 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_) select case(dupl_)
case(psb_dupl_ovwrt_) case(psb_dupl_ovwrt_)
do i = 1, m do i = 1, m
!loop over all val's rows !loop over all val's rows
! row actual block row ! row actual block row
glob_row=irw(i) loc_row = irl(i)
if ((glob_row>0).and.(glob_row <= mglob)) then if (loc_row > 0) then
! this row belongs to me
loc_row=desc_a%glob_to_loc(glob_row) ! copy i-th row of block val in x
if (loc_row.ge.1) then do j=1,n
! this row belongs to me x(loc_row,j) = val(i,j)
! copy i-th row of block val in x end do
do j=1,n
x(loc_row,j) = val(i,j)
end do
end if
end if end if
enddo enddo
@ -331,17 +341,13 @@ subroutine psb_zinsi(m,irw, val, x, desc_a, info, dupl)
!loop over all val's rows !loop over all val's rows
! row actual block row ! row actual block row
glob_row=irw(i) loc_row = irl(i)
if ((glob_row>0).and.(glob_row <= mglob)) then if (loc_row > 0) then
! this row belongs to me
loc_row=desc_a%glob_to_loc(glob_row) ! copy i-th row of block val in x
if (loc_row.ge.1) then do j=1,n
! this row belongs to me x(loc_row,j) = x(loc_row,j) + val(i,j)
! copy i-th row of block val in x end do
do j=1,n
x(loc_row,j) = x(loc_row,j) + val(i,j)
end do
end if
end if end if
enddo enddo
@ -350,6 +356,7 @@ subroutine psb_zinsi(m,irw, val, x, desc_a, info, dupl)
call psb_errpush(info,name) call psb_errpush(info,name)
goto 9999 goto 9999
end select end select
deallocate(irl)
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
return return

@ -56,7 +56,7 @@ subroutine psb_zspalloc(a, desc_a, info, nnz)
integer, optional, intent(in) :: nnz integer, optional, intent(in) :: nnz
!locals !locals
integer :: ictxt, dectype integer :: ictxt
integer :: np,me,loc_row,& integer :: np,me,loc_row,&
& length_ia1,length_ia2, err_act,m,n & length_ia1,length_ia2, err_act,m,n
integer :: int_err(5) integer :: int_err(5)
@ -69,7 +69,6 @@ subroutine psb_zspalloc(a, desc_a, info, nnz)
name = 'psb_zspalloc' name = 'psb_zspalloc'
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) call psb_info(ictxt, me, np)
! ....verify blacs grid correctness.. ! ....verify blacs grid correctness..
if (np == -1) then if (np == -1) then
@ -127,7 +126,7 @@ subroutine psb_zspalloc(a, desc_a, info, nnz)
a%infoa(psb_state_) = psb_spmat_bld_ a%infoa(psb_state_) = psb_spmat_bld_
if (debug) write(0,*) 'spall: ', & if (debug) write(0,*) 'spall: ', &
&psb_cd_get_dectype(desc_a),psb_desc_bld_ & psb_cd_get_dectype(desc_a),psb_desc_bld_
return return

@ -65,7 +65,7 @@ subroutine psb_zspasb(a,desc_a, info, afmt, upd, dupl)
integer :: int_err(5) integer :: int_err(5)
type(psb_zspmat_type) :: atemp type(psb_zspmat_type) :: atemp
integer :: np,me,n_col,iout, err_act integer :: np,me,n_col,iout, err_act
integer :: dscstate, spstate integer :: spstate
integer :: upd_, dupl_ integer :: upd_, dupl_
integer :: ictxt,n_row integer :: ictxt,n_row
logical, parameter :: debug=.false., debugwrt=.false. logical, parameter :: debug=.false., debugwrt=.false.
@ -76,8 +76,7 @@ subroutine psb_zspasb(a,desc_a, info, afmt, upd, dupl)
name = 'psb_spasb' name = 'psb_spasb'
call psb_erractionsave(err_act) 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_row = psb_cd_get_local_rows(desc_a)
n_col = psb_cd_get_local_cols(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 goto 9999
endif endif
if (.not.psb_is_asb_dec(dscstate)) then if (.not.psb_is_asb_desc(desc_a)) then
info = 600 info = 600
int_err(1) = dscstate int_err(1) = psb_cd_get_dectype(desc_a)
call psb_errpush(info,name) call psb_errpush(info,name)
goto 9999 goto 9999
endif endif
@ -100,7 +99,7 @@ subroutine psb_zspasb(a,desc_a, info, afmt, upd, dupl)
!check on errors encountered in psdspins !check on errors encountered in psdspins
spstate = a%infoa(psb_state_) spstate = a%infoa(psb_state_)
if (spstate == psb_spmat_bld_) then if (spstate == psb_spmat_bld_) then
! !
! First case: we come from a fresh build. ! First case: we come from a fresh build.

@ -122,8 +122,7 @@ subroutine psb_zspcnv(a,b,desc_a,info)
time(1) = mpi_wtime() 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) dectype = psb_cd_get_dectype(desc_a)
n_row = psb_cd_get_local_rows(desc_a) n_row = psb_cd_get_local_rows(desc_a)
n_col = psb_cd_get_local_cols(desc_a) n_col = psb_cd_get_local_cols(desc_a)

@ -61,11 +61,11 @@ subroutine psb_zspfree(a, desc_a,info)
call psb_erractionsave(err_act) call psb_erractionsave(err_act)
if (.not.allocated(desc_a%matrix_data)) then if (.not.allocated(desc_a%matrix_data)) then
info=295 info = 295
call psb_errpush(info,name) call psb_errpush(info,name)
return return
else else
ictxt=psb_cd_get_context(desc_a) ictxt = psb_cd_get_context(desc_a)
end if end if
!...deallocate a.... !...deallocate a....

@ -100,7 +100,8 @@ Subroutine psb_zsphalo(a,desc_a,blk,info,rwcnv,clcnv,outfmt)
outfmt_ = 'CSR' outfmt_ = 'CSR'
endif endif
ictxt=psb_cd_get_context(desc_a) ictxt = psb_cd_get_context(desc_a)
Call psb_info(ictxt, me, np) Call psb_info(ictxt, me, np)
t1 = mpi_wtime() t1 = mpi_wtime()

@ -52,6 +52,7 @@ subroutine psb_zspins(nz,ia,ja,val,a,desc_a,info,rebuild)
use psb_const_mod use psb_const_mod
use psb_error_mod use psb_error_mod
use psb_penv_mod use psb_penv_mod
use psb_tools_mod
implicit none implicit none
!....parameters... !....parameters...
@ -63,22 +64,34 @@ subroutine psb_zspins(nz,ia,ja,val,a,desc_a,info,rebuild)
logical, intent(in), optional :: rebuild logical, intent(in), optional :: rebuild
!locals..... !locals.....
integer :: nrow,err_act, dectype,mglob,ncol, spstate integer :: nrow, err_act,mglob,ncol, spstate
integer :: ictxt,np, me integer :: ictxt,np,me
logical, parameter :: debug=.false. logical, parameter :: debug=.false.
integer, parameter :: relocsz=200 integer, parameter :: relocsz=200
logical :: rebuild_ logical :: rebuild_
integer, allocatable :: ila(:),jla(:)
interface psb_cdins !!$ 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 !!$ use psb_descriptor_type
implicit none !!$ implicit none
type(psb_desc_type), intent(inout) :: desc_a !!$ type(psb_desc_type), intent(inout) :: desc_a
integer, intent(in) :: nz,ia(:),ja(:) !!$ integer, intent(in) :: nz,ia(:),ja(:)
integer, intent(out) :: info !!$ integer, intent(out) :: info
end subroutine psb_cdins !!$ integer, optional, intent(out) :: ila(:), jla(:)
end interface !!$ 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 character(len=20) :: name, ch_err
info = 0 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) 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) call psb_info(ictxt, me, np)
@ -98,7 +110,7 @@ subroutine psb_zspins(nz,ia,ja,val,a,desc_a,info,rebuild)
goto 9999 goto 9999
endif endif
if (nz <= 0) then if (nz < 0) then
info = 1111 info = 1111
call psb_errpush(info,name) call psb_errpush(info,name)
goto 9999 goto 9999
@ -119,6 +131,7 @@ subroutine psb_zspins(nz,ia,ja,val,a,desc_a,info,rebuild)
call psb_errpush(info,name) call psb_errpush(info,name)
goto 9999 goto 9999
end if end if
if (nz==0) return
if (present(rebuild)) then if (present(rebuild)) then
rebuild_ = rebuild rebuild_ = rebuild
@ -128,39 +141,102 @@ subroutine psb_zspins(nz,ia,ja,val,a,desc_a,info,rebuild)
spstate = a%infoa(psb_state_) spstate = a%infoa(psb_state_)
if (psb_is_bld_desc(desc_a)) then if (psb_is_bld_desc(desc_a)) then
call psb_cdins(nz,ia,ja,desc_a,info) if (psb_is_large_desc(desc_a)) then
if (info /= 0) then
info=4010 allocate(ila(nz),jla(nz),stat=info)
ch_err='psb_cdins' if (info /= 0) then
call psb_errpush(info,name,a_err=ch_err) ch_err='allocate'
goto 9999 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 end if
nrow = psb_cd_get_local_rows(desc_a)
ncol = psb_cd_get_local_cols(desc_a)
if (spstate == psb_spmat_bld_) then else if (psb_is_asb_desc(desc_a)) then
call psb_coins(nz,ia,ja,val,a,1,nrow,1,ncol,info,gtl=desc_a%glob_to_loc)
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 if (info /= 0) then
info=4010 info=4010
ch_err='psb_coins' ch_err='psb_coins'
call psb_errpush(info,name,a_err=ch_err) call psb_errpush(info,name,a_err=ch_err)
goto 9999 goto 9999
end if end if
else else
info = 1123 nrow = psb_cd_get_local_rows(desc_a)
call psb_errpush(info,name) ncol = psb_cd_get_local_cols(desc_a)
goto 9999 call psb_coins(nz,ia,ja,val,a,1,nrow,1,ncol,&
end if & info,gtl=desc_a%glob_to_loc,rebuild=rebuild_)
else if (psb_is_asb_desc(desc_a)) then if (info /= 0) then
nrow = psb_cd_get_local_rows(desc_a) info=4010
ncol = psb_cd_get_local_cols(desc_a) ch_err='psb_coins'
call psb_coins(nz,ia,ja,val,a,1,nrow,1,ncol,& call psb_errpush(info,name,a_err=ch_err)
& info,gtl=desc_a%glob_to_loc,rebuild=rebuild_) goto 9999
if (info /= 0) then end if
info=4010
ch_err='psb_coins'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if end if
else else
info = 1122 info = 1122

Loading…
Cancel
Save