diff --git a/mlprec/impl/aggregator/mld_c_dec_aggregator_tprol.f90 b/mlprec/impl/aggregator/mld_c_dec_aggregator_tprol.f90 index 7a62c5eb..da8fe18c 100644 --- a/mlprec/impl/aggregator/mld_c_dec_aggregator_tprol.f90 +++ b/mlprec/impl/aggregator/mld_c_dec_aggregator_tprol.f90 @@ -40,7 +40,7 @@ ! Subroutine: mld_c_dec_aggregator_tprol ! Version: complex ! -! This routine is mainly an interface to map_bld where the real work is performed. +! This routine is mainly an interface to soc_map_bld where the real work is performed. ! It takes care of some consistency checking, and calls map_to_tprol, which is ! refactored and shared among all the aggregation methods that produce a simple ! integer mapping. @@ -50,6 +50,8 @@ ! ag - type(mld_c_dec_aggregator_type), input/output. ! The aggregator object, carrying with itself the mapping algorithm. ! parms - The auxiliary parameters object +! ag_data - Auxiliary global aggregation parameters object +! ! ! a - type(psb_cspmat_type). ! The sparse matrix structure containing the local part of the @@ -72,13 +74,15 @@ ! info - integer, output. ! Error code. ! -subroutine mld_c_dec_aggregator_build_tprol(ag,parms,a,desc_a,ilaggr,nlaggr,op_prol,info) +subroutine mld_c_dec_aggregator_build_tprol(ag,parms,ag_data,& + & a,desc_a,ilaggr,nlaggr,op_prol,info) use psb_base_mod use mld_c_prec_type, mld_protect_name => mld_c_dec_aggregator_build_tprol use mld_c_inner_mod implicit none class(mld_c_dec_aggregator_type), target, intent(inout) :: ag type(mld_sml_parms), intent(inout) :: parms + type(mld_saggr_data), intent(in) :: ag_data type(psb_cspmat_type), intent(in) :: a type(psb_desc_type), intent(in) :: desc_a integer(psb_lpk_), allocatable, intent(out) :: ilaggr(:),nlaggr(:) @@ -91,6 +95,7 @@ subroutine mld_c_dec_aggregator_build_tprol(ag,parms,a,desc_a,ilaggr,nlaggr,op_ integer(psb_ipk_) :: err_act integer(psb_lpk_) :: ntaggr integer(psb_ipk_) :: debug_level, debug_unit + logical :: clean_zeros name='mld_c_dec_aggregator_tprol' call psb_erractionsave(err_act) @@ -111,13 +116,17 @@ subroutine mld_c_dec_aggregator_build_tprol(ag,parms,a,desc_a,ilaggr,nlaggr,op_ & mld_aggr_ord_nat_,is_legal_ml_aggr_ord) call mld_check_def(parms%aggr_thresh,'Aggr_Thresh',szero,is_legal_s_aggr_thrs) - - call ag%map_bld(parms%aggr_ord,parms%aggr_thresh,a,desc_a,nlaggr,ilaggr,info) + ! + ! The decoupled aggregator based on SOC measures ignores + ! ag_data except for clean_zeros; soc_map_bld is a procedure pointer. + ! + clean_zeros = ag%do_clean_zeros + call ag%soc_map_bld(parms%aggr_ord,parms%aggr_thresh,clean_zeros,a,desc_a,nlaggr,ilaggr,info) if (info==psb_success_) call mld_map_to_tprol(desc_a,ilaggr,nlaggr,op_prol,info) if (info /= psb_success_) then info=psb_err_from_subroutine_ - call psb_errpush(info,name,a_err='map_bld/map_to_tprol') + call psb_errpush(info,name,a_err='soc_map_bld/map_to_tprol') goto 9999 endif diff --git a/mlprec/impl/aggregator/mld_c_soc1_map_bld.f90 b/mlprec/impl/aggregator/mld_c_soc1_map_bld.f90 index 53bebfe0..576df2d1 100644 --- a/mlprec/impl/aggregator/mld_c_soc1_map_bld.f90 +++ b/mlprec/impl/aggregator/mld_c_soc1_map_bld.f90 @@ -67,7 +67,7 @@ ! ! ! -subroutine mld_c_soc1_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info) +subroutine mld_c_soc1_map_bld(iorder,theta,clean_zeros,a,desc_a,nlaggr,ilaggr,info) use psb_base_mod use mld_base_prec_type @@ -77,6 +77,7 @@ subroutine mld_c_soc1_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info) ! Arguments integer(psb_ipk_), intent(in) :: iorder + logical, intent(in) :: clean_zeros type(psb_cspmat_type), intent(in) :: a type(psb_desc_type), intent(in) :: desc_a real(psb_spk_), intent(in) :: theta @@ -109,8 +110,8 @@ subroutine mld_c_soc1_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info) ! ictxt=desc_a%get_context() call psb_info(ictxt,me,np) - nrow = desc_a%get_local_rows() - ncol = desc_a%get_local_cols() + nrow = desc_a%get_local_rows() + ncol = desc_a%get_local_cols() nrglob = desc_a%get_global_rows() nr = a%get_nrows() @@ -132,6 +133,7 @@ subroutine mld_c_soc1_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info) end if call a%cp_to(acsr) + if (clean_zeros) call acsr%clean_zeros(info) if (iorder == mld_aggr_ord_nat_) then do i=1, nr ilaggr(i) = -(nr+1) @@ -169,12 +171,6 @@ subroutine mld_c_soc1_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info) icol(1:nz) = acsr%ja(acsr%irp(i):acsr%irp(i+1)-1) val(1:nz) = acsr%val(acsr%irp(i):acsr%irp(i+1)-1) -!!$ call a%csget(i,i,nz,irow,icol,val,info,chksz=.false.) -!!$ if (info /= psb_success_) then -!!$ info=psb_err_from_subroutine_ -!!$ call psb_errpush(info,name,a_err='csget') -!!$ goto 9999 -!!$ end if ! ! Build the set of all strongly coupled nodes @@ -194,7 +190,7 @@ subroutine mld_c_soc1_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info) ! If the whole strongly coupled neighborhood of I is ! as yet unconnected, turn it into the next aggregate. ! Same if ip==0 (in which case, neighborhood only - ! contains I even if it does not look from matrix) + ! contains I even if it does not look like it from matrix) ! disjoint = all(ilaggr(icol(1:ip)) == -(nr+1)).or.(ip==0) if (disjoint) then @@ -222,14 +218,10 @@ subroutine mld_c_soc1_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info) if (ilaggr(i) == -(nr+1)) then nz = (acsr%irp(i+1)-acsr%irp(i)) + if (nz == 1) cycle step2 icol(1:nz) = acsr%ja(acsr%irp(i):acsr%irp(i+1)-1) val(1:nz) = acsr%val(acsr%irp(i):acsr%irp(i+1)-1) -!!$ call a%csget(i,i,nz,irow,icol,val,info,chksz=.false.) -!!$ if (info /= psb_success_) then -!!$ info=psb_err_from_subroutine_ -!!$ call psb_errpush(info,name,a_err='psb_sp_getrow') -!!$ goto 9999 -!!$ end if + ! ! Find the most strongly connected neighbour that is ! already aggregated, if any, and join its aggregate @@ -261,6 +253,7 @@ subroutine mld_c_soc1_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info) if (ilaggr(i) < 0) then nz = (acsr%irp(i+1)-acsr%irp(i)) + if (nz == 1) cycle step3 icol(1:nz) = acsr%ja(acsr%irp(i):acsr%irp(i+1)-1) val(1:nz) = acsr%val(acsr%irp(i):acsr%irp(i+1)-1) ! @@ -288,8 +281,8 @@ subroutine mld_c_soc1_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info) end do else ! - ! This should not happen: we did not even connect with ourselves. - ! Create an isolate anyway. + ! This should not happen: we did not even connect with ourselves, + ! but it's not a singleton. ! naggr = naggr + 1 ilaggr(i) = naggr @@ -297,7 +290,6 @@ subroutine mld_c_soc1_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info) end if end do step3 - ! Any leftovers? do i=1, nr if (ilaggr(i) < 0) then diff --git a/mlprec/impl/aggregator/mld_c_soc2_map_bld.f90 b/mlprec/impl/aggregator/mld_c_soc2_map_bld.f90 index 6b9bdc6e..81b43da5 100644 --- a/mlprec/impl/aggregator/mld_c_soc2_map_bld.f90 +++ b/mlprec/impl/aggregator/mld_c_soc2_map_bld.f90 @@ -66,7 +66,7 @@ ! ! ! -subroutine mld_c_soc2_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info) +subroutine mld_c_soc2_map_bld(iorder,theta,clean_zeros,a,desc_a,nlaggr,ilaggr,info) use psb_base_mod use mld_base_prec_type @@ -76,6 +76,7 @@ subroutine mld_c_soc2_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info) ! Arguments integer(psb_ipk_), intent(in) :: iorder + logical, intent(in) :: clean_zeros type(psb_cspmat_type), intent(in) :: a type(psb_desc_type), intent(in) :: desc_a real(psb_spk_), intent(in) :: theta @@ -134,13 +135,14 @@ subroutine mld_c_soc2_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info) ! Phase zero: compute muij ! call a%cp_to(muij) + if (clean_zeros) call muij%clean_zeros(info) do i=1, nr do k=muij%irp(i),muij%irp(i+1)-1 j = muij%ja(k) if (j<= nr) muij%val(k) = abs(muij%val(k))/sqrt(abs(diag(i)*diag(j))) end do end do - !write(*,*) 'murows/cols ',maxval(mu_rows),maxval(mu_cols) + ! ! Compute the 1-neigbour; mark strong links with +1, weak links with -1 ! diff --git a/mlprec/impl/aggregator/mld_c_symdec_aggregator_tprol.f90 b/mlprec/impl/aggregator/mld_c_symdec_aggregator_tprol.f90 index 256223f6..1ca19694 100644 --- a/mlprec/impl/aggregator/mld_c_symdec_aggregator_tprol.f90 +++ b/mlprec/impl/aggregator/mld_c_symdec_aggregator_tprol.f90 @@ -41,7 +41,7 @@ ! Version: complex ! ! -! This routine is mainly an interface to map_bld where the real work is performed. +! This routine is mainly an interface to soc_map_bld where the real work is performed. ! It takes care of some consistency checking, and calls map_to_tprol, which is ! refactored and shared among all the aggregation methods that produce a simple ! integer mapping. It also symmetrizes the pattern of the local matrix A. @@ -53,6 +53,7 @@ ! ag - type(mld_c_dec_aggregator_type), input/output. ! The aggregator object, carrying with itself the mapping algorithm. ! parms - The auxiliary parameters object +! ag_data - Auxiliary global aggregation parameters object ! ! a - type(psb_cspmat_type). ! The sparse matrix structure containing the local part of the @@ -75,7 +76,8 @@ ! info - integer, output. ! Error code. ! -subroutine mld_c_symdec_aggregator_build_tprol(ag,parms,a,desc_a,ilaggr,nlaggr,op_prol,info) +subroutine mld_c_symdec_aggregator_build_tprol(ag,parms,ag_data,& + & a,desc_a,ilaggr,nlaggr,op_prol,info) use psb_base_mod use mld_c_prec_type use mld_c_symdec_aggregator_mod, mld_protect_name => mld_c_symdec_aggregator_build_tprol @@ -83,6 +85,7 @@ subroutine mld_c_symdec_aggregator_build_tprol(ag,parms,a,desc_a,ilaggr,nlaggr, implicit none class(mld_c_symdec_aggregator_type), target, intent(inout) :: ag type(mld_sml_parms), intent(inout) :: parms + type(mld_saggr_data), intent(in) :: ag_data type(psb_cspmat_type), intent(in) :: a type(psb_desc_type), intent(in) :: desc_a integer(psb_lpk_), allocatable, intent(out) :: ilaggr(:),nlaggr(:) @@ -97,6 +100,7 @@ subroutine mld_c_symdec_aggregator_build_tprol(ag,parms,a,desc_a,ilaggr,nlaggr, integer(psb_ipk_) :: nr integer(psb_lpk_) :: ntaggr integer(psb_ipk_) :: debug_level, debug_unit + logical :: clean_zeros name='mld_c_symdec_aggregator_tprol' call psb_erractionsave(err_act) @@ -130,14 +134,20 @@ subroutine mld_c_symdec_aggregator_build_tprol(ag,parms,a,desc_a,ilaggr,nlaggr, if (info == psb_success_) call atrans%free() if (info == psb_success_) call atmp%cscnv(info,type='CSR') + ! + ! The decoupled aggregator based on SOC measures ignores + ! ag_data except for clean_zeros; soc_map_bld is a procedure pointer. + ! + clean_zeros = ag%do_clean_zeros if (info == psb_success_) & - & call ag%map_bld(parms%aggr_ord,parms%aggr_thresh,atmp,desc_a,nlaggr,ilaggr,info) + & call ag%soc_map_bld(parms%aggr_ord,parms%aggr_thresh,clean_zeros,atmp,& + & desc_a,nlaggr,ilaggr,info) if (info == psb_success_) call atmp%free() if (info == psb_success_) call mld_map_to_tprol(desc_a,ilaggr,nlaggr,op_prol,info) if (info /= psb_success_) then info=psb_err_from_subroutine_ - call psb_errpush(info,name,a_err='map_bld/map_to_tprol') + call psb_errpush(info,name,a_err='soc_map_bld/map_to_tprol') goto 9999 endif diff --git a/mlprec/impl/aggregator/mld_d_dec_aggregator_tprol.f90 b/mlprec/impl/aggregator/mld_d_dec_aggregator_tprol.f90 index d38e021d..4f6a88a8 100644 --- a/mlprec/impl/aggregator/mld_d_dec_aggregator_tprol.f90 +++ b/mlprec/impl/aggregator/mld_d_dec_aggregator_tprol.f90 @@ -40,7 +40,7 @@ ! Subroutine: mld_d_dec_aggregator_tprol ! Version: real ! -! This routine is mainly an interface to map_bld where the real work is performed. +! This routine is mainly an interface to soc_map_bld where the real work is performed. ! It takes care of some consistency checking, and calls map_to_tprol, which is ! refactored and shared among all the aggregation methods that produce a simple ! integer mapping. @@ -50,6 +50,8 @@ ! ag - type(mld_d_dec_aggregator_type), input/output. ! The aggregator object, carrying with itself the mapping algorithm. ! parms - The auxiliary parameters object +! ag_data - Auxiliary global aggregation parameters object +! ! ! a - type(psb_dspmat_type). ! The sparse matrix structure containing the local part of the @@ -72,13 +74,15 @@ ! info - integer, output. ! Error code. ! -subroutine mld_d_dec_aggregator_build_tprol(ag,parms,a,desc_a,ilaggr,nlaggr,op_prol,info) +subroutine mld_d_dec_aggregator_build_tprol(ag,parms,ag_data,& + & a,desc_a,ilaggr,nlaggr,op_prol,info) use psb_base_mod use mld_d_prec_type, mld_protect_name => mld_d_dec_aggregator_build_tprol use mld_d_inner_mod implicit none class(mld_d_dec_aggregator_type), target, intent(inout) :: ag type(mld_dml_parms), intent(inout) :: parms + type(mld_daggr_data), intent(in) :: ag_data type(psb_dspmat_type), intent(in) :: a type(psb_desc_type), intent(in) :: desc_a integer(psb_lpk_), allocatable, intent(out) :: ilaggr(:),nlaggr(:) @@ -91,6 +95,7 @@ subroutine mld_d_dec_aggregator_build_tprol(ag,parms,a,desc_a,ilaggr,nlaggr,op_ integer(psb_ipk_) :: err_act integer(psb_lpk_) :: ntaggr integer(psb_ipk_) :: debug_level, debug_unit + logical :: clean_zeros name='mld_d_dec_aggregator_tprol' call psb_erractionsave(err_act) @@ -111,13 +116,17 @@ subroutine mld_d_dec_aggregator_build_tprol(ag,parms,a,desc_a,ilaggr,nlaggr,op_ & mld_aggr_ord_nat_,is_legal_ml_aggr_ord) call mld_check_def(parms%aggr_thresh,'Aggr_Thresh',dzero,is_legal_d_aggr_thrs) - - call ag%map_bld(parms%aggr_ord,parms%aggr_thresh,a,desc_a,nlaggr,ilaggr,info) + ! + ! The decoupled aggregator based on SOC measures ignores + ! ag_data except for clean_zeros; soc_map_bld is a procedure pointer. + ! + clean_zeros = ag%do_clean_zeros + call ag%soc_map_bld(parms%aggr_ord,parms%aggr_thresh,clean_zeros,a,desc_a,nlaggr,ilaggr,info) if (info==psb_success_) call mld_map_to_tprol(desc_a,ilaggr,nlaggr,op_prol,info) if (info /= psb_success_) then info=psb_err_from_subroutine_ - call psb_errpush(info,name,a_err='map_bld/map_to_tprol') + call psb_errpush(info,name,a_err='soc_map_bld/map_to_tprol') goto 9999 endif diff --git a/mlprec/impl/aggregator/mld_d_soc1_map_bld.f90 b/mlprec/impl/aggregator/mld_d_soc1_map_bld.f90 index 24b83650..4da28284 100644 --- a/mlprec/impl/aggregator/mld_d_soc1_map_bld.f90 +++ b/mlprec/impl/aggregator/mld_d_soc1_map_bld.f90 @@ -67,7 +67,7 @@ ! ! ! -subroutine mld_d_soc1_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info) +subroutine mld_d_soc1_map_bld(iorder,theta,clean_zeros,a,desc_a,nlaggr,ilaggr,info) use psb_base_mod use mld_base_prec_type @@ -77,6 +77,7 @@ subroutine mld_d_soc1_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info) ! Arguments integer(psb_ipk_), intent(in) :: iorder + logical, intent(in) :: clean_zeros type(psb_dspmat_type), intent(in) :: a type(psb_desc_type), intent(in) :: desc_a real(psb_dpk_), intent(in) :: theta @@ -109,8 +110,8 @@ subroutine mld_d_soc1_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info) ! ictxt=desc_a%get_context() call psb_info(ictxt,me,np) - nrow = desc_a%get_local_rows() - ncol = desc_a%get_local_cols() + nrow = desc_a%get_local_rows() + ncol = desc_a%get_local_cols() nrglob = desc_a%get_global_rows() nr = a%get_nrows() @@ -132,6 +133,7 @@ subroutine mld_d_soc1_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info) end if call a%cp_to(acsr) + if (clean_zeros) call acsr%clean_zeros(info) if (iorder == mld_aggr_ord_nat_) then do i=1, nr ilaggr(i) = -(nr+1) @@ -169,12 +171,6 @@ subroutine mld_d_soc1_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info) icol(1:nz) = acsr%ja(acsr%irp(i):acsr%irp(i+1)-1) val(1:nz) = acsr%val(acsr%irp(i):acsr%irp(i+1)-1) -!!$ call a%csget(i,i,nz,irow,icol,val,info,chksz=.false.) -!!$ if (info /= psb_success_) then -!!$ info=psb_err_from_subroutine_ -!!$ call psb_errpush(info,name,a_err='csget') -!!$ goto 9999 -!!$ end if ! ! Build the set of all strongly coupled nodes @@ -194,7 +190,7 @@ subroutine mld_d_soc1_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info) ! If the whole strongly coupled neighborhood of I is ! as yet unconnected, turn it into the next aggregate. ! Same if ip==0 (in which case, neighborhood only - ! contains I even if it does not look from matrix) + ! contains I even if it does not look like it from matrix) ! disjoint = all(ilaggr(icol(1:ip)) == -(nr+1)).or.(ip==0) if (disjoint) then @@ -222,14 +218,10 @@ subroutine mld_d_soc1_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info) if (ilaggr(i) == -(nr+1)) then nz = (acsr%irp(i+1)-acsr%irp(i)) + if (nz == 1) cycle step2 icol(1:nz) = acsr%ja(acsr%irp(i):acsr%irp(i+1)-1) val(1:nz) = acsr%val(acsr%irp(i):acsr%irp(i+1)-1) -!!$ call a%csget(i,i,nz,irow,icol,val,info,chksz=.false.) -!!$ if (info /= psb_success_) then -!!$ info=psb_err_from_subroutine_ -!!$ call psb_errpush(info,name,a_err='psb_sp_getrow') -!!$ goto 9999 -!!$ end if + ! ! Find the most strongly connected neighbour that is ! already aggregated, if any, and join its aggregate @@ -261,6 +253,7 @@ subroutine mld_d_soc1_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info) if (ilaggr(i) < 0) then nz = (acsr%irp(i+1)-acsr%irp(i)) + if (nz == 1) cycle step3 icol(1:nz) = acsr%ja(acsr%irp(i):acsr%irp(i+1)-1) val(1:nz) = acsr%val(acsr%irp(i):acsr%irp(i+1)-1) ! @@ -288,8 +281,8 @@ subroutine mld_d_soc1_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info) end do else ! - ! This should not happen: we did not even connect with ourselves. - ! Create an isolate anyway. + ! This should not happen: we did not even connect with ourselves, + ! but it's not a singleton. ! naggr = naggr + 1 ilaggr(i) = naggr @@ -297,7 +290,6 @@ subroutine mld_d_soc1_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info) end if end do step3 - ! Any leftovers? do i=1, nr if (ilaggr(i) < 0) then diff --git a/mlprec/impl/aggregator/mld_d_soc2_map_bld.f90 b/mlprec/impl/aggregator/mld_d_soc2_map_bld.f90 index 0307ca73..ae8ceedc 100644 --- a/mlprec/impl/aggregator/mld_d_soc2_map_bld.f90 +++ b/mlprec/impl/aggregator/mld_d_soc2_map_bld.f90 @@ -66,7 +66,7 @@ ! ! ! -subroutine mld_d_soc2_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info) +subroutine mld_d_soc2_map_bld(iorder,theta,clean_zeros,a,desc_a,nlaggr,ilaggr,info) use psb_base_mod use mld_base_prec_type @@ -76,6 +76,7 @@ subroutine mld_d_soc2_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info) ! Arguments integer(psb_ipk_), intent(in) :: iorder + logical, intent(in) :: clean_zeros type(psb_dspmat_type), intent(in) :: a type(psb_desc_type), intent(in) :: desc_a real(psb_dpk_), intent(in) :: theta @@ -134,13 +135,14 @@ subroutine mld_d_soc2_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info) ! Phase zero: compute muij ! call a%cp_to(muij) + if (clean_zeros) call muij%clean_zeros(info) do i=1, nr do k=muij%irp(i),muij%irp(i+1)-1 j = muij%ja(k) if (j<= nr) muij%val(k) = abs(muij%val(k))/sqrt(abs(diag(i)*diag(j))) end do end do - !write(*,*) 'murows/cols ',maxval(mu_rows),maxval(mu_cols) + ! ! Compute the 1-neigbour; mark strong links with +1, weak links with -1 ! diff --git a/mlprec/impl/aggregator/mld_d_symdec_aggregator_tprol.f90 b/mlprec/impl/aggregator/mld_d_symdec_aggregator_tprol.f90 index 12466c00..f9f80452 100644 --- a/mlprec/impl/aggregator/mld_d_symdec_aggregator_tprol.f90 +++ b/mlprec/impl/aggregator/mld_d_symdec_aggregator_tprol.f90 @@ -41,7 +41,7 @@ ! Version: real ! ! -! This routine is mainly an interface to map_bld where the real work is performed. +! This routine is mainly an interface to soc_map_bld where the real work is performed. ! It takes care of some consistency checking, and calls map_to_tprol, which is ! refactored and shared among all the aggregation methods that produce a simple ! integer mapping. It also symmetrizes the pattern of the local matrix A. @@ -53,6 +53,7 @@ ! ag - type(mld_d_dec_aggregator_type), input/output. ! The aggregator object, carrying with itself the mapping algorithm. ! parms - The auxiliary parameters object +! ag_data - Auxiliary global aggregation parameters object ! ! a - type(psb_dspmat_type). ! The sparse matrix structure containing the local part of the @@ -75,7 +76,8 @@ ! info - integer, output. ! Error code. ! -subroutine mld_d_symdec_aggregator_build_tprol(ag,parms,a,desc_a,ilaggr,nlaggr,op_prol,info) +subroutine mld_d_symdec_aggregator_build_tprol(ag,parms,ag_data,& + & a,desc_a,ilaggr,nlaggr,op_prol,info) use psb_base_mod use mld_d_prec_type use mld_d_symdec_aggregator_mod, mld_protect_name => mld_d_symdec_aggregator_build_tprol @@ -83,6 +85,7 @@ subroutine mld_d_symdec_aggregator_build_tprol(ag,parms,a,desc_a,ilaggr,nlaggr, implicit none class(mld_d_symdec_aggregator_type), target, intent(inout) :: ag type(mld_dml_parms), intent(inout) :: parms + type(mld_daggr_data), intent(in) :: ag_data type(psb_dspmat_type), intent(in) :: a type(psb_desc_type), intent(in) :: desc_a integer(psb_lpk_), allocatable, intent(out) :: ilaggr(:),nlaggr(:) @@ -97,6 +100,7 @@ subroutine mld_d_symdec_aggregator_build_tprol(ag,parms,a,desc_a,ilaggr,nlaggr, integer(psb_ipk_) :: nr integer(psb_lpk_) :: ntaggr integer(psb_ipk_) :: debug_level, debug_unit + logical :: clean_zeros name='mld_d_symdec_aggregator_tprol' call psb_erractionsave(err_act) @@ -130,14 +134,20 @@ subroutine mld_d_symdec_aggregator_build_tprol(ag,parms,a,desc_a,ilaggr,nlaggr, if (info == psb_success_) call atrans%free() if (info == psb_success_) call atmp%cscnv(info,type='CSR') + ! + ! The decoupled aggregator based on SOC measures ignores + ! ag_data except for clean_zeros; soc_map_bld is a procedure pointer. + ! + clean_zeros = ag%do_clean_zeros if (info == psb_success_) & - & call ag%map_bld(parms%aggr_ord,parms%aggr_thresh,atmp,desc_a,nlaggr,ilaggr,info) + & call ag%soc_map_bld(parms%aggr_ord,parms%aggr_thresh,clean_zeros,atmp,& + & desc_a,nlaggr,ilaggr,info) if (info == psb_success_) call atmp%free() if (info == psb_success_) call mld_map_to_tprol(desc_a,ilaggr,nlaggr,op_prol,info) if (info /= psb_success_) then info=psb_err_from_subroutine_ - call psb_errpush(info,name,a_err='map_bld/map_to_tprol') + call psb_errpush(info,name,a_err='soc_map_bld/map_to_tprol') goto 9999 endif diff --git a/mlprec/impl/aggregator/mld_s_dec_aggregator_tprol.f90 b/mlprec/impl/aggregator/mld_s_dec_aggregator_tprol.f90 index 8ab80eb3..e6b43cf2 100644 --- a/mlprec/impl/aggregator/mld_s_dec_aggregator_tprol.f90 +++ b/mlprec/impl/aggregator/mld_s_dec_aggregator_tprol.f90 @@ -40,7 +40,7 @@ ! Subroutine: mld_s_dec_aggregator_tprol ! Version: real ! -! This routine is mainly an interface to map_bld where the real work is performed. +! This routine is mainly an interface to soc_map_bld where the real work is performed. ! It takes care of some consistency checking, and calls map_to_tprol, which is ! refactored and shared among all the aggregation methods that produce a simple ! integer mapping. @@ -50,6 +50,8 @@ ! ag - type(mld_s_dec_aggregator_type), input/output. ! The aggregator object, carrying with itself the mapping algorithm. ! parms - The auxiliary parameters object +! ag_data - Auxiliary global aggregation parameters object +! ! ! a - type(psb_sspmat_type). ! The sparse matrix structure containing the local part of the @@ -72,13 +74,15 @@ ! info - integer, output. ! Error code. ! -subroutine mld_s_dec_aggregator_build_tprol(ag,parms,a,desc_a,ilaggr,nlaggr,op_prol,info) +subroutine mld_s_dec_aggregator_build_tprol(ag,parms,ag_data,& + & a,desc_a,ilaggr,nlaggr,op_prol,info) use psb_base_mod use mld_s_prec_type, mld_protect_name => mld_s_dec_aggregator_build_tprol use mld_s_inner_mod implicit none class(mld_s_dec_aggregator_type), target, intent(inout) :: ag type(mld_sml_parms), intent(inout) :: parms + type(mld_saggr_data), intent(in) :: ag_data type(psb_sspmat_type), intent(in) :: a type(psb_desc_type), intent(in) :: desc_a integer(psb_lpk_), allocatable, intent(out) :: ilaggr(:),nlaggr(:) @@ -91,6 +95,7 @@ subroutine mld_s_dec_aggregator_build_tprol(ag,parms,a,desc_a,ilaggr,nlaggr,op_ integer(psb_ipk_) :: err_act integer(psb_lpk_) :: ntaggr integer(psb_ipk_) :: debug_level, debug_unit + logical :: clean_zeros name='mld_s_dec_aggregator_tprol' call psb_erractionsave(err_act) @@ -111,13 +116,17 @@ subroutine mld_s_dec_aggregator_build_tprol(ag,parms,a,desc_a,ilaggr,nlaggr,op_ & mld_aggr_ord_nat_,is_legal_ml_aggr_ord) call mld_check_def(parms%aggr_thresh,'Aggr_Thresh',szero,is_legal_s_aggr_thrs) - - call ag%map_bld(parms%aggr_ord,parms%aggr_thresh,a,desc_a,nlaggr,ilaggr,info) + ! + ! The decoupled aggregator based on SOC measures ignores + ! ag_data except for clean_zeros; soc_map_bld is a procedure pointer. + ! + clean_zeros = ag%do_clean_zeros + call ag%soc_map_bld(parms%aggr_ord,parms%aggr_thresh,clean_zeros,a,desc_a,nlaggr,ilaggr,info) if (info==psb_success_) call mld_map_to_tprol(desc_a,ilaggr,nlaggr,op_prol,info) if (info /= psb_success_) then info=psb_err_from_subroutine_ - call psb_errpush(info,name,a_err='map_bld/map_to_tprol') + call psb_errpush(info,name,a_err='soc_map_bld/map_to_tprol') goto 9999 endif diff --git a/mlprec/impl/aggregator/mld_s_soc1_map_bld.f90 b/mlprec/impl/aggregator/mld_s_soc1_map_bld.f90 index 4c197a1a..7992b425 100644 --- a/mlprec/impl/aggregator/mld_s_soc1_map_bld.f90 +++ b/mlprec/impl/aggregator/mld_s_soc1_map_bld.f90 @@ -67,7 +67,7 @@ ! ! ! -subroutine mld_s_soc1_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info) +subroutine mld_s_soc1_map_bld(iorder,theta,clean_zeros,a,desc_a,nlaggr,ilaggr,info) use psb_base_mod use mld_base_prec_type @@ -77,6 +77,7 @@ subroutine mld_s_soc1_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info) ! Arguments integer(psb_ipk_), intent(in) :: iorder + logical, intent(in) :: clean_zeros type(psb_sspmat_type), intent(in) :: a type(psb_desc_type), intent(in) :: desc_a real(psb_spk_), intent(in) :: theta @@ -109,8 +110,8 @@ subroutine mld_s_soc1_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info) ! ictxt=desc_a%get_context() call psb_info(ictxt,me,np) - nrow = desc_a%get_local_rows() - ncol = desc_a%get_local_cols() + nrow = desc_a%get_local_rows() + ncol = desc_a%get_local_cols() nrglob = desc_a%get_global_rows() nr = a%get_nrows() @@ -132,6 +133,7 @@ subroutine mld_s_soc1_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info) end if call a%cp_to(acsr) + if (clean_zeros) call acsr%clean_zeros(info) if (iorder == mld_aggr_ord_nat_) then do i=1, nr ilaggr(i) = -(nr+1) @@ -169,12 +171,6 @@ subroutine mld_s_soc1_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info) icol(1:nz) = acsr%ja(acsr%irp(i):acsr%irp(i+1)-1) val(1:nz) = acsr%val(acsr%irp(i):acsr%irp(i+1)-1) -!!$ call a%csget(i,i,nz,irow,icol,val,info,chksz=.false.) -!!$ if (info /= psb_success_) then -!!$ info=psb_err_from_subroutine_ -!!$ call psb_errpush(info,name,a_err='csget') -!!$ goto 9999 -!!$ end if ! ! Build the set of all strongly coupled nodes @@ -194,7 +190,7 @@ subroutine mld_s_soc1_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info) ! If the whole strongly coupled neighborhood of I is ! as yet unconnected, turn it into the next aggregate. ! Same if ip==0 (in which case, neighborhood only - ! contains I even if it does not look from matrix) + ! contains I even if it does not look like it from matrix) ! disjoint = all(ilaggr(icol(1:ip)) == -(nr+1)).or.(ip==0) if (disjoint) then @@ -222,14 +218,10 @@ subroutine mld_s_soc1_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info) if (ilaggr(i) == -(nr+1)) then nz = (acsr%irp(i+1)-acsr%irp(i)) + if (nz == 1) cycle step2 icol(1:nz) = acsr%ja(acsr%irp(i):acsr%irp(i+1)-1) val(1:nz) = acsr%val(acsr%irp(i):acsr%irp(i+1)-1) -!!$ call a%csget(i,i,nz,irow,icol,val,info,chksz=.false.) -!!$ if (info /= psb_success_) then -!!$ info=psb_err_from_subroutine_ -!!$ call psb_errpush(info,name,a_err='psb_sp_getrow') -!!$ goto 9999 -!!$ end if + ! ! Find the most strongly connected neighbour that is ! already aggregated, if any, and join its aggregate @@ -261,6 +253,7 @@ subroutine mld_s_soc1_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info) if (ilaggr(i) < 0) then nz = (acsr%irp(i+1)-acsr%irp(i)) + if (nz == 1) cycle step3 icol(1:nz) = acsr%ja(acsr%irp(i):acsr%irp(i+1)-1) val(1:nz) = acsr%val(acsr%irp(i):acsr%irp(i+1)-1) ! @@ -288,8 +281,8 @@ subroutine mld_s_soc1_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info) end do else ! - ! This should not happen: we did not even connect with ourselves. - ! Create an isolate anyway. + ! This should not happen: we did not even connect with ourselves, + ! but it's not a singleton. ! naggr = naggr + 1 ilaggr(i) = naggr @@ -297,7 +290,6 @@ subroutine mld_s_soc1_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info) end if end do step3 - ! Any leftovers? do i=1, nr if (ilaggr(i) < 0) then diff --git a/mlprec/impl/aggregator/mld_s_soc2_map_bld.f90 b/mlprec/impl/aggregator/mld_s_soc2_map_bld.f90 index 770c51ac..7fcce613 100644 --- a/mlprec/impl/aggregator/mld_s_soc2_map_bld.f90 +++ b/mlprec/impl/aggregator/mld_s_soc2_map_bld.f90 @@ -66,7 +66,7 @@ ! ! ! -subroutine mld_s_soc2_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info) +subroutine mld_s_soc2_map_bld(iorder,theta,clean_zeros,a,desc_a,nlaggr,ilaggr,info) use psb_base_mod use mld_base_prec_type @@ -76,6 +76,7 @@ subroutine mld_s_soc2_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info) ! Arguments integer(psb_ipk_), intent(in) :: iorder + logical, intent(in) :: clean_zeros type(psb_sspmat_type), intent(in) :: a type(psb_desc_type), intent(in) :: desc_a real(psb_spk_), intent(in) :: theta @@ -134,13 +135,14 @@ subroutine mld_s_soc2_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info) ! Phase zero: compute muij ! call a%cp_to(muij) + if (clean_zeros) call muij%clean_zeros(info) do i=1, nr do k=muij%irp(i),muij%irp(i+1)-1 j = muij%ja(k) if (j<= nr) muij%val(k) = abs(muij%val(k))/sqrt(abs(diag(i)*diag(j))) end do end do - !write(*,*) 'murows/cols ',maxval(mu_rows),maxval(mu_cols) + ! ! Compute the 1-neigbour; mark strong links with +1, weak links with -1 ! diff --git a/mlprec/impl/aggregator/mld_s_symdec_aggregator_tprol.f90 b/mlprec/impl/aggregator/mld_s_symdec_aggregator_tprol.f90 index 7f493c55..007250e9 100644 --- a/mlprec/impl/aggregator/mld_s_symdec_aggregator_tprol.f90 +++ b/mlprec/impl/aggregator/mld_s_symdec_aggregator_tprol.f90 @@ -41,7 +41,7 @@ ! Version: real ! ! -! This routine is mainly an interface to map_bld where the real work is performed. +! This routine is mainly an interface to soc_map_bld where the real work is performed. ! It takes care of some consistency checking, and calls map_to_tprol, which is ! refactored and shared among all the aggregation methods that produce a simple ! integer mapping. It also symmetrizes the pattern of the local matrix A. @@ -53,6 +53,7 @@ ! ag - type(mld_s_dec_aggregator_type), input/output. ! The aggregator object, carrying with itself the mapping algorithm. ! parms - The auxiliary parameters object +! ag_data - Auxiliary global aggregation parameters object ! ! a - type(psb_sspmat_type). ! The sparse matrix structure containing the local part of the @@ -75,7 +76,8 @@ ! info - integer, output. ! Error code. ! -subroutine mld_s_symdec_aggregator_build_tprol(ag,parms,a,desc_a,ilaggr,nlaggr,op_prol,info) +subroutine mld_s_symdec_aggregator_build_tprol(ag,parms,ag_data,& + & a,desc_a,ilaggr,nlaggr,op_prol,info) use psb_base_mod use mld_s_prec_type use mld_s_symdec_aggregator_mod, mld_protect_name => mld_s_symdec_aggregator_build_tprol @@ -83,6 +85,7 @@ subroutine mld_s_symdec_aggregator_build_tprol(ag,parms,a,desc_a,ilaggr,nlaggr, implicit none class(mld_s_symdec_aggregator_type), target, intent(inout) :: ag type(mld_sml_parms), intent(inout) :: parms + type(mld_saggr_data), intent(in) :: ag_data type(psb_sspmat_type), intent(in) :: a type(psb_desc_type), intent(in) :: desc_a integer(psb_lpk_), allocatable, intent(out) :: ilaggr(:),nlaggr(:) @@ -97,6 +100,7 @@ subroutine mld_s_symdec_aggregator_build_tprol(ag,parms,a,desc_a,ilaggr,nlaggr, integer(psb_ipk_) :: nr integer(psb_lpk_) :: ntaggr integer(psb_ipk_) :: debug_level, debug_unit + logical :: clean_zeros name='mld_s_symdec_aggregator_tprol' call psb_erractionsave(err_act) @@ -130,14 +134,20 @@ subroutine mld_s_symdec_aggregator_build_tprol(ag,parms,a,desc_a,ilaggr,nlaggr, if (info == psb_success_) call atrans%free() if (info == psb_success_) call atmp%cscnv(info,type='CSR') + ! + ! The decoupled aggregator based on SOC measures ignores + ! ag_data except for clean_zeros; soc_map_bld is a procedure pointer. + ! + clean_zeros = ag%do_clean_zeros if (info == psb_success_) & - & call ag%map_bld(parms%aggr_ord,parms%aggr_thresh,atmp,desc_a,nlaggr,ilaggr,info) + & call ag%soc_map_bld(parms%aggr_ord,parms%aggr_thresh,clean_zeros,atmp,& + & desc_a,nlaggr,ilaggr,info) if (info == psb_success_) call atmp%free() if (info == psb_success_) call mld_map_to_tprol(desc_a,ilaggr,nlaggr,op_prol,info) if (info /= psb_success_) then info=psb_err_from_subroutine_ - call psb_errpush(info,name,a_err='map_bld/map_to_tprol') + call psb_errpush(info,name,a_err='soc_map_bld/map_to_tprol') goto 9999 endif diff --git a/mlprec/impl/aggregator/mld_z_dec_aggregator_tprol.f90 b/mlprec/impl/aggregator/mld_z_dec_aggregator_tprol.f90 index ccd0dcbc..75f3b7da 100644 --- a/mlprec/impl/aggregator/mld_z_dec_aggregator_tprol.f90 +++ b/mlprec/impl/aggregator/mld_z_dec_aggregator_tprol.f90 @@ -40,7 +40,7 @@ ! Subroutine: mld_z_dec_aggregator_tprol ! Version: complex ! -! This routine is mainly an interface to map_bld where the real work is performed. +! This routine is mainly an interface to soc_map_bld where the real work is performed. ! It takes care of some consistency checking, and calls map_to_tprol, which is ! refactored and shared among all the aggregation methods that produce a simple ! integer mapping. @@ -50,6 +50,8 @@ ! ag - type(mld_z_dec_aggregator_type), input/output. ! The aggregator object, carrying with itself the mapping algorithm. ! parms - The auxiliary parameters object +! ag_data - Auxiliary global aggregation parameters object +! ! ! a - type(psb_zspmat_type). ! The sparse matrix structure containing the local part of the @@ -72,13 +74,15 @@ ! info - integer, output. ! Error code. ! -subroutine mld_z_dec_aggregator_build_tprol(ag,parms,a,desc_a,ilaggr,nlaggr,op_prol,info) +subroutine mld_z_dec_aggregator_build_tprol(ag,parms,ag_data,& + & a,desc_a,ilaggr,nlaggr,op_prol,info) use psb_base_mod use mld_z_prec_type, mld_protect_name => mld_z_dec_aggregator_build_tprol use mld_z_inner_mod implicit none class(mld_z_dec_aggregator_type), target, intent(inout) :: ag type(mld_dml_parms), intent(inout) :: parms + type(mld_daggr_data), intent(in) :: ag_data type(psb_zspmat_type), intent(in) :: a type(psb_desc_type), intent(in) :: desc_a integer(psb_lpk_), allocatable, intent(out) :: ilaggr(:),nlaggr(:) @@ -91,6 +95,7 @@ subroutine mld_z_dec_aggregator_build_tprol(ag,parms,a,desc_a,ilaggr,nlaggr,op_ integer(psb_ipk_) :: err_act integer(psb_lpk_) :: ntaggr integer(psb_ipk_) :: debug_level, debug_unit + logical :: clean_zeros name='mld_z_dec_aggregator_tprol' call psb_erractionsave(err_act) @@ -111,13 +116,17 @@ subroutine mld_z_dec_aggregator_build_tprol(ag,parms,a,desc_a,ilaggr,nlaggr,op_ & mld_aggr_ord_nat_,is_legal_ml_aggr_ord) call mld_check_def(parms%aggr_thresh,'Aggr_Thresh',dzero,is_legal_d_aggr_thrs) - - call ag%map_bld(parms%aggr_ord,parms%aggr_thresh,a,desc_a,nlaggr,ilaggr,info) + ! + ! The decoupled aggregator based on SOC measures ignores + ! ag_data except for clean_zeros; soc_map_bld is a procedure pointer. + ! + clean_zeros = ag%do_clean_zeros + call ag%soc_map_bld(parms%aggr_ord,parms%aggr_thresh,clean_zeros,a,desc_a,nlaggr,ilaggr,info) if (info==psb_success_) call mld_map_to_tprol(desc_a,ilaggr,nlaggr,op_prol,info) if (info /= psb_success_) then info=psb_err_from_subroutine_ - call psb_errpush(info,name,a_err='map_bld/map_to_tprol') + call psb_errpush(info,name,a_err='soc_map_bld/map_to_tprol') goto 9999 endif diff --git a/mlprec/impl/aggregator/mld_z_soc1_map_bld.f90 b/mlprec/impl/aggregator/mld_z_soc1_map_bld.f90 index 249a5ec0..48e44213 100644 --- a/mlprec/impl/aggregator/mld_z_soc1_map_bld.f90 +++ b/mlprec/impl/aggregator/mld_z_soc1_map_bld.f90 @@ -67,7 +67,7 @@ ! ! ! -subroutine mld_z_soc1_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info) +subroutine mld_z_soc1_map_bld(iorder,theta,clean_zeros,a,desc_a,nlaggr,ilaggr,info) use psb_base_mod use mld_base_prec_type @@ -77,6 +77,7 @@ subroutine mld_z_soc1_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info) ! Arguments integer(psb_ipk_), intent(in) :: iorder + logical, intent(in) :: clean_zeros type(psb_zspmat_type), intent(in) :: a type(psb_desc_type), intent(in) :: desc_a real(psb_dpk_), intent(in) :: theta @@ -109,8 +110,8 @@ subroutine mld_z_soc1_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info) ! ictxt=desc_a%get_context() call psb_info(ictxt,me,np) - nrow = desc_a%get_local_rows() - ncol = desc_a%get_local_cols() + nrow = desc_a%get_local_rows() + ncol = desc_a%get_local_cols() nrglob = desc_a%get_global_rows() nr = a%get_nrows() @@ -132,6 +133,7 @@ subroutine mld_z_soc1_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info) end if call a%cp_to(acsr) + if (clean_zeros) call acsr%clean_zeros(info) if (iorder == mld_aggr_ord_nat_) then do i=1, nr ilaggr(i) = -(nr+1) @@ -169,12 +171,6 @@ subroutine mld_z_soc1_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info) icol(1:nz) = acsr%ja(acsr%irp(i):acsr%irp(i+1)-1) val(1:nz) = acsr%val(acsr%irp(i):acsr%irp(i+1)-1) -!!$ call a%csget(i,i,nz,irow,icol,val,info,chksz=.false.) -!!$ if (info /= psb_success_) then -!!$ info=psb_err_from_subroutine_ -!!$ call psb_errpush(info,name,a_err='csget') -!!$ goto 9999 -!!$ end if ! ! Build the set of all strongly coupled nodes @@ -194,7 +190,7 @@ subroutine mld_z_soc1_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info) ! If the whole strongly coupled neighborhood of I is ! as yet unconnected, turn it into the next aggregate. ! Same if ip==0 (in which case, neighborhood only - ! contains I even if it does not look from matrix) + ! contains I even if it does not look like it from matrix) ! disjoint = all(ilaggr(icol(1:ip)) == -(nr+1)).or.(ip==0) if (disjoint) then @@ -222,14 +218,10 @@ subroutine mld_z_soc1_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info) if (ilaggr(i) == -(nr+1)) then nz = (acsr%irp(i+1)-acsr%irp(i)) + if (nz == 1) cycle step2 icol(1:nz) = acsr%ja(acsr%irp(i):acsr%irp(i+1)-1) val(1:nz) = acsr%val(acsr%irp(i):acsr%irp(i+1)-1) -!!$ call a%csget(i,i,nz,irow,icol,val,info,chksz=.false.) -!!$ if (info /= psb_success_) then -!!$ info=psb_err_from_subroutine_ -!!$ call psb_errpush(info,name,a_err='psb_sp_getrow') -!!$ goto 9999 -!!$ end if + ! ! Find the most strongly connected neighbour that is ! already aggregated, if any, and join its aggregate @@ -261,6 +253,7 @@ subroutine mld_z_soc1_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info) if (ilaggr(i) < 0) then nz = (acsr%irp(i+1)-acsr%irp(i)) + if (nz == 1) cycle step3 icol(1:nz) = acsr%ja(acsr%irp(i):acsr%irp(i+1)-1) val(1:nz) = acsr%val(acsr%irp(i):acsr%irp(i+1)-1) ! @@ -288,8 +281,8 @@ subroutine mld_z_soc1_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info) end do else ! - ! This should not happen: we did not even connect with ourselves. - ! Create an isolate anyway. + ! This should not happen: we did not even connect with ourselves, + ! but it's not a singleton. ! naggr = naggr + 1 ilaggr(i) = naggr @@ -297,7 +290,6 @@ subroutine mld_z_soc1_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info) end if end do step3 - ! Any leftovers? do i=1, nr if (ilaggr(i) < 0) then diff --git a/mlprec/impl/aggregator/mld_z_soc2_map_bld.f90 b/mlprec/impl/aggregator/mld_z_soc2_map_bld.f90 index 967ba393..7e1538e2 100644 --- a/mlprec/impl/aggregator/mld_z_soc2_map_bld.f90 +++ b/mlprec/impl/aggregator/mld_z_soc2_map_bld.f90 @@ -66,7 +66,7 @@ ! ! ! -subroutine mld_z_soc2_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info) +subroutine mld_z_soc2_map_bld(iorder,theta,clean_zeros,a,desc_a,nlaggr,ilaggr,info) use psb_base_mod use mld_base_prec_type @@ -76,6 +76,7 @@ subroutine mld_z_soc2_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info) ! Arguments integer(psb_ipk_), intent(in) :: iorder + logical, intent(in) :: clean_zeros type(psb_zspmat_type), intent(in) :: a type(psb_desc_type), intent(in) :: desc_a real(psb_dpk_), intent(in) :: theta @@ -134,13 +135,14 @@ subroutine mld_z_soc2_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info) ! Phase zero: compute muij ! call a%cp_to(muij) + if (clean_zeros) call muij%clean_zeros(info) do i=1, nr do k=muij%irp(i),muij%irp(i+1)-1 j = muij%ja(k) if (j<= nr) muij%val(k) = abs(muij%val(k))/sqrt(abs(diag(i)*diag(j))) end do end do - !write(*,*) 'murows/cols ',maxval(mu_rows),maxval(mu_cols) + ! ! Compute the 1-neigbour; mark strong links with +1, weak links with -1 ! diff --git a/mlprec/impl/aggregator/mld_z_symdec_aggregator_tprol.f90 b/mlprec/impl/aggregator/mld_z_symdec_aggregator_tprol.f90 index 17c5c962..dbf11de1 100644 --- a/mlprec/impl/aggregator/mld_z_symdec_aggregator_tprol.f90 +++ b/mlprec/impl/aggregator/mld_z_symdec_aggregator_tprol.f90 @@ -41,7 +41,7 @@ ! Version: complex ! ! -! This routine is mainly an interface to map_bld where the real work is performed. +! This routine is mainly an interface to soc_map_bld where the real work is performed. ! It takes care of some consistency checking, and calls map_to_tprol, which is ! refactored and shared among all the aggregation methods that produce a simple ! integer mapping. It also symmetrizes the pattern of the local matrix A. @@ -53,6 +53,7 @@ ! ag - type(mld_z_dec_aggregator_type), input/output. ! The aggregator object, carrying with itself the mapping algorithm. ! parms - The auxiliary parameters object +! ag_data - Auxiliary global aggregation parameters object ! ! a - type(psb_zspmat_type). ! The sparse matrix structure containing the local part of the @@ -75,7 +76,8 @@ ! info - integer, output. ! Error code. ! -subroutine mld_z_symdec_aggregator_build_tprol(ag,parms,a,desc_a,ilaggr,nlaggr,op_prol,info) +subroutine mld_z_symdec_aggregator_build_tprol(ag,parms,ag_data,& + & a,desc_a,ilaggr,nlaggr,op_prol,info) use psb_base_mod use mld_z_prec_type use mld_z_symdec_aggregator_mod, mld_protect_name => mld_z_symdec_aggregator_build_tprol @@ -83,6 +85,7 @@ subroutine mld_z_symdec_aggregator_build_tprol(ag,parms,a,desc_a,ilaggr,nlaggr, implicit none class(mld_z_symdec_aggregator_type), target, intent(inout) :: ag type(mld_dml_parms), intent(inout) :: parms + type(mld_daggr_data), intent(in) :: ag_data type(psb_zspmat_type), intent(in) :: a type(psb_desc_type), intent(in) :: desc_a integer(psb_lpk_), allocatable, intent(out) :: ilaggr(:),nlaggr(:) @@ -97,6 +100,7 @@ subroutine mld_z_symdec_aggregator_build_tprol(ag,parms,a,desc_a,ilaggr,nlaggr, integer(psb_ipk_) :: nr integer(psb_lpk_) :: ntaggr integer(psb_ipk_) :: debug_level, debug_unit + logical :: clean_zeros name='mld_z_symdec_aggregator_tprol' call psb_erractionsave(err_act) @@ -130,14 +134,20 @@ subroutine mld_z_symdec_aggregator_build_tprol(ag,parms,a,desc_a,ilaggr,nlaggr, if (info == psb_success_) call atrans%free() if (info == psb_success_) call atmp%cscnv(info,type='CSR') + ! + ! The decoupled aggregator based on SOC measures ignores + ! ag_data except for clean_zeros; soc_map_bld is a procedure pointer. + ! + clean_zeros = ag%do_clean_zeros if (info == psb_success_) & - & call ag%map_bld(parms%aggr_ord,parms%aggr_thresh,atmp,desc_a,nlaggr,ilaggr,info) + & call ag%soc_map_bld(parms%aggr_ord,parms%aggr_thresh,clean_zeros,atmp,& + & desc_a,nlaggr,ilaggr,info) if (info == psb_success_) call atmp%free() if (info == psb_success_) call mld_map_to_tprol(desc_a,ilaggr,nlaggr,op_prol,info) if (info /= psb_success_) then info=psb_err_from_subroutine_ - call psb_errpush(info,name,a_err='map_bld/map_to_tprol') + call psb_errpush(info,name,a_err='soc_map_bld/map_to_tprol') goto 9999 endif diff --git a/mlprec/impl/level/mld_c_base_onelev_cnv.f90 b/mlprec/impl/level/mld_c_base_onelev_cnv.f90 index e3d49b1b..cdc48dad 100644 --- a/mlprec/impl/level/mld_c_base_onelev_cnv.f90 +++ b/mlprec/impl/level/mld_c_base_onelev_cnv.f90 @@ -62,6 +62,6 @@ subroutine mld_c_base_onelev_cnv(lv,info,amold,vmold,imold) & call lv%ac%cscnv(info,mold=amold) if (info == psb_success_ .and. lv%desc_ac%is_ok() & & .and. present(imold)) call lv%desc_ac%cnv(imold) - call lv%map%cnv(info,mold=amold,imold=imold) + if (info == psb_success_) call lv%map%cnv(info,mold=amold,imold=imold) end if end subroutine mld_c_base_onelev_cnv diff --git a/mlprec/impl/level/mld_c_base_onelev_dump.f90 b/mlprec/impl/level/mld_c_base_onelev_dump.f90 index 5416dd63..1fa45bb2 100644 --- a/mlprec/impl/level/mld_c_base_onelev_dump.f90 +++ b/mlprec/impl/level/mld_c_base_onelev_dump.f90 @@ -36,7 +36,7 @@ ! ! subroutine mld_c_base_onelev_dump(lv,level,info,prefix,head,ac,rp,& - & smoother,solver,global_num) + & smoother,solver,tprol,global_num) use psb_base_mod use mld_c_onelev_mod, mld_protect_name => mld_c_base_onelev_dump @@ -45,13 +45,13 @@ subroutine mld_c_base_onelev_dump(lv,level,info,prefix,head,ac,rp,& integer(psb_ipk_), intent(in) :: level integer(psb_ipk_), intent(out) :: info character(len=*), intent(in), optional :: prefix, head - logical, optional, intent(in) :: ac, rp, smoother, solver,global_num + logical, optional, intent(in) :: ac, rp, smoother, solver, tprol, global_num ! Local variables integer(psb_ipk_) :: i, j, il1, iln, lname, lev integer(psb_ipk_) :: icontxt,iam, np character(len=80) :: prefix_ character(len=120) :: fname ! len should be at least 20 more than - logical :: ac_, rp_,global_num_ + logical :: ac_, rp_, tprol_, global_num_ integer(psb_ipk_), allocatable :: ivr(:), ivc(:) info = 0 @@ -79,6 +79,11 @@ subroutine mld_c_base_onelev_dump(lv,level,info,prefix,head,ac,rp,& else rp_ = .false. end if + if (present(tprol)) then + tprol_ = tprol + else + tprol_ = .false. + end if if (present(global_num)) then global_num_ = global_num else @@ -104,6 +109,15 @@ subroutine mld_c_base_onelev_dump(lv,level,info,prefix,head,ac,rp,& write(fname(lname+1:),'(a,i3.3,a)')'_l',level,'_p.mtx' call lv%map%map_Y2X%print(fname,head=head,ivr=ivr,ivc=ivc) end if + if (tprol_) then + ! Tentative prolongator is stored with column indices already + ! in global numbering, so only IVR is needed. + ivr = lv%map%p_desc_X%get_global_indices(owned=.false.) + write(fname(lname+1:),'(a,i3.3,a)')'_l',level,'_tprol.mtx' + ! + ! This is not implemented yet. + !call lv%tprol%print(fname,head=head,ivr=ivr) + end if end if else if (level >= 2) then @@ -117,6 +131,12 @@ subroutine mld_c_base_onelev_dump(lv,level,info,prefix,head,ac,rp,& write(fname(lname+1:),'(a,i3.3,a)')'_l',level,'_p.mtx' call lv%map%map_Y2X%print(fname,head=head) end if + if (tprol_) then + write(fname(lname+1:),'(a,i3.3,a)')'_l',level,'_tprol.mtx' + ! + ! This is not implemented yet. + !call lv%tprol%print(fname,head=head) + end if end if end if diff --git a/mlprec/impl/level/mld_c_base_onelev_setag.f90 b/mlprec/impl/level/mld_c_base_onelev_setag.f90 index 25da1459..37b0d1de 100644 --- a/mlprec/impl/level/mld_c_base_onelev_setag.f90 +++ b/mlprec/impl/level/mld_c_base_onelev_setag.f90 @@ -75,6 +75,7 @@ subroutine mld_c_base_onelev_setag(lv,val,info,pos) end if lv%parms%par_aggr_alg = mld_ext_aggr_ lv%parms%aggr_type = mld_noalg_ + call lv%aggr%default() end if end subroutine mld_c_base_onelev_setag diff --git a/mlprec/impl/level/mld_d_base_onelev_cnv.f90 b/mlprec/impl/level/mld_d_base_onelev_cnv.f90 index 23ff2c26..a43054a6 100644 --- a/mlprec/impl/level/mld_d_base_onelev_cnv.f90 +++ b/mlprec/impl/level/mld_d_base_onelev_cnv.f90 @@ -62,6 +62,6 @@ subroutine mld_d_base_onelev_cnv(lv,info,amold,vmold,imold) & call lv%ac%cscnv(info,mold=amold) if (info == psb_success_ .and. lv%desc_ac%is_ok() & & .and. present(imold)) call lv%desc_ac%cnv(imold) - call lv%map%cnv(info,mold=amold,imold=imold) + if (info == psb_success_) call lv%map%cnv(info,mold=amold,imold=imold) end if end subroutine mld_d_base_onelev_cnv diff --git a/mlprec/impl/level/mld_d_base_onelev_dump.f90 b/mlprec/impl/level/mld_d_base_onelev_dump.f90 index 36b59197..5522de4e 100644 --- a/mlprec/impl/level/mld_d_base_onelev_dump.f90 +++ b/mlprec/impl/level/mld_d_base_onelev_dump.f90 @@ -36,7 +36,7 @@ ! ! subroutine mld_d_base_onelev_dump(lv,level,info,prefix,head,ac,rp,& - & smoother,solver,global_num) + & smoother,solver,tprol,global_num) use psb_base_mod use mld_d_onelev_mod, mld_protect_name => mld_d_base_onelev_dump @@ -45,13 +45,13 @@ subroutine mld_d_base_onelev_dump(lv,level,info,prefix,head,ac,rp,& integer(psb_ipk_), intent(in) :: level integer(psb_ipk_), intent(out) :: info character(len=*), intent(in), optional :: prefix, head - logical, optional, intent(in) :: ac, rp, smoother, solver,global_num + logical, optional, intent(in) :: ac, rp, smoother, solver, tprol, global_num ! Local variables integer(psb_ipk_) :: i, j, il1, iln, lname, lev integer(psb_ipk_) :: icontxt,iam, np character(len=80) :: prefix_ character(len=120) :: fname ! len should be at least 20 more than - logical :: ac_, rp_,global_num_ + logical :: ac_, rp_, tprol_, global_num_ integer(psb_ipk_), allocatable :: ivr(:), ivc(:) info = 0 @@ -79,6 +79,11 @@ subroutine mld_d_base_onelev_dump(lv,level,info,prefix,head,ac,rp,& else rp_ = .false. end if + if (present(tprol)) then + tprol_ = tprol + else + tprol_ = .false. + end if if (present(global_num)) then global_num_ = global_num else @@ -104,6 +109,15 @@ subroutine mld_d_base_onelev_dump(lv,level,info,prefix,head,ac,rp,& write(fname(lname+1:),'(a,i3.3,a)')'_l',level,'_p.mtx' call lv%map%map_Y2X%print(fname,head=head,ivr=ivr,ivc=ivc) end if + if (tprol_) then + ! Tentative prolongator is stored with column indices already + ! in global numbering, so only IVR is needed. + ivr = lv%map%p_desc_X%get_global_indices(owned=.false.) + write(fname(lname+1:),'(a,i3.3,a)')'_l',level,'_tprol.mtx' + ! + ! This is not implemented yet. + !call lv%tprol%print(fname,head=head,ivr=ivr) + end if end if else if (level >= 2) then @@ -117,6 +131,12 @@ subroutine mld_d_base_onelev_dump(lv,level,info,prefix,head,ac,rp,& write(fname(lname+1:),'(a,i3.3,a)')'_l',level,'_p.mtx' call lv%map%map_Y2X%print(fname,head=head) end if + if (tprol_) then + write(fname(lname+1:),'(a,i3.3,a)')'_l',level,'_tprol.mtx' + ! + ! This is not implemented yet. + !call lv%tprol%print(fname,head=head) + end if end if end if diff --git a/mlprec/impl/level/mld_d_base_onelev_setag.f90 b/mlprec/impl/level/mld_d_base_onelev_setag.f90 index c4c26cca..83bd9f76 100644 --- a/mlprec/impl/level/mld_d_base_onelev_setag.f90 +++ b/mlprec/impl/level/mld_d_base_onelev_setag.f90 @@ -75,6 +75,7 @@ subroutine mld_d_base_onelev_setag(lv,val,info,pos) end if lv%parms%par_aggr_alg = mld_ext_aggr_ lv%parms%aggr_type = mld_noalg_ + call lv%aggr%default() end if end subroutine mld_d_base_onelev_setag diff --git a/mlprec/impl/level/mld_s_base_onelev_cnv.f90 b/mlprec/impl/level/mld_s_base_onelev_cnv.f90 index 89a3d087..732681b7 100644 --- a/mlprec/impl/level/mld_s_base_onelev_cnv.f90 +++ b/mlprec/impl/level/mld_s_base_onelev_cnv.f90 @@ -62,6 +62,6 @@ subroutine mld_s_base_onelev_cnv(lv,info,amold,vmold,imold) & call lv%ac%cscnv(info,mold=amold) if (info == psb_success_ .and. lv%desc_ac%is_ok() & & .and. present(imold)) call lv%desc_ac%cnv(imold) - call lv%map%cnv(info,mold=amold,imold=imold) + if (info == psb_success_) call lv%map%cnv(info,mold=amold,imold=imold) end if end subroutine mld_s_base_onelev_cnv diff --git a/mlprec/impl/level/mld_s_base_onelev_dump.f90 b/mlprec/impl/level/mld_s_base_onelev_dump.f90 index df93abe4..08792813 100644 --- a/mlprec/impl/level/mld_s_base_onelev_dump.f90 +++ b/mlprec/impl/level/mld_s_base_onelev_dump.f90 @@ -36,7 +36,7 @@ ! ! subroutine mld_s_base_onelev_dump(lv,level,info,prefix,head,ac,rp,& - & smoother,solver,global_num) + & smoother,solver,tprol,global_num) use psb_base_mod use mld_s_onelev_mod, mld_protect_name => mld_s_base_onelev_dump @@ -45,13 +45,13 @@ subroutine mld_s_base_onelev_dump(lv,level,info,prefix,head,ac,rp,& integer(psb_ipk_), intent(in) :: level integer(psb_ipk_), intent(out) :: info character(len=*), intent(in), optional :: prefix, head - logical, optional, intent(in) :: ac, rp, smoother, solver,global_num + logical, optional, intent(in) :: ac, rp, smoother, solver, tprol, global_num ! Local variables integer(psb_ipk_) :: i, j, il1, iln, lname, lev integer(psb_ipk_) :: icontxt,iam, np character(len=80) :: prefix_ character(len=120) :: fname ! len should be at least 20 more than - logical :: ac_, rp_,global_num_ + logical :: ac_, rp_, tprol_, global_num_ integer(psb_ipk_), allocatable :: ivr(:), ivc(:) info = 0 @@ -79,6 +79,11 @@ subroutine mld_s_base_onelev_dump(lv,level,info,prefix,head,ac,rp,& else rp_ = .false. end if + if (present(tprol)) then + tprol_ = tprol + else + tprol_ = .false. + end if if (present(global_num)) then global_num_ = global_num else @@ -104,6 +109,15 @@ subroutine mld_s_base_onelev_dump(lv,level,info,prefix,head,ac,rp,& write(fname(lname+1:),'(a,i3.3,a)')'_l',level,'_p.mtx' call lv%map%map_Y2X%print(fname,head=head,ivr=ivr,ivc=ivc) end if + if (tprol_) then + ! Tentative prolongator is stored with column indices already + ! in global numbering, so only IVR is needed. + ivr = lv%map%p_desc_X%get_global_indices(owned=.false.) + write(fname(lname+1:),'(a,i3.3,a)')'_l',level,'_tprol.mtx' + ! + ! This is not implemented yet. + !call lv%tprol%print(fname,head=head,ivr=ivr) + end if end if else if (level >= 2) then @@ -117,6 +131,12 @@ subroutine mld_s_base_onelev_dump(lv,level,info,prefix,head,ac,rp,& write(fname(lname+1:),'(a,i3.3,a)')'_l',level,'_p.mtx' call lv%map%map_Y2X%print(fname,head=head) end if + if (tprol_) then + write(fname(lname+1:),'(a,i3.3,a)')'_l',level,'_tprol.mtx' + ! + ! This is not implemented yet. + !call lv%tprol%print(fname,head=head) + end if end if end if diff --git a/mlprec/impl/level/mld_s_base_onelev_setag.f90 b/mlprec/impl/level/mld_s_base_onelev_setag.f90 index 36bf2871..e9de3cbc 100644 --- a/mlprec/impl/level/mld_s_base_onelev_setag.f90 +++ b/mlprec/impl/level/mld_s_base_onelev_setag.f90 @@ -75,6 +75,7 @@ subroutine mld_s_base_onelev_setag(lv,val,info,pos) end if lv%parms%par_aggr_alg = mld_ext_aggr_ lv%parms%aggr_type = mld_noalg_ + call lv%aggr%default() end if end subroutine mld_s_base_onelev_setag diff --git a/mlprec/impl/level/mld_z_base_onelev_cnv.f90 b/mlprec/impl/level/mld_z_base_onelev_cnv.f90 index 407fca8d..11614681 100644 --- a/mlprec/impl/level/mld_z_base_onelev_cnv.f90 +++ b/mlprec/impl/level/mld_z_base_onelev_cnv.f90 @@ -62,6 +62,6 @@ subroutine mld_z_base_onelev_cnv(lv,info,amold,vmold,imold) & call lv%ac%cscnv(info,mold=amold) if (info == psb_success_ .and. lv%desc_ac%is_ok() & & .and. present(imold)) call lv%desc_ac%cnv(imold) - call lv%map%cnv(info,mold=amold,imold=imold) + if (info == psb_success_) call lv%map%cnv(info,mold=amold,imold=imold) end if end subroutine mld_z_base_onelev_cnv diff --git a/mlprec/impl/level/mld_z_base_onelev_dump.f90 b/mlprec/impl/level/mld_z_base_onelev_dump.f90 index c1181e97..8e96b373 100644 --- a/mlprec/impl/level/mld_z_base_onelev_dump.f90 +++ b/mlprec/impl/level/mld_z_base_onelev_dump.f90 @@ -36,7 +36,7 @@ ! ! subroutine mld_z_base_onelev_dump(lv,level,info,prefix,head,ac,rp,& - & smoother,solver,global_num) + & smoother,solver,tprol,global_num) use psb_base_mod use mld_z_onelev_mod, mld_protect_name => mld_z_base_onelev_dump @@ -45,13 +45,13 @@ subroutine mld_z_base_onelev_dump(lv,level,info,prefix,head,ac,rp,& integer(psb_ipk_), intent(in) :: level integer(psb_ipk_), intent(out) :: info character(len=*), intent(in), optional :: prefix, head - logical, optional, intent(in) :: ac, rp, smoother, solver,global_num + logical, optional, intent(in) :: ac, rp, smoother, solver, tprol, global_num ! Local variables integer(psb_ipk_) :: i, j, il1, iln, lname, lev integer(psb_ipk_) :: icontxt,iam, np character(len=80) :: prefix_ character(len=120) :: fname ! len should be at least 20 more than - logical :: ac_, rp_,global_num_ + logical :: ac_, rp_, tprol_, global_num_ integer(psb_ipk_), allocatable :: ivr(:), ivc(:) info = 0 @@ -79,6 +79,11 @@ subroutine mld_z_base_onelev_dump(lv,level,info,prefix,head,ac,rp,& else rp_ = .false. end if + if (present(tprol)) then + tprol_ = tprol + else + tprol_ = .false. + end if if (present(global_num)) then global_num_ = global_num else @@ -104,6 +109,15 @@ subroutine mld_z_base_onelev_dump(lv,level,info,prefix,head,ac,rp,& write(fname(lname+1:),'(a,i3.3,a)')'_l',level,'_p.mtx' call lv%map%map_Y2X%print(fname,head=head,ivr=ivr,ivc=ivc) end if + if (tprol_) then + ! Tentative prolongator is stored with column indices already + ! in global numbering, so only IVR is needed. + ivr = lv%map%p_desc_X%get_global_indices(owned=.false.) + write(fname(lname+1:),'(a,i3.3,a)')'_l',level,'_tprol.mtx' + ! + ! This is not implemented yet. + !call lv%tprol%print(fname,head=head,ivr=ivr) + end if end if else if (level >= 2) then @@ -117,6 +131,12 @@ subroutine mld_z_base_onelev_dump(lv,level,info,prefix,head,ac,rp,& write(fname(lname+1:),'(a,i3.3,a)')'_l',level,'_p.mtx' call lv%map%map_Y2X%print(fname,head=head) end if + if (tprol_) then + write(fname(lname+1:),'(a,i3.3,a)')'_l',level,'_tprol.mtx' + ! + ! This is not implemented yet. + !call lv%tprol%print(fname,head=head) + end if end if end if diff --git a/mlprec/impl/level/mld_z_base_onelev_setag.f90 b/mlprec/impl/level/mld_z_base_onelev_setag.f90 index e6c7d51b..84bf26be 100644 --- a/mlprec/impl/level/mld_z_base_onelev_setag.f90 +++ b/mlprec/impl/level/mld_z_base_onelev_setag.f90 @@ -75,6 +75,7 @@ subroutine mld_z_base_onelev_setag(lv,val,info,pos) end if lv%parms%par_aggr_alg = mld_ext_aggr_ lv%parms%aggr_type = mld_noalg_ + call lv%aggr%default() end if end subroutine mld_z_base_onelev_setag diff --git a/mlprec/impl/mld_c_extprol_bld.F90 b/mlprec/impl/mld_c_extprol_bld.F90 index db934bfb..d52ce68a 100644 --- a/mlprec/impl/mld_c_extprol_bld.F90 +++ b/mlprec/impl/mld_c_extprol_bld.F90 @@ -160,9 +160,9 @@ subroutine mld_c_extprol_bld(a,desc_a,p,prolv,restrv,info,amold,vmold,imold) ! Check to ensure all procs have the same ! newsz = -1 - mxplevs = p%max_levs - mnaggratio = p%min_cr_ratio - casize = p%min_coarse_size + mxplevs = p%ag_data%max_levs + mnaggratio = p%ag_data%min_cr_ratio + casize = p%ag_data%min_coarse_size iszv = size(p%precv) nprolv = size(prolv) nrestrv = size(restrv) @@ -172,17 +172,17 @@ subroutine mld_c_extprol_bld(a,desc_a,p,prolv,restrv,info,amold,vmold,imold) call psb_bcast(ictxt,mnaggratio) call psb_bcast(ictxt,nprolv) call psb_bcast(ictxt,nrestrv) - if (casize /= p%min_coarse_size) then + if (casize /= p%ag_data%min_coarse_size) then info=psb_err_internal_error_ call psb_errpush(info,name,a_err='Inconsistent min_coarse_size') goto 9999 end if - if (mxplevs /= p%max_levs) then + if (mxplevs /= p%ag_data%max_levs) then info=psb_err_internal_error_ call psb_errpush(info,name,a_err='Inconsistent max_levs') goto 9999 end if - if (mnaggratio /= p%min_cr_ratio) then + if (mnaggratio /= p%ag_data%min_cr_ratio) then info=psb_err_internal_error_ call psb_errpush(info,name,a_err='Inconsistent min_cr_ratio') goto 9999 @@ -226,7 +226,7 @@ subroutine mld_c_extprol_bld(a,desc_a,p,prolv,restrv,info,amold,vmold,imold) ! nplevs = nrestrv + 1 - p%max_levs = nplevs + p%ag_data%max_levs = nplevs ! ! Fixed number of levels. diff --git a/mlprec/impl/mld_c_hierarchy_bld.f90 b/mlprec/impl/mld_c_hierarchy_bld.f90 index 08079ab2..56c3057c 100644 --- a/mlprec/impl/mld_c_hierarchy_bld.f90 +++ b/mlprec/impl/mld_c_hierarchy_bld.f90 @@ -124,25 +124,25 @@ subroutine mld_c_hierarchy_bld(a,desc_a,prec,info) ! Check to ensure all procs have the same ! newsz = -1 - mxplevs = prec%max_levs - mnaggratio = prec%min_cr_ratio - casize = prec%min_coarse_size + mxplevs = prec%ag_data%max_levs + mnaggratio = prec%ag_data%min_cr_ratio + casize = prec%ag_data%min_coarse_size iszv = size(prec%precv) call psb_bcast(ictxt,iszv) call psb_bcast(ictxt,casize) call psb_bcast(ictxt,mxplevs) call psb_bcast(ictxt,mnaggratio) - if (casize /= prec%min_coarse_size) then + if (casize /= prec%ag_data%min_coarse_size) then info=psb_err_internal_error_ call psb_errpush(info,name,a_err='Inconsistent min_coarse_size') goto 9999 end if - if (mxplevs /= prec%max_levs) then + if (mxplevs /= prec%ag_data%max_levs) then info=psb_err_internal_error_ call psb_errpush(info,name,a_err='Inconsistent max_levs') goto 9999 end if - if (mnaggratio /= prec%min_cr_ratio) then + if (mnaggratio /= prec%ag_data%min_cr_ratio) then info=psb_err_internal_error_ call psb_errpush(info,name,a_err='Inconsistent min_cr_ratio') goto 9999 @@ -196,7 +196,7 @@ subroutine mld_c_hierarchy_bld(a,desc_a,prec,info) casize = max(casize,ione) casize = casize*40_psb_ipk_ call psb_bcast(ictxt,casize) - prec%min_coarse_size = casize + prec%ag_data%min_coarse_size = casize end if nplevs = max(itwo,mxplevs) @@ -298,7 +298,7 @@ subroutine mld_c_hierarchy_bld(a,desc_a,prec,info) if (info == psb_success_)& & call prec%precv(i)%bld_tprol(prec%precv(i-1)%base_a,& & prec%precv(i-1)%base_desc,& - & ilaggr,nlaggr,op_prol,info) + & ilaggr,nlaggr,op_prol,prec%ag_data,info) if (info /= psb_success_) then call psb_errpush(psb_err_internal_error_,name,& diff --git a/mlprec/impl/mld_ccprecset.F90 b/mlprec/impl/mld_ccprecset.F90 index 04684638..bd4f0985 100644 --- a/mlprec/impl/mld_ccprecset.F90 +++ b/mlprec/impl/mld_ccprecset.F90 @@ -129,7 +129,7 @@ subroutine mld_ccprecseti(p,what,val,info,ilev,ilmax,pos,idx) end if else ilev_ = 1 - ilmax_ = ilev_ + ilmax_ = nlev_ end if if ((ilev_<1).or.(ilev_ > nlev_)) then info = -1 @@ -146,10 +146,10 @@ subroutine mld_ccprecseti(p,what,val,info,ilev,ilmax,pos,idx) select case(psb_toupper(what)) case ('MIN_COARSE_SIZE') - p%min_coarse_size = max(val,-1) + p%ag_data%min_coarse_size = max(val,-1) return case('MAX_LEVS') - p%max_levs = max(val,1) + p%ag_data%max_levs = max(val,1) return case ('OUTER_SWEEPS') p%outer_sweeps = max(val,1) @@ -502,7 +502,7 @@ subroutine mld_ccprecsetc(p,what,string,info,ilev,ilmax,pos,idx) end if else ilev_ = 1 - ilmax_ = ilev_ + ilmax_ = nlev_ end if if ((ilev_<1).or.(ilev_ > nlev_)) then info = -1 @@ -593,7 +593,7 @@ subroutine mld_ccprecsetr(p,what,val,info,ilev,ilmax,pos,idx) select case(psb_toupper(what)) case ('MIN_CR_RATIO') - p%min_cr_ratio = max(sone,val) + p%ag_data%min_cr_ratio = max(sone,val) return end select @@ -615,7 +615,7 @@ subroutine mld_ccprecsetr(p,what,val,info,ilev,ilmax,pos,idx) end if else ilev_ = 1 - ilmax_ = ilev_ + ilmax_ = nlev_ end if if ((ilev_<1).or.(ilev_ > nlev_)) then info = -1 diff --git a/mlprec/impl/mld_cprecinit.F90 b/mlprec/impl/mld_cprecinit.F90 index 9574ab0c..837b4260 100644 --- a/mlprec/impl/mld_cprecinit.F90 +++ b/mlprec/impl/mld_cprecinit.F90 @@ -119,7 +119,7 @@ subroutine mld_cprecinit(ictxt,prec,ptype,info) endif endif prec%ictxt = ictxt - prec%min_coarse_size = -1 + prec%ag_data%min_coarse_size = -1 select case(psb_toupper(ptype(1:len_trim(ptype)))) case ('NOPREC','NONE') @@ -186,7 +186,7 @@ subroutine mld_cprecinit(ictxt,prec,ptype,info) case ('ML') - nlev_ = prec%max_levs + nlev_ = prec%ag_data%max_levs ilev_ = 1 allocate(prec%precv(nlev_),stat=info) diff --git a/mlprec/impl/mld_d_extprol_bld.F90 b/mlprec/impl/mld_d_extprol_bld.F90 index a0be249e..aff2d843 100644 --- a/mlprec/impl/mld_d_extprol_bld.F90 +++ b/mlprec/impl/mld_d_extprol_bld.F90 @@ -160,9 +160,9 @@ subroutine mld_d_extprol_bld(a,desc_a,p,prolv,restrv,info,amold,vmold,imold) ! Check to ensure all procs have the same ! newsz = -1 - mxplevs = p%max_levs - mnaggratio = p%min_cr_ratio - casize = p%min_coarse_size + mxplevs = p%ag_data%max_levs + mnaggratio = p%ag_data%min_cr_ratio + casize = p%ag_data%min_coarse_size iszv = size(p%precv) nprolv = size(prolv) nrestrv = size(restrv) @@ -172,17 +172,17 @@ subroutine mld_d_extprol_bld(a,desc_a,p,prolv,restrv,info,amold,vmold,imold) call psb_bcast(ictxt,mnaggratio) call psb_bcast(ictxt,nprolv) call psb_bcast(ictxt,nrestrv) - if (casize /= p%min_coarse_size) then + if (casize /= p%ag_data%min_coarse_size) then info=psb_err_internal_error_ call psb_errpush(info,name,a_err='Inconsistent min_coarse_size') goto 9999 end if - if (mxplevs /= p%max_levs) then + if (mxplevs /= p%ag_data%max_levs) then info=psb_err_internal_error_ call psb_errpush(info,name,a_err='Inconsistent max_levs') goto 9999 end if - if (mnaggratio /= p%min_cr_ratio) then + if (mnaggratio /= p%ag_data%min_cr_ratio) then info=psb_err_internal_error_ call psb_errpush(info,name,a_err='Inconsistent min_cr_ratio') goto 9999 @@ -226,7 +226,7 @@ subroutine mld_d_extprol_bld(a,desc_a,p,prolv,restrv,info,amold,vmold,imold) ! nplevs = nrestrv + 1 - p%max_levs = nplevs + p%ag_data%max_levs = nplevs ! ! Fixed number of levels. diff --git a/mlprec/impl/mld_d_hierarchy_bld.f90 b/mlprec/impl/mld_d_hierarchy_bld.f90 index 655a2000..40e8f216 100644 --- a/mlprec/impl/mld_d_hierarchy_bld.f90 +++ b/mlprec/impl/mld_d_hierarchy_bld.f90 @@ -124,25 +124,25 @@ subroutine mld_d_hierarchy_bld(a,desc_a,prec,info) ! Check to ensure all procs have the same ! newsz = -1 - mxplevs = prec%max_levs - mnaggratio = prec%min_cr_ratio - casize = prec%min_coarse_size + mxplevs = prec%ag_data%max_levs + mnaggratio = prec%ag_data%min_cr_ratio + casize = prec%ag_data%min_coarse_size iszv = size(prec%precv) call psb_bcast(ictxt,iszv) call psb_bcast(ictxt,casize) call psb_bcast(ictxt,mxplevs) call psb_bcast(ictxt,mnaggratio) - if (casize /= prec%min_coarse_size) then + if (casize /= prec%ag_data%min_coarse_size) then info=psb_err_internal_error_ call psb_errpush(info,name,a_err='Inconsistent min_coarse_size') goto 9999 end if - if (mxplevs /= prec%max_levs) then + if (mxplevs /= prec%ag_data%max_levs) then info=psb_err_internal_error_ call psb_errpush(info,name,a_err='Inconsistent max_levs') goto 9999 end if - if (mnaggratio /= prec%min_cr_ratio) then + if (mnaggratio /= prec%ag_data%min_cr_ratio) then info=psb_err_internal_error_ call psb_errpush(info,name,a_err='Inconsistent min_cr_ratio') goto 9999 @@ -196,7 +196,7 @@ subroutine mld_d_hierarchy_bld(a,desc_a,prec,info) casize = max(casize,ione) casize = casize*40_psb_ipk_ call psb_bcast(ictxt,casize) - prec%min_coarse_size = casize + prec%ag_data%min_coarse_size = casize end if nplevs = max(itwo,mxplevs) @@ -298,7 +298,7 @@ subroutine mld_d_hierarchy_bld(a,desc_a,prec,info) if (info == psb_success_)& & call prec%precv(i)%bld_tprol(prec%precv(i-1)%base_a,& & prec%precv(i-1)%base_desc,& - & ilaggr,nlaggr,op_prol,info) + & ilaggr,nlaggr,op_prol,prec%ag_data,info) if (info /= psb_success_) then call psb_errpush(psb_err_internal_error_,name,& diff --git a/mlprec/impl/mld_dcprecset.F90 b/mlprec/impl/mld_dcprecset.F90 index 2532504f..0338d640 100644 --- a/mlprec/impl/mld_dcprecset.F90 +++ b/mlprec/impl/mld_dcprecset.F90 @@ -135,7 +135,7 @@ subroutine mld_dcprecseti(p,what,val,info,ilev,ilmax,pos,idx) end if else ilev_ = 1 - ilmax_ = ilev_ + ilmax_ = nlev_ end if if ((ilev_<1).or.(ilev_ > nlev_)) then info = -1 @@ -152,10 +152,10 @@ subroutine mld_dcprecseti(p,what,val,info,ilev,ilmax,pos,idx) select case(psb_toupper(what)) case ('MIN_COARSE_SIZE') - p%min_coarse_size = max(val,-1) + p%ag_data%min_coarse_size = max(val,-1) return case('MAX_LEVS') - p%max_levs = max(val,1) + p%ag_data%max_levs = max(val,1) return case ('OUTER_SWEEPS') p%outer_sweeps = max(val,1) @@ -536,7 +536,7 @@ subroutine mld_dcprecsetc(p,what,string,info,ilev,ilmax,pos,idx) end if else ilev_ = 1 - ilmax_ = ilev_ + ilmax_ = nlev_ end if if ((ilev_<1).or.(ilev_ > nlev_)) then info = -1 @@ -627,7 +627,7 @@ subroutine mld_dcprecsetr(p,what,val,info,ilev,ilmax,pos,idx) select case(psb_toupper(what)) case ('MIN_CR_RATIO') - p%min_cr_ratio = max(done,val) + p%ag_data%min_cr_ratio = max(done,val) return end select @@ -649,7 +649,7 @@ subroutine mld_dcprecsetr(p,what,val,info,ilev,ilmax,pos,idx) end if else ilev_ = 1 - ilmax_ = ilev_ + ilmax_ = nlev_ end if if ((ilev_<1).or.(ilev_ > nlev_)) then info = -1 diff --git a/mlprec/impl/mld_dprecinit.F90 b/mlprec/impl/mld_dprecinit.F90 index f9b24687..6703dddf 100644 --- a/mlprec/impl/mld_dprecinit.F90 +++ b/mlprec/impl/mld_dprecinit.F90 @@ -122,7 +122,7 @@ subroutine mld_dprecinit(ictxt,prec,ptype,info) endif endif prec%ictxt = ictxt - prec%min_coarse_size = -1 + prec%ag_data%min_coarse_size = -1 select case(psb_toupper(ptype(1:len_trim(ptype)))) case ('NOPREC','NONE') @@ -189,7 +189,7 @@ subroutine mld_dprecinit(ictxt,prec,ptype,info) case ('ML') - nlev_ = prec%max_levs + nlev_ = prec%ag_data%max_levs ilev_ = 1 allocate(prec%precv(nlev_),stat=info) diff --git a/mlprec/impl/mld_s_extprol_bld.F90 b/mlprec/impl/mld_s_extprol_bld.F90 index 96b4d8d8..9ee6755c 100644 --- a/mlprec/impl/mld_s_extprol_bld.F90 +++ b/mlprec/impl/mld_s_extprol_bld.F90 @@ -160,9 +160,9 @@ subroutine mld_s_extprol_bld(a,desc_a,p,prolv,restrv,info,amold,vmold,imold) ! Check to ensure all procs have the same ! newsz = -1 - mxplevs = p%max_levs - mnaggratio = p%min_cr_ratio - casize = p%min_coarse_size + mxplevs = p%ag_data%max_levs + mnaggratio = p%ag_data%min_cr_ratio + casize = p%ag_data%min_coarse_size iszv = size(p%precv) nprolv = size(prolv) nrestrv = size(restrv) @@ -172,17 +172,17 @@ subroutine mld_s_extprol_bld(a,desc_a,p,prolv,restrv,info,amold,vmold,imold) call psb_bcast(ictxt,mnaggratio) call psb_bcast(ictxt,nprolv) call psb_bcast(ictxt,nrestrv) - if (casize /= p%min_coarse_size) then + if (casize /= p%ag_data%min_coarse_size) then info=psb_err_internal_error_ call psb_errpush(info,name,a_err='Inconsistent min_coarse_size') goto 9999 end if - if (mxplevs /= p%max_levs) then + if (mxplevs /= p%ag_data%max_levs) then info=psb_err_internal_error_ call psb_errpush(info,name,a_err='Inconsistent max_levs') goto 9999 end if - if (mnaggratio /= p%min_cr_ratio) then + if (mnaggratio /= p%ag_data%min_cr_ratio) then info=psb_err_internal_error_ call psb_errpush(info,name,a_err='Inconsistent min_cr_ratio') goto 9999 @@ -226,7 +226,7 @@ subroutine mld_s_extprol_bld(a,desc_a,p,prolv,restrv,info,amold,vmold,imold) ! nplevs = nrestrv + 1 - p%max_levs = nplevs + p%ag_data%max_levs = nplevs ! ! Fixed number of levels. diff --git a/mlprec/impl/mld_s_hierarchy_bld.f90 b/mlprec/impl/mld_s_hierarchy_bld.f90 index a1f63975..e159b368 100644 --- a/mlprec/impl/mld_s_hierarchy_bld.f90 +++ b/mlprec/impl/mld_s_hierarchy_bld.f90 @@ -124,25 +124,25 @@ subroutine mld_s_hierarchy_bld(a,desc_a,prec,info) ! Check to ensure all procs have the same ! newsz = -1 - mxplevs = prec%max_levs - mnaggratio = prec%min_cr_ratio - casize = prec%min_coarse_size + mxplevs = prec%ag_data%max_levs + mnaggratio = prec%ag_data%min_cr_ratio + casize = prec%ag_data%min_coarse_size iszv = size(prec%precv) call psb_bcast(ictxt,iszv) call psb_bcast(ictxt,casize) call psb_bcast(ictxt,mxplevs) call psb_bcast(ictxt,mnaggratio) - if (casize /= prec%min_coarse_size) then + if (casize /= prec%ag_data%min_coarse_size) then info=psb_err_internal_error_ call psb_errpush(info,name,a_err='Inconsistent min_coarse_size') goto 9999 end if - if (mxplevs /= prec%max_levs) then + if (mxplevs /= prec%ag_data%max_levs) then info=psb_err_internal_error_ call psb_errpush(info,name,a_err='Inconsistent max_levs') goto 9999 end if - if (mnaggratio /= prec%min_cr_ratio) then + if (mnaggratio /= prec%ag_data%min_cr_ratio) then info=psb_err_internal_error_ call psb_errpush(info,name,a_err='Inconsistent min_cr_ratio') goto 9999 @@ -196,7 +196,7 @@ subroutine mld_s_hierarchy_bld(a,desc_a,prec,info) casize = max(casize,ione) casize = casize*40_psb_ipk_ call psb_bcast(ictxt,casize) - prec%min_coarse_size = casize + prec%ag_data%min_coarse_size = casize end if nplevs = max(itwo,mxplevs) @@ -298,7 +298,7 @@ subroutine mld_s_hierarchy_bld(a,desc_a,prec,info) if (info == psb_success_)& & call prec%precv(i)%bld_tprol(prec%precv(i-1)%base_a,& & prec%precv(i-1)%base_desc,& - & ilaggr,nlaggr,op_prol,info) + & ilaggr,nlaggr,op_prol,prec%ag_data,info) if (info /= psb_success_) then call psb_errpush(psb_err_internal_error_,name,& diff --git a/mlprec/impl/mld_scprecset.F90 b/mlprec/impl/mld_scprecset.F90 index 7af6c714..0032994b 100644 --- a/mlprec/impl/mld_scprecset.F90 +++ b/mlprec/impl/mld_scprecset.F90 @@ -129,7 +129,7 @@ subroutine mld_scprecseti(p,what,val,info,ilev,ilmax,pos,idx) end if else ilev_ = 1 - ilmax_ = ilev_ + ilmax_ = nlev_ end if if ((ilev_<1).or.(ilev_ > nlev_)) then info = -1 @@ -146,10 +146,10 @@ subroutine mld_scprecseti(p,what,val,info,ilev,ilmax,pos,idx) select case(psb_toupper(what)) case ('MIN_COARSE_SIZE') - p%min_coarse_size = max(val,-1) + p%ag_data%min_coarse_size = max(val,-1) return case('MAX_LEVS') - p%max_levs = max(val,1) + p%ag_data%max_levs = max(val,1) return case ('OUTER_SWEEPS') p%outer_sweeps = max(val,1) @@ -502,7 +502,7 @@ subroutine mld_scprecsetc(p,what,string,info,ilev,ilmax,pos,idx) end if else ilev_ = 1 - ilmax_ = ilev_ + ilmax_ = nlev_ end if if ((ilev_<1).or.(ilev_ > nlev_)) then info = -1 @@ -593,7 +593,7 @@ subroutine mld_scprecsetr(p,what,val,info,ilev,ilmax,pos,idx) select case(psb_toupper(what)) case ('MIN_CR_RATIO') - p%min_cr_ratio = max(sone,val) + p%ag_data%min_cr_ratio = max(sone,val) return end select @@ -615,7 +615,7 @@ subroutine mld_scprecsetr(p,what,val,info,ilev,ilmax,pos,idx) end if else ilev_ = 1 - ilmax_ = ilev_ + ilmax_ = nlev_ end if if ((ilev_<1).or.(ilev_ > nlev_)) then info = -1 diff --git a/mlprec/impl/mld_sprecinit.F90 b/mlprec/impl/mld_sprecinit.F90 index 6ed32381..c239d081 100644 --- a/mlprec/impl/mld_sprecinit.F90 +++ b/mlprec/impl/mld_sprecinit.F90 @@ -119,7 +119,7 @@ subroutine mld_sprecinit(ictxt,prec,ptype,info) endif endif prec%ictxt = ictxt - prec%min_coarse_size = -1 + prec%ag_data%min_coarse_size = -1 select case(psb_toupper(ptype(1:len_trim(ptype)))) case ('NOPREC','NONE') @@ -186,7 +186,7 @@ subroutine mld_sprecinit(ictxt,prec,ptype,info) case ('ML') - nlev_ = prec%max_levs + nlev_ = prec%ag_data%max_levs ilev_ = 1 allocate(prec%precv(nlev_),stat=info) diff --git a/mlprec/impl/mld_z_extprol_bld.F90 b/mlprec/impl/mld_z_extprol_bld.F90 index 028f3cbc..1fb7b727 100644 --- a/mlprec/impl/mld_z_extprol_bld.F90 +++ b/mlprec/impl/mld_z_extprol_bld.F90 @@ -160,9 +160,9 @@ subroutine mld_z_extprol_bld(a,desc_a,p,prolv,restrv,info,amold,vmold,imold) ! Check to ensure all procs have the same ! newsz = -1 - mxplevs = p%max_levs - mnaggratio = p%min_cr_ratio - casize = p%min_coarse_size + mxplevs = p%ag_data%max_levs + mnaggratio = p%ag_data%min_cr_ratio + casize = p%ag_data%min_coarse_size iszv = size(p%precv) nprolv = size(prolv) nrestrv = size(restrv) @@ -172,17 +172,17 @@ subroutine mld_z_extprol_bld(a,desc_a,p,prolv,restrv,info,amold,vmold,imold) call psb_bcast(ictxt,mnaggratio) call psb_bcast(ictxt,nprolv) call psb_bcast(ictxt,nrestrv) - if (casize /= p%min_coarse_size) then + if (casize /= p%ag_data%min_coarse_size) then info=psb_err_internal_error_ call psb_errpush(info,name,a_err='Inconsistent min_coarse_size') goto 9999 end if - if (mxplevs /= p%max_levs) then + if (mxplevs /= p%ag_data%max_levs) then info=psb_err_internal_error_ call psb_errpush(info,name,a_err='Inconsistent max_levs') goto 9999 end if - if (mnaggratio /= p%min_cr_ratio) then + if (mnaggratio /= p%ag_data%min_cr_ratio) then info=psb_err_internal_error_ call psb_errpush(info,name,a_err='Inconsistent min_cr_ratio') goto 9999 @@ -226,7 +226,7 @@ subroutine mld_z_extprol_bld(a,desc_a,p,prolv,restrv,info,amold,vmold,imold) ! nplevs = nrestrv + 1 - p%max_levs = nplevs + p%ag_data%max_levs = nplevs ! ! Fixed number of levels. diff --git a/mlprec/impl/mld_z_hierarchy_bld.f90 b/mlprec/impl/mld_z_hierarchy_bld.f90 index f33a76e7..c4dcac39 100644 --- a/mlprec/impl/mld_z_hierarchy_bld.f90 +++ b/mlprec/impl/mld_z_hierarchy_bld.f90 @@ -124,25 +124,25 @@ subroutine mld_z_hierarchy_bld(a,desc_a,prec,info) ! Check to ensure all procs have the same ! newsz = -1 - mxplevs = prec%max_levs - mnaggratio = prec%min_cr_ratio - casize = prec%min_coarse_size + mxplevs = prec%ag_data%max_levs + mnaggratio = prec%ag_data%min_cr_ratio + casize = prec%ag_data%min_coarse_size iszv = size(prec%precv) call psb_bcast(ictxt,iszv) call psb_bcast(ictxt,casize) call psb_bcast(ictxt,mxplevs) call psb_bcast(ictxt,mnaggratio) - if (casize /= prec%min_coarse_size) then + if (casize /= prec%ag_data%min_coarse_size) then info=psb_err_internal_error_ call psb_errpush(info,name,a_err='Inconsistent min_coarse_size') goto 9999 end if - if (mxplevs /= prec%max_levs) then + if (mxplevs /= prec%ag_data%max_levs) then info=psb_err_internal_error_ call psb_errpush(info,name,a_err='Inconsistent max_levs') goto 9999 end if - if (mnaggratio /= prec%min_cr_ratio) then + if (mnaggratio /= prec%ag_data%min_cr_ratio) then info=psb_err_internal_error_ call psb_errpush(info,name,a_err='Inconsistent min_cr_ratio') goto 9999 @@ -196,7 +196,7 @@ subroutine mld_z_hierarchy_bld(a,desc_a,prec,info) casize = max(casize,ione) casize = casize*40_psb_ipk_ call psb_bcast(ictxt,casize) - prec%min_coarse_size = casize + prec%ag_data%min_coarse_size = casize end if nplevs = max(itwo,mxplevs) @@ -298,7 +298,7 @@ subroutine mld_z_hierarchy_bld(a,desc_a,prec,info) if (info == psb_success_)& & call prec%precv(i)%bld_tprol(prec%precv(i-1)%base_a,& & prec%precv(i-1)%base_desc,& - & ilaggr,nlaggr,op_prol,info) + & ilaggr,nlaggr,op_prol,prec%ag_data,info) if (info /= psb_success_) then call psb_errpush(psb_err_internal_error_,name,& diff --git a/mlprec/impl/mld_zcprecset.F90 b/mlprec/impl/mld_zcprecset.F90 index 37b41d67..72de1558 100644 --- a/mlprec/impl/mld_zcprecset.F90 +++ b/mlprec/impl/mld_zcprecset.F90 @@ -135,7 +135,7 @@ subroutine mld_zcprecseti(p,what,val,info,ilev,ilmax,pos,idx) end if else ilev_ = 1 - ilmax_ = ilev_ + ilmax_ = nlev_ end if if ((ilev_<1).or.(ilev_ > nlev_)) then info = -1 @@ -152,10 +152,10 @@ subroutine mld_zcprecseti(p,what,val,info,ilev,ilmax,pos,idx) select case(psb_toupper(what)) case ('MIN_COARSE_SIZE') - p%min_coarse_size = max(val,-1) + p%ag_data%min_coarse_size = max(val,-1) return case('MAX_LEVS') - p%max_levs = max(val,1) + p%ag_data%max_levs = max(val,1) return case ('OUTER_SWEEPS') p%outer_sweeps = max(val,1) @@ -536,7 +536,7 @@ subroutine mld_zcprecsetc(p,what,string,info,ilev,ilmax,pos,idx) end if else ilev_ = 1 - ilmax_ = ilev_ + ilmax_ = nlev_ end if if ((ilev_<1).or.(ilev_ > nlev_)) then info = -1 @@ -627,7 +627,7 @@ subroutine mld_zcprecsetr(p,what,val,info,ilev,ilmax,pos,idx) select case(psb_toupper(what)) case ('MIN_CR_RATIO') - p%min_cr_ratio = max(done,val) + p%ag_data%min_cr_ratio = max(done,val) return end select @@ -649,7 +649,7 @@ subroutine mld_zcprecsetr(p,what,val,info,ilev,ilmax,pos,idx) end if else ilev_ = 1 - ilmax_ = ilev_ + ilmax_ = nlev_ end if if ((ilev_<1).or.(ilev_ > nlev_)) then info = -1 diff --git a/mlprec/impl/mld_zprecinit.F90 b/mlprec/impl/mld_zprecinit.F90 index 197b04b6..27299b57 100644 --- a/mlprec/impl/mld_zprecinit.F90 +++ b/mlprec/impl/mld_zprecinit.F90 @@ -122,7 +122,7 @@ subroutine mld_zprecinit(ictxt,prec,ptype,info) endif endif prec%ictxt = ictxt - prec%min_coarse_size = -1 + prec%ag_data%min_coarse_size = -1 select case(psb_toupper(ptype(1:len_trim(ptype)))) case ('NOPREC','NONE') @@ -189,7 +189,7 @@ subroutine mld_zprecinit(ictxt,prec,ptype,info) case ('ML') - nlev_ = prec%max_levs + nlev_ = prec%ag_data%max_levs ilev_ = 1 allocate(prec%precv(nlev_),stat=info) diff --git a/mlprec/mld_base_prec_type.F90 b/mlprec/mld_base_prec_type.F90 index 11179b6f..c6d3957e 100644 --- a/mlprec/mld_base_prec_type.F90 +++ b/mlprec/mld_base_prec_type.F90 @@ -120,6 +120,44 @@ module mld_base_prec_type procedure, pass(pm) :: printout => d_ml_parms_printout end type mld_dml_parms + type mld_saggr_data + ! + ! Aggregation data and defaults: + ! + ! 1. min_coarse_size = 0 Default target size will be computed as + ! 40*(N_fine)**(1./3.) + ! We are assuming that the coarse size fits in + ! integer range of psb_ipk_, but this is + ! not very restrictive + integer(psb_ipk_) :: min_coarse_size = izero + ! 2. maximum number of levels. Defaults to 20 + integer(psb_ipk_) :: max_levs = 20_psb_ipk_ + ! 3. min_cr_ratio = 1.5 + real(psb_spk_) :: min_cr_ratio = 1.5_psb_spk_ + real(psb_spk_) :: op_complexity = szero + real(psb_spk_) :: avg_cr = szero + end type mld_saggr_data + + type mld_daggr_data + ! + ! Aggregation data and defaults: + ! + ! + ! 1. min_coarse_size = 0 Default target size will be computed as + ! 40*(N_fine)**(1./3.) + ! We are assuming that the coarse size fits in + ! integer range of psb_ipk_, but this is + ! not very restrictive + integer(psb_ipk_) :: min_coarse_size = izero + ! 2. maximum number of levels. Defaults to 20 + integer(psb_ipk_) :: max_levs = 20_psb_ipk_ + ! 3. min_cr_ratio = 1.5 + real(psb_dpk_) :: min_cr_ratio = 1.5_psb_dpk_ + real(psb_dpk_) :: op_complexity = dzero + real(psb_dpk_) :: avg_cr = dzero + end type mld_daggr_data + + ! ! Entries in iprcparm diff --git a/mlprec/mld_c_base_aggregator_mod.f90 b/mlprec/mld_c_base_aggregator_mod.f90 index b089cccb..61b8e80a 100644 --- a/mlprec/mld_c_base_aggregator_mod.f90 +++ b/mlprec/mld_c_base_aggregator_mod.f90 @@ -41,11 +41,11 @@ ! module mld_c_base_aggregator_mod - use mld_base_prec_type, only : mld_sml_parms + use mld_base_prec_type, only : mld_sml_parms, mld_saggr_data use psb_base_mod, only : psb_cspmat_type, psb_lcspmat_type, psb_c_vect_type, & & psb_c_base_vect_type, psb_clinmap_type, psb_spk_, & & psb_ipk_, psb_epk_, psb_lpk_, psb_desc_type, psb_i_base_vect_type, & - & psb_erractionsave, psb_error_handler, psb_success_ + & psb_erractionsave, psb_error_handler, psb_success_, psb_toupper ! ! ! @@ -79,7 +79,8 @@ module mld_c_base_aggregator_mod !! cseti, csetr, csetc - Set internal parameters, if any ! type mld_c_base_aggregator_type - + ! Do we want to purge explicit zeros when aggregating? + logical :: do_clean_zeros contains procedure, pass(ag) :: bld_tprol => mld_c_base_aggregator_build_tprol procedure, pass(ag) :: mat_asb => mld_c_base_aggregator_mat_asb @@ -96,6 +97,19 @@ module mld_c_base_aggregator_mod generic, public :: set => cseti, csetr, csetc end type mld_c_base_aggregator_type + abstract interface + subroutine mld_c_soc_map_bld(iorder,theta,clean_zeros,a,desc_a,nlaggr,ilaggr,info) + import :: psb_cspmat_type, psb_desc_type, psb_spk_, psb_ipk_, psb_lpk_ + implicit none + integer(psb_ipk_), intent(in) :: iorder + logical, intent(in) :: clean_zeros + type(psb_cspmat_type), intent(in) :: a + type(psb_desc_type), intent(in) :: desc_a + real(psb_spk_), intent(in) :: theta + integer(psb_lpk_), allocatable, intent(out) :: ilaggr(:),nlaggr(:) + integer(psb_ipk_), intent(out) :: info + end subroutine mld_c_soc_map_bld + end interface contains @@ -137,7 +151,16 @@ contains character(len=*), intent(in) :: val integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(in), optional :: idx - ! Do nothing + ! Set clean zeros, or do nothing. + select case (psb_toupper(trim(what))) + case('AGGR_CLEAN_ZEROS') + select case (psb_toupper(trim(val))) + case('TRUE','T') + ag%do_clean_zeros = .true. + case('FALSE','F') + ag%do_clean_zeros = .false. + end select + end select info = 0 end subroutine mld_c_base_aggregator_csetc @@ -181,8 +204,8 @@ contains subroutine mld_c_base_aggregator_default(ag) implicit none class(mld_c_base_aggregator_type), intent(inout) :: ag - - ! Here we need do nothing + ! Only one default setting + ag%do_clean_zeros = .true. return end subroutine mld_c_base_aggregator_default @@ -228,9 +251,12 @@ contains !! will contribute, in global numbering. !! Many aggregation produce a binary tentative prolongator, but some !! do not, hence we also need the OP_PROL output. + !! AG_DATA is passed here just in case some of the + !! aggregators need it internally, most of them will ignore. !! !! \param ag The input aggregator object !! \param parms The auxiliary parameters object + !! \param ag_data Auxiliary global aggregation info !! \param a The local matrix part !! \param desc_a The descriptor !! \param ilaggr Output aggregation map @@ -239,11 +265,13 @@ contains !! \param info Return code !! ! - subroutine mld_c_base_aggregator_build_tprol(ag,parms,a,desc_a,ilaggr,nlaggr,op_prol,info) + subroutine mld_c_base_aggregator_build_tprol(ag,parms,ag_data,& + & a,desc_a,ilaggr,nlaggr,op_prol,info) use psb_base_mod implicit none class(mld_c_base_aggregator_type), target, intent(inout) :: ag - type(mld_sml_parms), intent(inout) :: parms + type(mld_sml_parms), intent(inout) :: parms + type(mld_saggr_data), intent(in) :: ag_data type(psb_cspmat_type), intent(in) :: a type(psb_desc_type), intent(in) :: desc_a integer(psb_lpk_), allocatable, intent(out) :: ilaggr(:),nlaggr(:) diff --git a/mlprec/mld_c_dec_aggregator_mod.f90 b/mlprec/mld_c_dec_aggregator_mod.f90 index 3d21c09c..40070838 100644 --- a/mlprec/mld_c_dec_aggregator_mod.f90 +++ b/mlprec/mld_c_dec_aggregator_mod.f90 @@ -48,7 +48,7 @@ module mld_c_dec_aggregator_mod !! \extends mld_c_base_aggregator_mod::mld_c_base_aggregator_type !! !! type, extends(mld_c_base_aggregator_type) :: mld_c_dec_aggregator_type - !! procedure(mld_c_map_bld), nopass, pointer :: map_bld => null() + !! procedure(mld_c_soc_map_bld), nopass, pointer :: soc_map_bld => null() !! end type !! !! This is the simplest aggregation method: starting from the @@ -71,12 +71,12 @@ module mld_c_dec_aggregator_mod !! PSBLAS-based parallel two-level Schwarz preconditioners, Appl. Num. Math. !! 57 (2007), 1181-1196. !! - !! The map_bld method is used inside the implementation of build_tprol + !! The soc_map_bld method is used inside the implementation of build_tprol !! ! ! type, extends(mld_c_base_aggregator_type) :: mld_c_dec_aggregator_type - procedure(mld_c_map_bld), nopass, pointer :: map_bld => null() + procedure(mld_c_soc_map_bld), nopass, pointer :: soc_map_bld => null() contains procedure, pass(ag) :: bld_tprol => mld_c_dec_aggregator_build_tprol @@ -88,28 +88,17 @@ module mld_c_dec_aggregator_mod end type mld_c_dec_aggregator_type - abstract interface - subroutine mld_c_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info) - import :: psb_cspmat_type, psb_desc_type, psb_spk_, psb_ipk_, psb_lpk_ - implicit none - integer(psb_ipk_), intent(in) :: iorder - type(psb_cspmat_type), intent(in) :: a - type(psb_desc_type), intent(in) :: desc_a - real(psb_spk_), intent(in) :: theta - integer(psb_lpk_), allocatable, intent(out) :: ilaggr(:),nlaggr(:) - integer(psb_ipk_), intent(out) :: info - end subroutine mld_c_map_bld - end interface - - procedure(mld_c_map_bld) :: mld_c_soc1_map_bld, mld_c_soc2_map_bld + procedure(mld_c_soc_map_bld) :: mld_c_soc1_map_bld, mld_c_soc2_map_bld interface - subroutine mld_c_dec_aggregator_build_tprol(ag,parms,a,desc_a,ilaggr,nlaggr,op_prol,info) + subroutine mld_c_dec_aggregator_build_tprol(ag,parms,ag_data,& + & a,desc_a,ilaggr,nlaggr,op_prol,info) import :: mld_c_dec_aggregator_type, psb_desc_type, psb_cspmat_type, psb_spk_, & - & psb_ipk_, psb_lpk_, psb_lcspmat_type, mld_sml_parms + & psb_ipk_, psb_lpk_, psb_lcspmat_type, mld_sml_parms, mld_saggr_data implicit none class(mld_c_dec_aggregator_type), target, intent(inout) :: ag type(mld_sml_parms), intent(inout) :: parms + type(mld_saggr_data), intent(in) :: ag_data type(psb_cspmat_type), intent(in) :: a type(psb_desc_type), intent(in) :: desc_a integer(psb_lpk_), allocatable, intent(out) :: ilaggr(:),nlaggr(:) @@ -147,14 +136,14 @@ contains select case(parms%aggr_type) case (mld_noalg_) - ag%map_bld => null() + ag%soc_map_bld => null() case (mld_soc1_) - ag%map_bld => mld_c_soc1_map_bld + ag%soc_map_bld => mld_c_soc1_map_bld case (mld_soc2_) - ag%map_bld => mld_c_soc2_map_bld + ag%soc_map_bld => mld_c_soc2_map_bld case default write(0,*) 'Unknown aggregation type, defaulting to SOC1' - ag%map_bld => mld_c_soc1_map_bld + ag%soc_map_bld => mld_c_soc1_map_bld end select return @@ -165,7 +154,8 @@ contains implicit none class(mld_c_dec_aggregator_type), intent(inout) :: ag - ag%map_bld => mld_c_soc1_map_bld + call ag%mld_c_base_aggregator_type%default() + ag%soc_map_bld => mld_c_soc1_map_bld return end subroutine mld_c_dec_aggregator_default diff --git a/mlprec/mld_c_inner_mod.f90 b/mlprec/mld_c_inner_mod.f90 index 2ac6caf8..9a141310 100644 --- a/mlprec/mld_c_inner_mod.f90 +++ b/mlprec/mld_c_inner_mod.f90 @@ -97,32 +97,6 @@ module mld_c_inner_mod end subroutine mld_cmlprec_aply_vect end interface mld_mlprec_aply - interface mld_aggrmap_bld - subroutine mld_c_lev_aggrmap_bld(p,a,desc_a,ilaggr,nlaggr,op_prol,info) - import :: psb_cspmat_type, psb_desc_type, psb_spk_, psb_ipk_, psb_lpk_, psb_lcspmat_type - import :: mld_c_onelev_type - implicit none - type(mld_c_onelev_type), intent(inout), target :: p - type(psb_cspmat_type), intent(in) :: a - type(psb_desc_type), intent(in) :: desc_a - integer(psb_lpk_), intent(out) :: ilaggr(:),nlaggr(:) - type(psb_lcspmat_type), intent(out) :: op_prol - integer(psb_ipk_), intent(out) :: info - end subroutine mld_c_lev_aggrmap_bld - subroutine mld_caggrmap_bld(aggr_type,iorder,theta,a,desc_a,ilaggr,nlaggr,op_prol,info) - import :: psb_cspmat_type, psb_desc_type, psb_spk_, psb_ipk_, psb_lpk_, psb_lcspmat_type - implicit none - integer(psb_ipk_), intent(in) :: iorder - integer(psb_ipk_), intent(in) :: aggr_type - real(psb_spk_), intent(in) :: theta - type(psb_cspmat_type), intent(in) :: a - type(psb_desc_type), intent(in) :: desc_a - integer(psb_lpk_), allocatable, intent(out) :: ilaggr(:),nlaggr(:) - type(psb_lcspmat_type), intent(out) :: op_prol - integer(psb_ipk_), intent(out) :: info - end subroutine mld_caggrmap_bld - end interface mld_aggrmap_bld - interface mld_map_to_tprol subroutine mld_c_map_to_tprol(desc_a,ilaggr,nlaggr,op_prol,info) use psb_base_mod, only : psb_cspmat_type, psb_desc_type, psb_spk_, psb_ipk_, psb_lpk_, psb_lcspmat_type diff --git a/mlprec/mld_c_onelev_mod.f90 b/mlprec/mld_c_onelev_mod.f90 index 475e91e1..2ab4253b 100644 --- a/mlprec/mld_c_onelev_mod.f90 +++ b/mlprec/mld_c_onelev_mod.f90 @@ -381,7 +381,7 @@ interface interface subroutine mld_c_base_onelev_dump(lv,level,info,prefix,head,ac,rp,smoother,& - & solver,global_num) + & solver,tprol,global_num) import :: psb_cspmat_type, psb_c_vect_type, psb_c_base_vect_type, & & psb_clinmap_type, psb_spk_, mld_c_onelev_type, & & psb_ipk_, psb_epk_, psb_desc_type @@ -390,7 +390,7 @@ interface integer(psb_ipk_), intent(in) :: level integer(psb_ipk_), intent(out) :: info character(len=*), intent(in), optional :: prefix, head - logical, optional, intent(in) :: ac, rp, smoother, solver, global_num + logical, optional, intent(in) :: ac, rp, smoother, solver, tprol, global_num end subroutine mld_c_base_onelev_dump end interface @@ -484,16 +484,18 @@ contains end subroutine c_base_onelev_default - subroutine c_base_onelev_bld_tprol(lv,a,desc_a,ilaggr,nlaggr,op_prol,info) + subroutine c_base_onelev_bld_tprol(lv,a,desc_a,& + & ilaggr,nlaggr,op_prol,ag_data,info) implicit none class(mld_c_onelev_type), intent(inout), target :: lv type(psb_cspmat_type), intent(in) :: a type(psb_desc_type), intent(in) :: desc_a integer(psb_lpk_), allocatable, intent(out) :: ilaggr(:),nlaggr(:) type(psb_lcspmat_type), intent(out) :: op_prol + type(mld_saggr_data), intent(in) :: ag_data integer(psb_ipk_), intent(out) :: info - call lv%aggr%bld_tprol(lv%parms,a,desc_a,ilaggr,nlaggr,op_prol,info) + call lv%aggr%bld_tprol(lv%parms,ag_data,a,desc_a,ilaggr,nlaggr,op_prol,info) end subroutine c_base_onelev_bld_tprol diff --git a/mlprec/mld_c_prec_type.f90 b/mlprec/mld_c_prec_type.f90 index 686defd5..96f660dd 100644 --- a/mlprec/mld_c_prec_type.f90 +++ b/mlprec/mld_c_prec_type.f90 @@ -85,18 +85,8 @@ module mld_c_prec_type integer, parameter, private :: wv_size_=4 type, extends(psb_cprec_type) :: mld_cprec_type - ! integer(psb_ipk_) :: ictxt ! Now it's in the PSBLAS prec. - ! - ! Aggregation defaults: - ! - ! 1. min_coarse_size = 0 Default target size will be computed as 40*(N_fine)**(1./3.) - integer(psb_ipk_) :: min_coarse_size = izero - ! 2. maximum number of levels. Defaults to 20 - integer(psb_ipk_) :: max_levs = 20_psb_ipk_ - ! 3. min_cr_ratio = 1.5 - real(psb_spk_) :: min_cr_ratio = 1.5_psb_spk_ - real(psb_spk_) :: op_complexity = szero - real(psb_spk_) :: avg_cr = szero + ! integer(psb_ipk_) :: ictxt ! Now it's in the PSBLAS prec. + type(mld_saggr_data) :: ag_data ! ! Number of outer sweeps. Sometimes 2 V-cycles may be better than 1 W-cycle. ! @@ -445,7 +435,7 @@ contains class(mld_cprec_type), intent(in) :: prec complex(psb_spk_) :: val - val = prec%op_complexity + val = prec%ag_data%op_complexity end function mld_c_get_compl @@ -480,7 +470,7 @@ contains call psb_sum(ictxt,num) call psb_sum(ictxt,den) end if - prec%op_complexity = num/den + prec%ag_data%op_complexity = num/den end subroutine mld_c_cmp_compl ! @@ -492,7 +482,7 @@ contains class(mld_cprec_type), intent(in) :: prec complex(psb_spk_) :: val - val = prec%avg_cr + val = prec%ag_data%avg_cr end function mld_c_get_avg_cr @@ -517,7 +507,7 @@ contains avgcr = avgcr / (nl-1) end if call psb_sum(ictxt,avgcr) - prec%avg_cr = avgcr/np + prec%ag_data%avg_cr = avgcr/np end subroutine mld_c_cmp_avg_cr ! @@ -722,14 +712,15 @@ contains end subroutine mld_c_apply1v - subroutine mld_c_dump(prec,info,istart,iend,prefix,head,ac,rp,smoother,solver,global_num) + subroutine mld_c_dump(prec,info,istart,iend,prefix,head,ac,rp,smoother,solver,tprol,& + & global_num) implicit none class(mld_cprec_type), intent(in) :: prec integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(in), optional :: istart, iend character(len=*), intent(in), optional :: prefix, head - logical, optional, intent(in) :: smoother, solver,ac, rp, global_num + logical, optional, intent(in) :: smoother, solver,ac, rp, tprol, global_num integer(psb_ipk_) :: i, j, il1, iln, lname, lev integer(psb_ipk_) :: icontxt,iam, np character(len=80) :: prefix_ @@ -750,7 +741,8 @@ contains do lev=il1, iln call prec%precv(lev)%dump(lev,info,prefix=prefix,head=head,& - & ac=ac,smoother=smoother,solver=solver,rp=rp,global_num=global_num) + & ac=ac,smoother=smoother,solver=solver,rp=rp,tprol=tprol, & + & global_num=global_num) end do end subroutine mld_c_dump @@ -801,13 +793,9 @@ contains info = psb_success_ select type(pout => precout) class is (mld_cprec_type) - pout%ictxt = prec%ictxt - pout%max_levs = prec%max_levs - pout%min_coarse_size = prec%min_coarse_size - pout%min_cr_ratio = prec%min_cr_ratio - pout%outer_sweeps = prec%outer_sweeps - pout%op_complexity = prec%op_complexity - pout%avg_cr = prec%avg_cr + pout%ictxt = prec%ictxt + pout%ag_data = prec%ag_data + pout%outer_sweeps = prec%outer_sweeps if (allocated(prec%precv)) then ln = size(prec%precv) allocate(pout%precv(ln),stat=info) @@ -853,7 +841,10 @@ contains !!$ return endif end if - + b%ictxt = prec%ictxt + b%ag_data = prec%ag_data + b%outer_sweeps = prec%outer_sweeps + call move_alloc(prec%precv,b%precv) ! Fix the pointers except on level 1. do i=2, size(b%precv) diff --git a/mlprec/mld_c_symdec_aggregator_mod.f90 b/mlprec/mld_c_symdec_aggregator_mod.f90 index f9436ddf..b96dfa8d 100644 --- a/mlprec/mld_c_symdec_aggregator_mod.f90 +++ b/mlprec/mld_c_symdec_aggregator_mod.f90 @@ -62,12 +62,14 @@ module mld_c_symdec_aggregator_mod interface - subroutine mld_c_symdec_aggregator_build_tprol(ag,parms,a,desc_a,ilaggr,nlaggr,op_prol,info) + subroutine mld_c_symdec_aggregator_build_tprol(ag,parms,ag_data,& + & a,desc_a,ilaggr,nlaggr,op_prol,info) import :: mld_c_symdec_aggregator_type, psb_desc_type, psb_cspmat_type, psb_spk_, & - & psb_ipk_, psb_lpk_, psb_lcspmat_type, mld_sml_parms + & psb_ipk_, psb_lpk_, psb_lcspmat_type, mld_sml_parms, mld_saggr_data implicit none class(mld_c_symdec_aggregator_type), target, intent(inout) :: ag type(mld_sml_parms), intent(inout) :: parms + type(mld_saggr_data), intent(in) :: ag_data type(psb_cspmat_type), intent(in) :: a type(psb_desc_type), intent(in) :: desc_a integer(psb_lpk_), allocatable, intent(out) :: ilaggr(:),nlaggr(:) diff --git a/mlprec/mld_d_base_aggregator_mod.f90 b/mlprec/mld_d_base_aggregator_mod.f90 index 5e154f80..fe92d7be 100644 --- a/mlprec/mld_d_base_aggregator_mod.f90 +++ b/mlprec/mld_d_base_aggregator_mod.f90 @@ -41,11 +41,11 @@ ! module mld_d_base_aggregator_mod - use mld_base_prec_type, only : mld_dml_parms + use mld_base_prec_type, only : mld_dml_parms, mld_daggr_data use psb_base_mod, only : psb_dspmat_type, psb_ldspmat_type, psb_d_vect_type, & & psb_d_base_vect_type, psb_dlinmap_type, psb_dpk_, & & psb_ipk_, psb_epk_, psb_lpk_, psb_desc_type, psb_i_base_vect_type, & - & psb_erractionsave, psb_error_handler, psb_success_ + & psb_erractionsave, psb_error_handler, psb_success_, psb_toupper ! ! ! @@ -79,7 +79,8 @@ module mld_d_base_aggregator_mod !! cseti, csetr, csetc - Set internal parameters, if any ! type mld_d_base_aggregator_type - + ! Do we want to purge explicit zeros when aggregating? + logical :: do_clean_zeros contains procedure, pass(ag) :: bld_tprol => mld_d_base_aggregator_build_tprol procedure, pass(ag) :: mat_asb => mld_d_base_aggregator_mat_asb @@ -96,6 +97,19 @@ module mld_d_base_aggregator_mod generic, public :: set => cseti, csetr, csetc end type mld_d_base_aggregator_type + abstract interface + subroutine mld_d_soc_map_bld(iorder,theta,clean_zeros,a,desc_a,nlaggr,ilaggr,info) + import :: psb_dspmat_type, psb_desc_type, psb_dpk_, psb_ipk_, psb_lpk_ + implicit none + integer(psb_ipk_), intent(in) :: iorder + logical, intent(in) :: clean_zeros + type(psb_dspmat_type), intent(in) :: a + type(psb_desc_type), intent(in) :: desc_a + real(psb_dpk_), intent(in) :: theta + integer(psb_lpk_), allocatable, intent(out) :: ilaggr(:),nlaggr(:) + integer(psb_ipk_), intent(out) :: info + end subroutine mld_d_soc_map_bld + end interface contains @@ -137,7 +151,16 @@ contains character(len=*), intent(in) :: val integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(in), optional :: idx - ! Do nothing + ! Set clean zeros, or do nothing. + select case (psb_toupper(trim(what))) + case('AGGR_CLEAN_ZEROS') + select case (psb_toupper(trim(val))) + case('TRUE','T') + ag%do_clean_zeros = .true. + case('FALSE','F') + ag%do_clean_zeros = .false. + end select + end select info = 0 end subroutine mld_d_base_aggregator_csetc @@ -181,8 +204,8 @@ contains subroutine mld_d_base_aggregator_default(ag) implicit none class(mld_d_base_aggregator_type), intent(inout) :: ag - - ! Here we need do nothing + ! Only one default setting + ag%do_clean_zeros = .true. return end subroutine mld_d_base_aggregator_default @@ -228,9 +251,12 @@ contains !! will contribute, in global numbering. !! Many aggregation produce a binary tentative prolongator, but some !! do not, hence we also need the OP_PROL output. + !! AG_DATA is passed here just in case some of the + !! aggregators need it internally, most of them will ignore. !! !! \param ag The input aggregator object !! \param parms The auxiliary parameters object + !! \param ag_data Auxiliary global aggregation info !! \param a The local matrix part !! \param desc_a The descriptor !! \param ilaggr Output aggregation map @@ -239,11 +265,13 @@ contains !! \param info Return code !! ! - subroutine mld_d_base_aggregator_build_tprol(ag,parms,a,desc_a,ilaggr,nlaggr,op_prol,info) + subroutine mld_d_base_aggregator_build_tprol(ag,parms,ag_data,& + & a,desc_a,ilaggr,nlaggr,op_prol,info) use psb_base_mod implicit none class(mld_d_base_aggregator_type), target, intent(inout) :: ag - type(mld_dml_parms), intent(inout) :: parms + type(mld_dml_parms), intent(inout) :: parms + type(mld_daggr_data), intent(in) :: ag_data type(psb_dspmat_type), intent(in) :: a type(psb_desc_type), intent(in) :: desc_a integer(psb_lpk_), allocatable, intent(out) :: ilaggr(:),nlaggr(:) diff --git a/mlprec/mld_d_dec_aggregator_mod.f90 b/mlprec/mld_d_dec_aggregator_mod.f90 index 5d9b63b7..71dfaf0a 100644 --- a/mlprec/mld_d_dec_aggregator_mod.f90 +++ b/mlprec/mld_d_dec_aggregator_mod.f90 @@ -48,7 +48,7 @@ module mld_d_dec_aggregator_mod !! \extends mld_d_base_aggregator_mod::mld_d_base_aggregator_type !! !! type, extends(mld_d_base_aggregator_type) :: mld_d_dec_aggregator_type - !! procedure(mld_d_map_bld), nopass, pointer :: map_bld => null() + !! procedure(mld_d_soc_map_bld), nopass, pointer :: soc_map_bld => null() !! end type !! !! This is the simplest aggregation method: starting from the @@ -71,12 +71,12 @@ module mld_d_dec_aggregator_mod !! PSBLAS-based parallel two-level Schwarz preconditioners, Appl. Num. Math. !! 57 (2007), 1181-1196. !! - !! The map_bld method is used inside the implementation of build_tprol + !! The soc_map_bld method is used inside the implementation of build_tprol !! ! ! type, extends(mld_d_base_aggregator_type) :: mld_d_dec_aggregator_type - procedure(mld_d_map_bld), nopass, pointer :: map_bld => null() + procedure(mld_d_soc_map_bld), nopass, pointer :: soc_map_bld => null() contains procedure, pass(ag) :: bld_tprol => mld_d_dec_aggregator_build_tprol @@ -88,28 +88,17 @@ module mld_d_dec_aggregator_mod end type mld_d_dec_aggregator_type - abstract interface - subroutine mld_d_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info) - import :: psb_dspmat_type, psb_desc_type, psb_dpk_, psb_ipk_, psb_lpk_ - implicit none - integer(psb_ipk_), intent(in) :: iorder - type(psb_dspmat_type), intent(in) :: a - type(psb_desc_type), intent(in) :: desc_a - real(psb_dpk_), intent(in) :: theta - integer(psb_lpk_), allocatable, intent(out) :: ilaggr(:),nlaggr(:) - integer(psb_ipk_), intent(out) :: info - end subroutine mld_d_map_bld - end interface - - procedure(mld_d_map_bld) :: mld_d_soc1_map_bld, mld_d_soc2_map_bld + procedure(mld_d_soc_map_bld) :: mld_d_soc1_map_bld, mld_d_soc2_map_bld interface - subroutine mld_d_dec_aggregator_build_tprol(ag,parms,a,desc_a,ilaggr,nlaggr,op_prol,info) + subroutine mld_d_dec_aggregator_build_tprol(ag,parms,ag_data,& + & a,desc_a,ilaggr,nlaggr,op_prol,info) import :: mld_d_dec_aggregator_type, psb_desc_type, psb_dspmat_type, psb_dpk_, & - & psb_ipk_, psb_lpk_, psb_ldspmat_type, mld_dml_parms + & psb_ipk_, psb_lpk_, psb_ldspmat_type, mld_dml_parms, mld_daggr_data implicit none class(mld_d_dec_aggregator_type), target, intent(inout) :: ag type(mld_dml_parms), intent(inout) :: parms + type(mld_daggr_data), intent(in) :: ag_data type(psb_dspmat_type), intent(in) :: a type(psb_desc_type), intent(in) :: desc_a integer(psb_lpk_), allocatable, intent(out) :: ilaggr(:),nlaggr(:) @@ -147,14 +136,14 @@ contains select case(parms%aggr_type) case (mld_noalg_) - ag%map_bld => null() + ag%soc_map_bld => null() case (mld_soc1_) - ag%map_bld => mld_d_soc1_map_bld + ag%soc_map_bld => mld_d_soc1_map_bld case (mld_soc2_) - ag%map_bld => mld_d_soc2_map_bld + ag%soc_map_bld => mld_d_soc2_map_bld case default write(0,*) 'Unknown aggregation type, defaulting to SOC1' - ag%map_bld => mld_d_soc1_map_bld + ag%soc_map_bld => mld_d_soc1_map_bld end select return @@ -165,7 +154,8 @@ contains implicit none class(mld_d_dec_aggregator_type), intent(inout) :: ag - ag%map_bld => mld_d_soc1_map_bld + call ag%mld_d_base_aggregator_type%default() + ag%soc_map_bld => mld_d_soc1_map_bld return end subroutine mld_d_dec_aggregator_default diff --git a/mlprec/mld_d_inner_mod.f90 b/mlprec/mld_d_inner_mod.f90 index 4d9b4138..83977dd9 100644 --- a/mlprec/mld_d_inner_mod.f90 +++ b/mlprec/mld_d_inner_mod.f90 @@ -97,32 +97,6 @@ module mld_d_inner_mod end subroutine mld_dmlprec_aply_vect end interface mld_mlprec_aply - interface mld_aggrmap_bld - subroutine mld_d_lev_aggrmap_bld(p,a,desc_a,ilaggr,nlaggr,op_prol,info) - import :: psb_dspmat_type, psb_desc_type, psb_dpk_, psb_ipk_, psb_lpk_, psb_ldspmat_type - import :: mld_d_onelev_type - implicit none - type(mld_d_onelev_type), intent(inout), target :: p - type(psb_dspmat_type), intent(in) :: a - type(psb_desc_type), intent(in) :: desc_a - integer(psb_lpk_), intent(out) :: ilaggr(:),nlaggr(:) - type(psb_ldspmat_type), intent(out) :: op_prol - integer(psb_ipk_), intent(out) :: info - end subroutine mld_d_lev_aggrmap_bld - subroutine mld_daggrmap_bld(aggr_type,iorder,theta,a,desc_a,ilaggr,nlaggr,op_prol,info) - import :: psb_dspmat_type, psb_desc_type, psb_dpk_, psb_ipk_, psb_lpk_, psb_ldspmat_type - implicit none - integer(psb_ipk_), intent(in) :: iorder - integer(psb_ipk_), intent(in) :: aggr_type - real(psb_dpk_), intent(in) :: theta - type(psb_dspmat_type), intent(in) :: a - type(psb_desc_type), intent(in) :: desc_a - integer(psb_lpk_), allocatable, intent(out) :: ilaggr(:),nlaggr(:) - type(psb_ldspmat_type), intent(out) :: op_prol - integer(psb_ipk_), intent(out) :: info - end subroutine mld_daggrmap_bld - end interface mld_aggrmap_bld - interface mld_map_to_tprol subroutine mld_d_map_to_tprol(desc_a,ilaggr,nlaggr,op_prol,info) use psb_base_mod, only : psb_dspmat_type, psb_desc_type, psb_dpk_, psb_ipk_, psb_lpk_, psb_ldspmat_type diff --git a/mlprec/mld_d_onelev_mod.f90 b/mlprec/mld_d_onelev_mod.f90 index 34f3be63..5c8c1c82 100644 --- a/mlprec/mld_d_onelev_mod.f90 +++ b/mlprec/mld_d_onelev_mod.f90 @@ -381,7 +381,7 @@ interface interface subroutine mld_d_base_onelev_dump(lv,level,info,prefix,head,ac,rp,smoother,& - & solver,global_num) + & solver,tprol,global_num) import :: psb_dspmat_type, psb_d_vect_type, psb_d_base_vect_type, & & psb_dlinmap_type, psb_dpk_, mld_d_onelev_type, & & psb_ipk_, psb_epk_, psb_desc_type @@ -390,7 +390,7 @@ interface integer(psb_ipk_), intent(in) :: level integer(psb_ipk_), intent(out) :: info character(len=*), intent(in), optional :: prefix, head - logical, optional, intent(in) :: ac, rp, smoother, solver, global_num + logical, optional, intent(in) :: ac, rp, smoother, solver, tprol, global_num end subroutine mld_d_base_onelev_dump end interface @@ -422,7 +422,7 @@ contains val = val + lv%desc_ac%sizeof() val = val + lv%ac%sizeof() val = val + lv%tprol%sizeof() - val = val + lv%map%sizeof() + val = val + lv%map%sizeof() if (allocated(lv%sm)) val = val + lv%sm%sizeof() if (allocated(lv%sm2a)) val = val + lv%sm2a%sizeof() end function d_base_onelev_sizeof @@ -484,16 +484,18 @@ contains end subroutine d_base_onelev_default - subroutine d_base_onelev_bld_tprol(lv,a,desc_a,ilaggr,nlaggr,op_prol,info) + subroutine d_base_onelev_bld_tprol(lv,a,desc_a,& + & ilaggr,nlaggr,op_prol,ag_data,info) implicit none class(mld_d_onelev_type), intent(inout), target :: lv type(psb_dspmat_type), intent(in) :: a type(psb_desc_type), intent(in) :: desc_a integer(psb_lpk_), allocatable, intent(out) :: ilaggr(:),nlaggr(:) type(psb_ldspmat_type), intent(out) :: op_prol + type(mld_daggr_data), intent(in) :: ag_data integer(psb_ipk_), intent(out) :: info - call lv%aggr%bld_tprol(lv%parms,a,desc_a,ilaggr,nlaggr,op_prol,info) + call lv%aggr%bld_tprol(lv%parms,ag_data,a,desc_a,ilaggr,nlaggr,op_prol,info) end subroutine d_base_onelev_bld_tprol diff --git a/mlprec/mld_d_prec_type.f90 b/mlprec/mld_d_prec_type.f90 index 4ce187f9..07a993be 100644 --- a/mlprec/mld_d_prec_type.f90 +++ b/mlprec/mld_d_prec_type.f90 @@ -85,18 +85,8 @@ module mld_d_prec_type integer, parameter, private :: wv_size_=4 type, extends(psb_dprec_type) :: mld_dprec_type - ! integer(psb_ipk_) :: ictxt ! Now it's in the PSBLAS prec. - ! - ! Aggregation defaults: - ! - ! 1. min_coarse_size = 0 Default target size will be computed as 40*(N_fine)**(1./3.) - integer(psb_ipk_) :: min_coarse_size = izero - ! 2. maximum number of levels. Defaults to 20 - integer(psb_ipk_) :: max_levs = 20_psb_ipk_ - ! 3. min_cr_ratio = 1.5 - real(psb_dpk_) :: min_cr_ratio = 1.5_psb_dpk_ - real(psb_dpk_) :: op_complexity = dzero - real(psb_dpk_) :: avg_cr = dzero + ! integer(psb_ipk_) :: ictxt ! Now it's in the PSBLAS prec. + type(mld_daggr_data) :: ag_data ! ! Number of outer sweeps. Sometimes 2 V-cycles may be better than 1 W-cycle. ! @@ -445,7 +435,7 @@ contains class(mld_dprec_type), intent(in) :: prec real(psb_dpk_) :: val - val = prec%op_complexity + val = prec%ag_data%op_complexity end function mld_d_get_compl @@ -480,7 +470,7 @@ contains call psb_sum(ictxt,num) call psb_sum(ictxt,den) end if - prec%op_complexity = num/den + prec%ag_data%op_complexity = num/den end subroutine mld_d_cmp_compl ! @@ -492,7 +482,7 @@ contains class(mld_dprec_type), intent(in) :: prec real(psb_dpk_) :: val - val = prec%avg_cr + val = prec%ag_data%avg_cr end function mld_d_get_avg_cr @@ -517,7 +507,7 @@ contains avgcr = avgcr / (nl-1) end if call psb_sum(ictxt,avgcr) - prec%avg_cr = avgcr/np + prec%ag_data%avg_cr = avgcr/np end subroutine mld_d_cmp_avg_cr ! @@ -722,14 +712,15 @@ contains end subroutine mld_d_apply1v - subroutine mld_d_dump(prec,info,istart,iend,prefix,head,ac,rp,smoother,solver,global_num) + subroutine mld_d_dump(prec,info,istart,iend,prefix,head,ac,rp,smoother,solver,tprol,& + & global_num) implicit none class(mld_dprec_type), intent(in) :: prec integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(in), optional :: istart, iend character(len=*), intent(in), optional :: prefix, head - logical, optional, intent(in) :: smoother, solver,ac, rp, global_num + logical, optional, intent(in) :: smoother, solver,ac, rp, tprol, global_num integer(psb_ipk_) :: i, j, il1, iln, lname, lev integer(psb_ipk_) :: icontxt,iam, np character(len=80) :: prefix_ @@ -750,7 +741,8 @@ contains do lev=il1, iln call prec%precv(lev)%dump(lev,info,prefix=prefix,head=head,& - & ac=ac,smoother=smoother,solver=solver,rp=rp,global_num=global_num) + & ac=ac,smoother=smoother,solver=solver,rp=rp,tprol=tprol, & + & global_num=global_num) end do end subroutine mld_d_dump @@ -801,13 +793,9 @@ contains info = psb_success_ select type(pout => precout) class is (mld_dprec_type) - pout%ictxt = prec%ictxt - pout%max_levs = prec%max_levs - pout%min_coarse_size = prec%min_coarse_size - pout%min_cr_ratio = prec%min_cr_ratio - pout%outer_sweeps = prec%outer_sweeps - pout%op_complexity = prec%op_complexity - pout%avg_cr = prec%avg_cr + pout%ictxt = prec%ictxt + pout%ag_data = prec%ag_data + pout%outer_sweeps = prec%outer_sweeps if (allocated(prec%precv)) then ln = size(prec%precv) allocate(pout%precv(ln),stat=info) @@ -853,7 +841,10 @@ contains !!$ return endif end if - + b%ictxt = prec%ictxt + b%ag_data = prec%ag_data + b%outer_sweeps = prec%outer_sweeps + call move_alloc(prec%precv,b%precv) ! Fix the pointers except on level 1. do i=2, size(b%precv) diff --git a/mlprec/mld_d_symdec_aggregator_mod.f90 b/mlprec/mld_d_symdec_aggregator_mod.f90 index 2ec8f933..be247126 100644 --- a/mlprec/mld_d_symdec_aggregator_mod.f90 +++ b/mlprec/mld_d_symdec_aggregator_mod.f90 @@ -62,12 +62,14 @@ module mld_d_symdec_aggregator_mod interface - subroutine mld_d_symdec_aggregator_build_tprol(ag,parms,a,desc_a,ilaggr,nlaggr,op_prol,info) + subroutine mld_d_symdec_aggregator_build_tprol(ag,parms,ag_data,& + & a,desc_a,ilaggr,nlaggr,op_prol,info) import :: mld_d_symdec_aggregator_type, psb_desc_type, psb_dspmat_type, psb_dpk_, & - & psb_ipk_, psb_lpk_, psb_ldspmat_type, mld_dml_parms + & psb_ipk_, psb_lpk_, psb_ldspmat_type, mld_dml_parms, mld_daggr_data implicit none class(mld_d_symdec_aggregator_type), target, intent(inout) :: ag type(mld_dml_parms), intent(inout) :: parms + type(mld_daggr_data), intent(in) :: ag_data type(psb_dspmat_type), intent(in) :: a type(psb_desc_type), intent(in) :: desc_a integer(psb_lpk_), allocatable, intent(out) :: ilaggr(:),nlaggr(:) diff --git a/mlprec/mld_s_base_aggregator_mod.f90 b/mlprec/mld_s_base_aggregator_mod.f90 index 27aafd3d..7209eeb1 100644 --- a/mlprec/mld_s_base_aggregator_mod.f90 +++ b/mlprec/mld_s_base_aggregator_mod.f90 @@ -41,11 +41,11 @@ ! module mld_s_base_aggregator_mod - use mld_base_prec_type, only : mld_sml_parms + use mld_base_prec_type, only : mld_sml_parms, mld_saggr_data use psb_base_mod, only : psb_sspmat_type, psb_lsspmat_type, psb_s_vect_type, & & psb_s_base_vect_type, psb_slinmap_type, psb_spk_, & & psb_ipk_, psb_epk_, psb_lpk_, psb_desc_type, psb_i_base_vect_type, & - & psb_erractionsave, psb_error_handler, psb_success_ + & psb_erractionsave, psb_error_handler, psb_success_, psb_toupper ! ! ! @@ -79,7 +79,8 @@ module mld_s_base_aggregator_mod !! cseti, csetr, csetc - Set internal parameters, if any ! type mld_s_base_aggregator_type - + ! Do we want to purge explicit zeros when aggregating? + logical :: do_clean_zeros contains procedure, pass(ag) :: bld_tprol => mld_s_base_aggregator_build_tprol procedure, pass(ag) :: mat_asb => mld_s_base_aggregator_mat_asb @@ -96,6 +97,19 @@ module mld_s_base_aggregator_mod generic, public :: set => cseti, csetr, csetc end type mld_s_base_aggregator_type + abstract interface + subroutine mld_s_soc_map_bld(iorder,theta,clean_zeros,a,desc_a,nlaggr,ilaggr,info) + import :: psb_sspmat_type, psb_desc_type, psb_spk_, psb_ipk_, psb_lpk_ + implicit none + integer(psb_ipk_), intent(in) :: iorder + logical, intent(in) :: clean_zeros + type(psb_sspmat_type), intent(in) :: a + type(psb_desc_type), intent(in) :: desc_a + real(psb_spk_), intent(in) :: theta + integer(psb_lpk_), allocatable, intent(out) :: ilaggr(:),nlaggr(:) + integer(psb_ipk_), intent(out) :: info + end subroutine mld_s_soc_map_bld + end interface contains @@ -137,7 +151,16 @@ contains character(len=*), intent(in) :: val integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(in), optional :: idx - ! Do nothing + ! Set clean zeros, or do nothing. + select case (psb_toupper(trim(what))) + case('AGGR_CLEAN_ZEROS') + select case (psb_toupper(trim(val))) + case('TRUE','T') + ag%do_clean_zeros = .true. + case('FALSE','F') + ag%do_clean_zeros = .false. + end select + end select info = 0 end subroutine mld_s_base_aggregator_csetc @@ -181,8 +204,8 @@ contains subroutine mld_s_base_aggregator_default(ag) implicit none class(mld_s_base_aggregator_type), intent(inout) :: ag - - ! Here we need do nothing + ! Only one default setting + ag%do_clean_zeros = .true. return end subroutine mld_s_base_aggregator_default @@ -228,9 +251,12 @@ contains !! will contribute, in global numbering. !! Many aggregation produce a binary tentative prolongator, but some !! do not, hence we also need the OP_PROL output. + !! AG_DATA is passed here just in case some of the + !! aggregators need it internally, most of them will ignore. !! !! \param ag The input aggregator object !! \param parms The auxiliary parameters object + !! \param ag_data Auxiliary global aggregation info !! \param a The local matrix part !! \param desc_a The descriptor !! \param ilaggr Output aggregation map @@ -239,11 +265,13 @@ contains !! \param info Return code !! ! - subroutine mld_s_base_aggregator_build_tprol(ag,parms,a,desc_a,ilaggr,nlaggr,op_prol,info) + subroutine mld_s_base_aggregator_build_tprol(ag,parms,ag_data,& + & a,desc_a,ilaggr,nlaggr,op_prol,info) use psb_base_mod implicit none class(mld_s_base_aggregator_type), target, intent(inout) :: ag - type(mld_sml_parms), intent(inout) :: parms + type(mld_sml_parms), intent(inout) :: parms + type(mld_saggr_data), intent(in) :: ag_data type(psb_sspmat_type), intent(in) :: a type(psb_desc_type), intent(in) :: desc_a integer(psb_lpk_), allocatable, intent(out) :: ilaggr(:),nlaggr(:) diff --git a/mlprec/mld_s_dec_aggregator_mod.f90 b/mlprec/mld_s_dec_aggregator_mod.f90 index a3944890..0e8380f1 100644 --- a/mlprec/mld_s_dec_aggregator_mod.f90 +++ b/mlprec/mld_s_dec_aggregator_mod.f90 @@ -48,7 +48,7 @@ module mld_s_dec_aggregator_mod !! \extends mld_s_base_aggregator_mod::mld_s_base_aggregator_type !! !! type, extends(mld_s_base_aggregator_type) :: mld_s_dec_aggregator_type - !! procedure(mld_s_map_bld), nopass, pointer :: map_bld => null() + !! procedure(mld_s_soc_map_bld), nopass, pointer :: soc_map_bld => null() !! end type !! !! This is the simplest aggregation method: starting from the @@ -71,12 +71,12 @@ module mld_s_dec_aggregator_mod !! PSBLAS-based parallel two-level Schwarz preconditioners, Appl. Num. Math. !! 57 (2007), 1181-1196. !! - !! The map_bld method is used inside the implementation of build_tprol + !! The soc_map_bld method is used inside the implementation of build_tprol !! ! ! type, extends(mld_s_base_aggregator_type) :: mld_s_dec_aggregator_type - procedure(mld_s_map_bld), nopass, pointer :: map_bld => null() + procedure(mld_s_soc_map_bld), nopass, pointer :: soc_map_bld => null() contains procedure, pass(ag) :: bld_tprol => mld_s_dec_aggregator_build_tprol @@ -88,28 +88,17 @@ module mld_s_dec_aggregator_mod end type mld_s_dec_aggregator_type - abstract interface - subroutine mld_s_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info) - import :: psb_sspmat_type, psb_desc_type, psb_spk_, psb_ipk_, psb_lpk_ - implicit none - integer(psb_ipk_), intent(in) :: iorder - type(psb_sspmat_type), intent(in) :: a - type(psb_desc_type), intent(in) :: desc_a - real(psb_spk_), intent(in) :: theta - integer(psb_lpk_), allocatable, intent(out) :: ilaggr(:),nlaggr(:) - integer(psb_ipk_), intent(out) :: info - end subroutine mld_s_map_bld - end interface - - procedure(mld_s_map_bld) :: mld_s_soc1_map_bld, mld_s_soc2_map_bld + procedure(mld_s_soc_map_bld) :: mld_s_soc1_map_bld, mld_s_soc2_map_bld interface - subroutine mld_s_dec_aggregator_build_tprol(ag,parms,a,desc_a,ilaggr,nlaggr,op_prol,info) + subroutine mld_s_dec_aggregator_build_tprol(ag,parms,ag_data,& + & a,desc_a,ilaggr,nlaggr,op_prol,info) import :: mld_s_dec_aggregator_type, psb_desc_type, psb_sspmat_type, psb_spk_, & - & psb_ipk_, psb_lpk_, psb_lsspmat_type, mld_sml_parms + & psb_ipk_, psb_lpk_, psb_lsspmat_type, mld_sml_parms, mld_saggr_data implicit none class(mld_s_dec_aggregator_type), target, intent(inout) :: ag type(mld_sml_parms), intent(inout) :: parms + type(mld_saggr_data), intent(in) :: ag_data type(psb_sspmat_type), intent(in) :: a type(psb_desc_type), intent(in) :: desc_a integer(psb_lpk_), allocatable, intent(out) :: ilaggr(:),nlaggr(:) @@ -147,14 +136,14 @@ contains select case(parms%aggr_type) case (mld_noalg_) - ag%map_bld => null() + ag%soc_map_bld => null() case (mld_soc1_) - ag%map_bld => mld_s_soc1_map_bld + ag%soc_map_bld => mld_s_soc1_map_bld case (mld_soc2_) - ag%map_bld => mld_s_soc2_map_bld + ag%soc_map_bld => mld_s_soc2_map_bld case default write(0,*) 'Unknown aggregation type, defaulting to SOC1' - ag%map_bld => mld_s_soc1_map_bld + ag%soc_map_bld => mld_s_soc1_map_bld end select return @@ -165,7 +154,8 @@ contains implicit none class(mld_s_dec_aggregator_type), intent(inout) :: ag - ag%map_bld => mld_s_soc1_map_bld + call ag%mld_s_base_aggregator_type%default() + ag%soc_map_bld => mld_s_soc1_map_bld return end subroutine mld_s_dec_aggregator_default diff --git a/mlprec/mld_s_inner_mod.f90 b/mlprec/mld_s_inner_mod.f90 index 4dd44d02..29e01457 100644 --- a/mlprec/mld_s_inner_mod.f90 +++ b/mlprec/mld_s_inner_mod.f90 @@ -97,32 +97,6 @@ module mld_s_inner_mod end subroutine mld_smlprec_aply_vect end interface mld_mlprec_aply - interface mld_aggrmap_bld - subroutine mld_s_lev_aggrmap_bld(p,a,desc_a,ilaggr,nlaggr,op_prol,info) - import :: psb_sspmat_type, psb_desc_type, psb_spk_, psb_ipk_, psb_lpk_, psb_lsspmat_type - import :: mld_s_onelev_type - implicit none - type(mld_s_onelev_type), intent(inout), target :: p - type(psb_sspmat_type), intent(in) :: a - type(psb_desc_type), intent(in) :: desc_a - integer(psb_lpk_), intent(out) :: ilaggr(:),nlaggr(:) - type(psb_lsspmat_type), intent(out) :: op_prol - integer(psb_ipk_), intent(out) :: info - end subroutine mld_s_lev_aggrmap_bld - subroutine mld_saggrmap_bld(aggr_type,iorder,theta,a,desc_a,ilaggr,nlaggr,op_prol,info) - import :: psb_sspmat_type, psb_desc_type, psb_spk_, psb_ipk_, psb_lpk_, psb_lsspmat_type - implicit none - integer(psb_ipk_), intent(in) :: iorder - integer(psb_ipk_), intent(in) :: aggr_type - real(psb_spk_), intent(in) :: theta - type(psb_sspmat_type), intent(in) :: a - type(psb_desc_type), intent(in) :: desc_a - integer(psb_lpk_), allocatable, intent(out) :: ilaggr(:),nlaggr(:) - type(psb_lsspmat_type), intent(out) :: op_prol - integer(psb_ipk_), intent(out) :: info - end subroutine mld_saggrmap_bld - end interface mld_aggrmap_bld - interface mld_map_to_tprol subroutine mld_s_map_to_tprol(desc_a,ilaggr,nlaggr,op_prol,info) use psb_base_mod, only : psb_sspmat_type, psb_desc_type, psb_spk_, psb_ipk_, psb_lpk_, psb_lsspmat_type diff --git a/mlprec/mld_s_onelev_mod.f90 b/mlprec/mld_s_onelev_mod.f90 index e9667402..d743b366 100644 --- a/mlprec/mld_s_onelev_mod.f90 +++ b/mlprec/mld_s_onelev_mod.f90 @@ -381,7 +381,7 @@ interface interface subroutine mld_s_base_onelev_dump(lv,level,info,prefix,head,ac,rp,smoother,& - & solver,global_num) + & solver,tprol,global_num) import :: psb_sspmat_type, psb_s_vect_type, psb_s_base_vect_type, & & psb_slinmap_type, psb_spk_, mld_s_onelev_type, & & psb_ipk_, psb_epk_, psb_desc_type @@ -390,7 +390,7 @@ interface integer(psb_ipk_), intent(in) :: level integer(psb_ipk_), intent(out) :: info character(len=*), intent(in), optional :: prefix, head - logical, optional, intent(in) :: ac, rp, smoother, solver, global_num + logical, optional, intent(in) :: ac, rp, smoother, solver, tprol, global_num end subroutine mld_s_base_onelev_dump end interface @@ -484,16 +484,18 @@ contains end subroutine s_base_onelev_default - subroutine s_base_onelev_bld_tprol(lv,a,desc_a,ilaggr,nlaggr,op_prol,info) + subroutine s_base_onelev_bld_tprol(lv,a,desc_a,& + & ilaggr,nlaggr,op_prol,ag_data,info) implicit none class(mld_s_onelev_type), intent(inout), target :: lv type(psb_sspmat_type), intent(in) :: a type(psb_desc_type), intent(in) :: desc_a integer(psb_lpk_), allocatable, intent(out) :: ilaggr(:),nlaggr(:) type(psb_lsspmat_type), intent(out) :: op_prol + type(mld_saggr_data), intent(in) :: ag_data integer(psb_ipk_), intent(out) :: info - call lv%aggr%bld_tprol(lv%parms,a,desc_a,ilaggr,nlaggr,op_prol,info) + call lv%aggr%bld_tprol(lv%parms,ag_data,a,desc_a,ilaggr,nlaggr,op_prol,info) end subroutine s_base_onelev_bld_tprol diff --git a/mlprec/mld_s_prec_type.f90 b/mlprec/mld_s_prec_type.f90 index b685433e..37c40eb3 100644 --- a/mlprec/mld_s_prec_type.f90 +++ b/mlprec/mld_s_prec_type.f90 @@ -85,18 +85,8 @@ module mld_s_prec_type integer, parameter, private :: wv_size_=4 type, extends(psb_sprec_type) :: mld_sprec_type - ! integer(psb_ipk_) :: ictxt ! Now it's in the PSBLAS prec. - ! - ! Aggregation defaults: - ! - ! 1. min_coarse_size = 0 Default target size will be computed as 40*(N_fine)**(1./3.) - integer(psb_ipk_) :: min_coarse_size = izero - ! 2. maximum number of levels. Defaults to 20 - integer(psb_ipk_) :: max_levs = 20_psb_ipk_ - ! 3. min_cr_ratio = 1.5 - real(psb_spk_) :: min_cr_ratio = 1.5_psb_spk_ - real(psb_spk_) :: op_complexity = szero - real(psb_spk_) :: avg_cr = szero + ! integer(psb_ipk_) :: ictxt ! Now it's in the PSBLAS prec. + type(mld_saggr_data) :: ag_data ! ! Number of outer sweeps. Sometimes 2 V-cycles may be better than 1 W-cycle. ! @@ -445,7 +435,7 @@ contains class(mld_sprec_type), intent(in) :: prec real(psb_spk_) :: val - val = prec%op_complexity + val = prec%ag_data%op_complexity end function mld_s_get_compl @@ -480,7 +470,7 @@ contains call psb_sum(ictxt,num) call psb_sum(ictxt,den) end if - prec%op_complexity = num/den + prec%ag_data%op_complexity = num/den end subroutine mld_s_cmp_compl ! @@ -492,7 +482,7 @@ contains class(mld_sprec_type), intent(in) :: prec real(psb_spk_) :: val - val = prec%avg_cr + val = prec%ag_data%avg_cr end function mld_s_get_avg_cr @@ -517,7 +507,7 @@ contains avgcr = avgcr / (nl-1) end if call psb_sum(ictxt,avgcr) - prec%avg_cr = avgcr/np + prec%ag_data%avg_cr = avgcr/np end subroutine mld_s_cmp_avg_cr ! @@ -722,14 +712,15 @@ contains end subroutine mld_s_apply1v - subroutine mld_s_dump(prec,info,istart,iend,prefix,head,ac,rp,smoother,solver,global_num) + subroutine mld_s_dump(prec,info,istart,iend,prefix,head,ac,rp,smoother,solver,tprol,& + & global_num) implicit none class(mld_sprec_type), intent(in) :: prec integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(in), optional :: istart, iend character(len=*), intent(in), optional :: prefix, head - logical, optional, intent(in) :: smoother, solver,ac, rp, global_num + logical, optional, intent(in) :: smoother, solver,ac, rp, tprol, global_num integer(psb_ipk_) :: i, j, il1, iln, lname, lev integer(psb_ipk_) :: icontxt,iam, np character(len=80) :: prefix_ @@ -750,7 +741,8 @@ contains do lev=il1, iln call prec%precv(lev)%dump(lev,info,prefix=prefix,head=head,& - & ac=ac,smoother=smoother,solver=solver,rp=rp,global_num=global_num) + & ac=ac,smoother=smoother,solver=solver,rp=rp,tprol=tprol, & + & global_num=global_num) end do end subroutine mld_s_dump @@ -801,13 +793,9 @@ contains info = psb_success_ select type(pout => precout) class is (mld_sprec_type) - pout%ictxt = prec%ictxt - pout%max_levs = prec%max_levs - pout%min_coarse_size = prec%min_coarse_size - pout%min_cr_ratio = prec%min_cr_ratio - pout%outer_sweeps = prec%outer_sweeps - pout%op_complexity = prec%op_complexity - pout%avg_cr = prec%avg_cr + pout%ictxt = prec%ictxt + pout%ag_data = prec%ag_data + pout%outer_sweeps = prec%outer_sweeps if (allocated(prec%precv)) then ln = size(prec%precv) allocate(pout%precv(ln),stat=info) @@ -853,7 +841,10 @@ contains !!$ return endif end if - + b%ictxt = prec%ictxt + b%ag_data = prec%ag_data + b%outer_sweeps = prec%outer_sweeps + call move_alloc(prec%precv,b%precv) ! Fix the pointers except on level 1. do i=2, size(b%precv) diff --git a/mlprec/mld_s_symdec_aggregator_mod.f90 b/mlprec/mld_s_symdec_aggregator_mod.f90 index 0a750908..026abbe6 100644 --- a/mlprec/mld_s_symdec_aggregator_mod.f90 +++ b/mlprec/mld_s_symdec_aggregator_mod.f90 @@ -62,12 +62,14 @@ module mld_s_symdec_aggregator_mod interface - subroutine mld_s_symdec_aggregator_build_tprol(ag,parms,a,desc_a,ilaggr,nlaggr,op_prol,info) + subroutine mld_s_symdec_aggregator_build_tprol(ag,parms,ag_data,& + & a,desc_a,ilaggr,nlaggr,op_prol,info) import :: mld_s_symdec_aggregator_type, psb_desc_type, psb_sspmat_type, psb_spk_, & - & psb_ipk_, psb_lpk_, psb_lsspmat_type, mld_sml_parms + & psb_ipk_, psb_lpk_, psb_lsspmat_type, mld_sml_parms, mld_saggr_data implicit none class(mld_s_symdec_aggregator_type), target, intent(inout) :: ag type(mld_sml_parms), intent(inout) :: parms + type(mld_saggr_data), intent(in) :: ag_data type(psb_sspmat_type), intent(in) :: a type(psb_desc_type), intent(in) :: desc_a integer(psb_lpk_), allocatable, intent(out) :: ilaggr(:),nlaggr(:) diff --git a/mlprec/mld_z_base_aggregator_mod.f90 b/mlprec/mld_z_base_aggregator_mod.f90 index 9e407514..c1cda281 100644 --- a/mlprec/mld_z_base_aggregator_mod.f90 +++ b/mlprec/mld_z_base_aggregator_mod.f90 @@ -41,11 +41,11 @@ ! module mld_z_base_aggregator_mod - use mld_base_prec_type, only : mld_dml_parms + use mld_base_prec_type, only : mld_dml_parms, mld_daggr_data use psb_base_mod, only : psb_zspmat_type, psb_lzspmat_type, psb_z_vect_type, & & psb_z_base_vect_type, psb_zlinmap_type, psb_dpk_, & & psb_ipk_, psb_epk_, psb_lpk_, psb_desc_type, psb_i_base_vect_type, & - & psb_erractionsave, psb_error_handler, psb_success_ + & psb_erractionsave, psb_error_handler, psb_success_, psb_toupper ! ! ! @@ -79,7 +79,8 @@ module mld_z_base_aggregator_mod !! cseti, csetr, csetc - Set internal parameters, if any ! type mld_z_base_aggregator_type - + ! Do we want to purge explicit zeros when aggregating? + logical :: do_clean_zeros contains procedure, pass(ag) :: bld_tprol => mld_z_base_aggregator_build_tprol procedure, pass(ag) :: mat_asb => mld_z_base_aggregator_mat_asb @@ -96,6 +97,19 @@ module mld_z_base_aggregator_mod generic, public :: set => cseti, csetr, csetc end type mld_z_base_aggregator_type + abstract interface + subroutine mld_z_soc_map_bld(iorder,theta,clean_zeros,a,desc_a,nlaggr,ilaggr,info) + import :: psb_zspmat_type, psb_desc_type, psb_dpk_, psb_ipk_, psb_lpk_ + implicit none + integer(psb_ipk_), intent(in) :: iorder + logical, intent(in) :: clean_zeros + type(psb_zspmat_type), intent(in) :: a + type(psb_desc_type), intent(in) :: desc_a + real(psb_dpk_), intent(in) :: theta + integer(psb_lpk_), allocatable, intent(out) :: ilaggr(:),nlaggr(:) + integer(psb_ipk_), intent(out) :: info + end subroutine mld_z_soc_map_bld + end interface contains @@ -137,7 +151,16 @@ contains character(len=*), intent(in) :: val integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(in), optional :: idx - ! Do nothing + ! Set clean zeros, or do nothing. + select case (psb_toupper(trim(what))) + case('AGGR_CLEAN_ZEROS') + select case (psb_toupper(trim(val))) + case('TRUE','T') + ag%do_clean_zeros = .true. + case('FALSE','F') + ag%do_clean_zeros = .false. + end select + end select info = 0 end subroutine mld_z_base_aggregator_csetc @@ -181,8 +204,8 @@ contains subroutine mld_z_base_aggregator_default(ag) implicit none class(mld_z_base_aggregator_type), intent(inout) :: ag - - ! Here we need do nothing + ! Only one default setting + ag%do_clean_zeros = .true. return end subroutine mld_z_base_aggregator_default @@ -228,9 +251,12 @@ contains !! will contribute, in global numbering. !! Many aggregation produce a binary tentative prolongator, but some !! do not, hence we also need the OP_PROL output. + !! AG_DATA is passed here just in case some of the + !! aggregators need it internally, most of them will ignore. !! !! \param ag The input aggregator object !! \param parms The auxiliary parameters object + !! \param ag_data Auxiliary global aggregation info !! \param a The local matrix part !! \param desc_a The descriptor !! \param ilaggr Output aggregation map @@ -239,11 +265,13 @@ contains !! \param info Return code !! ! - subroutine mld_z_base_aggregator_build_tprol(ag,parms,a,desc_a,ilaggr,nlaggr,op_prol,info) + subroutine mld_z_base_aggregator_build_tprol(ag,parms,ag_data,& + & a,desc_a,ilaggr,nlaggr,op_prol,info) use psb_base_mod implicit none class(mld_z_base_aggregator_type), target, intent(inout) :: ag - type(mld_dml_parms), intent(inout) :: parms + type(mld_dml_parms), intent(inout) :: parms + type(mld_daggr_data), intent(in) :: ag_data type(psb_zspmat_type), intent(in) :: a type(psb_desc_type), intent(in) :: desc_a integer(psb_lpk_), allocatable, intent(out) :: ilaggr(:),nlaggr(:) diff --git a/mlprec/mld_z_dec_aggregator_mod.f90 b/mlprec/mld_z_dec_aggregator_mod.f90 index d76345d2..307716b0 100644 --- a/mlprec/mld_z_dec_aggregator_mod.f90 +++ b/mlprec/mld_z_dec_aggregator_mod.f90 @@ -48,7 +48,7 @@ module mld_z_dec_aggregator_mod !! \extends mld_z_base_aggregator_mod::mld_z_base_aggregator_type !! !! type, extends(mld_z_base_aggregator_type) :: mld_z_dec_aggregator_type - !! procedure(mld_z_map_bld), nopass, pointer :: map_bld => null() + !! procedure(mld_z_soc_map_bld), nopass, pointer :: soc_map_bld => null() !! end type !! !! This is the simplest aggregation method: starting from the @@ -71,12 +71,12 @@ module mld_z_dec_aggregator_mod !! PSBLAS-based parallel two-level Schwarz preconditioners, Appl. Num. Math. !! 57 (2007), 1181-1196. !! - !! The map_bld method is used inside the implementation of build_tprol + !! The soc_map_bld method is used inside the implementation of build_tprol !! ! ! type, extends(mld_z_base_aggregator_type) :: mld_z_dec_aggregator_type - procedure(mld_z_map_bld), nopass, pointer :: map_bld => null() + procedure(mld_z_soc_map_bld), nopass, pointer :: soc_map_bld => null() contains procedure, pass(ag) :: bld_tprol => mld_z_dec_aggregator_build_tprol @@ -88,28 +88,17 @@ module mld_z_dec_aggregator_mod end type mld_z_dec_aggregator_type - abstract interface - subroutine mld_z_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info) - import :: psb_zspmat_type, psb_desc_type, psb_dpk_, psb_ipk_, psb_lpk_ - implicit none - integer(psb_ipk_), intent(in) :: iorder - type(psb_zspmat_type), intent(in) :: a - type(psb_desc_type), intent(in) :: desc_a - real(psb_dpk_), intent(in) :: theta - integer(psb_lpk_), allocatable, intent(out) :: ilaggr(:),nlaggr(:) - integer(psb_ipk_), intent(out) :: info - end subroutine mld_z_map_bld - end interface - - procedure(mld_z_map_bld) :: mld_z_soc1_map_bld, mld_z_soc2_map_bld + procedure(mld_z_soc_map_bld) :: mld_z_soc1_map_bld, mld_z_soc2_map_bld interface - subroutine mld_z_dec_aggregator_build_tprol(ag,parms,a,desc_a,ilaggr,nlaggr,op_prol,info) + subroutine mld_z_dec_aggregator_build_tprol(ag,parms,ag_data,& + & a,desc_a,ilaggr,nlaggr,op_prol,info) import :: mld_z_dec_aggregator_type, psb_desc_type, psb_zspmat_type, psb_dpk_, & - & psb_ipk_, psb_lpk_, psb_lzspmat_type, mld_dml_parms + & psb_ipk_, psb_lpk_, psb_lzspmat_type, mld_dml_parms, mld_daggr_data implicit none class(mld_z_dec_aggregator_type), target, intent(inout) :: ag type(mld_dml_parms), intent(inout) :: parms + type(mld_daggr_data), intent(in) :: ag_data type(psb_zspmat_type), intent(in) :: a type(psb_desc_type), intent(in) :: desc_a integer(psb_lpk_), allocatable, intent(out) :: ilaggr(:),nlaggr(:) @@ -147,14 +136,14 @@ contains select case(parms%aggr_type) case (mld_noalg_) - ag%map_bld => null() + ag%soc_map_bld => null() case (mld_soc1_) - ag%map_bld => mld_z_soc1_map_bld + ag%soc_map_bld => mld_z_soc1_map_bld case (mld_soc2_) - ag%map_bld => mld_z_soc2_map_bld + ag%soc_map_bld => mld_z_soc2_map_bld case default write(0,*) 'Unknown aggregation type, defaulting to SOC1' - ag%map_bld => mld_z_soc1_map_bld + ag%soc_map_bld => mld_z_soc1_map_bld end select return @@ -165,7 +154,8 @@ contains implicit none class(mld_z_dec_aggregator_type), intent(inout) :: ag - ag%map_bld => mld_z_soc1_map_bld + call ag%mld_z_base_aggregator_type%default() + ag%soc_map_bld => mld_z_soc1_map_bld return end subroutine mld_z_dec_aggregator_default diff --git a/mlprec/mld_z_inner_mod.f90 b/mlprec/mld_z_inner_mod.f90 index 4fe6c8e6..3e46b6f0 100644 --- a/mlprec/mld_z_inner_mod.f90 +++ b/mlprec/mld_z_inner_mod.f90 @@ -97,32 +97,6 @@ module mld_z_inner_mod end subroutine mld_zmlprec_aply_vect end interface mld_mlprec_aply - interface mld_aggrmap_bld - subroutine mld_z_lev_aggrmap_bld(p,a,desc_a,ilaggr,nlaggr,op_prol,info) - import :: psb_zspmat_type, psb_desc_type, psb_dpk_, psb_ipk_, psb_lpk_, psb_lzspmat_type - import :: mld_z_onelev_type - implicit none - type(mld_z_onelev_type), intent(inout), target :: p - type(psb_zspmat_type), intent(in) :: a - type(psb_desc_type), intent(in) :: desc_a - integer(psb_lpk_), intent(out) :: ilaggr(:),nlaggr(:) - type(psb_lzspmat_type), intent(out) :: op_prol - integer(psb_ipk_), intent(out) :: info - end subroutine mld_z_lev_aggrmap_bld - subroutine mld_zaggrmap_bld(aggr_type,iorder,theta,a,desc_a,ilaggr,nlaggr,op_prol,info) - import :: psb_zspmat_type, psb_desc_type, psb_dpk_, psb_ipk_, psb_lpk_, psb_lzspmat_type - implicit none - integer(psb_ipk_), intent(in) :: iorder - integer(psb_ipk_), intent(in) :: aggr_type - real(psb_dpk_), intent(in) :: theta - type(psb_zspmat_type), intent(in) :: a - type(psb_desc_type), intent(in) :: desc_a - integer(psb_lpk_), allocatable, intent(out) :: ilaggr(:),nlaggr(:) - type(psb_lzspmat_type), intent(out) :: op_prol - integer(psb_ipk_), intent(out) :: info - end subroutine mld_zaggrmap_bld - end interface mld_aggrmap_bld - interface mld_map_to_tprol subroutine mld_z_map_to_tprol(desc_a,ilaggr,nlaggr,op_prol,info) use psb_base_mod, only : psb_zspmat_type, psb_desc_type, psb_dpk_, psb_ipk_, psb_lpk_, psb_lzspmat_type diff --git a/mlprec/mld_z_onelev_mod.f90 b/mlprec/mld_z_onelev_mod.f90 index 124fc8d2..a560ed5f 100644 --- a/mlprec/mld_z_onelev_mod.f90 +++ b/mlprec/mld_z_onelev_mod.f90 @@ -381,7 +381,7 @@ interface interface subroutine mld_z_base_onelev_dump(lv,level,info,prefix,head,ac,rp,smoother,& - & solver,global_num) + & solver,tprol,global_num) import :: psb_zspmat_type, psb_z_vect_type, psb_z_base_vect_type, & & psb_zlinmap_type, psb_dpk_, mld_z_onelev_type, & & psb_ipk_, psb_epk_, psb_desc_type @@ -390,7 +390,7 @@ interface integer(psb_ipk_), intent(in) :: level integer(psb_ipk_), intent(out) :: info character(len=*), intent(in), optional :: prefix, head - logical, optional, intent(in) :: ac, rp, smoother, solver, global_num + logical, optional, intent(in) :: ac, rp, smoother, solver, tprol, global_num end subroutine mld_z_base_onelev_dump end interface @@ -484,16 +484,18 @@ contains end subroutine z_base_onelev_default - subroutine z_base_onelev_bld_tprol(lv,a,desc_a,ilaggr,nlaggr,op_prol,info) + subroutine z_base_onelev_bld_tprol(lv,a,desc_a,& + & ilaggr,nlaggr,op_prol,ag_data,info) implicit none class(mld_z_onelev_type), intent(inout), target :: lv type(psb_zspmat_type), intent(in) :: a type(psb_desc_type), intent(in) :: desc_a integer(psb_lpk_), allocatable, intent(out) :: ilaggr(:),nlaggr(:) type(psb_lzspmat_type), intent(out) :: op_prol + type(mld_daggr_data), intent(in) :: ag_data integer(psb_ipk_), intent(out) :: info - call lv%aggr%bld_tprol(lv%parms,a,desc_a,ilaggr,nlaggr,op_prol,info) + call lv%aggr%bld_tprol(lv%parms,ag_data,a,desc_a,ilaggr,nlaggr,op_prol,info) end subroutine z_base_onelev_bld_tprol diff --git a/mlprec/mld_z_prec_type.f90 b/mlprec/mld_z_prec_type.f90 index 714d6db3..fda7429a 100644 --- a/mlprec/mld_z_prec_type.f90 +++ b/mlprec/mld_z_prec_type.f90 @@ -85,18 +85,8 @@ module mld_z_prec_type integer, parameter, private :: wv_size_=4 type, extends(psb_zprec_type) :: mld_zprec_type - ! integer(psb_ipk_) :: ictxt ! Now it's in the PSBLAS prec. - ! - ! Aggregation defaults: - ! - ! 1. min_coarse_size = 0 Default target size will be computed as 40*(N_fine)**(1./3.) - integer(psb_ipk_) :: min_coarse_size = izero - ! 2. maximum number of levels. Defaults to 20 - integer(psb_ipk_) :: max_levs = 20_psb_ipk_ - ! 3. min_cr_ratio = 1.5 - real(psb_dpk_) :: min_cr_ratio = 1.5_psb_dpk_ - real(psb_dpk_) :: op_complexity = dzero - real(psb_dpk_) :: avg_cr = dzero + ! integer(psb_ipk_) :: ictxt ! Now it's in the PSBLAS prec. + type(mld_daggr_data) :: ag_data ! ! Number of outer sweeps. Sometimes 2 V-cycles may be better than 1 W-cycle. ! @@ -445,7 +435,7 @@ contains class(mld_zprec_type), intent(in) :: prec complex(psb_dpk_) :: val - val = prec%op_complexity + val = prec%ag_data%op_complexity end function mld_z_get_compl @@ -480,7 +470,7 @@ contains call psb_sum(ictxt,num) call psb_sum(ictxt,den) end if - prec%op_complexity = num/den + prec%ag_data%op_complexity = num/den end subroutine mld_z_cmp_compl ! @@ -492,7 +482,7 @@ contains class(mld_zprec_type), intent(in) :: prec complex(psb_dpk_) :: val - val = prec%avg_cr + val = prec%ag_data%avg_cr end function mld_z_get_avg_cr @@ -517,7 +507,7 @@ contains avgcr = avgcr / (nl-1) end if call psb_sum(ictxt,avgcr) - prec%avg_cr = avgcr/np + prec%ag_data%avg_cr = avgcr/np end subroutine mld_z_cmp_avg_cr ! @@ -722,14 +712,15 @@ contains end subroutine mld_z_apply1v - subroutine mld_z_dump(prec,info,istart,iend,prefix,head,ac,rp,smoother,solver,global_num) + subroutine mld_z_dump(prec,info,istart,iend,prefix,head,ac,rp,smoother,solver,tprol,& + & global_num) implicit none class(mld_zprec_type), intent(in) :: prec integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(in), optional :: istart, iend character(len=*), intent(in), optional :: prefix, head - logical, optional, intent(in) :: smoother, solver,ac, rp, global_num + logical, optional, intent(in) :: smoother, solver,ac, rp, tprol, global_num integer(psb_ipk_) :: i, j, il1, iln, lname, lev integer(psb_ipk_) :: icontxt,iam, np character(len=80) :: prefix_ @@ -750,7 +741,8 @@ contains do lev=il1, iln call prec%precv(lev)%dump(lev,info,prefix=prefix,head=head,& - & ac=ac,smoother=smoother,solver=solver,rp=rp,global_num=global_num) + & ac=ac,smoother=smoother,solver=solver,rp=rp,tprol=tprol, & + & global_num=global_num) end do end subroutine mld_z_dump @@ -801,13 +793,9 @@ contains info = psb_success_ select type(pout => precout) class is (mld_zprec_type) - pout%ictxt = prec%ictxt - pout%max_levs = prec%max_levs - pout%min_coarse_size = prec%min_coarse_size - pout%min_cr_ratio = prec%min_cr_ratio - pout%outer_sweeps = prec%outer_sweeps - pout%op_complexity = prec%op_complexity - pout%avg_cr = prec%avg_cr + pout%ictxt = prec%ictxt + pout%ag_data = prec%ag_data + pout%outer_sweeps = prec%outer_sweeps if (allocated(prec%precv)) then ln = size(prec%precv) allocate(pout%precv(ln),stat=info) @@ -853,7 +841,10 @@ contains !!$ return endif end if - + b%ictxt = prec%ictxt + b%ag_data = prec%ag_data + b%outer_sweeps = prec%outer_sweeps + call move_alloc(prec%precv,b%precv) ! Fix the pointers except on level 1. do i=2, size(b%precv) diff --git a/mlprec/mld_z_symdec_aggregator_mod.f90 b/mlprec/mld_z_symdec_aggregator_mod.f90 index 4da6fe6f..9fa7be0b 100644 --- a/mlprec/mld_z_symdec_aggregator_mod.f90 +++ b/mlprec/mld_z_symdec_aggregator_mod.f90 @@ -62,12 +62,14 @@ module mld_z_symdec_aggregator_mod interface - subroutine mld_z_symdec_aggregator_build_tprol(ag,parms,a,desc_a,ilaggr,nlaggr,op_prol,info) + subroutine mld_z_symdec_aggregator_build_tprol(ag,parms,ag_data,& + & a,desc_a,ilaggr,nlaggr,op_prol,info) import :: mld_z_symdec_aggregator_type, psb_desc_type, psb_zspmat_type, psb_dpk_, & - & psb_ipk_, psb_lpk_, psb_lzspmat_type, mld_dml_parms + & psb_ipk_, psb_lpk_, psb_lzspmat_type, mld_dml_parms, mld_daggr_data implicit none class(mld_z_symdec_aggregator_type), target, intent(inout) :: ag type(mld_dml_parms), intent(inout) :: parms + type(mld_daggr_data), intent(in) :: ag_data type(psb_zspmat_type), intent(in) :: a type(psb_desc_type), intent(in) :: desc_a integer(psb_lpk_), allocatable, intent(out) :: ilaggr(:),nlaggr(:)