diff --git a/base/internals/psi_extrct_dl.F90 b/base/internals/psi_extrct_dl.F90 index ef004dfc..f5fe3c58 100644 --- a/base/internals/psi_extrct_dl.F90 +++ b/base/internals/psi_extrct_dl.F90 @@ -126,6 +126,7 @@ subroutine psi_i_extract_dep_list(ictxt,is_bld,is_upd,desc_str,dep_list,& use psb_const_mod use psb_error_mod use psb_desc_mod + use psb_sort_mod implicit none #ifdef MPI_H include 'mpif.h' @@ -145,7 +146,8 @@ subroutine psi_i_extract_dep_list(ictxt,is_bld,is_upd,desc_str,dep_list,& integer(psb_ipk_) :: i,pointer_dep_list,proc,j,err_act integer(psb_ipk_) :: err integer(psb_ipk_) :: debug_level, debug_unit - integer(psb_mpk_) :: iictxt, icomm, me, np, dl_mpi, minfo + integer(psb_mpk_) :: iictxt, icomm, me, np, minfo + logical, parameter :: dist_symm_list=.false., print_dl=.false. character name*20 name='psi_extrct_dl' @@ -156,13 +158,14 @@ subroutine psi_i_extract_dep_list(ictxt,is_bld,is_upd,desc_str,dep_list,& info = psb_success_ call psb_info(iictxt,me,np) - allocate(itmp(np+1),length_dl(0:np),stat=info) + allocate(itmp(2*np+1),length_dl(0:np),stat=info) if (info /= psb_success_) then info=psb_err_alloc_dealloc_ goto 9999 end if do i=0,np length_dl(i) = 0 + itmp(i+1) = -1 enddo i=1 if (debug_level >= psb_debug_inner_)& @@ -249,32 +252,116 @@ subroutine psi_i_extract_dep_list(ictxt,is_bld,is_upd,desc_str,dep_list,& endif length_dl(me)=pointer_dep_list-1 - dl_lda = max(length_dl(me),1) - call psb_max(iictxt, dl_lda) - ! - ! This doubling of DL_LDA is not 100% safe, - ! but should work most of the time. - ! Will need to be improved later, perhaps move - ! from a 2D allocation (ellpack style) to - ! a 1D allocation (csr like). - ! - dl_lda = min(2*dl_lda,np+1) - allocate(dep_list(dl_lda,0:np),stat=info) - if (info /= psb_success_) then - call psb_errpush(psb_err_from_subroutine_,name,a_err='Allocate') - goto 9999 - end if + if (dist_symm_list) then + call psb_realloc(length_dl(me),itmp,info) + call psi_symm_dep_list(itmp,ictxt,info) + dl_lda = max(size(itmp),1) + call psb_max(iictxt, dl_lda) + call psb_realloc(dl_lda,itmp,info) + ! dl_lda = min(np,2*dl_lda) + allocate(dep_list(dl_lda,0:np),stat=info) + if (info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,a_err='Allocate') + goto 9999 + end if - call psb_sum(iictxt,length_dl(0:np)) - icomm = psb_get_mpicomm(iictxt) - call mpi_allgather(itmp,dl_lda,psb_mpi_ipk_,& - & dep_list,dl_lda,psb_mpi_ipk_,icomm,minfo) - info = minfo - if (info == 0) deallocate(itmp,stat=info) - if (info /= psb_success_) then - info=psb_err_alloc_dealloc_ - goto 9999 + call psb_sum(iictxt,length_dl(0:np)) + icomm = psb_get_mpicomm(iictxt) + call mpi_allgather(itmp,dl_lda,psb_mpi_ipk_,& + & dep_list,dl_lda,psb_mpi_ipk_,icomm,minfo) + info = minfo + if (info == 0) deallocate(itmp,stat=info) + if (info /= psb_success_) then + info=psb_err_alloc_dealloc_ + goto 9999 + endif + else + block + integer(psb_ipk_), allocatable :: list1(:,:), ldl2(:), list2(:,:) + integer(psb_ipk_) :: i,j,ip,dlsym, ldu, mdl, l1, l2 + + dl_lda = max(length_dl(me),1) + call psb_max(iictxt, dl_lda) + allocate(dep_list(dl_lda,0:np),stat=info) + if (info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,a_err='Allocate') + goto 9999 + end if + call psb_sum(iictxt,length_dl(0:np)) + icomm = psb_get_mpicomm(iictxt) + call mpi_allgather(itmp,dl_lda,psb_mpi_ipk_,& + & dep_list,dl_lda,psb_mpi_ipk_,icomm,minfo) + info = minfo + if (info /= psb_success_) then + info=psb_err_alloc_dealloc_ + goto 9999 + endif + allocate(ldl2(0:np),stat=info) + ldl2 = 0 + do j=0, np-1 + do i=1,length_dl(j) + ip = dep_list(i,j) + ldl2(ip) = ldl2(ip) + 1 + end do + end do + dlsym = maxval(ldl2) + allocate(list2(dlsym,0:np),stat=info) + call move_alloc(dep_list,list1) + if (info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,a_err='Allocate') + goto 9999 + end if + ldl2 = 0 + do j=0, np-1 + do i=1,length_dl(j) + ip = list1(i,j) + ldl2(ip) = ldl2(ip) + 1 + list2(ldl2(ip),ip) = j + end do + end do + mdl = 0 + do j = 0, np-1 + l1 = length_dl(j) + l2 = ldl2(j) + itmp(1:l1) = list1(1:l1,j) + itmp(l1+1:l1+l2) = list2(1:l2,j) + ldu = l1 + l2 + !if (me == 0) write(0,*) 'Iter ',j,':',l1,l2,':',itmp(1:l1),':',itmp(l1+1:l1+l2) + call psb_msort_unique(itmp(1:l1+l2),ldu) + mdl = max(mdl, ldu) + !if (me == 0) write(0,*) 'Iter ',j,':',ldu,':',itmp(1:ldu) + end do + dl_lda = mdl + allocate(dep_list(dl_lda,0:np),stat=info) + dep_list = -1 + if (info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,a_err='Allocate') + goto 9999 + end if + do j = 0, np-1 + l1 = length_dl(j) + l2 = ldl2(j) + itmp(1:l1) = list1(1:l1,j) + itmp(l1+1:l1+l2) = list2(1:l2,j) + ldu = l1 + l2 + call psb_msort_unique(itmp(1:l1+l2),ldu) + length_dl(j) = ldu + dep_list(1:ldu,j) = itmp(1:ldu) + end do + + end block endif + if (print_dl) then + if (me == 0) then + write(0,*) ' Dep_list ' + do i=0,np-1 + j = length_dl(i) + write(0,*) 'Proc ',i,':',dep_list(1:j,i) + end do + flush(0) + end if + call psb_barrier(ictxt) + end if call psb_erractionrestore(err_act) return diff --git a/base/internals/psi_symm_dep_list.F90 b/base/internals/psi_symm_dep_list.F90 index c25f61e5..c3c4db84 100644 --- a/base/internals/psi_symm_dep_list.F90 +++ b/base/internals/psi_symm_dep_list.F90 @@ -37,13 +37,13 @@ ! ! Arguments: ! -subroutine psi_symm_dep_list(rvsz,adj,idxmap,info,flag) +subroutine psi_symm_dep_list_inrv(rvsz,adj,ictxt,info) use psb_serial_mod use psb_const_mod use psb_error_mod use psb_penv_mod use psb_realloc_mod - use psb_indx_map_mod, psb_protect_name => psi_symm_dep_list + use psb_indx_map_mod, psb_protect_name => psi_symm_dep_list_inrv #ifdef MPI_MOD use mpi #endif @@ -54,30 +54,23 @@ subroutine psi_symm_dep_list(rvsz,adj,idxmap,info,flag) #endif integer(psb_mpk_), intent(inout) :: rvsz(0:) integer(psb_ipk_), allocatable, intent(inout) :: adj(:) - class(psb_indx_map), intent(in) :: idxmap + integer(psb_ipk_), intent(in) :: ictxt integer(psb_ipk_), intent(out) :: info - integer(psb_ipk_), intent(in), optional :: flag ! integer(psb_ipk_), allocatable :: ladj(:) - integer(psb_mpk_), allocatable :: sdsz(:) - integer(psb_mpk_) :: icomm, minfo, iictxt + integer(psb_mpk_) :: icomm, minfo integer(psb_ipk_) :: i,n_row,n_col,err_act,hsize,ip,& & last_ih, last_j, nidx, nrecv, nadj, flag_ integer(psb_lpk_) :: mglob, ih - integer(psb_ipk_) :: ictxt,np,me + integer(psb_ipk_) :: np,me character(len=20) :: name info = psb_success_ name = 'psi_symm_dep_list' call psb_erractionsave(err_act) - ictxt = idxmap%get_ctxt() - icomm = idxmap%get_mpic() - mglob = idxmap%get_gr() - n_row = idxmap%get_lr() - n_col = idxmap%get_lc() - iictxt = ictxt + icomm = psb_get_mpicomm(ictxt) call psb_info(ictxt, me, np) @@ -87,51 +80,15 @@ subroutine psi_symm_dep_list(rvsz,adj,idxmap,info,flag) goto 9999 endif - if (.not.(idxmap%is_valid())) then - call psb_errpush(psb_err_from_subroutine_,name,a_err='invalid idxmap') - goto 9999 - end if - - nadj = size(adj) + + ! + ! Am getting this from the outside, possibly from the wrapper call + ! if (size(rvsz) 0) nrecv = nrecv + 1 @@ -145,13 +102,87 @@ subroutine psi_symm_dep_list(rvsz,adj,idxmap,info,flag) do ip=0, np-1 if (rvsz(ip)>0) then nadj = nadj + 1 - ladj(nadj+1:nadj+nrecv) = ip + ladj(nadj) = ip end if end do call psb_msort_unique(ladj,nadj) call psb_realloc(nadj,adj,info) - adj(1:nadj) = ladj(1:nadj) + if (info == 0) adj(1:nadj) = ladj(1:nadj) + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(ictxt,err_act) + + return + +end subroutine psi_symm_dep_list_inrv + +subroutine psi_symm_dep_list_norv(adj,ictxt,info) + use psb_serial_mod + use psb_const_mod + use psb_error_mod + use psb_penv_mod + use psb_realloc_mod + use psb_indx_map_mod, psb_protect_name => psi_symm_dep_list_norv +#ifdef MPI_MOD + use mpi +#endif + + implicit none +#ifdef MPI_H + include 'mpif.h' +#endif + integer(psb_ipk_), allocatable, intent(inout) :: adj(:) + integer(psb_ipk_), intent(in) :: ictxt + integer(psb_ipk_), intent(out) :: info + ! + integer(psb_mpk_), allocatable :: rvsz(:), sdsz(:) + + integer(psb_mpk_) :: icomm, minfo + integer(psb_ipk_) :: i,n_row,n_col,err_act,hsize,ip,& + & last_ih, last_j, nidx, nrecv, nadj + integer(psb_ipk_) :: mglob, ih + integer(psb_ipk_) :: np,me + character(len=20) :: name + + info = psb_success_ + name = 'psi_symm_dep_list' + call psb_erractionsave(err_act) + + icomm = psb_get_mpicomm(ictxt) + + call psb_info(ictxt, me, np) + + if (np == -1) then + info = psb_err_context_error_ + call psb_errpush(info,name) + goto 9999 + endif + + nadj = size(adj) + + ! write(0,*) me,name,' Going through ',nidx,nadj + + Allocate(sdsz(0:np-1), rvsz(0:np-1), stat=info) + ! + ! First, send sizes according to adjcncy list + ! + sdsz = 0 + do i=1, nadj + sdsz(adj(i)) = 1 + end do + !write(0,*)me,' Check on sizes into a2a:',adj(:),nadj,':',sdsz(:) + + call mpi_alltoall(sdsz,1,psb_mpi_mpk_,& + & rvsz,1,psb_mpi_mpk_,icomm,minfo) + if (minfo == 0) call psi_symm_dep_list(rvsz,adj,ictxt,info) + if ((minfo /=0).or.(info /= 0)) then + call psb_errpush(psb_err_from_subroutine_,name,a_err='inner call symm_dep') + goto 9999 + end if + call psb_erractionrestore(err_act) return @@ -159,4 +190,4 @@ subroutine psi_symm_dep_list(rvsz,adj,idxmap,info,flag) return -end subroutine psi_symm_dep_list +end subroutine psi_symm_dep_list_norv diff --git a/base/modules/desc/psb_indx_map_mod.f90 b/base/modules/desc/psb_indx_map_mod.f90 index eed22819..378b38af 100644 --- a/base/modules/desc/psb_indx_map_mod.f90 +++ b/base/modules/desc/psb_indx_map_mod.f90 @@ -313,17 +313,23 @@ module psb_indx_map_mod integer, parameter :: psi_symm_flag_norv_ = 0 integer, parameter :: psi_symm_flag_inrv_ = 1 - interface - subroutine psi_symm_dep_list(rvsz,adj,idxmap,info,flag) + interface psi_symm_dep_list + subroutine psi_symm_dep_list_inrv(rvsz,adj,ictxt,info) import :: psb_indx_map, psb_ipk_, psb_lpk_, psb_mpk_ implicit none - integer(psb_mpk_), intent(inout) :: rvsz(:) - integer(psb_ipk_), intent(in) :: adj(:) - class(psb_indx_map), intent(in) :: idxmap - integer(psb_ipk_), intent(out) :: info - integer(psb_ipk_), intent(in), optional :: flag - end subroutine psi_symm_dep_list - end interface + integer(psb_mpk_), intent(inout) :: rvsz(0:) + integer(psb_ipk_), allocatable, intent(inout) :: adj(:) + integer(psb_ipk_), intent(in) :: ictxt + integer(psb_ipk_), intent(out) :: info + end subroutine psi_symm_dep_list_inrv + subroutine psi_symm_dep_list_norv(adj,ictxt,info) + import :: psb_indx_map, psb_ipk_, psb_lpk_, psb_mpk_ + implicit none + integer(psb_ipk_), allocatable, intent(inout) :: adj(:) + integer(psb_ipk_), intent(in) :: ictxt + integer(psb_ipk_), intent(out) :: info + end subroutine psi_symm_dep_list_norv + end interface psi_symm_dep_list contains