diff --git a/amgprec/Makefile b/amgprec/Makefile index c518fe8e..6b8035b1 100644 --- a/amgprec/Makefile +++ b/amgprec/Makefile @@ -17,7 +17,7 @@ DMODOBJS=amg_d_prec_type.o \ amg_d_ainv_solver.o amg_d_base_ainv_mod.o \ amg_d_invk_solver.o amg_d_invt_solver.o amg_d_krm_solver.o \ amg_d_matchboxp_mod.o amg_d_parmatch_aggregator_mod.o \ - amg_d_newmatch_aggregator_mod.o + amg_d_newmatch_aggregator_mod.o amg_d_decmatch_mod.o SMODOBJS=amg_s_prec_type.o amg_s_ilu_fact_mod.o \ amg_s_inner_mod.o amg_s_ilu_solver.o amg_s_diag_solver.o amg_s_jac_smoother.o amg_s_as_smoother.o \ @@ -127,7 +127,8 @@ amg_d_base_aggregator_mod.o: amg_base_prec_type.o amg_d_parmatch_aggregator_mod.o amg_d_dec_aggregator_mod.o: amg_d_base_aggregator_mod.o amg_d_hybrid_aggregator_mod.o amg_d_symdec_aggregator_mod.o: amg_d_dec_aggregator_mod.o amg_d_parmatch_aggregator_mod.o: amg_d_matchboxp_mod.o -amg_d_newmatch_aggregator_mod.o: amg_d_dec_aggregator_mod.o +amg_d_newmatch_aggregator_mod.o: amg_d_base_aggregator_mod.o +amg_d_newmatch_aggregator_mod.o: amg_d_decmatch_mod.o amg_c_base_aggregator_mod.o: amg_base_prec_type.o amg_c_parmatch_aggregator_mod.o amg_c_dec_aggregator_mod.o: amg_c_base_aggregator_mod.o diff --git a/amgprec/amg_d_newmatch_aggregator_mod.F90 b/amgprec/amg_d_newmatch_aggregator_mod.F90 index 59ca4c76..e8b1833d 100644 --- a/amgprec/amg_d_newmatch_aggregator_mod.F90 +++ b/amgprec/amg_d_newmatch_aggregator_mod.F90 @@ -73,16 +73,25 @@ module amg_d_newmatch_aggregator_mod ! the W argument with its update. Hence, copy it in w_nxt ! before passing it to the matching ! + integer(psb_ipk_) :: orig_aggr_size + integer(psb_ipk_) :: jacobi_sweeps real(psb_dpk_), allocatable :: w(:), w_nxt(:) + type(psb_dspmat_type), allocatable :: prol, restr + type(psb_dspmat_type), allocatable :: ac, base_a, rwa + type(psb_desc_type), allocatable :: desc_ac, desc_ax, base_desc, rwdesc type(nwm_Vector) :: w_c_nxt integer(psb_ipk_) :: max_csize integer(psb_ipk_) :: max_nlevels + logical :: reproducible_matching = .false. + logical :: need_symmetrize = .false. + logical :: unsmoothed_hierarchy = .true. contains procedure, pass(ag) :: bld_tprol => amg_d_newmatch_aggregator_build_tprol procedure, pass(ag) :: cseti => d_newmatch_aggr_cseti procedure, pass(ag) :: default => d_newmatch_aggr_set_default procedure, pass(ag) :: mat_asb => amg_d_newmatch_aggregator_mat_asb procedure, pass(ag) :: mat_bld => amg_d_newmatch_aggregator_mat_bld + procedure, pass(ag) :: inner_mat_asb => amg_d_newmatch_aggregator_inner_mat_asb procedure, pass(ag) :: update_next => d_newmatch_aggregator_update_next procedure, pass(ag) :: bld_wnxt => d_newmatch_bld_wnxt procedure, pass(ag) :: bld_default_w => d_bld_default_w @@ -179,6 +188,24 @@ module amg_d_newmatch_aggregator_mod integer(psb_ipk_), intent(out) :: info end subroutine amg_daggrmat_unsmth_spmm_asb end interface + + interface + subroutine amg_d_newmatch_aggregator_inner_mat_asb(ag,parms,a,desc_a,& + & ac,desc_ac, op_prol,op_restr,info) + import :: amg_d_newmatch_aggregator_type, psb_desc_type, psb_dspmat_type,& + & psb_ldspmat_type, psb_dpk_, psb_ipk_, psb_lpk_, amg_dml_parms, amg_daggr_data + implicit none + class(amg_d_newmatch_aggregator_type), target, intent(inout) :: ag + type(amg_dml_parms), intent(inout) :: parms + type(psb_dspmat_type), intent(in) :: a + type(psb_desc_type), intent(in) :: desc_a + type(psb_dspmat_type), intent(inout) :: op_prol,op_restr + type(psb_dspmat_type), intent(inout) :: ac + type(psb_desc_type), intent(inout) :: desc_ac + integer(psb_ipk_), intent(out) :: info + end subroutine amg_d_newmatch_aggregator_inner_mat_asb + end interface + !!$ interface !!$ subroutine amg_d_newmatch_unsmth_spmm_bld(a,desc_a,ilaggr,nlaggr,parms,& diff --git a/amgprec/impl/aggregator/amg_d_newmatch_aggregator_tprol.f90 b/amgprec/impl/aggregator/amg_d_newmatch_aggregator_tprol.F90 similarity index 89% rename from amgprec/impl/aggregator/amg_d_newmatch_aggregator_tprol.f90 rename to amgprec/impl/aggregator/amg_d_newmatch_aggregator_tprol.F90 index fa8b70eb..c2154fa9 100644 --- a/amgprec/impl/aggregator/amg_d_newmatch_aggregator_tprol.f90 +++ b/amgprec/impl/aggregator/amg_d_newmatch_aggregator_tprol.F90 @@ -12,6 +12,7 @@ subroutine amg_d_newmatch_aggregator_build_tprol(ag,parms,ag_data,& use psb_base_mod use amg_base_prec_type use amg_d_inner_mod + use amg_d_decmatch_mod #if defined(SERIAL_MPI) use amg_d_newmatch_aggregator_mod #else @@ -55,43 +56,6 @@ subroutine amg_d_newmatch_aggregator_build_tprol(ag,parms,ag_data,& integer(psb_ipk_), save :: idx_mboxp=-1, idx_spmmbld=-1, idx_sweeps_mult=-1 logical, parameter :: dump=.false., do_timings=.true., debug=.false., & & dump_prol_restr=.false. - interface - function bootCMatch(C,match_alg,n_sweeps,max_nlevels,max_csize,w)& - & bind(c,name='bootCMatch') result(P) - use iso_c_binding - import - implicit none - type(nwm_CSRMatrix) :: C, P - type(nwm_Vector) :: w - integer(c_int) :: match_alg - integer(c_int) :: n_sweeps - integer(c_int) :: max_nlevels - integer(c_int) :: max_csize - end function bootCMatch - end interface - - interface - function amg_bootCMatch_if(C,match_alg,n_sweeps,max_nlevels,max_csize,& - & w,isz,ilaggr,valaggr, num_cols) & - & bind(c,name='amg_bootCMatch_if') result(iret) - use iso_c_binding - import - implicit none - type(nwm_CSRMatrix) :: C, P - type(nwm_Vector) :: w - integer(c_int), value :: match_alg - integer(c_int), value :: n_sweeps - integer(c_int), value :: max_nlevels - integer(c_int), value :: max_csize - integer(c_int), value :: isz - integer(c_int) :: num_cols - integer(c_int) :: ilaggr(*) - real(c_double) :: valaggr(*) - integer(c_int) :: iret - end function amg_bootCMatch_if - end interface - - name='d_newmatch_tprol' ictxt = desc_a%get_context() call psb_info(ictxt,me,np) @@ -121,16 +85,6 @@ subroutine amg_d_newmatch_aggregator_build_tprol(ag,parms,ag_data,& & amg_aggr_ord_nat_,is_legal_ml_aggr_ord) call amg_check_def(parms%aggr_thresh,'Aggr_Thresh',dzero,is_legal_d_aggr_thrs) - !write(*,*) 'Build_tprol:',acsr%get_nrows(),acsr%get_ncols() - C%num_rows = acsr%get_nrows() - C%num_cols = acsr%get_ncols() - C%num_nonzeros = acsr%get_nzeros() - C%owns_data = 0 - acsr%irp = acsr%irp - 1 - acsr%ja = acsr%ja - 1 - C%i = c_loc(acsr%irp) - C%j = c_loc(acsr%ja) - C%data = c_loc(acsr%val) #if !defined(SERIAL_MPI) match_algorithm = ag%matching_alg @@ -166,8 +120,9 @@ subroutine amg_d_newmatch_aggregator_build_tprol(ag,parms,ag_data,& if (n_sweeps /= ag%n_sweeps) then write(0,*) me,' Inconsistent N_SWEEPS ',n_sweeps,ag%n_sweeps end if -!!$ if (me==0) write(0,*) 'Matching sweeps: ',n_sweeps n_sweeps = max(1,n_sweeps) + + if (debug) write(0,*) me,' Copies, with n_sweeps: ',n_sweeps,max_csize if (ag%unsmoothed_hierarchy.and.allocated(ag%base_a)) then call ag%base_a%cp_to(acsr) @@ -276,8 +231,8 @@ subroutine amg_d_newmatch_aggregator_build_tprol(ag,parms,ag_data,& ! if (debug) write(0,*) me,' Into matchbox_build_prol ',info if (do_timings) call psb_tic(idx_mboxp) -!!$ call amg_dmatchboxp_build_prol(tmpw,acv(i-1),desc_acv(i-1),ixaggr,nxaggr,tmp_prol,info,& -!!$ & symmetrize=ag%need_symmetrize,reproducible=ag%reproducible_matching) + call amg_ddecmatch_build_prol(tmpw,acv(i-1),desc_acv(i-1),ixaggr,nxaggr,tmp_prol,info,& + & symmetrize=ag%need_symmetrize,reproducible=ag%reproducible_matching) if (do_timings) call psb_toc(idx_mboxp) if (debug) write(0,*) me,' Out from matchbox_build_prol ',info if (psb_errstatus_fatal()) write(0,*)me,trim(name),'Error fatal on exit bld_tprol',info diff --git a/amgprec/impl/aggregator/newmatch_interface.cpp b/amgprec/impl/aggregator/newmatch_interface.cpp index d02dd398..7ec204b0 100644 --- a/amgprec/impl/aggregator/newmatch_interface.cpp +++ b/amgprec/impl/aggregator/newmatch_interface.cpp @@ -1,72 +1,76 @@ - #include #include +#include "psb_base_cbind.h" +#include "MatchingAlgorithms.h" -#include "bcm.h" - -bcm_CSRMatrix bootCMatch(bcm_CSRMatrix *C, int *match_algorithm, int *n_sweeps, int *max_nlevels, int *max_csize, bcm_Vector *w); -bcm_CSRMatrix bootCMatch(bcm_CSRMatrix *C, int *match_algorithm, int *n_sweeps, int *max_nlevels, int *max_csize, bcm_Vector *w){ - bcm_Vector *w_temp; - int info; - //double *w_inp; - //w_inp=bcm_VectorData(w); - - bcm_CSRMatrix *P; - bcm_CSRMatrix *Ac; - int ftcoarse=1; - int cr_it=0, cr_relax_type=0; - double cr_relax_weight=0.0; - // Here I am building Ac but I won't use it. - Ac=bcm_CSRMatchingAgg(C, &w, &P, *match_algorithm, *n_sweeps,*max_csize, *max_nlevels, &ftcoarse, - cr_it, cr_relax_type, cr_relax_weight); - //w_inp=bcm_VectorData(w); - bcm_CSRMatrixDestroy(Ac); - return *P; -} +#ifdef __cplusplus +extern "C" { +#endif -int mld_bootCMatch_if(bcm_CSRMatrix *C, int match_algorithm, int n_sweeps, - int max_nlevels, int max_csize, bcm_Vector *w, - int isz, int ilaggr[], double valaggr[], int *num_cols){ - bcm_Vector *w_temp; - int info; - //double *w_inp; - //w_inp=bcm_VectorData(w); +psb_i_t dnew_Match_If(psb_i_t nr, psb_i_t irp[], psb_i_t ja[], + psb_d_t val[], psb_d_t diag[], + psb_d_t w[], psb_i_t mate[]) - bcm_CSRMatrix *P; - bcm_CSRMatrix *Ac; - int *irp, *ja, nr, nz, nc,i,j; - double *val; - int ftcoarse=1; - int cr_it=0, cr_relax_type=0; - double cr_relax_weight=0.0; +#ifdef __cplusplus + } +#endif + +psb_i_t dnew_Match_If(psb_i_t nr, psb_i_t irp[], psb_i_t ja[], + psb_d_t val[], psb_d_t diag[], psb_d_t w[], + psb_i_t mate[]) +{ + psb_i_t info; + psb_i_t i,j; + psb_i_t ftcoarse=1; + psb_i_t cr_it=0, cr_relax_type=0; + psb_d_t cr_relax_weight=0.0; - // Sanity checks - nr = bcm_CSRMatrixNumRows(C); - nc = bcm_VectorSize(w); -// fprintf(stderr,"Sanity check: %d %d \n",nr,nc); + vector s; + vector t; + vector weights; + vector mateNode; + NODE_T u,v; + VAL_T weight; + psb_i_t preprocess = atoi(argv[2]); + psb_i_t romaInput = atoi(argv[3]); + VAL_T lambda = atof(argv[4]); + psb_d_t aii, ajj, aij, wii, wjj, tmp1, tmp2, minabs, edgnrm; + psb_i_t nt; + psb_d_t timeDiff; + MatchStat pstat; - // Here I am building Ac but I won't use it. - Ac=bcm_CSRMatchingAgg(C, &w, &P, match_algorithm, n_sweeps, max_csize, max_nlevels, &ftcoarse, - cr_it, cr_relax_type, cr_relax_weight); - irp = bcm_CSRMatrixI(P); - ja = bcm_CSRMatrixJ(P); - val = bcm_CSRMatrixData(P); - nr = bcm_CSRMatrixNumRows(P); - nc = bcm_CSRMatrixNumCols(P); - nz = bcm_CSRMatrixNumNonzeros(P); + // fprintf(stderr,"Sanity check: %d %d \n",nr,nc); + for (i=1; iu) { + // Define Ahat entry + aij = val[j-1]; + aii = diag[v]; + ajj = diag[u]; + wii = w[v]; + wjj = w[u]; + edgnrm = aii*(wii*wii) + ajj*(wjj*wjj); + if (edgnrm > eps) { + weight = 1.0 - (2*1.0*aij*wii*wjj)/(aii*(wii*wii) + ajj*(wjj*wjj)); + } else { + weight = 1e-16; + } + // + s.push_back(u); + t.push_back(v); + weights.push_back(weight); + } + } + } + runRomaWrapper(s,t,weights, nr, mateNode,preprocess,romaInput,lambda ,nt, pstat, timeDiff); if (isz < nr) return(-1); if (nz != nr) return(-2); /* loop here only makes sense when nr==nz */ for (i=0; i< nr; i++) { - for (j=irp[i]; j