diff --git a/mlprec/impl/mld_daggrmat_asb.f90 b/mlprec/impl/mld_daggrmat_asb.f90 index 82613c50..69965097 100644 --- a/mlprec/impl/mld_daggrmat_asb.f90 +++ b/mlprec/impl/mld_daggrmat_asb.f90 @@ -98,7 +98,7 @@ ! info - integer, output. ! Error code. ! -subroutine mld_daggrmat_asb(a,desc_a,ilaggr,nlaggr,p,info) +subroutine mld_daggrmat_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_restr,info) use psb_base_mod use mld_d_inner_mod, mld_protect_name => mld_daggrmat_asb @@ -109,11 +109,12 @@ subroutine mld_daggrmat_asb(a,desc_a,ilaggr,nlaggr,p,info) type(psb_dspmat_type), intent(in) :: a type(psb_desc_type), intent(in) :: desc_a integer(psb_ipk_), intent(inout) :: ilaggr(:), nlaggr(:) - type(mld_d_onelev_type), intent(inout), target :: p + type(mld_dml_parms), intent(inout) :: parms +!!$ type(mld_d_onelev_type), intent(inout), target :: p integer(psb_ipk_), intent(out) :: info + type(psb_dspmat_type), intent(inout) :: ac, op_prol,op_restr ! Local variables - type(psb_dspmat_type) :: ac, op_prol,op_restr type(psb_d_coo_sparse_mat) :: acoo, bcoo type(psb_d_csr_sparse_mat) :: acsr1 integer(psb_ipk_) :: nzl,ntaggr, err_act @@ -133,26 +134,26 @@ subroutine mld_daggrmat_asb(a,desc_a,ilaggr,nlaggr,p,info) call psb_info(ictxt, me, np) - select case (p%parms%aggr_kind) + select case (parms%aggr_kind) case (mld_no_smooth_) call mld_daggrmat_nosmth_asb(a,desc_a,ilaggr,nlaggr,& - & p%parms,ac,op_prol,op_restr,info) + & parms,ac,op_prol,op_restr,info) case(mld_smooth_prol_) call mld_daggrmat_smth_asb(a,desc_a,ilaggr,nlaggr, & - & p%parms,ac,op_prol,op_restr,info) + & parms,ac,op_prol,op_restr,info) case(mld_biz_prol_) call mld_daggrmat_biz_asb(a,desc_a,ilaggr,nlaggr, & - & p%parms,ac,op_prol,op_restr,info) + & parms,ac,op_prol,op_restr,info) case(mld_min_energy_) call mld_daggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr, & - & p%parms,ac,op_prol,op_restr,info) + & parms,ac,op_prol,op_restr,info) case default info = psb_err_internal_error_ @@ -166,113 +167,113 @@ subroutine mld_daggrmat_asb(a,desc_a,ilaggr,nlaggr,p,info) end if - - ntaggr = sum(nlaggr) - - select case(p%parms%coarse_mat) - - case(mld_distr_mat_) - - call ac%mv_to(bcoo) - if (p%parms%clean_zeros) call bcoo%clean_zeros(info) - nzl = bcoo%get_nzeros() - - if (info == psb_success_) call psb_cdall(ictxt,p%desc_ac,info,nl=nlaggr(me+1)) - if (info == psb_success_) call psb_cdins(nzl,bcoo%ia,bcoo%ja,p%desc_ac,info) - if (info == psb_success_) call psb_cdasb(p%desc_ac,info) - if (info == psb_success_) call psb_glob_to_loc(bcoo%ia(1:nzl),p%desc_ac,info,iact='I') - if (info == psb_success_) call psb_glob_to_loc(bcoo%ja(1:nzl),p%desc_ac,info,iact='I') - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Creating p%desc_ac and converting ac') - goto 9999 - end if - if (debug_level >= psb_debug_outer_) & - & write(debug_unit,*) me,' ',trim(name),& - & 'Assembld aux descr. distr.' - call p%ac%mv_from(bcoo) - - call p%ac%set_nrows(p%desc_ac%get_local_rows()) - call p%ac%set_ncols(p%desc_ac%get_local_cols()) - call p%ac%set_asb() - - if (info /= psb_success_) then - call psb_errpush(psb_err_from_subroutine_,name,a_err='psb_sp_free') - goto 9999 - end if - - if (np>1) then - call op_prol%mv_to(acsr1) - nzl = acsr1%get_nzeros() - call psb_glob_to_loc(acsr1%ja(1:nzl),p%desc_ac,info,'I') - if(info /= psb_success_) then - call psb_errpush(psb_err_from_subroutine_,name,a_err='psb_glob_to_loc') - goto 9999 - end if - call op_prol%mv_from(acsr1) - endif - call op_prol%set_ncols(p%desc_ac%get_local_cols()) - - if (np>1) then - call op_restr%cscnv(info,type='coo',dupl=psb_dupl_add_) - call op_restr%mv_to(acoo) - nzl = acoo%get_nzeros() - if (info == psb_success_) call psb_glob_to_loc(acoo%ia(1:nzl),p%desc_ac,info,'I') - call acoo%set_dupl(psb_dupl_add_) - if (info == psb_success_) call op_restr%mv_from(acoo) - if (info == psb_success_) call op_restr%cscnv(info,type='csr') - if(info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Converting op_restr to local') - goto 9999 - end if - end if - ! - ! Clip to local rows. - ! - call op_restr%set_nrows(p%desc_ac%get_local_rows()) - - if (debug_level >= psb_debug_outer_) & - & write(debug_unit,*) me,' ',trim(name),& - & 'Done ac ' - - case(mld_repl_mat_) - ! - ! - call psb_cdall(ictxt,p%desc_ac,info,mg=ntaggr,repl=.true.) - if (info == psb_success_) call psb_cdasb(p%desc_ac,info) - if ((info == psb_success_).and.p%parms%clean_zeros) call ac%clean_zeros(info) - if (info == psb_success_) & - & call psb_gather(p%ac,ac,p%desc_ac,info,dupl=psb_dupl_add_,keeploc=.false.) - - if (info /= psb_success_) goto 9999 - - case default - info = psb_err_internal_error_ - call psb_errpush(info,name,a_err='invalid mld_coarse_mat_') - goto 9999 - end select - - call p%ac%cscnv(info,type='csr',dupl=psb_dupl_add_) - if(info /= psb_success_) then - call psb_errpush(psb_err_from_subroutine_,name,a_err='spcnv') - goto 9999 - end if - - ! - ! Copy the prolongation/restriction matrices into the descriptor map. - ! op_restr => PR^T i.e. restriction operator - ! op_prol => PR i.e. prolongation operator - ! - - p%map = psb_linmap(psb_map_aggr_,desc_a,& - & p%desc_ac,op_restr,op_prol,ilaggr,nlaggr) - if (info == psb_success_) call op_prol%free() - if (info == psb_success_) call op_restr%free() - if(info /= psb_success_) then - call psb_errpush(psb_err_from_subroutine_,name,a_err='sp_Free') - goto 9999 - end if +!!$ +!!$ ntaggr = sum(nlaggr) +!!$ +!!$ select case(p%parms%coarse_mat) +!!$ +!!$ case(mld_distr_mat_) +!!$ +!!$ call ac%mv_to(bcoo) +!!$ if (p%parms%clean_zeros) call bcoo%clean_zeros(info) +!!$ nzl = bcoo%get_nzeros() +!!$ +!!$ if (info == psb_success_) call psb_cdall(ictxt,p%desc_ac,info,nl=nlaggr(me+1)) +!!$ if (info == psb_success_) call psb_cdins(nzl,bcoo%ia,bcoo%ja,p%desc_ac,info) +!!$ if (info == psb_success_) call psb_cdasb(p%desc_ac,info) +!!$ if (info == psb_success_) call psb_glob_to_loc(bcoo%ia(1:nzl),p%desc_ac,info,iact='I') +!!$ if (info == psb_success_) call psb_glob_to_loc(bcoo%ja(1:nzl),p%desc_ac,info,iact='I') +!!$ if (info /= psb_success_) then +!!$ call psb_errpush(psb_err_internal_error_,name,& +!!$ & a_err='Creating p%desc_ac and converting ac') +!!$ goto 9999 +!!$ end if +!!$ if (debug_level >= psb_debug_outer_) & +!!$ & write(debug_unit,*) me,' ',trim(name),& +!!$ & 'Assembld aux descr. distr.' +!!$ call p%ac%mv_from(bcoo) +!!$ +!!$ call p%ac%set_nrows(p%desc_ac%get_local_rows()) +!!$ call p%ac%set_ncols(p%desc_ac%get_local_cols()) +!!$ call p%ac%set_asb() +!!$ +!!$ if (info /= psb_success_) then +!!$ call psb_errpush(psb_err_from_subroutine_,name,a_err='psb_sp_free') +!!$ goto 9999 +!!$ end if +!!$ +!!$ if (np>1) then +!!$ call op_prol%mv_to(acsr1) +!!$ nzl = acsr1%get_nzeros() +!!$ call psb_glob_to_loc(acsr1%ja(1:nzl),p%desc_ac,info,'I') +!!$ if(info /= psb_success_) then +!!$ call psb_errpush(psb_err_from_subroutine_,name,a_err='psb_glob_to_loc') +!!$ goto 9999 +!!$ end if +!!$ call op_prol%mv_from(acsr1) +!!$ endif +!!$ call op_prol%set_ncols(p%desc_ac%get_local_cols()) +!!$ +!!$ if (np>1) then +!!$ call op_restr%cscnv(info,type='coo',dupl=psb_dupl_add_) +!!$ call op_restr%mv_to(acoo) +!!$ nzl = acoo%get_nzeros() +!!$ if (info == psb_success_) call psb_glob_to_loc(acoo%ia(1:nzl),p%desc_ac,info,'I') +!!$ call acoo%set_dupl(psb_dupl_add_) +!!$ if (info == psb_success_) call op_restr%mv_from(acoo) +!!$ if (info == psb_success_) call op_restr%cscnv(info,type='csr') +!!$ if(info /= psb_success_) then +!!$ call psb_errpush(psb_err_internal_error_,name,& +!!$ & a_err='Converting op_restr to local') +!!$ goto 9999 +!!$ end if +!!$ end if +!!$ ! +!!$ ! Clip to local rows. +!!$ ! +!!$ call op_restr%set_nrows(p%desc_ac%get_local_rows()) +!!$ +!!$ if (debug_level >= psb_debug_outer_) & +!!$ & write(debug_unit,*) me,' ',trim(name),& +!!$ & 'Done ac ' +!!$ +!!$ case(mld_repl_mat_) +!!$ ! +!!$ ! +!!$ call psb_cdall(ictxt,p%desc_ac,info,mg=ntaggr,repl=.true.) +!!$ if (info == psb_success_) call psb_cdasb(p%desc_ac,info) +!!$ if ((info == psb_success_).and.p%parms%clean_zeros) call ac%clean_zeros(info) +!!$ if (info == psb_success_) & +!!$ & call psb_gather(p%ac,ac,p%desc_ac,info,dupl=psb_dupl_add_,keeploc=.false.) +!!$ +!!$ if (info /= psb_success_) goto 9999 +!!$ +!!$ case default +!!$ info = psb_err_internal_error_ +!!$ call psb_errpush(info,name,a_err='invalid mld_coarse_mat_') +!!$ goto 9999 +!!$ end select +!!$ +!!$ call p%ac%cscnv(info,type='csr',dupl=psb_dupl_add_) +!!$ if(info /= psb_success_) then +!!$ call psb_errpush(psb_err_from_subroutine_,name,a_err='spcnv') +!!$ goto 9999 +!!$ end if +!!$ +!!$ ! +!!$ ! Copy the prolongation/restriction matrices into the descriptor map. +!!$ ! op_restr => PR^T i.e. restriction operator +!!$ ! op_prol => PR i.e. prolongation operator +!!$ ! +!!$ +!!$ p%map = psb_linmap(psb_map_aggr_,desc_a,& +!!$ & p%desc_ac,op_restr,op_prol,ilaggr,nlaggr) +!!$ if (info == psb_success_) call op_prol%free() +!!$ if (info == psb_success_) call op_restr%free() +!!$ if(info /= psb_success_) then +!!$ call psb_errpush(psb_err_from_subroutine_,name,a_err='sp_Free') +!!$ goto 9999 +!!$ end if call psb_erractionrestore(err_act) diff --git a/mlprec/impl/mld_dcoarse_bld.f90 b/mlprec/impl/mld_dcoarse_bld.f90 index 9cf4190e..11b368c6 100644 --- a/mlprec/impl/mld_dcoarse_bld.f90 +++ b/mlprec/impl/mld_dcoarse_bld.f90 @@ -83,11 +83,18 @@ subroutine mld_dcoarse_bld(a,desc_a,p,info) integer(psb_mpik_) :: ictxt, np, me integer(psb_ipk_) :: err_act integer(psb_ipk_), allocatable :: ilaggr(:), nlaggr(:) + type(psb_dspmat_type) :: ac, op_prol,op_restr + type(psb_d_coo_sparse_mat) :: acoo, bcoo + type(psb_d_csr_sparse_mat) :: acsr1 + integer(psb_ipk_) :: nzl, ntaggr + integer(psb_ipk_) :: debug_level, debug_unit name='mld_dcoarse_bld' if (psb_get_errstatus().ne.0) return call psb_erractionsave(err_act) - info = psb_success_ + debug_unit = psb_get_debug_unit() + debug_level = psb_get_debug_level() + info = psb_success_ ictxt = desc_a%get_context() call psb_info(ictxt,me,np) @@ -114,7 +121,7 @@ subroutine mld_dcoarse_bld(a,desc_a,p,info) select case(p%parms%aggr_alg) case (mld_dec_aggr_, mld_sym_dec_aggr_) - + ! ! Build a mapping between the row indices of the fine-level matrix ! and the row indices of the coarse-level matrix, according to a decoupled @@ -134,7 +141,7 @@ subroutine mld_dcoarse_bld(a,desc_a,p,info) ! the mapping defined by mld_aggrmap_bld and applying the aggregation ! algorithm specified by p%iprcparm(mld_aggr_kind_) ! - call mld_aggrmat_asb(a,desc_a,ilaggr,nlaggr,p,info) + call mld_daggrmat_asb(a,desc_a,ilaggr,nlaggr,p%parms,ac,op_prol,op_restr,info) if(info /= psb_success_) then call psb_errpush(psb_err_from_subroutine_,name,a_err='mld_aggrmat_asb') @@ -156,7 +163,116 @@ subroutine mld_dcoarse_bld(a,desc_a,p,info) goto 9999 end select - + + + ! Common code refactored here. + + ntaggr = sum(nlaggr) + + select case(p%parms%coarse_mat) + + case(mld_distr_mat_) + + call ac%mv_to(bcoo) + if (p%parms%clean_zeros) call bcoo%clean_zeros(info) + nzl = bcoo%get_nzeros() + + if (info == psb_success_) call psb_cdall(ictxt,p%desc_ac,info,nl=nlaggr(me+1)) + if (info == psb_success_) call psb_cdins(nzl,bcoo%ia,bcoo%ja,p%desc_ac,info) + if (info == psb_success_) call psb_cdasb(p%desc_ac,info) + if (info == psb_success_) call psb_glob_to_loc(bcoo%ia(1:nzl),p%desc_ac,info,iact='I') + if (info == psb_success_) call psb_glob_to_loc(bcoo%ja(1:nzl),p%desc_ac,info,iact='I') + if (info /= psb_success_) then + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Creating p%desc_ac and converting ac') + goto 9999 + end if + if (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),& + & 'Assembld aux descr. distr.' + call p%ac%mv_from(bcoo) + + call p%ac%set_nrows(p%desc_ac%get_local_rows()) + call p%ac%set_ncols(p%desc_ac%get_local_cols()) + call p%ac%set_asb() + + if (info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,a_err='psb_sp_free') + goto 9999 + end if + + if (np>1) then + call op_prol%mv_to(acsr1) + nzl = acsr1%get_nzeros() + call psb_glob_to_loc(acsr1%ja(1:nzl),p%desc_ac,info,'I') + if(info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,a_err='psb_glob_to_loc') + goto 9999 + end if + call op_prol%mv_from(acsr1) + endif + call op_prol%set_ncols(p%desc_ac%get_local_cols()) + + if (np>1) then + call op_restr%cscnv(info,type='coo',dupl=psb_dupl_add_) + call op_restr%mv_to(acoo) + nzl = acoo%get_nzeros() + if (info == psb_success_) call psb_glob_to_loc(acoo%ia(1:nzl),p%desc_ac,info,'I') + call acoo%set_dupl(psb_dupl_add_) + if (info == psb_success_) call op_restr%mv_from(acoo) + if (info == psb_success_) call op_restr%cscnv(info,type='csr') + if(info /= psb_success_) then + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Converting op_restr to local') + goto 9999 + end if + end if + ! + ! Clip to local rows. + ! + call op_restr%set_nrows(p%desc_ac%get_local_rows()) + + if (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),& + & 'Done ac ' + + case(mld_repl_mat_) + ! + ! + call psb_cdall(ictxt,p%desc_ac,info,mg=ntaggr,repl=.true.) + if (info == psb_success_) call psb_cdasb(p%desc_ac,info) + if ((info == psb_success_).and.p%parms%clean_zeros) call ac%clean_zeros(info) + if (info == psb_success_) & + & call psb_gather(p%ac,ac,p%desc_ac,info,dupl=psb_dupl_add_,keeploc=.false.) + + if (info /= psb_success_) goto 9999 + + case default + info = psb_err_internal_error_ + call psb_errpush(info,name,a_err='invalid mld_coarse_mat_') + goto 9999 + end select + + call p%ac%cscnv(info,type='csr',dupl=psb_dupl_add_) + if(info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,a_err='spcnv') + goto 9999 + end if + + ! + ! Copy the prolongation/restriction matrices into the descriptor map. + ! op_restr => PR^T i.e. restriction operator + ! op_prol => PR i.e. prolongation operator + ! + + p%map = psb_linmap(psb_map_aggr_,desc_a,& + & p%desc_ac,op_restr,op_prol,ilaggr,nlaggr) + if (info == psb_success_) call op_prol%free() + if (info == psb_success_) call op_restr%free() + if(info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,a_err='sp_Free') + goto 9999 + end if ! ! Fix the base_a and base_desc pointers for handling of residuals. ! This is correct because this routine is only called at levels >=2. diff --git a/mlprec/mld_d_inner_mod.f90 b/mlprec/mld_d_inner_mod.f90 index bf891db2..3c50afd7 100644 --- a/mlprec/mld_d_inner_mod.f90 +++ b/mlprec/mld_d_inner_mod.f90 @@ -165,18 +165,18 @@ module mld_d_inner_mod end interface mld_dec_map_bld - interface mld_aggrmat_asb - subroutine mld_daggrmat_asb(a,desc_a,ilaggr,nlaggr,p,info) - use psb_base_mod, only : psb_dspmat_type, psb_desc_type, psb_dpk_, psb_ipk_ - use mld_d_prec_type, only : mld_d_onelev_type - implicit none - type(psb_dspmat_type), intent(in) :: a - type(psb_desc_type), intent(in) :: desc_a - integer(psb_ipk_), intent(inout) :: ilaggr(:), nlaggr(:) - type(mld_d_onelev_type), intent(inout), target :: p - integer(psb_ipk_), intent(out) :: info - end subroutine mld_daggrmat_asb - end interface mld_aggrmat_asb +!!$ interface mld_aggrmat_asb +!!$ subroutine mld_daggrmat_asb(a,desc_a,ilaggr,nlaggr,p,info) +!!$ use psb_base_mod, only : psb_dspmat_type, psb_desc_type, psb_dpk_, psb_ipk_ +!!$ use mld_d_prec_type, only : mld_d_onelev_type +!!$ implicit none +!!$ type(psb_dspmat_type), intent(in) :: a +!!$ type(psb_desc_type), intent(in) :: desc_a +!!$ integer(psb_ipk_), intent(inout) :: ilaggr(:), nlaggr(:) +!!$ type(mld_d_onelev_type), intent(inout), target :: p +!!$ integer(psb_ipk_), intent(out) :: info +!!$ end subroutine mld_daggrmat_asb +!!$ end interface mld_aggrmat_asb @@ -197,7 +197,7 @@ module mld_d_inner_mod procedure(mld_daggrmat_var_asb) :: mld_daggrmat_nosmth_asb, & & mld_daggrmat_smth_asb, mld_daggrmat_minnrg_asb, & - & mld_daggrmat_biz_asb + & mld_daggrmat_biz_asb, mld_daggrmat_asb end module mld_d_inner_mod