diff --git a/amgprec/amg_d_decmatch_mod.f90 b/amgprec/amg_d_decmatch_mod.f90 index 2a34a599..d2d3d466 100644 --- a/amgprec/amg_d_decmatch_mod.f90 +++ b/amgprec/amg_d_decmatch_mod.f90 @@ -40,49 +40,30 @@ module amg_d_decmatch_mod use psb_base_cbind_mod interface new_Match_If - function dnew_Match_If(nr, irp, ja, val, diag, w, mate) & + function dnew_Match_If(ipar,matching,lambda,nr, irp, ja, val, diag, w, mate) & & bind(c,name="dnew_Match_If") result(res) use iso_c_binding import :: psb_c_ipk_, psb_c_lpk_, psb_c_mpk_, psb_c_epk_ implicit none integer(psb_c_ipk_) :: res - integer(psb_c_ipk_), value :: nr - + integer(psb_c_ipk_), value :: nr,ipar,matching + real(c_double), value :: lambda type(c_ptr), value :: irp, ja, mate type(c_ptr), value :: val, diag, w end function dnew_Match_If end interface new_Match_If -!!$ interface amg_i_aggr_assign -!!$ module procedure amg_i_daggr_assign -!!$ end interface amg_i_aggr_assign -!!$ interface amg_build_decmatch module procedure amg_dbuild_decmatch end interface amg_build_decmatch -!!$ -!!$ interface amg_build_ahat -!!$ module procedure amg_dbuild_ahat -!!$ end interface amg_build_ahat -!!$ -!!$ interface amg_gtranspose -!!$ module procedure amg_dgtranspose -!!$ end interface amg_gtranspose -!!$ -!!$ interface amg_htranspose -!!$ module procedure amg_dhtranspose -!!$ end interface amg_htranspose -!!$ -!!$ interface amg_PMatchBox -!!$ module procedure amg_dPMatchBox -!!$ end interface amg_PMatchBox logical, parameter, private :: print_statistics=.false. contains subroutine amg_ddecmatch_build_prol(w,a,desc_a,ilaggr,nlaggr,prol,info,& - & symmetrize,reproducible,display_inp, display_out, print_out) + & symmetrize,reproducible,display_inp, display_out, print_out, & + & parallel, matching,lambda) use psb_base_mod use psb_util_mod use iso_c_binding @@ -95,7 +76,9 @@ contains type(psb_ldspmat_type), intent(out) :: prol integer(psb_ipk_), intent(out) :: info logical, optional, intent(in) :: display_inp, display_out, reproducible - logical, optional, intent(in) :: symmetrize, print_out + logical, optional, intent(in) :: symmetrize, print_out, parallel + integer(psb_ipk_), optional, intent(in) :: matching + real(psb_dpk_), optional, intent(in) :: lambda ! ! @@ -111,7 +94,9 @@ contains integer(psb_ipk_), save :: cnt=1 character(len=256) :: aname type(psb_ld_coo_sparse_mat) :: tmpcoo - logical :: display_out_, print_out_, reproducible_ + logical :: display_out_, print_out_, reproducible_, parallel_ + integer(psb_ipk_) :: matching_ + real(psb_dpk_) :: lambda_ logical, parameter :: dump=.false., debug=.false., dump_mate=.false., & & debug_ilaggr=.false., debug_sync=.false. integer(psb_ipk_), save :: idx_bldmtc=-1, idx_phase1=-1, idx_phase2=-1, idx_phase3=-1 @@ -147,6 +132,24 @@ contains reproducible_ = .false. end if + if (present(parallel)) then + parallel_ = parallel + else + parallel_ = .true. + end if + + if (present(matching)) then + matching_ = matching + else + matching_ = 2 + end if + + if (present(lambda)) then + lambda_ = lambda + else + lambda_ = 2.0 + end if + allocate(nlaggr(0:np-1),stat=info) if (info /= 0) then return @@ -175,7 +178,7 @@ contains end if if (do_timings) call psb_toc(idx_phase1) if (do_timings) call psb_tic(idx_bldmtc) - call amg_dbuild_decmatch(w,a,desc_a,mate,info) + call amg_dbuild_decmatch(parallel_,matching_,lambda_,w,a,desc_a,mate,info) if (do_timings) call psb_toc(idx_bldmtc) if (debug) write(0,*) iam,' buildprol from buildmatching:',& & info @@ -248,41 +251,6 @@ contains wtemp(idx) = w(idx)/nrmagg end if nlpairs = nlpairs+1 -!!$ else if (idx <= nc) then -!!$ ! -!!$ ! This pair involves a non-local vertex. -!!$ ! Set wtemp, then apply a tie-breaking algorithm -!!$ wtemp(k) = w(k)/nrmagg -!!$ idxg = ilv(idx) -!!$ kg = ilv(k) -!!$ if (reproducible_) then -!!$ ! -!!$ ! Tie-break by assigning to the lowest-index process -!!$ ! so that the numbering is the same as with NP=1 -!!$ ! -!!$ if (kg < idxg) then -!!$ nlaggr(iam) = nlaggr(iam) + 1 -!!$ ilaggr(k) = nlaggr(iam) -!!$ nlpairs = nlpairs+1 -!!$ else -!!$ ilaggr(k) = -2 -!!$ end if -!!$ else -!!$ ! Use a statistically unbiased tie-breaking rule, -!!$ ! this will give an even spread. -!!$ ! Delegate to i_aggr_assign. -!!$ ! Should be a symmetric function. -!!$ ! -!!$ call desc_a%indxmap%qry_halo_owner(idx,iown,info) -!!$! !$ ip = amg_i_aggr_assign(iam, iown, kg, idxg, wk, widx, nrmagg) -!!$ if (iam == ip) then -!!$ nlaggr(iam) = nlaggr(iam) + 1 -!!$ ilaggr(k) = nlaggr(iam) -!!$ nlpairs = nlpairs+1 -!!$ else -!!$ ilaggr(k) = -2 -!!$ end if -!!$ end if else write(0,*) 'Really? mate(k) > nr? ',mate(k),nr end if @@ -486,56 +454,14 @@ contains end subroutine amg_ddecmatch_build_prol -!!$ function amg_i_daggr_assign(iam, iown, kg, idxg, wk, widx, nrmagg) & -!!$ & result(iproc) -!!$ ! -!!$ ! How to break ties? This -!!$ ! must be a symmetric function, i.e. -!!$ ! the result value iproc MUST be the same upon -!!$ ! swapping simultaneously all pairs -!!$ ! (iam,iown) (kg,idxg) (wk,widx) -!!$ ! -!!$ implicit none -!!$ integer(psb_ipk_) :: iproc -!!$ integer(psb_ipk_), intent(in) :: iam, iown -!!$ integer(psb_lpk_), intent(in) :: kg, idxg -!!$ real(psb_dpk_), intent(in) :: wk, widx, nrmagg -!!$ ! -!!$ integer(psb_lpk_) :: kg2, idxg2 -!!$ -!!$ idxg2 = mod(idxg,2) -!!$ kg2 = mod(kg,2) -!!$ -!!$ ! -!!$ ! In this particular tie-breaking rule we are ignoring WK, WIDX. -!!$ ! This should statistically entail an even spread of aggregates. -!!$ ! -!!$ if ((kg2/=0).and.(idxg2/=0)) then -!!$ ! If both global indices are odd, -!!$ ! assign to the higher number process -!!$ iproc = max(iam,iown) -!!$ -!!$ else if ((kg2==0).and.(idxg2==0)) then -!!$ ! If both global indices are even, -!!$ ! assign to the lower number process; -!!$ iproc = min(iam,iown) -!!$ else -!!$ ! If the global indices are mixed, -!!$ ! then assign to the owner of the odd index -!!$ if (kg2 /= 0) then -!!$ iproc = iam -!!$ else -!!$ iproc = iown -!!$ end if -!!$ end if -!!$ end function amg_i_daggr_assign - - - subroutine amg_dbuild_decmatch(w,a,desc_a,mate,info) + subroutine amg_dbuild_decmatch(parallel,matching,lambda,w,a,desc_a,mate,info) use psb_base_mod use psb_util_mod use iso_c_binding implicit none + logical, intent(in) :: parallel + integer(psb_ipk_), intent(in) :: matching + real(psb_dpk_), intent(in) :: lambda real(psb_dpk_), target :: w(:) type(psb_dspmat_type), intent(inout) :: a type(psb_desc_type) :: desc_a @@ -547,7 +473,7 @@ contains real(psb_dpk_) :: ph0t,ph1t,ph2t type(psb_ctxt_type) :: ictxt integer(psb_ipk_) :: iam, np - integer(psb_ipk_) :: nr, nc, nz, i, nunmatch + integer(psb_ipk_) :: nr, nc, nz, i, nunmatch, ipar integer(psb_ipk_), save :: cnt=2 logical, parameter :: debug=.false., dump_ahat=.false., debug_sync=.false. logical, parameter :: old_style=.false., sort_minp=.true. @@ -562,12 +488,17 @@ contains call a%cp_to(tcsr) call psb_realloc(nr,mate,info) diag = a%get_diag(info) + if (parallel) then + ipar = 2 + else + ipar = 1 + end if ! ! Now call matching! ! if (debug) write(0,*) iam,' buildmatching into PMatchBox:' if (do_timings) call psb_tic(idx_cmboxp) - info = dnew_Match_If(nr,c_loc(tcsr%irp),c_loc(tcsr%ja),& + info = dnew_Match_If(ipar,matching,lambda,nr,c_loc(tcsr%irp),c_loc(tcsr%ja),& & c_loc(tcsr%val),c_loc(diag),c_loc(w),c_loc(mate)) if (do_timings) call psb_toc(idx_cmboxp) if (debug) write(0,*) iam,' buildmatching from PMatchBox:', info @@ -592,691 +523,7 @@ contains 9999 continue call psb_error(ictxt) -!!$ contains -!!$ subroutine fix_order(n,ja,val,iret) -!!$ use psb_base_mod -!!$ implicit none -!!$ integer(psb_ipk_), intent(in) :: n -!!$ integer(psb_lpk_), intent(inout) :: ja(:) -!!$ real(psb_dpk_), intent(inout) :: val(:) -!!$ integer(psb_ipk_),intent(out) :: iret -!!$ integer(psb_lpk_), allocatable :: ix(:) -!!$ real(psb_dpk_), allocatable :: tmp(:) -!!$ -!!$ allocate(ix(n), tmp(n),stat=iret) -!!$ if (iret /= 0) return -!!$ call psb_msort(ja(1:n),ix=ix,dir=psb_sort_up_) -!!$ tmp(1:n) = val(ix(1:n)) -!!$ val(1:n) = tmp(1:n) -!!$ end subroutine fix_order end subroutine amg_dbuild_decmatch -!!$ subroutine amg_dbuild_ahat(w,a,ahat,desc_a,info,symmetrize) -!!$ use psb_base_mod -!!$ implicit none -!!$ real(psb_dpk_), intent(in) :: w(:) -!!$ type(psb_dspmat_type), intent(inout) :: a -!!$ type(psb_dspmat_type), intent(out) :: ahat -!!$ type(psb_desc_type) :: desc_a -!!$ integer(psb_ipk_), intent(out) :: info -!!$ logical, optional :: symmetrize -!!$ -!!$ type(psb_dspmat_type) :: atnd -!!$ type(psb_d_coo_sparse_mat) :: tcoo1, tcoo2, tcoo3 -!!$ real(psb_dpk_), allocatable :: diag(:) -!!$ integer(psb_lpk_), allocatable :: ilv(:) -!!$ logical, parameter :: debug=.false., dump=.false., dump_ahat=.false., dump_inp=.false. -!!$ logical :: symmetrize_ -!!$ real(psb_dpk_) :: aii, ajj, aij, wii, wjj, tmp1, tmp2, minabs, edgnrm -!!$ integer(psb_ipk_) :: nr, nc, nz, i, nz2, nrg, ii, jj, k, k2, ncols -!!$ integer(psb_lpk_) :: nzglob -!!$ type(psb_ctxt_type) :: ictxt -!!$ integer(psb_ipk_) :: me, np -!!$ character(len=80) :: aname -!!$ real(psb_dpk_), parameter :: eps=epsilon(1.d0) -!!$ integer(psb_ipk_), save :: idx_glbt=-1, idx_phase1=-1, idx_phase2=-1 -!!$ logical, parameter :: do_timings=.true. -!!$ logical, parameter :: debug_symmetry = .false., check_size=.false. -!!$ logical, parameter :: unroll_logtrans=.false. -!!$ -!!$ ictxt = desc_a%get_ctxt() -!!$ call psb_info(ictxt,me,np) -!!$ if (present(symmetrize)) then -!!$ symmetrize_ = symmetrize -!!$ else -!!$ symmetrize_ = .false. -!!$ end if -!!$ if ((do_timings).and.(idx_phase1==-1)) & -!!$ & idx_phase1 = psb_get_timer_idx("BLD_AHAT: phase1 ") -!!$ if ((do_timings).and.(idx_glbt==-1)) & -!!$ & idx_glbt = psb_get_timer_idx("BLD_AHAT: glob_transp") -!!$ if ((do_timings).and.(idx_phase2==-1)) & -!!$ & idx_phase2 = psb_get_timer_idx("BLD_AHAT: phase2 ") -!!$ if (do_timings) call psb_tic(idx_phase1) -!!$ if (dump_inp) then -!!$ block -!!$ type(psb_ldspmat_type) :: amglob -!!$ write(aname,'(a,i3.3,a)') 'a-bld-inp-',me,'.mtx' -!!$ call a%print(fname=aname,head='Test ') -!!$ call psb_gather(amglob,a,desc_a,info) -!!$ if (me==psb_root_) then -!!$ write(aname,'(a,i0,a)') 'a-bld-inp-g-',amglob%get_nrows(),'.mtx' -!!$ call amglob%print(fname=aname,head='Test ') -!!$ end if -!!$ end block -!!$ end if -!!$ -!!$ ! -!!$ ! Extract diagonal of A -!!$ ! -!!$ call a%cp_to(tcoo1) -!!$! !$ call a%triu(atnd,info) -!!$! !$ call atnd%mv_to(tcoo1) -!!$ nr = tcoo1%get_nrows() -!!$ nc = tcoo1%get_ncols() -!!$ nz = tcoo1%get_nzeros() -!!$ diag = a%get_diag(info) -!!$ call psb_realloc(nc,diag,info) -!!$ call psb_halo(diag,desc_a,info) -!!$ -!!$ ncols = desc_a%get_local_cols() -!!$ call psb_realloc(ncols,ilv,info) -!!$ do i=1, ncols -!!$ ilv(i) = i -!!$ end do -!!$ call desc_a%l2gip(ilv,info,owned=.false.) -!!$ nr = tcoo1%get_nrows() -!!$ nc = tcoo1%get_ncols() -!!$ nz = tcoo1%get_nzeros() -!!$ call tcoo2%allocate(nr,nc,int(1.25*nz)) -!!$ k2 = 0 -!!$ ! -!!$ ! Build the entries of \^A for matching -!!$ ! -!!$ minabs = huge(dzero) -!!$ do k = 1, nz -!!$ ii = tcoo1%ia(k) -!!$ jj = tcoo1%ja(k) -!!$ ! -!!$ ! Run over only the upper triangle, then transpose. -!!$ ! In theory, A was SPD; in practice, since we -!!$ ! are building the hierarchy with Galerkin products P^T A P, -!!$ ! A will be affected by round-off -!!$ ! -!!$ if (ilv(ii) eps) then -!!$ tcoo2%val(k2) = done - (2*done*aij*wii*wjj)/(aii*(wii**2) + ajj*(wjj**2)) -!!$ else -!!$ tcoo2%val(k2) = eps -!!$ end if -!!$ end if -!!$ !else -!!$ ! write(0,*) 'build_ahat: index error :',ii,jj,' : boundaries :',nr,nc -!!$ !end if -!!$ end do -!!$ call tcoo2%set_nzeros(k2) -!!$ -!!$ if (dump) then -!!$ block -!!$ type(psb_ldspmat_type) :: amglob -!!$ call ahat%cp_from(tcoo2) -!!$ call psb_gather(amglob,ahat,desc_a,info) -!!$ if (me==psb_root_) then -!!$ write(aname,'(a,i0,a)') 'ahat-gu-',desc_a%get_global_rows(),'.mtx' -!!$ call amglob%print(fname=aname,head='Test ') -!!$ end if -!!$ write(aname,'(a,i0,a,i0,a)') 'ahat-gu-',desc_a%get_global_rows(),'-p',me,'.mtx' -!!$ call ahat%print(fname=aname,head='Test ',iv=ilv) -!!$ write(aname,'(a,i0,a,i0,a)') 'ainp-ahat-',desc_a%get_global_rows(),'-p',me,'.mtx' -!!$ call a%print(fname=aname,head='Test ',iv=ilv) -!!$ call psb_gather(amglob,a,desc_a,info) -!!$ if (me==psb_root_) then -!!$ write(aname,'(a,i0,a)') 'ainp-ahat-g-',desc_a%get_global_rows(),'.mtx' -!!$ call amglob%print(fname=aname,head='Test ') -!!$ end if -!!$ write(0,*) 'Done build_ahat' -!!$ end block -!!$ end if -!!$ if (do_timings) call psb_toc(idx_phase1) -!!$ if (do_timings) call psb_tic(idx_glbt) -!!$ -!!$ call psb_glob_transpose(tcoo2,desc_a,info,atrans=tcoo1) -!!$ if (do_timings) call psb_toc(idx_glbt) -!!$ if (do_timings) call psb_tic(idx_phase2) -!!$ if (dump) then -!!$ block -!!$ type(psb_ldspmat_type) :: amglob -!!$ call atnd%cp_from(tcoo1) -!!$ call psb_gather(amglob,atnd,desc_a,info) -!!$ if (me==psb_root_) then -!!$ write(aname,'(a,i0,a)') 'ahat-gl-',desc_a%get_global_rows(),'.mtx' -!!$ call amglob%print(fname=aname,head='Test ') -!!$ end if -!!$ write(aname,'(a,i0,a,i0,a)') 'ahat-gl-',desc_a%get_global_rows(),'-p',me,'.mtx' -!!$ call ahat%print(fname=aname,head='Test ',iv=ilv) -!!$ write(0,*) 'Done build_ahat' -!!$ end block -!!$ end if -!!$ -!!$ nz = tcoo1%get_nzeros() -!!$ nz2 = tcoo2%get_nzeros() -!!$ call tcoo2%ensure_size(nz+nz2) -!!$ tcoo2%ia(nz2+1:nz2+nz) = tcoo1%ia(1:nz) -!!$ tcoo2%ja(nz2+1:nz2+nz) = tcoo1%ja(1:nz) -!!$ tcoo2%val(nz2+1:nz2+nz) = tcoo1%val(1:nz) -!!$ call tcoo2%set_nzeros(nz+nz2) -!!$ call tcoo2%fix(info) -!!$ call tcoo1%free() -!!$ nz = tcoo2%get_nzeros() -!!$ if (unroll_logtrans) then -!!$ minabs = huge(dzero) -!!$ do k = 1, nz -!!$ tcoo2%val(k) = abs(tcoo2%val(k)) -!!$ minabs = min(minabs,tcoo2%val(k)) -!!$ end do -!!$ else -!!$ minabs = minval(abs(tcoo2%val(1:nz))) -!!$ end if -!!$ call psb_min(ictxt,minabs) -!!$ if (minabs <= dzero) then -!!$ if (me == 0) write(0,*) me, 'Min value for log correction is <=zero! ' -!!$ minabs = done -!!$ end if -!!$ if (unroll_logtrans) then -!!$ do k = 1, nz -!!$ ! tcoo2%val has already been subject to abs() above. -!!$ tcoo2%val(k) = log(tcoo2%val(k)/(0.999*minabs)) -!!$ if (tcoo2%val(k)<0) write(0,*) me, 'Warning: negative log output!' -!!$ end do -!!$ else -!!$ tcoo2%val(1:nz) = log(abs(tcoo2%val(1:nz))/(0.999*minabs)) -!!$ if (any(tcoo2%val(1:nz)<0)) write(0,*) me, 'Warning: negative log output!' -!!$ end if -!!$ -!!$ call ahat%mv_from(tcoo2) -!!$ -!!$ if (do_timings) call psb_toc(idx_phase2) -!!$ -!!$ if (check_size) then -!!$ nzglob = ahat%get_nzeros() -!!$ call psb_sum(ictxt,nzglob) -!!$ if (me==0) write(0,*) 'From build_ahat: global nzeros ',desc_a%get_global_rows(),nzglob -!!$ end if -!!$ if (debug_symmetry) then -!!$ block -!!$ integer(psb_ipk_) :: nz1, nz2 -!!$ real(psb_dpk_) :: mxv -!!$ call ahat%cp_to(tcoo1) -!!$ call psb_glob_transpose(tcoo1,desc_a,info,atrans=tcoo2) -!!$ nz1 = tcoo1%get_nzeros() -!!$ nz2 = tcoo2%get_nzeros() -!!$ call tcoo1%reallocate(nz1+nz2) -!!$ tcoo1%ia(nz1+1:nz1+nz2) = tcoo2%ia(1:nz2) -!!$ tcoo1%ja(nz1+1:nz1+nz2) = tcoo2%ja(1:nz2) -!!$ tcoo1%val(nz1+1:nz1+nz2) = -tcoo2%val(1:nz2) -!!$ call tcoo1%set_nzeros(nz1+nz2) -!!$ call tcoo1%set_dupl(psb_dupl_add_) -!!$ call tcoo1%fix(info) -!!$ nz1 = tcoo1%get_nzeros() -!!$ mxv = maxval(abs(tcoo1%val(1:nz1))) -!!$ call psb_max(ictxt,mxv) -!!$ if (me==0) write(0,*) 'Maximum disp from symmetry:',mxv -!!$ end block -!!$ end if -!!$ -!!$ if (dump_ahat) then -!!$ block -!!$ type(psb_ldspmat_type) :: amglob -!!$ write(aname,'(a,i3.3,a)') 'ahat-hfa-l-',me,'.mtx' -!!$ call ahat%print(fname=aname,head='Test ') -!!$ call psb_gather(amglob,ahat,desc_a,info) -!!$ if (me==psb_root_) then -!!$ write(aname,'(a,i0,a)') 'ahat-hfa-g-',amglob%get_nrows(),'.mtx' -!!$ call amglob%print(fname=aname,head='Test ') -!!$ end if -!!$ write(0,*) 'Done build_ahat' -!!$ end block -!!$ end if -!!$ -!!$ end subroutine amg_dbuild_ahat -!!$ -!!$ subroutine amg_dgtranspose(ain,aout,desc_a,info) -!!$ use psb_base_mod -!!$ implicit none -!!$ type(psb_ldspmat_type), intent(in) :: ain -!!$ type(psb_ldspmat_type), intent(out) :: aout -!!$ type(psb_desc_type) :: desc_a -!!$ integer(psb_ipk_), intent(out) :: info -!!$ -!!$ ! -!!$ ! BEWARE: This routine works under the assumption -!!$ ! that the same DESC_A works for both A and A^T, which -!!$ ! essentially means that A has a symmetric pattern. -!!$ ! -!!$ type(psb_ldspmat_type) :: atmp, ahalo, aglb -!!$ type(psb_ld_coo_sparse_mat) :: tmpcoo -!!$ type(psb_ld_csr_sparse_mat) :: tmpcsr -!!$ type(psb_ctxt_type) :: ictxt -!!$ integer(psb_ipk_) :: me, np -!!$ integer(psb_lpk_) :: i, j, k, nrow, ncol -!!$ integer(psb_lpk_), allocatable :: ilv(:) -!!$ character(len=80) :: aname -!!$ logical, parameter :: debug=.false., dump=.false., debug_sync=.false. -!!$ -!!$ ictxt = desc_a%get_context() -!!$ call psb_info(ictxt,me,np) -!!$ -!!$ nrow = desc_a%get_local_rows() -!!$ ncol = desc_a%get_local_cols() -!!$ if (debug_sync) then -!!$ call psb_barrier(ictxt) -!!$ if (me == 0) write(0,*) 'Start gtranspose ' -!!$ end if -!!$ call ain%cscnv(tmpcsr,info) -!!$ -!!$ if (debug) then -!!$ ilv = [(i,i=1,ncol)] -!!$ call desc_a%l2gip(ilv,info,owned=.false.) -!!$ write(aname,'(a,i3.3,a)') 'atmp-preh-',me,'.mtx' -!!$ call ain%print(fname=aname,head='atmp before haloTest ',iv=ilv) -!!$ end if -!!$ if (dump) then -!!$ call ain%cscnv(atmp,info) -!!$ call psb_gather(aglb,atmp,desc_a,info) -!!$ if (me==psb_root_) then -!!$ write(aname,'(a,i3.3,a)') 'aglob-prehalo.mtx' -!!$ call aglb%print(fname=aname,head='Test ') -!!$ end if -!!$ end if -!!$ -!!$ !call psb_loc_to_glob(tmpcsr%ja,desc_a,info) -!!$ call atmp%mv_from(tmpcsr) -!!$ -!!$ if (debug) then -!!$ write(aname,'(a,i3.3,a)') 'tmpcsr-',me,'.mtx' -!!$ call atmp%print(fname=aname,head='tmpcsr ',iv=ilv) -!!$ !call psb_set_debug_level(9999) -!!$ end if -!!$ -!!$ ! FIXME THIS NEEDS REWORKING -!!$ if (debug) write(0,*) me,' Gtranspose into sphalo :',atmp%get_nrows(),atmp%get_ncols() -!!$ call psb_sphalo(atmp,desc_a,ahalo,info,rowscale=.true.) -!!$ if (debug) write(0,*) me,' Gtranspose from sphalo :',ahalo%get_nrows(),ahalo%get_ncols() -!!$ if (info == psb_success_) call psb_rwextd(ncol,atmp,info,b=ahalo) -!!$ -!!$ if (debug) then -!!$ write(aname,'(a,i3.3,a)') 'ahalo-',me,'.mtx' -!!$ call ahalo%print(fname=aname,head='ahalo after haloTest ',iv=ilv) -!!$ write(aname,'(a,i3.3,a)') 'atmp-h-',me,'.mtx' -!!$ call atmp%print(fname=aname,head='atmp after haloTest ',iv=ilv) -!!$ end if -!!$ -!!$ if (info == psb_success_) call ahalo%free() -!!$ -!!$ call atmp%cp_to(tmpcoo) -!!$ call tmpcoo%transp() -!!$ !call psb_glob_to_loc(tmpcoo%ia,desc_a,info,iact='I') -!!$ if (debug) write(0,*) 'Before cleanup:',tmpcoo%get_nzeros() -!!$ -!!$ j = 0 -!!$ do k=1, tmpcoo%get_nzeros() -!!$ if ((tmpcoo%ia(k) > 0).and.(tmpcoo%ja(k)>0)) then -!!$ j = j+1 -!!$ tmpcoo%ia(j) = tmpcoo%ia(k) -!!$ tmpcoo%ja(j) = tmpcoo%ja(k) -!!$ tmpcoo%val(j) = tmpcoo%val(k) -!!$ end if -!!$ end do -!!$ call tmpcoo%set_nzeros(j) -!!$ -!!$ if (debug) write(0,*) 'After cleanup:',tmpcoo%get_nzeros() -!!$ -!!$ call ahalo%mv_from(tmpcoo) -!!$ if (dump) then -!!$ call psb_gather(aglb,ahalo,desc_a,info) -!!$ if (me==psb_root_) then -!!$ write(aname,'(a,i3.3,a)') 'atran-preclip.mtx' -!!$ call aglb%print(fname=aname,head='Test ') -!!$ end if -!!$ end if -!!$ -!!$ -!!$ call ahalo%csclip(aout,info,imax=nrow) -!!$ -!!$ if (debug) write(0,*) 'After clip:',aout%get_nzeros() -!!$ -!!$ if (debug_sync) then -!!$ call psb_barrier(ictxt) -!!$ if (me == 0) write(0,*) 'End gtranspose ' -!!$ end if -!!$ !call aout%cscnv(info,type='csr') -!!$ -!!$ if (dump) then -!!$ write(aname,'(a,i3.3,a)') 'atran-',me,'.mtx' -!!$ call aout%print(fname=aname,head='atrans ',iv=ilv) -!!$ call psb_gather(aglb,aout,desc_a,info) -!!$ if (me==psb_root_) then -!!$ write(aname,'(a,i3.3,a)') 'atran.mtx' -!!$ call aglb%print(fname=aname,head='Test ') -!!$ end if -!!$ end if -!!$ -!!$ end subroutine amg_dgtranspose -!!$ -!!$ subroutine amg_dhtranspose(ain,aout,desc_a,info) -!!$ use psb_base_mod -!!$ implicit none -!!$ type(psb_ldspmat_type), intent(in) :: ain -!!$ type(psb_ldspmat_type), intent(out) :: aout -!!$ type(psb_desc_type) :: desc_a -!!$ integer(psb_ipk_), intent(out) :: info -!!$ -!!$ ! -!!$ ! BEWARE: This routine works under the assumption -!!$ ! that the same DESC_A works for both A and A^T, which -!!$ ! essentially means that A has a symmetric pattern. -!!$ ! -!!$ type(psb_ldspmat_type) :: atmp, ahalo, aglb -!!$ type(psb_ld_coo_sparse_mat) :: tmpcoo, tmpc1, tmpc2, tmpch -!!$ type(psb_ld_csr_sparse_mat) :: tmpcsr -!!$ integer(psb_ipk_) :: nz1, nz2, nzh, nz -!!$ type(psb_ctxt_type) :: ictxt -!!$ integer(psb_ipk_) :: me, np -!!$ integer(psb_lpk_) :: i, j, k, nrow, ncol, nlz -!!$ integer(psb_lpk_), allocatable :: ilv(:) -!!$ character(len=80) :: aname -!!$ logical, parameter :: debug=.false., dump=.false., debug_sync=.false. -!!$ -!!$ ictxt = desc_a%get_context() -!!$ call psb_info(ictxt,me,np) -!!$ -!!$ nrow = desc_a%get_local_rows() -!!$ ncol = desc_a%get_local_cols() -!!$ if (debug_sync) then -!!$ call psb_barrier(ictxt) -!!$ if (me == 0) write(0,*) 'Start htranspose ' -!!$ end if -!!$ call ain%cscnv(tmpcsr,info) -!!$ -!!$ if (debug) then -!!$ ilv = [(i,i=1,ncol)] -!!$ call desc_a%l2gip(ilv,info,owned=.false.) -!!$ write(aname,'(a,i3.3,a)') 'atmp-preh-',me,'.mtx' -!!$ call ain%print(fname=aname,head='atmp before haloTest ',iv=ilv) -!!$ end if -!!$ if (dump) then -!!$ call ain%cscnv(atmp,info) -!!$ call psb_gather(aglb,atmp,desc_a,info) -!!$ if (me==psb_root_) then -!!$ write(aname,'(a,i3.3,a)') 'aglob-prehalo.mtx' -!!$ call aglb%print(fname=aname,head='Test ') -!!$ end if -!!$ end if -!!$ -!!$ !call psb_loc_to_glob(tmpcsr%ja,desc_a,info) -!!$ call atmp%mv_from(tmpcsr) -!!$ -!!$ if (debug) then -!!$ write(aname,'(a,i3.3,a)') 'tmpcsr-',me,'.mtx' -!!$ call atmp%print(fname=aname,head='tmpcsr ',iv=ilv) -!!$ !call psb_set_debug_level(9999) -!!$ end if -!!$ -!!$ ! FIXME THIS NEEDS REWORKING -!!$ if (debug) write(0,*) me,' Htranspose into sphalo :',atmp%get_nrows(),atmp%get_ncols() -!!$ if (.true.) then -!!$ call psb_sphalo(atmp,desc_a,ahalo,info, outfmt='coo ') -!!$ call atmp%mv_to(tmpc1) -!!$ call ahalo%mv_to(tmpch) -!!$ nz1 = tmpc1%get_nzeros() -!!$ call psb_loc_to_glob(tmpc1%ia(1:nz1),desc_a,info,iact='I') -!!$ call psb_loc_to_glob(tmpc1%ja(1:nz1),desc_a,info,iact='I') -!!$ nzh = tmpch%get_nzeros() -!!$ call psb_loc_to_glob(tmpch%ia(1:nzh),desc_a,info,iact='I') -!!$ call psb_loc_to_glob(tmpch%ja(1:nzh),desc_a,info,iact='I') -!!$ nlz = nz1+nzh -!!$ call tmpcoo%allocate(ncol,ncol,nlz) -!!$ tmpcoo%ia(1:nz1) = tmpc1%ia(1:nz1) -!!$ tmpcoo%ja(1:nz1) = tmpc1%ja(1:nz1) -!!$ tmpcoo%val(1:nz1) = tmpc1%val(1:nz1) -!!$ tmpcoo%ia(nz1+1:nz1+nzh) = tmpch%ia(1:nzh) -!!$ tmpcoo%ja(nz1+1:nz1+nzh) = tmpch%ja(1:nzh) -!!$ tmpcoo%val(nz1+1:nz1+nzh) = tmpch%val(1:nzh) -!!$ call tmpcoo%set_nzeros(nlz) -!!$ call tmpcoo%transp() -!!$ nz = tmpcoo%get_nzeros() -!!$ call psb_glob_to_loc(tmpcoo%ia(1:nz),desc_a,info,iact='I') -!!$ call psb_glob_to_loc(tmpcoo%ja(1:nz),desc_a,info,iact='I') -!!$ if (.true.) then -!!$ call tmpcoo%clean_negidx(info) -!!$ else -!!$ j = 0 -!!$ do k=1, tmpcoo%get_nzeros() -!!$ if ((tmpcoo%ia(k) > 0).and.(tmpcoo%ja(k)>0)) then -!!$ j = j+1 -!!$ tmpcoo%ia(j) = tmpcoo%ia(k) -!!$ tmpcoo%ja(j) = tmpcoo%ja(k) -!!$ tmpcoo%val(j) = tmpcoo%val(k) -!!$ end if -!!$ end do -!!$ call tmpcoo%set_nzeros(j) -!!$ end if -!!$ call ahalo%mv_from(tmpcoo) -!!$ call ahalo%csclip(aout,info,imax=nrow) -!!$ -!!$ else -!!$ call psb_sphalo(atmp,desc_a,ahalo,info, rowscale=.true.) -!!$ if (debug) write(0,*) me,' Htranspose from sphalo :',ahalo%get_nrows(),ahalo%get_ncols() -!!$ if (info == psb_success_) call psb_rwextd(ncol,atmp,info,b=ahalo) -!!$ -!!$ if (debug) then -!!$ write(aname,'(a,i3.3,a)') 'ahalo-',me,'.mtx' -!!$ call ahalo%print(fname=aname,head='ahalo after haloTest ',iv=ilv) -!!$ write(aname,'(a,i3.3,a)') 'atmp-h-',me,'.mtx' -!!$ call atmp%print(fname=aname,head='atmp after haloTest ',iv=ilv) -!!$ end if -!!$ -!!$ if (info == psb_success_) call ahalo%free() -!!$ -!!$ call atmp%cp_to(tmpcoo) -!!$ call tmpcoo%transp() -!!$ if (debug) write(0,*) 'Before cleanup:',tmpcoo%get_nzeros() -!!$ if (.true.) then -!!$ call tmpcoo%clean_negidx(info) -!!$ else -!!$ -!!$ j = 0 -!!$ do k=1, tmpcoo%get_nzeros() -!!$ if ((tmpcoo%ia(k) > 0).and.(tmpcoo%ja(k)>0)) then -!!$ j = j+1 -!!$ tmpcoo%ia(j) = tmpcoo%ia(k) -!!$ tmpcoo%ja(j) = tmpcoo%ja(k) -!!$ tmpcoo%val(j) = tmpcoo%val(k) -!!$ end if -!!$ end do -!!$ call tmpcoo%set_nzeros(j) -!!$ end if -!!$ -!!$ if (debug) write(0,*) 'After cleanup:',tmpcoo%get_nzeros() -!!$ -!!$ call ahalo%mv_from(tmpcoo) -!!$ if (dump) then -!!$ call psb_gather(aglb,ahalo,desc_a,info) -!!$ if (me==psb_root_) then -!!$ write(aname,'(a,i3.3,a)') 'atran-preclip.mtx' -!!$ call aglb%print(fname=aname,head='Test ') -!!$ end if -!!$ end if -!!$ -!!$ -!!$ call ahalo%csclip(aout,info,imax=nrow) -!!$ end if -!!$ -!!$ if (debug) write(0,*) 'After clip:',aout%get_nzeros() -!!$ -!!$ if (debug_sync) then -!!$ call psb_barrier(ictxt) -!!$ if (me == 0) write(0,*) 'End htranspose ' -!!$ end if -!!$ !call aout%cscnv(info,type='csr') -!!$ -!!$ if (dump) then -!!$ write(aname,'(a,i3.3,a)') 'atran-',me,'.mtx' -!!$ call aout%print(fname=aname,head='atrans ',iv=ilv) -!!$ call psb_gather(aglb,aout,desc_a,info) -!!$ if (me==psb_root_) then -!!$ write(aname,'(a,i3.3,a)') 'atran.mtx' -!!$ call aglb%print(fname=aname,head='Test ') -!!$ end if -!!$ end if -!!$ -!!$ end subroutine amg_dhtranspose -!!$ -!!$ subroutine amg_dPMatchBox(nlver,nledge,verlocptr,verlocind,edgelocweight,& -!!$ & verdistance, mate, myrank, numprocs, ictxt,& -!!$ & msgindsent,msgactualsent,msgpercent,& -!!$ & ph0_time, ph1_time, ph2_time, ph1_card, ph2_card,info,display_inp) -!!$ use psb_base_mod -!!$ implicit none -!!$ type(psb_ctxt_type) :: ictxt -!!$ integer(psb_c_ipk_), value :: myrank, numprocs -!!$ integer(psb_c_lpk_), value :: nlver,nledge -!!$ integer(psb_c_lpk_) :: verlocptr(:),verlocind(:), verdistance(:) -!!$ integer(psb_c_lpk_) :: mate(:) -!!$ integer(psb_c_lpk_) :: msgindsent(*),msgactualsent(*) -!!$ real(c_double) :: ph0_time, ph1_time, ph2_time -!!$ integer(psb_c_lpk_) :: ph1_card(*),ph2_card(*) -!!$ real(c_double) :: edgelocweight(:) -!!$ real(c_double) :: msgpercent(*) -!!$ integer(psb_ipk_) :: info, me, np -!!$ integer(psb_c_mpk_) :: icomm, mrank, mnp -!!$ logical, optional :: display_inp -!!$ ! -!!$ logical, parameter :: debug=.false., debug_out=.false., debug_sync=.false., dump_input=.false. -!!$ logical :: display_ -!!$ integer(psb_lpk_) :: i,k -!!$ integer(psb_ipk_), save :: cnt=1 -!!$ -!!$ call psb_info(ictxt,me,np) -!!$ icomm = psb_get_mpi_comm(ictxt) -!!$ mrank = psb_get_mpi_rank(ictxt,me) -!!$ mnp = np -!!$ if (present(display_inp)) then -!!$ display_ = display_inp -!!$ else -!!$ display_ = .false. -!!$ end if -!!$ -!!$ verlocptr(:) = verlocptr(:) - 1 -!!$ verlocind(:) = verlocind(:) - 1 -!!$ verdistance(:) = verdistance(:) -1 -!!$ -!!$ if (dump_input) then -!!$ block -!!$ integer(psb_ipk_) :: iout=20,info,i,j,k,nr -!!$ character(len=80) :: fname -!!$ write(fname,'(a,i4.4,a,i4.4,a)') 'verlocptr-l',cnt,'-i',me,'.mtx' -!!$ open(iout,file=fname,iostat=info) -!!$ if (info == 0) then -!!$ write(iout,'(a)') '%verlocptr ' -!!$ write(iout,*) nlver -!!$ do i=1, nlver+1 -!!$ write(iout,*) verlocptr(i) -!!$ end do -!!$ close(iout) -!!$ else -!!$ write(psb_err_unit,*) 'Error: could not open ',fname,' for output' -!!$ end if -!!$ write(fname,'(a,i4.4,a,i4.4,a)') 'edges-l',cnt,'-i',me,'.mtx' -!!$ open(iout,file=fname,iostat=info) -!!$ if (info == 0) then -!!$ write(iout,'(a)') '%verlocind/edgeweight ' -!!$ write(iout,*) nlver -!!$ do i=1, nlver -!!$ do j=verlocptr(i)+1,verlocptr(i+1) -!!$ write(iout,*) i-1, verlocind(j),edgelocweight(j) -!!$ end do -!!$ end do -!!$ close(iout) -!!$ else -!!$ write(psb_err_unit,*) 'Error: could not open ',fname,' for output' -!!$ end if -!!$ if (me==0) then -!!$ write(fname,'(a,i4.4,a,i4.4,a)') 'verdistance-l',cnt,'-i',me,'.mtx' -!!$ open(iout,file=fname,iostat=info) -!!$ if (info == 0) then -!!$ write(iout,'(a)') '%verdistance' -!!$ write(iout,*) np -!!$ do i=1, np+1 -!!$ write(iout,*) verdistance(i) -!!$ end do -!!$ close(iout) -!!$ else -!!$ write(psb_err_unit,*) 'Error: could not open ',fname,' for output' -!!$ end if -!!$ end if -!!$ end block -!!$ cnt = cnt + 1 -!!$ end if -!!$ -!!$ if (debug.or.display_) then -!!$ do i=0,np-1 -!!$ if (me == i) then -!!$ write(6,*) 'Process: ',me,' : Input into matching: nlver ',nlver,' nledge ',nledge -!!$ write(6,*) 'Process: ',me,' : VERDISTANCE 0-base : ',verdistance(1:np+1) -!!$ write(6,*) 'Process: ',me,' : VERLOCPTR 0-base : ',verlocptr(1:nlver+1) -!!$ write(6,*) 'Process: ',me,' : VERLOCIND 0-base : ',verlocind(1:nledge) -!!$ write(6,*) 'Process: ',me,' : EDGELOCWEIGHT : ',edgelocweight(1:nledge) -!!$ write(6,*) 'Process: ',me,' : Initial MATE : ',mate(1:nlver) -!!$ flush(6) -!!$ end if -!!$ call psb_barrier(ictxt) -!!$ end do -!!$ end if -!!$ if (debug_sync) then -!!$ call psb_barrier(ictxt) -!!$ if (me == 0) write(0,*)' Calling Decmatch ' -!!$ end if -!!$ -!!$ call DecmatchC(nlver,nledge,verlocptr,verlocind,edgelocweight,& -!!$ & verdistance, mate, mrank, mnp, icomm,& -!!$ & msgindsent,msgactualsent,msgpercent,& -!!$ & ph0_time, ph1_time, ph2_time, ph1_card, ph2_card) -!!$ verlocptr(:) = verlocptr(:) + 1 -!!$ verlocind(:) = verlocind(:) + 1 -!!$ verdistance(:) = verdistance(:) + 1 -!!$ -!!$ if (debug_sync) then -!!$ call psb_barrier(ictxt) -!!$ if (me == 0) write(0,*)' Done Decmatch ' -!!$ end if -!!$ -!!$ if (debug_out) then -!!$ do k=0,np-1 -!!$ if (me == k) then -!!$ write(6,*) 'Process: ',me,' : from Matching (0-base): ',info -!!$ do i=1,nlver -!!$ write(6,*) '(',i,',',mate(i),')' -!!$ !mate(i) = mate(i) +1 -!!$ ! -!!$ end do -!!$ flush(6) -!!$ end if -!!$ call psb_barrier(ictxt) -!!$ end do -!!$ end if -!!$ where(mate>=0) mate = mate + 1 -!!$ -!!$ end subroutine amg_dPMatchBox - end module amg_d_decmatch_mod diff --git a/amgprec/amg_d_newmatch_aggregator_mod.F90 b/amgprec/amg_d_newmatch_aggregator_mod.F90 index 62a45f00..23bd5393 100644 --- a/amgprec/amg_d_newmatch_aggregator_mod.F90 +++ b/amgprec/amg_d_newmatch_aggregator_mod.F90 @@ -82,11 +82,15 @@ module amg_d_newmatch_aggregator_mod type(nwm_Vector) :: w_c_nxt integer(psb_ipk_) :: max_csize integer(psb_ipk_) :: max_nlevels + real(psb_dpk_) :: lambda logical :: reproducible_matching = .false. logical :: need_symmetrize = .false. logical :: unsmoothed_hierarchy = .true. + logical :: parallel_matching = .true. contains procedure, pass(ag) :: bld_tprol => amg_d_newmatch_aggregator_build_tprol + procedure, pass(ag) :: csetc => amg_d_newmatch_aggr_csetc + procedure, pass(ag) :: csetr => amg_d_newmatch_aggr_csetr 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 @@ -335,7 +339,6 @@ contains write(iout,*) 'NewMatch Aggregator' write(iout,*) ' Number of Matching sweeps: ',ag%n_sweeps write(iout,*) ' Matching algorithm : ',ag%matching_alg - !write(iout,*) ' 0: Preis 1: MC64 2: SPRAL ' write(iout,*) 'Aggregator object type: ',ag%fmt() call parms%mldescr(iout,info) @@ -399,6 +402,83 @@ contains info = 0 end subroutine d_newmatch_aggregator_update_next + + subroutine amg_d_newmatch_aggr_csetr(ag,what,val,info,idx) + + Implicit None + + ! Arguments + class(amg_d_newmatch_aggregator_type), intent(inout) :: ag + character(len=*), intent(in) :: what + real(psb_dpk_), intent(in) :: val + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_), intent(in), optional :: idx + integer(psb_ipk_) :: err_act, iwhat + character(len=20) :: name='d_newmatch_aggr_csetr' + info = psb_success_ + + ! For now we ignore IDX + + select case(psb_toupper(trim(what))) + case('NWM_LAMBDA') + ag%lambda = val + case default + ! Do nothing + end select + return + end subroutine amg_d_newmatch_aggr_csetr + + subroutine amg_d_newmatch_aggr_csetc(ag,what,val,info,idx) + + Implicit None + + ! Arguments + class(amg_d_newmatch_aggregator_type), intent(inout) :: ag + character(len=*), intent(in) :: what + character(len=*), intent(in) :: val + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_), intent(in), optional :: idx + integer(psb_ipk_) :: err_act, iwhat + character(len=20) :: name='d_newmatch_aggr_csetc' + info = psb_success_ + + ! For now we ignore IDX + + select case(psb_toupper(trim(what))) + case('NWM_PARALLEL_MATCHING') + select case(psb_toupper(trim(val))) + case('SEQUENTIAL','F','FALSE') + ag%parallel_matching = .false. + case('PARALLEL','TRUE','T') + ag%parallel_matching =.true. + end select + case('NWM_REPRODUCIBLE_MATCHING') + select case(psb_toupper(trim(val))) + case('F','FALSE') + ag%reproducible_matching = .false. + case('REPRODUCIBLE','TRUE','T') + ag%reproducible_matching =.true. + end select + case('NWM_NEED_SYMMETRIZE') + select case(psb_toupper(trim(val))) + case('FALSE','F') + ag%need_symmetrize = .false. + case('SYMMETRIZE','TRUE','T') + ag%need_symmetrize =.true. + end select + case('NWM_UNSMOOTHED_HIERARCHY') + select case(psb_toupper(trim(val))) + case('F','FALSE') + ag%unsmoothed_hierarchy = .false. + case('T','TRUE') + ag%unsmoothed_hierarchy =.true. + end select + case default + ! Do nothing + end select + return + end subroutine amg_d_newmatch_aggr_csetc + subroutine d_newmatch_aggr_cseti(ag,what,val,info,idx) Implicit None @@ -416,14 +496,14 @@ contains ! For now we ignore IDX select case(psb_toupper(trim(what))) - case('NWM_MATCH_ALG') - ag%matching_alg=val + case('NWM_MATCH_ALG','NWM_MATCHING_ALG') + ag%matching_alg = val case('NWM_SWEEPS') ag%n_sweeps=val case('NWM_MAX_CSIZE') - ag%max_csize=val + ag%max_csize = val case('NWM_MAX_NLEVELS') - ag%max_nlevels=val + ag%max_nlevels = val case('NWM_W_SIZE') call ag%bld_default_w(val) case('AGGR_SIZE') @@ -442,7 +522,7 @@ contains ! Arguments class(amg_d_newmatch_aggregator_type), intent(inout) :: ag character(len=20) :: name='d_newmatch_aggr_set_default' - ag%matching_alg = 0 + ag%matching_alg = 1 ag%n_sweeps = 1 ag%max_nlevels = 36 ag%max_csize = -1 diff --git a/amgprec/impl/aggregator/Makefile b/amgprec/impl/aggregator/Makefile index 839c4f3a..39bb57e8 100644 --- a/amgprec/impl/aggregator/Makefile +++ b/amgprec/impl/aggregator/Makefile @@ -59,14 +59,7 @@ amg_s_parmatch_spmm_bld.o \ amg_s_parmatch_spmm_bld_ov.o \ amg_s_parmatch_unsmth_bld.o \ amg_s_parmatch_smth_bld.o \ -amg_s_parmatch_spmm_bld_inner.o \ -amg_d_newmatch_aggregator_mat_asb.o \ -amg_d_newmatch_aggregator_inner_mat_asb.o \ -amg_d_newmatch_aggregator_mat_bld.o \ -amg_d_newmatch_aggregator_tprol.o \ -amg_d_newmatch_map_to_tprol.o \ -amg_d_newmatch_spmm_bld_inner.o \ -amg_d_newmatch_spmm_bld_ov.o +amg_s_parmatch_spmm_bld_inner.o MPCXXOBJS=MatchBoxPC.o \ algoDistEdgeApproxDomEdgesLinearSearchMesgBndlSmallMateC.o \ diff --git a/amgprec/impl/aggregator/amg_d_newmatch_aggregator_tprol.F90 b/amgprec/impl/aggregator/amg_d_newmatch_aggregator_tprol.F90 index 597b39f1..4d796f3d 100644 --- a/amgprec/impl/aggregator/amg_d_newmatch_aggregator_tprol.F90 +++ b/amgprec/impl/aggregator/amg_d_newmatch_aggregator_tprol.F90 @@ -122,7 +122,6 @@ subroutine amg_d_newmatch_aggregator_build_tprol(ag,parms,ag_data,& end if 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) @@ -154,7 +153,7 @@ subroutine amg_d_newmatch_aggregator_build_tprol(ag,parms,ag_data,& nr = acsr%get_nrows() if (psb_size(ag%w) < nr) call ag%bld_default_w(nr) isz = acsr%get_ncols() - + call psb_realloc(isz,ixaggr,info) if (info == psb_success_) & & allocate(acv(0:n_sweeps), desc_acv(0:n_sweeps),& @@ -172,7 +171,7 @@ subroutine amg_d_newmatch_aggregator_build_tprol(ag,parms,ag_data,& call acv(0)%mv_from(acsr) call desc_a%clone(desc_acv(0),info) end if - + nrac = desc_acv(0)%get_local_rows() ncac = desc_acv(0)%get_local_cols() if (debug) write(0,*) me,' On input to level: ',nrac, ncac @@ -232,7 +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_ddecmatch_build_prol(tmpw,acv(i-1),desc_acv(i-1),ixaggr,nxaggr,tmp_prol,info,& - & symmetrize=ag%need_symmetrize,reproducible=ag%reproducible_matching) + & symmetrize=ag%need_symmetrize,reproducible=ag%reproducible_matching,& + & parallel=ag%parallel_matching,matching=ag%matching_alg,lambda=ag%lambda) 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/samples/advanced/pdegen/amg_d_pde3d.f90 b/samples/advanced/pdegen/amg_d_pde3d.f90 index 9d3028a2..0efb58f1 100644 --- a/samples/advanced/pdegen/amg_d_pde3d.f90 +++ b/samples/advanced/pdegen/amg_d_pde3d.f90 @@ -133,6 +133,9 @@ program amg_d_pde3d character(len=16) :: aggr_ord ! ordering for aggregation: NATURAL, DEGREE character(len=16) :: aggr_filter ! filtering: FILTER, NO_FILTER real(psb_dpk_) :: mncrratio ! minimum aggregation ratio + integer(psb_ipk_) :: matching_alg ! For NEW matching 1 2 3 variant + real(psb_dpk_) :: lambda ! matching LAMBDA + real(psb_dpk_), allocatable :: athresv(:) ! smoothed aggregation threshold vector integer(psb_ipk_) :: thrvsz ! size of threshold vector real(psb_dpk_) :: athres ! smoothed aggregation threshold @@ -298,6 +301,11 @@ program amg_d_pde3d call prec%set('par_aggr_alg', p_choice%par_aggr_alg, info) call prec%set('aggr_type', p_choice%aggr_type, info) call prec%set('aggr_size', p_choice%aggr_size, info) + write(0,*) 'match variant ',p_choice%matching_alg + if (p_choice%matching_alg>0)& + & call prec%set('nwm_matching_alg', p_choice%matching_alg, info) + if (p_choice%lambda>0)& + & call prec%set('nwm_lambda', p_choice%lambda, info) call prec%set('aggr_ord', p_choice%aggr_ord, info) call prec%set('aggr_filter', p_choice%aggr_filter,info) @@ -569,6 +577,8 @@ contains call read_data(prec%aggr_ord,inp_unit) ! ordering for aggregation call read_data(prec%aggr_filter,inp_unit) ! filtering call read_data(prec%mncrratio,inp_unit) ! minimum aggregation ratio + call read_data(prec%matching_alg,inp_unit) ! matching variant + call read_data(prec%lambda,inp_unit) ! lambda call read_data(prec%thrvsz,inp_unit) ! size of aggr thresh vector if (prec%thrvsz > 0) then call psb_realloc(prec%thrvsz,prec%athresv,info) @@ -638,6 +648,10 @@ contains call psb_bcast(ctxt,prec%aggr_ord) call psb_bcast(ctxt,prec%aggr_filter) call psb_bcast(ctxt,prec%mncrratio) + call psb_bcast(ctxt,prec%matching_alg) + call psb_bcast(ctxt,prec%lambda) + + call psb_bcast(ctxt,prec%thrvsz) if (prec%thrvsz > 0) then if (iam /= psb_root_) call psb_realloc(prec%thrvsz,prec%athresv,info) diff --git a/samples/advanced/pdegen/runs/tnewmtc.inp b/samples/advanced/pdegen/runs/tnewmtc.inp new file mode 100644 index 00000000..fc6e5f8e --- /dev/null +++ b/samples/advanced/pdegen/runs/tnewmtc.inp @@ -0,0 +1,59 @@ +%%%%%%%%%%% General arguments % Lines starting with % are ignored. +CSR ! Storage format CSR COO JAD +0080 ! IDIM; domain size. Linear system size is IDIM**3 +CONST ! PDECOEFF: CONST, EXP, GAUSS Coefficients of the PDE +BICGSTAB ! Iterative method: BiCGSTAB BiCGSTABL BiCG CG CGS FCG GCR RGMRES +2 ! ISTOPC +00500 ! ITMAX +1 ! ITRACE +30 ! IRST (restart for RGMRES and BiCGSTABL) +1.d-6 ! EPS +%%%%%%%%%%% Main preconditioner choices %%%%%%%%%%%%%%%% +ML-VCYCLE-BJAC-D-BJAC ! Longer descriptive name for preconditioner (up to 20 chars) +ML ! Preconditioner type: NONE JACOBI GS FBGS BJAC AS ML +%%%%%%%%%%% First smoother (for all levels but coarsest) %%%%%%%%%%%%%%%% +FBGS ! Smoother type JACOBI FBGS GS BWGS BJAC AS. For 1-level, repeats previous. +1 ! Number of sweeps for smoother +0 ! Number of overlap layers for AS preconditioner +HALO ! AS restriction operator: NONE HALO +NONE ! AS prolongation operator: NONE SUM AVG +INVK ! Subdomain solver for BJAC/AS: JACOBI GS BGS ILU ILUT MILU MUMPS SLU UMF +LLK ! AINV variant +0 ! Fill level P for ILU(P) and ILU(T,P) +1 ! Inverse Fill level P for INVK +1.d-4 ! Threshold T for ILU(T,P) +%%%%%%%%%%% Second smoother, always ignored for non-ML %%%%%%%%%%%%%%%% +NONE ! Second (post) smoother, ignored if NONE +1 ! Number of sweeps for (post) smoother +0 ! Number of overlap layers for AS preconditioner +HALO ! AS restriction operator: NONE HALO +NONE ! AS prolongation operator: NONE SUM AVG +ILU ! Subdomain solver for BJAC/AS: JACOBI GS BGS ILU ILUT MILU MUMPS SLU UMF +LLK ! AINV variant +0 ! Fill level P for ILU(P) and ILU(T,P) +8 ! Inverse Fill level P for INVK +1.d-4 ! Threshold T for ILU(T,P) +%%%%%%%%%%% Multilevel parameters %%%%%%%%%%%%%%%% +VCYCLE ! Type of multilevel CYCLE: VCYCLE WCYCLE KCYCLE MULT ADD +1 ! Number of outer sweeps for ML +-3 ! Max Number of levels in a multilevel preconditioner; if <0, lib default +-3 ! Target coarse matrix size per process; if <0, lib default +SMOOTHED ! Type of aggregation: SMOOTHED UNSMOOTHED +NEWMTC ! Parallel aggregation: DEC, SYMDEC, COUPLED NEWMTC +MATCHBOXP ! aggregation measure SOC1, MATCHBOXP NEWMTC +8 ! Requested size of the aggregates for MATCHBOXP +NATURAL ! Ordering of aggregation NATURAL DEGREE +NOFILTER ! Filtering of matrix: FILTER NOFILTER +-1.5 ! Coarsening ratio, if < 0 use library default +2 ! MATCHING variant +8.0 ! LAMBDA +-2 ! Number of thresholds in vector, next line ignored if <= 0 +0.05 0.025 ! Thresholds +-0.0100d0 ! Smoothed aggregation threshold, ignored if < 0 +%%%%%%%%%%% Coarse level solver %%%%%%%%%%%%%%%% +BJAC ! Coarsest-level solver: MUMPS UMF SLU SLUDIST JACOBI GS BJAC +ILU ! Coarsest-level subsolver for BJAC: ILU ILUT MILU UMF MUMPS SLU +DIST ! Coarsest-level matrix distribution: DIST REPL +1 ! Coarsest-level fillin P for ILU(P) and ILU(T,P) +1.d-4 ! Coarsest-level threshold T for ILU(T,P) +1 ! Number of sweeps for JACOBI/GS/BJAC coarsest-level solver