diff --git a/mlprec/impl/mld_c_hierarchy_bld.f90 b/mlprec/impl/mld_c_hierarchy_bld.f90 index bad8ba82..836a6905 100644 --- a/mlprec/impl/mld_c_hierarchy_bld.f90 +++ b/mlprec/impl/mld_c_hierarchy_bld.f90 @@ -82,7 +82,7 @@ subroutine mld_c_hierarchy_bld(a,desc_a,p,info) integer(psb_ipk_) :: ictxt, me,np integer(psb_ipk_) :: err,i,k, err_act, iszv, newsz, casize, nplevs, mxplevs, iaggsize real(psb_spk_) :: mnaggratio, sizeratio - class(mld_c_base_smoother_type), allocatable :: coarse_sm, base_sm, med_sm + class(mld_c_base_smoother_type), allocatable :: coarse_sm, base_sm, med_sm, base_sm2, med_sm2, coarse_sm2 type(mld_sml_parms) :: baseparms, medparms, coarseparms integer(psb_ipk_), allocatable :: ilaggr(:), nlaggr(:) type(psb_cspmat_type) :: op_prol @@ -215,50 +215,40 @@ subroutine mld_c_hierarchy_bld(a,desc_a,p,info) baseparms = p%precv(1)%parms medparms = p%precv(2)%parms - allocate(coarse_sm, source=p%precv(iszv)%sm,stat=info) - if (info == psb_success_) & - & allocate(med_sm, source=p%precv(2)%sm,stat=info) - if (info == psb_success_) & - & allocate(base_sm, source=p%precv(1)%sm,stat=info) + call save_smoothers(p%precv(iszv),coarse_sm,coarse_sm2,info) + if (info == 0) call save_smoothers(p%precv(2),med_sm,med_sm2,info) + if (info == 0) call save_smoothers(p%precv(1),base_sm,base_sm2,info) if (info /= psb_success_) then write(0,*) 'Error in saving smoothers',info call psb_errpush(psb_err_internal_error_,name,a_err='Base level precbuild.') goto 9999 end if - ! ! First set desired number of levels ! if (iszv /= nplevs) then allocate(tprecv(nplevs),stat=info) - if (info /= psb_success_) then - call psb_errpush(psb_err_from_subroutine_,name,& - & a_err='prec reallocation') - goto 9999 - endif - tprecv(1)%parms = baseparms - allocate(tprecv(1)%sm,source=base_sm,stat=info) - if (info /= psb_success_) then - call psb_errpush(psb_err_from_subroutine_,name,& - & a_err='prec reallocation') - goto 9999 - endif - do i=2,nplevs-1 - tprecv(i)%parms = medparms - allocate(tprecv(i)%sm,source=med_sm,stat=info) - if (info /= psb_success_) then - call psb_errpush(psb_err_from_subroutine_,name,& - & a_err='prec reallocation') - goto 9999 - endif + ! First all existing levels + if (info == 0) tprecv(1)%parms = baseparms + if (info == 0) call restore_smoothers(tprecv(1),p%precv(1)%sm,p%precv(1)%sm2a,info) + do i=2, min(iszv,nplevs) - 1 + if (info == 0) tprecv(i)%parms = medparms + if (info == 0) call restore_smoothers(tprecv(i),p%precv(i)%sm,p%precv(i)%sm2a,info) + end do + ! Further intermediates, if any + do i=iszv-1, nplevs - 1 + if (info == 0) tprecv(i)%parms = medparms + if (info == 0) call restore_smoothers(tprecv(i),med_sm,med_sm2,info) end do - tprecv(nplevs)%parms = coarseparms - allocate(tprecv(nplevs)%sm,source=coarse_sm,stat=info) + ! Then coarse + if (info == 0) tprecv(nplevs)%parms = coarseparms + if (info == 0) call restore_smoothers(tprecv(nplevs),coarse_sm,coarse_sm2,info) if (info /= psb_success_) then call psb_errpush(psb_err_from_subroutine_,name,& & a_err='prec reallocation') goto 9999 endif + do i=1,iszv call p%precv(i)%free(info) end do @@ -330,7 +320,7 @@ subroutine mld_c_hierarchy_bld(a,desc_a,p,info) ! ! We are not gaining ! - newsz = newsz-1 + newsz = i-1 end if end if @@ -350,20 +340,25 @@ subroutine mld_c_hierarchy_bld(a,desc_a,p,info) end if end if call psb_bcast(ictxt,newsz) - if (newsz > 0) & - & call p%precv(i)%parms%get_coarse(p%precv(iszv)%parms) - - if (info == psb_success_) call mld_lev_mat_asb(p%precv(i),& - & p%precv(i-1)%base_a,p%precv(i-1)%base_desc,& - & ilaggr,nlaggr,op_prol,info) - + if (newsz > 0) then + if (info == 0) p%precv(newsz)%parms = coarseparms + if (info == 0) call restore_smoothers(p%precv(newsz),coarse_sm,coarse_sm2,info) + + if (info == psb_success_) call mld_lev_mat_asb(p%precv(newsz),& + & p%precv(newsz-1)%base_a,p%precv(newsz-1)%base_desc,& + & ilaggr,nlaggr,op_prol,info) + exit array_build_loop + else + if (info == psb_success_) call mld_lev_mat_asb(p%precv(i),& + & p%precv(i-1)%base_a,p%precv(i-1)%base_desc,& + & ilaggr,nlaggr,op_prol,info) + end if if (info /= psb_success_) then call psb_errpush(psb_err_internal_error_,name,& & a_err='Map build') goto 9999 endif - if (newsz > 0) exit array_build_loop end do array_build_loop if (newsz > 0) then @@ -399,7 +394,6 @@ subroutine mld_c_hierarchy_bld(a,desc_a,p,info) end do end if - if (info /= psb_success_) then call psb_errpush(psb_err_internal_error_,name,& & a_err='Internal hierarchy build' ) @@ -419,4 +413,58 @@ subroutine mld_c_hierarchy_bld(a,desc_a,p,info) return +contains + subroutine save_smoothers(level,save1, save2,info) + type(mld_c_onelev_type), intent(in) :: level + class(mld_c_base_smoother_type), allocatable , intent(inout) :: save1, save2 + integer(psb_ipk_), intent(out) :: info + + info = 0 + if (allocated(save1)) then + call save1%free(info) + if (info == 0) deallocate(save1,stat=info) + if (info /= 0) return + end if + if (allocated(save2)) then + call save2%free(info) + if (info == 0) deallocate(save2,stat=info) + if (info /= 0) return + end if + allocate(save1, source=level%sm,stat=info) + if ((info == 0).and.allocated(level%sm2a)) allocate(save2, source=level%sm2a,stat=info) + + return + end subroutine save_smoothers + + subroutine restore_smoothers(level,save1, save2,info) + type(mld_c_onelev_type), intent(inout), target :: level + class(mld_c_base_smoother_type), allocatable, intent(in) :: save1, save2 + integer(psb_ipk_), intent(out) :: info + + info = 0 + + if (allocated(level%sm)) then + if (info == 0) call level%sm%free(info) + if (info == 0) deallocate(level%sm,stat=info) + end if + if (allocated(save1)) then + if (info == 0) allocate(level%sm,source=save1,stat=info) + end if + + if (info /= 0) return + + if (allocated(level%sm2a)) then + if (info == 0) call level%sm2a%free(info) + if (info == 0) deallocate(level%sm2a,stat=info) + end if + if (allocated(save2)) then + if (info == 0) allocate(level%sm2a,source=save2,stat=info) + if (info == 0) level%sm2 => level%sm2a + else + if (allocated(level%sm)) level%sm2 => level%sm + end if + + return + end subroutine restore_smoothers + end subroutine mld_c_hierarchy_bld diff --git a/mlprec/impl/mld_ccprecset.F90 b/mlprec/impl/mld_ccprecset.F90 index 783b4582..8c39848f 100644 --- a/mlprec/impl/mld_ccprecset.F90 +++ b/mlprec/impl/mld_ccprecset.F90 @@ -519,6 +519,11 @@ subroutine mld_ccprecsetr(p,what,val,info,ilev,pos) ilev_=nlev_ call p%precv(ilev_)%set('SUB_ILUTHRS',val,info,pos=pos) + case('AGGR_SCALE') + do ilev_ = 2, nlev_ + call p%precv(ilev_)%set('AGGR_SCALE',val,info,pos=pos) + end do + case('AGGR_THRESH') thr = val do ilev_ = 2, nlev_ diff --git a/mlprec/impl/mld_cprecset.F90 b/mlprec/impl/mld_cprecset.F90 index 38c1d791..02f32538 100644 --- a/mlprec/impl/mld_cprecset.F90 +++ b/mlprec/impl/mld_cprecset.F90 @@ -621,6 +621,11 @@ subroutine mld_cprecsetr(p,what,val,info,ilev,pos) ilev_=nlev_ call p%precv(ilev_)%set(mld_sub_iluthrs_,val,info,pos=pos) + case(mld_aggr_scale_) + do ilev_ = 2, nlev_ + call p%precv(ilev_)%set(mld_aggr_scale_,val,info,pos=pos) + end do + case(mld_aggr_thresh_) thr = val do ilev_ = 2, nlev_ diff --git a/mlprec/impl/mld_d_hierarchy_bld.f90 b/mlprec/impl/mld_d_hierarchy_bld.f90 index e8af7fed..c6061345 100644 --- a/mlprec/impl/mld_d_hierarchy_bld.f90 +++ b/mlprec/impl/mld_d_hierarchy_bld.f90 @@ -82,7 +82,7 @@ subroutine mld_d_hierarchy_bld(a,desc_a,p,info) integer(psb_ipk_) :: ictxt, me,np integer(psb_ipk_) :: err,i,k, err_act, iszv, newsz, casize, nplevs, mxplevs, iaggsize real(psb_dpk_) :: mnaggratio, sizeratio - class(mld_d_base_smoother_type), allocatable :: coarse_sm, base_sm, med_sm + class(mld_d_base_smoother_type), allocatable :: coarse_sm, base_sm, med_sm, base_sm2, med_sm2, coarse_sm2 type(mld_dml_parms) :: baseparms, medparms, coarseparms integer(psb_ipk_), allocatable :: ilaggr(:), nlaggr(:) type(psb_dspmat_type) :: op_prol @@ -215,50 +215,40 @@ subroutine mld_d_hierarchy_bld(a,desc_a,p,info) baseparms = p%precv(1)%parms medparms = p%precv(2)%parms - allocate(coarse_sm, source=p%precv(iszv)%sm,stat=info) - if (info == psb_success_) & - & allocate(med_sm, source=p%precv(2)%sm,stat=info) - if (info == psb_success_) & - & allocate(base_sm, source=p%precv(1)%sm,stat=info) + call save_smoothers(p%precv(iszv),coarse_sm,coarse_sm2,info) + if (info == 0) call save_smoothers(p%precv(2),med_sm,med_sm2,info) + if (info == 0) call save_smoothers(p%precv(1),base_sm,base_sm2,info) if (info /= psb_success_) then write(0,*) 'Error in saving smoothers',info call psb_errpush(psb_err_internal_error_,name,a_err='Base level precbuild.') goto 9999 end if - ! ! First set desired number of levels ! if (iszv /= nplevs) then allocate(tprecv(nplevs),stat=info) - if (info /= psb_success_) then - call psb_errpush(psb_err_from_subroutine_,name,& - & a_err='prec reallocation') - goto 9999 - endif - tprecv(1)%parms = baseparms - allocate(tprecv(1)%sm,source=base_sm,stat=info) - if (info /= psb_success_) then - call psb_errpush(psb_err_from_subroutine_,name,& - & a_err='prec reallocation') - goto 9999 - endif - do i=2,nplevs-1 - tprecv(i)%parms = medparms - allocate(tprecv(i)%sm,source=med_sm,stat=info) - if (info /= psb_success_) then - call psb_errpush(psb_err_from_subroutine_,name,& - & a_err='prec reallocation') - goto 9999 - endif + ! First all existing levels + if (info == 0) tprecv(1)%parms = baseparms + if (info == 0) call restore_smoothers(tprecv(1),p%precv(1)%sm,p%precv(1)%sm2a,info) + do i=2, min(iszv,nplevs) - 1 + if (info == 0) tprecv(i)%parms = medparms + if (info == 0) call restore_smoothers(tprecv(i),p%precv(i)%sm,p%precv(i)%sm2a,info) + end do + ! Further intermediates, if any + do i=iszv-1, nplevs - 1 + if (info == 0) tprecv(i)%parms = medparms + if (info == 0) call restore_smoothers(tprecv(i),med_sm,med_sm2,info) end do - tprecv(nplevs)%parms = coarseparms - allocate(tprecv(nplevs)%sm,source=coarse_sm,stat=info) + ! Then coarse + if (info == 0) tprecv(nplevs)%parms = coarseparms + if (info == 0) call restore_smoothers(tprecv(nplevs),coarse_sm,coarse_sm2,info) if (info /= psb_success_) then call psb_errpush(psb_err_from_subroutine_,name,& & a_err='prec reallocation') goto 9999 endif + do i=1,iszv call p%precv(i)%free(info) end do @@ -330,7 +320,7 @@ subroutine mld_d_hierarchy_bld(a,desc_a,p,info) ! ! We are not gaining ! - newsz = newsz-1 + newsz = i-1 end if end if @@ -350,20 +340,25 @@ subroutine mld_d_hierarchy_bld(a,desc_a,p,info) end if end if call psb_bcast(ictxt,newsz) - if (newsz > 0) & - & call p%precv(i)%parms%get_coarse(p%precv(iszv)%parms) - - if (info == psb_success_) call mld_lev_mat_asb(p%precv(i),& - & p%precv(i-1)%base_a,p%precv(i-1)%base_desc,& - & ilaggr,nlaggr,op_prol,info) - + if (newsz > 0) then + if (info == 0) p%precv(newsz)%parms = coarseparms + if (info == 0) call restore_smoothers(p%precv(newsz),coarse_sm,coarse_sm2,info) + + if (info == psb_success_) call mld_lev_mat_asb(p%precv(newsz),& + & p%precv(newsz-1)%base_a,p%precv(newsz-1)%base_desc,& + & ilaggr,nlaggr,op_prol,info) + exit array_build_loop + else + if (info == psb_success_) call mld_lev_mat_asb(p%precv(i),& + & p%precv(i-1)%base_a,p%precv(i-1)%base_desc,& + & ilaggr,nlaggr,op_prol,info) + end if if (info /= psb_success_) then call psb_errpush(psb_err_internal_error_,name,& & a_err='Map build') goto 9999 endif - if (newsz > 0) exit array_build_loop end do array_build_loop if (newsz > 0) then @@ -399,7 +394,6 @@ subroutine mld_d_hierarchy_bld(a,desc_a,p,info) end do end if - if (info /= psb_success_) then call psb_errpush(psb_err_internal_error_,name,& & a_err='Internal hierarchy build' ) @@ -419,4 +413,58 @@ subroutine mld_d_hierarchy_bld(a,desc_a,p,info) return +contains + subroutine save_smoothers(level,save1, save2,info) + type(mld_d_onelev_type), intent(in) :: level + class(mld_d_base_smoother_type), allocatable , intent(inout) :: save1, save2 + integer(psb_ipk_), intent(out) :: info + + info = 0 + if (allocated(save1)) then + call save1%free(info) + if (info == 0) deallocate(save1,stat=info) + if (info /= 0) return + end if + if (allocated(save2)) then + call save2%free(info) + if (info == 0) deallocate(save2,stat=info) + if (info /= 0) return + end if + allocate(save1, source=level%sm,stat=info) + if ((info == 0).and.allocated(level%sm2a)) allocate(save2, source=level%sm2a,stat=info) + + return + end subroutine save_smoothers + + subroutine restore_smoothers(level,save1, save2,info) + type(mld_d_onelev_type), intent(inout), target :: level + class(mld_d_base_smoother_type), allocatable, intent(in) :: save1, save2 + integer(psb_ipk_), intent(out) :: info + + info = 0 + + if (allocated(level%sm)) then + if (info == 0) call level%sm%free(info) + if (info == 0) deallocate(level%sm,stat=info) + end if + if (allocated(save1)) then + if (info == 0) allocate(level%sm,source=save1,stat=info) + end if + + if (info /= 0) return + + if (allocated(level%sm2a)) then + if (info == 0) call level%sm2a%free(info) + if (info == 0) deallocate(level%sm2a,stat=info) + end if + if (allocated(save2)) then + if (info == 0) allocate(level%sm2a,source=save2,stat=info) + if (info == 0) level%sm2 => level%sm2a + else + if (allocated(level%sm)) level%sm2 => level%sm + end if + + return + end subroutine restore_smoothers + end subroutine mld_d_hierarchy_bld diff --git a/mlprec/impl/mld_dcprecset.F90 b/mlprec/impl/mld_dcprecset.F90 index bb61822a..005ecbba 100644 --- a/mlprec/impl/mld_dcprecset.F90 +++ b/mlprec/impl/mld_dcprecset.F90 @@ -529,6 +529,11 @@ subroutine mld_dcprecsetr(p,what,val,info,ilev,pos) ilev_=nlev_ call p%precv(ilev_)%set('SUB_ILUTHRS',val,info,pos=pos) + case('AGGR_SCALE') + do ilev_ = 2, nlev_ + call p%precv(ilev_)%set('AGGR_SCALE',val,info,pos=pos) + end do + case('AGGR_THRESH') thr = val do ilev_ = 2, nlev_ diff --git a/mlprec/impl/mld_dprecset.F90 b/mlprec/impl/mld_dprecset.F90 index 37236b03..6d868885 100644 --- a/mlprec/impl/mld_dprecset.F90 +++ b/mlprec/impl/mld_dprecset.F90 @@ -631,6 +631,11 @@ subroutine mld_dprecsetr(p,what,val,info,ilev,pos) ilev_=nlev_ call p%precv(ilev_)%set(mld_sub_iluthrs_,val,info,pos=pos) + case(mld_aggr_scale_) + do ilev_ = 2, nlev_ + call p%precv(ilev_)%set(mld_aggr_scale_,val,info,pos=pos) + end do + case(mld_aggr_thresh_) thr = val do ilev_ = 2, nlev_ diff --git a/mlprec/impl/mld_s_hierarchy_bld.f90 b/mlprec/impl/mld_s_hierarchy_bld.f90 index acd8aa53..b6c1f269 100644 --- a/mlprec/impl/mld_s_hierarchy_bld.f90 +++ b/mlprec/impl/mld_s_hierarchy_bld.f90 @@ -82,7 +82,7 @@ subroutine mld_s_hierarchy_bld(a,desc_a,p,info) integer(psb_ipk_) :: ictxt, me,np integer(psb_ipk_) :: err,i,k, err_act, iszv, newsz, casize, nplevs, mxplevs, iaggsize real(psb_spk_) :: mnaggratio, sizeratio - class(mld_s_base_smoother_type), allocatable :: coarse_sm, base_sm, med_sm + class(mld_s_base_smoother_type), allocatable :: coarse_sm, base_sm, med_sm, base_sm2, med_sm2, coarse_sm2 type(mld_sml_parms) :: baseparms, medparms, coarseparms integer(psb_ipk_), allocatable :: ilaggr(:), nlaggr(:) type(psb_sspmat_type) :: op_prol @@ -215,50 +215,40 @@ subroutine mld_s_hierarchy_bld(a,desc_a,p,info) baseparms = p%precv(1)%parms medparms = p%precv(2)%parms - allocate(coarse_sm, source=p%precv(iszv)%sm,stat=info) - if (info == psb_success_) & - & allocate(med_sm, source=p%precv(2)%sm,stat=info) - if (info == psb_success_) & - & allocate(base_sm, source=p%precv(1)%sm,stat=info) + call save_smoothers(p%precv(iszv),coarse_sm,coarse_sm2,info) + if (info == 0) call save_smoothers(p%precv(2),med_sm,med_sm2,info) + if (info == 0) call save_smoothers(p%precv(1),base_sm,base_sm2,info) if (info /= psb_success_) then write(0,*) 'Error in saving smoothers',info call psb_errpush(psb_err_internal_error_,name,a_err='Base level precbuild.') goto 9999 end if - ! ! First set desired number of levels ! if (iszv /= nplevs) then allocate(tprecv(nplevs),stat=info) - if (info /= psb_success_) then - call psb_errpush(psb_err_from_subroutine_,name,& - & a_err='prec reallocation') - goto 9999 - endif - tprecv(1)%parms = baseparms - allocate(tprecv(1)%sm,source=base_sm,stat=info) - if (info /= psb_success_) then - call psb_errpush(psb_err_from_subroutine_,name,& - & a_err='prec reallocation') - goto 9999 - endif - do i=2,nplevs-1 - tprecv(i)%parms = medparms - allocate(tprecv(i)%sm,source=med_sm,stat=info) - if (info /= psb_success_) then - call psb_errpush(psb_err_from_subroutine_,name,& - & a_err='prec reallocation') - goto 9999 - endif + ! First all existing levels + if (info == 0) tprecv(1)%parms = baseparms + if (info == 0) call restore_smoothers(tprecv(1),p%precv(1)%sm,p%precv(1)%sm2a,info) + do i=2, min(iszv,nplevs) - 1 + if (info == 0) tprecv(i)%parms = medparms + if (info == 0) call restore_smoothers(tprecv(i),p%precv(i)%sm,p%precv(i)%sm2a,info) + end do + ! Further intermediates, if any + do i=iszv-1, nplevs - 1 + if (info == 0) tprecv(i)%parms = medparms + if (info == 0) call restore_smoothers(tprecv(i),med_sm,med_sm2,info) end do - tprecv(nplevs)%parms = coarseparms - allocate(tprecv(nplevs)%sm,source=coarse_sm,stat=info) + ! Then coarse + if (info == 0) tprecv(nplevs)%parms = coarseparms + if (info == 0) call restore_smoothers(tprecv(nplevs),coarse_sm,coarse_sm2,info) if (info /= psb_success_) then call psb_errpush(psb_err_from_subroutine_,name,& & a_err='prec reallocation') goto 9999 endif + do i=1,iszv call p%precv(i)%free(info) end do @@ -330,7 +320,7 @@ subroutine mld_s_hierarchy_bld(a,desc_a,p,info) ! ! We are not gaining ! - newsz = newsz-1 + newsz = i-1 end if end if @@ -350,20 +340,25 @@ subroutine mld_s_hierarchy_bld(a,desc_a,p,info) end if end if call psb_bcast(ictxt,newsz) - if (newsz > 0) & - & call p%precv(i)%parms%get_coarse(p%precv(iszv)%parms) - - if (info == psb_success_) call mld_lev_mat_asb(p%precv(i),& - & p%precv(i-1)%base_a,p%precv(i-1)%base_desc,& - & ilaggr,nlaggr,op_prol,info) - + if (newsz > 0) then + if (info == 0) p%precv(newsz)%parms = coarseparms + if (info == 0) call restore_smoothers(p%precv(newsz),coarse_sm,coarse_sm2,info) + + if (info == psb_success_) call mld_lev_mat_asb(p%precv(newsz),& + & p%precv(newsz-1)%base_a,p%precv(newsz-1)%base_desc,& + & ilaggr,nlaggr,op_prol,info) + exit array_build_loop + else + if (info == psb_success_) call mld_lev_mat_asb(p%precv(i),& + & p%precv(i-1)%base_a,p%precv(i-1)%base_desc,& + & ilaggr,nlaggr,op_prol,info) + end if if (info /= psb_success_) then call psb_errpush(psb_err_internal_error_,name,& & a_err='Map build') goto 9999 endif - if (newsz > 0) exit array_build_loop end do array_build_loop if (newsz > 0) then @@ -399,7 +394,6 @@ subroutine mld_s_hierarchy_bld(a,desc_a,p,info) end do end if - if (info /= psb_success_) then call psb_errpush(psb_err_internal_error_,name,& & a_err='Internal hierarchy build' ) @@ -419,4 +413,58 @@ subroutine mld_s_hierarchy_bld(a,desc_a,p,info) return +contains + subroutine save_smoothers(level,save1, save2,info) + type(mld_s_onelev_type), intent(in) :: level + class(mld_s_base_smoother_type), allocatable , intent(inout) :: save1, save2 + integer(psb_ipk_), intent(out) :: info + + info = 0 + if (allocated(save1)) then + call save1%free(info) + if (info == 0) deallocate(save1,stat=info) + if (info /= 0) return + end if + if (allocated(save2)) then + call save2%free(info) + if (info == 0) deallocate(save2,stat=info) + if (info /= 0) return + end if + allocate(save1, source=level%sm,stat=info) + if ((info == 0).and.allocated(level%sm2a)) allocate(save2, source=level%sm2a,stat=info) + + return + end subroutine save_smoothers + + subroutine restore_smoothers(level,save1, save2,info) + type(mld_s_onelev_type), intent(inout), target :: level + class(mld_s_base_smoother_type), allocatable, intent(in) :: save1, save2 + integer(psb_ipk_), intent(out) :: info + + info = 0 + + if (allocated(level%sm)) then + if (info == 0) call level%sm%free(info) + if (info == 0) deallocate(level%sm,stat=info) + end if + if (allocated(save1)) then + if (info == 0) allocate(level%sm,source=save1,stat=info) + end if + + if (info /= 0) return + + if (allocated(level%sm2a)) then + if (info == 0) call level%sm2a%free(info) + if (info == 0) deallocate(level%sm2a,stat=info) + end if + if (allocated(save2)) then + if (info == 0) allocate(level%sm2a,source=save2,stat=info) + if (info == 0) level%sm2 => level%sm2a + else + if (allocated(level%sm)) level%sm2 => level%sm + end if + + return + end subroutine restore_smoothers + end subroutine mld_s_hierarchy_bld diff --git a/mlprec/impl/mld_scprecset.F90 b/mlprec/impl/mld_scprecset.F90 index af641414..ae70a83b 100644 --- a/mlprec/impl/mld_scprecset.F90 +++ b/mlprec/impl/mld_scprecset.F90 @@ -519,6 +519,11 @@ subroutine mld_scprecsetr(p,what,val,info,ilev,pos) ilev_=nlev_ call p%precv(ilev_)%set('SUB_ILUTHRS',val,info,pos=pos) + case('AGGR_SCALE') + do ilev_ = 2, nlev_ + call p%precv(ilev_)%set('AGGR_SCALE',val,info,pos=pos) + end do + case('AGGR_THRESH') thr = val do ilev_ = 2, nlev_ diff --git a/mlprec/impl/mld_sprecset.F90 b/mlprec/impl/mld_sprecset.F90 index a8b5c289..4bc463ca 100644 --- a/mlprec/impl/mld_sprecset.F90 +++ b/mlprec/impl/mld_sprecset.F90 @@ -621,6 +621,11 @@ subroutine mld_sprecsetr(p,what,val,info,ilev,pos) ilev_=nlev_ call p%precv(ilev_)%set(mld_sub_iluthrs_,val,info,pos=pos) + case(mld_aggr_scale_) + do ilev_ = 2, nlev_ + call p%precv(ilev_)%set(mld_aggr_scale_,val,info,pos=pos) + end do + case(mld_aggr_thresh_) thr = val do ilev_ = 2, nlev_ diff --git a/mlprec/impl/mld_z_hierarchy_bld.f90 b/mlprec/impl/mld_z_hierarchy_bld.f90 index 33982c77..6b6a306e 100644 --- a/mlprec/impl/mld_z_hierarchy_bld.f90 +++ b/mlprec/impl/mld_z_hierarchy_bld.f90 @@ -82,7 +82,7 @@ subroutine mld_z_hierarchy_bld(a,desc_a,p,info) integer(psb_ipk_) :: ictxt, me,np integer(psb_ipk_) :: err,i,k, err_act, iszv, newsz, casize, nplevs, mxplevs, iaggsize real(psb_dpk_) :: mnaggratio, sizeratio - class(mld_z_base_smoother_type), allocatable :: coarse_sm, base_sm, med_sm + class(mld_z_base_smoother_type), allocatable :: coarse_sm, base_sm, med_sm, base_sm2, med_sm2, coarse_sm2 type(mld_dml_parms) :: baseparms, medparms, coarseparms integer(psb_ipk_), allocatable :: ilaggr(:), nlaggr(:) type(psb_zspmat_type) :: op_prol @@ -215,50 +215,40 @@ subroutine mld_z_hierarchy_bld(a,desc_a,p,info) baseparms = p%precv(1)%parms medparms = p%precv(2)%parms - allocate(coarse_sm, source=p%precv(iszv)%sm,stat=info) - if (info == psb_success_) & - & allocate(med_sm, source=p%precv(2)%sm,stat=info) - if (info == psb_success_) & - & allocate(base_sm, source=p%precv(1)%sm,stat=info) + call save_smoothers(p%precv(iszv),coarse_sm,coarse_sm2,info) + if (info == 0) call save_smoothers(p%precv(2),med_sm,med_sm2,info) + if (info == 0) call save_smoothers(p%precv(1),base_sm,base_sm2,info) if (info /= psb_success_) then write(0,*) 'Error in saving smoothers',info call psb_errpush(psb_err_internal_error_,name,a_err='Base level precbuild.') goto 9999 end if - ! ! First set desired number of levels ! if (iszv /= nplevs) then allocate(tprecv(nplevs),stat=info) - if (info /= psb_success_) then - call psb_errpush(psb_err_from_subroutine_,name,& - & a_err='prec reallocation') - goto 9999 - endif - tprecv(1)%parms = baseparms - allocate(tprecv(1)%sm,source=base_sm,stat=info) - if (info /= psb_success_) then - call psb_errpush(psb_err_from_subroutine_,name,& - & a_err='prec reallocation') - goto 9999 - endif - do i=2,nplevs-1 - tprecv(i)%parms = medparms - allocate(tprecv(i)%sm,source=med_sm,stat=info) - if (info /= psb_success_) then - call psb_errpush(psb_err_from_subroutine_,name,& - & a_err='prec reallocation') - goto 9999 - endif + ! First all existing levels + if (info == 0) tprecv(1)%parms = baseparms + if (info == 0) call restore_smoothers(tprecv(1),p%precv(1)%sm,p%precv(1)%sm2a,info) + do i=2, min(iszv,nplevs) - 1 + if (info == 0) tprecv(i)%parms = medparms + if (info == 0) call restore_smoothers(tprecv(i),p%precv(i)%sm,p%precv(i)%sm2a,info) + end do + ! Further intermediates, if any + do i=iszv-1, nplevs - 1 + if (info == 0) tprecv(i)%parms = medparms + if (info == 0) call restore_smoothers(tprecv(i),med_sm,med_sm2,info) end do - tprecv(nplevs)%parms = coarseparms - allocate(tprecv(nplevs)%sm,source=coarse_sm,stat=info) + ! Then coarse + if (info == 0) tprecv(nplevs)%parms = coarseparms + if (info == 0) call restore_smoothers(tprecv(nplevs),coarse_sm,coarse_sm2,info) if (info /= psb_success_) then call psb_errpush(psb_err_from_subroutine_,name,& & a_err='prec reallocation') goto 9999 endif + do i=1,iszv call p%precv(i)%free(info) end do @@ -330,7 +320,7 @@ subroutine mld_z_hierarchy_bld(a,desc_a,p,info) ! ! We are not gaining ! - newsz = newsz-1 + newsz = i-1 end if end if @@ -350,20 +340,25 @@ subroutine mld_z_hierarchy_bld(a,desc_a,p,info) end if end if call psb_bcast(ictxt,newsz) - if (newsz > 0) & - & call p%precv(i)%parms%get_coarse(p%precv(iszv)%parms) - - if (info == psb_success_) call mld_lev_mat_asb(p%precv(i),& - & p%precv(i-1)%base_a,p%precv(i-1)%base_desc,& - & ilaggr,nlaggr,op_prol,info) - + if (newsz > 0) then + if (info == 0) p%precv(newsz)%parms = coarseparms + if (info == 0) call restore_smoothers(p%precv(newsz),coarse_sm,coarse_sm2,info) + + if (info == psb_success_) call mld_lev_mat_asb(p%precv(newsz),& + & p%precv(newsz-1)%base_a,p%precv(newsz-1)%base_desc,& + & ilaggr,nlaggr,op_prol,info) + exit array_build_loop + else + if (info == psb_success_) call mld_lev_mat_asb(p%precv(i),& + & p%precv(i-1)%base_a,p%precv(i-1)%base_desc,& + & ilaggr,nlaggr,op_prol,info) + end if if (info /= psb_success_) then call psb_errpush(psb_err_internal_error_,name,& & a_err='Map build') goto 9999 endif - if (newsz > 0) exit array_build_loop end do array_build_loop if (newsz > 0) then @@ -399,7 +394,6 @@ subroutine mld_z_hierarchy_bld(a,desc_a,p,info) end do end if - if (info /= psb_success_) then call psb_errpush(psb_err_internal_error_,name,& & a_err='Internal hierarchy build' ) @@ -419,4 +413,58 @@ subroutine mld_z_hierarchy_bld(a,desc_a,p,info) return +contains + subroutine save_smoothers(level,save1, save2,info) + type(mld_z_onelev_type), intent(in) :: level + class(mld_z_base_smoother_type), allocatable , intent(inout) :: save1, save2 + integer(psb_ipk_), intent(out) :: info + + info = 0 + if (allocated(save1)) then + call save1%free(info) + if (info == 0) deallocate(save1,stat=info) + if (info /= 0) return + end if + if (allocated(save2)) then + call save2%free(info) + if (info == 0) deallocate(save2,stat=info) + if (info /= 0) return + end if + allocate(save1, source=level%sm,stat=info) + if ((info == 0).and.allocated(level%sm2a)) allocate(save2, source=level%sm2a,stat=info) + + return + end subroutine save_smoothers + + subroutine restore_smoothers(level,save1, save2,info) + type(mld_z_onelev_type), intent(inout), target :: level + class(mld_z_base_smoother_type), allocatable, intent(in) :: save1, save2 + integer(psb_ipk_), intent(out) :: info + + info = 0 + + if (allocated(level%sm)) then + if (info == 0) call level%sm%free(info) + if (info == 0) deallocate(level%sm,stat=info) + end if + if (allocated(save1)) then + if (info == 0) allocate(level%sm,source=save1,stat=info) + end if + + if (info /= 0) return + + if (allocated(level%sm2a)) then + if (info == 0) call level%sm2a%free(info) + if (info == 0) deallocate(level%sm2a,stat=info) + end if + if (allocated(save2)) then + if (info == 0) allocate(level%sm2a,source=save2,stat=info) + if (info == 0) level%sm2 => level%sm2a + else + if (allocated(level%sm)) level%sm2 => level%sm + end if + + return + end subroutine restore_smoothers + end subroutine mld_z_hierarchy_bld diff --git a/mlprec/impl/mld_zcprecset.F90 b/mlprec/impl/mld_zcprecset.F90 index 338616fb..abdd925b 100644 --- a/mlprec/impl/mld_zcprecset.F90 +++ b/mlprec/impl/mld_zcprecset.F90 @@ -529,6 +529,11 @@ subroutine mld_zcprecsetr(p,what,val,info,ilev,pos) ilev_=nlev_ call p%precv(ilev_)%set('SUB_ILUTHRS',val,info,pos=pos) + case('AGGR_SCALE') + do ilev_ = 2, nlev_ + call p%precv(ilev_)%set('AGGR_SCALE',val,info,pos=pos) + end do + case('AGGR_THRESH') thr = val do ilev_ = 2, nlev_ diff --git a/mlprec/impl/mld_zprecset.F90 b/mlprec/impl/mld_zprecset.F90 index e57c4c38..ab1a3f9c 100644 --- a/mlprec/impl/mld_zprecset.F90 +++ b/mlprec/impl/mld_zprecset.F90 @@ -631,6 +631,11 @@ subroutine mld_zprecsetr(p,what,val,info,ilev,pos) ilev_=nlev_ call p%precv(ilev_)%set(mld_sub_iluthrs_,val,info,pos=pos) + case(mld_aggr_scale_) + do ilev_ = 2, nlev_ + call p%precv(ilev_)%set(mld_aggr_scale_,val,info,pos=pos) + end do + case(mld_aggr_thresh_) thr = val do ilev_ = 2, nlev_