diff --git a/amgprec/Makefile b/amgprec/Makefile index 79842c27..c6ccb4b7 100644 --- a/amgprec/Makefile +++ b/amgprec/Makefile @@ -21,7 +21,7 @@ DMODOBJS=amg_d_prec_type.o \ SMODOBJS=amg_s_prec_type.o amg_s_ilu_fact_mod.o \ amg_s_inner_mod.o amg_s_ilu_solver.o amg_s_diag_solver.o amg_s_jac_smoother.o amg_s_as_smoother.o \ - amg_s_slu_solver.o amg_s_id_solver.o\ + amg_s_poly_smoother.o amg_s_slu_solver.o amg_s_id_solver.o\ amg_s_base_solver_mod.o amg_s_base_smoother_mod.o amg_s_onelev_mod.o \ amg_s_gs_solver.o amg_s_mumps_solver.o amg_s_jac_solver.o \ amg_s_base_aggregator_mod.o \ @@ -166,6 +166,7 @@ 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_poly_coeff_mod.o +amg_s_poly_smoother.o: amg_s_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 5d75e274..3434d675 100644 --- a/amgprec/amg_base_prec_type.F90 +++ b/amgprec/amg_base_prec_type.F90 @@ -687,10 +687,10 @@ contains & ml_names(pm%ml_cycle) select case (pm%ml_cycle) case (amg_add_ml_) - write(iout,*) ' Number of smoother sweeps : ',& + write(iout,*) ' Number of smoother sweeps/degree : ',& & pm%sweeps_pre case (amg_mult_ml_,amg_vcycle_ml_, amg_wcycle_ml_, amg_kcycle_ml_, amg_kcyclesym_ml_) - write(iout,*) ' Number of smoother sweeps : pre: ',& + write(iout,*) ' Number of smoother sweeps/degree : pre: ',& & pm%sweeps_pre ,' post: ', pm%sweeps_post end select diff --git a/amgprec/impl/amg_cfile_prec_descr.f90 b/amgprec/impl/amg_cfile_prec_descr.f90 index 396a9467..bdaa7c41 100644 --- a/amgprec/impl/amg_cfile_prec_descr.f90 +++ b/amgprec/impl/amg_cfile_prec_descr.f90 @@ -170,7 +170,7 @@ subroutine amg_cfile_prec_descr(prec,info,iout,root, verbosity,prefix) call prec%precv(1)%sm%descr(info,iout=iout_,prefix=prefix) nswps = prec%precv(1)%parms%sweeps_pre end if - write(iout_,*) trim(prefix_), ' Number of sweeps : ',nswps + write(iout_,*) trim(prefix_), ' Number of sweeps/degree : ',nswps write(iout_,*) trim(prefix_) else if (nlev > 1) then diff --git a/amgprec/impl/amg_cprecinit.F90 b/amgprec/impl/amg_cprecinit.F90 index 2335281b..5c08b1d1 100644 --- a/amgprec/impl/amg_cprecinit.F90 +++ b/amgprec/impl/amg_cprecinit.F90 @@ -98,6 +98,7 @@ subroutine amg_cprecinit(ctxt,prec,ptype,info) use amg_c_diag_solver use amg_c_ilu_solver use amg_c_gs_solver + #if defined(HAVE_SLU_) use amg_c_slu_solver #endif @@ -152,7 +153,6 @@ subroutine amg_cprecinit(ctxt,prec,ptype,info) if (info /= psb_success_) return allocate(amg_c_diag_solver_type :: prec%precv(ilev_)%sm%sv, stat=info) call prec%precv(ilev_)%default() - case ('L1-DIAG','L1-JACOBI','L1_DIAG','L1_JACOBI') nlev_ = 1 ilev_ = 1 diff --git a/amgprec/impl/amg_dfile_prec_descr.f90 b/amgprec/impl/amg_dfile_prec_descr.f90 index 3213df29..4cdde58f 100644 --- a/amgprec/impl/amg_dfile_prec_descr.f90 +++ b/amgprec/impl/amg_dfile_prec_descr.f90 @@ -170,7 +170,7 @@ subroutine amg_dfile_prec_descr(prec,info,iout,root, verbosity,prefix) call prec%precv(1)%sm%descr(info,iout=iout_,prefix=prefix) nswps = prec%precv(1)%parms%sweeps_pre end if - write(iout_,*) trim(prefix_), ' Number of sweeps : ',nswps + write(iout_,*) trim(prefix_), ' Number of sweeps/degree : ',nswps write(iout_,*) trim(prefix_) else if (nlev > 1) then diff --git a/amgprec/impl/amg_dprecinit.F90 b/amgprec/impl/amg_dprecinit.F90 index 8f3c0cb6..176ceb0d 100644 --- a/amgprec/impl/amg_dprecinit.F90 +++ b/amgprec/impl/amg_dprecinit.F90 @@ -93,12 +93,13 @@ subroutine amg_dprecinit(ctxt,prec,ptype,info) use psb_base_mod use amg_d_prec_mod, amg_protect_name => amg_dprecinit use amg_d_jac_smoother - use amg_d_poly_smoother use amg_d_as_smoother use amg_d_id_solver use amg_d_diag_solver use amg_d_ilu_solver use amg_d_gs_solver + use amg_d_poly_smoother + #if defined(HAVE_UMF_) use amg_d_umf_solver #endif @@ -156,7 +157,6 @@ subroutine amg_dprecinit(ctxt,prec,ptype,info) if (info /= psb_success_) return allocate(amg_d_diag_solver_type :: prec%precv(ilev_)%sm%sv, stat=info) call prec%precv(ilev_)%default() - case ('POLY') nlev_ = 1 ilev_ = 1 @@ -165,7 +165,6 @@ subroutine amg_dprecinit(ctxt,prec,ptype,info) if (info /= psb_success_) return 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') nlev_ = 1 ilev_ = 1 diff --git a/amgprec/impl/amg_sfile_prec_descr.f90 b/amgprec/impl/amg_sfile_prec_descr.f90 index 5996e2a1..07bde2a4 100644 --- a/amgprec/impl/amg_sfile_prec_descr.f90 +++ b/amgprec/impl/amg_sfile_prec_descr.f90 @@ -170,7 +170,7 @@ subroutine amg_sfile_prec_descr(prec,info,iout,root, verbosity,prefix) call prec%precv(1)%sm%descr(info,iout=iout_,prefix=prefix) nswps = prec%precv(1)%parms%sweeps_pre end if - write(iout_,*) trim(prefix_), ' Number of sweeps : ',nswps + write(iout_,*) trim(prefix_), ' Number of sweeps/degree : ',nswps write(iout_,*) trim(prefix_) else if (nlev > 1) then diff --git a/amgprec/impl/amg_sprecinit.F90 b/amgprec/impl/amg_sprecinit.F90 index cd91708d..3ef58406 100644 --- a/amgprec/impl/amg_sprecinit.F90 +++ b/amgprec/impl/amg_sprecinit.F90 @@ -98,6 +98,8 @@ subroutine amg_sprecinit(ctxt,prec,ptype,info) use amg_s_diag_solver use amg_s_ilu_solver use amg_s_gs_solver + use amg_s_poly_smoother + #if defined(HAVE_SLU_) use amg_s_slu_solver #endif @@ -152,7 +154,14 @@ subroutine amg_sprecinit(ctxt,prec,ptype,info) if (info /= psb_success_) return allocate(amg_s_diag_solver_type :: prec%precv(ilev_)%sm%sv, stat=info) call prec%precv(ilev_)%default() - + case ('POLY') + nlev_ = 1 + ilev_ = 1 + allocate(prec%precv(nlev_),stat=info) + allocate(amg_s_poly_smoother_type :: prec%precv(ilev_)%sm, stat=info) + if (info /= psb_success_) return + allocate(amg_s_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') nlev_ = 1 ilev_ = 1 diff --git a/amgprec/impl/amg_zfile_prec_descr.f90 b/amgprec/impl/amg_zfile_prec_descr.f90 index f3002cfd..54e687bf 100644 --- a/amgprec/impl/amg_zfile_prec_descr.f90 +++ b/amgprec/impl/amg_zfile_prec_descr.f90 @@ -170,7 +170,7 @@ subroutine amg_zfile_prec_descr(prec,info,iout,root, verbosity,prefix) call prec%precv(1)%sm%descr(info,iout=iout_,prefix=prefix) nswps = prec%precv(1)%parms%sweeps_pre end if - write(iout_,*) trim(prefix_), ' Number of sweeps : ',nswps + write(iout_,*) trim(prefix_), ' Number of sweeps/degree : ',nswps write(iout_,*) trim(prefix_) else if (nlev > 1) then diff --git a/amgprec/impl/amg_zprecinit.F90 b/amgprec/impl/amg_zprecinit.F90 index 9bef4d1a..ab88c80c 100644 --- a/amgprec/impl/amg_zprecinit.F90 +++ b/amgprec/impl/amg_zprecinit.F90 @@ -98,6 +98,7 @@ subroutine amg_zprecinit(ctxt,prec,ptype,info) use amg_z_diag_solver use amg_z_ilu_solver use amg_z_gs_solver + #if defined(HAVE_UMF_) use amg_z_umf_solver #endif @@ -155,7 +156,6 @@ subroutine amg_zprecinit(ctxt,prec,ptype,info) if (info /= psb_success_) return allocate(amg_z_diag_solver_type :: prec%precv(ilev_)%sm%sv, stat=info) call prec%precv(ilev_)%default() - case ('L1-DIAG','L1-JACOBI','L1_DIAG','L1_JACOBI') nlev_ = 1 ilev_ = 1 diff --git a/amgprec/impl/level/amg_c_base_onelev_csetc.F90 b/amgprec/impl/level/amg_c_base_onelev_csetc.F90 index a372d972..12ed5ea6 100644 --- a/amgprec/impl/level/amg_c_base_onelev_csetc.F90 +++ b/amgprec/impl/level/amg_c_base_onelev_csetc.F90 @@ -188,16 +188,11 @@ subroutine amg_c_base_onelev_csetc(lv,what,val,info,pos,idx) case ('NONE','NOPREC','FACT_NONE') call lv%set(amg_c_id_solver_mold,info,pos=pos) - case ('DIAG') + case ('DIAG','JACOBI') call lv%set(amg_c_diag_solver_mold,info,pos=pos) - case ('JACOBI') - call lv%set(amg_c_jac_solver_mold,info,pos=pos) - - case ('L1-DIAG') + case ('L1-DIAG','L1-JACOBI') call lv%set(amg_c_l1_diag_solver_mold,info,pos=pos) - case ('L1-JACOBI') - call lv%set(amg_c_l1_jac_solver_mold,info,pos=pos) case ('GS','FGS','FWGS') call lv%set(amg_c_gs_solver_mold,info,pos=pos) diff --git a/amgprec/impl/level/amg_d_base_onelev_csetc.F90 b/amgprec/impl/level/amg_d_base_onelev_csetc.F90 index d4d37262..3b2e4f37 100644 --- a/amgprec/impl/level/amg_d_base_onelev_csetc.F90 +++ b/amgprec/impl/level/amg_d_base_onelev_csetc.F90 @@ -43,9 +43,9 @@ subroutine amg_d_base_onelev_csetc(lv,what,val,info,pos,idx) use amg_d_dec_aggregator_mod use amg_d_symdec_aggregator_mod use amg_d_parmatch_aggregator_mod + use amg_d_poly_smoother 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_jac_solver @@ -85,7 +85,6 @@ subroutine amg_d_base_onelev_csetc(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_jac_solver_type) :: amg_d_jac_solver_mold @@ -97,6 +96,7 @@ subroutine amg_d_base_onelev_csetc(lv,what,val,info,pos,idx) type(amg_d_ainv_solver_type) :: amg_d_ainv_solver_mold 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 #if defined(HAVE_UMF_) type(amg_d_umf_solver_type) :: amg_d_umf_solver_mold #endif @@ -161,7 +161,6 @@ 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 ('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_s_base_onelev_csetc.F90 b/amgprec/impl/level/amg_s_base_onelev_csetc.F90 index 60329291..f76d4b61 100644 --- a/amgprec/impl/level/amg_s_base_onelev_csetc.F90 +++ b/amgprec/impl/level/amg_s_base_onelev_csetc.F90 @@ -43,6 +43,7 @@ subroutine amg_s_base_onelev_csetc(lv,what,val,info,pos,idx) use amg_s_dec_aggregator_mod use amg_s_symdec_aggregator_mod use amg_s_parmatch_aggregator_mod + use amg_s_poly_smoother use amg_s_jac_smoother use amg_s_as_smoother use amg_s_diag_solver @@ -89,6 +90,7 @@ subroutine amg_s_base_onelev_csetc(lv,what,val,info,pos,idx) type(amg_s_ainv_solver_type) :: amg_s_ainv_solver_mold type(amg_s_invk_solver_type) :: amg_s_invk_solver_mold type(amg_s_invt_solver_type) :: amg_s_invt_solver_mold + type(amg_s_poly_smoother_type) :: amg_s_poly_smoother_mold #if defined(HAVE_SLU_) type(amg_s_slu_solver_type) :: amg_s_slu_solver_mold #endif @@ -144,6 +146,9 @@ subroutine amg_s_base_onelev_csetc(lv,what,val,info,pos,idx) call lv%set(amg_s_as_smoother_mold,info,pos=pos) if (info == 0) call lv%set(amg_s_ilu_solver_mold,info,pos=pos) + case ('POLY') + call lv%set(amg_s_poly_smoother_mold,info,pos=pos) + if (info == 0) call lv%set(amg_s_l1_diag_solver_mold,info,pos=pos) case ('GS','FWGS') call lv%set(amg_s_jac_smoother_mold,info,pos='pre') if (info == 0) call lv%set(amg_s_gs_solver_mold,info,pos='pre') @@ -189,16 +194,11 @@ subroutine amg_s_base_onelev_csetc(lv,what,val,info,pos,idx) case ('NONE','NOPREC','FACT_NONE') call lv%set(amg_s_id_solver_mold,info,pos=pos) - case ('DIAG') + case ('DIAG','JACOBI') call lv%set(amg_s_diag_solver_mold,info,pos=pos) - case ('JACOBI') - call lv%set(amg_s_jac_solver_mold,info,pos=pos) - - case ('L1-DIAG') + case ('L1-DIAG','L1-JACOBI') call lv%set(amg_s_l1_diag_solver_mold,info,pos=pos) - case ('L1-JACOBI') - call lv%set(amg_s_l1_jac_solver_mold,info,pos=pos) case ('GS','FGS','FWGS') call lv%set(amg_s_gs_solver_mold,info,pos=pos) diff --git a/amgprec/impl/level/amg_z_base_onelev_csetc.F90 b/amgprec/impl/level/amg_z_base_onelev_csetc.F90 index 5ca4233c..e0eddc4d 100644 --- a/amgprec/impl/level/amg_z_base_onelev_csetc.F90 +++ b/amgprec/impl/level/amg_z_base_onelev_csetc.F90 @@ -200,16 +200,11 @@ subroutine amg_z_base_onelev_csetc(lv,what,val,info,pos,idx) case ('NONE','NOPREC','FACT_NONE') call lv%set(amg_z_id_solver_mold,info,pos=pos) - case ('DIAG') + case ('DIAG','JACOBI') call lv%set(amg_z_diag_solver_mold,info,pos=pos) - case ('JACOBI') - call lv%set(amg_z_jac_solver_mold,info,pos=pos) - - case ('L1-DIAG') + case ('L1-DIAG','L1-JACOBI') call lv%set(amg_z_l1_diag_solver_mold,info,pos=pos) - case ('L1-JACOBI') - call lv%set(amg_z_l1_jac_solver_mold,info,pos=pos) case ('GS','FGS','FWGS') call lv%set(amg_z_gs_solver_mold,info,pos=pos) diff --git a/amgprec/impl/smoother/Makefile b/amgprec/impl/smoother/Makefile index 58884bc8..89a1906e 100644 --- a/amgprec/impl/smoother/Makefile +++ b/amgprec/impl/smoother/Makefile @@ -153,6 +153,17 @@ amg_s_jac_smoother_csetr.o \ amg_s_l1_jac_smoother_bld.o \ amg_s_l1_jac_smoother_descr.o \ amg_s_l1_jac_smoother_clone.o \ +amg_s_poly_smoother_apply_vect.o \ +amg_s_poly_smoother_bld.o \ +amg_s_poly_smoother_cnv.o \ +amg_s_poly_smoother_clone.o \ +amg_s_poly_smoother_clone_settings.o \ +amg_s_poly_smoother_clear_data.o \ +amg_s_poly_smoother_descr.o \ +amg_s_poly_smoother_dmp.o \ +amg_s_poly_smoother_csetc.o \ +amg_s_poly_smoother_cseti.o \ +amg_s_poly_smoother_csetr.o \ amg_z_as_smoother_apply.o \ amg_z_as_smoother_apply_vect.o \ amg_z_as_smoother_bld.o \ diff --git a/amgprec/impl/smoother/amg_c_jac_smoother_apply_vect.f90 b/amgprec/impl/smoother/amg_c_jac_smoother_apply_vect.f90 index 040f9fb3..a4238980 100644 --- a/amgprec/impl/smoother/amg_c_jac_smoother_apply_vect.f90 +++ b/amgprec/impl/smoother/amg_c_jac_smoother_apply_vect.f90 @@ -175,7 +175,7 @@ subroutine amg_c_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,& if (info /= psb_success_) exit if ( sm%checkres.and.(mod(i,sm%checkiter) == 0) ) then - call psb_geaxpby(cone,x,czero,r,r,desc_data,info) + call psb_geaxpby(cone,x,czero,r,desc_data,info) call psb_spmm(-cone,sm%pa,ty,cone,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_apply_vect.f90 b/amgprec/impl/smoother/amg_d_jac_smoother_apply_vect.f90 index 1c206c27..cafb6c8c 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(.true..or.sm%checkres) then + if(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,10 +159,7 @@ 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)), @@ -176,10 +173,6 @@ 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,desc_data,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 9b085f62..a7a4202f 100644 --- a/amgprec/impl/smoother/amg_d_poly_smoother_apply_vect.f90 +++ b/amgprec/impl/smoother/amg_d_poly_smoother_apply_vect.f90 @@ -120,7 +120,7 @@ subroutine amg_d_poly_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,& & a_err='invalid wv size in smoother_apply') goto 9999 end if - + sm%pdegree = sweeps associate(tx => wv(1), ty => wv(2), tz => wv(3), r => wv(4)) call psb_geaxpby(done,x,dzero,r,desc_data,info) @@ -135,7 +135,7 @@ subroutine amg_d_poly_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,& ! b == x ! x == tx ! - do i=1, sm%pdegree + do i=1, sweeps ! 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) @@ -163,7 +163,15 @@ subroutine amg_d_poly_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,& ! b == x ! x == tx ! - do i=1, sm%pdegree + if (allocated(sm%poly_beta)) then + if (size(sm%poly_beta) /= sweeps) deallocate(sm%poly_beta) + end if + if (.not.allocated(sm%poly_beta)) then + call psb_realloc(sweeps,sm%poly_beta,info) + sm%poly_beta(1:sweeps) = amg_d_poly_beta_mat(1:sweeps,sweeps) + end if + + do i=1, sweeps ! 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) @@ -190,6 +198,8 @@ subroutine amg_d_poly_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,& ! b == x ! x == tx ! + sm%cf_a = amg_d_poly_a_vect(sweeps) + theta = (done+sm%cf_a)/2 delta = (done-sm%cf_a)/2 sigma = theta/delta @@ -198,7 +208,7 @@ subroutine amg_d_poly_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,& call psb_geaxpby((done/sm%rho_ba),ty,dzero,r,desc_data,info) call psb_geaxpby((done/theta),r,dzero,tz,desc_data,info) ! tz == d - do i=1, sm%pdegree + do i=1, sweeps ! x_{k+1} = x_k + d_k call psb_geaxpby(done,tz,done,tx,desc_data,info) ! diff --git a/amgprec/impl/smoother/amg_d_poly_smoother_bld.f90 b/amgprec/impl/smoother/amg_d_poly_smoother_bld.f90 index a7f0a72c..d9d39c03 100644 --- a/amgprec/impl/smoother/amg_d_poly_smoother_bld.f90 +++ b/amgprec/impl/smoother/amg_d_poly_smoother_bld.f90 @@ -74,39 +74,39 @@ subroutine amg_d_poly_smoother_bld(a,desc_a,sm,info,amold,vmold,imold) n_col = desc_a%get_local_cols() nrow_a = a%get_nrows() nztota = a%get_nzeros() + if (.false.) then + select case(sm%variant) + case(amg_poly_lottes_) + ! do nothing + case(amg_poly_lottes_beta_) + if ((1<=sm%pdegree).and.(sm%pdegree<=30)) then + call psb_realloc(sm%pdegree,sm%poly_beta,info) + sm%poly_beta(1:sm%pdegree) = amg_d_poly_beta_mat(1:sm%pdegree,sm%pdegree) + else + info = psb_err_internal_error_ + call psb_errpush(info,name,& + & a_err='invalid sm%degree for poly_beta') + goto 9999 + end if + case(amg_poly_new_) - select case(sm%variant) - case(amg_poly_lottes_) - ! do nothing - case(amg_poly_lottes_beta_) - if ((1<=sm%pdegree).and.(sm%pdegree<=30)) then - call psb_realloc(sm%pdegree,sm%poly_beta,info) - sm%poly_beta(1:sm%pdegree) = amg_d_poly_beta_mat(1:sm%pdegree,sm%pdegree) - else - info = psb_err_internal_error_ - call psb_errpush(info,name,& - & a_err='invalid sm%degree for poly_beta') - goto 9999 - end if - case(amg_poly_new_) + if ((1<=sm%pdegree).and.(sm%pdegree<=30)) then + !Ok + sm%cf_a = amg_d_poly_a_vect(sm%pdegree) + else + info = psb_err_internal_error_ + call psb_errpush(info,name,& + & a_err='invalid sm%degree for poly_a') + goto 9999 + end if - if ((1<=sm%pdegree).and.(sm%pdegree<=30)) then - !Ok - sm%cf_a = amg_d_poly_a_vect(sm%pdegree) - else + case default info = psb_err_internal_error_ call psb_errpush(info,name,& - & a_err='invalid sm%degree for poly_a') + & a_err='invalid sm%variant') goto 9999 - end if - - case default - info = psb_err_internal_error_ - call psb_errpush(info,name,& - & a_err='invalid sm%variant') - goto 9999 - end select - + end select + end if sm%pa => a if (.not.allocated(sm%sv)) then info = psb_err_internal_error_ diff --git a/amgprec/impl/smoother/amg_d_poly_smoother_clear_data.f90 b/amgprec/impl/smoother/amg_d_poly_smoother_clear_data.f90 index ac526bca..a6df5486 100644 --- a/amgprec/impl/smoother/amg_d_poly_smoother_clear_data.f90 +++ b/amgprec/impl/smoother/amg_d_poly_smoother_clear_data.f90 @@ -44,7 +44,7 @@ subroutine amg_d_poly_smoother_clear_data(sm,info) class(amg_d_poly_smoother_type), intent(inout) :: sm integer(psb_ipk_), intent(out) :: info integer(psb_ipk_) :: err_act - character(len=20) :: name='d_poly_smoother_clear_data' + character(len=20) :: name='amg_d_poly_smoother_clear_data' call psb_erractionsave(err_act) diff --git a/amgprec/impl/smoother/amg_d_poly_smoother_clone_settings.f90 b/amgprec/impl/smoother/amg_d_poly_smoother_clone_settings.f90 index 1fbdac37..d72cce67 100644 --- a/amgprec/impl/smoother/amg_d_poly_smoother_clone_settings.f90 +++ b/amgprec/impl/smoother/amg_d_poly_smoother_clone_settings.f90 @@ -56,7 +56,13 @@ subroutine amg_d_poly_smoother_clone_settings(sm,smout,info) smout%pa => null() smout%pdegree = sm%pdegree - + smout%variant = sm%variant + smout%cf_a = sm%cf_a + smout%rho_ba = sm%rho_ba + smout%rho_estimate = sm%rho_estimate + smout%rho_estimate_iterations = sm%rho_estimate_iterations + smout%poly_beta = sm%poly_beta + if (allocated(smout%sv)) then if (.not.same_type_as(sm%sv,smout%sv)) then call smout%sv%free(info) diff --git a/amgprec/impl/smoother/amg_d_poly_smoother_csetc.f90 b/amgprec/impl/smoother/amg_d_poly_smoother_csetc.f90 index e61b09e3..3d0ac0fe 100644 --- a/amgprec/impl/smoother/amg_d_poly_smoother_csetc.f90 +++ b/amgprec/impl/smoother/amg_d_poly_smoother_csetc.f90 @@ -38,7 +38,7 @@ subroutine amg_d_poly_smoother_csetc(sm,what,val,info,idx) use psb_base_mod - use amg_d_poly_smoother, amg_protect_nam => amg_d_poly_smoother_csetc + use amg_d_poly_smoother, amg_protect_name => amg_d_poly_smoother_csetc Implicit None ! Arguments class(amg_d_poly_smoother_type), intent(inout) :: sm diff --git a/amgprec/impl/smoother/amg_d_poly_smoother_descr.f90 b/amgprec/impl/smoother/amg_d_poly_smoother_descr.f90 index 97521259..1535388c 100644 --- a/amgprec/impl/smoother/amg_d_poly_smoother_descr.f90 +++ b/amgprec/impl/smoother/amg_d_poly_smoother_descr.f90 @@ -40,8 +40,6 @@ subroutine amg_d_poly_smoother_descr(sm,info,iout,coarse,prefix) use psb_base_mod use amg_d_diag_solver use amg_d_poly_smoother, amg_protect_name => amg_d_poly_smoother_descr - use amg_d_diag_solver - use amg_d_gs_solver Implicit None @@ -81,18 +79,18 @@ subroutine amg_d_poly_smoother_descr(sm,info,iout,coarse,prefix) select case(sm%variant) case(amg_poly_lottes_) write(iout_,*) trim(prefix_), ' variant: ','POLY_LOTTES' - write(iout_,*) trim(prefix_), ' Degree: ',sm%pdegree + !write(iout_,*) trim(prefix_), ' Degree: ',sm%pdegree write(iout_,*) trim(prefix_), ' rho_ba: ',sm%rho_ba case(amg_poly_lottes_beta_) write(iout_,*) trim(prefix_), ' variant: ','POLY_LOTTES_BETA' - write(iout_,*) trim(prefix_), ' Degree: ',sm%pdegree + !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) + !if (allocated(sm%poly_beta)) write(iout_,*) trim(prefix_), ' Coefficients: ',sm%poly_beta(1:sm%pdegree) case(amg_poly_new_) write(iout_,*) trim(prefix_), ' variant: ','POLY_NEW' - write(iout_,*) trim(prefix_), ' Degree: ',sm%pdegree + !write(iout_,*) trim(prefix_), ' Degree: ',sm%pdegree write(iout_,*) trim(prefix_), ' rho_ba: ',sm%rho_ba - write(iout_,*) trim(prefix_), ' Coefficient: ',sm%cf_a + !write(iout_,*) trim(prefix_), ' Coefficient: ',sm%cf_a case default write(iout_,*) trim(prefix_), ' variant: ','UNKNOWN???' end select diff --git a/amgprec/impl/smoother/amg_s_jac_smoother_apply_vect.f90 b/amgprec/impl/smoother/amg_s_jac_smoother_apply_vect.f90 index 9fe3888a..fff7ac1e 100644 --- a/amgprec/impl/smoother/amg_s_jac_smoother_apply_vect.f90 +++ b/amgprec/impl/smoother/amg_s_jac_smoother_apply_vect.f90 @@ -175,7 +175,7 @@ subroutine amg_s_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,& if (info /= psb_success_) exit if ( sm%checkres.and.(mod(i,sm%checkiter) == 0) ) then - call psb_geaxpby(sone,x,szero,r,r,desc_data,info) + call psb_geaxpby(sone,x,szero,r,desc_data,info) call psb_spmm(-sone,sm%pa,ty,sone,r,desc_data,info) res = psb_genrm2(r,desc_data,info) if( sm%printres ) then diff --git a/amgprec/impl/smoother/amg_s_poly_smoother_apply_vect.f90 b/amgprec/impl/smoother/amg_s_poly_smoother_apply_vect.f90 new file mode 100644 index 00000000..76be3e99 --- /dev/null +++ b/amgprec/impl/smoother/amg_s_poly_smoother_apply_vect.f90 @@ -0,0 +1,453 @@ +! +! +! AMG4PSBLAS version 1.0 +! Algebraic Multigrid Package +! based on PSBLAS (Parallel Sparse BLAS version 3.7) +! +! (C) Copyright 2021 +! +! Salvatore Filippone +! Pasqua D'Ambra +! Daniela di Serafino +! +! Redistribution and use in source and binary forms, with or without +! modification, are permitted provided that the following conditions +! are met: +! 1. Redistributions of source code must retain the above copyright +! notice, this list of conditions and the following disclaimer. +! 2. Redistributions in binary form must reproduce the above copyright +! notice, this list of conditions, and the following disclaimer in the +! documentation and/or other materials provided with the distribution. +! 3. The name of the AMG4PSBLAS group or the names of its contributors may +! not be used to endorse or promote products derived from this +! software without specific written permission. +! +! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +! ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED +! TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR +! PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE AMG4PSBLAS GROUP OR ITS CONTRIBUTORS +! BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +! CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +! SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +! INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +! CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +! ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +! POSSIBILITY OF SUCH DAMAGE. +! +! +subroutine amg_s_poly_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,& + & sweeps,work,wv,info,init,initu) + + use psb_base_mod + use amg_s_diag_solver + use psb_base_krylov_conv_mod, only : log_conv + use amg_s_poly_smoother, amg_protect_name => amg_s_poly_smoother_apply_vect + implicit none + type(psb_desc_type), intent(in) :: desc_data + class(amg_s_poly_smoother_type), intent(inout) :: sm + type(psb_s_vect_type),intent(inout) :: x + type(psb_s_vect_type),intent(inout) :: y + real(psb_spk_),intent(in) :: alpha,beta + character(len=1),intent(in) :: trans + integer(psb_ipk_), intent(in) :: sweeps + real(psb_spk_),target, intent(inout) :: work(:) + type(psb_s_vect_type),intent(inout) :: wv(:) + integer(psb_ipk_), intent(out) :: info + character, intent(in), optional :: init + type(psb_s_vect_type),intent(inout), optional :: initu + ! + integer(psb_ipk_) :: n_row,n_col + type(psb_s_vect_type) :: tx, ty, tz, r + real(psb_spk_), pointer :: aux(:) + type(psb_ctxt_type) :: ctxt + integer(psb_ipk_) :: np, me, i, err_act + character :: trans_, init_ + real(psb_spk_) :: res, resdenum + character(len=20) :: name='d_poly_smoother_apply_v' + + call psb_erractionsave(err_act) + + info = psb_success_ + ctxt = desc_data%get_context() + call psb_info(ctxt,me,np) + + + if (present(init)) then + init_ = psb_toupper(init) + else + init_='Z' + end if + + trans_ = psb_toupper(trans) + select case(trans_) + case('N') + case('T','C') + case default + call psb_errpush(psb_err_iarg_invalid_i_,name) + goto 9999 + end select + + if (.not.allocated(sm%sv)) then + info = 1121 + call psb_errpush(info,name) + goto 9999 + end if + + n_row = desc_data%get_local_rows() + n_col = desc_data%get_local_cols() + + if (4*n_col <= size(work)) then + aux => work(:) + else + allocate(aux(4*n_col),stat=info) + if (info /= psb_success_) then + info=psb_err_alloc_request_ + call psb_errpush(info,name,& + & i_err=(/4*n_col,izero,izero,izero,izero/),& + & a_err='real(psb_spk_)') + goto 9999 + end if + endif +!!$ if (me == 0) write(0,*) name,' Unimplemented apply_vect ' +!!$ info =psb_err_internal_error_ +!!$ call psb_errpush(info,& +!!$ & name,a_err='Error in sub_aply Polynomial not implemented') +!!$ goto 9999 + + if (size(wv) < 4) then + info = psb_err_internal_error_ + call psb_errpush(info,name,& + & a_err='invalid wv size in smoother_apply') + goto 9999 + end if + sm%pdegree = sweeps + associate(tx => wv(1), ty => wv(2), tz => wv(3), r => wv(4)) + + call psb_geaxpby(sone,x,szero,r,desc_data,info) + call tx%zero() + call ty%zero() + call tz%zero() + + select case(sm%variant) + case(amg_poly_lottes_) + block + real(psb_spk_) :: cz, cr + ! b == x + ! x == tx + ! + do i=1, sweeps + ! B r_{k-1} + call sm%sv%apply(sone,r,szero,ty,desc_data,trans_,aux,wv(5:),info,init='Z') + cz = (2*i*sone-3)/(2*i*sone+sone) + cr = (8*i*sone-4)/((2*i*sone+sone)*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(sone,tz,sone,tx,desc_data,info) + if (.false.) then + call psb_geaxpby(sone,x,szero,r,desc_data,info) + call psb_spmm(-sone,sm%pa,tx,sone,r,desc_data,info,work=aux,trans=trans_) + else + call psb_spmm(-sone,sm%pa,tz,sone,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 + end block + + case(amg_poly_lottes_beta_) + + block + real(psb_spk_) :: cz, cr + ! b == x + ! x == tx + ! + if (allocated(sm%poly_beta)) then + if (size(sm%poly_beta) /= sweeps) deallocate(sm%poly_beta) + end if + if (.not.allocated(sm%poly_beta)) then + call psb_realloc(sweeps,sm%poly_beta,info) + sm%poly_beta(1:sweeps) = amg_d_poly_beta_mat(1:sweeps,sweeps) + end if + + do i=1, sweeps + ! B r_{k-1} + call sm%sv%apply(sone,r,szero,ty,desc_data,trans_,aux,wv(5:),info,init='Z') + cz = (2*i*sone-3)/(2*i*sone+sone) + cr = (8*i*sone-4)/((2*i*sone+sone)*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,sone,tx,desc_data,info) + if (.false.) then + call psb_geaxpby(sone,x,szero,r,desc_data,info) + call psb_spmm(-sone,sm%pa,tx,sone,r,desc_data,info,work=aux,trans=trans_) + else + call psb_spmm(-sone,sm%pa,tz,sone,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 + end block + + case(amg_poly_new_) + block + real(psb_spk_) :: sigma, theta, delta, rho_old, rho + ! b == x + ! x == tx + ! + sm%cf_a = amg_d_poly_a_vect(sweeps) + + theta = (sone+sm%cf_a)/2 + delta = (sone-sm%cf_a)/2 + sigma = theta/delta + rho_old = sone/sigma + call sm%sv%apply(sone,r,szero,ty,desc_data,trans_,aux,wv(5:),info,init='Z') + call psb_geaxpby((sone/sm%rho_ba),ty,szero,r,desc_data,info) + call psb_geaxpby((sone/theta),r,szero,tz,desc_data,info) + ! tz == d + do i=1, sweeps + ! x_{k+1} = x_k + d_k + call psb_geaxpby(sone,tz,sone,tx,desc_data,info) + ! + ! r_{k-1} = r_k - (1/rho(BA)) B A d_k + call psb_spmm(sone,sm%pa,tz,szero,ty,desc_data,info,work=aux,trans=trans_) + call sm%sv%apply(-sone,ty,sone,r,desc_data,trans_,aux,wv(5:),info,init='Z') + + ! + ! d_{k+1} = (rho rho_old) d_k + 2(rho/delta) r_{k+1} + rho = sone/(2*sigma - rho_old) + call psb_geaxpby((2*rho/delta),r,(rho*rho_old),tz,desc_data,info) +!!$ res = psb_genrm2(r,desc_data,info) +!!$ write(0,*) 'Polynomial smoother ',i,res + ! x_k = x_{k-1} + z_k + end do + end block + + + 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 + info=psb_err_internal_error_ + call psb_errpush(info,name,& + & a_err='polynomial smoother') + goto 9999 + end if + end associate + + + + +!!$ if ((.not.sm%sv%is_iterative()).and.((sweeps == 1).or.(sm%nd_nnz_tot==0))) then +!!$ ! if .not.sv%is_iterative, there's no need to pass init +!!$ call sm%sv%apply(alpha,x,beta,y,desc_data,trans_,aux,wv,info) +!!$ +!!$ if (info /= psb_success_) then +!!$ call psb_errpush(psb_err_internal_error_,& +!!$ & name,a_err='Error in sub_aply Jacobi Sweeps = 1') +!!$ goto 9999 +!!$ endif +!!$ +!!$ else if (sweeps >= 0) then +!!$ select type (smsv => sm%sv) +!!$ class is (amg_s_diag_solver_type) +!!$ ! +!!$ ! This means we are dealing with a pure Jacobi smoother/solver. +!!$ ! +!!$ associate(tx => wv(1), ty => wv(2)) +!!$ select case (init_) +!!$ case('Z') +!!$ +!!$ call sm%sv%apply(sone,x,szero,ty,desc_data,trans_,aux,wv(3:),info,init='Z') +!!$ +!!$ case('Y') +!!$ call psb_geaxpby(sone,x,szero,tx,desc_data,info) +!!$ call psb_geaxpby(sone,y,szero,ty,desc_data,info) +!!$ call psb_spmm(-sone,sm%pa,ty,sone,tx,desc_data,info,work=aux,trans=trans_) +!!$ call sm%sv%apply(sone,tx,szero,ty,desc_data,trans_,aux,wv(3:),info,init='Y') +!!$ +!!$ case('U') +!!$ if (.not.present(initu)) then +!!$ call psb_errpush(psb_err_internal_error_,name,& +!!$ & a_err='missing initu to smoother_apply') +!!$ goto 9999 +!!$ end if +!!$ call psb_geaxpby(sone,x,szero,tx,desc_data,info) +!!$ call psb_geaxpby(sone,initu,szero,ty,desc_data,info) +!!$ call psb_spmm(-sone,sm%pa,ty,sone,tx,desc_data,info,work=aux,trans=trans_) +!!$ call sm%sv%apply(sone,tx,szero,ty,desc_data,trans_,aux,wv(3:),info,init='Y') +!!$ +!!$ case default +!!$ call psb_errpush(psb_err_internal_error_,name,& +!!$ & a_err='wrong init to smoother_apply') +!!$ goto 9999 +!!$ end select +!!$ +!!$ do i=1, sweeps-1 +!!$ ! +!!$ ! Compute Y(j+1) = Y(j)+ D^(-1)*(X-A*Y(j)), +!!$ ! where is the diagonal and A the matrix. +!!$ ! +!!$ call psb_geaxpby(sone,x,szero,tx,desc_data,info) +!!$ call psb_spmm(-sone,sm%pa,ty,sone,tx,desc_data,info,work=aux,trans=trans_) +!!$ +!!$ if (info /= psb_success_) exit +!!$ +!!$ call sm%sv%apply(sone,tx,sone,ty,desc_data,trans_,aux,wv(3:),info,init='Y') +!!$ +!!$ if (info /= psb_success_) exit +!!$ +!!$ if ( sm%checkres.and.(mod(i,sm%checkiter) == 0) ) then +!!$ call psb_geaxpby(sone,x,szero,r,r,desc_data,info) +!!$ call psb_spmm(-sone,sm%pa,ty,sone,r,desc_data,info) +!!$ res = psb_genrm2(r,desc_data,info) +!!$ if( sm%printres ) then +!!$ call log_conv("BJAC",me,i,sm%printiter,res,resdenum,sm%tol) +!!$ end if +!!$ if ( res < sm%tol*resdenum ) then +!!$ if( (sm%printres).and.(mod(sm%printiter,sm%checkiter)/=0) ) & +!!$ & call log_conv("BJAC",me,i,1,res,resdenum,sm%tol) +!!$ exit +!!$ end if +!!$ end if +!!$ +!!$ end do +!!$ +!!$ +!!$ if (info == psb_success_) call psb_geaxpby(alpha,ty,beta,y,desc_data,info) +!!$ +!!$ if (info /= psb_success_) then +!!$ info=psb_err_internal_error_ +!!$ call psb_errpush(info,name,& +!!$ & a_err='subsolve with Jacobi sweeps > 1') +!!$ goto 9999 +!!$ end if +!!$ +!!$ +!!$ end associate +!!$ +!!$ class default +!!$ ! +!!$ ! +!!$ ! Apply multiple sweeps of a block-Jacobi solver +!!$ ! to compute an approximate solution of a linear system. +!!$ ! +!!$ ! +!!$ if (size(wv) < 2) then +!!$ info = psb_err_internal_error_ +!!$ call psb_errpush(info,name,& +!!$ & a_err='invalid wv size in smoother_apply') +!!$ goto 9999 +!!$ end if +!!$ associate(tx => wv(1), ty => wv(2)) +!!$ +!!$ ! +!!$ ! Unroll the first iteration and fold it inside SELECT CASE +!!$ ! this will save one AXPBY and one SPMM when INIT=Z, and will be +!!$ ! significant when sweeps=1 (a common case) +!!$ ! +!!$ select case (init_) +!!$ case('Z') +!!$ +!!$ call sm%sv%apply(sone,x,szero,ty,desc_data,trans_,aux,wv(3:),info,init='Z') +!!$ +!!$ case('Y') +!!$ call psb_geaxpby(sone,x,szero,tx,desc_data,info) +!!$ call psb_geaxpby(sone,y,szero,ty,desc_data,info) +!!$ call psb_spmm(-sone,sm%nd,ty,sone,tx,desc_data,info,work=aux,trans=trans_) +!!$ call sm%sv%apply(sone,tx,szero,ty,desc_data,trans_,aux,wv(3:),info,init='Y') +!!$ +!!$ case('U') +!!$ if (.not.present(initu)) then +!!$ call psb_errpush(psb_err_internal_error_,name,& +!!$ & a_err='missing initu to smoother_apply') +!!$ goto 9999 +!!$ end if +!!$ call psb_geaxpby(sone,x,szero,tx,desc_data,info) +!!$ call psb_geaxpby(sone,initu,szero,ty,desc_data,info) +!!$ call psb_spmm(-sone,sm%nd,ty,sone,tx,desc_data,info,work=aux,trans=trans_) +!!$ call sm%sv%apply(sone,tx,szero,ty,desc_data,trans_,aux,wv(3:),info,init='Y') +!!$ +!!$ case default +!!$ call psb_errpush(psb_err_internal_error_,name,& +!!$ & a_err='wrong init to smoother_apply') +!!$ goto 9999 +!!$ end select +!!$ +!!$ do i=1, sweeps-1 +!!$ ! +!!$ ! Compute Y(j+1) = D^(-1)*(X-ND*Y(j)), where D and ND are the +!!$ ! block diagonal part and the remaining part of the local matrix +!!$ ! and Y(j) is the approximate solution at sweep j. +!!$ ! +!!$ call psb_geaxpby(sone,x,szero,tx,desc_data,info) +!!$ call psb_spmm(-sone,sm%nd,ty,sone,tx,desc_data,info,work=aux,trans=trans_) +!!$ +!!$ if (info /= psb_success_) exit +!!$ +!!$ call sm%sv%apply(sone,tx,szero,ty,desc_data,trans_,aux,wv(3:),info,init='Y') +!!$ +!!$ if (info /= psb_success_) exit +!!$ +!!$ if ( sm%checkres.and.(mod(i,sm%checkiter) == 0) ) then +!!$ call psb_geaxpby(sone,x,szero,r,r,desc_data,info) +!!$ call psb_spmm(-sone,sm%pa,ty,sone,r,desc_data,info) +!!$ res = psb_genrm2(r,desc_data,info) +!!$ if( sm%printres ) then +!!$ call log_conv("BJAC",me,i,sm%printiter,res,resdenum,sm%tol) +!!$ end if +!!$ if (res < sm%tol*resdenum ) then +!!$ if( (sm%printres).and.( mod(sm%printiter,sm%checkiter) /=0 ) ) & +!!$ & call log_conv("BJAC",me,i,1,res,resdenum,sm%tol) +!!$ exit +!!$ end if +!!$ end if +!!$ +!!$ end do +!!$ +!!$ if (info == psb_success_) call psb_geaxpby(alpha,ty,beta,y,desc_data,info) +!!$ +!!$ if (info /= psb_success_) then +!!$ info=psb_err_internal_error_ +!!$ call psb_errpush(info,name,& +!!$ & a_err='subsolve with Jacobi sweeps > 1') +!!$ goto 9999 +!!$ end if +!!$ +!!$ end associate +!!$ end select +!!$ +!!$ else +!!$ +!!$ info = psb_err_iarg_neg_ +!!$ call psb_errpush(info,name,& +!!$ & i_err=(/itwo,sweeps,izero,izero,izero/)) +!!$ goto 9999 +!!$ +!!$ 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 + +9999 call psb_error_handler(err_act) + + return + +end subroutine amg_s_poly_smoother_apply_vect diff --git a/amgprec/impl/smoother/amg_s_poly_smoother_bld.f90 b/amgprec/impl/smoother/amg_s_poly_smoother_bld.f90 new file mode 100644 index 00000000..231136f1 --- /dev/null +++ b/amgprec/impl/smoother/amg_s_poly_smoother_bld.f90 @@ -0,0 +1,180 @@ +! +! +! AMG4PSBLAS version 1.0 +! Algebraic Multigrid Package +! based on PSBLAS (Parallel Sparse BLAS version 3.7) +! +! (C) Copyright 2021 +! +! Salvatore Filippone +! Pasqua D'Ambra +! Daniela di Serafino +! +! Redistribution and use in source and binary forms, with or without +! modification, are permitted provided that the following conditions +! are met: +! 1. Redistributions of source code must retain the above copyright +! notice, this list of conditions and the following disclaimer. +! 2. Redistributions in binary form must reproduce the above copyright +! notice, this list of conditions, and the following disclaimer in the +! documentation and/or other materials provided with the distribution. +! 3. The name of the AMG4PSBLAS group or the names of its contributors may +! not be used to endorse or promote products derived from this +! software without specific written permission. +! +! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +! ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED +! TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR +! PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE AMG4PSBLAS GROUP OR ITS CONTRIBUTORS +! BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +! CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +! SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +! INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +! CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +! ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +! POSSIBILITY OF SUCH DAMAGE. +! +! +subroutine amg_s_poly_smoother_bld(a,desc_a,sm,info,amold,vmold,imold) + + use psb_base_mod + use amg_s_diag_solver + use amg_s_l1_diag_solver + use amg_d_poly_coeff_mod + use amg_s_poly_smoother, amg_protect_name => amg_s_poly_smoother_bld + Implicit None + + ! Arguments + type(psb_sspmat_type), intent(in), target :: a + Type(psb_desc_type), Intent(inout) :: desc_a + class(amg_s_poly_smoother_type), intent(inout) :: sm + integer(psb_ipk_), intent(out) :: info + class(psb_s_base_sparse_mat), intent(in), optional :: amold + class(psb_s_base_vect_type), intent(in), optional :: vmold + class(psb_i_base_vect_type), intent(in), optional :: imold + ! Local variables + type(psb_sspmat_type) :: tmpa + integer(psb_ipk_) :: n_row,n_col, nrow_a, nztota, nzeros + type(psb_ctxt_type) :: ctxt + real(psb_spk_), 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 + + 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, me, np) + if (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),' start' + + + n_row = desc_a%get_local_rows() + n_col = desc_a%get_local_cols() + nrow_a = a%get_nrows() + nztota = a%get_nzeros() + if (.false.) then + select case(sm%variant) + case(amg_poly_lottes_) + ! do nothing + case(amg_poly_lottes_beta_) + if ((1<=sm%pdegree).and.(sm%pdegree<=30)) then + call psb_realloc(sm%pdegree,sm%poly_beta,info) + sm%poly_beta(1:sm%pdegree) = amg_d_poly_beta_mat(1:sm%pdegree,sm%pdegree) + else + info = psb_err_internal_error_ + call psb_errpush(info,name,& + & a_err='invalid sm%degree for poly_beta') + goto 9999 + end if + case(amg_poly_new_) + + if ((1<=sm%pdegree).and.(sm%pdegree<=30)) then + !Ok + sm%cf_a = amg_d_poly_a_vect(sm%pdegree) + else + info = psb_err_internal_error_ + call psb_errpush(info,name,& + & a_err='invalid sm%degree for poly_a') + goto 9999 + end if + + case default + info = psb_err_internal_error_ + call psb_errpush(info,name,& + & a_err='invalid sm%variant') + goto 9999 + end select + 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 (.false.) then +!!$ select type(ssv => sm%sv) +!!$ class is(amg_s_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 = sone +!!$ end select +!!$ else + if (sm%rho_ba <= szero) then + select case(sm%rho_estimate) + case(amg_poly_rho_est_power_) + block + type(psb_s_vect_type) :: tq, tt, tz,wv(2) + real(psb_spk_) :: znrm, lambda + real(psb_spk_),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(tt,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(sone) + call psb_geasb(tq,desc_a,info,mold=vmold) + call psb_spmm(sone,a,tq,szero,tt,desc_a,info) ! + call sm%sv%apply_v(sone,tt,szero,tz,desc_a,'NoTrans',work,wv,info) ! z_{k+1} = BA q_k + do i=1,sm%rho_estimate_iterations + znrm = psb_genrm2(tz,desc_a,info) ! znrm = |z_k|_2 + call psb_geaxpby((sone/znrm),tz,szero,tq,desc_a,info) ! q_k = z_k/znrm + call psb_spmm(sone,a,tq,szero,tt,desc_a,info) ! t_{k+1} = BA q_k + call sm%sv%apply_v(sone,tt,szero,tz,desc_a,'NoTrans',work,wv,info) ! z_{k+1} = B t_{k+1} + lambda = psb_gedot(tq,tz,desc_a,info) ! lambda = q_k^T z_{k+1} = q_k^T BA q_k + !write(0,*) 'BLD: lambda estimate ',i,lambda + end do + sm%rho_ba = lambda + end block + case default + write(0,*) ' Unknown algorithm for RHO(BA) estimate, defaulting to a value of 1.0 ' + sm%rho_ba = sone + end select + end if + + if (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),' end' + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + + return + +end subroutine amg_s_poly_smoother_bld diff --git a/amgprec/impl/smoother/amg_s_poly_smoother_clear_data.f90 b/amgprec/impl/smoother/amg_s_poly_smoother_clear_data.f90 new file mode 100644 index 00000000..5b0d88b7 --- /dev/null +++ b/amgprec/impl/smoother/amg_s_poly_smoother_clear_data.f90 @@ -0,0 +1,70 @@ +! +! +! AMG4PSBLAS version 1.0 +! Algebraic Multigrid Package +! based on PSBLAS (Parallel Sparse BLAS version 3.7) +! +! (C) Copyright 2021 +! +! Salvatore Filippone +! Pasqua D'Ambra +! Fabio Durastante +! +! Redistribution and use in source and binary forms, with or without +! modification, are permitted provided that the following conditions +! are met: +! 1. Redistributions of source code must retain the above copyright +! notice, this list of conditions and the following disclaimer. +! 2. Redistributions in binary form must reproduce the above copyright +! notice, this list of conditions, and the following disclaimer in the +! documentation and/or other materials provided with the distribution. +! 3. The name of the AMG4PSBLAS group or the names of its contributors may +! not be used to endorse or promote products derived from this +! software without specific written permission. +! +! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +! ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED +! TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR +! PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE AMG4PSBLAS GROUP OR ITS CONTRIBUTORS +! BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +! CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +! SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +! INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +! CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +! ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +! POSSIBILITY OF SUCH DAMAGE. +! +! +subroutine amg_s_poly_smoother_clear_data(sm,info) + + use psb_base_mod + use amg_s_poly_smoother, amg_protect_name => amg_s_poly_smoother_clear_data + Implicit None + ! Arguments + class(amg_s_poly_smoother_type), intent(inout) :: sm + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_) :: err_act + character(len=20) :: name='amg_s_poly_smoother_clear_data' + + call psb_erractionsave(err_act) + + info = 0 + sm%pdegree = 0 + if (allocated(sm%poly_beta)) deallocate(sm%poly_beta) + sm%pa => null() + if ((info==0).and.allocated(sm%sv)) then + call sm%sv%clear_data(info) + end if + if (info /= 0) then + info = psb_err_internal_error_ + call psb_errpush(info,name) + goto 9999 + end if + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + + return +end subroutine amg_s_poly_smoother_clear_data diff --git a/amgprec/impl/smoother/amg_s_poly_smoother_clone.f90 b/amgprec/impl/smoother/amg_s_poly_smoother_clone.f90 new file mode 100644 index 00000000..7bf3fbcc --- /dev/null +++ b/amgprec/impl/smoother/amg_s_poly_smoother_clone.f90 @@ -0,0 +1,90 @@ +! +! +! AMG4PSBLAS version 1.0 +! Algebraic Multigrid Package +! based on PSBLAS (Parallel Sparse BLAS version 3.7) +! +! (C) Copyright 2021 +! +! Salvatore Filippone +! Pasqua D'Ambra +! Daniela di Serafino +! +! Redistribution and use in source and binary forms, with or without +! modification, are permitted provided that the following conditions +! are met: +! 1. Redistributions of source code must retain the above copyright +! notice, this list of conditions and the following disclaimer. +! 2. Redistributions in binary form must reproduce the above copyright +! notice, this list of conditions, and the following disclaimer in the +! documentation and/or other materials provided with the distribution. +! 3. The name of the AMG4PSBLAS group or the names of its contributors may +! not be used to endorse or promote products derived from this +! software without specific written permission. +! +! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +! ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED +! TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR +! PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE AMG4PSBLAS GROUP OR ITS CONTRIBUTORS +! BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +! CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +! SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +! INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +! CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +! ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +! POSSIBILITY OF SUCH DAMAGE. +! +! +subroutine amg_s_poly_smoother_clone(sm,smout,info) + + use psb_base_mod + use amg_s_poly_smoother, amg_protect_name => amg_s_poly_smoother_clone + + Implicit None + + ! Arguments + class(amg_s_poly_smoother_type), intent(inout) :: sm + class(amg_s_base_smoother_type), allocatable, intent(inout) :: smout + integer(psb_ipk_), intent(out) :: info + ! Local variables + integer(psb_ipk_) :: err_act + + + info=psb_success_ + call psb_erractionsave(err_act) + + if (allocated(smout)) then + call smout%free(info) + if (info == psb_success_) deallocate(smout, stat=info) + end if + if (info == psb_success_) & + & allocate(amg_s_poly_smoother_type :: smout, stat=info) + if (info /= 0) then + info = psb_err_alloc_dealloc_ + goto 9999 + end if + + select type(smo => smout) + type is (amg_s_poly_smoother_type) + smo%pdegree = sm%pdegree + smo%rho_ba = sm%rho_ba + smo%poly_beta = sm%poly_beta + smo%pa => sm%pa + if ((info==psb_success_).and.(allocated(sm%sv))) then + allocate(smout%sv,mold=sm%sv,stat=info) + if (info == psb_success_) call sm%sv%clone(smo%sv,info) + end if + + class default + info = psb_err_internal_error_ + end select + + if (info /= 0) goto 9999 + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + + return +end subroutine amg_s_poly_smoother_clone diff --git a/amgprec/impl/smoother/amg_s_poly_smoother_clone_settings.f90 b/amgprec/impl/smoother/amg_s_poly_smoother_clone_settings.f90 new file mode 100644 index 00000000..ddbad88f --- /dev/null +++ b/amgprec/impl/smoother/amg_s_poly_smoother_clone_settings.f90 @@ -0,0 +1,102 @@ +! +! +! AMG4PSBLAS version 1.0 +! Algebraic Multigrid Package +! asd on PSBLAS (Parallel Sparse BLAS version 3.7) +! +! (C) Copyright 2021 +! +! Salvatore Filippone +! Pasqua D'Ambra +! Fabio Durastante +! +! Redistribution and use in source and binary forms, with or without +! modification, are permitted provided that the following conditions +! are met: +! 1. Redistributions of source code must retain the above copyright +! notice, this list of conditions and the following disclaimer. +! 2. Redistributions in binary form must reproduce the above copyright +! notice, this list of conditions, and the following disclaimer in the +! documentation and/or other materials provided with the distribution. +! 3. The name of the AMG4PSBLAS group or the names of its contributors may +! not be used to endorse or promote products derived from this +! software without specific written permission. +! +! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +! ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED +! TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR +! PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE AMG4PSBLAS GROUP OR ITS CONTRIBUTORS +! BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +! CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +! SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +! INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +! CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +! ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +! POSSIBILITY OF SUCH DAMAGE. +! +! +subroutine amg_s_poly_smoother_clone_settings(sm,smout,info) + + use psb_base_mod + use amg_s_poly_smoother, amg_protect_name => amg_s_poly_smoother_clone_settings + Implicit None + ! Arguments + class(amg_s_poly_smoother_type), intent(inout) :: sm + class(amg_s_base_smoother_type), intent(inout) :: smout + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_) :: err_act + character(len=20) :: name='d_poly_smoother_clone_settings' + + call psb_erractionsave(err_act) + + info = psb_success_ + + select type(smout) + class is(amg_s_poly_smoother_type) + + smout%pa => null() + smout%pdegree = sm%pdegree + smout%variant = sm%variant + smout%cf_a = sm%cf_a + smout%rho_ba = sm%rho_ba + smout%rho_estimate = sm%rho_estimate + smout%rho_estimate_iterations = sm%rho_estimate_iterations + smout%poly_beta = sm%poly_beta + + if (allocated(smout%sv)) then + if (.not.same_type_as(sm%sv,smout%sv)) then + call smout%sv%free(info) + if (info == 0) deallocate(smout%sv,stat=info) + end if + end if + if (info /= 0) then + info = psb_err_internal_error_ + else + if (allocated(smout%sv)) then + if (same_type_as(sm%sv,smout%sv)) then + call sm%sv%clone_settings(smout%sv,info) + else + info = psb_err_internal_error_ + end if + else + allocate(smout%sv,mold=sm%sv,stat=info) + if (info == 0) call sm%sv%clone_settings(smout%sv,info) + if (info /= 0) info = psb_err_internal_error_ + end if + end if + class default + info = psb_err_internal_error_ + end select + + if (info /= 0) then + call psb_errpush(info,name) + goto 9999 + end if + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + + return +end subroutine amg_s_poly_smoother_clone_settings diff --git a/amgprec/impl/smoother/amg_s_poly_smoother_cnv.f90 b/amgprec/impl/smoother/amg_s_poly_smoother_cnv.f90 new file mode 100644 index 00000000..810809fb --- /dev/null +++ b/amgprec/impl/smoother/amg_s_poly_smoother_cnv.f90 @@ -0,0 +1,77 @@ +! +! +! AMG4PSBLAS version 1.0 +! Algebraic Multigrid Package +! based on PSBLAS (Parallel Sparse BLAS version 3.7) +! +! (C) Copyright 2021 +! +! Salvatore Filippone +! Pasqua D'Ambra +! Fabio Durastante +! +! Redistribution and use in source and binary forms, with or without +! modification, are permitted provided that the following conditions +! are met: +! 1. Redistributions of source code must retain the above copyright +! notice, this list of conditions and the following disclaimer. +! 2. Redistributions in binary form must reproduce the above copyright +! notice, this list of conditions, and the following disclaimer in the +! documentation and/or other materials provided with the distribution. +! 3. The name of the AMG4PSBLAS group or the names of its contributors may +! not be used to endorse or promote products derived from this +! software without specific written permission. +! +! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +! ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED +! TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR +! PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE AMG4PSBLAS GROUP OR ITS CONTRIBUTORS +! BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +! CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +! SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +! INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +! CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +! ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +! POSSIBILITY OF SUCH DAMAGE. +! +! +subroutine amg_s_poly_smoother_cnv(sm,info,amold,vmold,imold) + + use psb_base_mod + use amg_s_diag_solver + use amg_s_poly_smoother, amg_protect_name => amg_s_poly_smoother_cnv + Implicit None + + ! Arguments + class(amg_s_poly_smoother_type), intent(inout) :: sm + integer(psb_ipk_), intent(out) :: info + class(psb_s_base_sparse_mat), intent(in), optional :: amold + class(psb_s_base_vect_type), intent(in), optional :: vmold + class(psb_i_base_vect_type), intent(in), optional :: imold + ! Local variables + integer(psb_ipk_) :: i, err_act, debug_unit, debug_level + character(len=20) :: name='d_poly_smoother_cnv', ch_err + + info=psb_success_ + call psb_erractionsave(err_act) + debug_unit = psb_get_debug_unit() + debug_level = psb_get_debug_level() + + + + if (allocated(sm%sv)) & + & call sm%sv%cnv(info,amold=amold,vmold=vmold,imold=imold) + if (info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,& + & a_err='solver cnv') + goto 9999 + end if + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + + return + +end subroutine amg_s_poly_smoother_cnv diff --git a/amgprec/impl/smoother/amg_s_poly_smoother_csetc.f90 b/amgprec/impl/smoother/amg_s_poly_smoother_csetc.f90 new file mode 100644 index 00000000..b00d1d0d --- /dev/null +++ b/amgprec/impl/smoother/amg_s_poly_smoother_csetc.f90 @@ -0,0 +1,76 @@ +! +! +! AMG4PSBLAS version 1.0 +! Algebraic Multigrid Package +! based on PSBLAS (Parallel Sparse BLAS version 3.7) +! +! (C) Copyright 2021 +! +! Salvatore Filippone +! Pasqua D'Ambra +! Daniela di Serafino +! +! Redistribution and use in source and binary forms, with or without +! modification, are permitted provided that the following conditions +! are met: +! 1. Redistributions of source code must retain the above copyright +! notice, this list of conditions and the following disclaimer. +! 2. Redistributions in binary form must reproduce the above copyright +! notice, this list of conditions, and the following disclaimer in the +! documentation and/or other materials provided with the distribution. +! 3. The name of the AMG4PSBLAS group or the names of its contributors may +! not be used to endorse or promote products derived from this +! software without specific written permission. +! +! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +! ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED +! TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR +! PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE AMG4PSBLAS GROUP OR ITS CONTRIBUTORS +! BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +! CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +! SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +! INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +! CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +! ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +! POSSIBILITY OF SUCH DAMAGE. +! +! +subroutine amg_s_poly_smoother_csetc(sm,what,val,info,idx) + + use psb_base_mod + use amg_s_poly_smoother, amg_protect_name => amg_s_poly_smoother_csetc + Implicit None + ! Arguments + class(amg_s_poly_smoother_type), intent(inout) :: sm + 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, ival + character(len=20) :: name='d_poly_smoother_csetc' + + info = psb_success_ + call psb_erractionsave(err_act) + + select case(psb_toupper(trim(what))) + case('POLY_VARIANT') + call sm%set(what,amg_stringval(val),info,idx=idx) + case('POLY_RHO_ESTIMATE') + call sm%set(what,amg_stringval(val),info,idx=idx) + case default + call sm%amg_s_base_smoother_type%set(what,val,info,idx=idx) + end select + + if (info /= psb_success_) then + info = psb_err_from_subroutine_ + call psb_errpush(info, name) + goto 9999 + end if + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + + return +end subroutine amg_s_poly_smoother_csetc diff --git a/amgprec/impl/smoother/amg_s_poly_smoother_cseti.f90 b/amgprec/impl/smoother/amg_s_poly_smoother_cseti.f90 new file mode 100644 index 00000000..78ce1b39 --- /dev/null +++ b/amgprec/impl/smoother/amg_s_poly_smoother_cseti.f90 @@ -0,0 +1,92 @@ +! +! +! AMG4PSBLAS version 1.0 +! Algebraic Multigrid Package +! based on PSBLAS (Parallel Sparse BLAS version 3.7) +! +! (C) Copyright 2021 +! +! Salvatore Filippone +! Pasqua D'Ambra +! Daniela di Serafino +! +! Redistribution and use in source and binary forms, with or without +! modification, are permitted provided that the following conditions +! are met: +! 1. Redistributions of source code must retain the above copyright +! notice, this list of conditions and the following disclaimer. +! 2. Redistributions in binary form must reproduce the above copyright +! notice, this list of conditions, and the following disclaimer in the +! documentation and/or other materials provided with the distribution. +! 3. The name of the AMG4PSBLAS group or the names of its contributors may +! not be used to endorse or promote products derived from this +! software without specific written permission. +! +! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +! ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED +! TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR +! PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE AMG4PSBLAS GROUP OR ITS CONTRIBUTORS +! BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +! CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +! SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +! INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +! CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +! ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +! POSSIBILITY OF SUCH DAMAGE. +! +! +subroutine amg_s_poly_smoother_cseti(sm,what,val,info,idx) + + use psb_base_mod + use amg_s_poly_smoother, amg_protect_nam => amg_s_poly_smoother_cseti + Implicit None + + ! Arguments + class(amg_s_poly_smoother_type), intent(inout) :: sm + character(len=*), intent(in) :: what + integer(psb_ipk_), intent(in) :: val + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_), intent(in), optional :: idx + integer(psb_ipk_) :: err_act + character(len=20) :: name='d_poly_smoother_cseti' + + info = psb_success_ + call psb_erractionsave(err_act) + + select case(psb_toupper(what)) + 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_',val + sm%variant = amg_poly_lottes_ + end select + case('POLY_RHO_ESTIMATE') + select case(val) + case (amg_poly_rho_est_power_) + sm%rho_estimate = val + case default + write(0,*) 'Invalid choice for POLY_RHO_ESTIMATE, defaulting to amg_poly_rho_power' + sm%variant = amg_poly_rho_est_power_ + end select + case('POLY_RHO_ESTIMATE_ITERATIONS') + if (val>0) then + sm%rho_estimate_iterations = val + else + write(0,*) 'Invalid choice for POLY_RHO_ESTIMATE_ITERATIONS, defaulting to 20' + sm%variant = 20 + end if + case default + call sm%amg_s_base_smoother_type%set(what,val,info,idx=idx) + end select + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + + return +end subroutine amg_s_poly_smoother_cseti diff --git a/amgprec/impl/smoother/amg_s_poly_smoother_csetr.f90 b/amgprec/impl/smoother/amg_s_poly_smoother_csetr.f90 new file mode 100644 index 00000000..c7ec74cd --- /dev/null +++ b/amgprec/impl/smoother/amg_s_poly_smoother_csetr.f90 @@ -0,0 +1,74 @@ +! +! +! AMG4PSBLAS version 1.0 +! Algebraic Multigrid Package +! based on PSBLAS (Parallel Sparse BLAS version 3.7) +! +! (C) Copyright 2021 +! +! Salvatore Filippone +! Pasqua D'Ambra +! Daniela di Serafino +! +! Redistribution and use in source and binary forms, with or without +! modification, are permitted provided that the following conditions +! are met: +! 1. Redistributions of source code must retain the above copyright +! notice, this list of conditions and the following disclaimer. +! 2. Redistributions in binary form must reproduce the above copyright +! notice, this list of conditions, and the following disclaimer in the +! documentation and/or other materials provided with the distribution. +! 3. The name of the AMG4PSBLAS group or the names of its contributors may +! not be used to endorse or promote products derived from this +! software without specific written permission. +! +! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +! ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED +! TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR +! PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE AMG4PSBLAS GROUP OR ITS CONTRIBUTORS +! BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +! CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +! SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +! INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +! CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +! ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +! POSSIBILITY OF SUCH DAMAGE. +! +! +subroutine amg_s_poly_smoother_csetr(sm,what,val,info,idx) + + use psb_base_mod + use amg_s_poly_smoother, amg_protect_nam => amg_s_poly_smoother_csetr + Implicit None + + ! Arguments + class(amg_s_poly_smoother_type), intent(inout) :: sm + character(len=*), intent(in) :: what + real(psb_spk_), intent(in) :: val + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_), intent(in), optional :: idx + integer(psb_ipk_) :: err_act + character(len=20) :: name='d_poly_smoother_csetr' + + info = psb_success_ + call psb_erractionsave(err_act) + + select case(psb_toupper(what)) + case('POLY_RHO_BA') + if ((szero amg_s_poly_smoother_descr + + Implicit None + + ! Arguments + class(amg_s_poly_smoother_type), intent(in) :: sm + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_), intent(in), optional :: iout + logical, intent(in), optional :: coarse + character(len=*), intent(in), optional :: prefix + + ! Local variables + integer(psb_ipk_) :: err_act + character(len=20), parameter :: name='amg_s_poly_smoother_descr' + integer(psb_ipk_) :: iout_ + logical :: coarse_ + character(1024) :: prefix_ + + call psb_erractionsave(err_act) + info = psb_success_ + if (present(coarse)) then + coarse_ = coarse + else + coarse_ = .false. + end if + if (present(iout)) then + iout_ = iout + else + iout_ = psb_out_unit + endif + if (present(prefix)) then + prefix_ = prefix + else + prefix_ = "" + end if + + write(iout_,*) trim(prefix_), ' Polynomial smoother ' + select case(sm%variant) + case(amg_poly_lottes_) + write(iout_,*) trim(prefix_), ' variant: ','POLY_LOTTES' + !write(iout_,*) trim(prefix_), ' Degree: ',sm%pdegree + write(iout_,*) trim(prefix_), ' rho_ba: ',sm%rho_ba + case(amg_poly_lottes_beta_) + write(iout_,*) trim(prefix_), ' variant: ','POLY_LOTTES_BETA' + !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) + case(amg_poly_new_) + write(iout_,*) trim(prefix_), ' variant: ','POLY_NEW' + !write(iout_,*) trim(prefix_), ' Degree: ',sm%pdegree + write(iout_,*) trim(prefix_), ' rho_ba: ',sm%rho_ba + !write(iout_,*) trim(prefix_), ' Coefficient: ',sm%cf_a + case default + write(iout_,*) trim(prefix_), ' variant: ','UNKNOWN???' + end select + if (allocated(sm%sv)) then + write(iout_,*) trim(prefix_), ' Local solver details:' + call sm%sv%descr(info,iout_,coarse=coarse,prefix=prefix) + end if + + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + return +end subroutine amg_s_poly_smoother_descr diff --git a/amgprec/impl/smoother/amg_s_poly_smoother_dmp.f90 b/amgprec/impl/smoother/amg_s_poly_smoother_dmp.f90 new file mode 100644 index 00000000..f6fa2f8a --- /dev/null +++ b/amgprec/impl/smoother/amg_s_poly_smoother_dmp.f90 @@ -0,0 +1,90 @@ +! +! +! AMG4PSBLAS version 1.0 +! Algebraic Multigrid Package +! based on PSBLAS (Parallel Sparse BLAS version 3.7) +! +! (C) Copyright 2021 +! +! Salvatore Filippone +! Pasqua D'Ambra +! Fabio Durastante +! +! Redistribution and use in source and binary forms, with or without +! modification, are permitted provided that the following conditions +! are met: +! 1. Redistributions of source code must retain the above copyright +! notice, this list of conditions and the following disclaimer. +! 2. Redistributions in binary form must reproduce the above copyright +! notice, this list of conditions, and the following disclaimer in the +! documentation and/or other materials provided with the distribution. +! 3. The name of the AMG4PSBLAS group or the names of its contributors may +! not be used to endorse or promote products derived from this +! software without specific written permission. +! +! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +! ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED +! TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR +! PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE AMG4PSBLAS GROUP OR ITS CONTRIBUTORS +! BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +! CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +! SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +! INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +! CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +! ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +! POSSIBILITY OF SUCH DAMAGE. +! +! +subroutine amg_s_poly_smoother_dmp(sm,desc,level,info,prefix,head,smoother,solver,global_num) + + use psb_base_mod + use amg_s_poly_smoother, amg_protect_nam => amg_s_poly_smoother_dmp + implicit none + class(amg_s_poly_smoother_type), intent(in) :: sm + type(psb_desc_type), intent(in) :: desc + integer(psb_ipk_), intent(in) :: level + integer(psb_ipk_), intent(out) :: info + character(len=*), intent(in), optional :: prefix, head + logical, optional, intent(in) :: smoother, solver, global_num + integer(psb_ipk_) :: i, j, il1, iln, lname, lev + type(psb_ctxt_type) :: ctxt + integer(psb_ipk_) :: iam, np + character(len=80) :: prefix_ + character(len=120) :: fname ! len should be at least 20 more than + integer(psb_lpk_), allocatable :: iv(:) + logical :: smoother_, global_num_ + ! len of prefix_ + + info = 0 + + if (present(prefix)) then + prefix_ = trim(prefix(1:min(len(prefix),len(prefix_)))) + else + prefix_ = "dump_smth_d" + end if + ctxt = desc%get_context() + call psb_info(ctxt,iam,np) + + if (present(smoother)) then + smoother_ = smoother + else + smoother_ = .false. + end if + if (present(global_num)) then + global_num_ = global_num + else + global_num_ = .false. + end if + lname = len_trim(prefix_) + fname = trim(prefix_) + write(fname(lname+1:lname+5),'(a,i3.3)') '_poly',iam + lname = lname + 8 + ! to be completed + + + + ! At base level do nothing for the smoother + if (allocated(sm%sv)) & + & call sm%sv%dump(desc,level,info,solver=solver,prefix=prefix,global_num=global_num) + +end subroutine amg_s_poly_smoother_dmp diff --git a/amgprec/impl/smoother/amg_z_jac_smoother_apply_vect.f90 b/amgprec/impl/smoother/amg_z_jac_smoother_apply_vect.f90 index ac3aaf3e..16d2a484 100644 --- a/amgprec/impl/smoother/amg_z_jac_smoother_apply_vect.f90 +++ b/amgprec/impl/smoother/amg_z_jac_smoother_apply_vect.f90 @@ -175,7 +175,7 @@ subroutine amg_z_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,& if (info /= psb_success_) exit if ( sm%checkres.and.(mod(i,sm%checkiter) == 0) ) then - call psb_geaxpby(zone,x,zzero,r,r,desc_data,info) + call psb_geaxpby(zone,x,zzero,r,desc_data,info) call psb_spmm(-zone,sm%pa,ty,zone,r,desc_data,info) res = psb_genrm2(r,desc_data,info) if( sm%printres ) then diff --git a/samples/advanced/pdegen/Makefile b/samples/advanced/pdegen/Makefile index b5092a22..fc9a7f21 100644 --- a/samples/advanced/pdegen/Makefile +++ b/samples/advanced/pdegen/Makefile @@ -8,23 +8,27 @@ FINCLUDES=$(FMFLAG). $(FMFLAG)$(AMGMODDIR) $(FMFLAG)$(AMGINCDIR) $(PSBLAS_INCLUD LINKOPT= EXEDIR=./runs +DGEN2D=amg_d_pde2d_base_mod.o amg_d_pde2d_exp_mod.o amg_d_pde2d_gauss_mod.o amg_d_pde2d_box_mod.o +DGEN3D=amg_d_pde3d_base_mod.o amg_d_pde3d_exp_mod.o amg_d_pde3d_gauss_mod.o amg_d_pde3d_box_mod.o +SGEN2D=amg_s_pde2d_base_mod.o amg_s_pde2d_exp_mod.o amg_s_pde2d_gauss_mod.o amg_s_pde2d_box_mod.o +SGEN3D=amg_s_pde3d_base_mod.o amg_s_pde3d_exp_mod.o amg_s_pde3d_gauss_mod.o amg_s_pde3d_box_mod.o all: amg_s_pde3d amg_d_pde3d amg_s_pde2d amg_d_pde2d -amg_d_pde3d: amg_d_pde3d.o amg_d_genpde_mod.o amg_d_pde3d_base_mod.o amg_d_pde3d_exp_mod.o amg_d_pde3d_gauss_mod.o data_input.o - $(FLINK) $(LINKOPT) amg_d_pde3d.o amg_d_genpde_mod.o amg_d_pde3d_base_mod.o amg_d_pde3d_exp_mod.o amg_d_pde3d_gauss_mod.o data_input.o -o amg_d_pde3d $(AMG_LIBS) $(PSBLAS_LIBS) $(LDLIBS) +amg_d_pde3d: amg_d_pde3d.o amg_d_genpde_mod.o $(DGEN3D) data_input.o + $(FLINK) $(LINKOPT) amg_d_pde3d.o amg_d_genpde_mod.o $(DGEN3D) data_input.o -o amg_d_pde3d $(AMG_LIBS) $(PSBLAS_LIBS) $(LDLIBS) /bin/mv amg_d_pde3d $(EXEDIR) -amg_s_pde3d: amg_s_pde3d.o amg_s_genpde_mod.o amg_s_pde3d_base_mod.o amg_s_pde3d_exp_mod.o amg_s_pde3d_gauss_mod.o data_input.o - $(FLINK) $(LINKOPT) amg_s_pde3d.o amg_s_genpde_mod.o amg_s_pde3d_base_mod.o amg_s_pde3d_exp_mod.o amg_s_pde3d_gauss_mod.o data_input.o -o amg_s_pde3d $(AMG_LIBS) $(PSBLAS_LIBS) $(LDLIBS) +amg_s_pde3d: amg_s_pde3d.o amg_s_genpde_mod.o $(SGEN3D) data_input.o + $(FLINK) $(LINKOPT) amg_s_pde3d.o amg_s_genpde_mod.o $(SGEN3D) data_input.o -o amg_s_pde3d $(AMG_LIBS) $(PSBLAS_LIBS) $(LDLIBS) /bin/mv amg_s_pde3d $(EXEDIR) -amg_d_pde2d: amg_d_pde2d.o amg_d_genpde_mod.o amg_d_pde2d_base_mod.o amg_d_pde2d_exp_mod.o amg_d_pde2d_box_mod.o data_input.o - $(FLINK) $(LINKOPT) amg_d_pde2d.o amg_d_genpde_mod.o amg_d_pde2d_base_mod.o amg_d_pde2d_exp_mod.o amg_d_pde2d_box_mod.o data_input.o -o amg_d_pde2d $(AMG_LIBS) $(PSBLAS_LIBS) $(LDLIBS) +amg_d_pde2d: amg_d_pde2d.o amg_d_genpde_mod.o $(DGEN2D) data_input.o + $(FLINK) $(LINKOPT) amg_d_pde2d.o amg_d_genpde_mod.o $(DGEN2D) data_input.o -o amg_d_pde2d $(AMG_LIBS) $(PSBLAS_LIBS) $(LDLIBS) /bin/mv amg_d_pde2d $(EXEDIR) -amg_s_pde2d: amg_s_pde2d.o amg_s_genpde_mod.o amg_s_pde2d_base_mod.o amg_s_pde2d_exp_mod.o amg_s_pde2d_box_mod.o data_input.o - $(FLINK) $(LINKOPT) amg_s_pde2d.o amg_s_genpde_mod.o amg_s_pde2d_base_mod.o amg_s_pde2d_exp_mod.o amg_s_pde2d_box_mod.o data_input.o -o amg_s_pde2d $(AMG_LIBS) $(PSBLAS_LIBS) $(LDLIBS) +amg_s_pde2d: amg_s_pde2d.o amg_s_genpde_mod.o $(SGEN2D) data_input.o + $(FLINK) $(LINKOPT) amg_s_pde2d.o amg_s_genpde_mod.o $(SGEN2D) data_input.o -o amg_s_pde2d $(AMG_LIBS) $(PSBLAS_LIBS) $(LDLIBS) /bin/mv amg_s_pde2d $(EXEDIR) amg_d_pde3d_rebld: amg_d_pde3d_rebld.o data_input.o @@ -33,10 +37,10 @@ amg_d_pde3d_rebld: amg_d_pde3d_rebld.o data_input.o amg_d_pde3d.o amg_s_pde3d.o amg_d_pde2d.o amg_s_pde2d.o: data_input.o -amg_d_pde3d.o: amg_d_genpde_mod.o amg_d_pde3d_base_mod.o amg_d_pde3d_exp_mod.o amg_d_pde3d_gauss_mod.o -amg_s_pde3d.o: amg_s_genpde_mod.o amg_s_pde3d_base_mod.o amg_s_pde3d_exp_mod.o amg_s_pde3d_gauss_mod.o -amg_d_pde2d.o: amg_d_genpde_mod.o amg_d_pde2d_base_mod.o amg_d_pde2d_exp_mod.o amg_d_pde2d_box_mod.o -amg_s_pde2d.o: amg_s_genpde_mod.o amg_s_pde2d_base_mod.o amg_s_pde2d_exp_mod.o amg_s_pde2d_box_mod.o +amg_d_pde3d.o: amg_d_genpde_mod.o $(DGEN3D) +amg_s_pde3d.o: amg_s_genpde_mod.o $(SGEN3D) +amg_d_pde2d.o: amg_d_genpde_mod.o $(DGEN2D) +amg_s_pde2d.o: amg_s_genpde_mod.o $(SGEN2D) check: all cd runs && ./amg_d_pde2d 0)& - & call prec%set('min_coarse_size_per_process', p_choice%csizepp, info) + & call prec%set('min_coarse_size_per_process', p_choice%csizepp, info) if (p_choice%mncrratio>1)& & call prec%set('min_cr_ratio', p_choice%mncrratio, info) if (p_choice%maxlevs>0)& @@ -332,7 +342,9 @@ program amg_d_pde2d call prec%set('smoother_type', p_choice%smther, info) call prec%set('smoother_sweeps', p_choice%jsweeps, info) - + call prec%set('poly_degree', p_choice%degree, info) + call prec%set('poly_variant', p_choice%pvariant, info) + select case (psb_toupper(p_choice%smther)) case ('GS','BWGS','FBGS','JACOBI','L1-JACOBI','L1-FBGS') ! do nothing @@ -362,6 +374,8 @@ program amg_d_pde2d 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('poly_degree', p_choice%degree2, info,pos='post') + call prec%set('poly_variant', p_choice%pvariant2, info,pos='post') select case (psb_toupper(p_choice%smther2)) case ('GS','BWGS','FBGS','JACOBI','L1-JACOBI','L1-FBGS') ! do nothing @@ -577,6 +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%pvariant,inp_unit) ! 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 @@ -589,11 +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%pvariant2,inp_unit) ! 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 @@ -659,6 +677,8 @@ contains ! broadcast first (pre-)smoother / 1-lev prec data call psb_bcast(ctxt,prec%smther) call psb_bcast(ctxt,prec%jsweeps) + call psb_bcast(ctxt,prec%degree) + call psb_bcast(ctxt,prec%pvariant) call psb_bcast(ctxt,prec%novr) call psb_bcast(ctxt,prec%restr) call psb_bcast(ctxt,prec%prol) @@ -671,6 +691,8 @@ contains ! broadcast second (post-)smoother call psb_bcast(ctxt,prec%smther2) call psb_bcast(ctxt,prec%jsweeps2) + call psb_bcast(ctxt,prec%degree2) + call psb_bcast(ctxt,prec%pvariant2) call psb_bcast(ctxt,prec%novr2) call psb_bcast(ctxt,prec%restr2) call psb_bcast(ctxt,prec%prol2) diff --git a/samples/advanced/pdegen/amg_d_pde2d_base_mod.f90 b/samples/advanced/pdegen/amg_d_pde2d_base_mod.f90 index a406e90e..e6613370 100644 --- a/samples/advanced/pdegen/amg_d_pde2d_base_mod.f90 +++ b/samples/advanced/pdegen/amg_d_pde2d_base_mod.f90 @@ -38,52 +38,52 @@ module amg_d_pde2d_base_mod use psb_base_mod, only : psb_dpk_, dzero, done real(psb_dpk_), save, private :: epsilon=done/80 contains - subroutine pde_set_parm(dat) + subroutine pde_set_parm2d_base(dat) real(psb_dpk_), intent(in) :: dat epsilon = dat - end subroutine pde_set_parm + end subroutine pde_set_parm2d_base ! ! functions parametrizing the differential equation ! - function b1(x,y) + function b1_base(x,y) implicit none - real(psb_dpk_) :: b1 + real(psb_dpk_) :: b1_base real(psb_dpk_), intent(in) :: x,y - b1 = dzero/1.414_psb_dpk_ - end function b1 - function b2(x,y) + b1_base = dzero/1.414_psb_dpk_ + end function b1_base + function b2_base(x,y) implicit none - real(psb_dpk_) :: b2 + real(psb_dpk_) :: b2_base real(psb_dpk_), intent(in) :: x,y - b2 = dzero/1.414_psb_dpk_ - end function b2 - function c(x,y) + b2_base = dzero/1.414_psb_dpk_ + end function b2_base + function c_base(x,y) implicit none - real(psb_dpk_) :: c + real(psb_dpk_) :: c_base real(psb_dpk_), intent(in) :: x,y - c = dzero - end function c - function a1(x,y) + c_base = dzero + end function c_base + function a1_base(x,y) implicit none - real(psb_dpk_) :: a1 + real(psb_dpk_) :: a1_base real(psb_dpk_), intent(in) :: x,y - a1=done*epsilon - end function a1 - function a2(x,y) + a1_base=done*epsilon + end function a1_base + function a2_base(x,y) implicit none - real(psb_dpk_) :: a2 + real(psb_dpk_) :: a2_base real(psb_dpk_), intent(in) :: x,y - a2=done*epsilon - end function a2 - function g(x,y) + a2_base=done*epsilon + end function a2_base + function g_base(x,y) implicit none - real(psb_dpk_) :: g + real(psb_dpk_) :: g_base real(psb_dpk_), intent(in) :: x,y - g = dzero + g_base = dzero if (x == done) then - g = done + g_base = done else if (x == dzero) then - g = done + g_base = done end if - end function g + end function g_base end module amg_d_pde2d_base_mod diff --git a/samples/advanced/pdegen/amg_d_pde2d_box_mod.f90 b/samples/advanced/pdegen/amg_d_pde2d_box_mod.f90 index db743633..11d5770b 100644 --- a/samples/advanced/pdegen/amg_d_pde2d_box_mod.f90 +++ b/samples/advanced/pdegen/amg_d_pde2d_box_mod.f90 @@ -38,10 +38,10 @@ module amg_d_pde2d_box_mod use psb_base_mod, only : psb_dpk_, dzero, done real(psb_dpk_), save, private :: epsilon=done/80 contains - subroutine pde_set_parm(dat) + subroutine pde_set_parm2d_box(dat) real(psb_dpk_), intent(in) :: dat epsilon = dat - end subroutine pde_set_parm + end subroutine pde_set_parm2d_box ! ! functions parametrizing the differential equation ! diff --git a/samples/advanced/pdegen/amg_d_pde2d_exp_mod.f90 b/samples/advanced/pdegen/amg_d_pde2d_exp_mod.f90 index 5dab37bc..76a733fc 100644 --- a/samples/advanced/pdegen/amg_d_pde2d_exp_mod.f90 +++ b/samples/advanced/pdegen/amg_d_pde2d_exp_mod.f90 @@ -38,10 +38,10 @@ module amg_d_pde2d_exp_mod use psb_base_mod, only : psb_dpk_, done, dzero real(psb_dpk_), save, private :: epsilon=done/80 contains - subroutine pde_set_parm(dat) + subroutine pde_set_parm2d_exp(dat) real(psb_dpk_), intent(in) :: dat epsilon = dat - end subroutine pde_set_parm + end subroutine pde_set_parm2d_exp ! ! functions parametrizing the differential equation ! diff --git a/samples/advanced/pdegen/amg_d_pde2d_gauss_mod.f90 b/samples/advanced/pdegen/amg_d_pde2d_gauss_mod.f90 new file mode 100644 index 00000000..93b911ea --- /dev/null +++ b/samples/advanced/pdegen/amg_d_pde2d_gauss_mod.f90 @@ -0,0 +1,89 @@ +! +! +! AMG4PSBLAS version 1.0 +! Algebraic Multigrid Package +! based on PSBLAS (Parallel Sparse BLAS version 3.7) +! +! (C) Copyright 2021 +! +! Salvatore Filippone +! Pasqua D'Ambra +! Fabio Durastante +! +! Redistribution and use in source and binary forms, with or without +! modification, are permitted provided that the following conditions +! are met: +! 1. Redistributions of source code must retain the above copyright +! notice, this list of conditions and the following disclaimer. +! 2. Redistributions in binary form must reproduce the above copyright +! notice, this list of conditions, and the following disclaimer in the +! documentation and/or other materials provided with the distribution. +! 3. The name of the AMG4PSBLAS group or the names of its contributors may +! not be used to endorse or promote products derived from this +! software without specific written permission. +! +! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +! ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED +! TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR +! PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE AMG4PSBLAS GROUP OR ITS CONTRIBUTORS +! BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +! CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +! SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +! INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +! CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +! ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +! POSSIBILITY OF SUCH DAMAGE. +! +module amg_d_pde2d_gauss_mod + use psb_base_mod, only : psb_dpk_, done, dzero + real(psb_dpk_), save, private :: epsilon=done/80 +contains + subroutine pde_set_parm2d_gauss(dat) + real(psb_dpk_), intent(in) :: dat + epsilon = dat + end subroutine pde_set_parm2d_gauss + ! + ! functions parametrizing the differential equation + ! + function b1_gauss(x,y) + implicit none + real(psb_dpk_) :: b1_gauss + real(psb_dpk_), intent(in) :: x,y + b1_gauss=done/sqrt(3.0_psb_dpk_)-2*x*exp(-(x**2+y**2)) + end function b1_gauss + function b2_gauss(x,y) + implicit none + real(psb_dpk_) :: b2_gauss + real(psb_dpk_), intent(in) :: x,y + b2_gauss=done/sqrt(3.0_psb_dpk_)-2*y*exp(-(x**2+y**2)) + end function b2_gauss + function c_gauss(x,y) + implicit none + real(psb_dpk_) :: c_gauss + real(psb_dpk_), intent(in) :: x,y + c_gauss=dzero + end function c_gauss + function a1_gauss(x,y) + implicit none + real(psb_dpk_) :: a1_gauss + real(psb_dpk_), intent(in) :: x,y + a1_gauss=epsilon*exp(-(x**2+y**2)) + end function a1_gauss + function a2_gauss(x,y) + implicit none + real(psb_dpk_) :: a2_gauss + real(psb_dpk_), intent(in) :: x,y + a2_gauss=epsilon*exp(-(x**2+y**2)) + end function a2_gauss + function g_gauss(x,y) + implicit none + real(psb_dpk_) :: g_gauss + real(psb_dpk_), intent(in) :: x,y + g_gauss = dzero + if (x == done) then + g_gauss = done + else if (x == dzero) then + g_gauss = done + end if + end function g_gauss +end module amg_d_pde2d_gauss_mod diff --git a/samples/advanced/pdegen/amg_d_pde3d.F90 b/samples/advanced/pdegen/amg_d_pde3d.F90 index 10fc96e0..8c4b7b6b 100644 --- a/samples/advanced/pdegen/amg_d_pde3d.F90 +++ b/samples/advanced/pdegen/amg_d_pde3d.F90 @@ -72,6 +72,7 @@ program amg_d_pde3d use data_input use amg_d_pde3d_base_mod use amg_d_pde3d_exp_mod + use amg_d_pde3d_box_mod use amg_d_pde3d_gauss_mod use amg_d_genpde_mod #if defined(OPENMP) @@ -201,7 +202,7 @@ program amg_d_pde3d ! other variables integer(psb_ipk_) :: info, i, k character(len=20) :: name,ch_err - type(psb_d_csr_sparse_mat) :: amold + info=psb_success_ @@ -247,10 +248,13 @@ program amg_d_pde3d select case(psb_toupper(trim(pdecoeff))) case("CONST") call amg_gen_pde3d(ctxt,idim,a,b,x,desc_a,afmt,& - & a1,a2,a3,b1,b2,b3,c,g,info) + & a1_base,a2_base,a3_base,b1_base,b2_base,b3_base,c_base,g_base,info) case("EXP") call amg_gen_pde3d(ctxt,idim,a,b,x,desc_a,afmt,& & a1_exp,a2_exp,a3_exp,b1_exp,b2_exp,b3_exp,c_exp,g_exp,info) + case("BOX") + call amg_gen_pde3d(ctxt,idim,a,b,x,desc_a,afmt,& + & a1_box,a2_box,a3_box,b1_box,b2_box,b3_box,c_box,g_box,info) case("GAUSS") call amg_gen_pde3d(ctxt,idim,a,b,x,desc_a,afmt,& & a1_gauss,a2_gauss,a3_gauss,b1_gauss,b2_gauss,b3_gauss,c_gauss,g_gauss,info) @@ -412,7 +416,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) @@ -425,7 +428,7 @@ program amg_d_pde3d end if call psb_barrier(ctxt) t1 = psb_wtime() - call prec%smoothers_build(a,desc_a,info,amold=amold) + call prec%smoothers_build(a,desc_a,info) tprec = psb_wtime()-t1 if (info /= psb_success_) then call psb_errpush(psb_err_from_subroutine_,name,a_err='amg_smoothers_bld') diff --git a/samples/advanced/pdegen/amg_d_pde3d_base_mod.f90 b/samples/advanced/pdegen/amg_d_pde3d_base_mod.f90 index 0eaf0a34..4e954654 100644 --- a/samples/advanced/pdegen/amg_d_pde3d_base_mod.f90 +++ b/samples/advanced/pdegen/amg_d_pde3d_base_mod.f90 @@ -38,64 +38,64 @@ module amg_d_pde3d_base_mod use psb_base_mod, only : psb_dpk_, done, dzero real(psb_dpk_), save, private :: epsilon=done/80 contains - subroutine pde_set_parm(dat) + subroutine pde_set_parm3d_base(dat) real(psb_dpk_), intent(in) :: dat epsilon = dat - end subroutine pde_set_parm + end subroutine pde_set_parm3d_base ! ! functions parametrizing the differential equation ! - function b1(x,y,z) + function b1_base(x,y,z) implicit none - real(psb_dpk_) :: b1 + real(psb_dpk_) :: b1_base real(psb_dpk_), intent(in) :: x,y,z - b1=dzero/sqrt(3.0_psb_dpk_) - end function b1 - function b2(x,y,z) + b1_base=dzero/sqrt(3.0_psb_dpk_) + end function b1_base + function b2_base(x,y,z) implicit none - real(psb_dpk_) :: b2 + real(psb_dpk_) :: b2_base real(psb_dpk_), intent(in) :: x,y,z - b2=dzero/sqrt(3.0_psb_dpk_) - end function b2 - function b3(x,y,z) + b2_base=dzero/sqrt(3.0_psb_dpk_) + end function b2_base + function b3_base(x,y,z) implicit none - real(psb_dpk_) :: b3 + real(psb_dpk_) :: b3_base real(psb_dpk_), intent(in) :: x,y,z - b3=dzero/sqrt(3.0_psb_dpk_) - end function b3 - function c(x,y,z) + b3_base=dzero/sqrt(3.0_psb_dpk_) + end function b3_base + function c_base(x,y,z) implicit none - real(psb_dpk_) :: c + real(psb_dpk_) :: c_base real(psb_dpk_), intent(in) :: x,y,z - c=dzero - end function c - function a1(x,y,z) + c_base=dzero + end function c_base + function a1_base(x,y,z) implicit none - real(psb_dpk_) :: a1 + real(psb_dpk_) :: a1_base real(psb_dpk_), intent(in) :: x,y,z - a1=epsilon - end function a1 - function a2(x,y,z) + a1_base=epsilon + end function a1_base + function a2_base(x,y,z) implicit none - real(psb_dpk_) :: a2 + real(psb_dpk_) :: a2_base real(psb_dpk_), intent(in) :: x,y,z - a2=epsilon - end function a2 - function a3(x,y,z) + a2_base=epsilon + end function a2_base + function a3_base(x,y,z) implicit none - real(psb_dpk_) :: a3 + real(psb_dpk_) :: a3_base real(psb_dpk_), intent(in) :: x,y,z - a3=epsilon - end function a3 - function g(x,y,z) + a3_base=epsilon + end function a3_base + function g_base(x,y,z) implicit none - real(psb_dpk_) :: g + real(psb_dpk_) :: g_base real(psb_dpk_), intent(in) :: x,y,z - g = dzero + g_base = dzero if (x == done) then - g = done + g_base = done else if (x == dzero) then - g = done + g_base = done end if - end function g + end function g_base end module amg_d_pde3d_base_mod diff --git a/samples/advanced/pdegen/amg_d_pde3d_box_mod.f90 b/samples/advanced/pdegen/amg_d_pde3d_box_mod.f90 new file mode 100644 index 00000000..7062ad27 --- /dev/null +++ b/samples/advanced/pdegen/amg_d_pde3d_box_mod.f90 @@ -0,0 +1,101 @@ +! +! +! AMG4PSBLAS version 1.0 +! Algebraic Multigrid Package +! based on PSBLAS (Parallel Sparse BLAS version 3.7) +! +! (C) Copyright 2021 +! +! Salvatore Filippone +! Pasqua D'Ambra +! Fabio Durastante +! +! Redistribution and use in source and binary forms, with or without +! modification, are permitted provided that the following conditions +! are met: +! 1. Redistributions of source code must retain the above copyright +! notice, this list of conditions and the following disclaimer. +! 2. Redistributions in binary form must reproduce the above copyright +! notice, this list of conditions, and the following disclaimer in the +! documentation and/or other materials provided with the distribution. +! 3. The name of the AMG4PSBLAS group or the names of its contributors may +! not be used to endorse or promote products derived from this +! software without specific written permission. +! +! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +! ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED +! TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR +! PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE AMG4PSBLAS GROUP OR ITS CONTRIBUTORS +! BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +! CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +! SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +! INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +! CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +! ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +! POSSIBILITY OF SUCH DAMAGE. +! +module amg_d_pde3d_box_mod + use psb_base_mod, only : psb_dpk_, done, dzero + real(psb_dpk_), save, private :: epsilon=done/80 +contains + subroutine pde_set_parm3d_box(dat) + real(psb_dpk_), intent(in) :: dat + epsilon = dat + end subroutine pde_set_parm3d_box + ! + ! functions parametrizing the differential equation + ! + function b1_box(x,y,z) + implicit none + real(psb_dpk_) :: b1_box + real(psb_dpk_), intent(in) :: x,y,z + b1_box=done/sqrt(3.0_psb_dpk_) + end function b1_box + function b2_box(x,y,z) + implicit none + real(psb_dpk_) :: b2_box + real(psb_dpk_), intent(in) :: x,y,z + b2_box=done/sqrt(3.0_psb_dpk_) + end function b2_box + function b3_box(x,y,z) + implicit none + real(psb_dpk_) :: b3_box + real(psb_dpk_), intent(in) :: x,y,z + b3_box=done/sqrt(3.0_psb_dpk_) + end function b3_box + function c_box(x,y,z) + implicit none + real(psb_dpk_) :: c_box + real(psb_dpk_), intent(in) :: x,y,z + c_box=dzero + end function c_box + function a1_box(x,y,z) + implicit none + real(psb_dpk_) :: a1_box + real(psb_dpk_), intent(in) :: x,y,z + a1_box=epsilon + end function a1_box + function a2_box(x,y,z) + implicit none + real(psb_dpk_) :: a2_box + real(psb_dpk_), intent(in) :: x,y,z + a2_box=epsilon + end function a2_box + function a3_box(x,y,z) + implicit none + real(psb_dpk_) :: a3_box + real(psb_dpk_), intent(in) :: x,y,z + a3_box=epsilon + end function a3_box + function g_box(x,y,z) + implicit none + real(psb_dpk_) :: g_box + real(psb_dpk_), intent(in) :: x,y,z + g_box= dzero + if (x == done) then + g_box = done + else if (x == dzero) then + g_box = done + end if + end function g_box +end module amg_d_pde3d_box_mod diff --git a/samples/advanced/pdegen/amg_d_pde3d_exp_mod.f90 b/samples/advanced/pdegen/amg_d_pde3d_exp_mod.f90 index e7bcf6ef..fea3b8a4 100644 --- a/samples/advanced/pdegen/amg_d_pde3d_exp_mod.f90 +++ b/samples/advanced/pdegen/amg_d_pde3d_exp_mod.f90 @@ -38,10 +38,10 @@ module amg_d_pde3d_exp_mod use psb_base_mod, only : psb_dpk_, done, dzero real(psb_dpk_), save, private :: epsilon=done/160 contains - subroutine pde_set_parm(dat) + subroutine pde_set_parm3d_exp(dat) real(psb_dpk_), intent(in) :: dat epsilon = dat - end subroutine pde_set_parm + end subroutine pde_set_parm3d_exp ! ! functions parametrizing the differential equation ! diff --git a/samples/advanced/pdegen/amg_d_pde3d_gauss_mod.f90 b/samples/advanced/pdegen/amg_d_pde3d_gauss_mod.f90 index 8dd5f71a..c2403131 100644 --- a/samples/advanced/pdegen/amg_d_pde3d_gauss_mod.f90 +++ b/samples/advanced/pdegen/amg_d_pde3d_gauss_mod.f90 @@ -38,10 +38,10 @@ module amg_d_pde3d_gauss_mod use psb_base_mod, only : psb_dpk_, done, dzero real(psb_dpk_), save, private :: epsilon=done/80 contains - subroutine pde_set_parm(dat) + subroutine pde_set_parm3d_gauss(dat) real(psb_dpk_), intent(in) :: dat epsilon = dat - end subroutine pde_set_parm + end subroutine pde_set_parm3d_gauss ! ! functions parametrizing the differential equation ! diff --git a/samples/advanced/pdegen/amg_s_pde2d.F90 b/samples/advanced/pdegen/amg_s_pde2d.F90 index f3842a65..bcc995ea 100644 --- a/samples/advanced/pdegen/amg_s_pde2d.F90 +++ b/samples/advanced/pdegen/amg_s_pde2d.F90 @@ -72,6 +72,7 @@ program amg_s_pde2d use amg_s_pde2d_base_mod use amg_s_pde2d_exp_mod use amg_s_pde2d_box_mod + use amg_s_pde2d_gauss_mod use amg_s_genpde_mod #if defined(OPENMP) use omp_lib @@ -124,60 +125,64 @@ program amg_s_pde2d integer(psb_ipk_) :: outer_sweeps ! number of outer sweeps: sweeps for 1-level, ! AMG cycles for ML ! general AMG data - character(len=16) :: mlcycle ! AMG cycle type - integer(psb_ipk_) :: maxlevs ! maximum number of levels in AMG preconditioner + character(len=32) :: mlcycle ! AMG cycle type + integer(psb_ipk_) :: maxlevs ! maximum number of levels in AMG preconditioner ! AMG aggregation - character(len=16) :: aggr_prol ! aggregation type: SMOOTHED, NONSMOOTHED - character(len=16) :: par_aggr_alg ! parallel aggregation algorithm: DEC, SYMDEC - character(len=16) :: aggr_type ! Type of aggregation SOC1, SOC2, MATCHBOXP - integer(psb_ipk_) :: aggr_size ! Requested size of the aggregates for MATCHBOXP - character(len=16) :: aggr_ord ! ordering for aggregation: NATURAL, DEGREE - character(len=16) :: aggr_filter ! filtering: FILTER, NO_FILTER - real(psb_spk_) :: mncrratio ! minimum aggregation ratio + character(len=32) :: aggr_prol ! aggregation type: SMOOTHED, NONSMOOTHED + character(len=32) :: par_aggr_alg ! parallel aggregation algorithm: DEC, SYMDEC + character(len=32) :: aggr_type ! Type of aggregation SOC1, SOC2, MATCHBOXP + integer(psb_ipk_) :: aggr_size ! Requested size of the aggregates for MATCHBOXP + character(len=32) :: aggr_ord ! ordering for aggregation: NATURAL, DEGREE + character(len=32) :: aggr_filter ! filtering: FILTER, NO_FILTER + real(psb_spk_) :: mncrratio ! minimum aggregation ratio real(psb_spk_), allocatable :: athresv(:) ! smoothed aggregation threshold vector - integer(psb_ipk_) :: thrvsz ! size of threshold vector - real(psb_spk_) :: athres ! smoothed aggregation threshold - integer(psb_ipk_) :: csizepp ! minimum size of coarsest matrix per process + integer(psb_ipk_) :: thrvsz ! size of threshold vector + real(psb_spk_) :: athres ! smoothed aggregation threshold + integer(psb_ipk_) :: csizepp ! minimum size of coarsest matrix per process ! AMG smoother or pre-smoother; also 1-lev preconditioner - character(len=16) :: smther ! (pre-)smoother type: BJAC, AS + character(len=32) :: 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=32) :: pvariant ! polynomial 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 - character(len=16) :: solve ! local subsolver type: ILU, MILU, ILUT, + character(len=32) :: restr ! restriction over application of AS + character(len=32) :: prol ! prolongation over application of AS + character(len=32) :: solve ! local subsolver type: ILU, MILU, ILUT, ! UMF, MUMPS, SLU, FWGS, BWGS, JAC integer(psb_ipk_) :: ssweeps ! inner solver sweeps - character(len=16) :: variant ! AINV variant: LLK, etc + character(len=32) :: 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_spk_) :: thr ! threshold for ILUT factorization ! AMG post-smoother; ignored by 1-lev preconditioner - character(len=16) :: smther2 ! post-smoother type: BJAC, AS + character(len=32) :: smther2 ! post-smoother type: BJAC, AS integer(psb_ipk_) :: jsweeps2 ! post-smoother sweeps + integer(psb_ipk_) :: degree2 ! degree for polynomial smoother + character(len=32) :: pvariant2 ! polynomial 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 - character(len=16) :: solve2 ! local subsolver type: ILU, MILU, ILUT, + character(len=32) :: restr2 ! restriction over application of AS + character(len=32) :: prol2 ! prolongation over application of AS + character(len=32) :: solve2 ! local subsolver type: ILU, MILU, ILUT, ! UMF, MUMPS, SLU, FWGS, BWGS, JAC integer(psb_ipk_) :: ssweeps2 ! inner solver sweeps - character(len=16) :: variant2 ! AINV variant: LLK, etc + character(len=32) :: 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_spk_) :: thr2 ! threshold for ILUT factorization ! coarsest-level solver - character(len=16) :: cmat ! coarsest matrix layout: REPL, DIST - character(len=16) :: csolve ! coarsest-lev solver: BJAC, SLUDIST (distr. - ! mat.); UMF, MUMPS, SLU, ILU, ILUT, MILU - ! (repl. mat.) - character(len=16) :: csbsolve ! coarsest-lev local subsolver: ILU, ILUT, - ! MILU, UMF, MUMPS, SLU - integer(psb_ipk_) :: cfill ! fill-in for incomplete LU factorization - real(psb_spk_) :: cthres ! threshold for ILUT factorization - integer(psb_ipk_) :: cjswp ! sweeps for GS or JAC coarsest-lev subsolver + character(len=32) :: cmat ! coarsest matrix layout: REPL, DIST + character(len=32) :: csolve ! coarsest-lev solver: BJAC, SLUDIST (distr. + ! mat.); UMF, MUMPS, SLU, ILU, ILUT, MILU + ! (repl. mat.) + character(len=32) :: csbsolve ! coarsest-lev local subsolver: ILU, ILUT, + ! MILU, UMF, MUMPS, SLU + integer(psb_ipk_) :: cfill ! fill-in for incomplete LU factorization + real(psb_spk_) :: cthres ! threshold for ILUT factorization + integer(psb_ipk_) :: cjswp ! sweeps for GS or JAC coarsest-lev subsolver ! Dump data logical :: dump = .false. @@ -241,13 +246,16 @@ program amg_s_pde2d select case(psb_toupper(trim(pdecoeff))) case("CONST") call amg_gen_pde2d(ctxt,idim,a,b,x,desc_a,afmt,& - & a1,a2,b1,b2,c,g,info) + & a1_base,a2_base,b1_base,b2_base,c_base,g_base,info) case("EXP") call amg_gen_pde2d(ctxt,idim,a,b,x,desc_a,afmt,& & a1_exp,a2_exp,b1_exp,b2_exp,c_exp,g_exp,info) case("BOX") call amg_gen_pde2d(ctxt,idim,a,b,x,desc_a,afmt,& & a1_box,a2_box,b1_box,b2_box,c_box,g_box,info) + case("GAUSS") + call amg_gen_pde2d(ctxt,idim,a,b,x,desc_a,afmt,& + & a1_gauss,a2_gauss,b1_gauss,b2_gauss,c_gauss,g_gauss,info) case default info=psb_err_from_subroutine_ ch_err='amg_gen_pdecoeff' @@ -281,10 +289,12 @@ program amg_s_pde2d ! 1-level sweeps from "outer_sweeps" call prec%set('smoother_sweeps', p_choice%jsweeps, info) - case ('BJAC') + case ('BJAC','POLY') 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('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) call prec%set('sub_fillin', p_choice%fill, info) @@ -308,7 +318,7 @@ program amg_s_pde2d call prec%set('ml_cycle', p_choice%mlcycle, info) call prec%set('outer_sweeps', p_choice%outer_sweeps,info) if (p_choice%csizepp>0)& - & call prec%set('min_coarse_size_per_process', p_choice%csizepp, info) + & call prec%set('min_coarse_size_per_process', p_choice%csizepp, info) if (p_choice%mncrratio>1)& & call prec%set('min_cr_ratio', p_choice%mncrratio, info) if (p_choice%maxlevs>0)& @@ -332,7 +342,9 @@ program amg_s_pde2d call prec%set('smoother_type', p_choice%smther, info) call prec%set('smoother_sweeps', p_choice%jsweeps, info) - + call prec%set('poly_degree', p_choice%degree, info) + call prec%set('poly_variant', p_choice%pvariant, info) + select case (psb_toupper(p_choice%smther)) case ('GS','BWGS','FBGS','JACOBI','L1-JACOBI','L1-FBGS') ! do nothing @@ -362,6 +374,8 @@ program amg_s_pde2d 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('poly_degree', p_choice%degree2, info,pos='post') + call prec%set('poly_variant', p_choice%pvariant2, info,pos='post') select case (psb_toupper(p_choice%smther2)) case ('GS','BWGS','FBGS','JACOBI','L1-JACOBI','L1-FBGS') ! do nothing @@ -577,6 +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%pvariant,inp_unit) ! 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 @@ -589,11 +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%pvariant2,inp_unit) ! 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 @@ -659,6 +677,8 @@ contains ! broadcast first (pre-)smoother / 1-lev prec data call psb_bcast(ctxt,prec%smther) call psb_bcast(ctxt,prec%jsweeps) + call psb_bcast(ctxt,prec%degree) + call psb_bcast(ctxt,prec%pvariant) call psb_bcast(ctxt,prec%novr) call psb_bcast(ctxt,prec%restr) call psb_bcast(ctxt,prec%prol) @@ -671,6 +691,8 @@ contains ! broadcast second (post-)smoother call psb_bcast(ctxt,prec%smther2) call psb_bcast(ctxt,prec%jsweeps2) + call psb_bcast(ctxt,prec%degree2) + call psb_bcast(ctxt,prec%pvariant2) call psb_bcast(ctxt,prec%novr2) call psb_bcast(ctxt,prec%restr2) call psb_bcast(ctxt,prec%prol2) diff --git a/samples/advanced/pdegen/amg_s_pde2d_base_mod.f90 b/samples/advanced/pdegen/amg_s_pde2d_base_mod.f90 index 462c6154..d5cbc6d0 100644 --- a/samples/advanced/pdegen/amg_s_pde2d_base_mod.f90 +++ b/samples/advanced/pdegen/amg_s_pde2d_base_mod.f90 @@ -38,52 +38,52 @@ module amg_s_pde2d_base_mod use psb_base_mod, only : psb_spk_, szero, sone real(psb_spk_), save, private :: epsilon=sone/80 contains - subroutine pde_set_parm(dat) + subroutine pde_set_parm2d_base(dat) real(psb_spk_), intent(in) :: dat epsilon = dat - end subroutine pde_set_parm + end subroutine pde_set_parm2d_base ! ! functions parametrizing the differential equation ! - function b1(x,y) + function b1_base(x,y) implicit none - real(psb_spk_) :: b1 + real(psb_spk_) :: b1_base real(psb_spk_), intent(in) :: x,y - b1 = szero/1.414_psb_spk_ - end function b1 - function b2(x,y) + b1_base = szero/1.414_psb_spk_ + end function b1_base + function b2_base(x,y) implicit none - real(psb_spk_) :: b2 + real(psb_spk_) :: b2_base real(psb_spk_), intent(in) :: x,y - b2 = szero/1.414_psb_spk_ - end function b2 - function c(x,y) + b2_base = szero/1.414_psb_spk_ + end function b2_base + function c_base(x,y) implicit none - real(psb_spk_) :: c + real(psb_spk_) :: c_base real(psb_spk_), intent(in) :: x,y - c = szero - end function c - function a1(x,y) + c_base = szero + end function c_base + function a1_base(x,y) implicit none - real(psb_spk_) :: a1 + real(psb_spk_) :: a1_base real(psb_spk_), intent(in) :: x,y - a1=sone*epsilon - end function a1 - function a2(x,y) + a1_base=sone*epsilon + end function a1_base + function a2_base(x,y) implicit none - real(psb_spk_) :: a2 + real(psb_spk_) :: a2_base real(psb_spk_), intent(in) :: x,y - a2=sone*epsilon - end function a2 - function g(x,y) + a2_base=sone*epsilon + end function a2_base + function g_base(x,y) implicit none - real(psb_spk_) :: g + real(psb_spk_) :: g_base real(psb_spk_), intent(in) :: x,y - g = szero + g_base = szero if (x == sone) then - g = sone + g_base = sone else if (x == szero) then - g = sone + g_base = sone end if - end function g + end function g_base end module amg_s_pde2d_base_mod diff --git a/samples/advanced/pdegen/amg_s_pde2d_box_mod.f90 b/samples/advanced/pdegen/amg_s_pde2d_box_mod.f90 index 9183521b..c36f8eb0 100644 --- a/samples/advanced/pdegen/amg_s_pde2d_box_mod.f90 +++ b/samples/advanced/pdegen/amg_s_pde2d_box_mod.f90 @@ -38,10 +38,10 @@ module amg_s_pde2d_box_mod use psb_base_mod, only : psb_spk_, szero, sone real(psb_spk_), save, private :: epsilon=sone/80 contains - subroutine pde_set_parm(dat) + subroutine pde_set_parm2d_box(dat) real(psb_spk_), intent(in) :: dat epsilon = dat - end subroutine pde_set_parm + end subroutine pde_set_parm2d_box ! ! functions parametrizing the differential equation ! diff --git a/samples/advanced/pdegen/amg_s_pde2d_exp_mod.f90 b/samples/advanced/pdegen/amg_s_pde2d_exp_mod.f90 index 3657546d..acc786b2 100644 --- a/samples/advanced/pdegen/amg_s_pde2d_exp_mod.f90 +++ b/samples/advanced/pdegen/amg_s_pde2d_exp_mod.f90 @@ -38,10 +38,10 @@ module amg_s_pde2d_exp_mod use psb_base_mod, only : psb_spk_, sone, szero real(psb_spk_), save, private :: epsilon=sone/80 contains - subroutine pde_set_parm(dat) + subroutine pde_set_parm2d_exp(dat) real(psb_spk_), intent(in) :: dat epsilon = dat - end subroutine pde_set_parm + end subroutine pde_set_parm2d_exp ! ! functions parametrizing the differential equation ! diff --git a/samples/advanced/pdegen/amg_s_pde2d_gauss_mod.f90 b/samples/advanced/pdegen/amg_s_pde2d_gauss_mod.f90 new file mode 100644 index 00000000..dd234d66 --- /dev/null +++ b/samples/advanced/pdegen/amg_s_pde2d_gauss_mod.f90 @@ -0,0 +1,89 @@ +! +! +! AMG4PSBLAS version 1.0 +! Algebraic Multigrid Package +! based on PSBLAS (Parallel Sparse BLAS version 3.7) +! +! (C) Copyright 2021 +! +! Salvatore Filippone +! Pasqua D'Ambra +! Fabio Durastante +! +! Redistribution and use in source and binary forms, with or without +! modification, are permitted provided that the following conditions +! are met: +! 1. Redistributions of source code must retain the above copyright +! notice, this list of conditions and the following disclaimer. +! 2. Redistributions in binary form must reproduce the above copyright +! notice, this list of conditions, and the following disclaimer in the +! documentation and/or other materials provided with the distribution. +! 3. The name of the AMG4PSBLAS group or the names of its contributors may +! not be used to endorse or promote products derived from this +! software without specific written permission. +! +! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +! ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED +! TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR +! PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE AMG4PSBLAS GROUP OR ITS CONTRIBUTORS +! BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +! CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +! SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +! INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +! CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +! ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +! POSSIBILITY OF SUCH DAMAGE. +! +module amg_s_pde2d_gauss_mod + use psb_base_mod, only : psb_spk_, sone, szero + real(psb_spk_), save, private :: epsilon=sone/80 +contains + subroutine pde_set_parm2d_gauss(dat) + real(psb_spk_), intent(in) :: dat + epsilon = dat + end subroutine pde_set_parm2d_gauss + ! + ! functions parametrizing the differential equation + ! + function b1_gauss(x,y) + implicit none + real(psb_spk_) :: b1_gauss + real(psb_spk_), intent(in) :: x,y + b1_gauss=sone/sqrt(3.0_psb_spk_)-2*x*exp(-(x**2+y**2)) + end function b1_gauss + function b2_gauss(x,y) + implicit none + real(psb_spk_) :: b2_gauss + real(psb_spk_), intent(in) :: x,y + b2_gauss=sone/sqrt(3.0_psb_spk_)-2*y*exp(-(x**2+y**2)) + end function b2_gauss + function c_gauss(x,y) + implicit none + real(psb_spk_) :: c_gauss + real(psb_spk_), intent(in) :: x,y + c_gauss=szero + end function c_gauss + function a1_gauss(x,y) + implicit none + real(psb_spk_) :: a1_gauss + real(psb_spk_), intent(in) :: x,y + a1_gauss=epsilon*exp(-(x**2+y**2)) + end function a1_gauss + function a2_gauss(x,y) + implicit none + real(psb_spk_) :: a2_gauss + real(psb_spk_), intent(in) :: x,y + a2_gauss=epsilon*exp(-(x**2+y**2)) + end function a2_gauss + function g_gauss(x,y) + implicit none + real(psb_spk_) :: g_gauss + real(psb_spk_), intent(in) :: x,y + g_gauss = szero + if (x == sone) then + g_gauss = sone + else if (x == szero) then + g_gauss = sone + end if + end function g_gauss +end module amg_s_pde2d_gauss_mod diff --git a/samples/advanced/pdegen/amg_s_pde3d.F90 b/samples/advanced/pdegen/amg_s_pde3d.F90 index 97abad28..fe53cd8b 100644 --- a/samples/advanced/pdegen/amg_s_pde3d.F90 +++ b/samples/advanced/pdegen/amg_s_pde3d.F90 @@ -72,6 +72,7 @@ program amg_s_pde3d use data_input use amg_s_pde3d_base_mod use amg_s_pde3d_exp_mod + use amg_s_pde3d_box_mod use amg_s_pde3d_gauss_mod use amg_s_genpde_mod #if defined(OPENMP) @@ -125,16 +126,16 @@ program amg_s_pde3d integer(psb_ipk_) :: outer_sweeps ! number of outer sweeps: sweeps for 1-level, ! AMG cycles for ML ! general AMG data - character(len=16) :: mlcycle ! AMG cycle type + character(len=32) :: mlcycle ! AMG cycle type integer(psb_ipk_) :: maxlevs ! maximum number of levels in AMG preconditioner ! AMG aggregation - character(len=16) :: aggr_prol ! aggregation type: SMOOTHED, NONSMOOTHED - character(len=16) :: par_aggr_alg ! parallel aggregation algorithm: DEC, SYMDEC - character(len=16) :: aggr_type ! Type of aggregation SOC1, SOC2, MATCHBOXP + character(len=32) :: aggr_prol ! aggregation type: SMOOTHED, NONSMOOTHED + character(len=32) :: par_aggr_alg ! parallel aggregation algorithm: DEC, SYMDEC + character(len=32) :: aggr_type ! Type of aggregation SOC1, SOC2, MATCHBOXP integer(psb_ipk_) :: aggr_size ! Requested size of the aggregates for MATCHBOXP - character(len=16) :: aggr_ord ! ordering for aggregation: NATURAL, DEGREE - character(len=16) :: aggr_filter ! filtering: FILTER, NO_FILTER + character(len=32) :: aggr_ord ! ordering for aggregation: NATURAL, DEGREE + character(len=32) :: aggr_filter ! filtering: FILTER, NO_FILTER real(psb_spk_) :: mncrratio ! minimum aggregation ratio real(psb_spk_), allocatable :: athresv(:) ! smoothed aggregation threshold vector integer(psb_ipk_) :: thrvsz ! size of threshold vector @@ -142,39 +143,43 @@ program amg_s_pde3d integer(psb_ipk_) :: csizepp ! minimum size of coarsest matrix per process ! AMG smoother or pre-smoother; also 1-lev preconditioner - character(len=16) :: smther ! (pre-)smoother type: BJAC, AS + character(len=32) :: 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=32) :: pvariant ! polynomial 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 - character(len=16) :: solve ! local subsolver type: ILU, MILU, ILUT, + character(len=32) :: restr ! restriction over application of AS + character(len=32) :: prol ! prolongation over application of AS + character(len=32) :: solve ! local subsolver type: ILU, MILU, ILUT, ! UMF, MUMPS, SLU, FWGS, BWGS, JAC integer(psb_ipk_) :: ssweeps ! inner solver sweeps - character(len=16) :: variant ! AINV variant: LLK, etc + character(len=32) :: 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_spk_) :: thr ! threshold for ILUT factorization ! AMG post-smoother; ignored by 1-lev preconditioner - character(len=16) :: smther2 ! post-smoother type: BJAC, AS + character(len=32) :: smther2 ! post-smoother type: BJAC, AS integer(psb_ipk_) :: jsweeps2 ! post-smoother sweeps + integer(psb_ipk_) :: degree2 ! degree for polynomial smoother + character(len=32) :: pvariant2 ! polynomial 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 - character(len=16) :: solve2 ! local subsolver type: ILU, MILU, ILUT, + character(len=32) :: restr2 ! restriction over application of AS + character(len=32) :: prol2 ! prolongation over application of AS + character(len=32) :: solve2 ! local subsolver type: ILU, MILU, ILUT, ! UMF, MUMPS, SLU, FWGS, BWGS, JAC integer(psb_ipk_) :: ssweeps2 ! inner solver sweeps - character(len=16) :: variant2 ! AINV variant: LLK, etc + character(len=32) :: 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_spk_) :: thr2 ! threshold for ILUT factorization ! coarsest-level solver - character(len=16) :: cmat ! coarsest matrix layout: REPL, DIST - character(len=16) :: csolve ! coarsest-lev solver: BJAC, SLUDIST (distr. + character(len=32) :: cmat ! coarsest matrix layout: REPL, DIST + character(len=32) :: csolve ! coarsest-lev solver: BJAC, SLUDIST (distr. ! mat.); UMF, MUMPS, SLU, ILU, ILUT, MILU ! (repl. mat.) - character(len=16) :: csbsolve ! coarsest-lev local subsolver: ILU, ILUT, + character(len=32) :: csbsolve ! coarsest-lev local subsolver: ILU, ILUT, ! MILU, UMF, MUMPS, SLU integer(psb_ipk_) :: cfill ! fill-in for incomplete LU factorization real(psb_spk_) :: cthres ! threshold for ILUT factorization @@ -197,7 +202,7 @@ program amg_s_pde3d ! other variables integer(psb_ipk_) :: info, i, k character(len=20) :: name,ch_err - type(psb_s_csr_sparse_mat) :: amold + info=psb_success_ @@ -243,10 +248,13 @@ program amg_s_pde3d select case(psb_toupper(trim(pdecoeff))) case("CONST") call amg_gen_pde3d(ctxt,idim,a,b,x,desc_a,afmt,& - & a1,a2,a3,b1,b2,b3,c,g,info) + & a1_base,a2_base,a3_base,b1_base,b2_base,b3_base,c_base,g_base,info) case("EXP") call amg_gen_pde3d(ctxt,idim,a,b,x,desc_a,afmt,& & a1_exp,a2_exp,a3_exp,b1_exp,b2_exp,b3_exp,c_exp,g_exp,info) + case("BOX") + call amg_gen_pde3d(ctxt,idim,a,b,x,desc_a,afmt,& + & a1_box,a2_box,a3_box,b1_box,b2_box,b3_box,c_box,g_box,info) case("GAUSS") call amg_gen_pde3d(ctxt,idim,a,b,x,desc_a,afmt,& & a1_gauss,a2_gauss,a3_gauss,b1_gauss,b2_gauss,b3_gauss,c_gauss,g_gauss,info) @@ -285,10 +293,12 @@ program amg_s_pde3d ! 1-level sweeps from "outer_sweeps" call prec%set('smoother_sweeps', p_choice%jsweeps, info) - case ('BJAC') + case ('BJAC','POLY') 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('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) call prec%set('sub_fillin', p_choice%fill, info) @@ -336,7 +346,9 @@ program amg_s_pde3d call prec%set('smoother_type', p_choice%smther, info) call prec%set('smoother_sweeps', p_choice%jsweeps, info) - + call prec%set('poly_degree', p_choice%degree, info) + call prec%set('poly_variant', p_choice%pvariant, info) + select case (psb_toupper(p_choice%smther)) case ('GS','BWGS','FBGS','JACOBI','L1-JACOBI','L1-FBGS') ! do nothing @@ -366,6 +378,8 @@ program amg_s_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('poly_degree', p_choice%degree2, info,pos='post') + call prec%set('poly_variant', p_choice%pvariant2, info,pos='post') select case (psb_toupper(p_choice%smther2)) case ('GS','BWGS','FBGS','JACOBI','L1-JACOBI','L1-FBGS') ! do nothing @@ -414,7 +428,7 @@ program amg_s_pde3d end if call psb_barrier(ctxt) t1 = psb_wtime() - call prec%smoothers_build(a,desc_a,info,amold=amold) + call prec%smoothers_build(a,desc_a,info) tprec = psb_wtime()-t1 if (info /= psb_success_) then call psb_errpush(psb_err_from_subroutine_,name,a_err='amg_smoothers_bld') @@ -581,6 +595,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%pvariant,inp_unit) ! 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 @@ -593,11 +609,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%pvariant2,inp_unit) ! 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 @@ -663,6 +681,8 @@ contains ! broadcast first (pre-)smoother / 1-lev prec data call psb_bcast(ctxt,prec%smther) call psb_bcast(ctxt,prec%jsweeps) + call psb_bcast(ctxt,prec%degree) + call psb_bcast(ctxt,prec%pvariant) call psb_bcast(ctxt,prec%novr) call psb_bcast(ctxt,prec%restr) call psb_bcast(ctxt,prec%prol) @@ -675,6 +695,8 @@ contains ! broadcast second (post-)smoother call psb_bcast(ctxt,prec%smther2) call psb_bcast(ctxt,prec%jsweeps2) + call psb_bcast(ctxt,prec%degree2) + call psb_bcast(ctxt,prec%pvariant2) call psb_bcast(ctxt,prec%novr2) call psb_bcast(ctxt,prec%restr2) call psb_bcast(ctxt,prec%prol2) diff --git a/samples/advanced/pdegen/amg_s_pde3d_base_mod.f90 b/samples/advanced/pdegen/amg_s_pde3d_base_mod.f90 index 0ce83989..3dbd039f 100644 --- a/samples/advanced/pdegen/amg_s_pde3d_base_mod.f90 +++ b/samples/advanced/pdegen/amg_s_pde3d_base_mod.f90 @@ -38,64 +38,64 @@ module amg_s_pde3d_base_mod use psb_base_mod, only : psb_spk_, sone, szero real(psb_spk_), save, private :: epsilon=sone/80 contains - subroutine pde_set_parm(dat) + subroutine pde_set_parm3d_base(dat) real(psb_spk_), intent(in) :: dat epsilon = dat - end subroutine pde_set_parm + end subroutine pde_set_parm3d_base ! ! functions parametrizing the differential equation ! - function b1(x,y,z) + function b1_base(x,y,z) implicit none - real(psb_spk_) :: b1 + real(psb_spk_) :: b1_base real(psb_spk_), intent(in) :: x,y,z - b1=sone/sqrt(3.0_psb_spk_) - end function b1 - function b2(x,y,z) + b1_base=szero/sqrt(3.0_psb_spk_) + end function b1_base + function b2_base(x,y,z) implicit none - real(psb_spk_) :: b2 + real(psb_spk_) :: b2_base real(psb_spk_), intent(in) :: x,y,z - b2=sone/sqrt(3.0_psb_spk_) - end function b2 - function b3(x,y,z) + b2_base=szero/sqrt(3.0_psb_spk_) + end function b2_base + function b3_base(x,y,z) implicit none - real(psb_spk_) :: b3 + real(psb_spk_) :: b3_base real(psb_spk_), intent(in) :: x,y,z - b3=sone/sqrt(3.0_psb_spk_) - end function b3 - function c(x,y,z) + b3_base=szero/sqrt(3.0_psb_spk_) + end function b3_base + function c_base(x,y,z) implicit none - real(psb_spk_) :: c + real(psb_spk_) :: c_base real(psb_spk_), intent(in) :: x,y,z - c=szero - end function c - function a1(x,y,z) + c_base=szero + end function c_base + function a1_base(x,y,z) implicit none - real(psb_spk_) :: a1 + real(psb_spk_) :: a1_base real(psb_spk_), intent(in) :: x,y,z - a1=epsilon - end function a1 - function a2(x,y,z) + a1_base=epsilon + end function a1_base + function a2_base(x,y,z) implicit none - real(psb_spk_) :: a2 + real(psb_spk_) :: a2_base real(psb_spk_), intent(in) :: x,y,z - a2=epsilon - end function a2 - function a3(x,y,z) + a2_base=epsilon + end function a2_base + function a3_base(x,y,z) implicit none - real(psb_spk_) :: a3 + real(psb_spk_) :: a3_base real(psb_spk_), intent(in) :: x,y,z - a3=epsilon - end function a3 - function g(x,y,z) + a3_base=epsilon + end function a3_base + function g_base(x,y,z) implicit none - real(psb_spk_) :: g + real(psb_spk_) :: g_base real(psb_spk_), intent(in) :: x,y,z - g = szero + g_base = szero if (x == sone) then - g = sone + g_base = sone else if (x == szero) then - g = sone + g_base = sone end if - end function g + end function g_base end module amg_s_pde3d_base_mod diff --git a/samples/advanced/pdegen/amg_s_pde3d_box_mod.f90 b/samples/advanced/pdegen/amg_s_pde3d_box_mod.f90 new file mode 100644 index 00000000..e0a1a5e3 --- /dev/null +++ b/samples/advanced/pdegen/amg_s_pde3d_box_mod.f90 @@ -0,0 +1,101 @@ +! +! +! AMG4PSBLAS version 1.0 +! Algebraic Multigrid Package +! based on PSBLAS (Parallel Sparse BLAS version 3.7) +! +! (C) Copyright 2021 +! +! Salvatore Filippone +! Pasqua D'Ambra +! Fabio Durastante +! +! Redistribution and use in source and binary forms, with or without +! modification, are permitted provided that the following conditions +! are met: +! 1. Redistributions of source code must retain the above copyright +! notice, this list of conditions and the following disclaimer. +! 2. Redistributions in binary form must reproduce the above copyright +! notice, this list of conditions, and the following disclaimer in the +! documentation and/or other materials provided with the distribution. +! 3. The name of the AMG4PSBLAS group or the names of its contributors may +! not be used to endorse or promote products derived from this +! software without specific written permission. +! +! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +! ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED +! TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR +! PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE AMG4PSBLAS GROUP OR ITS CONTRIBUTORS +! BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +! CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +! SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +! INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +! CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +! ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +! POSSIBILITY OF SUCH DAMAGE. +! +module amg_s_pde3d_box_mod + use psb_base_mod, only : psb_spk_, sone, szero + real(psb_spk_), save, private :: epsilon=sone/80 +contains + subroutine pde_set_parm3d_box(dat) + real(psb_spk_), intent(in) :: dat + epsilon = dat + end subroutine pde_set_parm3d_box + ! + ! functions parametrizing the differential equation + ! + function b1_box(x,y,z) + implicit none + real(psb_spk_) :: b1_box + real(psb_spk_), intent(in) :: x,y,z + b1_box=sone/sqrt(3.0_psb_spk_) + end function b1_box + function b2_box(x,y,z) + implicit none + real(psb_spk_) :: b2_box + real(psb_spk_), intent(in) :: x,y,z + b2_box=sone/sqrt(3.0_psb_spk_) + end function b2_box + function b3_box(x,y,z) + implicit none + real(psb_spk_) :: b3_box + real(psb_spk_), intent(in) :: x,y,z + b3_box=sone/sqrt(3.0_psb_spk_) + end function b3_box + function c_box(x,y,z) + implicit none + real(psb_spk_) :: c_box + real(psb_spk_), intent(in) :: x,y,z + c_box=szero + end function c_box + function a1_box(x,y,z) + implicit none + real(psb_spk_) :: a1_box + real(psb_spk_), intent(in) :: x,y,z + a1_box=epsilon + end function a1_box + function a2_box(x,y,z) + implicit none + real(psb_spk_) :: a2_box + real(psb_spk_), intent(in) :: x,y,z + a2_box=epsilon + end function a2_box + function a3_box(x,y,z) + implicit none + real(psb_spk_) :: a3_box + real(psb_spk_), intent(in) :: x,y,z + a3_box=epsilon + end function a3_box + function g_box(x,y,z) + implicit none + real(psb_spk_) :: g_box + real(psb_spk_), intent(in) :: x,y,z + g_box= szero + if (x == sone) then + g_box = sone + else if (x == szero) then + g_box = sone + end if + end function g_box +end module amg_s_pde3d_box_mod diff --git a/samples/advanced/pdegen/amg_s_pde3d_exp_mod.f90 b/samples/advanced/pdegen/amg_s_pde3d_exp_mod.f90 index 8ec96d00..d1f8fbea 100644 --- a/samples/advanced/pdegen/amg_s_pde3d_exp_mod.f90 +++ b/samples/advanced/pdegen/amg_s_pde3d_exp_mod.f90 @@ -38,10 +38,10 @@ module amg_s_pde3d_exp_mod use psb_base_mod, only : psb_spk_, sone, szero real(psb_spk_), save, private :: epsilon=sone/160 contains - subroutine pde_set_parm(dat) + subroutine pde_set_parm3d_exp(dat) real(psb_spk_), intent(in) :: dat epsilon = dat - end subroutine pde_set_parm + end subroutine pde_set_parm3d_exp ! ! functions parametrizing the differential equation ! diff --git a/samples/advanced/pdegen/amg_s_pde3d_gauss_mod.f90 b/samples/advanced/pdegen/amg_s_pde3d_gauss_mod.f90 index fa6362e0..d7bb81ab 100644 --- a/samples/advanced/pdegen/amg_s_pde3d_gauss_mod.f90 +++ b/samples/advanced/pdegen/amg_s_pde3d_gauss_mod.f90 @@ -38,10 +38,10 @@ module amg_s_pde3d_gauss_mod use psb_base_mod, only : psb_spk_, sone, szero real(psb_spk_), save, private :: epsilon=sone/80 contains - subroutine pde_set_parm(dat) + subroutine pde_set_parm3d_gauss(dat) real(psb_spk_), intent(in) :: dat epsilon = dat - end subroutine pde_set_parm + end subroutine pde_set_parm3d_gauss ! ! functions parametrizing the differential equation ! diff --git a/samples/advanced/pdegen/runs/amg_pde2d.inp b/samples/advanced/pdegen/runs/amg_pde2d.inp index 9f6a4c42..95f5f436 100644 --- a/samples/advanced/pdegen/runs/amg_pde2d.inp +++ b/samples/advanced/pdegen/runs/amg_pde2d.inp @@ -1,7 +1,7 @@ %%%%%%%%%%% General arguments % Lines starting with % are ignored. CSR ! Storage format CSR COO JAD 0200 ! IDIM; domain size. Linear system size is IDIM**2 -CONST ! PDECOEFF: CONST, EXP, BOX Coefficients of the PDE +CONST ! PDECOEFF: CONST, EXP, BOX, GAUSS Coefficients of the PDE CG ! Iterative method: BiCGSTAB BiCGSTABL BiCG CG CGS FCG GCR RGMRES 2 ! ISTOPC 00500 ! ITMAX @@ -9,23 +9,33 @@ CG ! Iterative method: BiCGSTAB BiCGSTABL BiCG CG CGS F 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 +ML-VBM-VCYCLE-FBGS-D-BJAC ! Longer descriptive name for preconditioner (up to 20 chars) +ML ! Preconditioner type: NONE JACOBI GS FBGS BJAC AS ML POLY +% %%%%%%%%%%% 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 POLY r 1-level, repeats previous. +6 ! Number of sweeps for smoother +1 ! degree for polynomial smoother +POLY_LOTTES_BETA ! Polynomial variant +% Fields to be added for POLY +% POLY_RHO_ESTIMATE Currently only POLY_RHO_EST_POWER +% POLY_RHO_ESTIMATE_ITERATIONS default = 20 +% POLY_RHO_BA set to value +% 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 -8 ! Inner solver sweeps (GS and JACOBI) -LLK ! AINV variant, ignored otherwise +L1-JACOBI ! 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) %%%%%%%%%%% Second smoother, always ignored for non-ML %%%%%%%%%%%%%%%% NONE ! Second (post) smoother, ignored if NONE -1 ! Number of sweeps for (post) smoother +6 ! Number of sweeps for (post) smoother +1 ! degree for polynomial smoother +POLY_LOTTES_BETA ! Polynomial variant 0 ! Number of overlap layers for AS preconditioner HALO ! AS restriction operator: NONE HALO NONE ! AS prolongation operator: NONE SUM AVG diff --git a/samples/advanced/pdegen/runs/amg_pde3d.inp b/samples/advanced/pdegen/runs/amg_pde3d.inp index 61e99a04..cb3b0819 100644 --- a/samples/advanced/pdegen/runs/amg_pde3d.inp +++ b/samples/advanced/pdegen/runs/amg_pde3d.inp @@ -1,7 +1,7 @@ %%%%%%%%%%% General arguments % Lines starting with % are ignored. CSR ! Storage format CSR COO JAD 0150 ! IDIM; domain size. Linear system size is IDIM**3 -CONST ! PDECOEFF: CONST, EXP, GAUSS Coefficients of the PDE +CONST ! PDECOEFF: CONST, EXP, BOX, GAUSS Coefficients of the PDE CG ! Iterative method: BiCGSTAB BiCGSTABL BiCG CG CGS FCG GCR RGMRES 2 ! ISTOPC 00500 ! ITMAX @@ -12,16 +12,16 @@ CG ! Iterative method: BiCGSTAB BiCGSTABL BiCG CG CGS F ML-VBM-VCYCLE-FBGS-D-BJAC ! Longer descriptive name for preconditioner (up to 20 chars) ML ! Preconditioner type: NONE JACOBI GS FBGS BJAC AS ML POLY % +%%%%%%%%%%% First smoother (for all levels but coarsest) %%%%%%%%%%%%%%%% +FBGS ! Smoother type JACOBI FBGS GS BWGS BJAC AS POLY r 1-level, repeats previous. +6 ! Number of sweeps for smoother +1 ! degree for polynomial smoother +POLY_LOTTES_BETA ! Polynomial variant % Fields to be added for POLY % POLY_RHO_ESTIMATE Currently only POLY_RHO_EST_POWER % POLY_RHO_ESTIMATE_ITERATIONS default = 20 % POLY_RHO_BA set to value % -%%%%%%%%%%% First smoother (for all levels but coarsest) %%%%%%%%%%%%%%%% -L1-JACOBI ! Smoother type JACOBI FBGS GS BWGS BJAC AS POLY r 1-level, repeats previous. -6 ! Number of sweeps for smoother -1 ! degree for polynomial smoother -POLY_LOTTES_BETA ! Polynomial variant 0 ! Number of overlap layers for AS preconditioner HALO ! AS restriction operator: NONE HALO NONE ! AS prolongation operator: NONE SUM AVG @@ -32,7 +32,7 @@ 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 %%%%%%%%%%%%%%%% -L1-JACOBI ! Second (post) smoother, ignored if NONE +NONE ! Second (post) smoother, ignored if NONE 6 ! Number of sweeps for (post) smoother 1 ! degree for polynomial smoother POLY_LOTTES_BETA ! Polynomial variant