From 42313bb843b7651f2708d17e1ab4725c154d79aa Mon Sep 17 00:00:00 2001 From: sfilippone Date: Sat, 11 Jul 2026 12:30:28 +0200 Subject: [PATCH] First D remapping version. Needs full testing an work on clone --- amgprec/amg_d_onelev_mod.f90 | 27 +- amgprec/amg_d_prec_type.f90 | 47 ++- amgprec/impl/amg_cmlprec_aply.f90 | 25 +- amgprec/impl/amg_d_extprol_bld.F90 | 8 +- amgprec/impl/amg_d_hierarchy_bld.F90 | 328 +++++++++--------- amgprec/impl/amg_d_hierarchy_rebld.f90 | 2 +- amgprec/impl/amg_d_smoothers_bld.f90 | 4 +- amgprec/impl/amg_dcprecset.F90 | 6 +- amgprec/impl/amg_dfile_prec_descr.f90 | 4 +- amgprec/impl/amg_dfile_prec_memory_use.f90 | 2 +- amgprec/impl/amg_dmlprec_aply.f90 | 65 ++-- amgprec/impl/amg_dmlprec_aply_a.f90 | 10 +- amgprec/impl/amg_dprecinit.F90 | 2 +- amgprec/impl/amg_smlprec_aply.f90 | 25 +- amgprec/impl/amg_zmlprec_aply.f90 | 25 +- .../impl/level/amg_c_base_onelev_map_prol.F90 | 12 +- .../impl/level/amg_c_base_onelev_map_rstr.F90 | 22 +- .../impl/level/amg_d_base_onelev_map_prol.F90 | 13 +- .../impl/level/amg_d_base_onelev_map_rstr.F90 | 46 ++- .../impl/level/amg_s_base_onelev_map_prol.F90 | 13 +- .../impl/level/amg_s_base_onelev_map_rstr.F90 | 22 +- .../impl/level/amg_z_base_onelev_map_prol.F90 | 13 +- .../impl/level/amg_z_base_onelev_map_rstr.F90 | 22 +- samples/advanced/remap/amg_d_pde3d.F90 | 12 + 24 files changed, 444 insertions(+), 311 deletions(-) diff --git a/amgprec/amg_d_onelev_mod.f90 b/amgprec/amg_d_onelev_mod.f90 index b73ee5f5..29b7f942 100644 --- a/amgprec/amg_d_onelev_mod.f90 +++ b/amgprec/amg_d_onelev_mod.f90 @@ -183,6 +183,7 @@ module amg_d_onelev_mod 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 end type amg_d_remap_data_type type amg_d_onelev_type @@ -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,7 +810,8 @@ 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 @@ -997,4 +1001,25 @@ contains call psb_safe_ab_cpy(rmp%nrsrc,remap_out%nrsrc,info) end subroutine d_remap_data_clone + subroutine d_remap_move_alloc(rmp, remap_out, info) + use psb_base_mod + implicit none + ! Arguments + class(amg_d_remap_data_type), target, intent(inout) :: rmp + class(amg_d_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 d_remap_move_alloc + end module amg_d_onelev_mod diff --git a/amgprec/amg_d_prec_type.f90 b/amgprec/amg_d_prec_type.f90 index 4ff41d15..40d53dd0 100644 --- a/amgprec/amg_d_prec_type.f90 +++ b/amgprec/amg_d_prec_type.f90 @@ -97,11 +97,11 @@ 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 ! type(amg_d_onelev_type), allocatable :: precv(:) + integer(psb_ipk_) :: nlevs contains procedure, pass(prec) :: psb_d_apply2_vect => amg_d_apply2_vect procedure, pass(prec) :: psb_d_apply1_vect => amg_d_apply1_vect @@ -119,6 +119,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) :: get_nlevs => amg_d_get_nlevs procedure, pass(prec) :: get_nzeros => amg_d_get_nzeros procedure, pass(prec) :: sizeof => amg_dprec_sizeof @@ -141,6 +142,7 @@ 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,& @@ -438,7 +440,17 @@ contains 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_d_get_nlevs + + subroutine amg_d_set_nlevs(prec,nl) + implicit none + class(amg_dprec_type), intent(inout) :: prec + integer(psb_ipk_) :: nl + prec%nlevs = nl + end subroutine amg_d_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 +521,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 +531,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 @@ -562,13 +577,11 @@ contains 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_d_cmp_avg_cr @@ -667,7 +680,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 @@ -857,7 +870,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 @@ -893,7 +906,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 @@ -930,8 +943,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 +953,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 +1033,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) @@ -1058,7 +1073,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_cmlprec_aply.f90 b/amgprec/impl/amg_cmlprec_aply.f90 index 4f724a46..c8dac94a 100644 --- a/amgprec/impl/amg_cmlprec_aply.f90 +++ b/amgprec/impl/amg_cmlprec_aply.f90 @@ -527,8 +527,7 @@ contains ! Apply the restriction call p%precv(level+1)%map_rstr(cone,vx2l,& & czero,p%precv(level+1)%wrk%vx2l,& - & info,work=work,& - & vtx=wv(1),vty=p%precv(level+1)%wrk%wv(1)) + & info,work=work,vtx=wv(1)) if (info /= psb_success_) then call psb_errpush(psb_err_internal_error_,name,& & a_err='Error during restriction') @@ -547,8 +546,7 @@ contains ! call p%precv(level+1)%map_prol(cone,& & p%precv(level+1)%wrk%vy2l, cone,vy2l,& - & info,work=work,& - & vtx=p%precv(level+1)%wrk%wv(1),vty=wv(1)) + & info,work=work, vty=wv(1)) if (info /= psb_success_) then call psb_errpush(psb_err_internal_error_,name,& & a_err='Error during prolongation') @@ -667,8 +665,7 @@ contains call p%precv(level+1)%map_rstr(cone,vty,& & czero,p%precv(level+1)%wrk%vx2l,& - & info,work=work,& - & vtx=wv(1),vty=p%precv(level+1)%wrk%wv(1)) + & info,work=work,vtx=wv(1)) if (info /= psb_success_) then call psb_errpush(psb_err_internal_error_,name,& & a_err='Error during restriction') @@ -678,8 +675,7 @@ contains ! Shortcut: just transfer x2l. call p%precv(level+1)%map_rstr(cone,vx2l,& & czero,p%precv(level+1)%wrk%vx2l,& - & info,work=work,& - & vtx=wv(1),vty=p%precv(level+1)%wrk%wv(1)) + & info,work=work,vtx=wv(1)) if (info /= psb_success_) then call psb_errpush(psb_err_internal_error_,name,& & a_err='Error during restriction') @@ -694,8 +690,7 @@ contains ! call p%precv(level+1)%map_prol(cone,& & p%precv(level+1)%wrk%vy2l,cone,vy2l,& - & info,work=work,& - & vtx=p%precv(level+1)%wrk%wv(1),vty=wv(1)) + & info,work=work,vty=wv(1)) if (info /= psb_success_) then call psb_errpush(psb_err_internal_error_,name,& & a_err='Error during prolongation') @@ -714,7 +709,7 @@ contains if (info == psb_success_) & & call p%precv(level+1)%map_rstr(cone,vty,& & czero,p%precv(level+1)%wrk%vx2l,info,work=work,& - & vtx=wv(1),vty=p%precv(level+1)%wrk%wv(1)) + & vtx=wv(1)) if (info /= psb_success_) then call psb_errpush(psb_err_internal_error_,name,& & a_err='Error during W-cycle restriction') @@ -725,8 +720,7 @@ contains if (info == psb_success_) call p%precv(level+1)%map_prol(cone, & & p%precv(level+1)%wrk%vy2l,cone,vy2l,& - & info,work=work,& - & vtx=p%precv(level+1)%wrk%wv(1),vty=wv(1)) + & info,work=work,vty=wv(1)) if (info /= psb_success_) then call psb_errpush(psb_err_internal_error_,name,& @@ -914,7 +908,7 @@ contains call p%precv(level + 1)%map_rstr(cone,vty,& & czero,p%precv(level + 1)%wrk%vx2l,& &info,work=work,& - & vtx=wv(1),vty=p%precv(level+1)%wrk%wv(1)) + & vtx=wv(1)) if (info /= psb_success_) then call psb_errpush(psb_err_internal_error_,name,& @@ -949,8 +943,7 @@ contains ! call p%precv(level+1)%map_prol(cone,& & p%precv(level+1)%wrk%vy2l,cone,vy2l,& - & info,work=work,& - & vtx=p%precv(level+1)%wrk%wv(1),vty=wv(1)) + & info,work=work,vty=wv(1)) if (info /= psb_success_) then call psb_errpush(psb_err_internal_error_,name,& diff --git a/amgprec/impl/amg_d_extprol_bld.F90 b/amgprec/impl/amg_d_extprol_bld.F90 index 68fa4726..e72e196d 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 = size(p%precv) + iszv = p%get_nlevs() 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 /= size(p%precv)) then + if (iszv /= p%get_nlevs()) 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 = size(p%precv) + iszv = p%get_nlevs() 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 = size(p%precv) + iszv = p%get_nlevs() 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 9f6c6df3..96f5f034 100644 --- a/amgprec/impl/amg_d_hierarchy_bld.F90 +++ b/amgprec/impl/amg_d_hierarchy_bld.F90 @@ -82,7 +82,7 @@ subroutine amg_d_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_d_base_smoother_type), allocatable :: coarse_sm, med_sm, & @@ -98,6 +98,7 @@ subroutine amg_d_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 +142,7 @@ subroutine amg_d_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 +168,7 @@ subroutine amg_d_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 +183,7 @@ subroutine amg_d_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 +231,6 @@ subroutine amg_d_hierarchy_bld(a,desc_a,prec,info,cpymat) casize = mncsize end if prec%ag_data%target_coarse_size = casize - nplevs = max(itwo,mxplevs) ! @@ -288,9 +289,9 @@ subroutine amg_d_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 - ! ! Finest level first; create a GEN_BLOCK ! copy of the descriptor. @@ -305,6 +306,7 @@ subroutine amg_d_hierarchy_bld(a,desc_a,prec,info,cpymat) prec%precv(1)%base_desc => prec%precv(1)%desc_ac newsz = 0 + stop_hierarchy_loop = .false. array_build_loop: do i=2, iszv ! ! Check on the iprcparm contents: they should be the same @@ -352,62 +354,23 @@ subroutine amg_d_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_hierarchy_bld_newsz(i,iszv,& + & desc_a%get_global_rows(),nlaggr,mnaggratio,sizeratio,newsz) else - sizeratio = sum(prec%precv(i-1)%linmap%naggr)/sizeratio + call amg_hierarchy_bld_newsz(i,iszv,& + & sum(prec%precv(i-1)%linmap%naggr),nlaggr,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 +404,166 @@ subroutine amg_d_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,& +!!$ & prec%precv(i)%remap_data%desc_ac_pre_remap%is_asb() +!!$ exit array_build_loop + 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 + ! + ! Do we want to remap onto a smaller subset of processes? + ! + + if (amg_policy_do_remap(level,sum(nlaggr))) then block type(psb_ctxt_type) :: lctxt integer(psb_ipk_) :: lme,lnp - lctxt = prec%precv(i)%desc_ac%get_ctxt() + lctxt = prec%precv(level)%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) + 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 + write(0,*) 'Second Check on out remapping ',level,& + & prec%precv(level)%remap_data%desc_ac_pre_remap%is_asb() 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) + 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) + if (stop_hierarchy_loop) then + exit array_build_loop + else + if (i0) 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 - do i=newsz+1, iszv - call prec%precv(i)%free(info) + 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 - 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 + end if + iszv = prec%get_nlevs() + call psb_barrier(ctxt) + + if (.false.) then + if (newsz > 0) then + ! + ! We exited early from the build loop, need to fix + ! the size. + ! + allocate(tprecv(newsz),stat=info) + if (info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,& + & a_err='prec reallocation') + 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 + end if + do i=2, iszv + prec%precv(i)%base_a => 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 if + else 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 end if - do i=2, iszv - prec%precv(i)%base_a => 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 if - write(0,*) ' Done reallocating precv',iszv,newsz,info,psb_errstatus_fatal() + 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) - !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 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 +703,45 @@ contains return end subroutine restore_smoothers #endif + + function amg_policy_do_remap(level,aggsize) result(res) + logical :: res + integer(psb_ipk_) :: level + integer(psb_lpk_) :: aggsize + res = amg_get_do_remap().and.(level>=2) +!!$ res = .false. + end function amg_policy_do_remap + + subroutine amg_hierarchy_bld_newsz(level,iszv,prevsize,& + & nlaggr,mnratio,sizeratio,newsz) + implicit none + integer(psb_ipk_) :: level,iszv,newsz + integer(psb_lpk_) :: nlaggr(:) + integer(psb_lpk_) :: prevsize + real(psb_dpk_) :: mnratio, sizeratio + ! ============================== + integer(psb_lpk_) :: iaggsize, casize + + newsz = 0 + iaggsize = sum(nlaggr) + sizeratio = prevsize + sizeratio = sizeratio/iaggsize + + if (iaggsize <= casize) newsz = level + if (level == iszv) newsz = level + + if (level>2) then + if (sizeratio < mnaggratio) then + if (sizeratio > 1) then + newsz = level + else + ! + ! We are not gaining + ! + newsz = level-1 + end if + end if + end if + end subroutine amg_hierarchy_bld_newsz + end subroutine amg_d_hierarchy_bld diff --git a/amgprec/impl/amg_d_hierarchy_rebld.f90 b/amgprec/impl/amg_d_hierarchy_rebld.f90 index 885f95fa..def6534e 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 = size(prec%precv) + iszv = prec%get_nlevs() do i=2, iszv call prec%precv(i-1)%base_a%cp_to(acsr) diff --git a/amgprec/impl/amg_d_smoothers_bld.f90 b/amgprec/impl/amg_d_smoothers_bld.f90 index 27ff39e0..4214d8b5 100644 --- a/amgprec/impl/amg_d_smoothers_bld.f90 +++ b/amgprec/impl/amg_d_smoothers_bld.f90 @@ -136,9 +136,9 @@ subroutine amg_d_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_dcprecset.F90 b/amgprec/impl/amg_dcprecset.F90 index 350ddffd..18e6f324 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_ = size(p%precv) + nlev_ = p%get_nlevs() if (present(ilev)) then ilev_ = ilev @@ -376,7 +376,7 @@ subroutine amg_dcprecsetc(p,what,string,info,ilev,ilmax,pos,idx) return endif - nlev_ = size(p%precv) + nlev_ = p%get_nlevs() 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_ = size(p%precv) + nlev_ = p%get_nlevs() if (present(ilev)) then ilev_ = ilev diff --git a/amgprec/impl/amg_dfile_prec_descr.f90 b/amgprec/impl/amg_dfile_prec_descr.f90 index 8d429643..b27d2c4c 100644 --- a/amgprec/impl/amg_dfile_prec_descr.f90 +++ b/amgprec/impl/amg_dfile_prec_descr.f90 @@ -129,7 +129,9 @@ subroutine amg_dfile_prec_descr(prec,info,iout,root, verbosity,prefix) ! ensured by amg_precbld). ! if (me == root_) then - nlev = size(prec%precv) + nlev = prec%get_nlevs() + write(0,*) 'From file_prec_descr: nlev ',nlev + flush(0) do ilev = 1, nlev if (.not.allocated(prec%precv(ilev)%sm)) then info = 3111 diff --git a/amgprec/impl/amg_dfile_prec_memory_use.f90 b/amgprec/impl/amg_dfile_prec_memory_use.f90 index 3757d6a4..2b689cea 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 = size(prec%precv) + nlev = prec%get_nlevs() 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 bcddbcc1..633ebd5b 100644 --- a/amgprec/impl/amg_dmlprec_aply.f90 +++ b/amgprec/impl/amg_dmlprec_aply.f90 @@ -244,10 +244,10 @@ subroutine amg_dmlprec_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(done,vx2l,dzero,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(done,& & vy2l,dzero,vty,& @@ -527,8 +528,7 @@ contains ! Apply the restriction call p%precv(level+1)%map_rstr(done,vx2l,& & dzero,p%precv(level+1)%wrk%vx2l,& - & info,work=work,& - & vtx=wv(1),vty=p%precv(level+1)%wrk%wv(1)) + & info,work=work,vtx=wv(1)) if (info /= psb_success_) then call psb_errpush(psb_err_internal_error_,name,& & a_err='Error during restriction') @@ -547,8 +547,7 @@ contains ! call p%precv(level+1)%map_prol(done,& & p%precv(level+1)%wrk%vy2l, done,vy2l,& - & info,work=work,& - & vtx=p%precv(level+1)%wrk%wv(1),vty=wv(1)) + & info,work=work, vty=wv(1)) if (info /= psb_success_) then call psb_errpush(psb_err_internal_error_,name,& & a_err='Error during prolongation') @@ -595,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') @@ -652,23 +651,23 @@ contains ! if (pre) then - 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 + 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 + + 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),vty=p%precv(level+1)%wrk%wv(1)) + & info,work=work,vtx=wv(1)) if (info /= psb_success_) then call psb_errpush(psb_err_internal_error_,name,& & a_err='Error during restriction') @@ -678,8 +677,7 @@ contains ! Shortcut: just transfer x2l. call p%precv(level+1)%map_rstr(done,vx2l,& & dzero,p%precv(level+1)%wrk%vx2l,& - & info,work=work,& - & vtx=wv(1),vty=p%precv(level+1)%wrk%wv(1)) + & info,work=work,vtx=wv(1)) if (info /= psb_success_) then call psb_errpush(psb_err_internal_error_,name,& & a_err='Error during restriction') @@ -694,8 +692,7 @@ contains ! call p%precv(level+1)%map_prol(done,& & p%precv(level+1)%wrk%vy2l,done,vy2l,& - & info,work=work,& - & vtx=p%precv(level+1)%wrk%wv(1),vty=wv(1)) + & info,work=work,vty=wv(1)) if (info /= psb_success_) then call psb_errpush(psb_err_internal_error_,name,& & a_err='Error during prolongation') @@ -714,7 +711,7 @@ contains if (info == psb_success_) & & call p%precv(level+1)%map_rstr(done,vty,& & dzero,p%precv(level+1)%wrk%vx2l,info,work=work,& - & vtx=wv(1),vty=p%precv(level+1)%wrk%wv(1)) + & vtx=wv(1)) if (info /= psb_success_) then call psb_errpush(psb_err_internal_error_,name,& & a_err='Error during W-cycle restriction') @@ -725,8 +722,7 @@ contains if (info == psb_success_) call p%precv(level+1)%map_prol(done, & & p%precv(level+1)%wrk%vy2l,done,vy2l,& - & info,work=work,& - & vtx=p%precv(level+1)%wrk%wv(1),vty=wv(1)) + & info,work=work,vty=wv(1)) if (info /= psb_success_) then call psb_errpush(psb_err_internal_error_,name,& @@ -837,7 +833,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') @@ -914,7 +910,7 @@ contains call p%precv(level + 1)%map_rstr(done,vty,& & dzero,p%precv(level + 1)%wrk%vx2l,& &info,work=work,& - & vtx=wv(1),vty=p%precv(level+1)%wrk%wv(1)) + & vtx=wv(1)) if (info /= psb_success_) then call psb_errpush(psb_err_internal_error_,name,& @@ -949,8 +945,7 @@ contains ! call p%precv(level+1)%map_prol(done,& & p%precv(level+1)%wrk%vy2l,done,vy2l,& - & info,work=work,& - & vtx=p%precv(level+1)%wrk%wv(1),vty=wv(1)) + & info,work=work,vty=wv(1)) if (info /= psb_success_) then call psb_errpush(psb_err_internal_error_,name,& diff --git a/amgprec/impl/amg_dmlprec_aply_a.f90 b/amgprec/impl/amg_dmlprec_aply_a.f90 index 832656b5..5ec271c5 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 ', size(p%precv) + & ' Entry ', p%get_nlevs() trans_ = psb_toupper(trans) - nlev = size(p%precv) + nlev = p%get_nlevs() 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 = 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') @@ -445,7 +445,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') @@ -550,7 +550,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') diff --git a/amgprec/impl/amg_dprecinit.F90 b/amgprec/impl/amg_dprecinit.F90 index f8c75f1d..62ad8d96 100644 --- a/amgprec/impl/amg_dprecinit.F90 +++ b/amgprec/impl/amg_dprecinit.F90 @@ -255,8 +255,8 @@ subroutine amg_dprecinit(ctxt,prec,ptype,info) write(psb_err_unit,*) name,& &': Warning: Unknown preconditioner type request "',ptype,'"' info = psb_err_pivot_too_small_ - end select + call prec%set_nlevs(nlev_) call psb_erractionrestore(err_act) return diff --git a/amgprec/impl/amg_smlprec_aply.f90 b/amgprec/impl/amg_smlprec_aply.f90 index 40091e60..439657bf 100644 --- a/amgprec/impl/amg_smlprec_aply.f90 +++ b/amgprec/impl/amg_smlprec_aply.f90 @@ -527,8 +527,7 @@ contains ! Apply the restriction call p%precv(level+1)%map_rstr(sone,vx2l,& & szero,p%precv(level+1)%wrk%vx2l,& - & info,work=work,& - & vtx=wv(1),vty=p%precv(level+1)%wrk%wv(1)) + & info,work=work,vtx=wv(1)) if (info /= psb_success_) then call psb_errpush(psb_err_internal_error_,name,& & a_err='Error during restriction') @@ -547,8 +546,7 @@ contains ! call p%precv(level+1)%map_prol(sone,& & p%precv(level+1)%wrk%vy2l, sone,vy2l,& - & info,work=work,& - & vtx=p%precv(level+1)%wrk%wv(1),vty=wv(1)) + & info,work=work, vty=wv(1)) if (info /= psb_success_) then call psb_errpush(psb_err_internal_error_,name,& & a_err='Error during prolongation') @@ -667,8 +665,7 @@ contains call p%precv(level+1)%map_rstr(sone,vty,& & szero,p%precv(level+1)%wrk%vx2l,& - & info,work=work,& - & vtx=wv(1),vty=p%precv(level+1)%wrk%wv(1)) + & info,work=work,vtx=wv(1)) if (info /= psb_success_) then call psb_errpush(psb_err_internal_error_,name,& & a_err='Error during restriction') @@ -678,8 +675,7 @@ contains ! Shortcut: just transfer x2l. call p%precv(level+1)%map_rstr(sone,vx2l,& & szero,p%precv(level+1)%wrk%vx2l,& - & info,work=work,& - & vtx=wv(1),vty=p%precv(level+1)%wrk%wv(1)) + & info,work=work,vtx=wv(1)) if (info /= psb_success_) then call psb_errpush(psb_err_internal_error_,name,& & a_err='Error during restriction') @@ -694,8 +690,7 @@ contains ! call p%precv(level+1)%map_prol(sone,& & p%precv(level+1)%wrk%vy2l,sone,vy2l,& - & info,work=work,& - & vtx=p%precv(level+1)%wrk%wv(1),vty=wv(1)) + & info,work=work,vty=wv(1)) if (info /= psb_success_) then call psb_errpush(psb_err_internal_error_,name,& & a_err='Error during prolongation') @@ -714,7 +709,7 @@ contains if (info == psb_success_) & & call p%precv(level+1)%map_rstr(sone,vty,& & szero,p%precv(level+1)%wrk%vx2l,info,work=work,& - & vtx=wv(1),vty=p%precv(level+1)%wrk%wv(1)) + & vtx=wv(1)) if (info /= psb_success_) then call psb_errpush(psb_err_internal_error_,name,& & a_err='Error during W-cycle restriction') @@ -725,8 +720,7 @@ contains if (info == psb_success_) call p%precv(level+1)%map_prol(sone, & & p%precv(level+1)%wrk%vy2l,sone,vy2l,& - & info,work=work,& - & vtx=p%precv(level+1)%wrk%wv(1),vty=wv(1)) + & info,work=work,vty=wv(1)) if (info /= psb_success_) then call psb_errpush(psb_err_internal_error_,name,& @@ -914,7 +908,7 @@ contains call p%precv(level + 1)%map_rstr(sone,vty,& & szero,p%precv(level + 1)%wrk%vx2l,& &info,work=work,& - & vtx=wv(1),vty=p%precv(level+1)%wrk%wv(1)) + & vtx=wv(1)) if (info /= psb_success_) then call psb_errpush(psb_err_internal_error_,name,& @@ -949,8 +943,7 @@ contains ! call p%precv(level+1)%map_prol(sone,& & p%precv(level+1)%wrk%vy2l,sone,vy2l,& - & info,work=work,& - & vtx=p%precv(level+1)%wrk%wv(1),vty=wv(1)) + & info,work=work,vty=wv(1)) if (info /= psb_success_) then call psb_errpush(psb_err_internal_error_,name,& diff --git a/amgprec/impl/amg_zmlprec_aply.f90 b/amgprec/impl/amg_zmlprec_aply.f90 index 8b22eab5..86e8850e 100644 --- a/amgprec/impl/amg_zmlprec_aply.f90 +++ b/amgprec/impl/amg_zmlprec_aply.f90 @@ -527,8 +527,7 @@ contains ! Apply the restriction call p%precv(level+1)%map_rstr(zone,vx2l,& & zzero,p%precv(level+1)%wrk%vx2l,& - & info,work=work,& - & vtx=wv(1),vty=p%precv(level+1)%wrk%wv(1)) + & info,work=work,vtx=wv(1)) if (info /= psb_success_) then call psb_errpush(psb_err_internal_error_,name,& & a_err='Error during restriction') @@ -547,8 +546,7 @@ contains ! call p%precv(level+1)%map_prol(zone,& & p%precv(level+1)%wrk%vy2l, zone,vy2l,& - & info,work=work,& - & vtx=p%precv(level+1)%wrk%wv(1),vty=wv(1)) + & info,work=work, vty=wv(1)) if (info /= psb_success_) then call psb_errpush(psb_err_internal_error_,name,& & a_err='Error during prolongation') @@ -667,8 +665,7 @@ contains call p%precv(level+1)%map_rstr(zone,vty,& & zzero,p%precv(level+1)%wrk%vx2l,& - & info,work=work,& - & vtx=wv(1),vty=p%precv(level+1)%wrk%wv(1)) + & info,work=work,vtx=wv(1)) if (info /= psb_success_) then call psb_errpush(psb_err_internal_error_,name,& & a_err='Error during restriction') @@ -678,8 +675,7 @@ contains ! Shortcut: just transfer x2l. call p%precv(level+1)%map_rstr(zone,vx2l,& & zzero,p%precv(level+1)%wrk%vx2l,& - & info,work=work,& - & vtx=wv(1),vty=p%precv(level+1)%wrk%wv(1)) + & info,work=work,vtx=wv(1)) if (info /= psb_success_) then call psb_errpush(psb_err_internal_error_,name,& & a_err='Error during restriction') @@ -694,8 +690,7 @@ contains ! call p%precv(level+1)%map_prol(zone,& & p%precv(level+1)%wrk%vy2l,zone,vy2l,& - & info,work=work,& - & vtx=p%precv(level+1)%wrk%wv(1),vty=wv(1)) + & info,work=work,vty=wv(1)) if (info /= psb_success_) then call psb_errpush(psb_err_internal_error_,name,& & a_err='Error during prolongation') @@ -714,7 +709,7 @@ contains if (info == psb_success_) & & call p%precv(level+1)%map_rstr(zone,vty,& & zzero,p%precv(level+1)%wrk%vx2l,info,work=work,& - & vtx=wv(1),vty=p%precv(level+1)%wrk%wv(1)) + & vtx=wv(1)) if (info /= psb_success_) then call psb_errpush(psb_err_internal_error_,name,& & a_err='Error during W-cycle restriction') @@ -725,8 +720,7 @@ contains if (info == psb_success_) call p%precv(level+1)%map_prol(zone, & & p%precv(level+1)%wrk%vy2l,zone,vy2l,& - & info,work=work,& - & vtx=p%precv(level+1)%wrk%wv(1),vty=wv(1)) + & info,work=work,vty=wv(1)) if (info /= psb_success_) then call psb_errpush(psb_err_internal_error_,name,& @@ -914,7 +908,7 @@ contains call p%precv(level + 1)%map_rstr(zone,vty,& & zzero,p%precv(level + 1)%wrk%vx2l,& &info,work=work,& - & vtx=wv(1),vty=p%precv(level+1)%wrk%wv(1)) + & vtx=wv(1)) if (info /= psb_success_) then call psb_errpush(psb_err_internal_error_,name,& @@ -949,8 +943,7 @@ contains ! call p%precv(level+1)%map_prol(zone,& & p%precv(level+1)%wrk%vy2l,zone,vy2l,& - & info,work=work,& - & vtx=p%precv(level+1)%wrk%wv(1),vty=wv(1)) + & info,work=work,vty=wv(1)) if (info /= psb_success_) then call psb_errpush(psb_err_internal_error_,name,& 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 08ee412f..99d46096 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,14 @@ 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() + if (present(vtx)) then + vtx_ => vtx + else + vtx_ => lv%wrk%wv(1) + end if !!$ write(0,*) 'New map_prol',lv%remap_data%ac_pre_remap%is_asb() if (lv%remap_data%ac_pre_remap%is_asb()) then @@ -98,7 +106,7 @@ subroutine amg_c_base_onelev_map_prol_v(lv,alpha,vect_v,beta,vect_u,info,work,vt call psb_rcv(ctxt,tv%v%v(1:nrl),idest) call tv%set_host() call lv%linmap%map_V2U(alpha,tv,beta,vect_u,info,& - & work=work,vtx=vtx,vty=vty) + & work=work,vtx=vtx_,vty=vty) end associate !!$ write(0,*) me, ' Prolongator with remap done ' !!$ flush(0) @@ -107,7 +115,7 @@ subroutine amg_c_base_onelev_map_prol_v(lv,alpha,vect_v,beta,vect_u,info,work,vt else ! Default transfer call lv%linmap%map_V2U(alpha,vect_v,beta,vect_u,info,& - & work=work,vtx=vtx,vty=vty) + & work=work,vtx=vtx_,vty=vty) end if end subroutine amg_c_base_onelev_map_prol_v 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 ad3b81e7..c49e8f50 100644 --- a/amgprec/impl/level/amg_c_base_onelev_map_rstr.F90 +++ b/amgprec/impl/level/amg_c_base_onelev_map_rstr.F90 @@ -47,8 +47,15 @@ subroutine amg_c_base_onelev_map_rstr_v(lv,alpha,vect_u,beta,vect_v,info,& 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 :: vty_ !!$ write(0,*) 'New map_rstr',lv%remap_data%ac_pre_remap%is_asb() + if (present(vty)) then + vty_ => vty + else + vty_ => lv%wrk%wv(1) + end if + if (lv%remap_data%ac_pre_remap%is_asb()) then ! ! Remap has happened, deal with it @@ -74,11 +81,11 @@ 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,' map_rstr calling U2V: ',me,np,rme,rnp,tv%get_nrows() + write(0,*) me,' remap map_rstr calling U2V: ',me,np,rme,rnp,tv%get_nrows() flush(0) call psb_barrier(ctxt) call lv%linmap%map_U2V(alpha,vect_u,beta,tv,info,& - & work=work,vtx=vtx,vty=vty) + & work=work,vtx=vtx,vty=vty_) call tv%sync() !rsnd = tv%get_vect() !call psb_snd(ctxt,rsnd(1:nrl),idest) @@ -103,8 +110,15 @@ subroutine amg_c_base_onelev_map_rstr_v(lv,alpha,vect_u,beta,vect_v,info,& else ! Default transfer - call lv%linmap%map_U2V(alpha,vect_u,beta,vect_v,info,& - & work=work,vtx=vtx,vty=vty) + 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 + call lv%linmap%map_U2V(alpha,vect_u,beta,vect_v,info,& + & work=work,vtx=vtx,vty=vty_) + end block end if end subroutine amg_c_base_onelev_map_rstr_v diff --git a/amgprec/impl/level/amg_d_base_onelev_map_prol.F90 b/amgprec/impl/level/amg_d_base_onelev_map_prol.F90 index dd1b0a9a..4dc7ea84 100644 --- a/amgprec/impl/level/amg_d_base_onelev_map_prol.F90 +++ b/amgprec/impl/level/amg_d_base_onelev_map_prol.F90 @@ -46,6 +46,15 @@ subroutine amg_d_base_onelev_map_prol_v(lv,alpha,vect_v,beta,vect_u,info,work,vt integer(psb_ipk_), intent(out) :: info real(psb_dpk_), optional :: work(:) type(psb_d_vect_type), optional, target, intent(inout) :: vtx,vty + + type(psb_d_vect_type), pointer :: vtx_ + +!!$ write(0,*) 'New map_rstr',lv%remap_data%ac_pre_remap%is_asb() + if (present(vtx)) then + vtx_ => vtx + else + vtx_ => lv%wrk%wv(1) + end if !!$ write(0,*) 'New map_prol',lv%remap_data%ac_pre_remap%is_asb() if (lv%remap_data%ac_pre_remap%is_asb()) then @@ -98,7 +107,7 @@ subroutine amg_d_base_onelev_map_prol_v(lv,alpha,vect_v,beta,vect_u,info,work,vt call psb_rcv(ctxt,tv%v%v(1:nrl),idest) call tv%set_host() call lv%linmap%map_V2U(alpha,tv,beta,vect_u,info,& - & work=work,vtx=vtx,vty=vty) + & work=work,vtx=vtx_,vty=vty) end associate !!$ write(0,*) me, ' Prolongator with remap done ' !!$ flush(0) @@ -107,7 +116,7 @@ subroutine amg_d_base_onelev_map_prol_v(lv,alpha,vect_v,beta,vect_u,info,work,vt else ! Default transfer call lv%linmap%map_V2U(alpha,vect_v,beta,vect_u,info,& - & work=work,vtx=vtx,vty=vty) + & work=work,vtx=vtx_,vty=vty) end if end subroutine amg_d_base_onelev_map_prol_v 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 ee17de60..dac4b3de 100644 --- a/amgprec/impl/level/amg_d_base_onelev_map_rstr.F90 +++ b/amgprec/impl/level/amg_d_base_onelev_map_rstr.F90 @@ -35,7 +35,6 @@ ! 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,17 +46,25 @@ subroutine amg_d_base_onelev_map_rstr_v(lv,alpha,vect_u,beta,vect_v,info,& integer(psb_ipk_), intent(out) :: 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 -!!$ 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 + vty_ => lv%wrk%wv(1) + end if + if (lv%remap_data%ac_pre_remap%is_asb()) then ! ! Remap has happened, deal with it ! -!!$ write(0,*) 'Remap handling not implemented yet ' + write(0,*) 'Level map_rstr with remapping ' 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_dpk_), allocatable :: rsnd(:), rrcv(:) type(psb_d_vect_type) :: tv @@ -65,24 +72,25 @@ subroutine amg_d_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 -!!$ if (rme >= 0) write(0,*) rme, ' Receiving data from ',isrc(:) + write(0,*) 'Should apply maps, then send data from ',me,' to ',idest + if (rme >= 0) write(0,*) rme, ' Receiving data from ',isrc(:) nsrc = size(isrc) 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,' map_rstr calling U2V: ',me,np,rme,rnp,tv%get_nrows() + 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() + & 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))) @@ -91,22 +99,28 @@ subroutine amg_d_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 - call lv%linmap%map_U2V(alpha,vect_u,beta,vect_v,info,& - & work=work,vtx=vtx,vty=vty) + block + type(psb_ctxt_type) :: ctxt, rctxt + ctxt = lv%linmap%p_desc_U%get_ctxt() + call psb_info(ctxt,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_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_s_base_onelev_map_prol.F90 b/amgprec/impl/level/amg_s_base_onelev_map_prol.F90 index 97ef0988..0ceaf150 100644 --- a/amgprec/impl/level/amg_s_base_onelev_map_prol.F90 +++ b/amgprec/impl/level/amg_s_base_onelev_map_prol.F90 @@ -46,6 +46,15 @@ subroutine amg_s_base_onelev_map_prol_v(lv,alpha,vect_v,beta,vect_u,info,work,vt integer(psb_ipk_), intent(out) :: info real(psb_spk_), optional :: work(:) type(psb_s_vect_type), optional, target, intent(inout) :: vtx,vty + + type(psb_s_vect_type), pointer :: vtx_ + +!!$ write(0,*) 'New map_rstr',lv%remap_data%ac_pre_remap%is_asb() + if (present(vtx)) then + vtx_ => vtx + else + vtx_ => lv%wrk%wv(1) + end if !!$ write(0,*) 'New map_prol',lv%remap_data%ac_pre_remap%is_asb() if (lv%remap_data%ac_pre_remap%is_asb()) then @@ -98,7 +107,7 @@ subroutine amg_s_base_onelev_map_prol_v(lv,alpha,vect_v,beta,vect_u,info,work,vt call psb_rcv(ctxt,tv%v%v(1:nrl),idest) call tv%set_host() call lv%linmap%map_V2U(alpha,tv,beta,vect_u,info,& - & work=work,vtx=vtx,vty=vty) + & work=work,vtx=vtx_,vty=vty) end associate !!$ write(0,*) me, ' Prolongator with remap done ' !!$ flush(0) @@ -107,7 +116,7 @@ subroutine amg_s_base_onelev_map_prol_v(lv,alpha,vect_v,beta,vect_u,info,work,vt else ! Default transfer call lv%linmap%map_V2U(alpha,vect_v,beta,vect_u,info,& - & work=work,vtx=vtx,vty=vty) + & work=work,vtx=vtx_,vty=vty) end if end subroutine amg_s_base_onelev_map_prol_v 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 e9e0cdb5..e81d5b31 100644 --- a/amgprec/impl/level/amg_s_base_onelev_map_rstr.F90 +++ b/amgprec/impl/level/amg_s_base_onelev_map_rstr.F90 @@ -47,8 +47,15 @@ subroutine amg_s_base_onelev_map_rstr_v(lv,alpha,vect_u,beta,vect_v,info,& integer(psb_ipk_), intent(out) :: info real(psb_spk_), optional :: work(:) type(psb_s_vect_type), optional, target, intent(inout) :: vtx,vty + type(psb_s_vect_type), pointer :: vty_ !!$ write(0,*) 'New map_rstr',lv%remap_data%ac_pre_remap%is_asb() + if (present(vty)) then + vty_ => vty + else + vty_ => lv%wrk%wv(1) + end if + if (lv%remap_data%ac_pre_remap%is_asb()) then ! ! Remap has happened, deal with it @@ -74,11 +81,11 @@ 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,' map_rstr calling U2V: ',me,np,rme,rnp,tv%get_nrows() + write(0,*) me,' remap map_rstr calling U2V: ',me,np,rme,rnp,tv%get_nrows() flush(0) call psb_barrier(ctxt) call lv%linmap%map_U2V(alpha,vect_u,beta,tv,info,& - & work=work,vtx=vtx,vty=vty) + & work=work,vtx=vtx,vty=vty_) call tv%sync() !rsnd = tv%get_vect() !call psb_snd(ctxt,rsnd(1:nrl),idest) @@ -103,8 +110,15 @@ subroutine amg_s_base_onelev_map_rstr_v(lv,alpha,vect_u,beta,vect_v,info,& else ! Default transfer - call lv%linmap%map_U2V(alpha,vect_u,beta,vect_v,info,& - & work=work,vtx=vtx,vty=vty) + 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 + call lv%linmap%map_U2V(alpha,vect_u,beta,vect_v,info,& + & work=work,vtx=vtx,vty=vty_) + end block end if end subroutine amg_s_base_onelev_map_rstr_v diff --git a/amgprec/impl/level/amg_z_base_onelev_map_prol.F90 b/amgprec/impl/level/amg_z_base_onelev_map_prol.F90 index 37c356ac..0a7d8029 100644 --- a/amgprec/impl/level/amg_z_base_onelev_map_prol.F90 +++ b/amgprec/impl/level/amg_z_base_onelev_map_prol.F90 @@ -46,6 +46,15 @@ subroutine amg_z_base_onelev_map_prol_v(lv,alpha,vect_v,beta,vect_u,info,work,vt integer(psb_ipk_), intent(out) :: info complex(psb_dpk_), optional :: work(:) type(psb_z_vect_type), optional, target, intent(inout) :: vtx,vty + + type(psb_z_vect_type), pointer :: vtx_ + +!!$ write(0,*) 'New map_rstr',lv%remap_data%ac_pre_remap%is_asb() + if (present(vtx)) then + vtx_ => vtx + else + vtx_ => lv%wrk%wv(1) + end if !!$ write(0,*) 'New map_prol',lv%remap_data%ac_pre_remap%is_asb() if (lv%remap_data%ac_pre_remap%is_asb()) then @@ -98,7 +107,7 @@ subroutine amg_z_base_onelev_map_prol_v(lv,alpha,vect_v,beta,vect_u,info,work,vt call psb_rcv(ctxt,tv%v%v(1:nrl),idest) call tv%set_host() call lv%linmap%map_V2U(alpha,tv,beta,vect_u,info,& - & work=work,vtx=vtx,vty=vty) + & work=work,vtx=vtx_,vty=vty) end associate !!$ write(0,*) me, ' Prolongator with remap done ' !!$ flush(0) @@ -107,7 +116,7 @@ subroutine amg_z_base_onelev_map_prol_v(lv,alpha,vect_v,beta,vect_u,info,work,vt else ! Default transfer call lv%linmap%map_V2U(alpha,vect_v,beta,vect_u,info,& - & work=work,vtx=vtx,vty=vty) + & work=work,vtx=vtx_,vty=vty) end if end subroutine amg_z_base_onelev_map_prol_v 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 c3a742dc..ffcd9a72 100644 --- a/amgprec/impl/level/amg_z_base_onelev_map_rstr.F90 +++ b/amgprec/impl/level/amg_z_base_onelev_map_rstr.F90 @@ -47,8 +47,15 @@ subroutine amg_z_base_onelev_map_rstr_v(lv,alpha,vect_u,beta,vect_v,info,& integer(psb_ipk_), intent(out) :: info complex(psb_dpk_), optional :: work(:) type(psb_z_vect_type), optional, target, intent(inout) :: vtx,vty + type(psb_z_vect_type), pointer :: vty_ !!$ write(0,*) 'New map_rstr',lv%remap_data%ac_pre_remap%is_asb() + if (present(vty)) then + vty_ => vty + else + vty_ => lv%wrk%wv(1) + end if + if (lv%remap_data%ac_pre_remap%is_asb()) then ! ! Remap has happened, deal with it @@ -74,11 +81,11 @@ 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,' map_rstr calling U2V: ',me,np,rme,rnp,tv%get_nrows() + write(0,*) me,' remap map_rstr calling U2V: ',me,np,rme,rnp,tv%get_nrows() flush(0) call psb_barrier(ctxt) call lv%linmap%map_U2V(alpha,vect_u,beta,tv,info,& - & work=work,vtx=vtx,vty=vty) + & work=work,vtx=vtx,vty=vty_) call tv%sync() !rsnd = tv%get_vect() !call psb_snd(ctxt,rsnd(1:nrl),idest) @@ -103,8 +110,15 @@ subroutine amg_z_base_onelev_map_rstr_v(lv,alpha,vect_u,beta,vect_v,info,& else ! Default transfer - call lv%linmap%map_U2V(alpha,vect_u,beta,vect_v,info,& - & work=work,vtx=vtx,vty=vty) + 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 + call lv%linmap%map_U2V(alpha,vect_u,beta,vect_v,info,& + & work=work,vtx=vtx,vty=vty_) + end block end if end subroutine amg_z_base_onelev_map_rstr_v diff --git a/samples/advanced/remap/amg_d_pde3d.F90 b/samples/advanced/remap/amg_d_pde3d.F90 index c092bf08..117e338f 100644 --- a/samples/advanced/remap/amg_d_pde3d.F90 +++ b/samples/advanced/remap/amg_d_pde3d.F90 @@ -472,6 +472,12 @@ program amg_d_pde3d call psb_errpush(psb_err_from_subroutine_,name,a_err='amg_hierarchy_bld') goto 9999 end if + + do i=2, size(prec%precv) + write(0,*) iam,'Between hier and smoothers_bld level',i,':',& + & prec%precv(i)%remap_data%desc_ac_pre_remap%is_asb() + end do + call psb_barrier(ctxt) t1 = psb_wtime() call prec%smoothers_build(a,desc_a,info) @@ -491,6 +497,12 @@ program amg_d_pde3d write(psb_out_unit,'(" ")') end if + do i=2, size(prec%precv) + write(0,*) iam,'After smoothers_bld level',i,':',& + & prec%precv(i)%remap_data%desc_ac_pre_remap%is_asb() + end do + + call prec%descr(info,iout=psb_out_unit) if (p_choice%dump) then call prec%dump(info,istart=p_choice%dlmin,iend=p_choice%dlmax,&