From f571b7aa9ab7f1a40f5801bc4dab97d8db5a2c64 Mon Sep 17 00:00:00 2001 From: Fabio Durastante Date: Tue, 19 May 2026 18:33:50 +0200 Subject: [PATCH] Added usage of distributed MUMPS inside Richardson iteration, and usage of distributed MUMPS from snapshot --- amgprec/amg_d_richards_smoother.f90 | 3 +- .../impl/level/amg_d_base_onelev_csetc.F90 | 7 + .../impl/level/amg_d_base_onelev_cseti.F90 | 6 + .../impl/solver/amg_d_mumps_solver_apply.F90 | 147 ++++++++++++++- .../impl/solver/amg_d_mumps_solver_bld.F90 | 168 +++++++++++++++++- samples/advanced/pdegen/amg_d_pde3d.F90 | 65 ++++++- samples/advanced/pdegen/runs/amg_pde3d.inp | 17 +- 7 files changed, 396 insertions(+), 17 deletions(-) diff --git a/amgprec/amg_d_richards_smoother.f90 b/amgprec/amg_d_richards_smoother.f90 index e31ce5c7..04e57cde 100644 --- a/amgprec/amg_d_richards_smoother.f90 +++ b/amgprec/amg_d_richards_smoother.f90 @@ -342,7 +342,8 @@ contains class(amg_d_richards_smoother_type), intent(inout) :: sm integer(psb_ipk_) :: val - val = 2 + ! apply_vect uses tx, ty, tz mapped to wv(1:3) + val = 3 if (allocated(sm%sv)) val = val + sm%sv%get_wrksz() end function d_richards_smoother_get_wrksize diff --git a/amgprec/impl/level/amg_d_base_onelev_csetc.F90 b/amgprec/impl/level/amg_d_base_onelev_csetc.F90 index 493f7842..fef99e57 100644 --- a/amgprec/impl/level/amg_d_base_onelev_csetc.F90 +++ b/amgprec/impl/level/amg_d_base_onelev_csetc.F90 @@ -46,6 +46,7 @@ subroutine amg_d_base_onelev_csetc(lv,what,val,info,pos,idx) use amg_d_poly_smoother use amg_d_jac_smoother use amg_d_as_smoother + use amg_d_richards_smoother use amg_d_diag_solver use amg_d_l1_diag_solver use amg_d_jac_solver @@ -97,6 +98,7 @@ subroutine amg_d_base_onelev_csetc(lv,what,val,info,pos,idx) type(amg_d_invk_solver_type) :: amg_d_invk_solver_mold type(amg_d_invt_solver_type) :: amg_d_invt_solver_mold type(amg_d_poly_smoother_type) :: amg_d_poly_smoother_mold + type(amg_d_richards_smoother_type) :: amg_d_richards_smoother_mold #if defined(AMG_HAVE_UMF) type(amg_d_umf_solver_type) :: amg_d_umf_solver_mold #endif @@ -161,6 +163,11 @@ subroutine amg_d_base_onelev_csetc(lv,what,val,info,pos,idx) case ('POLY') call lv%set(amg_d_poly_smoother_mold,info,pos=pos) if (info == 0) call lv%set(amg_d_l1_diag_solver_mold,info,pos=pos) + + case ('RICHARDS','RICHARDSON') + call lv%set(amg_d_richards_smoother_mold,info,pos=pos) + if (info == 0) call lv%set(amg_d_ilu_solver_mold,info,pos=pos) + case ('GS','FWGS') call lv%set(amg_d_jac_smoother_mold,info,pos='pre') if (info == 0) call lv%set(amg_d_gs_solver_mold,info,pos='pre') diff --git a/amgprec/impl/level/amg_d_base_onelev_cseti.F90 b/amgprec/impl/level/amg_d_base_onelev_cseti.F90 index a23ccf99..8a74b28a 100644 --- a/amgprec/impl/level/amg_d_base_onelev_cseti.F90 +++ b/amgprec/impl/level/amg_d_base_onelev_cseti.F90 @@ -45,6 +45,7 @@ subroutine amg_d_base_onelev_cseti(lv,what,val,info,pos,idx) use amg_d_parmatch_aggregator_mod use amg_d_jac_smoother use amg_d_as_smoother + use amg_d_richards_smoother use amg_d_diag_solver use amg_d_l1_diag_solver use amg_d_ilu_solver @@ -79,6 +80,7 @@ subroutine amg_d_base_onelev_cseti(lv,what,val,info,pos,idx) type(amg_d_jac_smoother_type) :: amg_d_jac_smoother_mold type(amg_d_l1_jac_smoother_type) :: amg_d_l1_jac_smoother_mold type(amg_d_as_smoother_type) :: amg_d_as_smoother_mold + type(amg_d_richards_smoother_type) :: amg_d_richards_smoother_mold type(amg_d_diag_solver_type) :: amg_d_diag_solver_mold type(amg_d_l1_diag_solver_type) :: amg_d_l1_diag_solver_mold type(amg_d_ilu_solver_type) :: amg_d_ilu_solver_mold @@ -141,6 +143,10 @@ subroutine amg_d_base_onelev_cseti(lv,what,val,info,pos,idx) call lv%set(amg_d_as_smoother_mold,info,pos=pos) if (info == 0) call lv%set(amg_d_ilu_solver_mold,info,pos=pos) + case (amg_richardson_) + call lv%set(amg_d_richards_smoother_mold,info,pos=pos) + if (info == 0) call lv%set(amg_d_ilu_solver_mold,info,pos=pos) + case (amg_fbgs_) call lv%set(amg_d_jac_smoother_mold,info,pos='pre') if (info == 0) call lv%set(amg_d_gs_solver_mold,info,pos='pre') diff --git a/amgprec/impl/solver/amg_d_mumps_solver_apply.F90 b/amgprec/impl/solver/amg_d_mumps_solver_apply.F90 index 9375b247..9b16e911 100644 --- a/amgprec/impl/solver/amg_d_mumps_solver_apply.F90 +++ b/amgprec/impl/solver/amg_d_mumps_solver_apply.F90 @@ -59,14 +59,18 @@ subroutine d_mumps_solver_apply(alpha,sv,x,beta,y,desc_data,& integer(psb_lpk_) :: nglob integer(psb_epk_) :: eng real(psb_dpk_), allocatable :: ww(:) + real(psb_dpk_), allocatable, target :: rhs_loc(:), sol_loc(:) real(psb_dpk_), allocatable, target :: gx(:) + integer(psb_lpk_), allocatable :: gidx(:) + integer, allocatable, target :: irhs_loc(:), isol_loc(:) integer(psb_ipk_) :: i, err_act character :: trans_ character(len=20) :: name='d_mumps_solver_apply' call psb_erractionsave(err_act) -#if defined(AMG_HAVE_MUMPS) +#if defined(AMG_HAVE_MUMPS) +#if defined(AMG_MUMPS_VERSION) && (AMG_MUMPS_VERSION <= 590) info = psb_success_ trans_ = psb_toupper(trans) select case(trans_) @@ -154,6 +158,147 @@ subroutine d_mumps_solver_apply(alpha,sv,x,beta,y,desc_data,& call psb_erractionrestore(err_act) return +#else + + ! Snapshot MUMPS branch (e.g. snapshot06052026 and newer snapshots). + info = psb_success_ + trans_ = psb_toupper(trans) + select case(trans_) + case('N') + case('T') + case default + call psb_errpush(psb_err_iarg_invalid_i_,name) + goto 9999 + end select + + nglob = desc_data%get_global_rows() + n_row = desc_data%get_local_rows() + n_col = desc_data%get_local_cols() + + if (sv%ipar(1) == amg_local_solver_ ) then + gx = x + sv%id%icntl(20) = 0 + sv%id%icntl(21) = 0 + else if (sv%ipar(1) == amg_global_solver_ ) then + + if (n_col <= size(work)) then + ww = work(1:n_col) + else + allocate(ww(n_col),stat=info) + if (info /= psb_success_) then + info=psb_err_alloc_request_ + call psb_errpush(info,name,i_err=(/n_col/),& + & a_err='real(psb_dpk_)') + goto 9999 + end if + end if + allocate(rhs_loc(max(n_row,1)),sol_loc(max(n_row,1)),& + & irhs_loc(max(n_row,1)),isol_loc(max(n_row,1)),stat=info) + if (info /= psb_success_) then + info=psb_err_alloc_request_; eng = max(n_row,1) + call psb_errpush(info,name,e_err=(/eng/),& + & a_err='real(psb_dpk_)') + goto 9999 + end if + ww = 0.0_psb_dpk_ + sol_loc = 0.0_psb_dpk_ + rhs_loc = 0.0_psb_dpk_ + gidx = desc_data%get_global_indices(owned=.true.) + if (size(gidx) /= n_row) then + info = psb_err_internal_error_ + call psb_errpush(info,name,& + & a_err='Invalid local row distribution in MUMPS') + goto 9999 + end if + if ((n_row > 0) .and. ((minval(gidx) < 1_psb_lpk_) .or. & + & (maxval(gidx) > int(huge(0),psb_lpk_)))) then + info = psb_err_internal_error_ + call psb_errpush(info,name,& + & a_err='Overflow in distributed MUMPS RHS indices') + goto 9999 + end if + if (n_row > 0) then + rhs_loc(1:n_row) = x(1:n_row) + irhs_loc(1:n_row) = int(gidx(1:n_row)) + isol_loc(1:n_row) = int(gidx(1:n_row)) + end if + else + info=psb_err_internal_error_ + call psb_errpush(info,name,& + & a_err='Invalid local/global solver in MUMPS') + goto 9999 + end if + + select case(trans_) + case('N') + sv%id%icntl(9) = 1 + case('T') + sv%id%icntl(9) = 2 + case default + call psb_errpush(psb_err_internal_error_,& + & name,a_err='Invalid TRANS in subsolve') + goto 9999 + end select + + sv%id%nrhs = 1 + if (sv%ipar(1) == amg_local_solver_ ) then + sv%id%nloc_rhs = 0 + sv%id%lrhs_loc = 0 + sv%id%nsol_loc = 0 + sv%id%lsol_loc = 0 + nullify(sv%id%rhs_loc) + nullify(sv%id%irhs_loc) + nullify(sv%id%sol_loc) + nullify(sv%id%isol_loc) + sv%id%rhs => gx + else + nullify(sv%id%rhs) + sv%id%icntl(20) = 10 + sv%id%nloc_rhs = n_row + sv%id%lrhs_loc = n_row + sv%id%rhs_loc => rhs_loc + sv%id%irhs_loc => irhs_loc + sv%id%icntl(21) = 2 + sv%id%nsol_loc = n_row + sv%id%lsol_loc = n_row + sv%id%sol_loc => sol_loc + sv%id%isol_loc => isol_loc + end if + ! Snapshot path: preserve stream settings coming from build, + ! only enforce silent mode if the user did not set them explicitly. + if (sv%id%icntl(1) == 0) sv%id%icntl(1) = -1 + if (sv%id%icntl(2) == 0) sv%id%icntl(2) = -1 + if (sv%id%icntl(3) == 0) sv%id%icntl(3) = -1 + if (sv%id%icntl(4) == 0) sv%id%icntl(4) = -1 + sv%id%job = 3 + call dmumps(sv%id) + + if (sv%ipar(1) == amg_local_solver_ ) then + call psb_geaxpby(alpha,gx,beta,y,desc_data,info) + else + if (n_row > 0) ww(1:n_row) = sol_loc(1:n_row) + call psb_geaxpby(alpha,ww,beta,y,desc_data,info) + end if + + nullify(sv%id%rhs) + nullify(sv%id%rhs_loc) + nullify(sv%id%irhs_loc) + nullify(sv%id%sol_loc) + nullify(sv%id%isol_loc) + + if (info /= psb_success_) then + call psb_errpush(psb_err_internal_error_,& + & name,a_err='Error in subsolve') + goto 9999 + endif + + if (allocated(ww)) deallocate(ww) + + call psb_erractionrestore(err_act) + return + +#endif + 9999 continue call psb_erractionrestore(err_act) if (err_act == psb_act_abort_) then diff --git a/amgprec/impl/solver/amg_d_mumps_solver_bld.F90 b/amgprec/impl/solver/amg_d_mumps_solver_bld.F90 index 98919497..3733a9d5 100644 --- a/amgprec/impl/solver/amg_d_mumps_solver_bld.F90 +++ b/amgprec/impl/solver/amg_d_mumps_solver_bld.F90 @@ -69,7 +69,8 @@ subroutine d_mumps_solver_bld(a,desc_a,sv,info,b,amold,vmold,imold) integer(psb_ipk_) :: np, iam, me, i, err_act, debug_unit, debug_level character(len=20) :: name='d_mumps_solver_bld', ch_err -#if defined(AMG_HAVE_MUMPS) +#if defined(AMG_HAVE_MUMPS) +#if defined(AMG_MUMPS_VERSION) && (AMG_MUMPS_VERSION <= 590) info=psb_success_ @@ -245,6 +246,171 @@ subroutine d_mumps_solver_bld(a,desc_a,sv,info,b,amold,vmold,imold) call psb_erractionrestore(err_act) return +#else + + ! Snapshot MUMPS branch (e.g. snapshot06052026 and newer snapshots). + + info=psb_success_ + + call psb_erractionsave(err_act) + debug_unit = psb_get_debug_unit() + debug_level = psb_get_debug_level() + ctxt = desc_a%get_context() + call psb_info(ctxt, iam, np) + if (sv%ipar(1) == amg_local_solver_ ) then + call psb_init(ctxt1,np=1,basectxt=ctxt,ids=(/iam/)) + else if (sv%ipar(1) == amg_global_solver_ ) then + call psb_init(ctxt1,basectxt=ctxt) + else + info = psb_err_internal_error_ + call psb_errpush(info,name,& + & a_err='Invalid local/global solver in MUMPS') + goto 9999 + end if + icomm = ctxt1%get_mpic() + sv%local_ctxt = ctxt1 + call psb_info(ctxt1, me, npr) + npc = 1 + if (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),' start' + + if(.not.allocated(sv%id)) then + allocate(sv%id,stat=info) + if (info /= psb_success_) then + info=psb_err_alloc_dealloc_ + call psb_errpush(info,name,a_err='amg_dmumps_default') + goto 9999 + end if + end if + + ! Snapshot path: bind directly to the local communicator. + sv%id%comm = icomm + + sv%id%job = -1 + sv%id%par = 1 + if (sv%ipar(3) == 2) then + sv%id%sym = 2 + else + sv%id%sym = 0 + end if + + call dmumps(sv%id) + + if (allocated(sv%icntl)) then + do i=1,amg_mumps_icntl_size + if (allocated(sv%icntl(i)%item)) then + sv%id%icntl(i) = sv%icntl(i)%item + end if + end do + end if + if (allocated(sv%rcntl)) then + do i=1,amg_mumps_rcntl_size + if (allocated(sv%rcntl(i)%item)) sv%id%cntl(i) = sv%rcntl(i)%item + end do + end if + sv%id%icntl(5)=0 + sv%id%icntl(3)=sv%ipar(2) + + nglob = desc_a%get_global_rows() + nrow_a = a%get_nrows() + if (sv%ipar(1) == amg_local_solver_ ) then + nglobrec=desc_a%get_local_rows() + if (sv%ipar(3) == 2) then + call a%triu(c,info,jmax=nrow_a) + call c%set_symmetric() + else + call a%csclip(c,info,jmax=nrow_a) + end if + call c%cp_to(acoo) + nglob = c%get_nrows() + if (nglobrec /= nglob) then + write(*,*)'WARNING: MUMPS solver does not allow overlap in AS yet. ' + write(*,*)'A zero-overlap is used instead' + end if + else + call a%cp_to(acoo) + end if + nza = acoo%get_nzeros() + + if ((sv%ipar(1) == amg_global_solver_) .and. (nza > 0)) then +#if defined(PSB_IPK4) && defined(PSB_LPK8) + if (nglob > huge(1_psb_ipk_)) then + write(0,*) iam,' ',trim(name),': Error: overflow of local indices ' + info=psb_err_internal_error_ + call psb_errpush(info,name) + goto 9999 + end if + + gia = acoo%ia(1:nza) + gja = acoo%ja(1:nza) + call psb_loc_to_glob(gia(1:nza), desc_a, info, iact='I') + call psb_loc_to_glob(gja(1:nza), desc_a, info, iact='I') + acoo%ia(1:nza) = gia(1:nza) + acoo%ja(1:nza) = gja(1:nza) +#else + call psb_loc_to_glob(acoo%ja(1:nza), desc_a, info, iact='I') + call psb_loc_to_glob(acoo%ia(1:nza), desc_a, info, iact='I') +#endif + + if (sv%ipar(3) == 2 ) then + block + integer(psb_ipk_) :: j,nz + nz = 0 + do j=1,nza + if (acoo%ja(j) >= acoo%ia(j)) then + nz = nz + 1 + acoo%ia(nz) = acoo%ia(j) + acoo%ja(nz) = acoo%ja(j) + acoo%val(nz) = acoo%val(j) + end if + end do + call acoo%set_nzeros(nz) + call acoo%set_triangle() + call acoo%set_upper() + call acoo%set_symmetric() + nza = nz + end block + end if + end if + + ! Snapshot input mode: distributed assembled matrix. + ! Each process contributes its own non-overlapping local triplets. + sv%id%irn_loc => acoo%ia + sv%id%jcn_loc => acoo%ja + sv%id%a_loc => acoo%val + sv%id%icntl(18) = 3 + sv%id%n = nglob + sv%id%nnz_loc = acoo%get_nzeros() + sv%id%nnz = acoo%get_nzeros() + sv%id%job = 4 + if (sv%ipar(1) == amg_global_solver_ ) then + call psb_sum(ctxt,sv%id%nnz) + end if + + call dmumps(sv%id) + info = sv%id%infog(1) + if (info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='amg_dmumps_fact ' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + + nullify(sv%id%irn) + nullify(sv%id%jcn) + nullify(sv%id%a) + + call acoo%free() + sv%built=.true. + + if (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) iam,' ',trim(name),' end' + + call psb_erractionrestore(err_act) + return + +#endif + 9999 continue call psb_erractionrestore(err_act) if (err_act == psb_act_abort_) then diff --git a/samples/advanced/pdegen/amg_d_pde3d.F90 b/samples/advanced/pdegen/amg_d_pde3d.F90 index 99590a72..ea92e36e 100644 --- a/samples/advanced/pdegen/amg_d_pde3d.F90 +++ b/samples/advanced/pdegen/amg_d_pde3d.F90 @@ -160,6 +160,11 @@ program amg_d_pde3d integer(psb_ipk_) :: fill ! fill-in for incomplete LU factorization integer(psb_ipk_) :: invfill ! Inverse fill-in for INVK real(psb_dpk_) :: thr ! threshold for ILUT factorization + integer(psb_ipk_) :: mumps_blr_icntl35 ! BLR activation/format (ICNTL(35)) + integer(psb_ipk_) :: mumps_blr_icntl36 ! BLR variant (ICNTL(36)) + integer(psb_ipk_) :: mumps_blr_icntl37 ! BLR block size (ICNTL(37)) + integer(psb_ipk_) :: mumps_blr_icntl38 ! BLR compression option (ICNTL(38)) + real(psb_dpk_) :: mumps_blr_cntl7 ! BLR drop parameter (CNTL(7)) ! AMG post-smoother; ignored by 1-lev preconditioner character(len=32) :: smther2 ! post-smoother type: BJAC, AS @@ -319,8 +324,14 @@ program amg_d_pde3d call prec%set('solver_sweeps', p_choice%ssweeps, info) call prec%set('poly_degree', p_choice%degree, info) call prec%set('poly_variant', p_choice%pvariant, info) - if (psb_toupper(p_choice%solve)=='MUMPS') & - & call prec%set('mumps_loc_glob','local_solver',info) + if (psb_toupper(p_choice%solve)=='MUMPS') then + call prec%set('mumps_loc_glob','local_solver',info) + call prec%set('mumps_ipar_entry',p_choice%mumps_blr_icntl35,info,idx=35_psb_ipk_) + call prec%set('mumps_ipar_entry',p_choice%mumps_blr_icntl36,info,idx=36_psb_ipk_) + call prec%set('mumps_ipar_entry',p_choice%mumps_blr_icntl37,info,idx=37_psb_ipk_) + call prec%set('mumps_ipar_entry',p_choice%mumps_blr_icntl38,info,idx=38_psb_ipk_) + call prec%set('mumps_rpar_entry',p_choice%mumps_blr_cntl7,info,idx=7_psb_ipk_) + end if call prec%set('sub_fillin', p_choice%fill, info) call prec%set('sub_iluthrs', p_choice%thr, info) @@ -331,8 +342,14 @@ program amg_d_pde3d call prec%set('sub_prol', p_choice%prol, info) call prec%set('sub_solve', p_choice%solve, info) call prec%set('solver_sweeps', p_choice%ssweeps, info) - if (psb_toupper(p_choice%solve)=='MUMPS') & - & call prec%set('mumps_loc_glob','local_solver',info) + if (psb_toupper(p_choice%solve)=='MUMPS') then + call prec%set('mumps_loc_glob','local_solver',info) + call prec%set('mumps_ipar_entry',p_choice%mumps_blr_icntl35,info,idx=35_psb_ipk_) + call prec%set('mumps_ipar_entry',p_choice%mumps_blr_icntl36,info,idx=36_psb_ipk_) + call prec%set('mumps_ipar_entry',p_choice%mumps_blr_icntl37,info,idx=37_psb_ipk_) + call prec%set('mumps_ipar_entry',p_choice%mumps_blr_icntl38,info,idx=38_psb_ipk_) + call prec%set('mumps_rpar_entry',p_choice%mumps_blr_cntl7,info,idx=7_psb_ipk_) + end if call prec%set('sub_fillin', p_choice%fill, info) call prec%set('sub_iluthrs', p_choice%thr, info) @@ -391,8 +408,19 @@ program amg_d_pde3d call prec%set('ainv_alg', p_choice%variant, info) case default call prec%set('sub_solve', p_choice%solve, info) - if (psb_toupper(p_choice%solve)=='MUMPS') & - & call prec%set('mumps_loc_glob','local_solver',info) + if (psb_toupper(p_choice%solve)=='MUMPS') then + if (psb_toupper(p_choice%smther) == 'RICHARDS' .or. & + & psb_toupper(p_choice%smther) == 'RICHARDSON') then + call prec%set('mumps_loc_glob','global_solver',info) + else + call prec%set('mumps_loc_glob','local_solver',info) + end if + call prec%set('mumps_ipar_entry',p_choice%mumps_blr_icntl35,info,idx=35_psb_ipk_) + call prec%set('mumps_ipar_entry',p_choice%mumps_blr_icntl36,info,idx=36_psb_ipk_) + call prec%set('mumps_ipar_entry',p_choice%mumps_blr_icntl37,info,idx=37_psb_ipk_) + call prec%set('mumps_ipar_entry',p_choice%mumps_blr_icntl38,info,idx=38_psb_ipk_) + call prec%set('mumps_rpar_entry',p_choice%mumps_blr_cntl7,info,idx=7_psb_ipk_) + end if end select call prec%set('solver_sweeps', p_choice%ssweeps, info) call prec%set('sub_fillin', p_choice%fill, info) @@ -428,8 +456,19 @@ program amg_d_pde3d call prec%set('ainv_alg', p_choice%variant2, info) case default call prec%set('sub_solve', p_choice%solve2, info, pos='post') - if (psb_toupper(p_choice%solve2)=='MUMPS') & - & call prec%set('mumps_loc_glob','local_solver',info) + if (psb_toupper(p_choice%solve2)=='MUMPS') then + if (psb_toupper(p_choice%smther2) == 'RICHARDS' .or. & + & psb_toupper(p_choice%smther2) == 'RICHARDSON') then + call prec%set('mumps_loc_glob','global_solver',info) + else + call prec%set('mumps_loc_glob','local_solver',info) + end if + call prec%set('mumps_ipar_entry',p_choice%mumps_blr_icntl35,info,idx=35_psb_ipk_) + call prec%set('mumps_ipar_entry',p_choice%mumps_blr_icntl36,info,idx=36_psb_ipk_) + call prec%set('mumps_ipar_entry',p_choice%mumps_blr_icntl37,info,idx=37_psb_ipk_) + call prec%set('mumps_ipar_entry',p_choice%mumps_blr_icntl38,info,idx=38_psb_ipk_) + call prec%set('mumps_rpar_entry',p_choice%mumps_blr_cntl7,info,idx=7_psb_ipk_) + end if end select call prec%set('solver_sweeps', p_choice%ssweeps2, info,pos='post') call prec%set('sub_fillin', p_choice%fill2, info,pos='post') @@ -674,6 +713,11 @@ contains call read_data(prec%fill,inp_unit) ! fill-in for incomplete LU call read_data(prec%invfill,inp_unit) !Inverse fill-in for INVK call read_data(prec%thr,inp_unit) ! threshold for ILUT + call read_data(prec%mumps_blr_icntl35,inp_unit) ! MUMPS ICNTL(35) + call read_data(prec%mumps_blr_icntl36,inp_unit) ! MUMPS ICNTL(36) + call read_data(prec%mumps_blr_icntl37,inp_unit) ! MUMPS ICNTL(37) + call read_data(prec%mumps_blr_icntl38,inp_unit) ! MUMPS ICNTL(38) + call read_data(prec%mumps_blr_cntl7,inp_unit) ! MUMPS CNTL(7) ! Second smoother/ AMG post-smoother (if NONE ignored in main) call read_data(prec%smther2,inp_unit) ! smoother type call read_data(prec%jsweeps2,inp_unit) ! (post-)smoother sweeps @@ -775,6 +819,11 @@ contains call psb_bcast(ctxt,prec%fill) call psb_bcast(ctxt,prec%invfill) call psb_bcast(ctxt,prec%thr) + call psb_bcast(ctxt,prec%mumps_blr_icntl35) + call psb_bcast(ctxt,prec%mumps_blr_icntl36) + call psb_bcast(ctxt,prec%mumps_blr_icntl37) + call psb_bcast(ctxt,prec%mumps_blr_icntl38) + call psb_bcast(ctxt,prec%mumps_blr_cntl7) ! broadcast second (post-)smoother call psb_bcast(ctxt,prec%smther2) call psb_bcast(ctxt,prec%jsweeps2) diff --git a/samples/advanced/pdegen/runs/amg_pde3d.inp b/samples/advanced/pdegen/runs/amg_pde3d.inp index 44929bb2..8d3c180b 100644 --- a/samples/advanced/pdegen/runs/amg_pde3d.inp +++ b/samples/advanced/pdegen/runs/amg_pde3d.inp @@ -1,12 +1,12 @@ %%%%%%%%%%% General arguments % Lines starting with % are ignored. CSR ! Storage format CSR COO JAD -0150 ! IDIM; domain size. Linear system size is IDIM**3 +00070 ! IDIM; domain size. Linear system size is IDIM**3 POISSON ! PDECOEFF: POISSON, EXP, BOX, GAUSS % ! Coefficients of the PDE CG ! Iterative method: % ! BiCGSTAB BiCGSTABL BiCG CG CGS FCG GCR RGMRES RICHARDSON 2 ! ISTOPC -00500 ! ITMAX +00100 ! ITMAX 1 ! ITRACE 30 ! IRST (restart for RGMRES and BiCGSTABL) 1.d-6 ! EPS @@ -15,8 +15,8 @@ ML-VBM-VCYC-POLY-L1-D-BJAC ! Longer descriptive name for preconditioner (u ML ! Preconditioner type: NONE JACOBI GS FBGS BJAC AS ML POLY % %%%%%%%%%%% First smoother (for all levels but coarsest) %%%%%%%%%%%%%%%% -POLY ! Smoother type JACOBI FBGS GS BWGS BJAC AS POLY r 1-level, repeats previous. -6 ! Number of sweeps for smoother +RICHARDS ! Smoother type JACOBI FBGS GS BWGS BJAC AS POLY RICHARDS r 1-level, repeats previous. +1 ! Number of sweeps for smoother 6 ! degree for polynomial smoother CHEB_4_OPT ! Polynomial variant % Fields to be added for POLY @@ -28,12 +28,17 @@ POLY_RHO_EST_POWER ! POLY_RHO_ESTIMATE Currently only POLY_RH 0 ! Number of overlap layers for AS preconditioner HALO ! AS restriction operator: NONE HALO NONE ! AS prolongation operator: NONE SUM AVG -L1-JACOBI ! Subdomain solver for BJAC/AS: JACOBI GS BGS ILU ILUT MILU MUMPS SLU UMF +MUMPS ! Subdomain solver for BJAC/AS: JACOBI GS BGS ILU ILUT MILU MUMPS SLU UMF 1 ! Inner solver sweeps (GS and JACOBI) 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) +1 ! MUMPS ICNTL(35) BLR activation/format +0 ! MUMPS ICNTL(36) BLR variant +256 ! MUMPS ICNTL(37) BLR block size +333 ! MUMPS ICNTL(38) BLR compression option +1.d-2 ! MUMPS CNTL(7) BLR dropping parameter %%%%%%%%%%% Second smoother, always ignored for non-ML %%%%%%%%%%%%%%%% NONE ! Second (post) smoother, ignored if NONE 6 ! Number of sweeps for (post) smoother @@ -68,7 +73,7 @@ FILTER ! Filtering of matrix: FILTER NOFILTER -2 ! Number of thresholds in vector, next line ignored if <= 0 0.05 0.025 ! Thresholds %%%%%%%%%%% Coarse level solver %%%%%%%%%%%%%%%% -BJAC ! Coarsest-level solver: MUMPS UMF SLU SLUDIST JACOBI GS BJAC KRM +MUMPS ! Coarsest-level solver: MUMPS UMF SLU SLUDIST JACOBI GS BJAC KRM 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)