diff --git a/mlprec/impl/mld_dmlprec_aply.f90 b/mlprec/impl/mld_dmlprec_aply.f90 index 5aeabe1b..a5cceced 100644 --- a/mlprec/impl/mld_dmlprec_aply.f90 +++ b/mlprec/impl/mld_dmlprec_aply.f90 @@ -406,7 +406,7 @@ contains ! Arguments integer(psb_ipk_) :: level - type(mld_dprec_type), intent(inout) :: p + type(mld_dprec_type), target, intent(inout) :: p type(mld_mlprec_wrk_type), intent(inout) :: mlprec_wrk(:) character, intent(in) :: trans real(psb_dpk_),target :: work(:) @@ -539,7 +539,6 @@ contains end if ! This is one step of post-smoothing - if (level < nlev) then call inner_ml_aply(level+1,p,mlprec_wrk,trans,work,info) if (info /= psb_success_) then @@ -572,7 +571,7 @@ contains end if sweeps = p%precv(level)%parms%sweeps_post - call p%precv(level)%sm%apply(done,& + call p%precv(level)%sm2%apply(done,& & mlprec_wrk(level)%x2l,done,mlprec_wrk(level)%y2l,& & p%precv(level)%base_desc, trans,& & sweeps,work,info) @@ -583,8 +582,9 @@ contains end if else + ! Here at coarse level sweeps = p%precv(level)%parms%sweeps - call p%precv(level)%sm%apply(done,& + call p%precv(level)%sm2%apply(done,& & mlprec_wrk(level)%x2l,dzero,mlprec_wrk(level)%y2l,& & p%precv(level)%base_desc, trans,& & sweeps,work,info) @@ -624,7 +624,7 @@ contains else sweeps = p%precv(level)%parms%sweeps end if - call p%precv(level)%sm%apply(done,& + call p%precv(level)%sm2%apply(done,& & mlprec_wrk(level)%x2l,dzero,mlprec_wrk(level)%y2l,& & p%precv(level)%base_desc, trans,& & sweeps,work,info) @@ -830,6 +830,13 @@ contains case(mld_twoside_smooth_) + ! CHECK + if (.not.(associated(p%precv(level)%sm2,p%precv(level)%sm2a))) then + write(0,*) 'inner_ml_aply: unassociated sm2 at level ',level + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Error during restriction') + goto 9999 + end if nc2l = p%precv(level)%base_desc%get_local_cols() nr2l = p%precv(level)%base_desc%get_local_rows() allocate(mlprec_wrk(level)%ty(nc2l), mlprec_wrk(level)%tx(nc2l), stat=info) @@ -866,10 +873,19 @@ contains else sweeps = p%precv(level)%parms%sweeps end if - if (info == psb_success_) call p%precv(level)%sm%apply(done,& - & mlprec_wrk(level)%x2l,dzero,mlprec_wrk(level)%y2l,& - & p%precv(level)%base_desc, trans,& - & sweeps,work,info) + + if (trans == 'N') then + if (info == psb_success_) call p%precv(level)%sm%apply(done,& + & mlprec_wrk(level)%x2l,dzero,mlprec_wrk(level)%y2l,& + & p%precv(level)%base_desc, trans,& + & sweeps,work,info) + else + if (info == psb_success_) call p%precv(level)%sm2%apply(done,& + & mlprec_wrk(level)%x2l,dzero,mlprec_wrk(level)%y2l,& + & p%precv(level)%base_desc, trans,& + & sweeps,work,info) + end if + if (info /= psb_success_) then call psb_errpush(psb_err_internal_error_,name,& & a_err='Error during smoother_apply') @@ -930,10 +946,18 @@ contains else sweeps = p%precv(level)%parms%sweeps_pre end if - if (info == psb_success_) call p%precv(level)%sm%apply(done,& - & mlprec_wrk(level)%tx,done,mlprec_wrk(level)%y2l,& - & p%precv(level)%base_desc, trans,& - & sweeps,work,info) + if (trans == 'N') then + if (info == psb_success_) call p%precv(level)%sm2%apply(done,& + & mlprec_wrk(level)%x2l,dzero,mlprec_wrk(level)%y2l,& + & p%precv(level)%base_desc, trans,& + & sweeps,work,info) + else + if (info == psb_success_) call p%precv(level)%sm%apply(done,& + & mlprec_wrk(level)%x2l,dzero,mlprec_wrk(level)%y2l,& + & p%precv(level)%base_desc, trans,& + & sweeps,work,info) + end if + if (info /= psb_success_) then call psb_errpush(psb_err_internal_error_,name,& & a_err='Error during smoother_apply') @@ -1043,7 +1067,7 @@ subroutine mld_dmlprec_aply_vect(alpha,p,x,beta,y,desc_data,trans,work,info) level = 1 call psb_geaxpby(done,x,dzero,mlprec_wrk(level)%vx2l,p%precv(level)%base_desc,info) - call mlprec_wrk(level)%vy2l%set(dzero) + call mlprec_wrk(level)%vy2l%zero() call inner_ml_aply(level,p,mlprec_wrk,trans_,work,info) @@ -1122,7 +1146,9 @@ contains nc2l = p%precv(level)%base_desc%get_local_cols() nr2l = p%precv(level)%base_desc%get_local_rows() - + if(debug_level > 1) then + write(debug_unit,*) me,' inner_ml_aply at level ',level + end if select case(p%precv(level)%parms%ml_type) @@ -1250,7 +1276,7 @@ contains sweeps = p%precv(level)%parms%sweeps_post - call p%precv(level)%sm%apply(done,& + call p%precv(level)%sm2%apply(done,& & mlprec_wrk(level)%vx2l,done,mlprec_wrk(level)%vy2l,& & p%precv(level)%base_desc, trans,& & sweeps,work,info) @@ -1262,7 +1288,7 @@ contains else sweeps = p%precv(level)%parms%sweeps - call p%precv(level)%sm%apply(done,& + call p%precv(level)%sm2%apply(done,& & mlprec_wrk(level)%vx2l,dzero,mlprec_wrk(level)%vy2l,& & p%precv(level)%base_desc, trans,& & sweeps,work,info) @@ -1302,7 +1328,7 @@ contains else sweeps = p%precv(level)%parms%sweeps end if - call p%precv(level)%sm%apply(done,& + call p%precv(level)%sm2%apply(done,& & mlprec_wrk(level)%vx2l,dzero,mlprec_wrk(level)%vy2l,& & p%precv(level)%base_desc, trans,& & sweeps,work,info) @@ -1536,13 +1562,20 @@ contains else sweeps = p%precv(level)%parms%sweeps end if - if (info == psb_success_) call p%precv(level)%sm%apply(done,& - & mlprec_wrk(level)%vx2l,dzero,mlprec_wrk(level)%vy2l,& - & p%precv(level)%base_desc, trans,& - & sweeps,work,info) + if (trans == 'N') then + if (info == psb_success_) call p%precv(level)%sm%apply(done,& + & mlprec_wrk(level)%vx2l,dzero,mlprec_wrk(level)%vy2l,& + & p%precv(level)%base_desc, trans,& + & sweeps,work,info) + else + if (info == psb_success_) call p%precv(level)%sm2%apply(done,& + & mlprec_wrk(level)%vx2l,dzero,mlprec_wrk(level)%vy2l,& + & p%precv(level)%base_desc, trans,& + & sweeps,work,info) + end if if (info /= psb_success_) then call psb_errpush(psb_err_internal_error_,name,& - & a_err='Error during smoother_apply') + & a_err='Error during 1st smoother_apply') goto 9999 end if @@ -1605,13 +1638,20 @@ contains else sweeps = p%precv(level)%parms%sweeps_pre end if - if (info == psb_success_) call p%precv(level)%sm%apply(done,& - & mlprec_wrk(level)%vtx,done,mlprec_wrk(level)%vy2l,& - & p%precv(level)%base_desc, trans,& - & sweeps,work,info) + if (trans == 'N') then + if (info == psb_success_) call p%precv(level)%sm2%apply(done,& + & mlprec_wrk(level)%vx2l,dzero,mlprec_wrk(level)%vy2l,& + & p%precv(level)%base_desc, trans,& + & sweeps,work,info) + else + if (info == psb_success_) call p%precv(level)%sm%apply(done,& + & mlprec_wrk(level)%vx2l,dzero,mlprec_wrk(level)%vy2l,& + & p%precv(level)%base_desc, trans,& + & sweeps,work,info) + end if if (info /= psb_success_) then call psb_errpush(psb_err_internal_error_,name,& - & a_err='Error during smoother_apply') + & a_err='Error during 2nd smoother_apply') goto 9999 end if diff --git a/mlprec/impl/mld_dmlprec_bld.f90 b/mlprec/impl/mld_dmlprec_bld.f90 index 11409363..3199ded8 100644 --- a/mlprec/impl/mld_dmlprec_bld.f90 +++ b/mlprec/impl/mld_dmlprec_bld.f90 @@ -495,7 +495,7 @@ subroutine mld_dmlprec_bld(a,desc_a,p,info,amold,vmold,imold) call p%precv(i)%sm%build(p%precv(i)%base_a,p%precv(i)%base_desc,& & 'F',info,amold=amold,vmold=vmold,imold=imold) - + p%precv(i)%sm2 => p%precv(i)%sm if ((info == psb_success_).and.(i>1)) then call p%precv(i)%cnv(info,amold=amold,vmold=vmold,imold=imold) end if diff --git a/mlprec/mld_d_gs_solver.f90 b/mlprec/mld_d_gs_solver.f90 index 6599d4ec..230c5d89 100644 --- a/mlprec/mld_d_gs_solver.f90 +++ b/mlprec/mld_d_gs_solver.f90 @@ -74,6 +74,16 @@ module mld_d_gs_solver procedure, nopass :: is_iterative => d_gs_solver_is_iterative end type mld_d_gs_solver_type + type, extends(mld_d_gs_solver_type) :: mld_d_bwgs_solver_type + contains +!!$ procedure, pass(sv) :: build => mld_d_bwgs_solver_bld +!!$ procedure, pass(sv) :: apply_v => mld_d_bwgs_solver_apply_vect +!!$ procedure, pass(sv) :: apply_a => mld_d_bwgs_solver_apply + procedure, nopass :: get_fmt => d_bwgs_solver_get_fmt + procedure, pass(sv) :: descr => d_bwgs_solver_descr + end type mld_d_bwgs_solver_type + + private :: d_gs_solver_bld, d_gs_solver_apply, & & d_gs_solver_free, d_gs_solver_seti, & @@ -82,7 +92,8 @@ module mld_d_gs_solver & d_gs_solver_default, d_gs_solver_dmp, & & d_gs_solver_apply_vect, d_gs_solver_get_nzeros, & & d_gs_solver_get_fmt, d_gs_solver_check,& - & d_gs_solver_is_iterative + & d_gs_solver_is_iterative, & + & d_bwgs_solver_get_fmt, d_bwgs_solver_descr interface @@ -452,10 +463,10 @@ contains endif if (sv%eps<=dzero) then - write(iout_,*) ' Gauss-Seidel iterative solver with ',& + write(iout_,*) ' Forward Gauss-Seidel iterative solver with ',& & sv%sweeps,' sweeps' else - write(iout_,*) ' Gauss-Seidel iterative solver with tolerance',& + write(iout_,*) ' Forward Gauss-Seidel iterative solver with tolerance',& & sv%eps,' and maxit', sv%sweeps end if @@ -500,7 +511,7 @@ contains implicit none character(len=32) :: val - val = "Gauss-Seidel solver" + val = "Forward Gauss-Seidel solver" end function d_gs_solver_get_fmt @@ -515,4 +526,51 @@ contains val = .true. end function d_gs_solver_is_iterative + + subroutine d_bwgs_solver_descr(sv,info,iout,coarse) + + Implicit None + + ! Arguments + class(mld_d_bwgs_solver_type), intent(in) :: sv + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_), intent(in), optional :: iout + logical, intent(in), optional :: coarse + + ! Local variables + integer(psb_ipk_) :: err_act + character(len=20), parameter :: name='mld_d_bwgs_solver_descr' + integer(psb_ipk_) :: iout_ + + call psb_erractionsave(err_act) + info = psb_success_ + if (present(iout)) then + iout_ = iout + else + iout_ = 6 + endif + + if (sv%eps<=dzero) then + write(iout_,*) ' Backward Gauss-Seidel iterative solver with ',& + & sv%sweeps,' sweeps' + else + write(iout_,*) ' Backward Gauss-Seidel iterative solver with tolerance',& + & sv%eps,' and maxit', sv%sweeps + end if + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + return + end subroutine d_bwgs_solver_descr + + function d_bwgs_solver_get_fmt() result(val) + implicit none + character(len=32) :: val + + val = "Backward Gauss-Seidel solver" + end function d_bwgs_solver_get_fmt + + end module mld_d_gs_solver diff --git a/mlprec/mld_d_onelev_mod.f90 b/mlprec/mld_d_onelev_mod.f90 index 3401115b..9773e4fb 100644 --- a/mlprec/mld_d_onelev_mod.f90 +++ b/mlprec/mld_d_onelev_mod.f90 @@ -121,7 +121,8 @@ module mld_d_onelev_mod ! ! type mld_d_onelev_type - class(mld_d_base_smoother_type), allocatable :: sm + class(mld_d_base_smoother_type), allocatable :: sm, sm2a + class(mld_d_base_smoother_type), pointer :: sm2 type(mld_dml_parms) :: parms type(psb_dspmat_type) :: ac integer(psb_ipk_) :: ac_nz_loc, ac_nz_tot @@ -331,6 +332,8 @@ contains val = 0 if (allocated(lv%sm)) & & val = lv%sm%get_nzeros() + if (allocated(lv%sm2a)) & + & val = val + lv%sm2a%get_nzeros() end function d_base_onelev_get_nzeros function d_base_onelev_sizeof(lv) result(val) @@ -343,6 +346,7 @@ contains val = val + lv%ac%sizeof() val = val + lv%map%sizeof() if (allocated(lv%sm)) val = val + lv%sm%sizeof() + if (allocated(lv%sm2a)) val = val + lv%sm2a%sizeof() end function d_base_onelev_sizeof @@ -353,7 +357,7 @@ contains nullify(lv%base_a) nullify(lv%base_desc) - + nullify(lv%sm2) end subroutine d_base_onelev_nullify ! @@ -371,7 +375,7 @@ contains Implicit None ! Arguments - class(mld_d_onelev_type), intent(inout) :: lv + class(mld_d_onelev_type), target, intent(inout) :: lv lv%parms%sweeps = 1 lv%parms%sweeps_pre = 1 @@ -388,6 +392,12 @@ contains lv%parms%aggr_thresh = dzero if (allocated(lv%sm)) call lv%sm%default() + if (allocated(lv%sm2a)) then + call lv%sm2a%default() + lv%sm2 => lv%sm2a + else + lv%sm2 => lv%sm + end if return @@ -401,7 +411,7 @@ contains ! Arguments class(mld_d_onelev_type), target, intent(inout) :: lv - class(mld_d_onelev_type), intent(inout) :: lvout + class(mld_d_onelev_type), target, intent(inout) :: lvout integer(psb_ipk_), intent(out) :: info info = psb_success_ @@ -413,6 +423,16 @@ contains if (info==psb_success_) deallocate(lvout%sm,stat=info) end if end if + if (allocated(lv%sm2a)) then + call lv%sm%clone(lvout%sm2a,info) + lvout%sm2 => lvout%sm2a + else + if (allocated(lvout%sm2a)) then + call lvout%sm2a%free(info) + if (info==psb_success_) deallocate(lvout%sm2a,stat=info) + end if + lvout%sm2 => lvout%sm + end if if (info == psb_success_) call lv%parms%clone(lvout%parms,info) if (info == psb_success_) call lv%ac%clone(lvout%ac,info) if (info == psb_success_) call lv%desc_ac%clone(lvout%desc_ac,info) @@ -428,12 +448,21 @@ contains subroutine mld_d_onelev_move_alloc(a, b,info) use psb_base_mod implicit none - type(mld_d_onelev_type), intent(inout) :: a, b + type(mld_d_onelev_type), target, intent(inout) :: a, b integer(psb_ipk_), intent(out) :: info call b%free(info) b%parms = a%parms - call move_alloc(a%sm,b%sm) + if (associated(a%sm2,a%sm2a)) then + call move_alloc(a%sm,b%sm) + call move_alloc(a%sm2a,b%sm2a) + b%sm2 =>b%sm2a + else + call move_alloc(a%sm,b%sm) + call move_alloc(a%sm2a,b%sm2a) + b%sm2 =>b%sm + end if + if (info == psb_success_) call psb_move_alloc(a%ac,b%ac,info) if (info == psb_success_) call psb_move_alloc(a%desc_ac,b%desc_ac,info) if (info == psb_success_) call psb_move_alloc(a%map,b%map,info) diff --git a/tests/pdegen/runs/ppde.inp b/tests/pdegen/runs/ppde.inp index 7e9993c4..d8bc10ce 100644 --- a/tests/pdegen/runs/ppde.inp +++ b/tests/pdegen/runs/ppde.inp @@ -3,7 +3,7 @@ CSR ! Storage format CSR COO JAD 0100 ! IDIM; domain size is idim**3 2 ! ISTOPC 0100 ! ITMAX -1 ! ITRACE +10 ! ITRACE 30 ! IRST (restart for RGMRES and BiCGSTABL) 1.d-6 ! EPS 3L-MUL-RAS-BJAC4-ILU ! Descriptive name for preconditioner (up to 40 chars) @@ -12,7 +12,7 @@ ML ! Preconditioner NONE JACOBI BJAC AS ML HALO ! Restriction operator NONE HALO NONE ! Prolongation operator NONE SUM AVG GS ! Subdomain solver DSCALE ILU MILU ILUT UMF SLU -1 ! sweeps for GS +2 ! sweeps for GS 0 ! Level-set N for ILU(N), and P for ILUT 1.d-4 ! Threshold T for ILU(T,P) 1 ! Smoother/Jacobi sweeps