diff --git a/amgprec/amg_c_base_aggregator_mod.f90 b/amgprec/amg_c_base_aggregator_mod.f90 index 26351dfe..376989b5 100644 --- a/amgprec/amg_c_base_aggregator_mod.f90 +++ b/amgprec/amg_c_base_aggregator_mod.f90 @@ -89,7 +89,7 @@ module amg_c_base_aggregator_mod procedure, pass(ag) :: bld_tprol => amg_c_base_aggregator_build_tprol procedure, pass(ag) :: mat_bld => amg_c_base_aggregator_mat_bld procedure, pass(ag) :: mat_asb => amg_c_base_aggregator_mat_asb - procedure, pass(ag) :: bld_map => amg_c_base_aggregator_bld_map + procedure, pass(ag) :: bld_linmap => amg_c_base_aggregator_bld_linmap procedure, pass(ag) :: update_next => amg_c_base_aggregator_update_next procedure, pass(ag) :: clone => amg_c_base_aggregator_clone procedure, pass(ag) :: free => amg_c_base_aggregator_free @@ -458,7 +458,7 @@ contains end subroutine amg_c_base_aggregator_mat_asb ! - !> Function bld_map + !> Function bld_linmap !! \memberof amg_c_base_aggregator_type !! \brief Build linear map between hierarchy levels !! @@ -473,7 +473,7 @@ contains !! \param map The output map !! \param info Return code !! - subroutine amg_c_base_aggregator_bld_map(ag,desc_a,desc_ac,ilaggr,nlaggr,& + subroutine amg_c_base_aggregator_bld_linmap(ag,desc_a,desc_ac,ilaggr,nlaggr,& & op_restr,op_prol,map,info) use psb_base_mod implicit none @@ -484,7 +484,7 @@ contains type(psb_clinmap_type), intent(out) :: map integer(psb_ipk_), intent(out) :: info integer(psb_ipk_) :: err_act - character(len=20) :: name='c_base_aggregator_bld_map' + character(len=20) :: name='c_base_aggregator_bld_linmap' info = psb_success_ call psb_erractionsave(err_act) @@ -508,7 +508,6 @@ contains 9999 call psb_error_handler(err_act) return - end subroutine amg_c_base_aggregator_bld_map - + end subroutine amg_c_base_aggregator_bld_linmap end module amg_c_base_aggregator_mod diff --git a/amgprec/amg_c_onelev_mod.f90 b/amgprec/amg_c_onelev_mod.f90 index f448a66c..ab044e74 100644 --- a/amgprec/amg_c_onelev_mod.f90 +++ b/amgprec/amg_c_onelev_mod.f90 @@ -181,7 +181,8 @@ module amg_c_onelev_mod integer(psb_ipk_) :: idest integer(psb_ipk_), allocatable :: isrc(:), nrsrc(:), naggr(:) contains - procedure, pass(rmp) :: clone => c_remap_data_clone + procedure, pass(rmp) :: clone => c_remap_data_clone + procedure, pass(rmp) :: move_alloc => c_remap_move_alloc end type amg_c_remap_data_type type amg_c_onelev_type @@ -227,7 +228,7 @@ module amg_c_onelev_mod procedure, pass(lv) :: get_wrksz => c_base_onelev_get_wrksize procedure, pass(lv) :: allocate_wrk => c_base_onelev_allocate_wrk procedure, pass(lv) :: free_wrk => c_base_onelev_free_wrk - procedure, nopass :: stringval => amg_stringval + procedure, nopass :: stringval => amg_stringval procedure, pass(lv) :: move_alloc => c_base_onelev_move_alloc @@ -639,7 +640,7 @@ contains ! Arguments class(amg_c_onelev_type), target, intent(inout) :: lv class(amg_c_onelev_type), target, intent(inout) :: lvout - integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_), intent(out) :: info info = psb_success_ if (allocated(lv%sm)) then @@ -705,6 +706,7 @@ contains if (info == psb_success_) call psb_move_alloc(lv%tprol,b%tprol,info) if (info == psb_success_) call psb_move_alloc(lv%desc_ac,b%desc_ac,info) if (info == psb_success_) call psb_move_alloc(lv%linmap,b%linmap,info) + if (info == psb_success_) call lv%remap_data%move_alloc(b%remap_data,info) b%base_a => lv%base_a b%base_desc => lv%base_desc @@ -759,6 +761,7 @@ contains info = psb_success_ nwv = lv%get_wrksz() if (.not.allocated(lv%wrk)) allocate(lv%wrk,stat=info) +!!$ write(0,*) 'From allocate_wrk :',lv%remap_data%desc_ac_pre_remap%is_asb() if (info == 0) then if (lv%remap_data%desc_ac_pre_remap%is_asb()) then ! @@ -806,11 +809,12 @@ contains info = psb_success_ call wk%free(info) - write(0,*) 'wrk_alloc D: "',trim(desc%get_fmt()),'"', present(desc2),desc%is_valid() +!!$ write(0,*) 'wrk_alloc D: "',trim(desc%get_fmt()),'"',& +!!$ & present(desc2),desc%is_valid() allocate(wk%wv(nwv),stat=info) if (present(desc2).and.(desc%is_valid())) then - write(0,*) 'wrk_alloc D2:',desc2%get_fmt(),desc2%is_asb() +!!$ write(0,*) 'wrk_alloc D2:',desc2%get_fmt(),desc2%is_asb() if (desc2%get_local_cols()>desc%get_local_cols()) then call inner_do_wrk_alloc(wk,nwv,desc2,vmold=vmold) else @@ -996,4 +1000,25 @@ contains call psb_safe_ab_cpy(rmp%nrsrc,remap_out%nrsrc,info) end subroutine c_remap_data_clone + subroutine c_remap_move_alloc(rmp, remap_out, info) + use psb_base_mod + implicit none + ! Arguments + class(amg_c_remap_data_type), target, intent(inout) :: rmp + class(amg_c_remap_data_type), target, intent(inout) :: remap_out + integer(psb_ipk_), intent(out) :: info + ! + integer(psb_ipk_) :: i + + info = psb_success_ + + call psb_move_alloc(rmp%ac_pre_remap,remap_out%ac_pre_remap,info) + if (info == psb_success_) & + & call psb_move_alloc(rmp%desc_ac_pre_remap,remap_out%desc_ac_pre_remap,info) + remap_out%idest = rmp%idest + call move_alloc(rmp%isrc,remap_out%isrc) + call move_alloc(rmp%nrsrc,remap_out%nrsrc) + call move_alloc(rmp%naggr,remap_out%naggr) + end subroutine c_remap_move_alloc + end module amg_c_onelev_mod diff --git a/amgprec/amg_c_prec_type.f90 b/amgprec/amg_c_prec_type.f90 index 6cb59939..3b0f62ce 100644 --- a/amgprec/amg_c_prec_type.f90 +++ b/amgprec/amg_c_prec_type.f90 @@ -102,6 +102,7 @@ module amg_c_prec_type ! The multilevel hierarchy ! type(amg_c_onelev_type), allocatable :: precv(:) + integer(psb_ipk_) :: nlevs contains procedure, pass(prec) :: psb_c_apply2_vect => amg_c_apply2_vect procedure, pass(prec) :: psb_c_apply1_vect => amg_c_apply1_vect @@ -119,6 +120,7 @@ module amg_c_prec_type procedure, pass(prec) :: cmp_complexity => amg_c_cmp_compl procedure, pass(prec) :: get_avg_cr => amg_c_get_avg_cr procedure, pass(prec) :: cmp_avg_cr => amg_c_cmp_avg_cr + procedure, pass(prec) :: set_nlevs => amg_c_set_nlevs procedure, pass(prec) :: get_nlevs => amg_c_get_nlevs procedure, pass(prec) :: get_nzeros => amg_c_get_nzeros procedure, pass(prec) :: sizeof => amg_cprec_sizeof @@ -435,10 +437,19 @@ contains class(amg_cprec_type), intent(in) :: prec integer(psb_ipk_) :: val val = 0 - if (allocated(prec%precv)) then - val = size(prec%precv) - end if +!!$ if (allocated(prec%precv)) then +!!$ val = size(prec%precv) +!!$ end if + val = prec%nlevs +!!$ write(0,*) ' NLEVS: ',prec%nlevs, val,size(prec%precv) end function amg_c_get_nlevs + + subroutine amg_c_set_nlevs(prec,nl) + implicit none + class(amg_cprec_type), intent(inout) :: prec + integer(psb_ipk_) :: nl + prec%nlevs = nl + end subroutine amg_c_set_nlevs ! ! Function returning the size of the amg_prec_type data structure ! in bytes or in number of nonzeros of the operator(s) involved. @@ -509,7 +520,7 @@ contains real(psb_spk_) :: num, den, nmin type(psb_ctxt_type) :: ctxt - integer(psb_ipk_) :: il + integer(psb_ipk_) :: il,nl num = -sone den = sone @@ -519,7 +530,10 @@ contains num = prec%precv(il)%base_a%get_nzeros() if (num >= szero) then den = num - do il=2,size(prec%precv) + nl = prec%get_nlevs() +!!$ write(0,*) 'Inside cmp_compl ',nl,size(prec%precv) + do il=2, nl +!!$ write(0,*) ' ',il,associated(prec%precv(il)%base_a) num = num + max(0,prec%precv(il)%base_a%get_nzeros()) end do end if @@ -550,7 +564,6 @@ contains end function amg_c_get_avg_cr subroutine amg_c_cmp_avg_cr(prec) - implicit none class(amg_cprec_type), intent(inout) :: prec @@ -558,17 +571,14 @@ contains type(psb_ctxt_type) :: ctxt integer(psb_ipk_) :: il, nl, iam, np - avgcr = szero ctxt = prec%ctxt call psb_info(ctxt,iam,np) - if (allocated(prec%precv)) then - nl = size(prec%precv) - do il=2,nl - avgcr = avgcr + max(szero,prec%precv(il)%szratio) - end do - avgcr = avgcr / (nl-1) - end if + nl = prec%get_nlevs() + do il=2,nl + avgcr = avgcr + max(szero,prec%precv(il)%szratio) + end do + avgcr = avgcr / (nl-1) call psb_sum(ctxt,avgcr) prec%ag_data%avg_cr = avgcr/np end subroutine amg_c_cmp_avg_cr @@ -586,9 +596,7 @@ contains ! error code. ! subroutine amg_cprecfree(p,info) - implicit none - ! Arguments type(amg_cprec_type), intent(inout) :: p integer(psb_ipk_), intent(out) :: info @@ -614,9 +622,7 @@ contains end subroutine amg_cprecfree subroutine amg_c_prec_free(prec,info) - implicit none - ! Arguments class(amg_cprec_type), intent(inout) :: prec integer(psb_ipk_), intent(out) :: info @@ -648,9 +654,7 @@ contains end subroutine amg_c_prec_free subroutine amg_c_smoothers_free(prec,info) - implicit none - ! Arguments class(amg_cprec_type), intent(inout) :: prec integer(psb_ipk_), intent(out) :: info @@ -667,7 +671,7 @@ contains end if if (allocated(prec%precv)) then - do i=1,size(prec%precv) + do i=1,prec%get_nlevs() call prec%precv(i)%free_smoothers(info) end do end if @@ -681,9 +685,7 @@ contains end subroutine amg_c_smoothers_free subroutine amg_c_hierarchy_free(prec,info) - implicit none - ! Arguments class(amg_cprec_type), intent(inout) :: prec integer(psb_ipk_), intent(out) :: info @@ -710,7 +712,6 @@ contains end subroutine amg_c_hierarchy_free - ! ! Top level methods. ! @@ -775,7 +776,6 @@ contains end subroutine amg_c_apply1_vect - subroutine amg_c_apply2v(prec,x,y,desc_data,info,trans,work) implicit none type(psb_desc_type),intent(in) :: desc_data @@ -840,7 +840,6 @@ contains subroutine amg_c_dump(prec,info,istart,iend,iproc,prefix,head,& & ac,rp,smoother,solver,tprol,& & global_num) - implicit none class(amg_cprec_type), intent(in) :: prec integer(psb_ipk_), intent(out) :: info @@ -857,7 +856,7 @@ contains info = 0 ctxt = prec%ctxt call psb_info(ctxt,iam,np) - iln = size(prec%precv) + iln = prec%get_nlevs() if (present(istart)) then il1 = max(1,istart) else @@ -881,7 +880,6 @@ contains end subroutine amg_c_dump subroutine amg_c_cnv(prec,info,amold,vmold,imold) - implicit none class(amg_cprec_type), intent(inout) :: prec integer(psb_ipk_), intent(out) :: info @@ -893,7 +891,7 @@ contains info = psb_success_ if (allocated(prec%precv)) then - do i=1,size(prec%precv) + do i=1,prec%get_nlevs() if (info == psb_success_ ) & & call prec%precv(i)%cnv(info,amold=amold,vmold=vmold,imold=imold) end do @@ -902,7 +900,6 @@ contains end subroutine amg_c_cnv subroutine amg_c_clone(prec,precout,info) - implicit none class(amg_cprec_type), intent(inout) :: prec class(psb_cprec_type), intent(inout) :: precout @@ -914,7 +911,6 @@ contains end subroutine amg_c_clone subroutine amg_c_inner_clone(prec,precout,info) - implicit none class(amg_cprec_type), intent(inout) :: prec class(psb_cprec_type), target, intent(inout) :: precout @@ -930,8 +926,9 @@ contains pout%ctxt = prec%ctxt pout%ag_data = prec%ag_data pout%outer_sweeps = prec%outer_sweeps + pout%nlevs = prec%nlevs if (allocated(prec%precv)) then - ln = size(prec%precv) + ln = prec%get_nlevs() allocate(pout%precv(ln),stat=info) if (info /= psb_success_) goto 9999 if (ln >= 1) then @@ -939,6 +936,7 @@ contains end if do lev=2, ln if (info /= psb_success_) exit +!!$ write(0,*) 'Inner_clone must be checked and reimplemented! ' call prec%precv(lev)%clone(pout%precv(lev),info) if (info == psb_success_) then pout%precv(lev)%base_a => pout%precv(lev)%ac @@ -1018,7 +1016,7 @@ contains if (psb_errstatus_fatal()) then info = psb_err_internal_error_; goto 9999 end if - nlev = size(prec%precv) + nlev = prec%get_nlevs() level = 1 do level = 1, nlev call prec%precv(level)%allocate_wrk(info,vmold=vmold) @@ -1041,7 +1039,6 @@ contains subroutine amg_c_free_wrk(prec,info) use psb_base_mod implicit none - ! Arguments class(amg_cprec_type), intent(inout) :: prec integer(psb_ipk_), intent(out) :: info @@ -1058,7 +1055,7 @@ contains end if if (allocated(prec%precv)) then - nlev = size(prec%precv) + nlev = prec%get_nlevs() do level = 1, nlev call prec%precv(level)%free_wrk(info) end do diff --git a/amgprec/amg_d_base_aggregator_mod.f90 b/amgprec/amg_d_base_aggregator_mod.f90 index c8fd2c4c..ba429532 100644 --- a/amgprec/amg_d_base_aggregator_mod.f90 +++ b/amgprec/amg_d_base_aggregator_mod.f90 @@ -89,7 +89,7 @@ module amg_d_base_aggregator_mod procedure, pass(ag) :: bld_tprol => amg_d_base_aggregator_build_tprol procedure, pass(ag) :: mat_bld => amg_d_base_aggregator_mat_bld procedure, pass(ag) :: mat_asb => amg_d_base_aggregator_mat_asb - procedure, pass(ag) :: bld_map => amg_d_base_aggregator_bld_map + procedure, pass(ag) :: bld_linmap => amg_d_base_aggregator_bld_linmap procedure, pass(ag) :: update_next => amg_d_base_aggregator_update_next procedure, pass(ag) :: clone => amg_d_base_aggregator_clone procedure, pass(ag) :: free => amg_d_base_aggregator_free @@ -458,7 +458,7 @@ contains end subroutine amg_d_base_aggregator_mat_asb ! - !> Function bld_map + !> Function bld_linmap !! \memberof amg_d_base_aggregator_type !! \brief Build linear map between hierarchy levels !! @@ -473,7 +473,7 @@ contains !! \param map The output map !! \param info Return code !! - subroutine amg_d_base_aggregator_bld_map(ag,desc_a,desc_ac,ilaggr,nlaggr,& + subroutine amg_d_base_aggregator_bld_linmap(ag,desc_a,desc_ac,ilaggr,nlaggr,& & op_restr,op_prol,map,info) use psb_base_mod implicit none @@ -484,7 +484,7 @@ contains type(psb_dlinmap_type), intent(out) :: map integer(psb_ipk_), intent(out) :: info integer(psb_ipk_) :: err_act - character(len=20) :: name='d_base_aggregator_bld_map' + character(len=20) :: name='d_base_aggregator_bld_linmap' info = psb_success_ call psb_erractionsave(err_act) @@ -508,7 +508,6 @@ contains 9999 call psb_error_handler(err_act) return - end subroutine amg_d_base_aggregator_bld_map - + end subroutine amg_d_base_aggregator_bld_linmap end module amg_d_base_aggregator_mod diff --git a/amgprec/amg_d_onelev_mod.f90 b/amgprec/amg_d_onelev_mod.f90 index 23150b4c..7cc43b7f 100644 --- a/amgprec/amg_d_onelev_mod.f90 +++ b/amgprec/amg_d_onelev_mod.f90 @@ -182,8 +182,8 @@ module amg_d_onelev_mod integer(psb_ipk_) :: idest integer(psb_ipk_), allocatable :: isrc(:), nrsrc(:), naggr(:) contains - procedure, pass(rmp) :: clone => d_remap_data_clone - procedure, pass(rmp) :: move_alloc => d_remap_move_alloc + procedure, pass(rmp) :: clone => d_remap_data_clone + procedure, pass(rmp) :: move_alloc => d_remap_move_alloc end type amg_d_remap_data_type type amg_d_onelev_type @@ -229,7 +229,7 @@ module amg_d_onelev_mod procedure, pass(lv) :: get_wrksz => d_base_onelev_get_wrksize procedure, pass(lv) :: allocate_wrk => d_base_onelev_allocate_wrk procedure, pass(lv) :: free_wrk => d_base_onelev_free_wrk - procedure, nopass :: stringval => amg_stringval + procedure, nopass :: stringval => amg_stringval procedure, pass(lv) :: move_alloc => d_base_onelev_move_alloc @@ -641,7 +641,7 @@ contains ! Arguments class(amg_d_onelev_type), target, intent(inout) :: lv class(amg_d_onelev_type), target, intent(inout) :: lvout - integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_), intent(out) :: info info = psb_success_ if (allocated(lv%sm)) then @@ -1021,5 +1021,5 @@ contains call move_alloc(rmp%nrsrc,remap_out%nrsrc) call move_alloc(rmp%naggr,remap_out%naggr) end subroutine d_remap_move_alloc - + end module amg_d_onelev_mod diff --git a/amgprec/amg_d_parmatch_aggregator_mod.F90 b/amgprec/amg_d_parmatch_aggregator_mod.F90 index 9fe5a053..cf58dff1 100644 --- a/amgprec/amg_d_parmatch_aggregator_mod.F90 +++ b/amgprec/amg_d_parmatch_aggregator_mod.F90 @@ -136,7 +136,7 @@ module amg_d_parmatch_aggregator_mod procedure, pass(ag) :: mat_bld => amg_d_parmatch_aggregator_mat_bld procedure, pass(ag) :: mat_asb => amg_d_parmatch_aggregator_mat_asb procedure, pass(ag) :: inner_mat_asb => amg_d_parmatch_aggregator_inner_mat_asb - procedure, pass(ag) :: bld_map => amg_d_parmatch_aggregator_bld_map + procedure, pass(ag) :: bld_linmap => amg_d_parmatch_aggregator_bld_linmap procedure, pass(ag) :: csetc => amg_d_parmatch_aggr_csetc procedure, pass(ag) :: cseti => amg_d_parmatch_aggr_cseti procedure, pass(ag) :: default => amg_d_parmatch_aggr_set_default @@ -643,7 +643,7 @@ contains end select end subroutine amg_d_parmatch_aggregator_clone - subroutine amg_d_parmatch_aggregator_bld_map(ag,desc_a,desc_ac,ilaggr,nlaggr,& + subroutine amg_d_parmatch_aggregator_bld_linmap(ag,desc_a,desc_ac,ilaggr,nlaggr,& & op_restr,op_prol,map,info) use psb_base_mod implicit none @@ -654,7 +654,7 @@ contains type(psb_dlinmap_type), intent(out) :: map integer(psb_ipk_), intent(out) :: info integer(psb_ipk_) :: err_act - character(len=20) :: name='d_parmatch_aggregator_bld_map' + character(len=20) :: name='d_parmatch_aggregator_bld_linmap' info = psb_success_ call psb_erractionsave(err_act) @@ -680,5 +680,5 @@ contains 9999 call psb_error_handler(err_act) return - end subroutine amg_d_parmatch_aggregator_bld_map + end subroutine amg_d_parmatch_aggregator_bld_linmap end module amg_d_parmatch_aggregator_mod diff --git a/amgprec/amg_d_prec_type.f90 b/amgprec/amg_d_prec_type.f90 index 29335283..9e5eeda5 100644 --- a/amgprec/amg_d_prec_type.f90 +++ b/amgprec/amg_d_prec_type.f90 @@ -97,6 +97,7 @@ module amg_d_prec_type ! to keep track against what is put later in the multilevel array ! integer(psb_ipk_) :: coarse_solver = -1 + ! ! The multilevel hierarchy ! @@ -119,7 +120,7 @@ module amg_d_prec_type procedure, pass(prec) :: cmp_complexity => amg_d_cmp_compl procedure, pass(prec) :: get_avg_cr => amg_d_get_avg_cr procedure, pass(prec) :: cmp_avg_cr => amg_d_cmp_avg_cr - procedure, pass(prec) :: set_nlevs => amg_d_set_nlevs + procedure, pass(prec) :: set_nlevs => amg_d_set_nlevs procedure, pass(prec) :: get_nlevs => amg_d_get_nlevs procedure, pass(prec) :: get_nzeros => amg_d_get_nzeros procedure, pass(prec) :: sizeof => amg_dprec_sizeof @@ -142,7 +143,6 @@ module amg_d_prec_type procedure, pass(prec) :: smoothers_free => amg_d_smoothers_free procedure, pass(prec) :: descr => amg_dfile_prec_descr procedure, pass(prec) :: memory_use => amg_dfile_prec_memory_use - end type amg_dprec_type private :: amg_d_dump, amg_d_get_compl, amg_d_cmp_compl,& @@ -442,7 +442,6 @@ contains !!$ end if val = prec%nlevs !!$ write(0,*) ' NLEVS: ',prec%nlevs, val,size(prec%precv) - end function amg_d_get_nlevs subroutine amg_d_set_nlevs(prec,nl) @@ -565,7 +564,6 @@ contains end function amg_d_get_avg_cr subroutine amg_d_cmp_avg_cr(prec) - implicit none class(amg_dprec_type), intent(inout) :: prec @@ -573,7 +571,6 @@ contains type(psb_ctxt_type) :: ctxt integer(psb_ipk_) :: il, nl, iam, np - avgcr = dzero ctxt = prec%ctxt call psb_info(ctxt,iam,np) @@ -599,9 +596,7 @@ contains ! error code. ! subroutine amg_dprecfree(p,info) - implicit none - ! Arguments type(amg_dprec_type), intent(inout) :: p integer(psb_ipk_), intent(out) :: info @@ -690,9 +685,7 @@ contains end subroutine amg_d_smoothers_free subroutine amg_d_hierarchy_free(prec,info) - implicit none - ! Arguments class(amg_dprec_type), intent(inout) :: prec integer(psb_ipk_), intent(out) :: info @@ -719,7 +712,6 @@ contains end subroutine amg_d_hierarchy_free - ! ! Top level methods. ! @@ -784,7 +776,6 @@ contains end subroutine amg_d_apply1_vect - subroutine amg_d_apply2v(prec,x,y,desc_data,info,trans,work) implicit none type(psb_desc_type),intent(in) :: desc_data @@ -849,7 +840,6 @@ contains subroutine amg_d_dump(prec,info,istart,iend,iproc,prefix,head,& & ac,rp,smoother,solver,tprol,& & global_num) - implicit none class(amg_dprec_type), intent(in) :: prec integer(psb_ipk_), intent(out) :: info @@ -890,7 +880,6 @@ contains end subroutine amg_d_dump subroutine amg_d_cnv(prec,info,amold,vmold,imold) - implicit none class(amg_dprec_type), intent(inout) :: prec integer(psb_ipk_), intent(out) :: info @@ -911,7 +900,6 @@ contains end subroutine amg_d_cnv subroutine amg_d_clone(prec,precout,info) - implicit none class(amg_dprec_type), intent(inout) :: prec class(psb_dprec_type), intent(inout) :: precout @@ -948,7 +936,7 @@ contains end if do lev=2, ln if (info /= psb_success_) exit - write(0,*) 'Inner_clone must be checked and reimplemented! ' +!!$ write(0,*) 'Inner_clone must be checked and reimplemented! ' call prec%precv(lev)%clone(pout%precv(lev),info) if (info == psb_success_) then pout%precv(lev)%base_a => pout%precv(lev)%ac @@ -1007,6 +995,7 @@ contains subroutine amg_d_allocate_wrk(prec,info,vmold,desc) use psb_base_mod implicit none + ! Arguments class(amg_dprec_type), intent(inout) :: prec integer(psb_ipk_), intent(out) :: info @@ -1050,7 +1039,6 @@ contains subroutine amg_d_free_wrk(prec,info) use psb_base_mod implicit none - ! Arguments class(amg_dprec_type), intent(inout) :: prec integer(psb_ipk_), intent(out) :: info @@ -1084,6 +1072,7 @@ contains function amg_d_is_allocated_wrk(prec) result(res) use psb_base_mod implicit none + ! Arguments class(amg_dprec_type), intent(in) :: prec logical :: res diff --git a/amgprec/amg_s_base_aggregator_mod.f90 b/amgprec/amg_s_base_aggregator_mod.f90 index 44eef5da..07f855e4 100644 --- a/amgprec/amg_s_base_aggregator_mod.f90 +++ b/amgprec/amg_s_base_aggregator_mod.f90 @@ -89,7 +89,7 @@ module amg_s_base_aggregator_mod procedure, pass(ag) :: bld_tprol => amg_s_base_aggregator_build_tprol procedure, pass(ag) :: mat_bld => amg_s_base_aggregator_mat_bld procedure, pass(ag) :: mat_asb => amg_s_base_aggregator_mat_asb - procedure, pass(ag) :: bld_map => amg_s_base_aggregator_bld_map + procedure, pass(ag) :: bld_linmap => amg_s_base_aggregator_bld_linmap procedure, pass(ag) :: update_next => amg_s_base_aggregator_update_next procedure, pass(ag) :: clone => amg_s_base_aggregator_clone procedure, pass(ag) :: free => amg_s_base_aggregator_free @@ -458,7 +458,7 @@ contains end subroutine amg_s_base_aggregator_mat_asb ! - !> Function bld_map + !> Function bld_linmap !! \memberof amg_s_base_aggregator_type !! \brief Build linear map between hierarchy levels !! @@ -473,7 +473,7 @@ contains !! \param map The output map !! \param info Return code !! - subroutine amg_s_base_aggregator_bld_map(ag,desc_a,desc_ac,ilaggr,nlaggr,& + subroutine amg_s_base_aggregator_bld_linmap(ag,desc_a,desc_ac,ilaggr,nlaggr,& & op_restr,op_prol,map,info) use psb_base_mod implicit none @@ -484,7 +484,7 @@ contains type(psb_slinmap_type), intent(out) :: map integer(psb_ipk_), intent(out) :: info integer(psb_ipk_) :: err_act - character(len=20) :: name='s_base_aggregator_bld_map' + character(len=20) :: name='s_base_aggregator_bld_linmap' info = psb_success_ call psb_erractionsave(err_act) @@ -508,7 +508,6 @@ contains 9999 call psb_error_handler(err_act) return - end subroutine amg_s_base_aggregator_bld_map - + end subroutine amg_s_base_aggregator_bld_linmap end module amg_s_base_aggregator_mod diff --git a/amgprec/amg_s_onelev_mod.f90 b/amgprec/amg_s_onelev_mod.f90 index dd5cff75..8f8b12b7 100644 --- a/amgprec/amg_s_onelev_mod.f90 +++ b/amgprec/amg_s_onelev_mod.f90 @@ -182,7 +182,8 @@ module amg_s_onelev_mod integer(psb_ipk_) :: idest integer(psb_ipk_), allocatable :: isrc(:), nrsrc(:), naggr(:) contains - procedure, pass(rmp) :: clone => s_remap_data_clone + procedure, pass(rmp) :: clone => s_remap_data_clone + procedure, pass(rmp) :: move_alloc => s_remap_move_alloc end type amg_s_remap_data_type type amg_s_onelev_type @@ -228,7 +229,7 @@ module amg_s_onelev_mod procedure, pass(lv) :: get_wrksz => s_base_onelev_get_wrksize procedure, pass(lv) :: allocate_wrk => s_base_onelev_allocate_wrk procedure, pass(lv) :: free_wrk => s_base_onelev_free_wrk - procedure, nopass :: stringval => amg_stringval + procedure, nopass :: stringval => amg_stringval procedure, pass(lv) :: move_alloc => s_base_onelev_move_alloc @@ -640,7 +641,7 @@ contains ! Arguments class(amg_s_onelev_type), target, intent(inout) :: lv class(amg_s_onelev_type), target, intent(inout) :: lvout - integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_), intent(out) :: info info = psb_success_ if (allocated(lv%sm)) then @@ -706,6 +707,7 @@ contains if (info == psb_success_) call psb_move_alloc(lv%tprol,b%tprol,info) if (info == psb_success_) call psb_move_alloc(lv%desc_ac,b%desc_ac,info) if (info == psb_success_) call psb_move_alloc(lv%linmap,b%linmap,info) + if (info == psb_success_) call lv%remap_data%move_alloc(b%remap_data,info) b%base_a => lv%base_a b%base_desc => lv%base_desc @@ -760,6 +762,7 @@ contains info = psb_success_ nwv = lv%get_wrksz() if (.not.allocated(lv%wrk)) allocate(lv%wrk,stat=info) +!!$ write(0,*) 'From allocate_wrk :',lv%remap_data%desc_ac_pre_remap%is_asb() if (info == 0) then if (lv%remap_data%desc_ac_pre_remap%is_asb()) then ! @@ -807,11 +810,12 @@ contains info = psb_success_ call wk%free(info) - write(0,*) 'wrk_alloc D: "',trim(desc%get_fmt()),'"', present(desc2),desc%is_valid() +!!$ write(0,*) 'wrk_alloc D: "',trim(desc%get_fmt()),'"',& +!!$ & present(desc2),desc%is_valid() allocate(wk%wv(nwv),stat=info) if (present(desc2).and.(desc%is_valid())) then - write(0,*) 'wrk_alloc D2:',desc2%get_fmt(),desc2%is_asb() +!!$ write(0,*) 'wrk_alloc D2:',desc2%get_fmt(),desc2%is_asb() if (desc2%get_local_cols()>desc%get_local_cols()) then call inner_do_wrk_alloc(wk,nwv,desc2,vmold=vmold) else @@ -997,4 +1001,25 @@ contains call psb_safe_ab_cpy(rmp%nrsrc,remap_out%nrsrc,info) end subroutine s_remap_data_clone + subroutine s_remap_move_alloc(rmp, remap_out, info) + use psb_base_mod + implicit none + ! Arguments + class(amg_s_remap_data_type), target, intent(inout) :: rmp + class(amg_s_remap_data_type), target, intent(inout) :: remap_out + integer(psb_ipk_), intent(out) :: info + ! + integer(psb_ipk_) :: i + + info = psb_success_ + + call psb_move_alloc(rmp%ac_pre_remap,remap_out%ac_pre_remap,info) + if (info == psb_success_) & + & call psb_move_alloc(rmp%desc_ac_pre_remap,remap_out%desc_ac_pre_remap,info) + remap_out%idest = rmp%idest + call move_alloc(rmp%isrc,remap_out%isrc) + call move_alloc(rmp%nrsrc,remap_out%nrsrc) + call move_alloc(rmp%naggr,remap_out%naggr) + end subroutine s_remap_move_alloc + end module amg_s_onelev_mod diff --git a/amgprec/amg_s_parmatch_aggregator_mod.F90 b/amgprec/amg_s_parmatch_aggregator_mod.F90 index 08c48cbc..bedb204f 100644 --- a/amgprec/amg_s_parmatch_aggregator_mod.F90 +++ b/amgprec/amg_s_parmatch_aggregator_mod.F90 @@ -136,7 +136,7 @@ module amg_s_parmatch_aggregator_mod procedure, pass(ag) :: mat_bld => amg_s_parmatch_aggregator_mat_bld procedure, pass(ag) :: mat_asb => amg_s_parmatch_aggregator_mat_asb procedure, pass(ag) :: inner_mat_asb => amg_s_parmatch_aggregator_inner_mat_asb - procedure, pass(ag) :: bld_map => amg_s_parmatch_aggregator_bld_map + procedure, pass(ag) :: bld_linmap => amg_s_parmatch_aggregator_bld_linmap procedure, pass(ag) :: csetc => amg_s_parmatch_aggr_csetc procedure, pass(ag) :: cseti => amg_s_parmatch_aggr_cseti procedure, pass(ag) :: default => amg_s_parmatch_aggr_set_default @@ -643,7 +643,7 @@ contains end select end subroutine amg_s_parmatch_aggregator_clone - subroutine amg_s_parmatch_aggregator_bld_map(ag,desc_a,desc_ac,ilaggr,nlaggr,& + subroutine amg_s_parmatch_aggregator_bld_linmap(ag,desc_a,desc_ac,ilaggr,nlaggr,& & op_restr,op_prol,map,info) use psb_base_mod implicit none @@ -654,7 +654,7 @@ contains type(psb_slinmap_type), intent(out) :: map integer(psb_ipk_), intent(out) :: info integer(psb_ipk_) :: err_act - character(len=20) :: name='s_parmatch_aggregator_bld_map' + character(len=20) :: name='s_parmatch_aggregator_bld_linmap' info = psb_success_ call psb_erractionsave(err_act) @@ -680,5 +680,5 @@ contains 9999 call psb_error_handler(err_act) return - end subroutine amg_s_parmatch_aggregator_bld_map + end subroutine amg_s_parmatch_aggregator_bld_linmap end module amg_s_parmatch_aggregator_mod diff --git a/amgprec/amg_s_prec_type.f90 b/amgprec/amg_s_prec_type.f90 index 68ff409b..6574fe12 100644 --- a/amgprec/amg_s_prec_type.f90 +++ b/amgprec/amg_s_prec_type.f90 @@ -102,6 +102,7 @@ module amg_s_prec_type ! The multilevel hierarchy ! type(amg_s_onelev_type), allocatable :: precv(:) + integer(psb_ipk_) :: nlevs contains procedure, pass(prec) :: psb_s_apply2_vect => amg_s_apply2_vect procedure, pass(prec) :: psb_s_apply1_vect => amg_s_apply1_vect @@ -119,6 +120,7 @@ module amg_s_prec_type procedure, pass(prec) :: cmp_complexity => amg_s_cmp_compl procedure, pass(prec) :: get_avg_cr => amg_s_get_avg_cr procedure, pass(prec) :: cmp_avg_cr => amg_s_cmp_avg_cr + procedure, pass(prec) :: set_nlevs => amg_s_set_nlevs procedure, pass(prec) :: get_nlevs => amg_s_get_nlevs procedure, pass(prec) :: get_nzeros => amg_s_get_nzeros procedure, pass(prec) :: sizeof => amg_sprec_sizeof @@ -435,10 +437,19 @@ contains class(amg_sprec_type), intent(in) :: prec integer(psb_ipk_) :: val val = 0 - if (allocated(prec%precv)) then - val = size(prec%precv) - end if +!!$ if (allocated(prec%precv)) then +!!$ val = size(prec%precv) +!!$ end if + val = prec%nlevs +!!$ write(0,*) ' NLEVS: ',prec%nlevs, val,size(prec%precv) end function amg_s_get_nlevs + + subroutine amg_s_set_nlevs(prec,nl) + implicit none + class(amg_sprec_type), intent(inout) :: prec + integer(psb_ipk_) :: nl + prec%nlevs = nl + end subroutine amg_s_set_nlevs ! ! Function returning the size of the amg_prec_type data structure ! in bytes or in number of nonzeros of the operator(s) involved. @@ -509,7 +520,7 @@ contains real(psb_spk_) :: num, den, nmin type(psb_ctxt_type) :: ctxt - integer(psb_ipk_) :: il + integer(psb_ipk_) :: il,nl num = -sone den = sone @@ -519,7 +530,10 @@ contains num = prec%precv(il)%base_a%get_nzeros() if (num >= szero) then den = num - do il=2,size(prec%precv) + nl = prec%get_nlevs() +!!$ write(0,*) 'Inside cmp_compl ',nl,size(prec%precv) + do il=2, nl +!!$ write(0,*) ' ',il,associated(prec%precv(il)%base_a) num = num + max(0,prec%precv(il)%base_a%get_nzeros()) end do end if @@ -550,7 +564,6 @@ contains end function amg_s_get_avg_cr subroutine amg_s_cmp_avg_cr(prec) - implicit none class(amg_sprec_type), intent(inout) :: prec @@ -558,17 +571,14 @@ contains type(psb_ctxt_type) :: ctxt integer(psb_ipk_) :: il, nl, iam, np - avgcr = szero ctxt = prec%ctxt call psb_info(ctxt,iam,np) - if (allocated(prec%precv)) then - nl = size(prec%precv) - do il=2,nl - avgcr = avgcr + max(szero,prec%precv(il)%szratio) - end do - avgcr = avgcr / (nl-1) - end if + nl = prec%get_nlevs() + do il=2,nl + avgcr = avgcr + max(szero,prec%precv(il)%szratio) + end do + avgcr = avgcr / (nl-1) call psb_sum(ctxt,avgcr) prec%ag_data%avg_cr = avgcr/np end subroutine amg_s_cmp_avg_cr @@ -586,9 +596,7 @@ contains ! error code. ! subroutine amg_sprecfree(p,info) - implicit none - ! Arguments type(amg_sprec_type), intent(inout) :: p integer(psb_ipk_), intent(out) :: info @@ -614,9 +622,7 @@ contains end subroutine amg_sprecfree subroutine amg_s_prec_free(prec,info) - implicit none - ! Arguments class(amg_sprec_type), intent(inout) :: prec integer(psb_ipk_), intent(out) :: info @@ -648,9 +654,7 @@ contains end subroutine amg_s_prec_free subroutine amg_s_smoothers_free(prec,info) - implicit none - ! Arguments class(amg_sprec_type), intent(inout) :: prec integer(psb_ipk_), intent(out) :: info @@ -667,7 +671,7 @@ contains end if if (allocated(prec%precv)) then - do i=1,size(prec%precv) + do i=1,prec%get_nlevs() call prec%precv(i)%free_smoothers(info) end do end if @@ -681,9 +685,7 @@ contains end subroutine amg_s_smoothers_free subroutine amg_s_hierarchy_free(prec,info) - implicit none - ! Arguments class(amg_sprec_type), intent(inout) :: prec integer(psb_ipk_), intent(out) :: info @@ -710,7 +712,6 @@ contains end subroutine amg_s_hierarchy_free - ! ! Top level methods. ! @@ -775,7 +776,6 @@ contains end subroutine amg_s_apply1_vect - subroutine amg_s_apply2v(prec,x,y,desc_data,info,trans,work) implicit none type(psb_desc_type),intent(in) :: desc_data @@ -840,7 +840,6 @@ contains subroutine amg_s_dump(prec,info,istart,iend,iproc,prefix,head,& & ac,rp,smoother,solver,tprol,& & global_num) - implicit none class(amg_sprec_type), intent(in) :: prec integer(psb_ipk_), intent(out) :: info @@ -857,7 +856,7 @@ contains info = 0 ctxt = prec%ctxt call psb_info(ctxt,iam,np) - iln = size(prec%precv) + iln = prec%get_nlevs() if (present(istart)) then il1 = max(1,istart) else @@ -881,7 +880,6 @@ contains end subroutine amg_s_dump subroutine amg_s_cnv(prec,info,amold,vmold,imold) - implicit none class(amg_sprec_type), intent(inout) :: prec integer(psb_ipk_), intent(out) :: info @@ -893,7 +891,7 @@ contains info = psb_success_ if (allocated(prec%precv)) then - do i=1,size(prec%precv) + do i=1,prec%get_nlevs() if (info == psb_success_ ) & & call prec%precv(i)%cnv(info,amold=amold,vmold=vmold,imold=imold) end do @@ -902,7 +900,6 @@ contains end subroutine amg_s_cnv subroutine amg_s_clone(prec,precout,info) - implicit none class(amg_sprec_type), intent(inout) :: prec class(psb_sprec_type), intent(inout) :: precout @@ -914,7 +911,6 @@ contains end subroutine amg_s_clone subroutine amg_s_inner_clone(prec,precout,info) - implicit none class(amg_sprec_type), intent(inout) :: prec class(psb_sprec_type), target, intent(inout) :: precout @@ -930,8 +926,9 @@ contains pout%ctxt = prec%ctxt pout%ag_data = prec%ag_data pout%outer_sweeps = prec%outer_sweeps + pout%nlevs = prec%nlevs if (allocated(prec%precv)) then - ln = size(prec%precv) + ln = prec%get_nlevs() allocate(pout%precv(ln),stat=info) if (info /= psb_success_) goto 9999 if (ln >= 1) then @@ -939,6 +936,7 @@ contains end if do lev=2, ln if (info /= psb_success_) exit +!!$ write(0,*) 'Inner_clone must be checked and reimplemented! ' call prec%precv(lev)%clone(pout%precv(lev),info) if (info == psb_success_) then pout%precv(lev)%base_a => pout%precv(lev)%ac @@ -1018,7 +1016,7 @@ contains if (psb_errstatus_fatal()) then info = psb_err_internal_error_; goto 9999 end if - nlev = size(prec%precv) + nlev = prec%get_nlevs() level = 1 do level = 1, nlev call prec%precv(level)%allocate_wrk(info,vmold=vmold) @@ -1041,7 +1039,6 @@ contains subroutine amg_s_free_wrk(prec,info) use psb_base_mod implicit none - ! Arguments class(amg_sprec_type), intent(inout) :: prec integer(psb_ipk_), intent(out) :: info @@ -1058,7 +1055,7 @@ contains end if if (allocated(prec%precv)) then - nlev = size(prec%precv) + nlev = prec%get_nlevs() do level = 1, nlev call prec%precv(level)%free_wrk(info) end do diff --git a/amgprec/amg_z_base_aggregator_mod.f90 b/amgprec/amg_z_base_aggregator_mod.f90 index 6157e660..fbf33c6f 100644 --- a/amgprec/amg_z_base_aggregator_mod.f90 +++ b/amgprec/amg_z_base_aggregator_mod.f90 @@ -89,7 +89,7 @@ module amg_z_base_aggregator_mod procedure, pass(ag) :: bld_tprol => amg_z_base_aggregator_build_tprol procedure, pass(ag) :: mat_bld => amg_z_base_aggregator_mat_bld procedure, pass(ag) :: mat_asb => amg_z_base_aggregator_mat_asb - procedure, pass(ag) :: bld_map => amg_z_base_aggregator_bld_map + procedure, pass(ag) :: bld_linmap => amg_z_base_aggregator_bld_linmap procedure, pass(ag) :: update_next => amg_z_base_aggregator_update_next procedure, pass(ag) :: clone => amg_z_base_aggregator_clone procedure, pass(ag) :: free => amg_z_base_aggregator_free @@ -458,7 +458,7 @@ contains end subroutine amg_z_base_aggregator_mat_asb ! - !> Function bld_map + !> Function bld_linmap !! \memberof amg_z_base_aggregator_type !! \brief Build linear map between hierarchy levels !! @@ -473,7 +473,7 @@ contains !! \param map The output map !! \param info Return code !! - subroutine amg_z_base_aggregator_bld_map(ag,desc_a,desc_ac,ilaggr,nlaggr,& + subroutine amg_z_base_aggregator_bld_linmap(ag,desc_a,desc_ac,ilaggr,nlaggr,& & op_restr,op_prol,map,info) use psb_base_mod implicit none @@ -484,7 +484,7 @@ contains type(psb_zlinmap_type), intent(out) :: map integer(psb_ipk_), intent(out) :: info integer(psb_ipk_) :: err_act - character(len=20) :: name='z_base_aggregator_bld_map' + character(len=20) :: name='z_base_aggregator_bld_linmap' info = psb_success_ call psb_erractionsave(err_act) @@ -508,7 +508,6 @@ contains 9999 call psb_error_handler(err_act) return - end subroutine amg_z_base_aggregator_bld_map - + end subroutine amg_z_base_aggregator_bld_linmap end module amg_z_base_aggregator_mod diff --git a/amgprec/amg_z_onelev_mod.f90 b/amgprec/amg_z_onelev_mod.f90 index 7c54f113..309bb6a4 100644 --- a/amgprec/amg_z_onelev_mod.f90 +++ b/amgprec/amg_z_onelev_mod.f90 @@ -181,7 +181,8 @@ module amg_z_onelev_mod integer(psb_ipk_) :: idest integer(psb_ipk_), allocatable :: isrc(:), nrsrc(:), naggr(:) contains - procedure, pass(rmp) :: clone => z_remap_data_clone + procedure, pass(rmp) :: clone => z_remap_data_clone + procedure, pass(rmp) :: move_alloc => z_remap_move_alloc end type amg_z_remap_data_type type amg_z_onelev_type @@ -227,7 +228,7 @@ module amg_z_onelev_mod procedure, pass(lv) :: get_wrksz => z_base_onelev_get_wrksize procedure, pass(lv) :: allocate_wrk => z_base_onelev_allocate_wrk procedure, pass(lv) :: free_wrk => z_base_onelev_free_wrk - procedure, nopass :: stringval => amg_stringval + procedure, nopass :: stringval => amg_stringval procedure, pass(lv) :: move_alloc => z_base_onelev_move_alloc @@ -639,7 +640,7 @@ contains ! Arguments class(amg_z_onelev_type), target, intent(inout) :: lv class(amg_z_onelev_type), target, intent(inout) :: lvout - integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_), intent(out) :: info info = psb_success_ if (allocated(lv%sm)) then @@ -705,6 +706,7 @@ contains if (info == psb_success_) call psb_move_alloc(lv%tprol,b%tprol,info) if (info == psb_success_) call psb_move_alloc(lv%desc_ac,b%desc_ac,info) if (info == psb_success_) call psb_move_alloc(lv%linmap,b%linmap,info) + if (info == psb_success_) call lv%remap_data%move_alloc(b%remap_data,info) b%base_a => lv%base_a b%base_desc => lv%base_desc @@ -759,6 +761,7 @@ contains info = psb_success_ nwv = lv%get_wrksz() if (.not.allocated(lv%wrk)) allocate(lv%wrk,stat=info) +!!$ write(0,*) 'From allocate_wrk :',lv%remap_data%desc_ac_pre_remap%is_asb() if (info == 0) then if (lv%remap_data%desc_ac_pre_remap%is_asb()) then ! @@ -806,11 +809,12 @@ contains info = psb_success_ call wk%free(info) - write(0,*) 'wrk_alloc D: "',trim(desc%get_fmt()),'"', present(desc2),desc%is_valid() +!!$ write(0,*) 'wrk_alloc D: "',trim(desc%get_fmt()),'"',& +!!$ & present(desc2),desc%is_valid() allocate(wk%wv(nwv),stat=info) if (present(desc2).and.(desc%is_valid())) then - write(0,*) 'wrk_alloc D2:',desc2%get_fmt(),desc2%is_asb() +!!$ write(0,*) 'wrk_alloc D2:',desc2%get_fmt(),desc2%is_asb() if (desc2%get_local_cols()>desc%get_local_cols()) then call inner_do_wrk_alloc(wk,nwv,desc2,vmold=vmold) else @@ -996,4 +1000,25 @@ contains call psb_safe_ab_cpy(rmp%nrsrc,remap_out%nrsrc,info) end subroutine z_remap_data_clone + subroutine z_remap_move_alloc(rmp, remap_out, info) + use psb_base_mod + implicit none + ! Arguments + class(amg_z_remap_data_type), target, intent(inout) :: rmp + class(amg_z_remap_data_type), target, intent(inout) :: remap_out + integer(psb_ipk_), intent(out) :: info + ! + integer(psb_ipk_) :: i + + info = psb_success_ + + call psb_move_alloc(rmp%ac_pre_remap,remap_out%ac_pre_remap,info) + if (info == psb_success_) & + & call psb_move_alloc(rmp%desc_ac_pre_remap,remap_out%desc_ac_pre_remap,info) + remap_out%idest = rmp%idest + call move_alloc(rmp%isrc,remap_out%isrc) + call move_alloc(rmp%nrsrc,remap_out%nrsrc) + call move_alloc(rmp%naggr,remap_out%naggr) + end subroutine z_remap_move_alloc + end module amg_z_onelev_mod diff --git a/amgprec/amg_z_prec_type.f90 b/amgprec/amg_z_prec_type.f90 index 7bc73992..2d8b362a 100644 --- a/amgprec/amg_z_prec_type.f90 +++ b/amgprec/amg_z_prec_type.f90 @@ -102,6 +102,7 @@ module amg_z_prec_type ! The multilevel hierarchy ! type(amg_z_onelev_type), allocatable :: precv(:) + integer(psb_ipk_) :: nlevs contains procedure, pass(prec) :: psb_z_apply2_vect => amg_z_apply2_vect procedure, pass(prec) :: psb_z_apply1_vect => amg_z_apply1_vect @@ -119,6 +120,7 @@ module amg_z_prec_type procedure, pass(prec) :: cmp_complexity => amg_z_cmp_compl procedure, pass(prec) :: get_avg_cr => amg_z_get_avg_cr procedure, pass(prec) :: cmp_avg_cr => amg_z_cmp_avg_cr + procedure, pass(prec) :: set_nlevs => amg_z_set_nlevs procedure, pass(prec) :: get_nlevs => amg_z_get_nlevs procedure, pass(prec) :: get_nzeros => amg_z_get_nzeros procedure, pass(prec) :: sizeof => amg_zprec_sizeof @@ -435,10 +437,19 @@ contains class(amg_zprec_type), intent(in) :: prec integer(psb_ipk_) :: val val = 0 - if (allocated(prec%precv)) then - val = size(prec%precv) - end if +!!$ if (allocated(prec%precv)) then +!!$ val = size(prec%precv) +!!$ end if + val = prec%nlevs +!!$ write(0,*) ' NLEVS: ',prec%nlevs, val,size(prec%precv) end function amg_z_get_nlevs + + subroutine amg_z_set_nlevs(prec,nl) + implicit none + class(amg_zprec_type), intent(inout) :: prec + integer(psb_ipk_) :: nl + prec%nlevs = nl + end subroutine amg_z_set_nlevs ! ! Function returning the size of the amg_prec_type data structure ! in bytes or in number of nonzeros of the operator(s) involved. @@ -509,7 +520,7 @@ contains real(psb_dpk_) :: num, den, nmin type(psb_ctxt_type) :: ctxt - integer(psb_ipk_) :: il + integer(psb_ipk_) :: il,nl num = -done den = done @@ -519,7 +530,10 @@ contains num = prec%precv(il)%base_a%get_nzeros() if (num >= dzero) then den = num - do il=2,size(prec%precv) + nl = prec%get_nlevs() +!!$ write(0,*) 'Inside cmp_compl ',nl,size(prec%precv) + do il=2, nl +!!$ write(0,*) ' ',il,associated(prec%precv(il)%base_a) num = num + max(0,prec%precv(il)%base_a%get_nzeros()) end do end if @@ -550,7 +564,6 @@ contains end function amg_z_get_avg_cr subroutine amg_z_cmp_avg_cr(prec) - implicit none class(amg_zprec_type), intent(inout) :: prec @@ -558,17 +571,14 @@ contains type(psb_ctxt_type) :: ctxt integer(psb_ipk_) :: il, nl, iam, np - avgcr = dzero ctxt = prec%ctxt call psb_info(ctxt,iam,np) - if (allocated(prec%precv)) then - nl = size(prec%precv) - do il=2,nl - avgcr = avgcr + max(dzero,prec%precv(il)%szratio) - end do - avgcr = avgcr / (nl-1) - end if + nl = prec%get_nlevs() + do il=2,nl + avgcr = avgcr + max(dzero,prec%precv(il)%szratio) + end do + avgcr = avgcr / (nl-1) call psb_sum(ctxt,avgcr) prec%ag_data%avg_cr = avgcr/np end subroutine amg_z_cmp_avg_cr @@ -586,9 +596,7 @@ contains ! error code. ! subroutine amg_zprecfree(p,info) - implicit none - ! Arguments type(amg_zprec_type), intent(inout) :: p integer(psb_ipk_), intent(out) :: info @@ -614,9 +622,7 @@ contains end subroutine amg_zprecfree subroutine amg_z_prec_free(prec,info) - implicit none - ! Arguments class(amg_zprec_type), intent(inout) :: prec integer(psb_ipk_), intent(out) :: info @@ -648,9 +654,7 @@ contains end subroutine amg_z_prec_free subroutine amg_z_smoothers_free(prec,info) - implicit none - ! Arguments class(amg_zprec_type), intent(inout) :: prec integer(psb_ipk_), intent(out) :: info @@ -667,7 +671,7 @@ contains end if if (allocated(prec%precv)) then - do i=1,size(prec%precv) + do i=1,prec%get_nlevs() call prec%precv(i)%free_smoothers(info) end do end if @@ -681,9 +685,7 @@ contains end subroutine amg_z_smoothers_free subroutine amg_z_hierarchy_free(prec,info) - implicit none - ! Arguments class(amg_zprec_type), intent(inout) :: prec integer(psb_ipk_), intent(out) :: info @@ -710,7 +712,6 @@ contains end subroutine amg_z_hierarchy_free - ! ! Top level methods. ! @@ -775,7 +776,6 @@ contains end subroutine amg_z_apply1_vect - subroutine amg_z_apply2v(prec,x,y,desc_data,info,trans,work) implicit none type(psb_desc_type),intent(in) :: desc_data @@ -840,7 +840,6 @@ contains subroutine amg_z_dump(prec,info,istart,iend,iproc,prefix,head,& & ac,rp,smoother,solver,tprol,& & global_num) - implicit none class(amg_zprec_type), intent(in) :: prec integer(psb_ipk_), intent(out) :: info @@ -857,7 +856,7 @@ contains info = 0 ctxt = prec%ctxt call psb_info(ctxt,iam,np) - iln = size(prec%precv) + iln = prec%get_nlevs() if (present(istart)) then il1 = max(1,istart) else @@ -881,7 +880,6 @@ contains end subroutine amg_z_dump subroutine amg_z_cnv(prec,info,amold,vmold,imold) - implicit none class(amg_zprec_type), intent(inout) :: prec integer(psb_ipk_), intent(out) :: info @@ -893,7 +891,7 @@ contains info = psb_success_ if (allocated(prec%precv)) then - do i=1,size(prec%precv) + do i=1,prec%get_nlevs() if (info == psb_success_ ) & & call prec%precv(i)%cnv(info,amold=amold,vmold=vmold,imold=imold) end do @@ -902,7 +900,6 @@ contains end subroutine amg_z_cnv subroutine amg_z_clone(prec,precout,info) - implicit none class(amg_zprec_type), intent(inout) :: prec class(psb_zprec_type), intent(inout) :: precout @@ -914,7 +911,6 @@ contains end subroutine amg_z_clone subroutine amg_z_inner_clone(prec,precout,info) - implicit none class(amg_zprec_type), intent(inout) :: prec class(psb_zprec_type), target, intent(inout) :: precout @@ -930,8 +926,9 @@ contains pout%ctxt = prec%ctxt pout%ag_data = prec%ag_data pout%outer_sweeps = prec%outer_sweeps + pout%nlevs = prec%nlevs if (allocated(prec%precv)) then - ln = size(prec%precv) + ln = prec%get_nlevs() allocate(pout%precv(ln),stat=info) if (info /= psb_success_) goto 9999 if (ln >= 1) then @@ -939,6 +936,7 @@ contains end if do lev=2, ln if (info /= psb_success_) exit +!!$ write(0,*) 'Inner_clone must be checked and reimplemented! ' call prec%precv(lev)%clone(pout%precv(lev),info) if (info == psb_success_) then pout%precv(lev)%base_a => pout%precv(lev)%ac @@ -1018,7 +1016,7 @@ contains if (psb_errstatus_fatal()) then info = psb_err_internal_error_; goto 9999 end if - nlev = size(prec%precv) + nlev = prec%get_nlevs() level = 1 do level = 1, nlev call prec%precv(level)%allocate_wrk(info,vmold=vmold) @@ -1041,7 +1039,6 @@ contains subroutine amg_z_free_wrk(prec,info) use psb_base_mod implicit none - ! Arguments class(amg_zprec_type), intent(inout) :: prec integer(psb_ipk_), intent(out) :: info @@ -1058,7 +1055,7 @@ contains end if if (allocated(prec%precv)) then - nlev = size(prec%precv) + nlev = prec%get_nlevs() do level = 1, nlev call prec%precv(level)%free_wrk(info) end do diff --git a/amgprec/impl/amg_c_hierarchy_bld.F90 b/amgprec/impl/amg_c_hierarchy_bld.F90 index 88702bba..b2d99585 100644 --- a/amgprec/impl/amg_c_hierarchy_bld.F90 +++ b/amgprec/impl/amg_c_hierarchy_bld.F90 @@ -64,11 +64,9 @@ ! Error code. ! subroutine amg_c_hierarchy_bld(a,desc_a,prec,info,cpymat) - use psb_base_mod use amg_c_inner_mod use amg_c_prec_mod, amg_protect_name => amg_c_hierarchy_bld - Implicit None ! Arguments @@ -82,7 +80,7 @@ subroutine amg_c_hierarchy_bld(a,desc_a,prec,info,cpymat) type(psb_ctxt_type) :: ctxt integer(psb_ipk_) :: me,np integer(psb_ipk_) :: err,i,k, err_act, iszv, newsz,& - & nplevs, mxplevs + & nplevs, mxplevs, level integer(psb_lpk_) :: iaggsize, casize, mncsize, mncszpp real(psb_spk_) :: mnaggratio, sizeratio, athresh, aomega class(amg_c_base_smoother_type), allocatable :: coarse_sm, med_sm, & @@ -98,6 +96,7 @@ subroutine amg_c_hierarchy_bld(a,desc_a,prec,info,cpymat) character(len=40) :: ch_err integer(psb_ipk_), save :: idx_bldtp=-1, idx_matasb=-1 logical, parameter :: do_timings=.false. + logical :: stop_hierarchy_loop type(psb_ctxt_type) :: lctxt integer(psb_ipk_) :: lme,lnp @@ -141,7 +140,7 @@ subroutine amg_c_hierarchy_bld(a,desc_a,prec,info,cpymat) mnaggratio = prec%ag_data%min_cr_ratio mncsize = prec%ag_data%min_coarse_size mncszpp = prec%ag_data%min_coarse_size_per_process - iszv = size(prec%precv) + iszv = prec%get_nlevs() call psb_bcast(ctxt,iszv) call psb_bcast(ctxt,mncsize) call psb_bcast(ctxt,mncszpp) @@ -167,7 +166,7 @@ subroutine amg_c_hierarchy_bld(a,desc_a,prec,info,cpymat) call psb_errpush(info,name,a_err='Inconsistent min_cr_ratio') goto 9999 end if - if (iszv /= size(prec%precv)) then + if (iszv /= prec%get_nlevs()) then info=psb_err_internal_error_ call psb_errpush(info,name,a_err='Inconsistent size of precv') goto 9999 @@ -182,6 +181,7 @@ subroutine amg_c_hierarchy_bld(a,desc_a,prec,info,cpymat) call psb_errpush(info,name,a_err=ch_err) goto 9999 endif + if (iszv == 1) then ! ! This is OK, since it may be called by the user even if there @@ -229,7 +229,6 @@ subroutine amg_c_hierarchy_bld(a,desc_a,prec,info,cpymat) casize = mncsize end if prec%ag_data%target_coarse_size = casize - nplevs = max(itwo,mxplevs) ! @@ -242,7 +241,7 @@ subroutine amg_c_hierarchy_bld(a,desc_a,prec,info,cpymat) goto 9999 end if ! - ! First set desired number of levels + ! First set desired number of levels if different from default. ! if (iszv /= nplevs) then allocate(tprecv(nplevs),stat=info) @@ -288,7 +287,8 @@ subroutine amg_c_hierarchy_bld(a,desc_a,prec,info,cpymat) call prec%precv(i)%free(info) end do call move_alloc(tprecv,prec%precv) - iszv = size(prec%precv) + call prec%set_nlevs(nplevs) + iszv = prec%get_nlevs() end if ! @@ -304,16 +304,23 @@ subroutine amg_c_hierarchy_bld(a,desc_a,prec,info,cpymat) call psb_cd_renum_block(desc_a,prec%precv(1)%desc_ac,info) prec%precv(1)%base_desc => prec%precv(1)%desc_ac + ! + ! Main build loop + ! newsz = 0 + stop_hierarchy_loop = .false. array_build_loop: do i=2, iszv ! ! Check on the iprcparm contents: they should be the same ! on all processes. ! call psb_bcast(ctxt,prec%precv(i)%parms) + ! + ! Get current context: might have performed remapping + ! lctxt = prec%precv(i-1)%base_desc%get_ctxt() call psb_info(lctxt,lme,lnp) - write(0,*) 'Check at level',i,lme,lnp +!!$ write(0,*) 'Check at level',i,lme,lnp ! ! Sanity checks on the parameters ! @@ -329,8 +336,8 @@ subroutine amg_c_hierarchy_bld(a,desc_a,prec,info,cpymat) & write(debug_unit,*) me,' ',trim(name),& & 'Calling mlprcbld at level ',i ! - ! Build the mapping between levels i-1 and i and the matrix - ! at level i + ! Build the tentative mapping between levels i-1 and i + ! and the matrix at level i ! if (do_timings) call psb_tic(idx_bldtp) if (info == psb_success_)& @@ -352,62 +359,26 @@ subroutine amg_c_hierarchy_bld(a,desc_a,prec,info,cpymat) ! Save op_prol just in case ! call op_prol%clone(prec%precv(i)%tprol,info) + ! ! Check for early termination of aggregation loop. - ! - iaggsize = sum(nlaggr) - - sizeratio = iaggsize - if (i==2) then - sizeratio = desc_a%get_global_rows()/sizeratio + ! + if (i == 2) then + call amg_c_hierarchy_bld_cmp_newsz(i,iszv,& + & desc_a%get_global_rows(),& + & nlaggr,casize,mnaggratio,sizeratio,newsz) else - sizeratio = sum(prec%precv(i-1)%linmap%naggr)/sizeratio + call amg_c_hierarchy_bld_cmp_newsz(i,iszv,& + & sum(prec%precv(i-1)%linmap%naggr),& + & nlaggr,casize,mnaggratio,sizeratio,newsz) end if prec%precv(i)%szratio = sizeratio - - if (iaggsize <= casize) newsz = i - if (i == iszv) newsz = i - - if (i>2) then - if (sizeratio < mnaggratio) then - if (sizeratio > 1) then - newsz = i - else - ! - ! We are not gaining - ! - newsz = i-1 - end if - end if - - ! - ! TO BE REVIWED - ! - write(0,*) me,lme,allocated(prec%precv(i-1)%linmap%naggr) - if (allocated(prec%precv(i-1)%linmap%naggr)) then - write(0,*) me,lme,size(prec%precv(i-1)%linmap%naggr), size(nlaggr) - end if - if (.false.) then - if (lme >=0) then - if (all(nlaggr == prec%precv(i-1)%linmap%naggr)) then - newsz=i-1 - if (me == 0) then - write(debug_unit,*) trim(name),& - &': Warning: aggregates from level ',& - & newsz - write(debug_unit,*) trim(name),& - &': to level ',& - & iszv,' coincide.' - write(debug_unit,*) trim(name),& - &': Number of levels actually used :',newsz - write(debug_unit,*) - end if - end if - end if - end if - end if call psb_bcast(ctxt,newsz) + ! + ! Handle reallocation, if needed, and then mat_asb to polish off the + ! construction + ! if (newsz > 0) then ! ! This is awkward, we are saving the aggregation parms, for the sake @@ -441,156 +412,112 @@ subroutine amg_c_hierarchy_bld(a,desc_a,prec,info,cpymat) & a_err=ch_err) goto 9999 endif - if (amg_get_do_remap().and.(i>=2)) then - block - type(psb_ctxt_type) :: lctxt - integer(psb_ipk_) :: lme,lnp - lctxt = prec%precv(i)%desc_ac%get_ctxt() - call psb_info(lctxt,lme,lnp) - write(0,*) ' Context on remapping ',lme,lnp - if ((lme >=0).and.(lnp>=2)) then - associate(lv=>prec%precv(i), rmp => prec%precv(i)%remap_data) - call lv%desc_ac%clone(rmp%desc_ac_pre_remap,info) - call lv%ac%clone(rmp%ac_pre_remap,info) - write(0,*) ' Doing remapping ',lnp, lnp/2 - call psb_remap(lnp/2,rmp%desc_ac_pre_remap,rmp%ac_pre_remap,& - & rmp%idest,rmp%isrc,rmp%nrsrc,rmp%naggr,lv%desc_ac,lv%ac,info) -!!$ write(0,*) me,' Out of remapping ',rmp%desc_ac_pre_remap%get_fmt(),' ',& -!!$ & lv%desc_ac%get_fmt(),sum(lv%linmap%naggr),sum(rmp%naggr) - write(0,*) 'Assignment ',size(lv%linmap%naggr),size(rmp%naggr) - lv%linmap%naggr(:) = rmp%naggr(:) - lv%linmap%p_desc_V => rmp%desc_ac_pre_remap - lv%base_a => lv%ac - lv%base_desc => lv%desc_ac - end associate - end if - end block - end if - exit array_build_loop +!!$ write(0,*) ' Early exit of array_build_loop',i,iszv,info,& + level = newsz + stop_hierarchy_loop = .true. else - if (do_timings) call psb_tic(idx_matasb) + if (do_timings) call psb_tic(idx_matasb) if (info == psb_success_) call prec%precv(i)%mat_asb(& & prec%precv(i-1)%base_a,prec%precv(i-1)%base_desc,& & ilaggr,nlaggr,op_prol,info) if (do_timings) call psb_toc(idx_matasb) + level = i end if - if (info /= psb_success_) then - write(ch_err,'(a,i7)') 'Mat asb fail @ level ',i - call psb_errpush(psb_err_internal_error_,name,& - & a_err=ch_err) - goto 9999 - endif - if (i=2)) then - block - type(psb_ctxt_type) :: lctxt - integer(psb_ipk_) :: lme,lnp - lctxt = prec%precv(i)%desc_ac%get_ctxt() - call psb_info(lctxt,lme,lnp) - write(0,*) ' Context on remapping ',lme,lnp + + ! + ! Do we want to remap onto a smaller subset of processes? + ! Will need a more sophisticated policy + ! + block + type(psb_ctxt_type) :: lctxt + integer(psb_ipk_) :: lme,lnp + lctxt = prec%precv(level)%desc_ac%get_ctxt() + call psb_info(lctxt,lme,lnp) + if (amg_c_policy_do_remap(lctxt,level,sum(nlaggr))) then +!!$ write(0,*) ' Context on remapping ',lme,lnp if ((lme >=0).and.(lnp>=2)) then - associate(lv=>prec%precv(i), rmp => prec%precv(i)%remap_data) + associate(lv=>prec%precv(level), rmp => prec%precv(level)%remap_data) call lv%desc_ac%clone(rmp%desc_ac_pre_remap,info) call lv%ac%clone(rmp%ac_pre_remap,info) - write(0,*) ' Doing remapping ',lnp, lnp/2 +!!$ write(0,*) 'During first remapping desc_ac:',lv%desc_ac%is_asb(),& +!!$ & rmp%desc_ac_pre_remap%is_asb() +!!$ write(0,*) ' First Doing remapping ',lnp, lnp/2 call psb_remap(lnp/2,rmp%desc_ac_pre_remap,rmp%ac_pre_remap,& & rmp%idest,rmp%isrc,rmp%nrsrc,rmp%naggr,lv%desc_ac,lv%ac,info) !!$ write(0,*) me,' Out of remapping ',rmp%desc_ac_pre_remap%get_fmt(),' ',& !!$ & lv%desc_ac%get_fmt(),sum(lv%linmap%naggr),sum(rmp%naggr) - write(0,*) 'Assignment ',size(lv%linmap%naggr),size(rmp%naggr) +!!$ write(0,*) 'First Assignment ',size(lv%linmap%naggr),size(rmp%naggr) lv%linmap%naggr(:) = rmp%naggr(:) lv%linmap%p_desc_V => rmp%desc_ac_pre_remap lv%base_a => lv%ac lv%base_desc => lv%desc_ac + block + integer(psb_ipk_) :: meu,npu,mev,npv + type(psb_ctxt_type) :: ct + ct = lv%linmap%p_desc_U%get_ctxt() + call psb_info(ct,meu,npu) + ct = lv%linmap%p_desc_V%get_ctxt() + call psb_info(ct,mev,npv) +!!$ write(0,*) 'First Check on out remapping ',i,& +!!$ & rmp%desc_ac_pre_remap%is_asb(),& +!!$ & ':',meu,npu,mev,npv + end block end associate end if - end block - end if - write(0,*) ' End of array_build_loop',i,iszv,info - end do array_build_loop - write(0,*) ' Done array_build_loop',iszv,newsz,info,psb_errstatus_fatal() - call psb_barrier(ctxt) - if (newsz > 0) then - ! - ! We exited early from the build loop, need to fix - ! the size. - ! - allocate(tprecv(newsz),stat=info) +!!$ write(0,*) 'Second Check on out remapping ',level,& +!!$ & prec%precv(level)%remap_data%desc_ac_pre_remap%is_asb(), newsz + end if + end block + if (info /= psb_success_) then - call psb_errpush(psb_err_from_subroutine_,name,& - & a_err='prec reallocation') + write(ch_err,'(a,i7)') 'Mat asb fail @ level ',i + call psb_errpush(psb_err_internal_error_,name,& + & a_err=ch_err) goto 9999 endif - do i=1,newsz - call prec%precv(i)%move_alloc(tprecv(i),info) - end do - do i=newsz+1, iszv - call prec%precv(i)%free(info) - end do - call move_alloc(tprecv,prec%precv) - ! Ignore errors from transfer - info = psb_success_ - ! - ! Restart - iszv = newsz - ! Fix the pointers, but the level 1 should - ! be treated differently - if (.not.associated(prec%precv(1)%base_a,a)) then - prec%precv(1)%base_a => prec%precv(1)%ac - end if - if (.not.associated(prec%precv(1)%base_desc,desc_a)) then - prec%precv(1)%base_desc => prec%precv(1)%desc_ac + if (stop_hierarchy_loop) then + exit array_build_loop + else + if (i prec%precv(i)%ac - prec%precv(i)%base_desc => prec%precv(i)%desc_ac - ! This is needed when the linmap object has been built - ! reusing the base_desc descriptor through a pointer. - ! With PSBLAS 4 we will have a better solution - if (associated(prec%precv(i)%linmap%p_desc_U)) & - & prec%precv(i)%linmap%p_desc_U => prec%precv(i-1)%base_desc - if (associated(prec%precv(i)%linmap%p_desc_V))& - & prec%precv(i)%linmap%p_desc_V => prec%precv(i)%base_desc - end do + end do array_build_loop + +!!$ write(0,*) ' Done array_build_loop',iszv,newsz,info,psb_errstatus_fatal() + + if (newsz>0) then +!!$ do i=2,newsz +!!$ write(0,*) me,'Newsz Out of array_build_loop ',i,':',& +!!$ & prec%precv(i)%remap_data%desc_ac_pre_remap%is_asb() +!!$ end do +!!$ write(0,*) 'Calling set_nlevs ',newsz + call prec%set_nlevs(newsz) + else +!!$ do i=2, iszv +!!$ write(0,*) me,'Out of array_build_loop ',i,':',& +!!$ & prec%precv(i)%remap_data%desc_ac_pre_remap%is_asb() +!!$ end do end if - write(0,*) ' Done reallocating precv',iszv,newsz,info,psb_errstatus_fatal() + iszv = prec%get_nlevs() call psb_barrier(ctxt) - !write(0,*) 'Should we remap? ' - if (.false.) then - if (amg_get_do_remap().and.(np>=4)) then -!!$ write(0,*) 'Going for remapping ' - if (.true.) then - associate(lv=>prec%precv(iszv), rmp => prec%precv(iszv)%remap_data) - call lv%desc_ac%clone(rmp%desc_ac_pre_remap,info) - call lv%ac%clone(rmp%ac_pre_remap,info) - if (np >= 8) then - call psb_remap(np/4,rmp%desc_ac_pre_remap,rmp%ac_pre_remap,& - & rmp%idest,rmp%isrc,rmp%nrsrc,rmp%naggr,lv%desc_ac,lv%ac,info) - else - call psb_remap(np/2,rmp%desc_ac_pre_remap,rmp%ac_pre_remap,& - & rmp%idest,rmp%isrc,rmp%nrsrc,rmp%naggr,lv%desc_ac,lv%ac,info) - end if -!!$ write(0,*) me,' Out of remapping ',rmp%desc_ac_pre_remap%get_fmt(),' ',& -!!$ & lv%desc_ac%get_fmt(),sum(lv%linmap%naggr),sum(rmp%naggr) - lv%linmap%naggr(:) = rmp%naggr(:) - lv%linmap%p_desc_V => rmp%desc_ac_pre_remap - lv%base_a => lv%ac - lv%base_desc => lv%desc_ac - end associate - end if - end if - end if +!!$ write(0,*) ' Done reallocating precv',iszv,newsz,info +!!$ +!!$ do i=2, iszv +!!$ write(0,*) me,'At end of hierarchy_bld level',i,':',& +!!$ & prec%precv(i)%remap_data%desc_ac_pre_remap%is_asb() +!!$ end do + call psb_barrier(ctxt) + if (info /= psb_success_) then call psb_errpush(psb_err_internal_error_,name,& & a_err='Internal hierarchy build' ) goto 9999 endif - iszv = size(prec%precv) - + iszv = prec%get_nlevs() +!!$ write(0,*) 'Going for cmp_complexity ',& +!!$ & allocated(prec%precv),iszv,size(prec%precv) call prec%cmp_complexity() call prec%cmp_avg_cr() @@ -730,4 +657,49 @@ contains return end subroutine restore_smoothers #endif + + function amg_c_policy_do_remap(ctxt,level,aggsize) result(res) + logical :: res + type(psb_ctxt_type) :: ctxt + integer(psb_ipk_) :: level + integer(psb_lpk_) :: aggsize + res = amg_get_do_remap().and.(level>=2) +!!$ res = .false. + end function amg_c_policy_do_remap + + subroutine amg_c_hierarchy_bld_cmp_newsz(level,iszv,prevsize,& + & nlaggr,casize,mnratio,sizeratio,newsz) + implicit none + integer(psb_ipk_) :: level,iszv,newsz + integer(psb_lpk_) :: nlaggr(:) + integer(psb_lpk_) :: prevsize, casize + real(psb_spk_) :: mnratio, sizeratio + ! ============================== + integer(psb_lpk_) :: iaggsize + + newsz = 0 + iaggsize = sum(nlaggr) + sizeratio = prevsize + sizeratio = sizeratio/iaggsize +!!$ write(0,*) 'From cmp_newsz: ',iaggsize,casize,& +!!$ & sizeratio,mnratio, level + + if (iaggsize <= casize) newsz = level + if (level == iszv) newsz = level + + if (level>2) then + if (sizeratio < mnratio) then + if (sizeratio > 1) then + newsz = level + else + ! + ! We are not gaining + ! + newsz = level-1 + end if + end if + end if +!!$ write(0,*) 'At end of cmp_newsz ',newsz + end subroutine amg_c_hierarchy_bld_cmp_newsz + end subroutine amg_c_hierarchy_bld diff --git a/amgprec/impl/amg_c_smoothers_bld.f90 b/amgprec/impl/amg_c_smoothers_bld.f90 index 3f425349..a21995ee 100644 --- a/amgprec/impl/amg_c_smoothers_bld.f90 +++ b/amgprec/impl/amg_c_smoothers_bld.f90 @@ -136,9 +136,9 @@ subroutine amg_c_smoothers_bld(a,desc_a,prec,info,amold,vmold,imold) ! ! Check to ensure all procs have the same ! - iszv = size(prec%precv) + iszv = prec%get_nlevs() call psb_bcast(ctxt,iszv) - if (iszv /= size(prec%precv)) then + if (iszv /= prec%get_nlevs()) then info=psb_err_internal_error_ call psb_errpush(info,name,a_err='Inconsistent size of precv') goto 9999 diff --git a/amgprec/impl/amg_cfile_prec_descr.f90 b/amgprec/impl/amg_cfile_prec_descr.f90 index 399e719b..d6ec343b 100644 --- a/amgprec/impl/amg_cfile_prec_descr.f90 +++ b/amgprec/impl/amg_cfile_prec_descr.f90 @@ -129,7 +129,7 @@ subroutine amg_cfile_prec_descr(prec,info,iout,root, verbosity,prefix) ! ensured by amg_precbld). ! if (me == root_) then - nlev = size(prec%precv) + nlev = prec%get_nlevs() do ilev = 1, nlev if (.not.allocated(prec%precv(ilev)%sm)) then info = 3111 @@ -141,6 +141,7 @@ subroutine amg_cfile_prec_descr(prec,info,iout,root, verbosity,prefix) write(iout_,*) write(iout_,'(a,1x,a)') trim(prefix_),'Preconditioner description' + write(iout_,*) 'At level :',1,' we have ',np,' processes' if (nlev == 1) then ! diff --git a/amgprec/impl/amg_cmlprec_aply.f90 b/amgprec/impl/amg_cmlprec_aply.f90 index c8dac94a..7d8f5df1 100644 --- a/amgprec/impl/amg_cmlprec_aply.f90 +++ b/amgprec/impl/amg_cmlprec_aply.f90 @@ -244,10 +244,10 @@ subroutine amg_cmlprec_aply_vect(alpha,p,x,beta,y,desc_data,trans,work,info) if (debug_level >= psb_debug_inner_) & & write(debug_unit,*) me,' ',trim(name),& - & ' Entry ', size(p%precv) + & ' Entry ', p%get_nlevs() trans_ = psb_toupper(trans) - nlev = size(p%precv) + nlev = p%get_nlevs() do_alloc_wrk = .not.allocated(p%precv(1)%wrk) @@ -382,7 +382,7 @@ contains call psb_erractionsave(err_act) debug_unit = psb_get_debug_unit() debug_level = psb_get_debug_level() - nlev = size(p%precv) + nlev = p%get_nlevs() if ((level < 1) .or. (level > nlev)) then call psb_errpush(psb_err_internal_error_,name,& & a_err='wrong call level to inner_ml') @@ -468,7 +468,7 @@ contains call psb_erractionsave(err_act) debug_unit = psb_get_debug_unit() debug_level = psb_get_debug_level() - nlev = size(p%precv) + nlev = p%get_nlevs() if ((level < 1) .or. (level > nlev)) then call psb_errpush(psb_err_internal_error_,name,& & a_err='wrong call level to inner_add') @@ -497,7 +497,8 @@ contains if (allocated(p%precv(level)%sm2a)) then call psb_geaxpby(cone,vx2l,czero,vy2l,base_desc,info) - sweeps = max(p%precv(level)%parms%sweeps_pre,p%precv(level)%parms%sweeps_post) + sweeps = max(p%precv(level)%parms%sweeps_pre,& + & p%precv(level)%parms%sweeps_post) do k=1, sweeps call p%precv(level)%sm%apply(cone,& & vy2l,czero,vty,& @@ -593,7 +594,7 @@ contains call psb_erractionsave(err_act) debug_unit = psb_get_debug_unit() debug_level = psb_get_debug_level() - nlev = size(p%precv) + nlev = p%get_nlevs() if ((level < 1) .or. (level > nlev)) then call psb_errpush(psb_err_internal_error_,name,& & a_err='wrong call level to inner_mult') @@ -604,7 +605,7 @@ contains if(debug_level > 1) then write(debug_unit,*) me,' inner_mult at level ',level end if - write(debug_unit,*) me,' inner_mult at level (1):',level,np +!!$ write(debug_unit,*) me,' inner_mult at level (1):',level,np sweeps_post = p%precv(level)%parms%sweeps_post sweeps_pre = p%precv(level)%parms%sweeps_pre pre = ((sweeps_pre>0).and.(trans=='N')).or.((sweeps_post>0).and.(trans/='N')) @@ -614,8 +615,8 @@ contains & vtx => p%precv(level)%wrk%vtx,vty => p%precv(level)%wrk%vty,& & base_a => p%precv(level)%base_a, base_desc=>p%precv(level)%base_desc,& & wv => p%precv(level)%wrk%wv) - write(0,*) 'Inner mult at level (2):',level,' :',me,np,':',& - & size(p%precv(level)%wrk%wv), allocated(p%precv(level)%wrk%wv) +!!$ write(0,*) 'Inner mult at level (2):',level,' :',me,np,':',& +!!$ & size(p%precv(level)%wrk%wv), allocated(p%precv(level)%wrk%wv) if (me >=0) then if (level < nlev) then @@ -831,7 +832,7 @@ contains call psb_erractionsave(err_act) debug_unit = psb_get_debug_unit() debug_level = psb_get_debug_level() - nlev = size(p%precv) + nlev = p%get_nlevs() if ((level < 1) .or. (level > nlev)) then call psb_errpush(psb_err_internal_error_,name,& & a_err='wrong call level to inner_add') diff --git a/amgprec/impl/amg_cprecinit.F90 b/amgprec/impl/amg_cprecinit.F90 index 6422aaac..68287b69 100644 --- a/amgprec/impl/amg_cprecinit.F90 +++ b/amgprec/impl/amg_cprecinit.F90 @@ -214,9 +214,7 @@ subroutine amg_cprecinit(ctxt,prec,ptype,info) allocate(amg_c_ilu_solver_type :: prec%precv(ilev_)%sm%sv, stat=info) call prec%precv(ilev_)%default() - case ('ML') - nlev_ = prec%ag_data%max_levs ilev_ = 1 allocate(prec%precv(nlev_),stat=info) @@ -228,6 +226,8 @@ subroutine amg_cprecinit(ctxt,prec,ptype,info) do ilev_ = 1, nlev_ call prec%precv(ilev_)%default() end do + call prec%set_nlevs(nlev_) + call prec%set('ML_CYCLE','VCYCLE',info) call prec%set('SMOOTHER_TYPE','FBGS',info) #if defined(AMG_HAVE_MUMPS) @@ -241,7 +241,6 @@ subroutine amg_cprecinit(ctxt,prec,ptype,info) write(psb_err_unit,*) name,& &': Warning: Unknown preconditioner type request "',ptype,'"' info = psb_err_pivot_too_small_ - end select call psb_erractionrestore(err_act) diff --git a/amgprec/impl/amg_d_extprol_bld.F90 b/amgprec/impl/amg_d_extprol_bld.F90 index e72e196d..68fa4726 100644 --- a/amgprec/impl/amg_d_extprol_bld.F90 +++ b/amgprec/impl/amg_d_extprol_bld.F90 @@ -164,7 +164,7 @@ subroutine amg_d_extprol_bld(a,desc_a,p,prolv,restrv,info,amold,vmold,imold) mxplevs = p%ag_data%max_levs mnaggratio = p%ag_data%min_cr_ratio casize = p%ag_data%min_coarse_size - iszv = p%get_nlevs() + iszv = size(p%precv) nprolv = size(prolv) nrestrv = size(restrv) call psb_bcast(ctxt,iszv) @@ -188,7 +188,7 @@ subroutine amg_d_extprol_bld(a,desc_a,p,prolv,restrv,info,amold,vmold,imold) call psb_errpush(info,name,a_err='Inconsistent min_cr_ratio') goto 9999 end if - if (iszv /= p%get_nlevs()) then + if (iszv /= size(p%precv)) then info=psb_err_internal_error_ call psb_errpush(info,name,a_err='Inconsistent size of precv') goto 9999 @@ -283,7 +283,7 @@ subroutine amg_d_extprol_bld(a,desc_a,p,prolv,restrv,info,amold,vmold,imold) call p%precv(i)%free(info) end do call move_alloc(tprecv,p%precv) - iszv = p%get_nlevs() + iszv = size(p%precv) end if ! ! Finest level first; remember to fix base_a and base_desc @@ -325,7 +325,7 @@ subroutine amg_d_extprol_bld(a,desc_a,p,prolv,restrv,info,amold,vmold,imold) goto 9999 endif - iszv = p%get_nlevs() + iszv = size(p%precv) if (debug_level >= psb_debug_outer_) & & write(debug_unit,*) me,' ',trim(name),& diff --git a/amgprec/impl/amg_d_hierarchy_bld.F90 b/amgprec/impl/amg_d_hierarchy_bld.F90 index e8540a51..47b6a510 100644 --- a/amgprec/impl/amg_d_hierarchy_bld.F90 +++ b/amgprec/impl/amg_d_hierarchy_bld.F90 @@ -64,11 +64,9 @@ ! Error code. ! subroutine amg_d_hierarchy_bld(a,desc_a,prec,info,cpymat) - use psb_base_mod use amg_d_inner_mod use amg_d_prec_mod, amg_protect_name => amg_d_hierarchy_bld - Implicit None ! Arguments @@ -306,7 +304,6 @@ subroutine amg_d_hierarchy_bld(a,desc_a,prec,info,cpymat) call psb_cd_renum_block(desc_a,prec%precv(1)%desc_ac,info) prec%precv(1)%base_desc => prec%precv(1)%desc_ac - ! ! Main build loop ! @@ -340,7 +337,7 @@ subroutine amg_d_hierarchy_bld(a,desc_a,prec,info,cpymat) & 'Calling mlprcbld at level ',i ! ! Build the tentative mapping between levels i-1 and i - ! and the matrixat level i + ! and the matrix at level i ! if (do_timings) call psb_tic(idx_bldtp) if (info == psb_success_)& @@ -512,7 +509,6 @@ subroutine amg_d_hierarchy_bld(a,desc_a,prec,info,cpymat) !!$ end do call psb_barrier(ctxt) - if (info /= psb_success_) then call psb_errpush(psb_err_internal_error_,name,& & a_err='Internal hierarchy build' ) diff --git a/amgprec/impl/amg_d_hierarchy_rebld.f90 b/amgprec/impl/amg_d_hierarchy_rebld.f90 index def6534e..885f95fa 100644 --- a/amgprec/impl/amg_d_hierarchy_rebld.f90 +++ b/amgprec/impl/amg_d_hierarchy_rebld.f90 @@ -122,7 +122,7 @@ subroutine amg_d_hierarchy_rebld(a,desc_a,prec,info) end if - iszv = prec%get_nlevs() + iszv = size(prec%precv) do i=2, iszv call prec%precv(i-1)%base_a%cp_to(acsr) diff --git a/amgprec/impl/amg_dcprecset.F90 b/amgprec/impl/amg_dcprecset.F90 index 18e6f324..350ddffd 100644 --- a/amgprec/impl/amg_dcprecset.F90 +++ b/amgprec/impl/amg_dcprecset.F90 @@ -129,7 +129,7 @@ subroutine amg_dcprecseti(p,what,val,info,ilev,ilmax,pos,idx) return endif - nlev_ = p%get_nlevs() + nlev_ = size(p%precv) if (present(ilev)) then ilev_ = ilev @@ -376,7 +376,7 @@ subroutine amg_dcprecsetc(p,what,string,info,ilev,ilmax,pos,idx) return endif - nlev_ = p%get_nlevs() + nlev_ = size(p%precv) if (present(ilev)) then ilev_ = ilev @@ -1057,7 +1057,7 @@ subroutine amg_dcprecsetr(p,what,val,info,ilev,ilmax,pos,idx) info = 3111 return endif - nlev_ = p%get_nlevs() + nlev_ = size(p%precv) if (present(ilev)) then ilev_ = ilev diff --git a/amgprec/impl/amg_dfile_prec_memory_use.f90 b/amgprec/impl/amg_dfile_prec_memory_use.f90 index 2b689cea..3757d6a4 100644 --- a/amgprec/impl/amg_dfile_prec_memory_use.f90 +++ b/amgprec/impl/amg_dfile_prec_memory_use.f90 @@ -135,7 +135,7 @@ subroutine amg_dfile_prec_memory_use(prec,info,iout,root, verbosity,prefix,globa write(iout_,*) write(iout_,'(a,1x,a)') trim(prefix_),'Preconditioner memory usage' end if - nlev = prec%get_nlevs() + nlev = size(prec%precv) do ilev=1,nlev call prec%precv(ilev)%memory_use(ilev,nlev,ilmin,info, & & iout=iout_,verbosity=verbosity_,prefix=trim(prefix_),global=global) diff --git a/amgprec/impl/amg_dmlprec_aply.f90 b/amgprec/impl/amg_dmlprec_aply.f90 index d0b3e0ad..2beda0a1 100644 --- a/amgprec/impl/amg_dmlprec_aply.f90 +++ b/amgprec/impl/amg_dmlprec_aply.f90 @@ -651,20 +651,19 @@ contains ! if (pre) then - call psb_geaxpby(done,vx2l,& - & dzero,vty,& - & base_desc,info) + call psb_geaxpby(done,vx2l,& + & dzero,vty,& + & base_desc,info) + + if (info == psb_success_) call psb_spmm(-done,base_a,& + & vy2l,done,vty,& + & base_desc,info,work=work,trans=trans) + if (info /= psb_success_) then + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Error during residue') + goto 9999 + end if - if (info == psb_success_) call psb_spmm(-done,base_a,& - & vy2l,done,vty,& - & base_desc,info,work=work,trans=trans) - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Error during residue') - goto 9999 - end if - -!!$ write(0,*) me,' Map_rstr from ', level,' to ',level+1 call p%precv(level+1)%map_rstr(done,vty,& & dzero,p%precv(level+1)%wrk%vx2l,& & info,work=work,vtx=wv(1)) diff --git a/amgprec/impl/amg_dmlprec_aply_a.f90 b/amgprec/impl/amg_dmlprec_aply_a.f90 index 5ec271c5..832656b5 100644 --- a/amgprec/impl/amg_dmlprec_aply_a.f90 +++ b/amgprec/impl/amg_dmlprec_aply_a.f90 @@ -249,11 +249,11 @@ subroutine amg_dmlprec_aply_a(alpha,p,x,beta,y,desc_data,trans,work,info) if (debug_level >= psb_debug_inner_) & & write(debug_unit,*) me,' ',trim(name),& - & ' Entry ', p%get_nlevs() + & ' Entry ', size(p%precv) trans_ = psb_toupper(trans) - nlev = p%get_nlevs() + nlev = size(p%precv) allocate(mlwrk(nlev),stat=info) if (info /= psb_success_) then call psb_errpush(psb_err_from_subroutine_,name,a_err='Allocate') @@ -362,7 +362,7 @@ contains call psb_erractionsave(err_act) debug_unit = psb_get_debug_unit() debug_level = psb_get_debug_level() - nlev = p%get_nlevs() + nlev = size(p%precv) if ((level < 1) .or. (level > nlev)) then call psb_errpush(psb_err_internal_error_,name,& & a_err='wrong call level to inner_ml') @@ -445,7 +445,7 @@ contains call psb_erractionsave(err_act) debug_unit = psb_get_debug_unit() debug_level = psb_get_debug_level() - nlev = p%get_nlevs() + nlev = size(p%precv) if ((level < 1) .or. (level > nlev)) then call psb_errpush(psb_err_internal_error_,name,& & a_err='wrong call level to inner_add') @@ -550,7 +550,7 @@ contains call psb_erractionsave(err_act) debug_unit = psb_get_debug_unit() debug_level = psb_get_debug_level() - nlev = p%get_nlevs() + nlev = size(p%precv) if ((level < 1) .or. (level > nlev)) then call psb_errpush(psb_err_internal_error_,name,& & a_err='wrong call level to inner_mult') diff --git a/amgprec/impl/amg_dprecinit.F90 b/amgprec/impl/amg_dprecinit.F90 index c8cf13da..e0d439e2 100644 --- a/amgprec/impl/amg_dprecinit.F90 +++ b/amgprec/impl/amg_dprecinit.F90 @@ -226,9 +226,7 @@ subroutine amg_dprecinit(ctxt,prec,ptype,info) allocate(amg_d_ilu_solver_type :: prec%precv(ilev_)%sm%sv, stat=info) call prec%precv(ilev_)%default() - case ('ML') - nlev_ = prec%ag_data%max_levs ilev_ = 1 allocate(prec%precv(nlev_),stat=info) diff --git a/amgprec/impl/amg_s_hierarchy_bld.F90 b/amgprec/impl/amg_s_hierarchy_bld.F90 index f54d4379..78b46dbf 100644 --- a/amgprec/impl/amg_s_hierarchy_bld.F90 +++ b/amgprec/impl/amg_s_hierarchy_bld.F90 @@ -64,11 +64,9 @@ ! Error code. ! subroutine amg_s_hierarchy_bld(a,desc_a,prec,info,cpymat) - use psb_base_mod use amg_s_inner_mod use amg_s_prec_mod, amg_protect_name => amg_s_hierarchy_bld - Implicit None ! Arguments @@ -82,7 +80,7 @@ subroutine amg_s_hierarchy_bld(a,desc_a,prec,info,cpymat) type(psb_ctxt_type) :: ctxt integer(psb_ipk_) :: me,np integer(psb_ipk_) :: err,i,k, err_act, iszv, newsz,& - & nplevs, mxplevs + & nplevs, mxplevs, level integer(psb_lpk_) :: iaggsize, casize, mncsize, mncszpp real(psb_spk_) :: mnaggratio, sizeratio, athresh, aomega class(amg_s_base_smoother_type), allocatable :: coarse_sm, med_sm, & @@ -98,6 +96,7 @@ subroutine amg_s_hierarchy_bld(a,desc_a,prec,info,cpymat) character(len=40) :: ch_err integer(psb_ipk_), save :: idx_bldtp=-1, idx_matasb=-1 logical, parameter :: do_timings=.false. + logical :: stop_hierarchy_loop type(psb_ctxt_type) :: lctxt integer(psb_ipk_) :: lme,lnp @@ -141,7 +140,7 @@ subroutine amg_s_hierarchy_bld(a,desc_a,prec,info,cpymat) mnaggratio = prec%ag_data%min_cr_ratio mncsize = prec%ag_data%min_coarse_size mncszpp = prec%ag_data%min_coarse_size_per_process - iszv = size(prec%precv) + iszv = prec%get_nlevs() call psb_bcast(ctxt,iszv) call psb_bcast(ctxt,mncsize) call psb_bcast(ctxt,mncszpp) @@ -167,7 +166,7 @@ subroutine amg_s_hierarchy_bld(a,desc_a,prec,info,cpymat) call psb_errpush(info,name,a_err='Inconsistent min_cr_ratio') goto 9999 end if - if (iszv /= size(prec%precv)) then + if (iszv /= prec%get_nlevs()) then info=psb_err_internal_error_ call psb_errpush(info,name,a_err='Inconsistent size of precv') goto 9999 @@ -182,6 +181,7 @@ subroutine amg_s_hierarchy_bld(a,desc_a,prec,info,cpymat) call psb_errpush(info,name,a_err=ch_err) goto 9999 endif + if (iszv == 1) then ! ! This is OK, since it may be called by the user even if there @@ -229,7 +229,6 @@ subroutine amg_s_hierarchy_bld(a,desc_a,prec,info,cpymat) casize = mncsize end if prec%ag_data%target_coarse_size = casize - nplevs = max(itwo,mxplevs) ! @@ -242,7 +241,7 @@ subroutine amg_s_hierarchy_bld(a,desc_a,prec,info,cpymat) goto 9999 end if ! - ! First set desired number of levels + ! First set desired number of levels if different from default. ! if (iszv /= nplevs) then allocate(tprecv(nplevs),stat=info) @@ -288,7 +287,8 @@ subroutine amg_s_hierarchy_bld(a,desc_a,prec,info,cpymat) call prec%precv(i)%free(info) end do call move_alloc(tprecv,prec%precv) - iszv = size(prec%precv) + call prec%set_nlevs(nplevs) + iszv = prec%get_nlevs() end if ! @@ -304,16 +304,23 @@ subroutine amg_s_hierarchy_bld(a,desc_a,prec,info,cpymat) call psb_cd_renum_block(desc_a,prec%precv(1)%desc_ac,info) prec%precv(1)%base_desc => prec%precv(1)%desc_ac + ! + ! Main build loop + ! newsz = 0 + stop_hierarchy_loop = .false. array_build_loop: do i=2, iszv ! ! Check on the iprcparm contents: they should be the same ! on all processes. ! call psb_bcast(ctxt,prec%precv(i)%parms) + ! + ! Get current context: might have performed remapping + ! lctxt = prec%precv(i-1)%base_desc%get_ctxt() call psb_info(lctxt,lme,lnp) - write(0,*) 'Check at level',i,lme,lnp +!!$ write(0,*) 'Check at level',i,lme,lnp ! ! Sanity checks on the parameters ! @@ -329,8 +336,8 @@ subroutine amg_s_hierarchy_bld(a,desc_a,prec,info,cpymat) & write(debug_unit,*) me,' ',trim(name),& & 'Calling mlprcbld at level ',i ! - ! Build the mapping between levels i-1 and i and the matrix - ! at level i + ! Build the tentative mapping between levels i-1 and i + ! and the matrix at level i ! if (do_timings) call psb_tic(idx_bldtp) if (info == psb_success_)& @@ -352,62 +359,26 @@ subroutine amg_s_hierarchy_bld(a,desc_a,prec,info,cpymat) ! Save op_prol just in case ! call op_prol%clone(prec%precv(i)%tprol,info) + ! ! Check for early termination of aggregation loop. - ! - iaggsize = sum(nlaggr) - - sizeratio = iaggsize - if (i==2) then - sizeratio = desc_a%get_global_rows()/sizeratio + ! + if (i == 2) then + call amg_s_hierarchy_bld_cmp_newsz(i,iszv,& + & desc_a%get_global_rows(),& + & nlaggr,casize,mnaggratio,sizeratio,newsz) else - sizeratio = sum(prec%precv(i-1)%linmap%naggr)/sizeratio + call amg_s_hierarchy_bld_cmp_newsz(i,iszv,& + & sum(prec%precv(i-1)%linmap%naggr),& + & nlaggr,casize,mnaggratio,sizeratio,newsz) end if prec%precv(i)%szratio = sizeratio - - if (iaggsize <= casize) newsz = i - if (i == iszv) newsz = i - - if (i>2) then - if (sizeratio < mnaggratio) then - if (sizeratio > 1) then - newsz = i - else - ! - ! We are not gaining - ! - newsz = i-1 - end if - end if - - ! - ! TO BE REVIWED - ! - write(0,*) me,lme,allocated(prec%precv(i-1)%linmap%naggr) - if (allocated(prec%precv(i-1)%linmap%naggr)) then - write(0,*) me,lme,size(prec%precv(i-1)%linmap%naggr), size(nlaggr) - end if - if (.false.) then - if (lme >=0) then - if (all(nlaggr == prec%precv(i-1)%linmap%naggr)) then - newsz=i-1 - if (me == 0) then - write(debug_unit,*) trim(name),& - &': Warning: aggregates from level ',& - & newsz - write(debug_unit,*) trim(name),& - &': to level ',& - & iszv,' coincide.' - write(debug_unit,*) trim(name),& - &': Number of levels actually used :',newsz - write(debug_unit,*) - end if - end if - end if - end if - end if call psb_bcast(ctxt,newsz) + ! + ! Handle reallocation, if needed, and then mat_asb to polish off the + ! construction + ! if (newsz > 0) then ! ! This is awkward, we are saving the aggregation parms, for the sake @@ -441,156 +412,112 @@ subroutine amg_s_hierarchy_bld(a,desc_a,prec,info,cpymat) & a_err=ch_err) goto 9999 endif - if (amg_get_do_remap().and.(i>=2)) then - block - type(psb_ctxt_type) :: lctxt - integer(psb_ipk_) :: lme,lnp - lctxt = prec%precv(i)%desc_ac%get_ctxt() - call psb_info(lctxt,lme,lnp) - write(0,*) ' Context on remapping ',lme,lnp - if ((lme >=0).and.(lnp>=2)) then - associate(lv=>prec%precv(i), rmp => prec%precv(i)%remap_data) - call lv%desc_ac%clone(rmp%desc_ac_pre_remap,info) - call lv%ac%clone(rmp%ac_pre_remap,info) - write(0,*) ' Doing remapping ',lnp, lnp/2 - call psb_remap(lnp/2,rmp%desc_ac_pre_remap,rmp%ac_pre_remap,& - & rmp%idest,rmp%isrc,rmp%nrsrc,rmp%naggr,lv%desc_ac,lv%ac,info) -!!$ write(0,*) me,' Out of remapping ',rmp%desc_ac_pre_remap%get_fmt(),' ',& -!!$ & lv%desc_ac%get_fmt(),sum(lv%linmap%naggr),sum(rmp%naggr) - write(0,*) 'Assignment ',size(lv%linmap%naggr),size(rmp%naggr) - lv%linmap%naggr(:) = rmp%naggr(:) - lv%linmap%p_desc_V => rmp%desc_ac_pre_remap - lv%base_a => lv%ac - lv%base_desc => lv%desc_ac - end associate - end if - end block - end if - exit array_build_loop +!!$ write(0,*) ' Early exit of array_build_loop',i,iszv,info,& + level = newsz + stop_hierarchy_loop = .true. else - if (do_timings) call psb_tic(idx_matasb) + if (do_timings) call psb_tic(idx_matasb) if (info == psb_success_) call prec%precv(i)%mat_asb(& & prec%precv(i-1)%base_a,prec%precv(i-1)%base_desc,& & ilaggr,nlaggr,op_prol,info) if (do_timings) call psb_toc(idx_matasb) + level = i end if - if (info /= psb_success_) then - write(ch_err,'(a,i7)') 'Mat asb fail @ level ',i - call psb_errpush(psb_err_internal_error_,name,& - & a_err=ch_err) - goto 9999 - endif - if (i=2)) then - block - type(psb_ctxt_type) :: lctxt - integer(psb_ipk_) :: lme,lnp - lctxt = prec%precv(i)%desc_ac%get_ctxt() - call psb_info(lctxt,lme,lnp) - write(0,*) ' Context on remapping ',lme,lnp + + ! + ! Do we want to remap onto a smaller subset of processes? + ! Will need a more sophisticated policy + ! + block + type(psb_ctxt_type) :: lctxt + integer(psb_ipk_) :: lme,lnp + lctxt = prec%precv(level)%desc_ac%get_ctxt() + call psb_info(lctxt,lme,lnp) + if (amg_s_policy_do_remap(lctxt,level,sum(nlaggr))) then +!!$ write(0,*) ' Context on remapping ',lme,lnp if ((lme >=0).and.(lnp>=2)) then - associate(lv=>prec%precv(i), rmp => prec%precv(i)%remap_data) + associate(lv=>prec%precv(level), rmp => prec%precv(level)%remap_data) call lv%desc_ac%clone(rmp%desc_ac_pre_remap,info) call lv%ac%clone(rmp%ac_pre_remap,info) - write(0,*) ' Doing remapping ',lnp, lnp/2 +!!$ write(0,*) 'During first remapping desc_ac:',lv%desc_ac%is_asb(),& +!!$ & rmp%desc_ac_pre_remap%is_asb() +!!$ write(0,*) ' First Doing remapping ',lnp, lnp/2 call psb_remap(lnp/2,rmp%desc_ac_pre_remap,rmp%ac_pre_remap,& & rmp%idest,rmp%isrc,rmp%nrsrc,rmp%naggr,lv%desc_ac,lv%ac,info) !!$ write(0,*) me,' Out of remapping ',rmp%desc_ac_pre_remap%get_fmt(),' ',& !!$ & lv%desc_ac%get_fmt(),sum(lv%linmap%naggr),sum(rmp%naggr) - write(0,*) 'Assignment ',size(lv%linmap%naggr),size(rmp%naggr) +!!$ write(0,*) 'First Assignment ',size(lv%linmap%naggr),size(rmp%naggr) lv%linmap%naggr(:) = rmp%naggr(:) lv%linmap%p_desc_V => rmp%desc_ac_pre_remap lv%base_a => lv%ac lv%base_desc => lv%desc_ac + block + integer(psb_ipk_) :: meu,npu,mev,npv + type(psb_ctxt_type) :: ct + ct = lv%linmap%p_desc_U%get_ctxt() + call psb_info(ct,meu,npu) + ct = lv%linmap%p_desc_V%get_ctxt() + call psb_info(ct,mev,npv) +!!$ write(0,*) 'First Check on out remapping ',i,& +!!$ & rmp%desc_ac_pre_remap%is_asb(),& +!!$ & ':',meu,npu,mev,npv + end block end associate end if - end block - end if - write(0,*) ' End of array_build_loop',i,iszv,info - end do array_build_loop - write(0,*) ' Done array_build_loop',iszv,newsz,info,psb_errstatus_fatal() - call psb_barrier(ctxt) - if (newsz > 0) then - ! - ! We exited early from the build loop, need to fix - ! the size. - ! - allocate(tprecv(newsz),stat=info) +!!$ write(0,*) 'Second Check on out remapping ',level,& +!!$ & prec%precv(level)%remap_data%desc_ac_pre_remap%is_asb(), newsz + end if + end block + if (info /= psb_success_) then - call psb_errpush(psb_err_from_subroutine_,name,& - & a_err='prec reallocation') + write(ch_err,'(a,i7)') 'Mat asb fail @ level ',i + call psb_errpush(psb_err_internal_error_,name,& + & a_err=ch_err) goto 9999 endif - do i=1,newsz - call prec%precv(i)%move_alloc(tprecv(i),info) - end do - do i=newsz+1, iszv - call prec%precv(i)%free(info) - end do - call move_alloc(tprecv,prec%precv) - ! Ignore errors from transfer - info = psb_success_ - ! - ! Restart - iszv = newsz - ! Fix the pointers, but the level 1 should - ! be treated differently - if (.not.associated(prec%precv(1)%base_a,a)) then - prec%precv(1)%base_a => prec%precv(1)%ac - end if - if (.not.associated(prec%precv(1)%base_desc,desc_a)) then - prec%precv(1)%base_desc => prec%precv(1)%desc_ac + if (stop_hierarchy_loop) then + exit array_build_loop + else + if (i prec%precv(i)%ac - prec%precv(i)%base_desc => prec%precv(i)%desc_ac - ! This is needed when the linmap object has been built - ! reusing the base_desc descriptor through a pointer. - ! With PSBLAS 4 we will have a better solution - if (associated(prec%precv(i)%linmap%p_desc_U)) & - & prec%precv(i)%linmap%p_desc_U => prec%precv(i-1)%base_desc - if (associated(prec%precv(i)%linmap%p_desc_V))& - & prec%precv(i)%linmap%p_desc_V => prec%precv(i)%base_desc - end do + end do array_build_loop + +!!$ write(0,*) ' Done array_build_loop',iszv,newsz,info,psb_errstatus_fatal() + + if (newsz>0) then +!!$ do i=2,newsz +!!$ write(0,*) me,'Newsz Out of array_build_loop ',i,':',& +!!$ & prec%precv(i)%remap_data%desc_ac_pre_remap%is_asb() +!!$ end do +!!$ write(0,*) 'Calling set_nlevs ',newsz + call prec%set_nlevs(newsz) + else +!!$ do i=2, iszv +!!$ write(0,*) me,'Out of array_build_loop ',i,':',& +!!$ & prec%precv(i)%remap_data%desc_ac_pre_remap%is_asb() +!!$ end do end if - write(0,*) ' Done reallocating precv',iszv,newsz,info,psb_errstatus_fatal() + iszv = prec%get_nlevs() call psb_barrier(ctxt) - !write(0,*) 'Should we remap? ' - if (.false.) then - if (amg_get_do_remap().and.(np>=4)) then -!!$ write(0,*) 'Going for remapping ' - if (.true.) then - associate(lv=>prec%precv(iszv), rmp => prec%precv(iszv)%remap_data) - call lv%desc_ac%clone(rmp%desc_ac_pre_remap,info) - call lv%ac%clone(rmp%ac_pre_remap,info) - if (np >= 8) then - call psb_remap(np/4,rmp%desc_ac_pre_remap,rmp%ac_pre_remap,& - & rmp%idest,rmp%isrc,rmp%nrsrc,rmp%naggr,lv%desc_ac,lv%ac,info) - else - call psb_remap(np/2,rmp%desc_ac_pre_remap,rmp%ac_pre_remap,& - & rmp%idest,rmp%isrc,rmp%nrsrc,rmp%naggr,lv%desc_ac,lv%ac,info) - end if -!!$ write(0,*) me,' Out of remapping ',rmp%desc_ac_pre_remap%get_fmt(),' ',& -!!$ & lv%desc_ac%get_fmt(),sum(lv%linmap%naggr),sum(rmp%naggr) - lv%linmap%naggr(:) = rmp%naggr(:) - lv%linmap%p_desc_V => rmp%desc_ac_pre_remap - lv%base_a => lv%ac - lv%base_desc => lv%desc_ac - end associate - end if - end if - end if +!!$ write(0,*) ' Done reallocating precv',iszv,newsz,info +!!$ +!!$ do i=2, iszv +!!$ write(0,*) me,'At end of hierarchy_bld level',i,':',& +!!$ & prec%precv(i)%remap_data%desc_ac_pre_remap%is_asb() +!!$ end do + call psb_barrier(ctxt) + if (info /= psb_success_) then call psb_errpush(psb_err_internal_error_,name,& & a_err='Internal hierarchy build' ) goto 9999 endif - iszv = size(prec%precv) - + iszv = prec%get_nlevs() +!!$ write(0,*) 'Going for cmp_complexity ',& +!!$ & allocated(prec%precv),iszv,size(prec%precv) call prec%cmp_complexity() call prec%cmp_avg_cr() @@ -730,4 +657,49 @@ contains return end subroutine restore_smoothers #endif + + function amg_s_policy_do_remap(ctxt,level,aggsize) result(res) + logical :: res + type(psb_ctxt_type) :: ctxt + integer(psb_ipk_) :: level + integer(psb_lpk_) :: aggsize + res = amg_get_do_remap().and.(level>=2) +!!$ res = .false. + end function amg_s_policy_do_remap + + subroutine amg_s_hierarchy_bld_cmp_newsz(level,iszv,prevsize,& + & nlaggr,casize,mnratio,sizeratio,newsz) + implicit none + integer(psb_ipk_) :: level,iszv,newsz + integer(psb_lpk_) :: nlaggr(:) + integer(psb_lpk_) :: prevsize, casize + real(psb_spk_) :: mnratio, sizeratio + ! ============================== + integer(psb_lpk_) :: iaggsize + + newsz = 0 + iaggsize = sum(nlaggr) + sizeratio = prevsize + sizeratio = sizeratio/iaggsize +!!$ write(0,*) 'From cmp_newsz: ',iaggsize,casize,& +!!$ & sizeratio,mnratio, level + + if (iaggsize <= casize) newsz = level + if (level == iszv) newsz = level + + if (level>2) then + if (sizeratio < mnratio) then + if (sizeratio > 1) then + newsz = level + else + ! + ! We are not gaining + ! + newsz = level-1 + end if + end if + end if +!!$ write(0,*) 'At end of cmp_newsz ',newsz + end subroutine amg_s_hierarchy_bld_cmp_newsz + end subroutine amg_s_hierarchy_bld diff --git a/amgprec/impl/amg_s_smoothers_bld.f90 b/amgprec/impl/amg_s_smoothers_bld.f90 index e274e68a..2f91de5f 100644 --- a/amgprec/impl/amg_s_smoothers_bld.f90 +++ b/amgprec/impl/amg_s_smoothers_bld.f90 @@ -136,9 +136,9 @@ subroutine amg_s_smoothers_bld(a,desc_a,prec,info,amold,vmold,imold) ! ! Check to ensure all procs have the same ! - iszv = size(prec%precv) + iszv = prec%get_nlevs() call psb_bcast(ctxt,iszv) - if (iszv /= size(prec%precv)) then + if (iszv /= prec%get_nlevs()) then info=psb_err_internal_error_ call psb_errpush(info,name,a_err='Inconsistent size of precv') goto 9999 diff --git a/amgprec/impl/amg_sfile_prec_descr.f90 b/amgprec/impl/amg_sfile_prec_descr.f90 index 08a4032e..d005fff1 100644 --- a/amgprec/impl/amg_sfile_prec_descr.f90 +++ b/amgprec/impl/amg_sfile_prec_descr.f90 @@ -129,7 +129,7 @@ subroutine amg_sfile_prec_descr(prec,info,iout,root, verbosity,prefix) ! ensured by amg_precbld). ! if (me == root_) then - nlev = size(prec%precv) + nlev = prec%get_nlevs() do ilev = 1, nlev if (.not.allocated(prec%precv(ilev)%sm)) then info = 3111 @@ -141,6 +141,7 @@ subroutine amg_sfile_prec_descr(prec,info,iout,root, verbosity,prefix) write(iout_,*) write(iout_,'(a,1x,a)') trim(prefix_),'Preconditioner description' + write(iout_,*) 'At level :',1,' we have ',np,' processes' if (nlev == 1) then ! diff --git a/amgprec/impl/amg_smlprec_aply.f90 b/amgprec/impl/amg_smlprec_aply.f90 index 439657bf..f7e48267 100644 --- a/amgprec/impl/amg_smlprec_aply.f90 +++ b/amgprec/impl/amg_smlprec_aply.f90 @@ -244,10 +244,10 @@ subroutine amg_smlprec_aply_vect(alpha,p,x,beta,y,desc_data,trans,work,info) if (debug_level >= psb_debug_inner_) & & write(debug_unit,*) me,' ',trim(name),& - & ' Entry ', size(p%precv) + & ' Entry ', p%get_nlevs() trans_ = psb_toupper(trans) - nlev = size(p%precv) + nlev = p%get_nlevs() do_alloc_wrk = .not.allocated(p%precv(1)%wrk) @@ -382,7 +382,7 @@ contains call psb_erractionsave(err_act) debug_unit = psb_get_debug_unit() debug_level = psb_get_debug_level() - nlev = size(p%precv) + nlev = p%get_nlevs() if ((level < 1) .or. (level > nlev)) then call psb_errpush(psb_err_internal_error_,name,& & a_err='wrong call level to inner_ml') @@ -468,7 +468,7 @@ contains call psb_erractionsave(err_act) debug_unit = psb_get_debug_unit() debug_level = psb_get_debug_level() - nlev = size(p%precv) + nlev = p%get_nlevs() if ((level < 1) .or. (level > nlev)) then call psb_errpush(psb_err_internal_error_,name,& & a_err='wrong call level to inner_add') @@ -497,7 +497,8 @@ contains if (allocated(p%precv(level)%sm2a)) then call psb_geaxpby(sone,vx2l,szero,vy2l,base_desc,info) - sweeps = max(p%precv(level)%parms%sweeps_pre,p%precv(level)%parms%sweeps_post) + sweeps = max(p%precv(level)%parms%sweeps_pre,& + & p%precv(level)%parms%sweeps_post) do k=1, sweeps call p%precv(level)%sm%apply(sone,& & vy2l,szero,vty,& @@ -593,7 +594,7 @@ contains call psb_erractionsave(err_act) debug_unit = psb_get_debug_unit() debug_level = psb_get_debug_level() - nlev = size(p%precv) + nlev = p%get_nlevs() if ((level < 1) .or. (level > nlev)) then call psb_errpush(psb_err_internal_error_,name,& & a_err='wrong call level to inner_mult') @@ -604,7 +605,7 @@ contains if(debug_level > 1) then write(debug_unit,*) me,' inner_mult at level ',level end if - write(debug_unit,*) me,' inner_mult at level (1):',level,np +!!$ write(debug_unit,*) me,' inner_mult at level (1):',level,np sweeps_post = p%precv(level)%parms%sweeps_post sweeps_pre = p%precv(level)%parms%sweeps_pre pre = ((sweeps_pre>0).and.(trans=='N')).or.((sweeps_post>0).and.(trans/='N')) @@ -614,8 +615,8 @@ contains & vtx => p%precv(level)%wrk%vtx,vty => p%precv(level)%wrk%vty,& & base_a => p%precv(level)%base_a, base_desc=>p%precv(level)%base_desc,& & wv => p%precv(level)%wrk%wv) - write(0,*) 'Inner mult at level (2):',level,' :',me,np,':',& - & size(p%precv(level)%wrk%wv), allocated(p%precv(level)%wrk%wv) +!!$ write(0,*) 'Inner mult at level (2):',level,' :',me,np,':',& +!!$ & size(p%precv(level)%wrk%wv), allocated(p%precv(level)%wrk%wv) if (me >=0) then if (level < nlev) then @@ -831,7 +832,7 @@ contains call psb_erractionsave(err_act) debug_unit = psb_get_debug_unit() debug_level = psb_get_debug_level() - nlev = size(p%precv) + nlev = p%get_nlevs() if ((level < 1) .or. (level > nlev)) then call psb_errpush(psb_err_internal_error_,name,& & a_err='wrong call level to inner_add') diff --git a/amgprec/impl/amg_sprecinit.F90 b/amgprec/impl/amg_sprecinit.F90 index cd67b094..7b577e3e 100644 --- a/amgprec/impl/amg_sprecinit.F90 +++ b/amgprec/impl/amg_sprecinit.F90 @@ -223,9 +223,7 @@ subroutine amg_sprecinit(ctxt,prec,ptype,info) allocate(amg_s_ilu_solver_type :: prec%precv(ilev_)%sm%sv, stat=info) call prec%precv(ilev_)%default() - case ('ML') - nlev_ = prec%ag_data%max_levs ilev_ = 1 allocate(prec%precv(nlev_),stat=info) @@ -237,6 +235,8 @@ subroutine amg_sprecinit(ctxt,prec,ptype,info) do ilev_ = 1, nlev_ call prec%precv(ilev_)%default() end do + call prec%set_nlevs(nlev_) + call prec%set('ML_CYCLE','VCYCLE',info) call prec%set('SMOOTHER_TYPE','FBGS',info) #if defined(AMG_HAVE_MUMPS) @@ -250,7 +250,6 @@ subroutine amg_sprecinit(ctxt,prec,ptype,info) write(psb_err_unit,*) name,& &': Warning: Unknown preconditioner type request "',ptype,'"' info = psb_err_pivot_too_small_ - end select call psb_erractionrestore(err_act) diff --git a/amgprec/impl/amg_z_hierarchy_bld.F90 b/amgprec/impl/amg_z_hierarchy_bld.F90 index 6a2cd661..223e9422 100644 --- a/amgprec/impl/amg_z_hierarchy_bld.F90 +++ b/amgprec/impl/amg_z_hierarchy_bld.F90 @@ -64,11 +64,9 @@ ! Error code. ! subroutine amg_z_hierarchy_bld(a,desc_a,prec,info,cpymat) - use psb_base_mod use amg_z_inner_mod use amg_z_prec_mod, amg_protect_name => amg_z_hierarchy_bld - Implicit None ! Arguments @@ -82,7 +80,7 @@ subroutine amg_z_hierarchy_bld(a,desc_a,prec,info,cpymat) type(psb_ctxt_type) :: ctxt integer(psb_ipk_) :: me,np integer(psb_ipk_) :: err,i,k, err_act, iszv, newsz,& - & nplevs, mxplevs + & nplevs, mxplevs, level integer(psb_lpk_) :: iaggsize, casize, mncsize, mncszpp real(psb_dpk_) :: mnaggratio, sizeratio, athresh, aomega class(amg_z_base_smoother_type), allocatable :: coarse_sm, med_sm, & @@ -98,6 +96,7 @@ subroutine amg_z_hierarchy_bld(a,desc_a,prec,info,cpymat) character(len=40) :: ch_err integer(psb_ipk_), save :: idx_bldtp=-1, idx_matasb=-1 logical, parameter :: do_timings=.false. + logical :: stop_hierarchy_loop type(psb_ctxt_type) :: lctxt integer(psb_ipk_) :: lme,lnp @@ -141,7 +140,7 @@ subroutine amg_z_hierarchy_bld(a,desc_a,prec,info,cpymat) mnaggratio = prec%ag_data%min_cr_ratio mncsize = prec%ag_data%min_coarse_size mncszpp = prec%ag_data%min_coarse_size_per_process - iszv = size(prec%precv) + iszv = prec%get_nlevs() call psb_bcast(ctxt,iszv) call psb_bcast(ctxt,mncsize) call psb_bcast(ctxt,mncszpp) @@ -167,7 +166,7 @@ subroutine amg_z_hierarchy_bld(a,desc_a,prec,info,cpymat) call psb_errpush(info,name,a_err='Inconsistent min_cr_ratio') goto 9999 end if - if (iszv /= size(prec%precv)) then + if (iszv /= prec%get_nlevs()) then info=psb_err_internal_error_ call psb_errpush(info,name,a_err='Inconsistent size of precv') goto 9999 @@ -182,6 +181,7 @@ subroutine amg_z_hierarchy_bld(a,desc_a,prec,info,cpymat) call psb_errpush(info,name,a_err=ch_err) goto 9999 endif + if (iszv == 1) then ! ! This is OK, since it may be called by the user even if there @@ -229,7 +229,6 @@ subroutine amg_z_hierarchy_bld(a,desc_a,prec,info,cpymat) casize = mncsize end if prec%ag_data%target_coarse_size = casize - nplevs = max(itwo,mxplevs) ! @@ -242,7 +241,7 @@ subroutine amg_z_hierarchy_bld(a,desc_a,prec,info,cpymat) goto 9999 end if ! - ! First set desired number of levels + ! First set desired number of levels if different from default. ! if (iszv /= nplevs) then allocate(tprecv(nplevs),stat=info) @@ -288,7 +287,8 @@ subroutine amg_z_hierarchy_bld(a,desc_a,prec,info,cpymat) call prec%precv(i)%free(info) end do call move_alloc(tprecv,prec%precv) - iszv = size(prec%precv) + call prec%set_nlevs(nplevs) + iszv = prec%get_nlevs() end if ! @@ -304,16 +304,23 @@ subroutine amg_z_hierarchy_bld(a,desc_a,prec,info,cpymat) call psb_cd_renum_block(desc_a,prec%precv(1)%desc_ac,info) prec%precv(1)%base_desc => prec%precv(1)%desc_ac + ! + ! Main build loop + ! newsz = 0 + stop_hierarchy_loop = .false. array_build_loop: do i=2, iszv ! ! Check on the iprcparm contents: they should be the same ! on all processes. ! call psb_bcast(ctxt,prec%precv(i)%parms) + ! + ! Get current context: might have performed remapping + ! lctxt = prec%precv(i-1)%base_desc%get_ctxt() call psb_info(lctxt,lme,lnp) - write(0,*) 'Check at level',i,lme,lnp +!!$ write(0,*) 'Check at level',i,lme,lnp ! ! Sanity checks on the parameters ! @@ -329,8 +336,8 @@ subroutine amg_z_hierarchy_bld(a,desc_a,prec,info,cpymat) & write(debug_unit,*) me,' ',trim(name),& & 'Calling mlprcbld at level ',i ! - ! Build the mapping between levels i-1 and i and the matrix - ! at level i + ! Build the tentative mapping between levels i-1 and i + ! and the matrix at level i ! if (do_timings) call psb_tic(idx_bldtp) if (info == psb_success_)& @@ -352,62 +359,26 @@ subroutine amg_z_hierarchy_bld(a,desc_a,prec,info,cpymat) ! Save op_prol just in case ! call op_prol%clone(prec%precv(i)%tprol,info) + ! ! Check for early termination of aggregation loop. - ! - iaggsize = sum(nlaggr) - - sizeratio = iaggsize - if (i==2) then - sizeratio = desc_a%get_global_rows()/sizeratio + ! + if (i == 2) then + call amg_z_hierarchy_bld_cmp_newsz(i,iszv,& + & desc_a%get_global_rows(),& + & nlaggr,casize,mnaggratio,sizeratio,newsz) else - sizeratio = sum(prec%precv(i-1)%linmap%naggr)/sizeratio + call amg_z_hierarchy_bld_cmp_newsz(i,iszv,& + & sum(prec%precv(i-1)%linmap%naggr),& + & nlaggr,casize,mnaggratio,sizeratio,newsz) end if prec%precv(i)%szratio = sizeratio - - if (iaggsize <= casize) newsz = i - if (i == iszv) newsz = i - - if (i>2) then - if (sizeratio < mnaggratio) then - if (sizeratio > 1) then - newsz = i - else - ! - ! We are not gaining - ! - newsz = i-1 - end if - end if - - ! - ! TO BE REVIWED - ! - write(0,*) me,lme,allocated(prec%precv(i-1)%linmap%naggr) - if (allocated(prec%precv(i-1)%linmap%naggr)) then - write(0,*) me,lme,size(prec%precv(i-1)%linmap%naggr), size(nlaggr) - end if - if (.false.) then - if (lme >=0) then - if (all(nlaggr == prec%precv(i-1)%linmap%naggr)) then - newsz=i-1 - if (me == 0) then - write(debug_unit,*) trim(name),& - &': Warning: aggregates from level ',& - & newsz - write(debug_unit,*) trim(name),& - &': to level ',& - & iszv,' coincide.' - write(debug_unit,*) trim(name),& - &': Number of levels actually used :',newsz - write(debug_unit,*) - end if - end if - end if - end if - end if call psb_bcast(ctxt,newsz) + ! + ! Handle reallocation, if needed, and then mat_asb to polish off the + ! construction + ! if (newsz > 0) then ! ! This is awkward, we are saving the aggregation parms, for the sake @@ -441,156 +412,112 @@ subroutine amg_z_hierarchy_bld(a,desc_a,prec,info,cpymat) & a_err=ch_err) goto 9999 endif - if (amg_get_do_remap().and.(i>=2)) then - block - type(psb_ctxt_type) :: lctxt - integer(psb_ipk_) :: lme,lnp - lctxt = prec%precv(i)%desc_ac%get_ctxt() - call psb_info(lctxt,lme,lnp) - write(0,*) ' Context on remapping ',lme,lnp - if ((lme >=0).and.(lnp>=2)) then - associate(lv=>prec%precv(i), rmp => prec%precv(i)%remap_data) - call lv%desc_ac%clone(rmp%desc_ac_pre_remap,info) - call lv%ac%clone(rmp%ac_pre_remap,info) - write(0,*) ' Doing remapping ',lnp, lnp/2 - call psb_remap(lnp/2,rmp%desc_ac_pre_remap,rmp%ac_pre_remap,& - & rmp%idest,rmp%isrc,rmp%nrsrc,rmp%naggr,lv%desc_ac,lv%ac,info) -!!$ write(0,*) me,' Out of remapping ',rmp%desc_ac_pre_remap%get_fmt(),' ',& -!!$ & lv%desc_ac%get_fmt(),sum(lv%linmap%naggr),sum(rmp%naggr) - write(0,*) 'Assignment ',size(lv%linmap%naggr),size(rmp%naggr) - lv%linmap%naggr(:) = rmp%naggr(:) - lv%linmap%p_desc_V => rmp%desc_ac_pre_remap - lv%base_a => lv%ac - lv%base_desc => lv%desc_ac - end associate - end if - end block - end if - exit array_build_loop +!!$ write(0,*) ' Early exit of array_build_loop',i,iszv,info,& + level = newsz + stop_hierarchy_loop = .true. else - if (do_timings) call psb_tic(idx_matasb) + if (do_timings) call psb_tic(idx_matasb) if (info == psb_success_) call prec%precv(i)%mat_asb(& & prec%precv(i-1)%base_a,prec%precv(i-1)%base_desc,& & ilaggr,nlaggr,op_prol,info) if (do_timings) call psb_toc(idx_matasb) + level = i end if - if (info /= psb_success_) then - write(ch_err,'(a,i7)') 'Mat asb fail @ level ',i - call psb_errpush(psb_err_internal_error_,name,& - & a_err=ch_err) - goto 9999 - endif - if (i=2)) then - block - type(psb_ctxt_type) :: lctxt - integer(psb_ipk_) :: lme,lnp - lctxt = prec%precv(i)%desc_ac%get_ctxt() - call psb_info(lctxt,lme,lnp) - write(0,*) ' Context on remapping ',lme,lnp + + ! + ! Do we want to remap onto a smaller subset of processes? + ! Will need a more sophisticated policy + ! + block + type(psb_ctxt_type) :: lctxt + integer(psb_ipk_) :: lme,lnp + lctxt = prec%precv(level)%desc_ac%get_ctxt() + call psb_info(lctxt,lme,lnp) + if (amg_z_policy_do_remap(lctxt,level,sum(nlaggr))) then +!!$ write(0,*) ' Context on remapping ',lme,lnp if ((lme >=0).and.(lnp>=2)) then - associate(lv=>prec%precv(i), rmp => prec%precv(i)%remap_data) + associate(lv=>prec%precv(level), rmp => prec%precv(level)%remap_data) call lv%desc_ac%clone(rmp%desc_ac_pre_remap,info) call lv%ac%clone(rmp%ac_pre_remap,info) - write(0,*) ' Doing remapping ',lnp, lnp/2 +!!$ write(0,*) 'During first remapping desc_ac:',lv%desc_ac%is_asb(),& +!!$ & rmp%desc_ac_pre_remap%is_asb() +!!$ write(0,*) ' First Doing remapping ',lnp, lnp/2 call psb_remap(lnp/2,rmp%desc_ac_pre_remap,rmp%ac_pre_remap,& & rmp%idest,rmp%isrc,rmp%nrsrc,rmp%naggr,lv%desc_ac,lv%ac,info) !!$ write(0,*) me,' Out of remapping ',rmp%desc_ac_pre_remap%get_fmt(),' ',& !!$ & lv%desc_ac%get_fmt(),sum(lv%linmap%naggr),sum(rmp%naggr) - write(0,*) 'Assignment ',size(lv%linmap%naggr),size(rmp%naggr) +!!$ write(0,*) 'First Assignment ',size(lv%linmap%naggr),size(rmp%naggr) lv%linmap%naggr(:) = rmp%naggr(:) lv%linmap%p_desc_V => rmp%desc_ac_pre_remap lv%base_a => lv%ac lv%base_desc => lv%desc_ac + block + integer(psb_ipk_) :: meu,npu,mev,npv + type(psb_ctxt_type) :: ct + ct = lv%linmap%p_desc_U%get_ctxt() + call psb_info(ct,meu,npu) + ct = lv%linmap%p_desc_V%get_ctxt() + call psb_info(ct,mev,npv) +!!$ write(0,*) 'First Check on out remapping ',i,& +!!$ & rmp%desc_ac_pre_remap%is_asb(),& +!!$ & ':',meu,npu,mev,npv + end block end associate end if - end block - end if - write(0,*) ' End of array_build_loop',i,iszv,info - end do array_build_loop - write(0,*) ' Done array_build_loop',iszv,newsz,info,psb_errstatus_fatal() - call psb_barrier(ctxt) - if (newsz > 0) then - ! - ! We exited early from the build loop, need to fix - ! the size. - ! - allocate(tprecv(newsz),stat=info) +!!$ write(0,*) 'Second Check on out remapping ',level,& +!!$ & prec%precv(level)%remap_data%desc_ac_pre_remap%is_asb(), newsz + end if + end block + if (info /= psb_success_) then - call psb_errpush(psb_err_from_subroutine_,name,& - & a_err='prec reallocation') + write(ch_err,'(a,i7)') 'Mat asb fail @ level ',i + call psb_errpush(psb_err_internal_error_,name,& + & a_err=ch_err) goto 9999 endif - do i=1,newsz - call prec%precv(i)%move_alloc(tprecv(i),info) - end do - do i=newsz+1, iszv - call prec%precv(i)%free(info) - end do - call move_alloc(tprecv,prec%precv) - ! Ignore errors from transfer - info = psb_success_ - ! - ! Restart - iszv = newsz - ! Fix the pointers, but the level 1 should - ! be treated differently - if (.not.associated(prec%precv(1)%base_a,a)) then - prec%precv(1)%base_a => prec%precv(1)%ac - end if - if (.not.associated(prec%precv(1)%base_desc,desc_a)) then - prec%precv(1)%base_desc => prec%precv(1)%desc_ac + if (stop_hierarchy_loop) then + exit array_build_loop + else + if (i prec%precv(i)%ac - prec%precv(i)%base_desc => prec%precv(i)%desc_ac - ! This is needed when the linmap object has been built - ! reusing the base_desc descriptor through a pointer. - ! With PSBLAS 4 we will have a better solution - if (associated(prec%precv(i)%linmap%p_desc_U)) & - & prec%precv(i)%linmap%p_desc_U => prec%precv(i-1)%base_desc - if (associated(prec%precv(i)%linmap%p_desc_V))& - & prec%precv(i)%linmap%p_desc_V => prec%precv(i)%base_desc - end do + end do array_build_loop + +!!$ write(0,*) ' Done array_build_loop',iszv,newsz,info,psb_errstatus_fatal() + + if (newsz>0) then +!!$ do i=2,newsz +!!$ write(0,*) me,'Newsz Out of array_build_loop ',i,':',& +!!$ & prec%precv(i)%remap_data%desc_ac_pre_remap%is_asb() +!!$ end do +!!$ write(0,*) 'Calling set_nlevs ',newsz + call prec%set_nlevs(newsz) + else +!!$ do i=2, iszv +!!$ write(0,*) me,'Out of array_build_loop ',i,':',& +!!$ & prec%precv(i)%remap_data%desc_ac_pre_remap%is_asb() +!!$ end do end if - write(0,*) ' Done reallocating precv',iszv,newsz,info,psb_errstatus_fatal() + iszv = prec%get_nlevs() call psb_barrier(ctxt) - !write(0,*) 'Should we remap? ' - if (.false.) then - if (amg_get_do_remap().and.(np>=4)) then -!!$ write(0,*) 'Going for remapping ' - if (.true.) then - associate(lv=>prec%precv(iszv), rmp => prec%precv(iszv)%remap_data) - call lv%desc_ac%clone(rmp%desc_ac_pre_remap,info) - call lv%ac%clone(rmp%ac_pre_remap,info) - if (np >= 8) then - call psb_remap(np/4,rmp%desc_ac_pre_remap,rmp%ac_pre_remap,& - & rmp%idest,rmp%isrc,rmp%nrsrc,rmp%naggr,lv%desc_ac,lv%ac,info) - else - call psb_remap(np/2,rmp%desc_ac_pre_remap,rmp%ac_pre_remap,& - & rmp%idest,rmp%isrc,rmp%nrsrc,rmp%naggr,lv%desc_ac,lv%ac,info) - end if -!!$ write(0,*) me,' Out of remapping ',rmp%desc_ac_pre_remap%get_fmt(),' ',& -!!$ & lv%desc_ac%get_fmt(),sum(lv%linmap%naggr),sum(rmp%naggr) - lv%linmap%naggr(:) = rmp%naggr(:) - lv%linmap%p_desc_V => rmp%desc_ac_pre_remap - lv%base_a => lv%ac - lv%base_desc => lv%desc_ac - end associate - end if - end if - end if +!!$ write(0,*) ' Done reallocating precv',iszv,newsz,info +!!$ +!!$ do i=2, iszv +!!$ write(0,*) me,'At end of hierarchy_bld level',i,':',& +!!$ & prec%precv(i)%remap_data%desc_ac_pre_remap%is_asb() +!!$ end do + call psb_barrier(ctxt) + if (info /= psb_success_) then call psb_errpush(psb_err_internal_error_,name,& & a_err='Internal hierarchy build' ) goto 9999 endif - iszv = size(prec%precv) - + iszv = prec%get_nlevs() +!!$ write(0,*) 'Going for cmp_complexity ',& +!!$ & allocated(prec%precv),iszv,size(prec%precv) call prec%cmp_complexity() call prec%cmp_avg_cr() @@ -730,4 +657,49 @@ contains return end subroutine restore_smoothers #endif + + function amg_z_policy_do_remap(ctxt,level,aggsize) result(res) + logical :: res + type(psb_ctxt_type) :: ctxt + integer(psb_ipk_) :: level + integer(psb_lpk_) :: aggsize + res = amg_get_do_remap().and.(level>=2) +!!$ res = .false. + end function amg_z_policy_do_remap + + subroutine amg_z_hierarchy_bld_cmp_newsz(level,iszv,prevsize,& + & nlaggr,casize,mnratio,sizeratio,newsz) + implicit none + integer(psb_ipk_) :: level,iszv,newsz + integer(psb_lpk_) :: nlaggr(:) + integer(psb_lpk_) :: prevsize, casize + real(psb_dpk_) :: mnratio, sizeratio + ! ============================== + integer(psb_lpk_) :: iaggsize + + newsz = 0 + iaggsize = sum(nlaggr) + sizeratio = prevsize + sizeratio = sizeratio/iaggsize +!!$ write(0,*) 'From cmp_newsz: ',iaggsize,casize,& +!!$ & sizeratio,mnratio, level + + if (iaggsize <= casize) newsz = level + if (level == iszv) newsz = level + + if (level>2) then + if (sizeratio < mnratio) then + if (sizeratio > 1) then + newsz = level + else + ! + ! We are not gaining + ! + newsz = level-1 + end if + end if + end if +!!$ write(0,*) 'At end of cmp_newsz ',newsz + end subroutine amg_z_hierarchy_bld_cmp_newsz + end subroutine amg_z_hierarchy_bld diff --git a/amgprec/impl/amg_z_smoothers_bld.f90 b/amgprec/impl/amg_z_smoothers_bld.f90 index 3a081409..e2d0c273 100644 --- a/amgprec/impl/amg_z_smoothers_bld.f90 +++ b/amgprec/impl/amg_z_smoothers_bld.f90 @@ -136,9 +136,9 @@ subroutine amg_z_smoothers_bld(a,desc_a,prec,info,amold,vmold,imold) ! ! Check to ensure all procs have the same ! - iszv = size(prec%precv) + iszv = prec%get_nlevs() call psb_bcast(ctxt,iszv) - if (iszv /= size(prec%precv)) then + if (iszv /= prec%get_nlevs()) then info=psb_err_internal_error_ call psb_errpush(info,name,a_err='Inconsistent size of precv') goto 9999 diff --git a/amgprec/impl/amg_zfile_prec_descr.f90 b/amgprec/impl/amg_zfile_prec_descr.f90 index 2fd88026..8285544b 100644 --- a/amgprec/impl/amg_zfile_prec_descr.f90 +++ b/amgprec/impl/amg_zfile_prec_descr.f90 @@ -129,7 +129,7 @@ subroutine amg_zfile_prec_descr(prec,info,iout,root, verbosity,prefix) ! ensured by amg_precbld). ! if (me == root_) then - nlev = size(prec%precv) + nlev = prec%get_nlevs() do ilev = 1, nlev if (.not.allocated(prec%precv(ilev)%sm)) then info = 3111 @@ -141,6 +141,7 @@ subroutine amg_zfile_prec_descr(prec,info,iout,root, verbosity,prefix) write(iout_,*) write(iout_,'(a,1x,a)') trim(prefix_),'Preconditioner description' + write(iout_,*) 'At level :',1,' we have ',np,' processes' if (nlev == 1) then ! diff --git a/amgprec/impl/amg_zmlprec_aply.f90 b/amgprec/impl/amg_zmlprec_aply.f90 index 86e8850e..50ea8ce0 100644 --- a/amgprec/impl/amg_zmlprec_aply.f90 +++ b/amgprec/impl/amg_zmlprec_aply.f90 @@ -244,10 +244,10 @@ subroutine amg_zmlprec_aply_vect(alpha,p,x,beta,y,desc_data,trans,work,info) if (debug_level >= psb_debug_inner_) & & write(debug_unit,*) me,' ',trim(name),& - & ' Entry ', size(p%precv) + & ' Entry ', p%get_nlevs() trans_ = psb_toupper(trans) - nlev = size(p%precv) + nlev = p%get_nlevs() do_alloc_wrk = .not.allocated(p%precv(1)%wrk) @@ -382,7 +382,7 @@ contains call psb_erractionsave(err_act) debug_unit = psb_get_debug_unit() debug_level = psb_get_debug_level() - nlev = size(p%precv) + nlev = p%get_nlevs() if ((level < 1) .or. (level > nlev)) then call psb_errpush(psb_err_internal_error_,name,& & a_err='wrong call level to inner_ml') @@ -468,7 +468,7 @@ contains call psb_erractionsave(err_act) debug_unit = psb_get_debug_unit() debug_level = psb_get_debug_level() - nlev = size(p%precv) + nlev = p%get_nlevs() if ((level < 1) .or. (level > nlev)) then call psb_errpush(psb_err_internal_error_,name,& & a_err='wrong call level to inner_add') @@ -497,7 +497,8 @@ contains if (allocated(p%precv(level)%sm2a)) then call psb_geaxpby(zone,vx2l,zzero,vy2l,base_desc,info) - sweeps = max(p%precv(level)%parms%sweeps_pre,p%precv(level)%parms%sweeps_post) + sweeps = max(p%precv(level)%parms%sweeps_pre,& + & p%precv(level)%parms%sweeps_post) do k=1, sweeps call p%precv(level)%sm%apply(zone,& & vy2l,zzero,vty,& @@ -593,7 +594,7 @@ contains call psb_erractionsave(err_act) debug_unit = psb_get_debug_unit() debug_level = psb_get_debug_level() - nlev = size(p%precv) + nlev = p%get_nlevs() if ((level < 1) .or. (level > nlev)) then call psb_errpush(psb_err_internal_error_,name,& & a_err='wrong call level to inner_mult') @@ -604,7 +605,7 @@ contains if(debug_level > 1) then write(debug_unit,*) me,' inner_mult at level ',level end if - write(debug_unit,*) me,' inner_mult at level (1):',level,np +!!$ write(debug_unit,*) me,' inner_mult at level (1):',level,np sweeps_post = p%precv(level)%parms%sweeps_post sweeps_pre = p%precv(level)%parms%sweeps_pre pre = ((sweeps_pre>0).and.(trans=='N')).or.((sweeps_post>0).and.(trans/='N')) @@ -614,8 +615,8 @@ contains & vtx => p%precv(level)%wrk%vtx,vty => p%precv(level)%wrk%vty,& & base_a => p%precv(level)%base_a, base_desc=>p%precv(level)%base_desc,& & wv => p%precv(level)%wrk%wv) - write(0,*) 'Inner mult at level (2):',level,' :',me,np,':',& - & size(p%precv(level)%wrk%wv), allocated(p%precv(level)%wrk%wv) +!!$ write(0,*) 'Inner mult at level (2):',level,' :',me,np,':',& +!!$ & size(p%precv(level)%wrk%wv), allocated(p%precv(level)%wrk%wv) if (me >=0) then if (level < nlev) then @@ -831,7 +832,7 @@ contains call psb_erractionsave(err_act) debug_unit = psb_get_debug_unit() debug_level = psb_get_debug_level() - nlev = size(p%precv) + nlev = p%get_nlevs() if ((level < 1) .or. (level > nlev)) then call psb_errpush(psb_err_internal_error_,name,& & a_err='wrong call level to inner_add') diff --git a/amgprec/impl/amg_zprecinit.F90 b/amgprec/impl/amg_zprecinit.F90 index b0f316b8..99b49a1d 100644 --- a/amgprec/impl/amg_zprecinit.F90 +++ b/amgprec/impl/amg_zprecinit.F90 @@ -217,9 +217,7 @@ subroutine amg_zprecinit(ctxt,prec,ptype,info) allocate(amg_z_ilu_solver_type :: prec%precv(ilev_)%sm%sv, stat=info) call prec%precv(ilev_)%default() - case ('ML') - nlev_ = prec%ag_data%max_levs ilev_ = 1 allocate(prec%precv(nlev_),stat=info) @@ -231,6 +229,8 @@ subroutine amg_zprecinit(ctxt,prec,ptype,info) do ilev_ = 1, nlev_ call prec%precv(ilev_)%default() end do + call prec%set_nlevs(nlev_) + call prec%set('ML_CYCLE','VCYCLE',info) call prec%set('SMOOTHER_TYPE','FBGS',info) #if defined(AMG_HAVE_UMF) @@ -246,7 +246,6 @@ subroutine amg_zprecinit(ctxt,prec,ptype,info) write(psb_err_unit,*) name,& &': Warning: Unknown preconditioner type request "',ptype,'"' info = psb_err_pivot_too_small_ - end select call psb_erractionrestore(err_act) diff --git a/amgprec/impl/level/amg_c_base_onelev_map_prol.F90 b/amgprec/impl/level/amg_c_base_onelev_map_prol.F90 index 99d46096..05580489 100644 --- a/amgprec/impl/level/amg_c_base_onelev_map_prol.F90 +++ b/amgprec/impl/level/amg_c_base_onelev_map_prol.F90 @@ -46,6 +46,7 @@ subroutine amg_c_base_onelev_map_prol_v(lv,alpha,vect_v,beta,vect_u,info,work,vt integer(psb_ipk_), intent(out) :: info complex(psb_spk_), optional :: work(:) type(psb_c_vect_type), optional, target, intent(inout) :: vtx,vty + type(psb_c_vect_type), pointer :: vtx_ !!$ write(0,*) 'New map_rstr',lv%remap_data%ac_pre_remap%is_asb() diff --git a/amgprec/impl/level/amg_c_base_onelev_map_rstr.F90 b/amgprec/impl/level/amg_c_base_onelev_map_rstr.F90 index c49e8f50..6beb0e0b 100644 --- a/amgprec/impl/level/amg_c_base_onelev_map_rstr.F90 +++ b/amgprec/impl/level/amg_c_base_onelev_map_rstr.F90 @@ -48,6 +48,7 @@ subroutine amg_c_base_onelev_map_rstr_v(lv,alpha,vect_u,beta,vect_v,info,& complex(psb_spk_), optional :: work(:) type(psb_c_vect_type), optional, target, intent(inout) :: vtx,vty type(psb_c_vect_type), pointer :: vty_ + integer(psb_mpk_) :: me, np !!$ write(0,*) 'New map_rstr',lv%remap_data%ac_pre_remap%is_asb() if (present(vty)) then @@ -64,7 +65,7 @@ subroutine amg_c_base_onelev_map_rstr_v(lv,alpha,vect_u,beta,vect_v,info,& block type(psb_ctxt_type) :: ctxt, rctxt integer(psb_mpk_) :: i,j,ip, idest, nsrc, nrl, kp - integer(psb_mpk_) :: me, np, rme, rnp + integer(psb_mpk_) :: rme, rnp complex(psb_spk_), allocatable :: rsnd(:), rrcv(:) type(psb_c_vect_type) :: tv @@ -72,7 +73,7 @@ subroutine amg_c_base_onelev_map_rstr_v(lv,alpha,vect_u,beta,vect_v,info,& call psb_info(ctxt,me,np) rctxt = lv%desc_ac%get_ctxt() call psb_info(rctxt,rme,rnp) -!!$ write(0,*) 'New context ',rme,rnp +!!$ write(0,*) 'New context map rstr',rme,rnp,me,np idest = lv%remap_data%idest associate(isrc => lv%remap_data%isrc, nrsrc => lv%remap_data%nrsrc) !!$ write(0,*) 'Should apply maps, then send data from ',me,' to ',idest @@ -81,15 +82,16 @@ subroutine amg_c_base_onelev_map_rstr_v(lv,alpha,vect_u,beta,vect_v,info,& nrl = lv%remap_data%desc_ac_pre_remap%get_local_rows() call psb_geall(tv,lv%remap_data%desc_ac_pre_remap,info) call psb_geasb(tv,lv%remap_data%desc_ac_pre_remap,info,mold=vect_u%v) - write(0,*) me,' remap map_rstr calling U2V: ',me,np,rme,rnp,tv%get_nrows() - flush(0) +!!$ write(0,*) me,' remap map_rstr calling U2V: ',me,np,rme,rnp,tv%get_nrows(),& +!!$ & psb_errstatus_fatal() +!!$ flush(0) call psb_barrier(ctxt) call lv%linmap%map_U2V(alpha,vect_u,beta,tv,info,& & work=work,vtx=vtx,vty=vty_) call tv%sync() !rsnd = tv%get_vect() !call psb_snd(ctxt,rsnd(1:nrl),idest) - write(0,*) me,' map_rstr sending ',me,idest +!!$ write(0,*) me,' map_rstr sending ',me,idest,psb_errstatus_fatal() call psb_snd(ctxt,tv%v%v(1:nrl),idest) if (rme >=0) then allocate(rrcv(sum(nrsrc))) @@ -98,29 +100,28 @@ subroutine amg_c_base_onelev_map_rstr_v(lv,alpha,vect_u,beta,vect_v,info,& do i = 1,size(isrc) ip = isrc(i) nrl = nrsrc(i) - write(0,*) me,' map_rstr receiving',rme,i,ip +!!$ write(0,*) me,' map_rstr receiving',rme,ip,psb_errstatus_fatal() call psb_rcv(ctxt,rrcv(kp+1:kp+nrl),ip) kp = kp + nrl end do call vect_v%set_vect(rrcv) end if end associate -!!$ write(0,*) me, ' Restrictor with remap done ' +!!$ write(0,*) me, ' Restrictor with remap done ',psb_errstatus_fatal() end block else ! Default transfer block type(psb_ctxt_type) :: ctxt, rctxt - integer(psb_mpk_) :: me, np ctxt = lv%linmap%p_desc_U%get_ctxt() call psb_info(ctxt,me,np) - write(0,*) me,' map_rstr calling U2V: ',me,np +!!$ write(0,*) me,' map_rstr calling U2V: ',me,np call lv%linmap%map_U2V(alpha,vect_u,beta,vect_v,info,& & work=work,vtx=vtx,vty=vty_) end block end if - +!!$ write(0,*) me, 'End of restriction ',info,psb_errstatus_fatal() end subroutine amg_c_base_onelev_map_rstr_v subroutine amg_c_base_onelev_map_rstr_a(lv,alpha,u,beta,v,info,work) diff --git a/amgprec/impl/level/amg_c_base_onelev_mat_asb.f90 b/amgprec/impl/level/amg_c_base_onelev_mat_asb.f90 index f5b3b738..1a762987 100644 --- a/amgprec/impl/level/amg_c_base_onelev_mat_asb.f90 +++ b/amgprec/impl/level/amg_c_base_onelev_mat_asb.f90 @@ -168,7 +168,7 @@ subroutine amg_c_base_onelev_mat_asb(lv,a,desc_a,ilaggr,nlaggr,t_prol,info) if (do_timings) call psb_tic(idx_mapbld) if (info == psb_success_) call lv%ac%cscnv(info,type='csr',dupl=psb_dupl_add_) - if (info == psb_success_) call lv%aggr%bld_map(desc_a, lv%desc_ac,& + if (info == psb_success_) call lv%aggr%bld_linmap(desc_a, lv%desc_ac,& & ilaggr,nlaggr,op_restr,op_prol,lv%linmap,info) if (do_timings) call psb_toc(idx_mapbld) if(info /= psb_success_) then diff --git a/amgprec/impl/level/amg_d_base_onelev_map_rstr.F90 b/amgprec/impl/level/amg_d_base_onelev_map_rstr.F90 index 74e3d29f..ecb9d4cb 100644 --- a/amgprec/impl/level/amg_d_base_onelev_map_rstr.F90 +++ b/amgprec/impl/level/amg_d_base_onelev_map_rstr.F90 @@ -35,6 +35,7 @@ ! POSSIBILITY OF SUCH DAMAGE. ! ! + subroutine amg_d_base_onelev_map_rstr_v(lv,alpha,vect_u,beta,vect_v,info,& & work,vtx,vty) use psb_base_mod @@ -47,9 +48,9 @@ subroutine amg_d_base_onelev_map_rstr_v(lv,alpha,vect_u,beta,vect_v,info,& real(psb_dpk_), optional :: work(:) type(psb_d_vect_type), optional, target, intent(inout) :: vtx,vty type(psb_d_vect_type), pointer :: vty_ - integer(psb_mpk_) :: me, np + integer(psb_mpk_) :: me, np -!!$ write(0,*) 'New map_rstr ',lv%remap_data%ac_pre_remap%is_asb() +!!$ write(0,*) 'New map_rstr',lv%remap_data%ac_pre_remap%is_asb() if (present(vty)) then vty_ => vty else @@ -60,7 +61,7 @@ subroutine amg_d_base_onelev_map_rstr_v(lv,alpha,vect_u,beta,vect_v,info,& ! ! Remap has happened, deal with it ! -!!$ write(0,*) 'Level map_rstr with remapping ' +!!$ write(0,*) 'Remap handling not implemented yet ' block type(psb_ctxt_type) :: ctxt, rctxt integer(psb_mpk_) :: i,j,ip, idest, nsrc, nrl, kp @@ -87,7 +88,7 @@ subroutine amg_d_base_onelev_map_rstr_v(lv,alpha,vect_u,beta,vect_v,info,& call psb_barrier(ctxt) call lv%linmap%map_U2V(alpha,vect_u,beta,tv,info,& & work=work,vtx=vtx,vty=vty_) - !call tv%sync() + call tv%sync() !rsnd = tv%get_vect() !call psb_snd(ctxt,rsnd(1:nrl),idest) !!$ write(0,*) me,' map_rstr sending ',me,idest,psb_errstatus_fatal() @@ -120,7 +121,7 @@ subroutine amg_d_base_onelev_map_rstr_v(lv,alpha,vect_u,beta,vect_v,info,& & work=work,vtx=vtx,vty=vty_) end block end if -!!$ write(0,*) me, 'End of restriction ',info,psb_errstatus_fatal() +!!$ write(0,*) me, 'End of restriction ',info,psb_errstatus_fatal() end subroutine amg_d_base_onelev_map_rstr_v subroutine amg_d_base_onelev_map_rstr_a(lv,alpha,u,beta,v,info,work) diff --git a/amgprec/impl/level/amg_d_base_onelev_mat_asb.f90 b/amgprec/impl/level/amg_d_base_onelev_mat_asb.f90 index 4584bc4f..7f94407a 100644 --- a/amgprec/impl/level/amg_d_base_onelev_mat_asb.f90 +++ b/amgprec/impl/level/amg_d_base_onelev_mat_asb.f90 @@ -168,7 +168,7 @@ subroutine amg_d_base_onelev_mat_asb(lv,a,desc_a,ilaggr,nlaggr,t_prol,info) if (do_timings) call psb_tic(idx_mapbld) if (info == psb_success_) call lv%ac%cscnv(info,type='csr',dupl=psb_dupl_add_) - if (info == psb_success_) call lv%aggr%bld_map(desc_a, lv%desc_ac,& + if (info == psb_success_) call lv%aggr%bld_linmap(desc_a, lv%desc_ac,& & ilaggr,nlaggr,op_restr,op_prol,lv%linmap,info) if (do_timings) call psb_toc(idx_mapbld) if(info /= psb_success_) then diff --git a/amgprec/impl/level/amg_s_base_onelev_map_rstr.F90 b/amgprec/impl/level/amg_s_base_onelev_map_rstr.F90 index e81d5b31..ce8cc74a 100644 --- a/amgprec/impl/level/amg_s_base_onelev_map_rstr.F90 +++ b/amgprec/impl/level/amg_s_base_onelev_map_rstr.F90 @@ -48,6 +48,7 @@ subroutine amg_s_base_onelev_map_rstr_v(lv,alpha,vect_u,beta,vect_v,info,& real(psb_spk_), optional :: work(:) type(psb_s_vect_type), optional, target, intent(inout) :: vtx,vty type(psb_s_vect_type), pointer :: vty_ + integer(psb_mpk_) :: me, np !!$ write(0,*) 'New map_rstr',lv%remap_data%ac_pre_remap%is_asb() if (present(vty)) then @@ -64,7 +65,7 @@ subroutine amg_s_base_onelev_map_rstr_v(lv,alpha,vect_u,beta,vect_v,info,& block type(psb_ctxt_type) :: ctxt, rctxt integer(psb_mpk_) :: i,j,ip, idest, nsrc, nrl, kp - integer(psb_mpk_) :: me, np, rme, rnp + integer(psb_mpk_) :: rme, rnp real(psb_spk_), allocatable :: rsnd(:), rrcv(:) type(psb_s_vect_type) :: tv @@ -72,7 +73,7 @@ subroutine amg_s_base_onelev_map_rstr_v(lv,alpha,vect_u,beta,vect_v,info,& call psb_info(ctxt,me,np) rctxt = lv%desc_ac%get_ctxt() call psb_info(rctxt,rme,rnp) -!!$ write(0,*) 'New context ',rme,rnp +!!$ write(0,*) 'New context map rstr',rme,rnp,me,np idest = lv%remap_data%idest associate(isrc => lv%remap_data%isrc, nrsrc => lv%remap_data%nrsrc) !!$ write(0,*) 'Should apply maps, then send data from ',me,' to ',idest @@ -81,15 +82,16 @@ subroutine amg_s_base_onelev_map_rstr_v(lv,alpha,vect_u,beta,vect_v,info,& nrl = lv%remap_data%desc_ac_pre_remap%get_local_rows() call psb_geall(tv,lv%remap_data%desc_ac_pre_remap,info) call psb_geasb(tv,lv%remap_data%desc_ac_pre_remap,info,mold=vect_u%v) - write(0,*) me,' remap map_rstr calling U2V: ',me,np,rme,rnp,tv%get_nrows() - flush(0) +!!$ write(0,*) me,' remap map_rstr calling U2V: ',me,np,rme,rnp,tv%get_nrows(),& +!!$ & psb_errstatus_fatal() +!!$ flush(0) call psb_barrier(ctxt) call lv%linmap%map_U2V(alpha,vect_u,beta,tv,info,& & work=work,vtx=vtx,vty=vty_) call tv%sync() !rsnd = tv%get_vect() !call psb_snd(ctxt,rsnd(1:nrl),idest) - write(0,*) me,' map_rstr sending ',me,idest +!!$ write(0,*) me,' map_rstr sending ',me,idest,psb_errstatus_fatal() call psb_snd(ctxt,tv%v%v(1:nrl),idest) if (rme >=0) then allocate(rrcv(sum(nrsrc))) @@ -98,29 +100,28 @@ subroutine amg_s_base_onelev_map_rstr_v(lv,alpha,vect_u,beta,vect_v,info,& do i = 1,size(isrc) ip = isrc(i) nrl = nrsrc(i) - write(0,*) me,' map_rstr receiving',rme,i,ip +!!$ write(0,*) me,' map_rstr receiving',rme,ip,psb_errstatus_fatal() call psb_rcv(ctxt,rrcv(kp+1:kp+nrl),ip) kp = kp + nrl end do call vect_v%set_vect(rrcv) end if end associate -!!$ write(0,*) me, ' Restrictor with remap done ' +!!$ write(0,*) me, ' Restrictor with remap done ',psb_errstatus_fatal() end block else ! Default transfer block type(psb_ctxt_type) :: ctxt, rctxt - integer(psb_mpk_) :: me, np ctxt = lv%linmap%p_desc_U%get_ctxt() call psb_info(ctxt,me,np) - write(0,*) me,' map_rstr calling U2V: ',me,np +!!$ write(0,*) me,' map_rstr calling U2V: ',me,np call lv%linmap%map_U2V(alpha,vect_u,beta,vect_v,info,& & work=work,vtx=vtx,vty=vty_) end block end if - +!!$ write(0,*) me, 'End of restriction ',info,psb_errstatus_fatal() end subroutine amg_s_base_onelev_map_rstr_v subroutine amg_s_base_onelev_map_rstr_a(lv,alpha,u,beta,v,info,work) diff --git a/amgprec/impl/level/amg_s_base_onelev_mat_asb.f90 b/amgprec/impl/level/amg_s_base_onelev_mat_asb.f90 index 5ecf7ac8..b613651b 100644 --- a/amgprec/impl/level/amg_s_base_onelev_mat_asb.f90 +++ b/amgprec/impl/level/amg_s_base_onelev_mat_asb.f90 @@ -168,7 +168,7 @@ subroutine amg_s_base_onelev_mat_asb(lv,a,desc_a,ilaggr,nlaggr,t_prol,info) if (do_timings) call psb_tic(idx_mapbld) if (info == psb_success_) call lv%ac%cscnv(info,type='csr',dupl=psb_dupl_add_) - if (info == psb_success_) call lv%aggr%bld_map(desc_a, lv%desc_ac,& + if (info == psb_success_) call lv%aggr%bld_linmap(desc_a, lv%desc_ac,& & ilaggr,nlaggr,op_restr,op_prol,lv%linmap,info) if (do_timings) call psb_toc(idx_mapbld) if(info /= psb_success_) then diff --git a/amgprec/impl/level/amg_z_base_onelev_map_rstr.F90 b/amgprec/impl/level/amg_z_base_onelev_map_rstr.F90 index ffcd9a72..808752fd 100644 --- a/amgprec/impl/level/amg_z_base_onelev_map_rstr.F90 +++ b/amgprec/impl/level/amg_z_base_onelev_map_rstr.F90 @@ -48,6 +48,7 @@ subroutine amg_z_base_onelev_map_rstr_v(lv,alpha,vect_u,beta,vect_v,info,& complex(psb_dpk_), optional :: work(:) type(psb_z_vect_type), optional, target, intent(inout) :: vtx,vty type(psb_z_vect_type), pointer :: vty_ + integer(psb_mpk_) :: me, np !!$ write(0,*) 'New map_rstr',lv%remap_data%ac_pre_remap%is_asb() if (present(vty)) then @@ -64,7 +65,7 @@ subroutine amg_z_base_onelev_map_rstr_v(lv,alpha,vect_u,beta,vect_v,info,& block type(psb_ctxt_type) :: ctxt, rctxt integer(psb_mpk_) :: i,j,ip, idest, nsrc, nrl, kp - integer(psb_mpk_) :: me, np, rme, rnp + integer(psb_mpk_) :: rme, rnp complex(psb_dpk_), allocatable :: rsnd(:), rrcv(:) type(psb_z_vect_type) :: tv @@ -72,7 +73,7 @@ subroutine amg_z_base_onelev_map_rstr_v(lv,alpha,vect_u,beta,vect_v,info,& call psb_info(ctxt,me,np) rctxt = lv%desc_ac%get_ctxt() call psb_info(rctxt,rme,rnp) -!!$ write(0,*) 'New context ',rme,rnp +!!$ write(0,*) 'New context map rstr',rme,rnp,me,np idest = lv%remap_data%idest associate(isrc => lv%remap_data%isrc, nrsrc => lv%remap_data%nrsrc) !!$ write(0,*) 'Should apply maps, then send data from ',me,' to ',idest @@ -81,15 +82,16 @@ subroutine amg_z_base_onelev_map_rstr_v(lv,alpha,vect_u,beta,vect_v,info,& nrl = lv%remap_data%desc_ac_pre_remap%get_local_rows() call psb_geall(tv,lv%remap_data%desc_ac_pre_remap,info) call psb_geasb(tv,lv%remap_data%desc_ac_pre_remap,info,mold=vect_u%v) - write(0,*) me,' remap map_rstr calling U2V: ',me,np,rme,rnp,tv%get_nrows() - flush(0) +!!$ write(0,*) me,' remap map_rstr calling U2V: ',me,np,rme,rnp,tv%get_nrows(),& +!!$ & psb_errstatus_fatal() +!!$ flush(0) call psb_barrier(ctxt) call lv%linmap%map_U2V(alpha,vect_u,beta,tv,info,& & work=work,vtx=vtx,vty=vty_) call tv%sync() !rsnd = tv%get_vect() !call psb_snd(ctxt,rsnd(1:nrl),idest) - write(0,*) me,' map_rstr sending ',me,idest +!!$ write(0,*) me,' map_rstr sending ',me,idest,psb_errstatus_fatal() call psb_snd(ctxt,tv%v%v(1:nrl),idest) if (rme >=0) then allocate(rrcv(sum(nrsrc))) @@ -98,29 +100,28 @@ subroutine amg_z_base_onelev_map_rstr_v(lv,alpha,vect_u,beta,vect_v,info,& do i = 1,size(isrc) ip = isrc(i) nrl = nrsrc(i) - write(0,*) me,' map_rstr receiving',rme,i,ip +!!$ write(0,*) me,' map_rstr receiving',rme,ip,psb_errstatus_fatal() call psb_rcv(ctxt,rrcv(kp+1:kp+nrl),ip) kp = kp + nrl end do call vect_v%set_vect(rrcv) end if end associate -!!$ write(0,*) me, ' Restrictor with remap done ' +!!$ write(0,*) me, ' Restrictor with remap done ',psb_errstatus_fatal() end block else ! Default transfer block type(psb_ctxt_type) :: ctxt, rctxt - integer(psb_mpk_) :: me, np ctxt = lv%linmap%p_desc_U%get_ctxt() call psb_info(ctxt,me,np) - write(0,*) me,' map_rstr calling U2V: ',me,np +!!$ write(0,*) me,' map_rstr calling U2V: ',me,np call lv%linmap%map_U2V(alpha,vect_u,beta,vect_v,info,& & work=work,vtx=vtx,vty=vty_) end block end if - +!!$ write(0,*) me, 'End of restriction ',info,psb_errstatus_fatal() end subroutine amg_z_base_onelev_map_rstr_v subroutine amg_z_base_onelev_map_rstr_a(lv,alpha,u,beta,v,info,work) diff --git a/amgprec/impl/level/amg_z_base_onelev_mat_asb.f90 b/amgprec/impl/level/amg_z_base_onelev_mat_asb.f90 index a2e93387..9cadc259 100644 --- a/amgprec/impl/level/amg_z_base_onelev_mat_asb.f90 +++ b/amgprec/impl/level/amg_z_base_onelev_mat_asb.f90 @@ -168,7 +168,7 @@ subroutine amg_z_base_onelev_mat_asb(lv,a,desc_a,ilaggr,nlaggr,t_prol,info) if (do_timings) call psb_tic(idx_mapbld) if (info == psb_success_) call lv%ac%cscnv(info,type='csr',dupl=psb_dupl_add_) - if (info == psb_success_) call lv%aggr%bld_map(desc_a, lv%desc_ac,& + if (info == psb_success_) call lv%aggr%bld_linmap(desc_a, lv%desc_ac,& & ilaggr,nlaggr,op_restr,op_prol,lv%linmap,info) if (do_timings) call psb_toc(idx_mapbld) if(info /= psb_success_) then