diff --git a/amgprec/Makefile b/amgprec/Makefile index 0442f5ac..79842c27 100644 --- a/amgprec/Makefile +++ b/amgprec/Makefile @@ -9,7 +9,7 @@ FINCLUDES=$(FMFLAG)$(HERE) $(FMFLAG)$(INCDIR) $(PSBLAS_INCLUDES) DMODOBJS=amg_d_prec_type.o \ amg_d_inner_mod.o amg_d_ilu_solver.o amg_d_diag_solver.o amg_d_jac_smoother.o amg_d_as_smoother.o \ - amg_d_poly_smoother.o amg_d_beta_coeff_mod.o\ + amg_d_poly_smoother.o amg_d_poly_coeff_mod.o\ amg_d_umf_solver.o amg_d_slu_solver.o amg_d_sludist_solver.o amg_d_id_solver.o\ amg_d_base_solver_mod.o amg_d_base_smoother_mod.o amg_d_onelev_mod.o \ amg_d_gs_solver.o amg_d_mumps_solver.o amg_d_jac_solver.o \ @@ -165,7 +165,7 @@ amg_d_jac_smoother.o: amg_d_diag_solver.o amg_dprecinit.o amg_dprecset.o: amg_d_diag_solver.o amg_d_ilu_solver.o \ amg_d_umf_solver.o amg_d_as_smoother.o amg_d_jac_smoother.o \ amg_d_id_solver.o amg_d_slu_solver.o amg_d_sludist_solver.o -amg_d_poly_smoother.o: amg_d_base_smoother_mod.o amg_d_beta_coeff_mod.o +amg_d_poly_smoother.o: amg_d_base_smoother_mod.o amg_d_poly_coeff_mod.o amg_s_mumps_solver.o amg_s_gs_solver.o amg_s_id_solver.o amg_s_slu_solver.o \ amg_s_diag_solver.o amg_s_ilu_solver.o amg_s_jac_solver.o: amg_s_base_solver_mod.o amg_s_prec_type.o diff --git a/amgprec/amg_base_prec_type.F90 b/amgprec/amg_base_prec_type.F90 index 63f446c2..0bd44079 100644 --- a/amgprec/amg_base_prec_type.F90 +++ b/amgprec/amg_base_prec_type.F90 @@ -321,6 +321,12 @@ module amg_base_prec_type integer(psb_ipk_), parameter :: amg_repl_mat_ = 1 integer(psb_ipk_), parameter :: amg_max_coarse_mat_ = amg_repl_mat_ ! + ! Legal values for entry: amg_poly_variant_ + ! + integer(psb_ipk_), parameter :: amg_poly_lottes_ = 0 + integer(psb_ipk_), parameter :: amg_poly_lottes_beta_ = 1 + integer(psb_ipk_), parameter :: amg_poly_new_ = 2 + ! ! Legal values for entry: amg_prec_status_ ! integer(psb_ipk_), parameter :: amg_prec_built_ = 98765 @@ -560,6 +566,12 @@ contains val = amg_as_ case('POLY') val = amg_poly_ + case('POLY_LOTTES') + val = amg_poly_lottes_ + case('POLY_LOTTES_BETA') + val = amg_poly_lottes_beta_ + case('POLY_NEW') + val = amg_poly_new_ case('A_NORMI') val = amg_max_norm_ case('USER_CHOICE') diff --git a/amgprec/amg_d_poly_coeff_mod.f90 b/amgprec/amg_d_poly_coeff_mod.f90 index 1bbeb876..8a7d8ad3 100644 --- a/amgprec/amg_d_poly_coeff_mod.f90 +++ b/amgprec/amg_d_poly_coeff_mod.f90 @@ -49,7 +49,7 @@ ! main diagonal block), so that it becomes possible to implement ! a pure Jacobi or L1-Jacobi global solver. ! -module amg_d_beta_coeff_mod +module amg_d_poly_coeff_mod use psb_base_mod real(psb_dpk_), parameter :: amg_d_beta_vect(900) = [ & @@ -513,4 +513,4 @@ module amg_d_beta_coeff_mod real(psb_dpk_), parameter :: amg_d_beta_mat(30,30)=reshape(amg_d_beta_vect,[30,30]) -end module amg_d_beta_coeff_mod +end module amg_d_poly_coeff_mod diff --git a/amgprec/amg_d_poly_smoother.f90 b/amgprec/amg_d_poly_smoother.f90 index aa1dd5d1..7e444550 100644 --- a/amgprec/amg_d_poly_smoother.f90 +++ b/amgprec/amg_d_poly_smoother.f90 @@ -51,14 +51,14 @@ ! module amg_d_poly_smoother use amg_d_base_smoother_mod - use amg_d_beta_coeff_mod + use amg_d_poly_coeff_mod type, extends(amg_d_base_smoother_type) :: amg_d_poly_smoother_type ! The local solver component is inherited from the ! parent type. ! class(amg_d_base_solver_type), allocatable :: sv ! - integer(psb_ipk_) :: pdegree + integer(psb_ipk_) :: pdegree, variant type(psb_dspmat_type), pointer :: pa => null() real(psb_dpk_), allocatable :: poly_beta(:) real(psb_dpk_) :: rho_ba @@ -319,6 +319,7 @@ contains ! sm%pdegree = 1 sm%rho_ba = dzero + sm%variant = amg_poly_lottes_ if (allocated(sm%sv)) then call sm%sv%default() diff --git a/amgprec/amg_d_prec_mod.f90 b/amgprec/amg_d_prec_mod.f90 index f86133d8..96e2f974 100644 --- a/amgprec/amg_d_prec_mod.f90 +++ b/amgprec/amg_d_prec_mod.f90 @@ -47,7 +47,6 @@ module amg_d_prec_mod use amg_d_prec_type use amg_d_jac_smoother use amg_d_as_smoother - use amg_d_poly_smoother use amg_d_id_solver use amg_d_diag_solver use amg_d_l1_diag_solver diff --git a/amgprec/impl/amg_dcprecset.F90 b/amgprec/impl/amg_dcprecset.F90 index c46f5224..83589e17 100644 --- a/amgprec/impl/amg_dcprecset.F90 +++ b/amgprec/impl/amg_dcprecset.F90 @@ -81,7 +81,6 @@ subroutine amg_dcprecseti(p,what,val,info,ilev,ilmax,pos,idx) use amg_d_prec_mod, amg_protect_name => amg_dcprecseti use amg_d_jac_smoother use amg_d_as_smoother - use amg_d_poly_smoother use amg_d_diag_solver use amg_d_l1_diag_solver use amg_d_ilu_solver @@ -126,7 +125,7 @@ subroutine amg_dcprecseti(p,what,val,info,ilev,ilmax,pos,idx) info = 3111 write(psb_err_unit,*) name,& & ': Error: uninitialized preconditioner,',& - & ' should call amg_PRECINIT' + &' should call amg_PRECINIT' return endif @@ -313,7 +312,6 @@ subroutine amg_dcprecsetc(p,what,string,info,ilev,ilmax,pos,idx) use amg_d_prec_mod, amg_protect_name => amg_dcprecsetc use amg_d_jac_smoother use amg_d_as_smoother - use amg_d_poly_smoother use amg_d_diag_solver use amg_d_l1_diag_solver use amg_d_ilu_solver @@ -404,10 +402,6 @@ subroutine amg_dcprecsetc(p,what,string,info,ilev,ilmax,pos,idx) do il=ilev_, ilmax_ call p%precv(il)%set('SMOOTHER_TYPE',amg_bjac_,info,pos=pos) end do - case ('POLY') - do il=ilev_, ilmax_ - call p%precv(il)%set('SMOOTHER_TYPE',amg_poly_,info,pos=pos) - end do case ('L1-BJAC') do il=ilev_, ilmax_ call p%precv(il)%set('SMOOTHER_TYPE',amg_l1_bjac_,info,pos=pos) diff --git a/amgprec/impl/amg_dprecinit.F90 b/amgprec/impl/amg_dprecinit.F90 index 491d8eca..8f3c0cb6 100644 --- a/amgprec/impl/amg_dprecinit.F90 +++ b/amgprec/impl/amg_dprecinit.F90 @@ -163,7 +163,7 @@ subroutine amg_dprecinit(ctxt,prec,ptype,info) allocate(prec%precv(nlev_),stat=info) allocate(amg_d_poly_smoother_type :: prec%precv(ilev_)%sm, stat=info) if (info /= psb_success_) return - allocate(amg_d_diag_solver_type :: prec%precv(ilev_)%sm%sv, stat=info) + allocate(amg_d_l1_diag_solver_type :: prec%precv(ilev_)%sm%sv, stat=info) call prec%precv(ilev_)%default() case ('L1-DIAG','L1-JACOBI','L1_DIAG','L1_JACOBI') diff --git a/amgprec/impl/level/amg_d_base_onelev_csetc.F90 b/amgprec/impl/level/amg_d_base_onelev_csetc.F90 index c8a8941d..d4d37262 100644 --- a/amgprec/impl/level/amg_d_base_onelev_csetc.F90 +++ b/amgprec/impl/level/amg_d_base_onelev_csetc.F90 @@ -207,16 +207,11 @@ subroutine amg_d_base_onelev_csetc(lv,what,val,info,pos,idx) case ('NONE','NOPREC','FACT_NONE') call lv%set(amg_d_id_solver_mold,info,pos=pos) - case ('DIAG') + case ('DIAG','JACOBI') call lv%set(amg_d_diag_solver_mold,info,pos=pos) - case ('JACOBI') - call lv%set(amg_d_jac_solver_mold,info,pos=pos) - - case ('L1-DIAG') + case ('L1-DIAG','L1-JACOBI') call lv%set(amg_d_l1_diag_solver_mold,info,pos=pos) - case ('L1-JACOBI') - call lv%set(amg_d_l1_jac_solver_mold,info,pos=pos) case ('GS','FGS','FWGS') call lv%set(amg_d_gs_solver_mold,info,pos=pos) diff --git a/amgprec/impl/level/amg_d_base_onelev_cseti.F90 b/amgprec/impl/level/amg_d_base_onelev_cseti.F90 index 195fbed7..c60ff895 100644 --- a/amgprec/impl/level/amg_d_base_onelev_cseti.F90 +++ b/amgprec/impl/level/amg_d_base_onelev_cseti.F90 @@ -45,7 +45,6 @@ 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_poly_smoother use amg_d_diag_solver use amg_d_l1_diag_solver use amg_d_ilu_solver @@ -80,7 +79,6 @@ 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_poly_smoother_type) :: amg_d_poly_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 @@ -143,10 +141,6 @@ 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_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 (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/smoother/amg_d_jac_smoother_apply_vect.f90 b/amgprec/impl/smoother/amg_d_jac_smoother_apply_vect.f90 index 7f91b358..1c206c27 100644 --- a/amgprec/impl/smoother/amg_d_jac_smoother_apply_vect.f90 +++ b/amgprec/impl/smoother/amg_d_jac_smoother_apply_vect.f90 @@ -109,7 +109,7 @@ subroutine amg_d_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,& end if endif - if(sm%checkres) then + if(.true..or.sm%checkres) then call psb_geall(r,desc_data,info) call psb_geasb(r,desc_data,info) resdenum = psb_genrm2(x,desc_data,info) @@ -159,7 +159,10 @@ subroutine amg_d_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,& & a_err='wrong init to smoother_apply') goto 9999 end select - +!!$ call psb_geaxpby(done,x,dzero,r,desc_data,info) +!!$ call psb_spmm(-done,sm%pa,ty,done,r,desc_data,info) +!!$ res = psb_genrm2(r,desc_data,info) +!!$ write(0,*) 'Jacobi smoother ',1,res do i=1, sweeps-1 ! ! Compute Y(j+1) = Y(j)+ D^(-1)*(X-A*Y(j)), @@ -173,9 +176,13 @@ subroutine amg_d_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,& call sm%sv%apply(done,tx,done,ty,desc_data,trans_,aux,wv(3:),info,init='Y') if (info /= psb_success_) exit +!!$ call psb_geaxpby(done,x,dzero,r,desc_data,info) +!!$ call psb_spmm(-done,sm%pa,ty,done,r,desc_data,info) +!!$ res = psb_genrm2(r,desc_data,info) +!!$ write(0,*) 'Jacobi smoother ',i+1,res if ( sm%checkres.and.(mod(i,sm%checkiter) == 0) ) then - call psb_geaxpby(done,x,dzero,r,r,desc_data,info) + call psb_geaxpby(done,x,dzero,r,desc_data,info) call psb_spmm(-done,sm%pa,ty,done,r,desc_data,info) res = psb_genrm2(r,desc_data,info) if( sm%printres ) then diff --git a/amgprec/impl/smoother/amg_d_jac_smoother_clone.f90 b/amgprec/impl/smoother/amg_d_jac_smoother_clone.f90 index 1a084f87..21b6da06 100644 --- a/amgprec/impl/smoother/amg_d_jac_smoother_clone.f90 +++ b/amgprec/impl/smoother/amg_d_jac_smoother_clone.f90 @@ -72,7 +72,6 @@ subroutine amg_d_jac_smoother_clone(sm,smout,info) smo%checkiter = sm%checkiter smo%printiter = sm%printiter smo%tol = sm%tol - smo%pa => sm%pa call sm%nd%clone(smo%nd,info) if ((info==psb_success_).and.(allocated(sm%sv))) then allocate(smout%sv,mold=sm%sv,stat=info) diff --git a/amgprec/impl/smoother/amg_d_poly_smoother_apply_vect.f90 b/amgprec/impl/smoother/amg_d_poly_smoother_apply_vect.f90 index 2fc3115c..044c4fb5 100644 --- a/amgprec/impl/smoother/amg_d_poly_smoother_apply_vect.f90 +++ b/amgprec/impl/smoother/amg_d_poly_smoother_apply_vect.f90 @@ -129,15 +129,64 @@ subroutine amg_d_poly_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,& call ty%zero() call tz%zero() - do i=1, sm%pdegree - call sm%sv%apply(done,r,dzero,ty,desc_data,trans_,aux,wv(5:),info,init='Y') - cz = (2*i*done-3)/(2*i*done+done) - cr = (8*i*done-4)/((2*i*done+done)*sm%rho_ba) - call psb_geaxpby(cr,ty,cz,tz,desc_data,info) - call psb_geaxpby(sm%poly_beta(i),tz,done,tx,desc_data,info) - call psb_spmm(-done,sm%pa,tz,done,r,desc_data,info,work=aux,trans=trans_) - end do + select case(sm%variant) + case(amg_poly_lottes_) + ! b == x + ! x == tx + ! + do i=1, sm%pdegree + ! B r_{k-1} + call sm%sv%apply(done,r,dzero,ty,desc_data,trans_,aux,wv(5:),info,init='Z') + cz = (2*i*done-3)/(2*i*done+done) + cr = (8*i*done-4)/((2*i*done+done)*sm%rho_ba) + ! z_k = cz z_{k-1} + cr ty = cz z_{k-1} + cr Br_{k-1} + call psb_geaxpby(cr,ty,cz,tz,desc_data,info) + ! r_k = b-Ax_k = x -A tx + call psb_geaxpby(done,tz,done,tx,desc_data,info) + if (.false.) then + call psb_geaxpby(done,x,dzero,r,desc_data,info) + call psb_spmm(-done,sm%pa,tx,done,r,desc_data,info,work=aux,trans=trans_) + else + call psb_spmm(-done,sm%pa,tz,done,r,desc_data,info,work=aux,trans=trans_) + end if +!!$ res = psb_genrm2(r,desc_data,info) +!!$ write(0,*) 'Polynomial smoother ',i,res + ! x_k = x_{k-1} + z_k + end do + case(amg_poly_lottes_beta_) + + ! b == x + ! x == tx + ! + do i=1, sm%pdegree + ! B r_{k-1} + call sm%sv%apply(done,r,dzero,ty,desc_data,trans_,aux,wv(5:),info,init='Z') + cz = (2*i*done-3)/(2*i*done+done) + cr = (8*i*done-4)/((2*i*done+done)*sm%rho_ba) + ! z_k = cz z_{k-1} + cr ty = cz z_{k-1} + cr Br_{k-1} + call psb_geaxpby(cr,ty,cz,tz,desc_data,info) + ! r_k = b-Ax_k = x -A tx + call psb_geaxpby(sm%poly_beta(i),tz,done,tx,desc_data,info) + if (.false.) then + call psb_geaxpby(done,x,dzero,r,desc_data,info) + call psb_spmm(-done,sm%pa,tx,done,r,desc_data,info,work=aux,trans=trans_) + else + call psb_spmm(-done,sm%pa,tz,done,r,desc_data,info,work=aux,trans=trans_) + end if +!!$ res = psb_genrm2(r,desc_data,info) +!!$ write(0,*) 'Polynomial smoother ',i,res + ! x_k = x_{k-1} + z_k + end do + + + case default + info=psb_err_internal_error_ + call psb_errpush(info,name,& + & a_err='wrong polynomial variant') + goto 9999 + end select + if (info == psb_success_) call psb_geaxpby(alpha,tx,beta,y,desc_data,info) if (info /= psb_success_) then @@ -339,19 +388,19 @@ subroutine amg_d_poly_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,& !!$ !!$ endif !!$ - if (.not.(4*n_col <= size(work))) then - deallocate(aux) - endif + if (.not.(4*n_col <= size(work))) then + deallocate(aux) + endif !!$ if(sm%checkres) then !!$ call psb_gefree(r,desc_data,info) !!$ end if - call psb_erractionrestore(err_act) - return + call psb_erractionrestore(err_act) + return 9999 call psb_error_handler(err_act) - return + return - end subroutine amg_d_poly_smoother_apply_vect +end subroutine amg_d_poly_smoother_apply_vect diff --git a/amgprec/impl/smoother/amg_d_poly_smoother_bld.f90 b/amgprec/impl/smoother/amg_d_poly_smoother_bld.f90 index be83a0ae..2858f9ba 100644 --- a/amgprec/impl/smoother/amg_d_poly_smoother_bld.f90 +++ b/amgprec/impl/smoother/amg_d_poly_smoother_bld.f90 @@ -39,7 +39,8 @@ subroutine amg_d_poly_smoother_bld(a,desc_a,sm,info,amold,vmold,imold) use psb_base_mod use amg_d_diag_solver - use amg_d_beta_coeff_mod + use amg_d_l1_diag_solver + use amg_d_poly_coeff_mod use amg_d_poly_smoother, amg_protect_name => amg_d_poly_smoother_bld Implicit None @@ -55,6 +56,7 @@ subroutine amg_d_poly_smoother_bld(a,desc_a,sm,info,amold,vmold,imold) type(psb_dspmat_type) :: tmpa integer(psb_ipk_) :: n_row,n_col, nrow_a, nztota, nzeros type(psb_ctxt_type) :: ctxt + real(psb_dpk_), allocatable :: da(:), dsv(:) integer(psb_ipk_) :: np, me, i, err_act, debug_unit, debug_level character(len=20) :: name='d_poly_smoother_bld', ch_err @@ -83,15 +85,54 @@ subroutine amg_d_poly_smoother_bld(a,desc_a,sm,info,amold,vmold,imold) goto 9999 end if sm%pa => a + if (.not.allocated(sm%sv)) then + info = psb_err_internal_error_ + call psb_errpush(info,name,& + & a_err='unallocated sm%sv') + goto 9999 + end if call sm%sv%build(a,desc_a,info,amold=amold,vmold=vmold) if (info /= psb_success_) then call psb_errpush(psb_err_from_subroutine_,name,& & a_err='sv%build') goto 9999 end if - - if (sm%rho_ba <= dzero) then - sm%rho_ba = psb_dspnrm1(a,desc_a,info) + + if (.true.) then + select type(ssv => sm%sv) + class is(amg_d_l1_diag_solver_type) + da = a%arwsum(info) + dsv = ssv%dv%get_vect() + sm%rho_ba = maxval(da(1:n_row)*dsv(1:n_row)) + class default + write(0,*) 'PolySmoother BUILD: only L1-Jacobi/L1-DIAG for now ',ssv%get_fmt() + sm%rho_ba = done + end select + else + block + type(psb_d_vect_type) :: tq, tz,wv(2) + real(psb_dpk_) :: qnrm, lambda + real(psb_dpk_),allocatable :: work(:) + integer(psb_ipk_) :: i, n_cols + n_cols = desc_a%get_local_cols() + allocate(work(4*n_cols)) + call psb_geasb(tz,desc_a,info,mold=vmold,scratch=.true.) + call psb_geasb(wv(1),desc_a,info,mold=vmold,scratch=.true.) + call psb_geasb(wv(2),desc_a,info,mold=vmold,scratch=.true.) + call psb_geall(tq,desc_a,info) + call tq%set(done) + call psb_geasb(tq,desc_a,info,mold=vmold) + call psb_spmm(done,a,tq,dzero,tz,desc_a,info) ! z_1 = A q_0 + do i=1,10 + call sm%sv%apply_v(done,tz,dzero,tq,desc_a,'NoTrans',work,wv,info) ! q_k = M^{-1} q_k + qnrm = psb_genrmi(tq,desc_a,info) ! qnrm = |q_k|_inf + call tq%scal((done/qnrm)) ! q_k = q_k/qnrm + call psb_spmm(done,a,tq,dzero,tz,desc_a,info) ! z_{k=1} = A q_k + lambda = psb_gedot(tq,tz,desc_a,info) ! lambda = q_k^T z_{k+1} = q_k^T A q_k + write(0,*) 'BLD: lambda estimate ',i,lambda + end do + sm%rho_ba = lambda + end block end if if (debug_level >= psb_debug_outer_) & diff --git a/amgprec/impl/smoother/amg_d_poly_smoother_csetc.f90 b/amgprec/impl/smoother/amg_d_poly_smoother_csetc.f90 index 3c1fef00..0cc786ed 100644 --- a/amgprec/impl/smoother/amg_d_poly_smoother_csetc.f90 +++ b/amgprec/impl/smoother/amg_d_poly_smoother_csetc.f90 @@ -53,8 +53,10 @@ subroutine amg_d_poly_smoother_csetc(sm,what,val,info,idx) call psb_erractionsave(err_act) select case(psb_toupper(trim(what))) - case default - call sm%amg_d_base_smoother_type%set(what,val,info,idx=idx) + case('POLY-VARIANT') + call sm%set(what,amg_stringval(val),info,idx=idx) + case default + call sm%amg_d_base_smoother_type%set(what,val,info,idx=idx) end select if (info /= psb_success_) then diff --git a/amgprec/impl/smoother/amg_d_poly_smoother_cseti.f90 b/amgprec/impl/smoother/amg_d_poly_smoother_cseti.f90 index 8bc48724..d3db3891 100644 --- a/amgprec/impl/smoother/amg_d_poly_smoother_cseti.f90 +++ b/amgprec/impl/smoother/amg_d_poly_smoother_cseti.f90 @@ -54,8 +54,16 @@ subroutine amg_d_poly_smoother_cseti(sm,what,val,info,idx) call psb_erractionsave(err_act) select case(psb_toupper(what)) - case('SMOOTHER_DEGREE') + case('POLY_DEGREE') sm%pdegree = val + case('POLY_VARIANT') + select case(val) + case(amg_poly_lottes_,amg_poly_lottes_beta_,amg_poly_new_) + sm%variant = val + case default + write(0,*) 'Invalid choice for POLY_VARIANT, defaulting to amg_poly_lottes_' + sm%variant = amg_poly_lottes_ + end select case default call sm%amg_d_base_smoother_type%set(what,val,info,idx=idx) end select diff --git a/amgprec/impl/smoother/amg_d_poly_smoother_descr.f90 b/amgprec/impl/smoother/amg_d_poly_smoother_descr.f90 index d17c5aa5..51c198b1 100644 --- a/amgprec/impl/smoother/amg_d_poly_smoother_descr.f90 +++ b/amgprec/impl/smoother/amg_d_poly_smoother_descr.f90 @@ -78,6 +78,14 @@ subroutine amg_d_poly_smoother_descr(sm,info,iout,coarse,prefix) end if write(iout_,*) trim(prefix_), ' Polynomial smoother ' + select case(sm%variant) + case(amg_poly_lottes_) + write(iout_,*) trim(prefix_), ' variant: ','POLY_LOTTES' + case(amg_poly_lottes_beta_) + write(iout_,*) trim(prefix_), ' variant: ','POLY_LOTTES_BETA' + case default + write(iout_,*) trim(prefix_), ' variant: ','UNKNOWN???' + end select write(iout_,*) trim(prefix_), ' Degree: ',sm%pdegree write(iout_,*) trim(prefix_), ' rho_ba: ',sm%rho_ba if (allocated(sm%poly_beta)) write(iout_,*) trim(prefix_), ' Coefficients: ',sm%poly_beta(1:sm%pdegree) diff --git a/amgprec/impl/solver/amg_c_ilu_solver_bld.f90 b/amgprec/impl/solver/amg_c_ilu_solver_bld.f90 index 9f6109db..a348fcea 100644 --- a/amgprec/impl/solver/amg_c_ilu_solver_bld.f90 +++ b/amgprec/impl/solver/amg_c_ilu_solver_bld.f90 @@ -99,11 +99,11 @@ subroutine amg_c_ilu_solver_bld(a,desc_a,sv,info,b,amold,vmold,imold) select case(sv%fact_type) case (amg_ilu_n_) - psb_fctype = amg_ilu_n_ + psb_fctype = psb_ilu_n_ case (amg_milu_n_) - psb_fctype = amg_milu_n_ + psb_fctype = psb_milu_n_ case (amg_ilu_t_) - psb_fctype = amg_ilu_t_ + psb_fctype = psb_ilu_t_ case default ! If we end up here, something was wrong up in the call chain. info = psb_err_input_value_invalid_i_ diff --git a/amgprec/impl/solver/amg_d_ilu_solver_bld.f90 b/amgprec/impl/solver/amg_d_ilu_solver_bld.f90 index c4ab246a..7a49e47e 100644 --- a/amgprec/impl/solver/amg_d_ilu_solver_bld.f90 +++ b/amgprec/impl/solver/amg_d_ilu_solver_bld.f90 @@ -99,11 +99,11 @@ subroutine amg_d_ilu_solver_bld(a,desc_a,sv,info,b,amold,vmold,imold) select case(sv%fact_type) case (amg_ilu_n_) - psb_fctype = amg_ilu_n_ + psb_fctype = psb_ilu_n_ case (amg_milu_n_) - psb_fctype = amg_milu_n_ + psb_fctype = psb_milu_n_ case (amg_ilu_t_) - psb_fctype = amg_ilu_t_ + psb_fctype = psb_ilu_t_ case default ! If we end up here, something was wrong up in the call chain. info = psb_err_input_value_invalid_i_ diff --git a/amgprec/impl/solver/amg_s_ilu_solver_bld.f90 b/amgprec/impl/solver/amg_s_ilu_solver_bld.f90 index 8cfadb3b..6c36bec2 100644 --- a/amgprec/impl/solver/amg_s_ilu_solver_bld.f90 +++ b/amgprec/impl/solver/amg_s_ilu_solver_bld.f90 @@ -99,11 +99,11 @@ subroutine amg_s_ilu_solver_bld(a,desc_a,sv,info,b,amold,vmold,imold) select case(sv%fact_type) case (amg_ilu_n_) - psb_fctype = amg_ilu_n_ + psb_fctype = psb_ilu_n_ case (amg_milu_n_) - psb_fctype = amg_milu_n_ + psb_fctype = psb_milu_n_ case (amg_ilu_t_) - psb_fctype = amg_ilu_t_ + psb_fctype = psb_ilu_t_ case default ! If we end up here, something was wrong up in the call chain. info = psb_err_input_value_invalid_i_ diff --git a/amgprec/impl/solver/amg_z_ilu_solver_bld.f90 b/amgprec/impl/solver/amg_z_ilu_solver_bld.f90 index 3afa82f7..36c91ad8 100644 --- a/amgprec/impl/solver/amg_z_ilu_solver_bld.f90 +++ b/amgprec/impl/solver/amg_z_ilu_solver_bld.f90 @@ -99,11 +99,11 @@ subroutine amg_z_ilu_solver_bld(a,desc_a,sv,info,b,amold,vmold,imold) select case(sv%fact_type) case (amg_ilu_n_) - psb_fctype = amg_ilu_n_ + psb_fctype = psb_ilu_n_ case (amg_milu_n_) - psb_fctype = amg_milu_n_ + psb_fctype = psb_milu_n_ case (amg_ilu_t_) - psb_fctype = amg_ilu_t_ + psb_fctype = psb_ilu_t_ case default ! If we end up here, something was wrong up in the call chain. info = psb_err_input_value_invalid_i_ diff --git a/samples/advanced/pdegen/amg_d_pde3d.F90 b/samples/advanced/pdegen/amg_d_pde3d.F90 index 410b6e01..927fb943 100644 --- a/samples/advanced/pdegen/amg_d_pde3d.F90 +++ b/samples/advanced/pdegen/amg_d_pde3d.F90 @@ -145,6 +145,7 @@ program amg_d_pde3d character(len=16) :: smther ! (pre-)smoother type: BJAC, AS integer(psb_ipk_) :: jsweeps ! (pre-)smoother / 1-lev prec. sweeps integer(psb_ipk_) :: degree ! degree for polynomial smoother + character(len=16) :: pvariant ! Polynomial smoother variant integer(psb_ipk_) :: novr ! number of overlap layers character(len=16) :: restr ! restriction over application of AS character(len=16) :: prol ! prolongation over application of AS @@ -154,12 +155,13 @@ program amg_d_pde3d character(len=16) :: variant ! AINV variant: LLK, etc 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 + real(psb_dpk_) :: thr ! threshold for ILUT factorization ! AMG post-smoother; ignored by 1-lev preconditioner character(len=16) :: smther2 ! post-smoother type: BJAC, AS integer(psb_ipk_) :: jsweeps2 ! post-smoother sweeps integer(psb_ipk_) :: degree2 ! degree for polynomial smoother + character(len=16) :: pvariant2 ! Polynomial smoother variant integer(psb_ipk_) :: novr2 ! number of overlap layers character(len=16) :: restr2 ! restriction over application of AS character(len=16) :: prol2 ! prolongation over application of AS @@ -169,7 +171,7 @@ program amg_d_pde3d character(len=16) :: variant2 ! AINV variant: LLK, etc integer(psb_ipk_) :: fill2 ! fill-in for incomplete LU factorization integer(psb_ipk_) :: invfill2 ! Inverse fill-in for INVK - real(psb_dpk_) :: thr2 ! threshold for ILUT factorization + real(psb_dpk_) :: thr2 ! threshold for ILUT factorization ! coarsest-level solver character(len=16) :: cmat ! coarsest matrix layout: REPL, DIST @@ -291,7 +293,8 @@ program amg_d_pde3d call prec%set('smoother_sweeps', p_choice%jsweeps, info) call prec%set('sub_solve', p_choice%solve, info) call prec%set('solver_sweeps', p_choice%ssweeps, info) - call prec%set('smoother_degree', p_choice%degree, info) + call prec%set('poly_degree', p_choice%degree, info) + call prec%set('poly_variant', p_choice%variant, info) if (psb_toupper(p_choice%solve)=='MUMPS') & & call prec%set('mumps_loc_glob','local_solver',info) call prec%set('sub_fillin', p_choice%fill, info) @@ -339,7 +342,8 @@ program amg_d_pde3d call prec%set('smoother_type', p_choice%smther, info) call prec%set('smoother_sweeps', p_choice%jsweeps, info) - call prec%set('smoother_degree', p_choice%degree, info) + call prec%set('poly_degree', p_choice%degree, info) + call prec%set('poly_variant', p_choice%variant, info) select case (psb_toupper(p_choice%smther)) case ('GS','BWGS','FBGS','JACOBI','L1-JACOBI','L1-FBGS') @@ -370,7 +374,8 @@ program amg_d_pde3d if (psb_toupper(p_choice%smther2) /= 'NONE') then call prec%set('smoother_type', p_choice%smther2, info,pos='post') call prec%set('smoother_sweeps', p_choice%jsweeps2, info,pos='post') - call prec%set('smoother_degree', p_choice%degree2, info,pos='post') + call prec%set('poly_degree', p_choice%degree2, info,pos='post') + call prec%set('poly_variant', p_choice%variant2, info,pos='post') select case (psb_toupper(p_choice%smther2)) case ('GS','BWGS','FBGS','JACOBI','L1-JACOBI','L1-FBGS') ! do nothing @@ -407,7 +412,6 @@ program amg_d_pde3d call prec%set('coarse_sweeps', p_choice%cjswp, info) end select -!!$ call prec%descr(info,iout=psb_out_unit) ! build the preconditioner call psb_barrier(ctxt) @@ -587,7 +591,8 @@ contains ! First smoother / 1-lev preconditioner call read_data(prec%smther,inp_unit) ! smoother type call read_data(prec%jsweeps,inp_unit) ! (pre-)smoother / 1-lev prec sweeps - call read_data(prec%degree,inp_unit) ! (pre-)smoother / 1-lev prec sweeps + call read_data(prec%degree,inp_unit) ! Degree of Polynomial smoother + call read_data(prec%variant,inp_unit) ! variant for Polynomial call read_data(prec%novr,inp_unit) ! number of overlap layers call read_data(prec%restr,inp_unit) ! restriction over application of AS call read_data(prec%prol,inp_unit) ! prolongation over application of AS @@ -600,12 +605,13 @@ contains ! 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 - call read_data(prec%degree2,inp_unit) ! (post-)smoother sweeps + call read_data(prec%degree2,inp_unit) ! Degree of Polynomial smoother + call read_data(prec%variant2,inp_unit) ! Polynomial smoother variant call read_data(prec%novr2,inp_unit) ! number of overlap layers call read_data(prec%restr2,inp_unit) ! restriction over application of AS call read_data(prec%prol2,inp_unit) ! prolongation over application of AS call read_data(prec%solve2,inp_unit) ! local subsolver - call read_data(prec%ssweeps2,inp_unit) ! inner solver sweeps + call read_data(prec%ssweeps2,inp_unit) ! inner solver sweeps call read_data(prec%variant2,inp_unit) ! AINV variant call read_data(prec%fill2,inp_unit) ! fill-in for incomplete LU call read_data(prec%invfill2,inp_unit) !Inverse fill-in for INVK @@ -672,6 +678,7 @@ contains call psb_bcast(ctxt,prec%smther) call psb_bcast(ctxt,prec%jsweeps) call psb_bcast(ctxt,prec%degree) + call psb_bcast(ctxt,prec%variant) call psb_bcast(ctxt,prec%novr) call psb_bcast(ctxt,prec%restr) call psb_bcast(ctxt,prec%prol) @@ -685,6 +692,7 @@ contains call psb_bcast(ctxt,prec%smther2) call psb_bcast(ctxt,prec%jsweeps2) call psb_bcast(ctxt,prec%degree2) + call psb_bcast(ctxt,prec%variant2) call psb_bcast(ctxt,prec%novr2) call psb_bcast(ctxt,prec%restr2) call psb_bcast(ctxt,prec%prol2) diff --git a/samples/advanced/pdegen/runs/amg_pde3d.inp b/samples/advanced/pdegen/runs/amg_pde3d.inp index ac39c4af..96b1b592 100644 --- a/samples/advanced/pdegen/runs/amg_pde3d.inp +++ b/samples/advanced/pdegen/runs/amg_pde3d.inp @@ -4,16 +4,17 @@ CSR ! Storage format CSR COO JAD CONST ! PDECOEFF: CONST, EXP, GAUSS Coefficients of the PDE BICGSTAB ! Iterative method: BiCGSTAB BiCGSTABL BiCG CG CGS FCG GCR RGMRES 2 ! ISTOPC -00500 ! ITMAX +00050 ! ITMAX 1 ! ITRACE 30 ! IRST (restart for RGMRES and BiCGSTABL) 1.d-6 ! EPS %%%%%%%%%%% Main preconditioner choices %%%%%%%%%%%%%%%% ML-VBM-VCYCLE-FBGS-D-BJAC ! Longer descriptive name for preconditioner (up to 20 chars) -BJAC ! Preconditioner type: NONE JACOBI GS FBGS BJAC AS ML +ML ! Preconditioner type: NONE JACOBI GS FBGS BJAC AS ML %%%%%%%%%%% First smoother (for all levels but coarsest) %%%%%%%%%%%%%%%% -BJAC ! Smoother type JACOBI FBGS GS BWGS BJAC AS. For 1-level, repeats previous. -1 ! Number of sweeps for smoother +FBGS ! Smoother type JACOBI FBGS GS BWGS BJAC AS. For 1-level, repeats previous. +4 ! Number of sweeps for smoother +4 ! degree for polynomial smoother 0 ! Number of overlap layers for AS preconditioner HALO ! AS restriction operator: NONE HALO NONE ! AS prolongation operator: NONE SUM AVG @@ -24,8 +25,9 @@ LLK ! AINV variant 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 +NONE ! Second (post) smoother, ignored if NONE +4 ! Number of sweeps for (post) smoother +4 ! degree for polynomial smoother 0 ! Number of overlap layers for AS preconditioner HALO ! AS restriction operator: NONE HALO NONE ! AS prolongation operator: NONE SUM AVG