From 788d394f58a1d23c36065ded464769834abc1a93 Mon Sep 17 00:00:00 2001 From: Salvatore Filippone Date: Wed, 11 Apr 2012 15:03:56 +0000 Subject: [PATCH] mld2p4-NewNL: mlprec/impl/mld_caggrmat_asb.f90 mlprec/impl/mld_caggrmat_biz_asb.f90 mlprec/impl/mld_caggrmat_minnrg_asb.f90 mlprec/impl/mld_caggrmat_nosmth_asb.f90 mlprec/impl/mld_caggrmat_smth_asb.f90 mlprec/impl/mld_daggrmat_asb.f90 mlprec/impl/mld_daggrmat_biz_asb.f90 mlprec/impl/mld_daggrmat_minnrg_asb.f90 mlprec/impl/mld_daggrmat_nosmth_asb.f90 mlprec/impl/mld_daggrmat_smth_asb.f90 mlprec/impl/mld_saggrmat_asb.f90 mlprec/impl/mld_saggrmat_biz_asb.f90 mlprec/impl/mld_saggrmat_minnrg_asb.f90 mlprec/impl/mld_saggrmat_nosmth_asb.f90 mlprec/impl/mld_saggrmat_smth_asb.f90 mlprec/impl/mld_zaggrmat_asb.f90 mlprec/impl/mld_zaggrmat_biz_asb.f90 mlprec/impl/mld_zaggrmat_minnrg_asb.f90 mlprec/impl/mld_zaggrmat_nosmth_asb.f90 mlprec/impl/mld_zaggrmat_smth_asb.f90 mlprec/mld_c_inner_mod.f90 mlprec/mld_d_inner_mod.f90 mlprec/mld_s_inner_mod.f90 mlprec/mld_z_inner_mod.f90 A bit of internal cleanup. --- mlprec/impl/mld_caggrmat_asb.f90 | 145 ++++++++++++++++--- mlprec/impl/mld_caggrmat_biz_asb.f90 | 165 ++++----------------- mlprec/impl/mld_caggrmat_minnrg_asb.f90 | 184 +++-------------------- mlprec/impl/mld_caggrmat_nosmth_asb.f90 | 144 +++--------------- mlprec/impl/mld_caggrmat_smth_asb.f90 | 185 ++++++------------------ mlprec/impl/mld_daggrmat_asb.f90 | 145 ++++++++++++++++--- mlprec/impl/mld_daggrmat_biz_asb.f90 | 165 ++++----------------- mlprec/impl/mld_daggrmat_minnrg_asb.f90 | 184 +++-------------------- mlprec/impl/mld_daggrmat_nosmth_asb.f90 | 144 +++--------------- mlprec/impl/mld_daggrmat_smth_asb.f90 | 185 ++++++------------------ mlprec/impl/mld_saggrmat_asb.f90 | 145 ++++++++++++++++--- mlprec/impl/mld_saggrmat_biz_asb.f90 | 165 ++++----------------- mlprec/impl/mld_saggrmat_minnrg_asb.f90 | 184 +++-------------------- mlprec/impl/mld_saggrmat_nosmth_asb.f90 | 144 +++--------------- mlprec/impl/mld_saggrmat_smth_asb.f90 | 185 ++++++------------------ mlprec/impl/mld_zaggrmat_asb.f90 | 145 ++++++++++++++++--- mlprec/impl/mld_zaggrmat_biz_asb.f90 | 165 ++++----------------- mlprec/impl/mld_zaggrmat_minnrg_asb.f90 | 184 +++-------------------- mlprec/impl/mld_zaggrmat_nosmth_asb.f90 | 144 +++--------------- mlprec/impl/mld_zaggrmat_smth_asb.f90 | 185 ++++++------------------ mlprec/mld_c_inner_mod.f90 | 7 +- mlprec/mld_d_inner_mod.f90 | 7 +- mlprec/mld_s_inner_mod.f90 | 7 +- mlprec/mld_z_inner_mod.f90 | 7 +- 24 files changed, 972 insertions(+), 2348 deletions(-) diff --git a/mlprec/impl/mld_caggrmat_asb.f90 b/mlprec/impl/mld_caggrmat_asb.f90 index 69cbc9c4..bfc5ef91 100644 --- a/mlprec/impl/mld_caggrmat_asb.f90 +++ b/mlprec/impl/mld_caggrmat_asb.f90 @@ -113,7 +113,11 @@ subroutine mld_caggrmat_asb(a,desc_a,ilaggr,nlaggr,p,info) integer, intent(out) :: info ! Local variables - type(psb_cspmat_type) :: b, op_prol,op_restr + type(psb_cspmat_type) :: ac, op_prol,op_restr + type(psb_c_coo_sparse_mat) :: acoo, bcoo + type(psb_c_csr_sparse_mat) :: acsr1 + integer :: nzl,ntaggr + integer :: debug_level, debug_unit integer :: ictxt,np,me, err_act character(len=20) :: name @@ -121,6 +125,9 @@ subroutine mld_caggrmat_asb(a,desc_a,ilaggr,nlaggr,p,info) if(psb_get_errstatus().ne.0) return info=psb_success_ call psb_erractionsave(err_act) + debug_unit = psb_get_debug_unit() + debug_level = psb_get_debug_level() + ictxt = desc_a%get_context() @@ -129,43 +136,139 @@ subroutine mld_caggrmat_asb(a,desc_a,ilaggr,nlaggr,p,info) select case (p%parms%aggr_kind) case (mld_no_smooth_) - call mld_caggrmat_nosmth_asb(a,desc_a,ilaggr,nlaggr,p,info) - if(info /= psb_success_) then - call psb_errpush(psb_err_from_subroutine_,name,a_err='mld_aggrmat_nosmth_asb') - goto 9999 - end if + call mld_caggrmat_nosmth_asb(a,desc_a,ilaggr,nlaggr,& + & p%parms,ac,op_prol,op_restr,info) case(mld_smooth_prol_) - call mld_caggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info) - if(info /= psb_success_) then - call psb_errpush(psb_err_from_subroutine_,name,a_err='mld_aggrmat_smth_asb') - goto 9999 - end if + call mld_caggrmat_smth_asb(a,desc_a,ilaggr,nlaggr, & + & p%parms,ac,op_prol,op_restr,info) case(mld_biz_prol_) - call mld_caggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,p,info) - if(info /= psb_success_) then - call psb_errpush(psb_err_from_subroutine_,name,a_err='mld_aggrmat_smth_asb') + call mld_caggrmat_biz_asb(a,desc_a,ilaggr,nlaggr, & + & p%parms,ac,op_prol,op_restr,info) + + case(mld_min_energy_) + + call mld_caggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr, & + & p%parms,ac,op_prol,op_restr,info) + + case default + info = psb_err_internal_error_ + call psb_errpush(info,name,a_err='Invalid aggr kind') + goto 9999 + + end select + if (info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,a_err='Inner aggrmat asb') + goto 9999 + end if + + + + ntaggr = sum(nlaggr) + + select case(p%parms%coarse_mat) + + case(mld_distr_mat_) + + nzl = ac%get_nzeros() + call ac%mv_to(bcoo) + + 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) - case(mld_min_energy_) + 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() - call mld_caggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info) - if(info /= psb_success_) then - call psb_errpush(psb_err_from_subroutine_,name,a_err='mld_aggrmat_smth_asb') + if (info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,a_err='psb_sp_free') goto 9999 end if - case default + 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()) - call psb_errpush(psb_err_internal_error_,name,a_err='Invalid aggr kind') - goto 9999 + 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 + call op_restr%set_nrows(p%desc_ac%get_local_cols()) + + 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_) & + & 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) return diff --git a/mlprec/impl/mld_caggrmat_biz_asb.f90 b/mlprec/impl/mld_caggrmat_biz_asb.f90 index 97d1ee4a..d278e4bf 100644 --- a/mlprec/impl/mld_caggrmat_biz_asb.f90 +++ b/mlprec/impl/mld_caggrmat_biz_asb.f90 @@ -78,7 +78,7 @@ ! info - integer, output. ! Error code. ! -subroutine mld_caggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,p,info) +subroutine mld_caggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_restr,info) use psb_base_mod use mld_c_inner_mod, mld_protect_name => mld_caggrmat_biz_asb @@ -88,9 +88,9 @@ subroutine mld_caggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,p,info) type(psb_cspmat_type), intent(in) :: a type(psb_desc_type), intent(in) :: desc_a integer, intent(inout) :: ilaggr(:), nlaggr(:) - type(mld_c_onelev_type), intent(inout), target :: p + type(mld_sml_parms), intent(inout) :: parms + type(psb_cspmat_type), intent(out) :: ac,op_prol,op_restr integer, intent(out) :: info - type(psb_cspmat_type) :: op_prol, op_restr, b ! Local variables integer :: nrow, nglob, ncol, ntaggr, ip, ndx,& @@ -98,7 +98,7 @@ subroutine mld_caggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,p,info) integer ::ictxt, np, me, err_act character(len=20) :: name type(psb_cspmat_type) :: am3, am4 - type(psb_c_coo_sparse_mat) :: acoo, acoof, bcoo + type(psb_c_coo_sparse_mat) :: tmpcoo type(psb_c_csr_sparse_mat) :: acsr1, acsr2, acsr3, acsrf, ptilde complex(psb_spk_), allocatable :: adiag(:) integer(psb_ipk_) :: ierr(5) @@ -123,7 +123,7 @@ subroutine mld_caggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,p,info) nrow = desc_a%get_local_rows() ncol = desc_a%get_local_cols() - theta = p%parms%aggr_thresh + theta = parms%aggr_thresh naggr = nlaggr(me+1) ntaggr = sum(nlaggr) @@ -131,7 +131,7 @@ subroutine mld_caggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,p,info) naggrm1 = sum(nlaggr(1:me)) naggrp1 = sum(nlaggr(1:me+1)) - filter_mat = (p%parms%aggr_filter == mld_filter_mat_) + filter_mat = (parms%aggr_filter == mld_filter_mat_) ilaggr(1:nrow) = ilaggr(1:nrow) + naggrm1 call psb_halo(ilaggr,desc_a,info) @@ -163,16 +163,16 @@ subroutine mld_caggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,p,info) end if ! 1. Allocate Ptilde in sparse matrix form - call acoo%allocate(ncol,naggr,ncol) + call tmpcoo%allocate(ncol,naggr,ncol) do i=1,nrow - acoo%val(i) = cone - acoo%ia(i) = i - acoo%ja(i) = ilaggr(i) + tmpcoo%val(i) = cone + tmpcoo%ia(i) = i + tmpcoo%ja(i) = ilaggr(i) end do - call acoo%set_nzeros(nrow) - call acoo%set_dupl(psb_dupl_add_) + call tmpcoo%set_nzeros(nrow) + call tmpcoo%set_dupl(psb_dupl_add_) - call ptilde%mv_from_coo(acoo,info) + call ptilde%mv_from_coo(tmpcoo,info) if (info == psb_success_) call a%cscnv(acsr3,info,dupl=psb_dupl_add_) if (debug_level >= psb_debug_outer_) & @@ -203,19 +203,19 @@ subroutine mld_caggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,p,info) end if enddo ! Take out zeroed terms - call acsrf%mv_to_coo(acoof,info) + call acsrf%mv_to_coo(tmpcoo,info) k = 0 - do j=1,acoof%get_nzeros() - if ((acoof%val(j) /= czero) .or. (acoof%ia(j) == acoof%ja(j))) then + do j=1,tmpcoo%get_nzeros() + if ((tmpcoo%val(j) /= czero) .or. (tmpcoo%ia(j) == tmpcoo%ja(j))) then k = k + 1 - acoof%val(k) = acoof%val(j) - acoof%ia(k) = acoof%ia(j) - acoof%ja(k) = acoof%ja(j) + tmpcoo%val(k) = tmpcoo%val(j) + tmpcoo%ia(k) = tmpcoo%ia(j) + tmpcoo%ja(k) = tmpcoo%ja(j) end if end do - call acoof%set_nzeros(k) - call acoof%set_dupl(psb_dupl_add_) - call acsrf%mv_from_coo(acoof,info) + call tmpcoo%set_nzeros(k) + call tmpcoo%set_dupl(psb_dupl_add_) + call acsrf%mv_from_coo(tmpcoo,info) end if @@ -232,9 +232,9 @@ subroutine mld_caggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,p,info) if (info /= psb_success_) goto 9999 - if (p%parms%aggr_omega_alg == mld_eig_est_) then + if (parms%aggr_omega_alg == mld_eig_est_) then - if (p%parms%aggr_eig == mld_max_norm_) then + if (parms%aggr_eig == mld_max_norm_) then ! ! This only works with CSR @@ -261,7 +261,7 @@ subroutine mld_caggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,p,info) goto 9999 end if omega = 4.d0/(3.d0*anorm) - p%parms%aggr_omega_val = omega + parms%aggr_omega_val = omega else info = psb_err_internal_error_ @@ -269,11 +269,11 @@ subroutine mld_caggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,p,info) goto 9999 end if - else if (p%parms%aggr_omega_alg == mld_user_choice_) then + else if (parms%aggr_omega_alg == mld_user_choice_) then - omega = p%parms%aggr_omega_val + omega = parms%aggr_omega_val - else if (p%parms%aggr_omega_alg /= mld_user_choice_) then + else if (parms%aggr_omega_alg /= mld_user_choice_) then info = psb_err_internal_error_ call psb_errpush(info,name,a_err='invalid mld_aggr_omega_alg_') goto 9999 @@ -372,7 +372,7 @@ subroutine mld_caggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,p,info) call psb_numbmm(a,op_prol,am3) if (debug_level >= psb_debug_outer_) & & write(debug_unit,*) me,' ',trim(name),& - & 'Done NUMBMM 2',p%parms%aggr_kind, mld_smooth_prol_ + & 'Done NUMBMM 2',parms%aggr_kind, mld_smooth_prol_ call op_prol%transp(op_restr) if (debug_level >= psb_debug_outer_) & @@ -389,10 +389,10 @@ subroutine mld_caggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,p,info) if (debug_level >= psb_debug_outer_) & & write(debug_unit,*) me,' ',trim(name),& & 'starting symbmm 3' - call psb_symbmm(op_restr,am3,b,info) - if (info == psb_success_) call psb_numbmm(op_restr,am3,b) + call psb_symbmm(op_restr,am3,ac,info) + if (info == psb_success_) call psb_numbmm(op_restr,am3,ac) if (info == psb_success_) call am3%free() - if (info == psb_success_) call b%cscnv(info,type='coo',dupl=psb_dupl_add_) + if (info == psb_success_) call ac%cscnv(info,type='coo',dupl=psb_dupl_add_) if (info /= psb_success_) then call psb_errpush(psb_err_internal_error_,name,a_err='Build b = op_restr x am3') goto 9999 @@ -401,107 +401,6 @@ subroutine mld_caggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,p,info) - select case(p%parms%coarse_mat) - - case(mld_distr_mat_) - - nzl = b%get_nzeros() - call b%mv_to(bcoo) - - 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 - call op_restr%set_nrows(p%desc_ac%get_local_cols()) - - 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_) & - & call psb_gather(p%ac,b,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 - if (debug_level >= psb_debug_outer_) & & write(debug_unit,*) me,' ',trim(name),& diff --git a/mlprec/impl/mld_caggrmat_minnrg_asb.f90 b/mlprec/impl/mld_caggrmat_minnrg_asb.f90 index 0573e6f5..6344b4a7 100644 --- a/mlprec/impl/mld_caggrmat_minnrg_asb.f90 +++ b/mlprec/impl/mld_caggrmat_minnrg_asb.f90 @@ -98,7 +98,7 @@ ! info - integer, output. ! Error code. ! -subroutine mld_caggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info) +subroutine mld_caggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_restr,info) use psb_base_mod use mld_c_inner_mod, mld_protect_name => mld_caggrmat_minnrg_asb @@ -108,9 +108,9 @@ subroutine mld_caggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info) type(psb_cspmat_type), intent(in) :: a type(psb_desc_type), intent(in) :: desc_a integer, intent(inout) :: ilaggr(:), nlaggr(:) - type(mld_c_onelev_type), intent(inout), target :: p + type(mld_sml_parms), intent(inout) :: parms + type(psb_cspmat_type), intent(out) :: ac,op_prol,op_restr integer, intent(out) :: info - type(psb_cspmat_type) :: op_prol,op_restr, b ! Local variables integer(psb_mpik_), allocatable :: nzbr(:), idisp(:) @@ -121,8 +121,8 @@ subroutine mld_caggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info) type(psb_cspmat_type) :: af, ptilde, rtilde, atran, atp, atdatp type(psb_cspmat_type) :: am3,am4, ap, adap,atmp,rada, ra, atmp2, dap, dadap, da type(psb_cspmat_type) :: dat, datp, datdatp, atmp3 - type(psb_c_coo_sparse_mat) :: acoo, acoof, bcoo, tmpcoo - type(psb_c_csr_sparse_mat) :: acsr1, acsr2, acsr3, bcsr, acsr, acsrf + type(psb_c_coo_sparse_mat) :: tmpcoo + type(psb_c_csr_sparse_mat) :: acsr1, acsr2, acsr3, acsr, acsrf type(psb_c_csc_sparse_mat) :: csc_dap, csc_dadap, csc_datp, csc_datdatp, acsc complex(psb_spk_), allocatable :: adiag(:), adinv(:) complex(psb_spk_), allocatable :: omf(:), omp(:), omi(:), oden(:) @@ -150,7 +150,7 @@ subroutine mld_caggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info) nrow = desc_a%get_local_rows() ncol = desc_a%get_local_cols() - theta = p%parms%aggr_thresh + theta = parms%aggr_thresh naggr = nlaggr(me+1) ntaggr = sum(nlaggr) @@ -165,7 +165,7 @@ subroutine mld_caggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info) naggrm1 = sum(nlaggr(1:me)) naggrp1 = sum(nlaggr(1:me+1)) - filter_mat = (p%parms%aggr_filter == mld_filter_mat_) + filter_mat = (parms%aggr_filter == mld_filter_mat_) ilaggr(1:nrow) = ilaggr(1:nrow) + naggrm1 call psb_halo(ilaggr,desc_a,info) @@ -207,16 +207,16 @@ subroutine mld_caggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info) ! 1. Allocate Ptilde in sparse matrix form - call acoo%allocate(ncol,ntaggr,ncol) + call tmpcoo%allocate(ncol,ntaggr,ncol) do i=1,ncol - acoo%val(i) = cone - acoo%ia(i) = i - acoo%ja(i) = ilaggr(i) + tmpcoo%val(i) = cone + tmpcoo%ia(i) = i + tmpcoo%ja(i) = ilaggr(i) end do - call acoo%set_nzeros(ncol) - call acoo%set_dupl(psb_dupl_add_) - call acoo%set_asb() - call ptilde%mv_from(acoo) + call tmpcoo%set_nzeros(ncol) + call tmpcoo%set_dupl(psb_dupl_add_) + call tmpcoo%set_asb() + call ptilde%mv_from(tmpcoo) call ptilde%cscnv(info,type='csr') !!$ call local_dump(me,ptilde,'csr-ptilde','Ptilde-1') @@ -570,164 +570,18 @@ subroutine mld_caggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info) & write(debug_unit,*) me,' ',trim(name),& & 'Done sphalo/ rwxtd' - call psb_symbmm(op_restr,am3,b,info) - if (info == psb_success_) call psb_numbmm(op_restr,am3,b) + call psb_symbmm(op_restr,am3,ac,info) + if (info == psb_success_) call psb_numbmm(op_restr,am3,ac) if (info == psb_success_) call am3%free() + if (info == psb_success_) call ac%cscnv(info,type='coo',dupl=psb_dupl_add_) if (info /= psb_success_) then call psb_errpush(psb_err_internal_error_,name,& - &a_err='Build b = op_restr x am3') + &a_err='Build ac = op_restr x am3') goto 9999 end if - if (debug_level >= psb_debug_outer_) & - & write(debug_unit,*) me,' ',trim(name),& - & 'Done mv_to_coo' - - - select case(p%parms%coarse_mat) - - case(mld_distr_mat_) - - call b%mv_to(bcoo) - nzl = bcoo%get_nzeros() - - if (debug_level >= psb_debug_outer_) & - & write(debug_unit,*) me,' ',trim(name),& - & ' B matrix nzl: ',nzl - - 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 bcoo%set_nrows(p%desc_ac%get_local_rows()) - call bcoo%set_ncols(p%desc_ac%get_local_cols()) - call bcoo%fix(info) - call p%ac%mv_from(bcoo) - call p%ac%set_asb() - - call p%ac%cscnv(info,type='csr') - - if (np>1) then - call op_prol%mv_to(acsr) - nzl = acsr%get_nzeros() - call psb_glob_to_loc(acsr%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(acsr) - endif - call op_prol%set_ncols(p%desc_ac%get_local_cols()) - - if (np>1) then - call op_restr%mv_to(tmpcoo) - nzl = tmpcoo%get_nzeros() - if (info == psb_success_) call psb_glob_to_loc(tmpcoo%ia(1:nzl),p%desc_ac,info,'I') - if(info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,a_err='Converting op_restr to local') - goto 9999 - end if - call op_restr%mv_from(tmpcoo) - call op_restr%cscnv(info,type='csr') - end if - call op_restr%set_nrows(p%desc_ac%get_local_cols()) - - 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_) then - call psb_errpush(psb_err_from_subroutine_,name,a_err='psb_cdall') - goto 9999 - end if - call psb_gather(p%ac,b,p%desc_ac,info,dupl=psb_dupl_add_,keeploc=.false.) - if(info /= psb_success_) goto 9999 -!!$ call psb_cdall(ictxt,p%desc_ac,info,mg=ntaggr,repl=.true.) -!!$ nzbr(:) = 0 -!!$ nzbr(me+1) = bcoo%get_nzeros() -!!$ -!!$ call psb_sum(ictxt,nzbr(1:np)) -!!$ nzac = sum(nzbr) -!!$ if (info == psb_success_) call tmpcoo%allocate(ntaggr,ntaggr,nzac) -!!$ if (info /= psb_success_) goto 9999 -!!$ -!!$ do ip=1,np -!!$ idispip) = sum(nzbr(1:ip-1)) -!!$ enddo -!!$ ndx = nzbr(me+1) -!!$ -!!$ call mpi_allgatherv(bcoo%val,ndx,mpi_complex,tmpcoo%val,nzbr,idisp,& -!!$ & mpi_complex,icomm,info) -!!$ if (info == psb_success_)& -!!$ & call mpi_allgatherv(bcoo%ia,ndx,psb_mpi_ipk_integer,tmpcoo%ia,nzbr,idisp,& -!!$ & psb_mpi_ipk_integer,icomm,info) -!!$ if (info == psb_success_)& -!!$ & call mpi_allgatherv(bcoo%ja,ndx,psb_mpi_ipk_integer,tmpcoo%ja,nzbr,idisp,& -!!$ & psb_mpi_ipk_integer,icomm,info) -!!$ -!!$ if (info /= psb_success_) then -!!$ call psb_errpush(psb_err_internal_error_,name,& -!!$ & a_err=' from mpi_allgatherv') -!!$ goto 9999 -!!$ end if -!!$ -!!$ call bcoo%free() -!!$ call tmpcoo%fix(info) -!!$ call p%ac%mv_from(tmpcoo) -!!$ call p%ac%cscnv(info,type='csr') -!!$ if(info /= psb_success_) goto 9999 -!!$ -!!$ deallocate(nzbr,idisp,stat=info) -!!$ if (info /= psb_success_) then -!!$ info = psb_err_alloc_dealloc_ -!!$ call psb_errpush(info,name) -!!$ goto 9999 -!!$ end if - 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') - 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 => R i.e. restriction operator - ! op_prol => P 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 if (debug_level >= psb_debug_outer_) & & write(debug_unit,*) me,' ',trim(name),& diff --git a/mlprec/impl/mld_caggrmat_nosmth_asb.f90 b/mlprec/impl/mld_caggrmat_nosmth_asb.f90 index 4f4ad6af..02294e6a 100644 --- a/mlprec/impl/mld_caggrmat_nosmth_asb.f90 +++ b/mlprec/impl/mld_caggrmat_nosmth_asb.f90 @@ -81,7 +81,7 @@ ! info - integer, output. ! Error code. ! -subroutine mld_caggrmat_nosmth_asb(a,desc_a,ilaggr,nlaggr,p,info) +subroutine mld_caggrmat_nosmth_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_restr,info) use psb_base_mod use mld_c_inner_mod, mld_protect_name => mld_caggrmat_nosmth_asb @@ -91,16 +91,16 @@ subroutine mld_caggrmat_nosmth_asb(a,desc_a,ilaggr,nlaggr,p,info) type(psb_cspmat_type), intent(in) :: a type(psb_desc_type), intent(in) :: desc_a integer, intent(inout) :: ilaggr(:), nlaggr(:) - type(mld_c_onelev_type), intent(inout), target :: p + type(mld_sml_parms), intent(inout) :: parms + type(psb_cspmat_type), intent(out) :: ac,op_prol,op_restr integer, intent(out) :: info - type(psb_cspmat_type) :: b, op_prol,op_restr ! Local variables integer :: ictxt,np,me, err_act integer(psb_mpik_) :: icomm, ndx, minfo character(len=20) :: name integer(psb_ipk_) :: ierr(5) - type(psb_c_coo_sparse_mat) :: acoo1, acoo2, bcoo, ac_coo, acoo + type(psb_c_coo_sparse_mat) :: ac_coo, acoo type(psb_c_csr_sparse_mat) :: acsr1, acsr2 integer :: debug_level, debug_unit integer :: nrow, nglob, ncol, ntaggr, nzl, ip, & @@ -134,136 +134,36 @@ subroutine mld_caggrmat_nosmth_asb(a,desc_a,ilaggr,nlaggr,p,info) goto 9999 end if - call acoo1%allocate(ncol,ntaggr,ncol) + call acoo%allocate(ncol,ntaggr,ncol) do i=1,nrow - acoo1%val(i) = cone - acoo1%ia(i) = i - acoo1%ja(i) = ilaggr(i) + acoo%val(i) = cone + acoo%ia(i) = i + acoo%ja(i) = ilaggr(i) end do - call acoo1%set_dupl(psb_dupl_add_) - call acoo1%set_nzeros(nrow) - call acoo1%set_asb() - call acoo1%fix(info) + call acoo%set_dupl(psb_dupl_add_) + call acoo%set_nzeros(nrow) + call acoo%set_asb() + call acoo%fix(info) - call op_prol%mv_from(acoo1) + call op_prol%mv_from(acoo) call op_prol%cscnv(info,type='csr',dupl=psb_dupl_add_) if (info == psb_success_) call op_prol%transp(op_restr) - call a%csclip(bcoo,info,jmax=nrow) + call a%csclip(ac_coo,info,jmax=nrow) - nzt = bcoo%get_nzeros() + nzt = ac_coo%get_nzeros() do i=1, nzt - bcoo%ia(i) = ilaggr(bcoo%ia(i)) - bcoo%ja(i) = ilaggr(bcoo%ja(i)) + ac_coo%ia(i) = ilaggr(ac_coo%ia(i)) + ac_coo%ja(i) = ilaggr(ac_coo%ja(i)) enddo - call bcoo%set_nrows(naggr) - call bcoo%set_ncols(naggr) - call bcoo%set_dupl(psb_dupl_add_) - call bcoo%fix(info) - - - call b%mv_from(bcoo) - - if (p%parms%coarse_mat == mld_repl_mat_) then - - 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_) then - call psb_errpush(psb_err_from_subroutine_,name,a_err='psb_cdall') - goto 9999 - end if - call psb_gather(p%ac,b,p%desc_ac,info,dupl=psb_dupl_add_,keeploc=.false.) - if(info /= psb_success_) goto 9999 - - else if (p%parms%coarse_mat == mld_distr_mat_) then - - nzl = b%get_nzeros() - call b%mv_to(bcoo) - - 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 - call op_restr%set_nrows(p%desc_ac%get_local_cols()) - - if (debug_level >= psb_debug_outer_) & - & write(debug_unit,*) me,' ',trim(name),& - & 'Done ac ' - else - info = psb_err_internal_error_ - call psb_errpush(psb_err_internal_error_,name,a_err='invalid mld_coarse_mat_') - goto 9999 - end if - - - 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='cscnv') - 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 - ! - if (info == psb_success_) & - & 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='linmap build') - goto 9999 - end if + call ac_coo%set_nrows(naggr) + call ac_coo%set_ncols(naggr) + call ac_coo%set_dupl(psb_dupl_add_) + call ac_coo%fix(info) + call ac%mv_from(ac_coo) call psb_erractionrestore(err_act) diff --git a/mlprec/impl/mld_caggrmat_smth_asb.f90 b/mlprec/impl/mld_caggrmat_smth_asb.f90 index ede23f6d..d29321b7 100644 --- a/mlprec/impl/mld_caggrmat_smth_asb.f90 +++ b/mlprec/impl/mld_caggrmat_smth_asb.f90 @@ -93,7 +93,7 @@ ! info - integer, output. ! Error code. ! -subroutine mld_caggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info) +subroutine mld_caggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_restr,info) use psb_base_mod use mld_c_inner_mod, mld_protect_name => mld_caggrmat_smth_asb @@ -103,9 +103,9 @@ subroutine mld_caggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info) type(psb_cspmat_type), intent(in) :: a type(psb_desc_type), intent(in) :: desc_a integer, intent(inout) :: ilaggr(:), nlaggr(:) - type(mld_c_onelev_type), intent(inout), target :: p + type(mld_sml_parms), intent(inout) :: parms + type(psb_cspmat_type), intent(out) :: ac,op_prol,op_restr integer, intent(out) :: info - type(psb_cspmat_type) :: op_prol, op_restr, b ! Local variables integer :: nrow, nglob, ncol, ntaggr, ip, ndx,& @@ -113,7 +113,7 @@ subroutine mld_caggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info) integer ::ictxt, np, me, err_act character(len=20) :: name type(psb_cspmat_type) :: am3, am4 - type(psb_c_coo_sparse_mat) :: acoo, acoof, bcoo + type(psb_c_coo_sparse_mat) :: tmpcoo type(psb_c_csr_sparse_mat) :: acsr1, acsr2, acsr3, acsrf, ptilde complex(psb_spk_), allocatable :: adiag(:) integer(psb_ipk_) :: ierr(5) @@ -138,14 +138,14 @@ subroutine mld_caggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info) nrow = desc_a%get_local_rows() ncol = desc_a%get_local_cols() - theta = p%parms%aggr_thresh + theta = parms%aggr_thresh naggr = nlaggr(me+1) ntaggr = sum(nlaggr) naggrm1 = sum(nlaggr(1:me)) naggrp1 = sum(nlaggr(1:me+1)) - filter_mat = (p%parms%aggr_filter == mld_filter_mat_) + filter_mat = (parms%aggr_filter == mld_filter_mat_) ilaggr(1:nrow) = ilaggr(1:nrow) + naggrm1 call psb_halo(ilaggr,desc_a,info) @@ -177,16 +177,16 @@ subroutine mld_caggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info) end if ! 1. Allocate Ptilde in sparse matrix form - call acoo%allocate(ncol,ntaggr,ncol) + call tmpcoo%allocate(ncol,ntaggr,ncol) do i=1,ncol - acoo%val(i) = cone - acoo%ia(i) = i - acoo%ja(i) = ilaggr(i) + tmpcoo%val(i) = cone + tmpcoo%ia(i) = i + tmpcoo%ja(i) = ilaggr(i) end do - call acoo%set_nzeros(ncol) - call acoo%set_dupl(psb_dupl_add_) + call tmpcoo%set_nzeros(ncol) + call tmpcoo%set_dupl(psb_dupl_add_) - call ptilde%mv_from_coo(acoo,info) + call ptilde%mv_from_coo(tmpcoo,info) if (info == psb_success_) call a%cscnv(acsr3,info,dupl=psb_dupl_add_) if (debug_level >= psb_debug_outer_) & @@ -217,19 +217,19 @@ subroutine mld_caggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info) end if enddo ! Take out zeroed terms - call acsrf%mv_to_coo(acoof,info) + call acsrf%mv_to_coo(tmpcoo,info) k = 0 - do j=1,acoof%get_nzeros() - if ((acoof%val(j) /= czero) .or. (acoof%ia(j) == acoof%ja(j))) then + do j=1,tmpcoo%get_nzeros() + if ((tmpcoo%val(j) /= czero) .or. (tmpcoo%ia(j) == tmpcoo%ja(j))) then k = k + 1 - acoof%val(k) = acoof%val(j) - acoof%ia(k) = acoof%ia(j) - acoof%ja(k) = acoof%ja(j) + tmpcoo%val(k) = tmpcoo%val(j) + tmpcoo%ia(k) = tmpcoo%ia(j) + tmpcoo%ja(k) = tmpcoo%ja(j) end if end do - call acoof%set_nzeros(k) - call acoof%set_dupl(psb_dupl_add_) - call acsrf%mv_from_coo(acoof,info) + call tmpcoo%set_nzeros(k) + call tmpcoo%set_dupl(psb_dupl_add_) + call acsrf%mv_from_coo(tmpcoo,info) end if @@ -246,13 +246,13 @@ subroutine mld_caggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info) if (info /= psb_success_) goto 9999 - if (p%parms%aggr_omega_alg == mld_eig_est_) then + if (parms%aggr_omega_alg == mld_eig_est_) then - if (p%parms%aggr_eig == mld_max_norm_) then + if (parms%aggr_eig == mld_max_norm_) then anorm = acsr3%csnmi() omega = 4.d0/(3.d0*anorm) - p%parms%aggr_omega_val = omega + parms%aggr_omega_val = omega else info = psb_err_internal_error_ @@ -260,11 +260,11 @@ subroutine mld_caggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info) goto 9999 end if - else if (p%parms%aggr_omega_alg == mld_user_choice_) then + else if (parms%aggr_omega_alg == mld_user_choice_) then - omega = p%parms%aggr_omega_val + omega = parms%aggr_omega_val - else if (p%parms%aggr_omega_alg /= mld_user_choice_) then + else if (parms%aggr_omega_alg /= mld_user_choice_) then info = psb_err_internal_error_ call psb_errpush(info,name,a_err='invalid mld_aggr_omega_alg_') goto 9999 @@ -369,27 +369,27 @@ subroutine mld_caggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info) call psb_numbmm(a,op_prol,am3) if (debug_level >= psb_debug_outer_) & & write(debug_unit,*) me,' ',trim(name),& - & 'Done NUMBMM 2',p%parms%aggr_kind, mld_smooth_prol_ + & 'Done NUMBMM 2',parms%aggr_kind, mld_smooth_prol_ call op_prol%transp(op_restr) - call op_restr%mv_to(acoo) - nzl = acoo%get_nzeros() + call op_restr%mv_to(tmpcoo) + nzl = tmpcoo%get_nzeros() i=0 ! ! Now we have to fix this. The only rows of B that are correct ! are those corresponding to "local" aggregates, i.e. indices in ilaggr(:) ! do k=1, nzl - if ((naggrm1 < acoo%ia(k)) .and.(acoo%ia(k) <= naggrp1)) then + if ((naggrm1 < tmpcoo%ia(k)) .and.(tmpcoo%ia(k) <= naggrp1)) then i = i+1 - acoo%val(i) = acoo%val(k) - acoo%ia(i) = acoo%ia(k) - acoo%ja(i) = acoo%ja(k) + tmpcoo%val(i) = tmpcoo%val(k) + tmpcoo%ia(i) = tmpcoo%ia(k) + tmpcoo%ja(i) = tmpcoo%ja(k) end if end do - call acoo%set_nzeros(i) - call acoo%trim() - call op_restr%mv_from(acoo) + call tmpcoo%set_nzeros(i) + call tmpcoo%trim() + call op_restr%mv_from(tmpcoo) call op_restr%cscnv(info,type='csr',dupl=psb_dupl_add_) if (info /= psb_success_) then call psb_errpush(psb_err_from_subroutine_,name,a_err='spcnv op_restr') @@ -413,113 +413,12 @@ subroutine mld_caggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info) if (debug_level >= psb_debug_outer_) & & write(debug_unit,*) me,' ',trim(name),& & 'starting symbmm 3' - call psb_symbmm(op_restr,am3,b,info) - if (info == psb_success_) call psb_numbmm(op_restr,am3,b) + call psb_symbmm(op_restr,am3,ac,info) + if (info == psb_success_) call psb_numbmm(op_restr,am3,ac) if (info == psb_success_) call am3%free() - if (info == psb_success_) call b%cscnv(info,type='coo',dupl=psb_dupl_add_) + if (info == psb_success_) call ac%cscnv(info,type='coo',dupl=psb_dupl_add_) if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,a_err='Build b = op_restr x am3') - goto 9999 - end if - - - - select case(p%parms%coarse_mat) - - case(mld_distr_mat_) - - nzl = b%get_nzeros() - call b%mv_to(bcoo) - - 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 - call op_restr%set_nrows(p%desc_ac%get_local_cols()) - - 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_) & - & call psb_gather(p%ac,b,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') + call psb_errpush(psb_err_internal_error_,name,a_err='Build ac = op_restr x am3') goto 9999 end if diff --git a/mlprec/impl/mld_daggrmat_asb.f90 b/mlprec/impl/mld_daggrmat_asb.f90 index 7024b969..6dbb942b 100644 --- a/mlprec/impl/mld_daggrmat_asb.f90 +++ b/mlprec/impl/mld_daggrmat_asb.f90 @@ -113,7 +113,11 @@ subroutine mld_daggrmat_asb(a,desc_a,ilaggr,nlaggr,p,info) integer, intent(out) :: info ! Local variables - type(psb_dspmat_type) :: b, op_prol,op_restr + 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 :: nzl,ntaggr + integer :: debug_level, debug_unit integer :: ictxt,np,me, err_act character(len=20) :: name @@ -121,6 +125,9 @@ subroutine mld_daggrmat_asb(a,desc_a,ilaggr,nlaggr,p,info) if(psb_get_errstatus().ne.0) return info=psb_success_ call psb_erractionsave(err_act) + debug_unit = psb_get_debug_unit() + debug_level = psb_get_debug_level() + ictxt = desc_a%get_context() @@ -129,43 +136,139 @@ subroutine mld_daggrmat_asb(a,desc_a,ilaggr,nlaggr,p,info) select case (p%parms%aggr_kind) case (mld_no_smooth_) - call mld_daggrmat_nosmth_asb(a,desc_a,ilaggr,nlaggr,p,info) - if(info /= psb_success_) then - call psb_errpush(psb_err_from_subroutine_,name,a_err='mld_aggrmat_nosmth_asb') - goto 9999 - end if + call mld_daggrmat_nosmth_asb(a,desc_a,ilaggr,nlaggr,& + & p%parms,ac,op_prol,op_restr,info) case(mld_smooth_prol_) - call mld_daggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info) - if(info /= psb_success_) then - call psb_errpush(psb_err_from_subroutine_,name,a_err='mld_aggrmat_smth_asb') - goto 9999 - end if + call mld_daggrmat_smth_asb(a,desc_a,ilaggr,nlaggr, & + & p%parms,ac,op_prol,op_restr,info) case(mld_biz_prol_) - call mld_daggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,p,info) - if(info /= psb_success_) then - call psb_errpush(psb_err_from_subroutine_,name,a_err='mld_aggrmat_smth_asb') + call mld_daggrmat_biz_asb(a,desc_a,ilaggr,nlaggr, & + & p%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) + + case default + info = psb_err_internal_error_ + call psb_errpush(info,name,a_err='Invalid aggr kind') + goto 9999 + + end select + if (info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,a_err='Inner aggrmat asb') + goto 9999 + end if + + + + ntaggr = sum(nlaggr) + + select case(p%parms%coarse_mat) + + case(mld_distr_mat_) + + nzl = ac%get_nzeros() + call ac%mv_to(bcoo) + + 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) - case(mld_min_energy_) + 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() - call mld_daggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info) - if(info /= psb_success_) then - call psb_errpush(psb_err_from_subroutine_,name,a_err='mld_aggrmat_smth_asb') + if (info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,a_err='psb_sp_free') goto 9999 end if - case default + 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()) - call psb_errpush(psb_err_internal_error_,name,a_err='Invalid aggr kind') - goto 9999 + 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 + call op_restr%set_nrows(p%desc_ac%get_local_cols()) + + 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_) & + & 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) return diff --git a/mlprec/impl/mld_daggrmat_biz_asb.f90 b/mlprec/impl/mld_daggrmat_biz_asb.f90 index 91e0178d..21ca67b6 100644 --- a/mlprec/impl/mld_daggrmat_biz_asb.f90 +++ b/mlprec/impl/mld_daggrmat_biz_asb.f90 @@ -78,7 +78,7 @@ ! info - integer, output. ! Error code. ! -subroutine mld_daggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,p,info) +subroutine mld_daggrmat_biz_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_biz_asb @@ -88,9 +88,9 @@ subroutine mld_daggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,p,info) type(psb_dspmat_type), intent(in) :: a type(psb_desc_type), intent(in) :: desc_a integer, intent(inout) :: ilaggr(:), nlaggr(:) - type(mld_d_onelev_type), intent(inout), target :: p + type(mld_dml_parms), intent(inout) :: parms + type(psb_dspmat_type), intent(out) :: ac,op_prol,op_restr integer, intent(out) :: info - type(psb_dspmat_type) :: op_prol, op_restr, b ! Local variables integer :: nrow, nglob, ncol, ntaggr, ip, ndx,& @@ -98,7 +98,7 @@ subroutine mld_daggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,p,info) integer ::ictxt, np, me, err_act character(len=20) :: name type(psb_dspmat_type) :: am3, am4 - type(psb_d_coo_sparse_mat) :: acoo, acoof, bcoo + type(psb_d_coo_sparse_mat) :: tmpcoo type(psb_d_csr_sparse_mat) :: acsr1, acsr2, acsr3, acsrf, ptilde real(psb_dpk_), allocatable :: adiag(:) integer(psb_ipk_) :: ierr(5) @@ -123,7 +123,7 @@ subroutine mld_daggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,p,info) nrow = desc_a%get_local_rows() ncol = desc_a%get_local_cols() - theta = p%parms%aggr_thresh + theta = parms%aggr_thresh naggr = nlaggr(me+1) ntaggr = sum(nlaggr) @@ -131,7 +131,7 @@ subroutine mld_daggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,p,info) naggrm1 = sum(nlaggr(1:me)) naggrp1 = sum(nlaggr(1:me+1)) - filter_mat = (p%parms%aggr_filter == mld_filter_mat_) + filter_mat = (parms%aggr_filter == mld_filter_mat_) ilaggr(1:nrow) = ilaggr(1:nrow) + naggrm1 call psb_halo(ilaggr,desc_a,info) @@ -163,16 +163,16 @@ subroutine mld_daggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,p,info) end if ! 1. Allocate Ptilde in sparse matrix form - call acoo%allocate(ncol,naggr,ncol) + call tmpcoo%allocate(ncol,naggr,ncol) do i=1,nrow - acoo%val(i) = done - acoo%ia(i) = i - acoo%ja(i) = ilaggr(i) + tmpcoo%val(i) = done + tmpcoo%ia(i) = i + tmpcoo%ja(i) = ilaggr(i) end do - call acoo%set_nzeros(nrow) - call acoo%set_dupl(psb_dupl_add_) + call tmpcoo%set_nzeros(nrow) + call tmpcoo%set_dupl(psb_dupl_add_) - call ptilde%mv_from_coo(acoo,info) + call ptilde%mv_from_coo(tmpcoo,info) if (info == psb_success_) call a%cscnv(acsr3,info,dupl=psb_dupl_add_) if (debug_level >= psb_debug_outer_) & @@ -203,19 +203,19 @@ subroutine mld_daggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,p,info) end if enddo ! Take out zeroed terms - call acsrf%mv_to_coo(acoof,info) + call acsrf%mv_to_coo(tmpcoo,info) k = 0 - do j=1,acoof%get_nzeros() - if ((acoof%val(j) /= dzero) .or. (acoof%ia(j) == acoof%ja(j))) then + do j=1,tmpcoo%get_nzeros() + if ((tmpcoo%val(j) /= dzero) .or. (tmpcoo%ia(j) == tmpcoo%ja(j))) then k = k + 1 - acoof%val(k) = acoof%val(j) - acoof%ia(k) = acoof%ia(j) - acoof%ja(k) = acoof%ja(j) + tmpcoo%val(k) = tmpcoo%val(j) + tmpcoo%ia(k) = tmpcoo%ia(j) + tmpcoo%ja(k) = tmpcoo%ja(j) end if end do - call acoof%set_nzeros(k) - call acoof%set_dupl(psb_dupl_add_) - call acsrf%mv_from_coo(acoof,info) + call tmpcoo%set_nzeros(k) + call tmpcoo%set_dupl(psb_dupl_add_) + call acsrf%mv_from_coo(tmpcoo,info) end if @@ -232,9 +232,9 @@ subroutine mld_daggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,p,info) if (info /= psb_success_) goto 9999 - if (p%parms%aggr_omega_alg == mld_eig_est_) then + if (parms%aggr_omega_alg == mld_eig_est_) then - if (p%parms%aggr_eig == mld_max_norm_) then + if (parms%aggr_eig == mld_max_norm_) then ! ! This only works with CSR @@ -261,7 +261,7 @@ subroutine mld_daggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,p,info) goto 9999 end if omega = 4.d0/(3.d0*anorm) - p%parms%aggr_omega_val = omega + parms%aggr_omega_val = omega else info = psb_err_internal_error_ @@ -269,11 +269,11 @@ subroutine mld_daggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,p,info) goto 9999 end if - else if (p%parms%aggr_omega_alg == mld_user_choice_) then + else if (parms%aggr_omega_alg == mld_user_choice_) then - omega = p%parms%aggr_omega_val + omega = parms%aggr_omega_val - else if (p%parms%aggr_omega_alg /= mld_user_choice_) then + else if (parms%aggr_omega_alg /= mld_user_choice_) then info = psb_err_internal_error_ call psb_errpush(info,name,a_err='invalid mld_aggr_omega_alg_') goto 9999 @@ -372,7 +372,7 @@ subroutine mld_daggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,p,info) call psb_numbmm(a,op_prol,am3) if (debug_level >= psb_debug_outer_) & & write(debug_unit,*) me,' ',trim(name),& - & 'Done NUMBMM 2',p%parms%aggr_kind, mld_smooth_prol_ + & 'Done NUMBMM 2',parms%aggr_kind, mld_smooth_prol_ call op_prol%transp(op_restr) if (debug_level >= psb_debug_outer_) & @@ -389,10 +389,10 @@ subroutine mld_daggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,p,info) if (debug_level >= psb_debug_outer_) & & write(debug_unit,*) me,' ',trim(name),& & 'starting symbmm 3' - call psb_symbmm(op_restr,am3,b,info) - if (info == psb_success_) call psb_numbmm(op_restr,am3,b) + call psb_symbmm(op_restr,am3,ac,info) + if (info == psb_success_) call psb_numbmm(op_restr,am3,ac) if (info == psb_success_) call am3%free() - if (info == psb_success_) call b%cscnv(info,type='coo',dupl=psb_dupl_add_) + if (info == psb_success_) call ac%cscnv(info,type='coo',dupl=psb_dupl_add_) if (info /= psb_success_) then call psb_errpush(psb_err_internal_error_,name,a_err='Build b = op_restr x am3') goto 9999 @@ -401,107 +401,6 @@ subroutine mld_daggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,p,info) - select case(p%parms%coarse_mat) - - case(mld_distr_mat_) - - nzl = b%get_nzeros() - call b%mv_to(bcoo) - - 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 - call op_restr%set_nrows(p%desc_ac%get_local_cols()) - - 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_) & - & call psb_gather(p%ac,b,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 - if (debug_level >= psb_debug_outer_) & & write(debug_unit,*) me,' ',trim(name),& diff --git a/mlprec/impl/mld_daggrmat_minnrg_asb.f90 b/mlprec/impl/mld_daggrmat_minnrg_asb.f90 index 7e46a943..89e8e015 100644 --- a/mlprec/impl/mld_daggrmat_minnrg_asb.f90 +++ b/mlprec/impl/mld_daggrmat_minnrg_asb.f90 @@ -98,7 +98,7 @@ ! info - integer, output. ! Error code. ! -subroutine mld_daggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info) +subroutine mld_daggrmat_minnrg_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_minnrg_asb @@ -108,9 +108,9 @@ subroutine mld_daggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info) type(psb_dspmat_type), intent(in) :: a type(psb_desc_type), intent(in) :: desc_a integer, intent(inout) :: ilaggr(:), nlaggr(:) - type(mld_d_onelev_type), intent(inout), target :: p + type(mld_dml_parms), intent(inout) :: parms + type(psb_dspmat_type), intent(out) :: ac,op_prol,op_restr integer, intent(out) :: info - type(psb_dspmat_type) :: op_prol,op_restr, b ! Local variables integer(psb_mpik_), allocatable :: nzbr(:), idisp(:) @@ -121,8 +121,8 @@ subroutine mld_daggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info) type(psb_dspmat_type) :: af, ptilde, rtilde, atran, atp, atdatp type(psb_dspmat_type) :: am3,am4, ap, adap,atmp,rada, ra, atmp2, dap, dadap, da type(psb_dspmat_type) :: dat, datp, datdatp, atmp3 - type(psb_d_coo_sparse_mat) :: acoo, acoof, bcoo, tmpcoo - type(psb_d_csr_sparse_mat) :: acsr1, acsr2, acsr3, bcsr, acsr, acsrf + type(psb_d_coo_sparse_mat) :: tmpcoo + type(psb_d_csr_sparse_mat) :: acsr1, acsr2, acsr3, acsr, acsrf type(psb_d_csc_sparse_mat) :: csc_dap, csc_dadap, csc_datp, csc_datdatp, acsc real(psb_dpk_), allocatable :: adiag(:), adinv(:) real(psb_dpk_), allocatable :: omf(:), omp(:), omi(:), oden(:) @@ -150,7 +150,7 @@ subroutine mld_daggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info) nrow = desc_a%get_local_rows() ncol = desc_a%get_local_cols() - theta = p%parms%aggr_thresh + theta = parms%aggr_thresh naggr = nlaggr(me+1) ntaggr = sum(nlaggr) @@ -165,7 +165,7 @@ subroutine mld_daggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info) naggrm1 = sum(nlaggr(1:me)) naggrp1 = sum(nlaggr(1:me+1)) - filter_mat = (p%parms%aggr_filter == mld_filter_mat_) + filter_mat = (parms%aggr_filter == mld_filter_mat_) ilaggr(1:nrow) = ilaggr(1:nrow) + naggrm1 call psb_halo(ilaggr,desc_a,info) @@ -207,16 +207,16 @@ subroutine mld_daggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info) ! 1. Allocate Ptilde in sparse matrix form - call acoo%allocate(ncol,ntaggr,ncol) + call tmpcoo%allocate(ncol,ntaggr,ncol) do i=1,ncol - acoo%val(i) = done - acoo%ia(i) = i - acoo%ja(i) = ilaggr(i) + tmpcoo%val(i) = done + tmpcoo%ia(i) = i + tmpcoo%ja(i) = ilaggr(i) end do - call acoo%set_nzeros(ncol) - call acoo%set_dupl(psb_dupl_add_) - call acoo%set_asb() - call ptilde%mv_from(acoo) + call tmpcoo%set_nzeros(ncol) + call tmpcoo%set_dupl(psb_dupl_add_) + call tmpcoo%set_asb() + call ptilde%mv_from(tmpcoo) call ptilde%cscnv(info,type='csr') !!$ call local_dump(me,ptilde,'csr-ptilde','Ptilde-1') @@ -570,164 +570,18 @@ subroutine mld_daggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info) & write(debug_unit,*) me,' ',trim(name),& & 'Done sphalo/ rwxtd' - call psb_symbmm(op_restr,am3,b,info) - if (info == psb_success_) call psb_numbmm(op_restr,am3,b) + call psb_symbmm(op_restr,am3,ac,info) + if (info == psb_success_) call psb_numbmm(op_restr,am3,ac) if (info == psb_success_) call am3%free() + if (info == psb_success_) call ac%cscnv(info,type='coo',dupl=psb_dupl_add_) if (info /= psb_success_) then call psb_errpush(psb_err_internal_error_,name,& - &a_err='Build b = op_restr x am3') + &a_err='Build ac = op_restr x am3') goto 9999 end if - if (debug_level >= psb_debug_outer_) & - & write(debug_unit,*) me,' ',trim(name),& - & 'Done mv_to_coo' - - - select case(p%parms%coarse_mat) - - case(mld_distr_mat_) - - call b%mv_to(bcoo) - nzl = bcoo%get_nzeros() - - if (debug_level >= psb_debug_outer_) & - & write(debug_unit,*) me,' ',trim(name),& - & ' B matrix nzl: ',nzl - - 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 bcoo%set_nrows(p%desc_ac%get_local_rows()) - call bcoo%set_ncols(p%desc_ac%get_local_cols()) - call bcoo%fix(info) - call p%ac%mv_from(bcoo) - call p%ac%set_asb() - - call p%ac%cscnv(info,type='csr') - - if (np>1) then - call op_prol%mv_to(acsr) - nzl = acsr%get_nzeros() - call psb_glob_to_loc(acsr%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(acsr) - endif - call op_prol%set_ncols(p%desc_ac%get_local_cols()) - - if (np>1) then - call op_restr%mv_to(tmpcoo) - nzl = tmpcoo%get_nzeros() - if (info == psb_success_) call psb_glob_to_loc(tmpcoo%ia(1:nzl),p%desc_ac,info,'I') - if(info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,a_err='Converting op_restr to local') - goto 9999 - end if - call op_restr%mv_from(tmpcoo) - call op_restr%cscnv(info,type='csr') - end if - call op_restr%set_nrows(p%desc_ac%get_local_cols()) - - 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_) then - call psb_errpush(psb_err_from_subroutine_,name,a_err='psb_cdall') - goto 9999 - end if - call psb_gather(p%ac,b,p%desc_ac,info,dupl=psb_dupl_add_,keeploc=.false.) - if(info /= psb_success_) goto 9999 -!!$ call psb_cdall(ictxt,p%desc_ac,info,mg=ntaggr,repl=.true.) -!!$ nzbr(:) = 0 -!!$ nzbr(me+1) = bcoo%get_nzeros() -!!$ -!!$ call psb_sum(ictxt,nzbr(1:np)) -!!$ nzac = sum(nzbr) -!!$ if (info == psb_success_) call tmpcoo%allocate(ntaggr,ntaggr,nzac) -!!$ if (info /= psb_success_) goto 9999 -!!$ -!!$ do ip=1,np -!!$ idispip) = sum(nzbr(1:ip-1)) -!!$ enddo -!!$ ndx = nzbr(me+1) -!!$ -!!$ call mpi_allgatherv(bcoo%val,ndx,mpi_double_precision,tmpcoo%val,nzbr,idisp,& -!!$ & mpi_double_precision,icomm,info) -!!$ if (info == psb_success_)& -!!$ & call mpi_allgatherv(bcoo%ia,ndx,psb_mpi_ipk_integer,tmpcoo%ia,nzbr,idisp,& -!!$ & psb_mpi_ipk_integer,icomm,info) -!!$ if (info == psb_success_)& -!!$ & call mpi_allgatherv(bcoo%ja,ndx,psb_mpi_ipk_integer,tmpcoo%ja,nzbr,idisp,& -!!$ & psb_mpi_ipk_integer,icomm,info) -!!$ -!!$ if (info /= psb_success_) then -!!$ call psb_errpush(psb_err_internal_error_,name,& -!!$ & a_err=' from mpi_allgatherv') -!!$ goto 9999 -!!$ end if -!!$ -!!$ call bcoo%free() -!!$ call tmpcoo%fix(info) -!!$ call p%ac%mv_from(tmpcoo) -!!$ call p%ac%cscnv(info,type='csr') -!!$ if(info /= psb_success_) goto 9999 -!!$ -!!$ deallocate(nzbr,idisp,stat=info) -!!$ if (info /= psb_success_) then -!!$ info = psb_err_alloc_dealloc_ -!!$ call psb_errpush(info,name) -!!$ goto 9999 -!!$ end if - 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') - 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 => R i.e. restriction operator - ! op_prol => P 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 if (debug_level >= psb_debug_outer_) & & write(debug_unit,*) me,' ',trim(name),& diff --git a/mlprec/impl/mld_daggrmat_nosmth_asb.f90 b/mlprec/impl/mld_daggrmat_nosmth_asb.f90 index b054fc78..32f076f7 100644 --- a/mlprec/impl/mld_daggrmat_nosmth_asb.f90 +++ b/mlprec/impl/mld_daggrmat_nosmth_asb.f90 @@ -81,7 +81,7 @@ ! info - integer, output. ! Error code. ! -subroutine mld_daggrmat_nosmth_asb(a,desc_a,ilaggr,nlaggr,p,info) +subroutine mld_daggrmat_nosmth_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_nosmth_asb @@ -91,16 +91,16 @@ subroutine mld_daggrmat_nosmth_asb(a,desc_a,ilaggr,nlaggr,p,info) type(psb_dspmat_type), intent(in) :: a type(psb_desc_type), intent(in) :: desc_a integer, intent(inout) :: ilaggr(:), nlaggr(:) - type(mld_d_onelev_type), intent(inout), target :: p + type(mld_dml_parms), intent(inout) :: parms + type(psb_dspmat_type), intent(out) :: ac,op_prol,op_restr integer, intent(out) :: info - type(psb_dspmat_type) :: b, op_prol,op_restr ! Local variables integer :: ictxt,np,me, err_act integer(psb_mpik_) :: icomm, ndx, minfo character(len=20) :: name integer(psb_ipk_) :: ierr(5) - type(psb_d_coo_sparse_mat) :: acoo1, acoo2, bcoo, ac_coo, acoo + type(psb_d_coo_sparse_mat) :: ac_coo, acoo type(psb_d_csr_sparse_mat) :: acsr1, acsr2 integer :: debug_level, debug_unit integer :: nrow, nglob, ncol, ntaggr, nzl, ip, & @@ -134,136 +134,36 @@ subroutine mld_daggrmat_nosmth_asb(a,desc_a,ilaggr,nlaggr,p,info) goto 9999 end if - call acoo1%allocate(ncol,ntaggr,ncol) + call acoo%allocate(ncol,ntaggr,ncol) do i=1,nrow - acoo1%val(i) = done - acoo1%ia(i) = i - acoo1%ja(i) = ilaggr(i) + acoo%val(i) = done + acoo%ia(i) = i + acoo%ja(i) = ilaggr(i) end do - call acoo1%set_dupl(psb_dupl_add_) - call acoo1%set_nzeros(nrow) - call acoo1%set_asb() - call acoo1%fix(info) + call acoo%set_dupl(psb_dupl_add_) + call acoo%set_nzeros(nrow) + call acoo%set_asb() + call acoo%fix(info) - call op_prol%mv_from(acoo1) + call op_prol%mv_from(acoo) call op_prol%cscnv(info,type='csr',dupl=psb_dupl_add_) if (info == psb_success_) call op_prol%transp(op_restr) - call a%csclip(bcoo,info,jmax=nrow) + call a%csclip(ac_coo,info,jmax=nrow) - nzt = bcoo%get_nzeros() + nzt = ac_coo%get_nzeros() do i=1, nzt - bcoo%ia(i) = ilaggr(bcoo%ia(i)) - bcoo%ja(i) = ilaggr(bcoo%ja(i)) + ac_coo%ia(i) = ilaggr(ac_coo%ia(i)) + ac_coo%ja(i) = ilaggr(ac_coo%ja(i)) enddo - call bcoo%set_nrows(naggr) - call bcoo%set_ncols(naggr) - call bcoo%set_dupl(psb_dupl_add_) - call bcoo%fix(info) - - - call b%mv_from(bcoo) - - if (p%parms%coarse_mat == mld_repl_mat_) then - - 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_) then - call psb_errpush(psb_err_from_subroutine_,name,a_err='psb_cdall') - goto 9999 - end if - call psb_gather(p%ac,b,p%desc_ac,info,dupl=psb_dupl_add_,keeploc=.false.) - if(info /= psb_success_) goto 9999 - - else if (p%parms%coarse_mat == mld_distr_mat_) then - - nzl = b%get_nzeros() - call b%mv_to(bcoo) - - 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 - call op_restr%set_nrows(p%desc_ac%get_local_cols()) - - if (debug_level >= psb_debug_outer_) & - & write(debug_unit,*) me,' ',trim(name),& - & 'Done ac ' - else - info = psb_err_internal_error_ - call psb_errpush(psb_err_internal_error_,name,a_err='invalid mld_coarse_mat_') - goto 9999 - end if - - - 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='cscnv') - 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 - ! - if (info == psb_success_) & - & 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='linmap build') - goto 9999 - end if + call ac_coo%set_nrows(naggr) + call ac_coo%set_ncols(naggr) + call ac_coo%set_dupl(psb_dupl_add_) + call ac_coo%fix(info) + call ac%mv_from(ac_coo) call psb_erractionrestore(err_act) diff --git a/mlprec/impl/mld_daggrmat_smth_asb.f90 b/mlprec/impl/mld_daggrmat_smth_asb.f90 index f1c3126e..af616ece 100644 --- a/mlprec/impl/mld_daggrmat_smth_asb.f90 +++ b/mlprec/impl/mld_daggrmat_smth_asb.f90 @@ -93,7 +93,7 @@ ! info - integer, output. ! Error code. ! -subroutine mld_daggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info) +subroutine mld_daggrmat_smth_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_smth_asb @@ -103,9 +103,9 @@ subroutine mld_daggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info) type(psb_dspmat_type), intent(in) :: a type(psb_desc_type), intent(in) :: desc_a integer, intent(inout) :: ilaggr(:), nlaggr(:) - type(mld_d_onelev_type), intent(inout), target :: p + type(mld_dml_parms), intent(inout) :: parms + type(psb_dspmat_type), intent(out) :: ac,op_prol,op_restr integer, intent(out) :: info - type(psb_dspmat_type) :: op_prol, op_restr, b ! Local variables integer :: nrow, nglob, ncol, ntaggr, ip, ndx,& @@ -113,7 +113,7 @@ subroutine mld_daggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info) integer ::ictxt, np, me, err_act character(len=20) :: name type(psb_dspmat_type) :: am3, am4 - type(psb_d_coo_sparse_mat) :: acoo, acoof, bcoo + type(psb_d_coo_sparse_mat) :: tmpcoo type(psb_d_csr_sparse_mat) :: acsr1, acsr2, acsr3, acsrf, ptilde real(psb_dpk_), allocatable :: adiag(:) integer(psb_ipk_) :: ierr(5) @@ -138,14 +138,14 @@ subroutine mld_daggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info) nrow = desc_a%get_local_rows() ncol = desc_a%get_local_cols() - theta = p%parms%aggr_thresh + theta = parms%aggr_thresh naggr = nlaggr(me+1) ntaggr = sum(nlaggr) naggrm1 = sum(nlaggr(1:me)) naggrp1 = sum(nlaggr(1:me+1)) - filter_mat = (p%parms%aggr_filter == mld_filter_mat_) + filter_mat = (parms%aggr_filter == mld_filter_mat_) ilaggr(1:nrow) = ilaggr(1:nrow) + naggrm1 call psb_halo(ilaggr,desc_a,info) @@ -177,16 +177,16 @@ subroutine mld_daggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info) end if ! 1. Allocate Ptilde in sparse matrix form - call acoo%allocate(ncol,ntaggr,ncol) + call tmpcoo%allocate(ncol,ntaggr,ncol) do i=1,ncol - acoo%val(i) = done - acoo%ia(i) = i - acoo%ja(i) = ilaggr(i) + tmpcoo%val(i) = done + tmpcoo%ia(i) = i + tmpcoo%ja(i) = ilaggr(i) end do - call acoo%set_nzeros(ncol) - call acoo%set_dupl(psb_dupl_add_) + call tmpcoo%set_nzeros(ncol) + call tmpcoo%set_dupl(psb_dupl_add_) - call ptilde%mv_from_coo(acoo,info) + call ptilde%mv_from_coo(tmpcoo,info) if (info == psb_success_) call a%cscnv(acsr3,info,dupl=psb_dupl_add_) if (debug_level >= psb_debug_outer_) & @@ -217,19 +217,19 @@ subroutine mld_daggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info) end if enddo ! Take out zeroed terms - call acsrf%mv_to_coo(acoof,info) + call acsrf%mv_to_coo(tmpcoo,info) k = 0 - do j=1,acoof%get_nzeros() - if ((acoof%val(j) /= dzero) .or. (acoof%ia(j) == acoof%ja(j))) then + do j=1,tmpcoo%get_nzeros() + if ((tmpcoo%val(j) /= dzero) .or. (tmpcoo%ia(j) == tmpcoo%ja(j))) then k = k + 1 - acoof%val(k) = acoof%val(j) - acoof%ia(k) = acoof%ia(j) - acoof%ja(k) = acoof%ja(j) + tmpcoo%val(k) = tmpcoo%val(j) + tmpcoo%ia(k) = tmpcoo%ia(j) + tmpcoo%ja(k) = tmpcoo%ja(j) end if end do - call acoof%set_nzeros(k) - call acoof%set_dupl(psb_dupl_add_) - call acsrf%mv_from_coo(acoof,info) + call tmpcoo%set_nzeros(k) + call tmpcoo%set_dupl(psb_dupl_add_) + call acsrf%mv_from_coo(tmpcoo,info) end if @@ -246,13 +246,13 @@ subroutine mld_daggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info) if (info /= psb_success_) goto 9999 - if (p%parms%aggr_omega_alg == mld_eig_est_) then + if (parms%aggr_omega_alg == mld_eig_est_) then - if (p%parms%aggr_eig == mld_max_norm_) then + if (parms%aggr_eig == mld_max_norm_) then anorm = acsr3%csnmi() omega = 4.d0/(3.d0*anorm) - p%parms%aggr_omega_val = omega + parms%aggr_omega_val = omega else info = psb_err_internal_error_ @@ -260,11 +260,11 @@ subroutine mld_daggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info) goto 9999 end if - else if (p%parms%aggr_omega_alg == mld_user_choice_) then + else if (parms%aggr_omega_alg == mld_user_choice_) then - omega = p%parms%aggr_omega_val + omega = parms%aggr_omega_val - else if (p%parms%aggr_omega_alg /= mld_user_choice_) then + else if (parms%aggr_omega_alg /= mld_user_choice_) then info = psb_err_internal_error_ call psb_errpush(info,name,a_err='invalid mld_aggr_omega_alg_') goto 9999 @@ -369,27 +369,27 @@ subroutine mld_daggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info) call psb_numbmm(a,op_prol,am3) if (debug_level >= psb_debug_outer_) & & write(debug_unit,*) me,' ',trim(name),& - & 'Done NUMBMM 2',p%parms%aggr_kind, mld_smooth_prol_ + & 'Done NUMBMM 2',parms%aggr_kind, mld_smooth_prol_ call op_prol%transp(op_restr) - call op_restr%mv_to(acoo) - nzl = acoo%get_nzeros() + call op_restr%mv_to(tmpcoo) + nzl = tmpcoo%get_nzeros() i=0 ! ! Now we have to fix this. The only rows of B that are correct ! are those corresponding to "local" aggregates, i.e. indices in ilaggr(:) ! do k=1, nzl - if ((naggrm1 < acoo%ia(k)) .and.(acoo%ia(k) <= naggrp1)) then + if ((naggrm1 < tmpcoo%ia(k)) .and.(tmpcoo%ia(k) <= naggrp1)) then i = i+1 - acoo%val(i) = acoo%val(k) - acoo%ia(i) = acoo%ia(k) - acoo%ja(i) = acoo%ja(k) + tmpcoo%val(i) = tmpcoo%val(k) + tmpcoo%ia(i) = tmpcoo%ia(k) + tmpcoo%ja(i) = tmpcoo%ja(k) end if end do - call acoo%set_nzeros(i) - call acoo%trim() - call op_restr%mv_from(acoo) + call tmpcoo%set_nzeros(i) + call tmpcoo%trim() + call op_restr%mv_from(tmpcoo) call op_restr%cscnv(info,type='csr',dupl=psb_dupl_add_) if (info /= psb_success_) then call psb_errpush(psb_err_from_subroutine_,name,a_err='spcnv op_restr') @@ -413,113 +413,12 @@ subroutine mld_daggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info) if (debug_level >= psb_debug_outer_) & & write(debug_unit,*) me,' ',trim(name),& & 'starting symbmm 3' - call psb_symbmm(op_restr,am3,b,info) - if (info == psb_success_) call psb_numbmm(op_restr,am3,b) + call psb_symbmm(op_restr,am3,ac,info) + if (info == psb_success_) call psb_numbmm(op_restr,am3,ac) if (info == psb_success_) call am3%free() - if (info == psb_success_) call b%cscnv(info,type='coo',dupl=psb_dupl_add_) + if (info == psb_success_) call ac%cscnv(info,type='coo',dupl=psb_dupl_add_) if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,a_err='Build b = op_restr x am3') - goto 9999 - end if - - - - select case(p%parms%coarse_mat) - - case(mld_distr_mat_) - - nzl = b%get_nzeros() - call b%mv_to(bcoo) - - 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 - call op_restr%set_nrows(p%desc_ac%get_local_cols()) - - 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_) & - & call psb_gather(p%ac,b,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') + call psb_errpush(psb_err_internal_error_,name,a_err='Build ac = op_restr x am3') goto 9999 end if diff --git a/mlprec/impl/mld_saggrmat_asb.f90 b/mlprec/impl/mld_saggrmat_asb.f90 index 9ce3b701..5a84f116 100644 --- a/mlprec/impl/mld_saggrmat_asb.f90 +++ b/mlprec/impl/mld_saggrmat_asb.f90 @@ -113,7 +113,11 @@ subroutine mld_saggrmat_asb(a,desc_a,ilaggr,nlaggr,p,info) integer, intent(out) :: info ! Local variables - type(psb_sspmat_type) :: b, op_prol,op_restr + type(psb_sspmat_type) :: ac, op_prol,op_restr + type(psb_s_coo_sparse_mat) :: acoo, bcoo + type(psb_s_csr_sparse_mat) :: acsr1 + integer :: nzl,ntaggr + integer :: debug_level, debug_unit integer :: ictxt,np,me, err_act character(len=20) :: name @@ -121,6 +125,9 @@ subroutine mld_saggrmat_asb(a,desc_a,ilaggr,nlaggr,p,info) if(psb_get_errstatus().ne.0) return info=psb_success_ call psb_erractionsave(err_act) + debug_unit = psb_get_debug_unit() + debug_level = psb_get_debug_level() + ictxt = desc_a%get_context() @@ -129,43 +136,139 @@ subroutine mld_saggrmat_asb(a,desc_a,ilaggr,nlaggr,p,info) select case (p%parms%aggr_kind) case (mld_no_smooth_) - call mld_saggrmat_nosmth_asb(a,desc_a,ilaggr,nlaggr,p,info) - if(info /= psb_success_) then - call psb_errpush(psb_err_from_subroutine_,name,a_err='mld_aggrmat_nosmth_asb') - goto 9999 - end if + call mld_saggrmat_nosmth_asb(a,desc_a,ilaggr,nlaggr,& + & p%parms,ac,op_prol,op_restr,info) case(mld_smooth_prol_) - call mld_saggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info) - if(info /= psb_success_) then - call psb_errpush(psb_err_from_subroutine_,name,a_err='mld_aggrmat_smth_asb') - goto 9999 - end if + call mld_saggrmat_smth_asb(a,desc_a,ilaggr,nlaggr, & + & p%parms,ac,op_prol,op_restr,info) case(mld_biz_prol_) - call mld_saggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,p,info) - if(info /= psb_success_) then - call psb_errpush(psb_err_from_subroutine_,name,a_err='mld_aggrmat_smth_asb') + call mld_saggrmat_biz_asb(a,desc_a,ilaggr,nlaggr, & + & p%parms,ac,op_prol,op_restr,info) + + case(mld_min_energy_) + + call mld_saggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr, & + & p%parms,ac,op_prol,op_restr,info) + + case default + info = psb_err_internal_error_ + call psb_errpush(info,name,a_err='Invalid aggr kind') + goto 9999 + + end select + if (info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,a_err='Inner aggrmat asb') + goto 9999 + end if + + + + ntaggr = sum(nlaggr) + + select case(p%parms%coarse_mat) + + case(mld_distr_mat_) + + nzl = ac%get_nzeros() + call ac%mv_to(bcoo) + + 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) - case(mld_min_energy_) + 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() - call mld_saggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info) - if(info /= psb_success_) then - call psb_errpush(psb_err_from_subroutine_,name,a_err='mld_aggrmat_smth_asb') + if (info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,a_err='psb_sp_free') goto 9999 end if - case default + 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()) - call psb_errpush(psb_err_internal_error_,name,a_err='Invalid aggr kind') - goto 9999 + 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 + call op_restr%set_nrows(p%desc_ac%get_local_cols()) + + 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_) & + & 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) return diff --git a/mlprec/impl/mld_saggrmat_biz_asb.f90 b/mlprec/impl/mld_saggrmat_biz_asb.f90 index e1bf0552..df43f098 100644 --- a/mlprec/impl/mld_saggrmat_biz_asb.f90 +++ b/mlprec/impl/mld_saggrmat_biz_asb.f90 @@ -78,7 +78,7 @@ ! info - integer, output. ! Error code. ! -subroutine mld_saggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,p,info) +subroutine mld_saggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_restr,info) use psb_base_mod use mld_s_inner_mod, mld_protect_name => mld_saggrmat_biz_asb @@ -88,9 +88,9 @@ subroutine mld_saggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,p,info) type(psb_sspmat_type), intent(in) :: a type(psb_desc_type), intent(in) :: desc_a integer, intent(inout) :: ilaggr(:), nlaggr(:) - type(mld_s_onelev_type), intent(inout), target :: p + type(mld_sml_parms), intent(inout) :: parms + type(psb_sspmat_type), intent(out) :: ac,op_prol,op_restr integer, intent(out) :: info - type(psb_sspmat_type) :: op_prol, op_restr, b ! Local variables integer :: nrow, nglob, ncol, ntaggr, ip, ndx,& @@ -98,7 +98,7 @@ subroutine mld_saggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,p,info) integer ::ictxt, np, me, err_act character(len=20) :: name type(psb_sspmat_type) :: am3, am4 - type(psb_s_coo_sparse_mat) :: acoo, acoof, bcoo + type(psb_s_coo_sparse_mat) :: tmpcoo type(psb_s_csr_sparse_mat) :: acsr1, acsr2, acsr3, acsrf, ptilde real(psb_spk_), allocatable :: adiag(:) integer(psb_ipk_) :: ierr(5) @@ -123,7 +123,7 @@ subroutine mld_saggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,p,info) nrow = desc_a%get_local_rows() ncol = desc_a%get_local_cols() - theta = p%parms%aggr_thresh + theta = parms%aggr_thresh naggr = nlaggr(me+1) ntaggr = sum(nlaggr) @@ -131,7 +131,7 @@ subroutine mld_saggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,p,info) naggrm1 = sum(nlaggr(1:me)) naggrp1 = sum(nlaggr(1:me+1)) - filter_mat = (p%parms%aggr_filter == mld_filter_mat_) + filter_mat = (parms%aggr_filter == mld_filter_mat_) ilaggr(1:nrow) = ilaggr(1:nrow) + naggrm1 call psb_halo(ilaggr,desc_a,info) @@ -163,16 +163,16 @@ subroutine mld_saggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,p,info) end if ! 1. Allocate Ptilde in sparse matrix form - call acoo%allocate(ncol,naggr,ncol) + call tmpcoo%allocate(ncol,naggr,ncol) do i=1,nrow - acoo%val(i) = sone - acoo%ia(i) = i - acoo%ja(i) = ilaggr(i) + tmpcoo%val(i) = sone + tmpcoo%ia(i) = i + tmpcoo%ja(i) = ilaggr(i) end do - call acoo%set_nzeros(nrow) - call acoo%set_dupl(psb_dupl_add_) + call tmpcoo%set_nzeros(nrow) + call tmpcoo%set_dupl(psb_dupl_add_) - call ptilde%mv_from_coo(acoo,info) + call ptilde%mv_from_coo(tmpcoo,info) if (info == psb_success_) call a%cscnv(acsr3,info,dupl=psb_dupl_add_) if (debug_level >= psb_debug_outer_) & @@ -203,19 +203,19 @@ subroutine mld_saggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,p,info) end if enddo ! Take out zeroed terms - call acsrf%mv_to_coo(acoof,info) + call acsrf%mv_to_coo(tmpcoo,info) k = 0 - do j=1,acoof%get_nzeros() - if ((acoof%val(j) /= szero) .or. (acoof%ia(j) == acoof%ja(j))) then + do j=1,tmpcoo%get_nzeros() + if ((tmpcoo%val(j) /= szero) .or. (tmpcoo%ia(j) == tmpcoo%ja(j))) then k = k + 1 - acoof%val(k) = acoof%val(j) - acoof%ia(k) = acoof%ia(j) - acoof%ja(k) = acoof%ja(j) + tmpcoo%val(k) = tmpcoo%val(j) + tmpcoo%ia(k) = tmpcoo%ia(j) + tmpcoo%ja(k) = tmpcoo%ja(j) end if end do - call acoof%set_nzeros(k) - call acoof%set_dupl(psb_dupl_add_) - call acsrf%mv_from_coo(acoof,info) + call tmpcoo%set_nzeros(k) + call tmpcoo%set_dupl(psb_dupl_add_) + call acsrf%mv_from_coo(tmpcoo,info) end if @@ -232,9 +232,9 @@ subroutine mld_saggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,p,info) if (info /= psb_success_) goto 9999 - if (p%parms%aggr_omega_alg == mld_eig_est_) then + if (parms%aggr_omega_alg == mld_eig_est_) then - if (p%parms%aggr_eig == mld_max_norm_) then + if (parms%aggr_eig == mld_max_norm_) then ! ! This only works with CSR @@ -261,7 +261,7 @@ subroutine mld_saggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,p,info) goto 9999 end if omega = 4.d0/(3.d0*anorm) - p%parms%aggr_omega_val = omega + parms%aggr_omega_val = omega else info = psb_err_internal_error_ @@ -269,11 +269,11 @@ subroutine mld_saggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,p,info) goto 9999 end if - else if (p%parms%aggr_omega_alg == mld_user_choice_) then + else if (parms%aggr_omega_alg == mld_user_choice_) then - omega = p%parms%aggr_omega_val + omega = parms%aggr_omega_val - else if (p%parms%aggr_omega_alg /= mld_user_choice_) then + else if (parms%aggr_omega_alg /= mld_user_choice_) then info = psb_err_internal_error_ call psb_errpush(info,name,a_err='invalid mld_aggr_omega_alg_') goto 9999 @@ -372,7 +372,7 @@ subroutine mld_saggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,p,info) call psb_numbmm(a,op_prol,am3) if (debug_level >= psb_debug_outer_) & & write(debug_unit,*) me,' ',trim(name),& - & 'Done NUMBMM 2',p%parms%aggr_kind, mld_smooth_prol_ + & 'Done NUMBMM 2',parms%aggr_kind, mld_smooth_prol_ call op_prol%transp(op_restr) if (debug_level >= psb_debug_outer_) & @@ -389,10 +389,10 @@ subroutine mld_saggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,p,info) if (debug_level >= psb_debug_outer_) & & write(debug_unit,*) me,' ',trim(name),& & 'starting symbmm 3' - call psb_symbmm(op_restr,am3,b,info) - if (info == psb_success_) call psb_numbmm(op_restr,am3,b) + call psb_symbmm(op_restr,am3,ac,info) + if (info == psb_success_) call psb_numbmm(op_restr,am3,ac) if (info == psb_success_) call am3%free() - if (info == psb_success_) call b%cscnv(info,type='coo',dupl=psb_dupl_add_) + if (info == psb_success_) call ac%cscnv(info,type='coo',dupl=psb_dupl_add_) if (info /= psb_success_) then call psb_errpush(psb_err_internal_error_,name,a_err='Build b = op_restr x am3') goto 9999 @@ -401,107 +401,6 @@ subroutine mld_saggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,p,info) - select case(p%parms%coarse_mat) - - case(mld_distr_mat_) - - nzl = b%get_nzeros() - call b%mv_to(bcoo) - - 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 - call op_restr%set_nrows(p%desc_ac%get_local_cols()) - - 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_) & - & call psb_gather(p%ac,b,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 - if (debug_level >= psb_debug_outer_) & & write(debug_unit,*) me,' ',trim(name),& diff --git a/mlprec/impl/mld_saggrmat_minnrg_asb.f90 b/mlprec/impl/mld_saggrmat_minnrg_asb.f90 index e9b46654..4e5ed449 100644 --- a/mlprec/impl/mld_saggrmat_minnrg_asb.f90 +++ b/mlprec/impl/mld_saggrmat_minnrg_asb.f90 @@ -98,7 +98,7 @@ ! info - integer, output. ! Error code. ! -subroutine mld_saggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info) +subroutine mld_saggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_restr,info) use psb_base_mod use mld_s_inner_mod, mld_protect_name => mld_saggrmat_minnrg_asb @@ -108,9 +108,9 @@ subroutine mld_saggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info) type(psb_sspmat_type), intent(in) :: a type(psb_desc_type), intent(in) :: desc_a integer, intent(inout) :: ilaggr(:), nlaggr(:) - type(mld_s_onelev_type), intent(inout), target :: p + type(mld_sml_parms), intent(inout) :: parms + type(psb_sspmat_type), intent(out) :: ac,op_prol,op_restr integer, intent(out) :: info - type(psb_sspmat_type) :: op_prol,op_restr, b ! Local variables integer(psb_mpik_), allocatable :: nzbr(:), idisp(:) @@ -121,8 +121,8 @@ subroutine mld_saggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info) type(psb_sspmat_type) :: af, ptilde, rtilde, atran, atp, atdatp type(psb_sspmat_type) :: am3,am4, ap, adap,atmp,rada, ra, atmp2, dap, dadap, da type(psb_sspmat_type) :: dat, datp, datdatp, atmp3 - type(psb_s_coo_sparse_mat) :: acoo, acoof, bcoo, tmpcoo - type(psb_s_csr_sparse_mat) :: acsr1, acsr2, acsr3, bcsr, acsr, acsrf + type(psb_s_coo_sparse_mat) :: tmpcoo + type(psb_s_csr_sparse_mat) :: acsr1, acsr2, acsr3, acsr, acsrf type(psb_s_csc_sparse_mat) :: csc_dap, csc_dadap, csc_datp, csc_datdatp, acsc real(psb_spk_), allocatable :: adiag(:), adinv(:) real(psb_spk_), allocatable :: omf(:), omp(:), omi(:), oden(:) @@ -150,7 +150,7 @@ subroutine mld_saggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info) nrow = desc_a%get_local_rows() ncol = desc_a%get_local_cols() - theta = p%parms%aggr_thresh + theta = parms%aggr_thresh naggr = nlaggr(me+1) ntaggr = sum(nlaggr) @@ -165,7 +165,7 @@ subroutine mld_saggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info) naggrm1 = sum(nlaggr(1:me)) naggrp1 = sum(nlaggr(1:me+1)) - filter_mat = (p%parms%aggr_filter == mld_filter_mat_) + filter_mat = (parms%aggr_filter == mld_filter_mat_) ilaggr(1:nrow) = ilaggr(1:nrow) + naggrm1 call psb_halo(ilaggr,desc_a,info) @@ -207,16 +207,16 @@ subroutine mld_saggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info) ! 1. Allocate Ptilde in sparse matrix form - call acoo%allocate(ncol,ntaggr,ncol) + call tmpcoo%allocate(ncol,ntaggr,ncol) do i=1,ncol - acoo%val(i) = sone - acoo%ia(i) = i - acoo%ja(i) = ilaggr(i) + tmpcoo%val(i) = sone + tmpcoo%ia(i) = i + tmpcoo%ja(i) = ilaggr(i) end do - call acoo%set_nzeros(ncol) - call acoo%set_dupl(psb_dupl_add_) - call acoo%set_asb() - call ptilde%mv_from(acoo) + call tmpcoo%set_nzeros(ncol) + call tmpcoo%set_dupl(psb_dupl_add_) + call tmpcoo%set_asb() + call ptilde%mv_from(tmpcoo) call ptilde%cscnv(info,type='csr') !!$ call local_dump(me,ptilde,'csr-ptilde','Ptilde-1') @@ -570,164 +570,18 @@ subroutine mld_saggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info) & write(debug_unit,*) me,' ',trim(name),& & 'Done sphalo/ rwxtd' - call psb_symbmm(op_restr,am3,b,info) - if (info == psb_success_) call psb_numbmm(op_restr,am3,b) + call psb_symbmm(op_restr,am3,ac,info) + if (info == psb_success_) call psb_numbmm(op_restr,am3,ac) if (info == psb_success_) call am3%free() + if (info == psb_success_) call ac%cscnv(info,type='coo',dupl=psb_dupl_add_) if (info /= psb_success_) then call psb_errpush(psb_err_internal_error_,name,& - &a_err='Build b = op_restr x am3') + &a_err='Build ac = op_restr x am3') goto 9999 end if - if (debug_level >= psb_debug_outer_) & - & write(debug_unit,*) me,' ',trim(name),& - & 'Done mv_to_coo' - - - select case(p%parms%coarse_mat) - - case(mld_distr_mat_) - - call b%mv_to(bcoo) - nzl = bcoo%get_nzeros() - - if (debug_level >= psb_debug_outer_) & - & write(debug_unit,*) me,' ',trim(name),& - & ' B matrix nzl: ',nzl - - 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 bcoo%set_nrows(p%desc_ac%get_local_rows()) - call bcoo%set_ncols(p%desc_ac%get_local_cols()) - call bcoo%fix(info) - call p%ac%mv_from(bcoo) - call p%ac%set_asb() - - call p%ac%cscnv(info,type='csr') - - if (np>1) then - call op_prol%mv_to(acsr) - nzl = acsr%get_nzeros() - call psb_glob_to_loc(acsr%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(acsr) - endif - call op_prol%set_ncols(p%desc_ac%get_local_cols()) - - if (np>1) then - call op_restr%mv_to(tmpcoo) - nzl = tmpcoo%get_nzeros() - if (info == psb_success_) call psb_glob_to_loc(tmpcoo%ia(1:nzl),p%desc_ac,info,'I') - if(info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,a_err='Converting op_restr to local') - goto 9999 - end if - call op_restr%mv_from(tmpcoo) - call op_restr%cscnv(info,type='csr') - end if - call op_restr%set_nrows(p%desc_ac%get_local_cols()) - - 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_) then - call psb_errpush(psb_err_from_subroutine_,name,a_err='psb_cdall') - goto 9999 - end if - call psb_gather(p%ac,b,p%desc_ac,info,dupl=psb_dupl_add_,keeploc=.false.) - if(info /= psb_success_) goto 9999 -!!$ call psb_cdall(ictxt,p%desc_ac,info,mg=ntaggr,repl=.true.) -!!$ nzbr(:) = 0 -!!$ nzbr(me+1) = bcoo%get_nzeros() -!!$ -!!$ call psb_sum(ictxt,nzbr(1:np)) -!!$ nzac = sum(nzbr) -!!$ if (info == psb_success_) call tmpcoo%allocate(ntaggr,ntaggr,nzac) -!!$ if (info /= psb_success_) goto 9999 -!!$ -!!$ do ip=1,np -!!$ idispip) = sum(nzbr(1:ip-1)) -!!$ enddo -!!$ ndx = nzbr(me+1) -!!$ -!!$ call mpi_allgatherv(bcoo%val,ndx,mpi_real,tmpcoo%val,nzbr,idisp,& -!!$ & mpi_real,icomm,info) -!!$ if (info == psb_success_)& -!!$ & call mpi_allgatherv(bcoo%ia,ndx,psb_mpi_ipk_integer,tmpcoo%ia,nzbr,idisp,& -!!$ & psb_mpi_ipk_integer,icomm,info) -!!$ if (info == psb_success_)& -!!$ & call mpi_allgatherv(bcoo%ja,ndx,psb_mpi_ipk_integer,tmpcoo%ja,nzbr,idisp,& -!!$ & psb_mpi_ipk_integer,icomm,info) -!!$ -!!$ if (info /= psb_success_) then -!!$ call psb_errpush(psb_err_internal_error_,name,& -!!$ & a_err=' from mpi_allgatherv') -!!$ goto 9999 -!!$ end if -!!$ -!!$ call bcoo%free() -!!$ call tmpcoo%fix(info) -!!$ call p%ac%mv_from(tmpcoo) -!!$ call p%ac%cscnv(info,type='csr') -!!$ if(info /= psb_success_) goto 9999 -!!$ -!!$ deallocate(nzbr,idisp,stat=info) -!!$ if (info /= psb_success_) then -!!$ info = psb_err_alloc_dealloc_ -!!$ call psb_errpush(info,name) -!!$ goto 9999 -!!$ end if - 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') - 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 => R i.e. restriction operator - ! op_prol => P 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 if (debug_level >= psb_debug_outer_) & & write(debug_unit,*) me,' ',trim(name),& diff --git a/mlprec/impl/mld_saggrmat_nosmth_asb.f90 b/mlprec/impl/mld_saggrmat_nosmth_asb.f90 index e7c26afe..6a10da69 100644 --- a/mlprec/impl/mld_saggrmat_nosmth_asb.f90 +++ b/mlprec/impl/mld_saggrmat_nosmth_asb.f90 @@ -81,7 +81,7 @@ ! info - integer, output. ! Error code. ! -subroutine mld_saggrmat_nosmth_asb(a,desc_a,ilaggr,nlaggr,p,info) +subroutine mld_saggrmat_nosmth_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_restr,info) use psb_base_mod use mld_s_inner_mod, mld_protect_name => mld_saggrmat_nosmth_asb @@ -91,16 +91,16 @@ subroutine mld_saggrmat_nosmth_asb(a,desc_a,ilaggr,nlaggr,p,info) type(psb_sspmat_type), intent(in) :: a type(psb_desc_type), intent(in) :: desc_a integer, intent(inout) :: ilaggr(:), nlaggr(:) - type(mld_s_onelev_type), intent(inout), target :: p + type(mld_sml_parms), intent(inout) :: parms + type(psb_sspmat_type), intent(out) :: ac,op_prol,op_restr integer, intent(out) :: info - type(psb_sspmat_type) :: b, op_prol,op_restr ! Local variables integer :: ictxt,np,me, err_act integer(psb_mpik_) :: icomm, ndx, minfo character(len=20) :: name integer(psb_ipk_) :: ierr(5) - type(psb_s_coo_sparse_mat) :: acoo1, acoo2, bcoo, ac_coo, acoo + type(psb_s_coo_sparse_mat) :: ac_coo, acoo type(psb_s_csr_sparse_mat) :: acsr1, acsr2 integer :: debug_level, debug_unit integer :: nrow, nglob, ncol, ntaggr, nzl, ip, & @@ -134,136 +134,36 @@ subroutine mld_saggrmat_nosmth_asb(a,desc_a,ilaggr,nlaggr,p,info) goto 9999 end if - call acoo1%allocate(ncol,ntaggr,ncol) + call acoo%allocate(ncol,ntaggr,ncol) do i=1,nrow - acoo1%val(i) = sone - acoo1%ia(i) = i - acoo1%ja(i) = ilaggr(i) + acoo%val(i) = sone + acoo%ia(i) = i + acoo%ja(i) = ilaggr(i) end do - call acoo1%set_dupl(psb_dupl_add_) - call acoo1%set_nzeros(nrow) - call acoo1%set_asb() - call acoo1%fix(info) + call acoo%set_dupl(psb_dupl_add_) + call acoo%set_nzeros(nrow) + call acoo%set_asb() + call acoo%fix(info) - call op_prol%mv_from(acoo1) + call op_prol%mv_from(acoo) call op_prol%cscnv(info,type='csr',dupl=psb_dupl_add_) if (info == psb_success_) call op_prol%transp(op_restr) - call a%csclip(bcoo,info,jmax=nrow) + call a%csclip(ac_coo,info,jmax=nrow) - nzt = bcoo%get_nzeros() + nzt = ac_coo%get_nzeros() do i=1, nzt - bcoo%ia(i) = ilaggr(bcoo%ia(i)) - bcoo%ja(i) = ilaggr(bcoo%ja(i)) + ac_coo%ia(i) = ilaggr(ac_coo%ia(i)) + ac_coo%ja(i) = ilaggr(ac_coo%ja(i)) enddo - call bcoo%set_nrows(naggr) - call bcoo%set_ncols(naggr) - call bcoo%set_dupl(psb_dupl_add_) - call bcoo%fix(info) - - - call b%mv_from(bcoo) - - if (p%parms%coarse_mat == mld_repl_mat_) then - - 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_) then - call psb_errpush(psb_err_from_subroutine_,name,a_err='psb_cdall') - goto 9999 - end if - call psb_gather(p%ac,b,p%desc_ac,info,dupl=psb_dupl_add_,keeploc=.false.) - if(info /= psb_success_) goto 9999 - - else if (p%parms%coarse_mat == mld_distr_mat_) then - - nzl = b%get_nzeros() - call b%mv_to(bcoo) - - 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 - call op_restr%set_nrows(p%desc_ac%get_local_cols()) - - if (debug_level >= psb_debug_outer_) & - & write(debug_unit,*) me,' ',trim(name),& - & 'Done ac ' - else - info = psb_err_internal_error_ - call psb_errpush(psb_err_internal_error_,name,a_err='invalid mld_coarse_mat_') - goto 9999 - end if - - - 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='cscnv') - 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 - ! - if (info == psb_success_) & - & 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='linmap build') - goto 9999 - end if + call ac_coo%set_nrows(naggr) + call ac_coo%set_ncols(naggr) + call ac_coo%set_dupl(psb_dupl_add_) + call ac_coo%fix(info) + call ac%mv_from(ac_coo) call psb_erractionrestore(err_act) diff --git a/mlprec/impl/mld_saggrmat_smth_asb.f90 b/mlprec/impl/mld_saggrmat_smth_asb.f90 index 532d7a2a..c38d9ec6 100644 --- a/mlprec/impl/mld_saggrmat_smth_asb.f90 +++ b/mlprec/impl/mld_saggrmat_smth_asb.f90 @@ -93,7 +93,7 @@ ! info - integer, output. ! Error code. ! -subroutine mld_saggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info) +subroutine mld_saggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_restr,info) use psb_base_mod use mld_s_inner_mod, mld_protect_name => mld_saggrmat_smth_asb @@ -103,9 +103,9 @@ subroutine mld_saggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info) type(psb_sspmat_type), intent(in) :: a type(psb_desc_type), intent(in) :: desc_a integer, intent(inout) :: ilaggr(:), nlaggr(:) - type(mld_s_onelev_type), intent(inout), target :: p + type(mld_sml_parms), intent(inout) :: parms + type(psb_sspmat_type), intent(out) :: ac,op_prol,op_restr integer, intent(out) :: info - type(psb_sspmat_type) :: op_prol, op_restr, b ! Local variables integer :: nrow, nglob, ncol, ntaggr, ip, ndx,& @@ -113,7 +113,7 @@ subroutine mld_saggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info) integer ::ictxt, np, me, err_act character(len=20) :: name type(psb_sspmat_type) :: am3, am4 - type(psb_s_coo_sparse_mat) :: acoo, acoof, bcoo + type(psb_s_coo_sparse_mat) :: tmpcoo type(psb_s_csr_sparse_mat) :: acsr1, acsr2, acsr3, acsrf, ptilde real(psb_spk_), allocatable :: adiag(:) integer(psb_ipk_) :: ierr(5) @@ -138,14 +138,14 @@ subroutine mld_saggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info) nrow = desc_a%get_local_rows() ncol = desc_a%get_local_cols() - theta = p%parms%aggr_thresh + theta = parms%aggr_thresh naggr = nlaggr(me+1) ntaggr = sum(nlaggr) naggrm1 = sum(nlaggr(1:me)) naggrp1 = sum(nlaggr(1:me+1)) - filter_mat = (p%parms%aggr_filter == mld_filter_mat_) + filter_mat = (parms%aggr_filter == mld_filter_mat_) ilaggr(1:nrow) = ilaggr(1:nrow) + naggrm1 call psb_halo(ilaggr,desc_a,info) @@ -177,16 +177,16 @@ subroutine mld_saggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info) end if ! 1. Allocate Ptilde in sparse matrix form - call acoo%allocate(ncol,ntaggr,ncol) + call tmpcoo%allocate(ncol,ntaggr,ncol) do i=1,ncol - acoo%val(i) = sone - acoo%ia(i) = i - acoo%ja(i) = ilaggr(i) + tmpcoo%val(i) = sone + tmpcoo%ia(i) = i + tmpcoo%ja(i) = ilaggr(i) end do - call acoo%set_nzeros(ncol) - call acoo%set_dupl(psb_dupl_add_) + call tmpcoo%set_nzeros(ncol) + call tmpcoo%set_dupl(psb_dupl_add_) - call ptilde%mv_from_coo(acoo,info) + call ptilde%mv_from_coo(tmpcoo,info) if (info == psb_success_) call a%cscnv(acsr3,info,dupl=psb_dupl_add_) if (debug_level >= psb_debug_outer_) & @@ -217,19 +217,19 @@ subroutine mld_saggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info) end if enddo ! Take out zeroed terms - call acsrf%mv_to_coo(acoof,info) + call acsrf%mv_to_coo(tmpcoo,info) k = 0 - do j=1,acoof%get_nzeros() - if ((acoof%val(j) /= szero) .or. (acoof%ia(j) == acoof%ja(j))) then + do j=1,tmpcoo%get_nzeros() + if ((tmpcoo%val(j) /= szero) .or. (tmpcoo%ia(j) == tmpcoo%ja(j))) then k = k + 1 - acoof%val(k) = acoof%val(j) - acoof%ia(k) = acoof%ia(j) - acoof%ja(k) = acoof%ja(j) + tmpcoo%val(k) = tmpcoo%val(j) + tmpcoo%ia(k) = tmpcoo%ia(j) + tmpcoo%ja(k) = tmpcoo%ja(j) end if end do - call acoof%set_nzeros(k) - call acoof%set_dupl(psb_dupl_add_) - call acsrf%mv_from_coo(acoof,info) + call tmpcoo%set_nzeros(k) + call tmpcoo%set_dupl(psb_dupl_add_) + call acsrf%mv_from_coo(tmpcoo,info) end if @@ -246,13 +246,13 @@ subroutine mld_saggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info) if (info /= psb_success_) goto 9999 - if (p%parms%aggr_omega_alg == mld_eig_est_) then + if (parms%aggr_omega_alg == mld_eig_est_) then - if (p%parms%aggr_eig == mld_max_norm_) then + if (parms%aggr_eig == mld_max_norm_) then anorm = acsr3%csnmi() omega = 4.d0/(3.d0*anorm) - p%parms%aggr_omega_val = omega + parms%aggr_omega_val = omega else info = psb_err_internal_error_ @@ -260,11 +260,11 @@ subroutine mld_saggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info) goto 9999 end if - else if (p%parms%aggr_omega_alg == mld_user_choice_) then + else if (parms%aggr_omega_alg == mld_user_choice_) then - omega = p%parms%aggr_omega_val + omega = parms%aggr_omega_val - else if (p%parms%aggr_omega_alg /= mld_user_choice_) then + else if (parms%aggr_omega_alg /= mld_user_choice_) then info = psb_err_internal_error_ call psb_errpush(info,name,a_err='invalid mld_aggr_omega_alg_') goto 9999 @@ -369,27 +369,27 @@ subroutine mld_saggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info) call psb_numbmm(a,op_prol,am3) if (debug_level >= psb_debug_outer_) & & write(debug_unit,*) me,' ',trim(name),& - & 'Done NUMBMM 2',p%parms%aggr_kind, mld_smooth_prol_ + & 'Done NUMBMM 2',parms%aggr_kind, mld_smooth_prol_ call op_prol%transp(op_restr) - call op_restr%mv_to(acoo) - nzl = acoo%get_nzeros() + call op_restr%mv_to(tmpcoo) + nzl = tmpcoo%get_nzeros() i=0 ! ! Now we have to fix this. The only rows of B that are correct ! are those corresponding to "local" aggregates, i.e. indices in ilaggr(:) ! do k=1, nzl - if ((naggrm1 < acoo%ia(k)) .and.(acoo%ia(k) <= naggrp1)) then + if ((naggrm1 < tmpcoo%ia(k)) .and.(tmpcoo%ia(k) <= naggrp1)) then i = i+1 - acoo%val(i) = acoo%val(k) - acoo%ia(i) = acoo%ia(k) - acoo%ja(i) = acoo%ja(k) + tmpcoo%val(i) = tmpcoo%val(k) + tmpcoo%ia(i) = tmpcoo%ia(k) + tmpcoo%ja(i) = tmpcoo%ja(k) end if end do - call acoo%set_nzeros(i) - call acoo%trim() - call op_restr%mv_from(acoo) + call tmpcoo%set_nzeros(i) + call tmpcoo%trim() + call op_restr%mv_from(tmpcoo) call op_restr%cscnv(info,type='csr',dupl=psb_dupl_add_) if (info /= psb_success_) then call psb_errpush(psb_err_from_subroutine_,name,a_err='spcnv op_restr') @@ -413,113 +413,12 @@ subroutine mld_saggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info) if (debug_level >= psb_debug_outer_) & & write(debug_unit,*) me,' ',trim(name),& & 'starting symbmm 3' - call psb_symbmm(op_restr,am3,b,info) - if (info == psb_success_) call psb_numbmm(op_restr,am3,b) + call psb_symbmm(op_restr,am3,ac,info) + if (info == psb_success_) call psb_numbmm(op_restr,am3,ac) if (info == psb_success_) call am3%free() - if (info == psb_success_) call b%cscnv(info,type='coo',dupl=psb_dupl_add_) + if (info == psb_success_) call ac%cscnv(info,type='coo',dupl=psb_dupl_add_) if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,a_err='Build b = op_restr x am3') - goto 9999 - end if - - - - select case(p%parms%coarse_mat) - - case(mld_distr_mat_) - - nzl = b%get_nzeros() - call b%mv_to(bcoo) - - 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 - call op_restr%set_nrows(p%desc_ac%get_local_cols()) - - 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_) & - & call psb_gather(p%ac,b,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') + call psb_errpush(psb_err_internal_error_,name,a_err='Build ac = op_restr x am3') goto 9999 end if diff --git a/mlprec/impl/mld_zaggrmat_asb.f90 b/mlprec/impl/mld_zaggrmat_asb.f90 index 383e10a6..f0d0e55a 100644 --- a/mlprec/impl/mld_zaggrmat_asb.f90 +++ b/mlprec/impl/mld_zaggrmat_asb.f90 @@ -113,7 +113,11 @@ subroutine mld_zaggrmat_asb(a,desc_a,ilaggr,nlaggr,p,info) integer, intent(out) :: info ! Local variables - type(psb_zspmat_type) :: b, op_prol,op_restr + type(psb_zspmat_type) :: ac, op_prol,op_restr + type(psb_z_coo_sparse_mat) :: acoo, bcoo + type(psb_z_csr_sparse_mat) :: acsr1 + integer :: nzl,ntaggr + integer :: debug_level, debug_unit integer :: ictxt,np,me, err_act character(len=20) :: name @@ -121,6 +125,9 @@ subroutine mld_zaggrmat_asb(a,desc_a,ilaggr,nlaggr,p,info) if(psb_get_errstatus().ne.0) return info=psb_success_ call psb_erractionsave(err_act) + debug_unit = psb_get_debug_unit() + debug_level = psb_get_debug_level() + ictxt = desc_a%get_context() @@ -129,43 +136,139 @@ subroutine mld_zaggrmat_asb(a,desc_a,ilaggr,nlaggr,p,info) select case (p%parms%aggr_kind) case (mld_no_smooth_) - call mld_zaggrmat_nosmth_asb(a,desc_a,ilaggr,nlaggr,p,info) - if(info /= psb_success_) then - call psb_errpush(psb_err_from_subroutine_,name,a_err='mld_aggrmat_nosmth_asb') - goto 9999 - end if + call mld_zaggrmat_nosmth_asb(a,desc_a,ilaggr,nlaggr,& + & p%parms,ac,op_prol,op_restr,info) case(mld_smooth_prol_) - call mld_zaggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info) - if(info /= psb_success_) then - call psb_errpush(psb_err_from_subroutine_,name,a_err='mld_aggrmat_smth_asb') - goto 9999 - end if + call mld_zaggrmat_smth_asb(a,desc_a,ilaggr,nlaggr, & + & p%parms,ac,op_prol,op_restr,info) case(mld_biz_prol_) - call mld_zaggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,p,info) - if(info /= psb_success_) then - call psb_errpush(psb_err_from_subroutine_,name,a_err='mld_aggrmat_smth_asb') + call mld_zaggrmat_biz_asb(a,desc_a,ilaggr,nlaggr, & + & p%parms,ac,op_prol,op_restr,info) + + case(mld_min_energy_) + + call mld_zaggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr, & + & p%parms,ac,op_prol,op_restr,info) + + case default + info = psb_err_internal_error_ + call psb_errpush(info,name,a_err='Invalid aggr kind') + goto 9999 + + end select + if (info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,a_err='Inner aggrmat asb') + goto 9999 + end if + + + + ntaggr = sum(nlaggr) + + select case(p%parms%coarse_mat) + + case(mld_distr_mat_) + + nzl = ac%get_nzeros() + call ac%mv_to(bcoo) + + 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) - case(mld_min_energy_) + 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() - call mld_zaggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info) - if(info /= psb_success_) then - call psb_errpush(psb_err_from_subroutine_,name,a_err='mld_aggrmat_smth_asb') + if (info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,a_err='psb_sp_free') goto 9999 end if - case default + 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()) - call psb_errpush(psb_err_internal_error_,name,a_err='Invalid aggr kind') - goto 9999 + 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 + call op_restr%set_nrows(p%desc_ac%get_local_cols()) + + 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_) & + & 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) return diff --git a/mlprec/impl/mld_zaggrmat_biz_asb.f90 b/mlprec/impl/mld_zaggrmat_biz_asb.f90 index 3f43e05c..f8d945f7 100644 --- a/mlprec/impl/mld_zaggrmat_biz_asb.f90 +++ b/mlprec/impl/mld_zaggrmat_biz_asb.f90 @@ -78,7 +78,7 @@ ! info - integer, output. ! Error code. ! -subroutine mld_zaggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,p,info) +subroutine mld_zaggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_restr,info) use psb_base_mod use mld_z_inner_mod, mld_protect_name => mld_zaggrmat_biz_asb @@ -88,9 +88,9 @@ subroutine mld_zaggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,p,info) type(psb_zspmat_type), intent(in) :: a type(psb_desc_type), intent(in) :: desc_a integer, intent(inout) :: ilaggr(:), nlaggr(:) - type(mld_z_onelev_type), intent(inout), target :: p + type(mld_dml_parms), intent(inout) :: parms + type(psb_zspmat_type), intent(out) :: ac,op_prol,op_restr integer, intent(out) :: info - type(psb_zspmat_type) :: op_prol, op_restr, b ! Local variables integer :: nrow, nglob, ncol, ntaggr, ip, ndx,& @@ -98,7 +98,7 @@ subroutine mld_zaggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,p,info) integer ::ictxt, np, me, err_act character(len=20) :: name type(psb_zspmat_type) :: am3, am4 - type(psb_z_coo_sparse_mat) :: acoo, acoof, bcoo + type(psb_z_coo_sparse_mat) :: tmpcoo type(psb_z_csr_sparse_mat) :: acsr1, acsr2, acsr3, acsrf, ptilde complex(psb_dpk_), allocatable :: adiag(:) integer(psb_ipk_) :: ierr(5) @@ -123,7 +123,7 @@ subroutine mld_zaggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,p,info) nrow = desc_a%get_local_rows() ncol = desc_a%get_local_cols() - theta = p%parms%aggr_thresh + theta = parms%aggr_thresh naggr = nlaggr(me+1) ntaggr = sum(nlaggr) @@ -131,7 +131,7 @@ subroutine mld_zaggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,p,info) naggrm1 = sum(nlaggr(1:me)) naggrp1 = sum(nlaggr(1:me+1)) - filter_mat = (p%parms%aggr_filter == mld_filter_mat_) + filter_mat = (parms%aggr_filter == mld_filter_mat_) ilaggr(1:nrow) = ilaggr(1:nrow) + naggrm1 call psb_halo(ilaggr,desc_a,info) @@ -163,16 +163,16 @@ subroutine mld_zaggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,p,info) end if ! 1. Allocate Ptilde in sparse matrix form - call acoo%allocate(ncol,naggr,ncol) + call tmpcoo%allocate(ncol,naggr,ncol) do i=1,nrow - acoo%val(i) = zone - acoo%ia(i) = i - acoo%ja(i) = ilaggr(i) + tmpcoo%val(i) = zone + tmpcoo%ia(i) = i + tmpcoo%ja(i) = ilaggr(i) end do - call acoo%set_nzeros(nrow) - call acoo%set_dupl(psb_dupl_add_) + call tmpcoo%set_nzeros(nrow) + call tmpcoo%set_dupl(psb_dupl_add_) - call ptilde%mv_from_coo(acoo,info) + call ptilde%mv_from_coo(tmpcoo,info) if (info == psb_success_) call a%cscnv(acsr3,info,dupl=psb_dupl_add_) if (debug_level >= psb_debug_outer_) & @@ -203,19 +203,19 @@ subroutine mld_zaggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,p,info) end if enddo ! Take out zeroed terms - call acsrf%mv_to_coo(acoof,info) + call acsrf%mv_to_coo(tmpcoo,info) k = 0 - do j=1,acoof%get_nzeros() - if ((acoof%val(j) /= zzero) .or. (acoof%ia(j) == acoof%ja(j))) then + do j=1,tmpcoo%get_nzeros() + if ((tmpcoo%val(j) /= zzero) .or. (tmpcoo%ia(j) == tmpcoo%ja(j))) then k = k + 1 - acoof%val(k) = acoof%val(j) - acoof%ia(k) = acoof%ia(j) - acoof%ja(k) = acoof%ja(j) + tmpcoo%val(k) = tmpcoo%val(j) + tmpcoo%ia(k) = tmpcoo%ia(j) + tmpcoo%ja(k) = tmpcoo%ja(j) end if end do - call acoof%set_nzeros(k) - call acoof%set_dupl(psb_dupl_add_) - call acsrf%mv_from_coo(acoof,info) + call tmpcoo%set_nzeros(k) + call tmpcoo%set_dupl(psb_dupl_add_) + call acsrf%mv_from_coo(tmpcoo,info) end if @@ -232,9 +232,9 @@ subroutine mld_zaggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,p,info) if (info /= psb_success_) goto 9999 - if (p%parms%aggr_omega_alg == mld_eig_est_) then + if (parms%aggr_omega_alg == mld_eig_est_) then - if (p%parms%aggr_eig == mld_max_norm_) then + if (parms%aggr_eig == mld_max_norm_) then ! ! This only works with CSR @@ -261,7 +261,7 @@ subroutine mld_zaggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,p,info) goto 9999 end if omega = 4.d0/(3.d0*anorm) - p%parms%aggr_omega_val = omega + parms%aggr_omega_val = omega else info = psb_err_internal_error_ @@ -269,11 +269,11 @@ subroutine mld_zaggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,p,info) goto 9999 end if - else if (p%parms%aggr_omega_alg == mld_user_choice_) then + else if (parms%aggr_omega_alg == mld_user_choice_) then - omega = p%parms%aggr_omega_val + omega = parms%aggr_omega_val - else if (p%parms%aggr_omega_alg /= mld_user_choice_) then + else if (parms%aggr_omega_alg /= mld_user_choice_) then info = psb_err_internal_error_ call psb_errpush(info,name,a_err='invalid mld_aggr_omega_alg_') goto 9999 @@ -372,7 +372,7 @@ subroutine mld_zaggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,p,info) call psb_numbmm(a,op_prol,am3) if (debug_level >= psb_debug_outer_) & & write(debug_unit,*) me,' ',trim(name),& - & 'Done NUMBMM 2',p%parms%aggr_kind, mld_smooth_prol_ + & 'Done NUMBMM 2',parms%aggr_kind, mld_smooth_prol_ call op_prol%transp(op_restr) if (debug_level >= psb_debug_outer_) & @@ -389,10 +389,10 @@ subroutine mld_zaggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,p,info) if (debug_level >= psb_debug_outer_) & & write(debug_unit,*) me,' ',trim(name),& & 'starting symbmm 3' - call psb_symbmm(op_restr,am3,b,info) - if (info == psb_success_) call psb_numbmm(op_restr,am3,b) + call psb_symbmm(op_restr,am3,ac,info) + if (info == psb_success_) call psb_numbmm(op_restr,am3,ac) if (info == psb_success_) call am3%free() - if (info == psb_success_) call b%cscnv(info,type='coo',dupl=psb_dupl_add_) + if (info == psb_success_) call ac%cscnv(info,type='coo',dupl=psb_dupl_add_) if (info /= psb_success_) then call psb_errpush(psb_err_internal_error_,name,a_err='Build b = op_restr x am3') goto 9999 @@ -401,107 +401,6 @@ subroutine mld_zaggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,p,info) - select case(p%parms%coarse_mat) - - case(mld_distr_mat_) - - nzl = b%get_nzeros() - call b%mv_to(bcoo) - - 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 - call op_restr%set_nrows(p%desc_ac%get_local_cols()) - - 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_) & - & call psb_gather(p%ac,b,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 - if (debug_level >= psb_debug_outer_) & & write(debug_unit,*) me,' ',trim(name),& diff --git a/mlprec/impl/mld_zaggrmat_minnrg_asb.f90 b/mlprec/impl/mld_zaggrmat_minnrg_asb.f90 index 834a7526..fc76615a 100644 --- a/mlprec/impl/mld_zaggrmat_minnrg_asb.f90 +++ b/mlprec/impl/mld_zaggrmat_minnrg_asb.f90 @@ -98,7 +98,7 @@ ! info - integer, output. ! Error code. ! -subroutine mld_zaggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info) +subroutine mld_zaggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_restr,info) use psb_base_mod use mld_z_inner_mod, mld_protect_name => mld_zaggrmat_minnrg_asb @@ -108,9 +108,9 @@ subroutine mld_zaggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info) type(psb_zspmat_type), intent(in) :: a type(psb_desc_type), intent(in) :: desc_a integer, intent(inout) :: ilaggr(:), nlaggr(:) - type(mld_z_onelev_type), intent(inout), target :: p + type(mld_dml_parms), intent(inout) :: parms + type(psb_zspmat_type), intent(out) :: ac,op_prol,op_restr integer, intent(out) :: info - type(psb_zspmat_type) :: op_prol,op_restr, b ! Local variables integer(psb_mpik_), allocatable :: nzbr(:), idisp(:) @@ -121,8 +121,8 @@ subroutine mld_zaggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info) type(psb_zspmat_type) :: af, ptilde, rtilde, atran, atp, atdatp type(psb_zspmat_type) :: am3,am4, ap, adap,atmp,rada, ra, atmp2, dap, dadap, da type(psb_zspmat_type) :: dat, datp, datdatp, atmp3 - type(psb_z_coo_sparse_mat) :: acoo, acoof, bcoo, tmpcoo - type(psb_z_csr_sparse_mat) :: acsr1, acsr2, acsr3, bcsr, acsr, acsrf + type(psb_z_coo_sparse_mat) :: tmpcoo + type(psb_z_csr_sparse_mat) :: acsr1, acsr2, acsr3, acsr, acsrf type(psb_z_csc_sparse_mat) :: csc_dap, csc_dadap, csc_datp, csc_datdatp, acsc complex(psb_dpk_), allocatable :: adiag(:), adinv(:) complex(psb_dpk_), allocatable :: omf(:), omp(:), omi(:), oden(:) @@ -150,7 +150,7 @@ subroutine mld_zaggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info) nrow = desc_a%get_local_rows() ncol = desc_a%get_local_cols() - theta = p%parms%aggr_thresh + theta = parms%aggr_thresh naggr = nlaggr(me+1) ntaggr = sum(nlaggr) @@ -165,7 +165,7 @@ subroutine mld_zaggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info) naggrm1 = sum(nlaggr(1:me)) naggrp1 = sum(nlaggr(1:me+1)) - filter_mat = (p%parms%aggr_filter == mld_filter_mat_) + filter_mat = (parms%aggr_filter == mld_filter_mat_) ilaggr(1:nrow) = ilaggr(1:nrow) + naggrm1 call psb_halo(ilaggr,desc_a,info) @@ -207,16 +207,16 @@ subroutine mld_zaggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info) ! 1. Allocate Ptilde in sparse matrix form - call acoo%allocate(ncol,ntaggr,ncol) + call tmpcoo%allocate(ncol,ntaggr,ncol) do i=1,ncol - acoo%val(i) = zone - acoo%ia(i) = i - acoo%ja(i) = ilaggr(i) + tmpcoo%val(i) = zone + tmpcoo%ia(i) = i + tmpcoo%ja(i) = ilaggr(i) end do - call acoo%set_nzeros(ncol) - call acoo%set_dupl(psb_dupl_add_) - call acoo%set_asb() - call ptilde%mv_from(acoo) + call tmpcoo%set_nzeros(ncol) + call tmpcoo%set_dupl(psb_dupl_add_) + call tmpcoo%set_asb() + call ptilde%mv_from(tmpcoo) call ptilde%cscnv(info,type='csr') !!$ call local_dump(me,ptilde,'csr-ptilde','Ptilde-1') @@ -570,164 +570,18 @@ subroutine mld_zaggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info) & write(debug_unit,*) me,' ',trim(name),& & 'Done sphalo/ rwxtd' - call psb_symbmm(op_restr,am3,b,info) - if (info == psb_success_) call psb_numbmm(op_restr,am3,b) + call psb_symbmm(op_restr,am3,ac,info) + if (info == psb_success_) call psb_numbmm(op_restr,am3,ac) if (info == psb_success_) call am3%free() + if (info == psb_success_) call ac%cscnv(info,type='coo',dupl=psb_dupl_add_) if (info /= psb_success_) then call psb_errpush(psb_err_internal_error_,name,& - &a_err='Build b = op_restr x am3') + &a_err='Build ac = op_restr x am3') goto 9999 end if - if (debug_level >= psb_debug_outer_) & - & write(debug_unit,*) me,' ',trim(name),& - & 'Done mv_to_coo' - - - select case(p%parms%coarse_mat) - - case(mld_distr_mat_) - - call b%mv_to(bcoo) - nzl = bcoo%get_nzeros() - - if (debug_level >= psb_debug_outer_) & - & write(debug_unit,*) me,' ',trim(name),& - & ' B matrix nzl: ',nzl - - 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 bcoo%set_nrows(p%desc_ac%get_local_rows()) - call bcoo%set_ncols(p%desc_ac%get_local_cols()) - call bcoo%fix(info) - call p%ac%mv_from(bcoo) - call p%ac%set_asb() - - call p%ac%cscnv(info,type='csr') - - if (np>1) then - call op_prol%mv_to(acsr) - nzl = acsr%get_nzeros() - call psb_glob_to_loc(acsr%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(acsr) - endif - call op_prol%set_ncols(p%desc_ac%get_local_cols()) - - if (np>1) then - call op_restr%mv_to(tmpcoo) - nzl = tmpcoo%get_nzeros() - if (info == psb_success_) call psb_glob_to_loc(tmpcoo%ia(1:nzl),p%desc_ac,info,'I') - if(info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,a_err='Converting op_restr to local') - goto 9999 - end if - call op_restr%mv_from(tmpcoo) - call op_restr%cscnv(info,type='csr') - end if - call op_restr%set_nrows(p%desc_ac%get_local_cols()) - - 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_) then - call psb_errpush(psb_err_from_subroutine_,name,a_err='psb_cdall') - goto 9999 - end if - call psb_gather(p%ac,b,p%desc_ac,info,dupl=psb_dupl_add_,keeploc=.false.) - if(info /= psb_success_) goto 9999 -!!$ call psb_cdall(ictxt,p%desc_ac,info,mg=ntaggr,repl=.true.) -!!$ nzbr(:) = 0 -!!$ nzbr(me+1) = bcoo%get_nzeros() -!!$ -!!$ call psb_sum(ictxt,nzbr(1:np)) -!!$ nzac = sum(nzbr) -!!$ if (info == psb_success_) call tmpcoo%allocate(ntaggr,ntaggr,nzac) -!!$ if (info /= psb_success_) goto 9999 -!!$ -!!$ do ip=1,np -!!$ idispip) = sum(nzbr(1:ip-1)) -!!$ enddo -!!$ ndx = nzbr(me+1) -!!$ -!!$ call mpi_allgatherv(bcoo%val,ndx,mpi_double_complex,tmpcoo%val,nzbr,idisp,& -!!$ & mpi_double_complex,icomm,info) -!!$ if (info == psb_success_)& -!!$ & call mpi_allgatherv(bcoo%ia,ndx,psb_mpi_ipk_integer,tmpcoo%ia,nzbr,idisp,& -!!$ & psb_mpi_ipk_integer,icomm,info) -!!$ if (info == psb_success_)& -!!$ & call mpi_allgatherv(bcoo%ja,ndx,psb_mpi_ipk_integer,tmpcoo%ja,nzbr,idisp,& -!!$ & psb_mpi_ipk_integer,icomm,info) -!!$ -!!$ if (info /= psb_success_) then -!!$ call psb_errpush(psb_err_internal_error_,name,& -!!$ & a_err=' from mpi_allgatherv') -!!$ goto 9999 -!!$ end if -!!$ -!!$ call bcoo%free() -!!$ call tmpcoo%fix(info) -!!$ call p%ac%mv_from(tmpcoo) -!!$ call p%ac%cscnv(info,type='csr') -!!$ if(info /= psb_success_) goto 9999 -!!$ -!!$ deallocate(nzbr,idisp,stat=info) -!!$ if (info /= psb_success_) then -!!$ info = psb_err_alloc_dealloc_ -!!$ call psb_errpush(info,name) -!!$ goto 9999 -!!$ end if - 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') - 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 => R i.e. restriction operator - ! op_prol => P 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 if (debug_level >= psb_debug_outer_) & & write(debug_unit,*) me,' ',trim(name),& diff --git a/mlprec/impl/mld_zaggrmat_nosmth_asb.f90 b/mlprec/impl/mld_zaggrmat_nosmth_asb.f90 index e83cc323..85239315 100644 --- a/mlprec/impl/mld_zaggrmat_nosmth_asb.f90 +++ b/mlprec/impl/mld_zaggrmat_nosmth_asb.f90 @@ -81,7 +81,7 @@ ! info - integer, output. ! Error code. ! -subroutine mld_zaggrmat_nosmth_asb(a,desc_a,ilaggr,nlaggr,p,info) +subroutine mld_zaggrmat_nosmth_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_restr,info) use psb_base_mod use mld_z_inner_mod, mld_protect_name => mld_zaggrmat_nosmth_asb @@ -91,16 +91,16 @@ subroutine mld_zaggrmat_nosmth_asb(a,desc_a,ilaggr,nlaggr,p,info) type(psb_zspmat_type), intent(in) :: a type(psb_desc_type), intent(in) :: desc_a integer, intent(inout) :: ilaggr(:), nlaggr(:) - type(mld_z_onelev_type), intent(inout), target :: p + type(mld_dml_parms), intent(inout) :: parms + type(psb_zspmat_type), intent(out) :: ac,op_prol,op_restr integer, intent(out) :: info - type(psb_zspmat_type) :: b, op_prol,op_restr ! Local variables integer :: ictxt,np,me, err_act integer(psb_mpik_) :: icomm, ndx, minfo character(len=20) :: name integer(psb_ipk_) :: ierr(5) - type(psb_z_coo_sparse_mat) :: acoo1, acoo2, bcoo, ac_coo, acoo + type(psb_z_coo_sparse_mat) :: ac_coo, acoo type(psb_z_csr_sparse_mat) :: acsr1, acsr2 integer :: debug_level, debug_unit integer :: nrow, nglob, ncol, ntaggr, nzl, ip, & @@ -134,136 +134,36 @@ subroutine mld_zaggrmat_nosmth_asb(a,desc_a,ilaggr,nlaggr,p,info) goto 9999 end if - call acoo1%allocate(ncol,ntaggr,ncol) + call acoo%allocate(ncol,ntaggr,ncol) do i=1,nrow - acoo1%val(i) = zone - acoo1%ia(i) = i - acoo1%ja(i) = ilaggr(i) + acoo%val(i) = zone + acoo%ia(i) = i + acoo%ja(i) = ilaggr(i) end do - call acoo1%set_dupl(psb_dupl_add_) - call acoo1%set_nzeros(nrow) - call acoo1%set_asb() - call acoo1%fix(info) + call acoo%set_dupl(psb_dupl_add_) + call acoo%set_nzeros(nrow) + call acoo%set_asb() + call acoo%fix(info) - call op_prol%mv_from(acoo1) + call op_prol%mv_from(acoo) call op_prol%cscnv(info,type='csr',dupl=psb_dupl_add_) if (info == psb_success_) call op_prol%transp(op_restr) - call a%csclip(bcoo,info,jmax=nrow) + call a%csclip(ac_coo,info,jmax=nrow) - nzt = bcoo%get_nzeros() + nzt = ac_coo%get_nzeros() do i=1, nzt - bcoo%ia(i) = ilaggr(bcoo%ia(i)) - bcoo%ja(i) = ilaggr(bcoo%ja(i)) + ac_coo%ia(i) = ilaggr(ac_coo%ia(i)) + ac_coo%ja(i) = ilaggr(ac_coo%ja(i)) enddo - call bcoo%set_nrows(naggr) - call bcoo%set_ncols(naggr) - call bcoo%set_dupl(psb_dupl_add_) - call bcoo%fix(info) - - - call b%mv_from(bcoo) - - if (p%parms%coarse_mat == mld_repl_mat_) then - - 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_) then - call psb_errpush(psb_err_from_subroutine_,name,a_err='psb_cdall') - goto 9999 - end if - call psb_gather(p%ac,b,p%desc_ac,info,dupl=psb_dupl_add_,keeploc=.false.) - if(info /= psb_success_) goto 9999 - - else if (p%parms%coarse_mat == mld_distr_mat_) then - - nzl = b%get_nzeros() - call b%mv_to(bcoo) - - 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 - call op_restr%set_nrows(p%desc_ac%get_local_cols()) - - if (debug_level >= psb_debug_outer_) & - & write(debug_unit,*) me,' ',trim(name),& - & 'Done ac ' - else - info = psb_err_internal_error_ - call psb_errpush(psb_err_internal_error_,name,a_err='invalid mld_coarse_mat_') - goto 9999 - end if - - - 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='cscnv') - 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 - ! - if (info == psb_success_) & - & 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='linmap build') - goto 9999 - end if + call ac_coo%set_nrows(naggr) + call ac_coo%set_ncols(naggr) + call ac_coo%set_dupl(psb_dupl_add_) + call ac_coo%fix(info) + call ac%mv_from(ac_coo) call psb_erractionrestore(err_act) diff --git a/mlprec/impl/mld_zaggrmat_smth_asb.f90 b/mlprec/impl/mld_zaggrmat_smth_asb.f90 index e5f72174..b33be1ce 100644 --- a/mlprec/impl/mld_zaggrmat_smth_asb.f90 +++ b/mlprec/impl/mld_zaggrmat_smth_asb.f90 @@ -93,7 +93,7 @@ ! info - integer, output. ! Error code. ! -subroutine mld_zaggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info) +subroutine mld_zaggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_restr,info) use psb_base_mod use mld_z_inner_mod, mld_protect_name => mld_zaggrmat_smth_asb @@ -103,9 +103,9 @@ subroutine mld_zaggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info) type(psb_zspmat_type), intent(in) :: a type(psb_desc_type), intent(in) :: desc_a integer, intent(inout) :: ilaggr(:), nlaggr(:) - type(mld_z_onelev_type), intent(inout), target :: p + type(mld_dml_parms), intent(inout) :: parms + type(psb_zspmat_type), intent(out) :: ac,op_prol,op_restr integer, intent(out) :: info - type(psb_zspmat_type) :: op_prol, op_restr, b ! Local variables integer :: nrow, nglob, ncol, ntaggr, ip, ndx,& @@ -113,7 +113,7 @@ subroutine mld_zaggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info) integer ::ictxt, np, me, err_act character(len=20) :: name type(psb_zspmat_type) :: am3, am4 - type(psb_z_coo_sparse_mat) :: acoo, acoof, bcoo + type(psb_z_coo_sparse_mat) :: tmpcoo type(psb_z_csr_sparse_mat) :: acsr1, acsr2, acsr3, acsrf, ptilde complex(psb_dpk_), allocatable :: adiag(:) integer(psb_ipk_) :: ierr(5) @@ -138,14 +138,14 @@ subroutine mld_zaggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info) nrow = desc_a%get_local_rows() ncol = desc_a%get_local_cols() - theta = p%parms%aggr_thresh + theta = parms%aggr_thresh naggr = nlaggr(me+1) ntaggr = sum(nlaggr) naggrm1 = sum(nlaggr(1:me)) naggrp1 = sum(nlaggr(1:me+1)) - filter_mat = (p%parms%aggr_filter == mld_filter_mat_) + filter_mat = (parms%aggr_filter == mld_filter_mat_) ilaggr(1:nrow) = ilaggr(1:nrow) + naggrm1 call psb_halo(ilaggr,desc_a,info) @@ -177,16 +177,16 @@ subroutine mld_zaggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info) end if ! 1. Allocate Ptilde in sparse matrix form - call acoo%allocate(ncol,ntaggr,ncol) + call tmpcoo%allocate(ncol,ntaggr,ncol) do i=1,ncol - acoo%val(i) = zone - acoo%ia(i) = i - acoo%ja(i) = ilaggr(i) + tmpcoo%val(i) = zone + tmpcoo%ia(i) = i + tmpcoo%ja(i) = ilaggr(i) end do - call acoo%set_nzeros(ncol) - call acoo%set_dupl(psb_dupl_add_) + call tmpcoo%set_nzeros(ncol) + call tmpcoo%set_dupl(psb_dupl_add_) - call ptilde%mv_from_coo(acoo,info) + call ptilde%mv_from_coo(tmpcoo,info) if (info == psb_success_) call a%cscnv(acsr3,info,dupl=psb_dupl_add_) if (debug_level >= psb_debug_outer_) & @@ -217,19 +217,19 @@ subroutine mld_zaggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info) end if enddo ! Take out zeroed terms - call acsrf%mv_to_coo(acoof,info) + call acsrf%mv_to_coo(tmpcoo,info) k = 0 - do j=1,acoof%get_nzeros() - if ((acoof%val(j) /= zzero) .or. (acoof%ia(j) == acoof%ja(j))) then + do j=1,tmpcoo%get_nzeros() + if ((tmpcoo%val(j) /= zzero) .or. (tmpcoo%ia(j) == tmpcoo%ja(j))) then k = k + 1 - acoof%val(k) = acoof%val(j) - acoof%ia(k) = acoof%ia(j) - acoof%ja(k) = acoof%ja(j) + tmpcoo%val(k) = tmpcoo%val(j) + tmpcoo%ia(k) = tmpcoo%ia(j) + tmpcoo%ja(k) = tmpcoo%ja(j) end if end do - call acoof%set_nzeros(k) - call acoof%set_dupl(psb_dupl_add_) - call acsrf%mv_from_coo(acoof,info) + call tmpcoo%set_nzeros(k) + call tmpcoo%set_dupl(psb_dupl_add_) + call acsrf%mv_from_coo(tmpcoo,info) end if @@ -246,13 +246,13 @@ subroutine mld_zaggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info) if (info /= psb_success_) goto 9999 - if (p%parms%aggr_omega_alg == mld_eig_est_) then + if (parms%aggr_omega_alg == mld_eig_est_) then - if (p%parms%aggr_eig == mld_max_norm_) then + if (parms%aggr_eig == mld_max_norm_) then anorm = acsr3%csnmi() omega = 4.d0/(3.d0*anorm) - p%parms%aggr_omega_val = omega + parms%aggr_omega_val = omega else info = psb_err_internal_error_ @@ -260,11 +260,11 @@ subroutine mld_zaggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info) goto 9999 end if - else if (p%parms%aggr_omega_alg == mld_user_choice_) then + else if (parms%aggr_omega_alg == mld_user_choice_) then - omega = p%parms%aggr_omega_val + omega = parms%aggr_omega_val - else if (p%parms%aggr_omega_alg /= mld_user_choice_) then + else if (parms%aggr_omega_alg /= mld_user_choice_) then info = psb_err_internal_error_ call psb_errpush(info,name,a_err='invalid mld_aggr_omega_alg_') goto 9999 @@ -369,27 +369,27 @@ subroutine mld_zaggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info) call psb_numbmm(a,op_prol,am3) if (debug_level >= psb_debug_outer_) & & write(debug_unit,*) me,' ',trim(name),& - & 'Done NUMBMM 2',p%parms%aggr_kind, mld_smooth_prol_ + & 'Done NUMBMM 2',parms%aggr_kind, mld_smooth_prol_ call op_prol%transp(op_restr) - call op_restr%mv_to(acoo) - nzl = acoo%get_nzeros() + call op_restr%mv_to(tmpcoo) + nzl = tmpcoo%get_nzeros() i=0 ! ! Now we have to fix this. The only rows of B that are correct ! are those corresponding to "local" aggregates, i.e. indices in ilaggr(:) ! do k=1, nzl - if ((naggrm1 < acoo%ia(k)) .and.(acoo%ia(k) <= naggrp1)) then + if ((naggrm1 < tmpcoo%ia(k)) .and.(tmpcoo%ia(k) <= naggrp1)) then i = i+1 - acoo%val(i) = acoo%val(k) - acoo%ia(i) = acoo%ia(k) - acoo%ja(i) = acoo%ja(k) + tmpcoo%val(i) = tmpcoo%val(k) + tmpcoo%ia(i) = tmpcoo%ia(k) + tmpcoo%ja(i) = tmpcoo%ja(k) end if end do - call acoo%set_nzeros(i) - call acoo%trim() - call op_restr%mv_from(acoo) + call tmpcoo%set_nzeros(i) + call tmpcoo%trim() + call op_restr%mv_from(tmpcoo) call op_restr%cscnv(info,type='csr',dupl=psb_dupl_add_) if (info /= psb_success_) then call psb_errpush(psb_err_from_subroutine_,name,a_err='spcnv op_restr') @@ -413,113 +413,12 @@ subroutine mld_zaggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info) if (debug_level >= psb_debug_outer_) & & write(debug_unit,*) me,' ',trim(name),& & 'starting symbmm 3' - call psb_symbmm(op_restr,am3,b,info) - if (info == psb_success_) call psb_numbmm(op_restr,am3,b) + call psb_symbmm(op_restr,am3,ac,info) + if (info == psb_success_) call psb_numbmm(op_restr,am3,ac) if (info == psb_success_) call am3%free() - if (info == psb_success_) call b%cscnv(info,type='coo',dupl=psb_dupl_add_) + if (info == psb_success_) call ac%cscnv(info,type='coo',dupl=psb_dupl_add_) if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,a_err='Build b = op_restr x am3') - goto 9999 - end if - - - - select case(p%parms%coarse_mat) - - case(mld_distr_mat_) - - nzl = b%get_nzeros() - call b%mv_to(bcoo) - - 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 - call op_restr%set_nrows(p%desc_ac%get_local_cols()) - - 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_) & - & call psb_gather(p%ac,b,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') + call psb_errpush(psb_err_internal_error_,name,a_err='Build ac = op_restr x am3') goto 9999 end if diff --git a/mlprec/mld_c_inner_mod.f90 b/mlprec/mld_c_inner_mod.f90 index 27f89f04..566ab38f 100644 --- a/mlprec/mld_c_inner_mod.f90 +++ b/mlprec/mld_c_inner_mod.f90 @@ -146,13 +146,14 @@ module mld_c_inner_mod abstract interface - subroutine mld_caggrmat_var_asb(a,desc_a,ilaggr,nlaggr,p,info) + subroutine mld_caggrmat_var_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_restr,info) use psb_base_mod, only : psb_cspmat_type, psb_desc_type, psb_spk_ - use mld_c_prec_type, only : mld_c_onelev_type + use mld_c_prec_type, only : mld_c_onelev_type, mld_sml_parms type(psb_cspmat_type), intent(in) :: a type(psb_desc_type), intent(in) :: desc_a integer, intent(inout) :: ilaggr(:), nlaggr(:) - type(mld_c_onelev_type), intent(inout), target :: p + type(mld_sml_parms), intent(inout) :: parms + type(psb_cspmat_type), intent(out) :: ac,op_prol,op_restr integer, intent(out) :: info end subroutine mld_caggrmat_var_asb end interface diff --git a/mlprec/mld_d_inner_mod.f90 b/mlprec/mld_d_inner_mod.f90 index b32bf8c7..5278ec15 100644 --- a/mlprec/mld_d_inner_mod.f90 +++ b/mlprec/mld_d_inner_mod.f90 @@ -146,13 +146,14 @@ module mld_d_inner_mod abstract interface - subroutine mld_daggrmat_var_asb(a,desc_a,ilaggr,nlaggr,p,info) + subroutine mld_daggrmat_var_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_restr,info) use psb_base_mod, only : psb_dspmat_type, psb_desc_type, psb_dpk_ - use mld_d_prec_type, only : mld_d_onelev_type + use mld_d_prec_type, only : mld_d_onelev_type, mld_dml_parms type(psb_dspmat_type), intent(in) :: a type(psb_desc_type), intent(in) :: desc_a integer, intent(inout) :: ilaggr(:), nlaggr(:) - type(mld_d_onelev_type), intent(inout), target :: p + type(mld_dml_parms), intent(inout) :: parms + type(psb_dspmat_type), intent(out) :: ac,op_prol,op_restr integer, intent(out) :: info end subroutine mld_daggrmat_var_asb end interface diff --git a/mlprec/mld_s_inner_mod.f90 b/mlprec/mld_s_inner_mod.f90 index 9cc9461c..d6b3badf 100644 --- a/mlprec/mld_s_inner_mod.f90 +++ b/mlprec/mld_s_inner_mod.f90 @@ -146,13 +146,14 @@ module mld_s_inner_mod abstract interface - subroutine mld_saggrmat_var_asb(a,desc_a,ilaggr,nlaggr,p,info) + subroutine mld_saggrmat_var_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_restr,info) use psb_base_mod, only : psb_sspmat_type, psb_desc_type, psb_spk_ - use mld_s_prec_type, only : mld_s_onelev_type + use mld_s_prec_type, only : mld_s_onelev_type, mld_sml_parms type(psb_sspmat_type), intent(in) :: a type(psb_desc_type), intent(in) :: desc_a integer, intent(inout) :: ilaggr(:), nlaggr(:) - type(mld_s_onelev_type), intent(inout), target :: p + type(mld_sml_parms), intent(inout) :: parms + type(psb_sspmat_type), intent(out) :: ac,op_prol,op_restr integer, intent(out) :: info end subroutine mld_saggrmat_var_asb end interface diff --git a/mlprec/mld_z_inner_mod.f90 b/mlprec/mld_z_inner_mod.f90 index c25f8683..4b630ae0 100644 --- a/mlprec/mld_z_inner_mod.f90 +++ b/mlprec/mld_z_inner_mod.f90 @@ -146,13 +146,14 @@ module mld_z_inner_mod abstract interface - subroutine mld_zaggrmat_var_asb(a,desc_a,ilaggr,nlaggr,p,info) + subroutine mld_zaggrmat_var_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_restr,info) use psb_base_mod, only : psb_zspmat_type, psb_desc_type, psb_dpk_ - use mld_z_prec_type, only : mld_z_onelev_type + use mld_z_prec_type, only : mld_z_onelev_type, mld_dml_parms type(psb_zspmat_type), intent(in) :: a type(psb_desc_type), intent(in) :: desc_a integer, intent(inout) :: ilaggr(:), nlaggr(:) - type(mld_z_onelev_type), intent(inout), target :: p + type(mld_dml_parms), intent(inout) :: parms + type(psb_zspmat_type), intent(out) :: ac,op_prol,op_restr integer, intent(out) :: info end subroutine mld_zaggrmat_var_asb end interface