From bb262275a19b4350f2f5e3a74863b4ddffe00976 Mon Sep 17 00:00:00 2001 From: sfilippone Date: Fri, 17 Nov 2023 14:45:18 +0100 Subject: [PATCH] Temporary checkpoint, working version, to be investigated further. --- amgprec/amg_base_prec_type.F90 | 6 +- .../impl/smoother/amg_d_poly_smoother_bld.f90 | 27 ++++---- .../smoother/amg_d_poly_smoother_csetc.f90 | 2 +- .../smoother/amg_d_poly_smoother_cseti.f90 | 2 +- .../smoother/amg_d_poly_smoother_descr.f90 | 10 +-- samples/advanced/pdegen/amg_d_pde3d.F90 | 63 ++++++++++--------- samples/advanced/pdegen/runs/amg_pde3d.inp | 10 +-- 7 files changed, 63 insertions(+), 57 deletions(-) diff --git a/amgprec/amg_base_prec_type.F90 b/amgprec/amg_base_prec_type.F90 index 0bd44079..59b1acce 100644 --- a/amgprec/amg_base_prec_type.F90 +++ b/amgprec/amg_base_prec_type.F90 @@ -464,12 +464,12 @@ contains character(len=*), parameter :: name='amg_stringval' ! Local variable integer :: index_tab - character(len=15) ::string2 + character(len=128) ::string2 index_tab=index(string,char(9)) if (index_tab.NE.0) then - string2=string(1:index_tab-1) + string2=string(1:index_tab-1) else - string2=string + string2=string endif select case(psb_toupper(trim(string2))) case('NONE') diff --git a/amgprec/impl/smoother/amg_d_poly_smoother_bld.f90 b/amgprec/impl/smoother/amg_d_poly_smoother_bld.f90 index 2858f9ba..b307fef5 100644 --- a/amgprec/impl/smoother/amg_d_poly_smoother_bld.f90 +++ b/amgprec/impl/smoother/amg_d_poly_smoother_bld.f90 @@ -98,7 +98,7 @@ subroutine amg_d_poly_smoother_bld(a,desc_a,sm,info,amold,vmold,imold) goto 9999 end if - if (.true.) then + if (.false.) then select type(ssv => sm%sv) class is(amg_d_l1_diag_solver_type) da = a%arwsum(info) @@ -110,28 +110,31 @@ subroutine amg_d_poly_smoother_bld(a,desc_a,sm,info,amold,vmold,imold) end select else block - type(psb_d_vect_type) :: tq, tz,wv(2) - real(psb_dpk_) :: qnrm, lambda + type(psb_d_vect_type) :: tq, tt, tz,wv(2) + real(psb_dpk_) :: znrm, lambda real(psb_dpk_),allocatable :: work(:) integer(psb_ipk_) :: i, n_cols n_cols = desc_a%get_local_cols() allocate(work(4*n_cols)) call psb_geasb(tz,desc_a,info,mold=vmold,scratch=.true.) + call psb_geasb(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(done) - call psb_geasb(tq,desc_a,info,mold=vmold) - call psb_spmm(done,a,tq,dzero,tz,desc_a,info) ! z_1 = A q_0 - do i=1,10 - call sm%sv%apply_v(done,tz,dzero,tq,desc_a,'NoTrans',work,wv,info) ! q_k = M^{-1} q_k - qnrm = psb_genrmi(tq,desc_a,info) ! qnrm = |q_k|_inf - call tq%scal((done/qnrm)) ! q_k = q_k/qnrm - call psb_spmm(done,a,tq,dzero,tz,desc_a,info) ! z_{k=1} = A q_k - lambda = psb_gedot(tq,tz,desc_a,info) ! lambda = q_k^T z_{k+1} = q_k^T A q_k - write(0,*) 'BLD: lambda estimate ',i,lambda + call psb_geasb(tq,desc_a,info,mold=vmold) + call psb_spmm(done,a,tq,dzero,tt,desc_a,info) ! + call sm%sv%apply_v(done,tt,dzero,tz,desc_a,'NoTrans',work,wv,info) ! z_{k+1} = BA q_k + do i=1,20 + znrm = psb_genrm2(tz,desc_a,info) ! znrm = |z_k|_2 + call psb_geaxpby((done/znrm),tz,dzero,tq,desc_a,info) ! q_k = z_k/znrm + call psb_spmm(done,a,tq,dzero,tt,desc_a,info) ! t_{k+1} = BA q_k + call sm%sv%apply_v(done,tt,dzero,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 + sm%rho_ba = done end block end if diff --git a/amgprec/impl/smoother/amg_d_poly_smoother_csetc.f90 b/amgprec/impl/smoother/amg_d_poly_smoother_csetc.f90 index 0cc786ed..0daa387b 100644 --- a/amgprec/impl/smoother/amg_d_poly_smoother_csetc.f90 +++ b/amgprec/impl/smoother/amg_d_poly_smoother_csetc.f90 @@ -53,7 +53,7 @@ subroutine amg_d_poly_smoother_csetc(sm,what,val,info,idx) call psb_erractionsave(err_act) select case(psb_toupper(trim(what))) - case('POLY-VARIANT') + case('POLY_VARIANT') call sm%set(what,amg_stringval(val),info,idx=idx) case default call sm%amg_d_base_smoother_type%set(what,val,info,idx=idx) diff --git a/amgprec/impl/smoother/amg_d_poly_smoother_cseti.f90 b/amgprec/impl/smoother/amg_d_poly_smoother_cseti.f90 index d3db3891..0b116deb 100644 --- a/amgprec/impl/smoother/amg_d_poly_smoother_cseti.f90 +++ b/amgprec/impl/smoother/amg_d_poly_smoother_cseti.f90 @@ -61,7 +61,7 @@ subroutine amg_d_poly_smoother_cseti(sm,what,val,info,idx) 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_' + write(0,*) 'Invalid choice for POLY_VARIANT, defaulting to amg_poly_lottes_',val sm%variant = amg_poly_lottes_ end select case default diff --git a/amgprec/impl/smoother/amg_d_poly_smoother_descr.f90 b/amgprec/impl/smoother/amg_d_poly_smoother_descr.f90 index 51c198b1..0607064d 100644 --- a/amgprec/impl/smoother/amg_d_poly_smoother_descr.f90 +++ b/amgprec/impl/smoother/amg_d_poly_smoother_descr.f90 @@ -80,15 +80,17 @@ subroutine amg_d_poly_smoother_descr(sm,info,iout,coarse,prefix) 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_), ' 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 default write(iout_,*) trim(prefix_), ' variant: ','UNKNOWN???' end select - write(iout_,*) trim(prefix_), ' Degree: ',sm%pdegree - write(iout_,*) trim(prefix_), ' rho_ba: ',sm%rho_ba - if (allocated(sm%poly_beta)) write(iout_,*) trim(prefix_), ' Coefficients: ',sm%poly_beta(1:sm%pdegree) if (allocated(sm%sv)) then write(iout_,*) trim(prefix_), ' Local solver details:' call sm%sv%descr(info,iout_,coarse=coarse,prefix=prefix) diff --git a/samples/advanced/pdegen/amg_d_pde3d.F90 b/samples/advanced/pdegen/amg_d_pde3d.F90 index 927fb943..11494e64 100644 --- a/samples/advanced/pdegen/amg_d_pde3d.F90 +++ b/samples/advanced/pdegen/amg_d_pde3d.F90 @@ -80,7 +80,7 @@ program amg_d_pde3d implicit none ! input parameters - character(len=20) :: kmethd, ptype + character(len=24) :: kmethd, ptype character(len=5) :: afmt, pdecoeff integer(psb_ipk_) :: idim integer(psb_epk_) :: system_size @@ -120,21 +120,21 @@ program amg_d_pde3d ! preconditioner type character(len=40) :: descr ! verbose description of the prec - character(len=10) :: ptype ! preconditioner type + character(len=24) :: ptype ! preconditioner type 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=24) :: 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=24) :: aggr_prol ! aggregation type: SMOOTHED, NONSMOOTHED + character(len=24) :: par_aggr_alg ! parallel aggregation algorithm: DEC, SYMDEC + character(len=24) :: 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=24) :: aggr_ord ! ordering for aggregation: NATURAL, DEGREE + character(len=24) :: aggr_filter ! filtering: FILTER, NO_FILTER real(psb_dpk_) :: mncrratio ! minimum aggregation ratio real(psb_dpk_), allocatable :: athresv(:) ! smoothed aggregation threshold vector integer(psb_ipk_) :: thrvsz ! size of threshold vector @@ -142,43 +142,43 @@ program amg_d_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=24) :: smther ! (pre-)smoother type: BJAC, AS integer(psb_ipk_) :: jsweeps ! (pre-)smoother / 1-lev prec. sweeps integer(psb_ipk_) :: degree ! degree for polynomial smoother - character(len=16) :: pvariant ! Polynomial smoother variant + character(len=24) :: pvariant ! Polynomial smoother variant integer(psb_ipk_) :: novr ! number of overlap layers - character(len=16) :: restr ! restriction over application of AS - character(len=16) :: prol ! prolongation over application of AS - character(len=16) :: solve ! local subsolver type: ILU, MILU, ILUT, + character(len=24) :: restr ! restriction over application of AS + character(len=24) :: prol ! prolongation over application of AS + character(len=24) :: 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=24) :: variant ! AINV variant: LLK, etc integer(psb_ipk_) :: fill ! fill-in for incomplete LU factorization integer(psb_ipk_) :: invfill ! Inverse fill-in for INVK real(psb_dpk_) :: thr ! threshold for ILUT factorization ! AMG post-smoother; ignored by 1-lev preconditioner - character(len=16) :: smther2 ! post-smoother type: BJAC, AS + character(len=24) :: smther2 ! post-smoother type: BJAC, AS integer(psb_ipk_) :: jsweeps2 ! post-smoother sweeps integer(psb_ipk_) :: degree2 ! degree for polynomial smoother - character(len=16) :: pvariant2 ! Polynomial smoother variant + character(len=24) :: pvariant2 ! Polynomial smoother variant integer(psb_ipk_) :: novr2 ! number of overlap layers - character(len=16) :: restr2 ! restriction over application of AS - character(len=16) :: prol2 ! prolongation over application of AS - character(len=16) :: solve2 ! local subsolver type: ILU, MILU, ILUT, + character(len=24) :: restr2 ! restriction over application of AS + character(len=24) :: prol2 ! prolongation over application of AS + character(len=24) :: 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=24) :: variant2 ! AINV variant: LLK, etc integer(psb_ipk_) :: fill2 ! fill-in for incomplete LU factorization integer(psb_ipk_) :: invfill2 ! Inverse fill-in for INVK real(psb_dpk_) :: thr2 ! threshold for ILUT factorization ! coarsest-level solver - character(len=16) :: cmat ! coarsest matrix layout: REPL, DIST - character(len=16) :: csolve ! coarsest-lev solver: BJAC, SLUDIST (distr. + character(len=24) :: cmat ! coarsest matrix layout: REPL, DIST + character(len=24) :: 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=24) :: csbsolve ! coarsest-lev local subsolver: ILU, ILUT, ! MILU, UMF, MUMPS, SLU integer(psb_ipk_) :: cfill ! fill-in for incomplete LU factorization real(psb_dpk_) :: cthres ! threshold for ILUT factorization @@ -200,7 +200,7 @@ program amg_d_pde3d ! other variables integer(psb_ipk_) :: info, i, k - character(len=20) :: name,ch_err + character(len=24) :: name,ch_err type(psb_d_csr_sparse_mat) :: amold info=psb_success_ @@ -294,7 +294,7 @@ program amg_d_pde3d 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%variant, 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) @@ -343,7 +343,8 @@ program amg_d_pde3d call prec%set('smoother_type', p_choice%smther, info) call prec%set('smoother_sweeps', p_choice%jsweeps, info) call prec%set('poly_degree', p_choice%degree, info) - call prec%set('poly_variant', p_choice%variant, info) + write(0,*) 'pvariant :',p_choice%pvariant + 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') @@ -375,7 +376,7 @@ program amg_d_pde3d 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%variant2, 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 @@ -592,7 +593,7 @@ contains 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) ! Degree of Polynomial smoother - call read_data(prec%variant,inp_unit) ! variant for Polynomial + call read_data(prec%pvariant,inp_unit) ! variant for Polynomial call read_data(prec%novr,inp_unit) ! number of overlap layers call read_data(prec%restr,inp_unit) ! restriction over application of AS call read_data(prec%prol,inp_unit) ! prolongation over application of AS @@ -606,7 +607,7 @@ contains 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) ! Degree of Polynomial smoother - call read_data(prec%variant2,inp_unit) ! Polynomial smoother variant + call read_data(prec%pvariant2,inp_unit) ! Polynomial smoother variant call read_data(prec%novr2,inp_unit) ! number of overlap layers call read_data(prec%restr2,inp_unit) ! restriction over application of AS call read_data(prec%prol2,inp_unit) ! prolongation over application of AS @@ -678,7 +679,7 @@ contains call psb_bcast(ctxt,prec%smther) call psb_bcast(ctxt,prec%jsweeps) call psb_bcast(ctxt,prec%degree) - call psb_bcast(ctxt,prec%variant) + call psb_bcast(ctxt,prec%pvariant) call psb_bcast(ctxt,prec%novr) call psb_bcast(ctxt,prec%restr) call psb_bcast(ctxt,prec%prol) @@ -692,7 +693,7 @@ contains call psb_bcast(ctxt,prec%smther2) call psb_bcast(ctxt,prec%jsweeps2) call psb_bcast(ctxt,prec%degree2) - call psb_bcast(ctxt,prec%variant2) + call psb_bcast(ctxt,prec%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/runs/amg_pde3d.inp b/samples/advanced/pdegen/runs/amg_pde3d.inp index 4a0e605c..95de34a6 100644 --- a/samples/advanced/pdegen/runs/amg_pde3d.inp +++ b/samples/advanced/pdegen/runs/amg_pde3d.inp @@ -1,6 +1,6 @@ %%%%%%%%%%% General arguments % Lines starting with % are ignored. CSR ! Storage format CSR COO JAD -0050 ! IDIM; domain size. Linear system size is IDIM**3 +0150 ! IDIM; domain size. Linear system size is IDIM**3 CONST ! PDECOEFF: CONST, EXP, GAUSS Coefficients of the PDE BICGSTAB ! Iterative method: BiCGSTAB BiCGSTABL BiCG CG CGS FCG GCR RGMRES 2 ! ISTOPC @@ -12,8 +12,8 @@ BICGSTAB ! 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 %%%%%%%%%%% First smoother (for all levels but coarsest) %%%%%%%%%%%%%%%% -FBGS ! Smoother type JACOBI FBGS GS BWGS BJAC AS. For 1-level, repeats previous. -1 ! Number of sweeps for smoother +L1-JACOBI ! Smoother type JACOBI FBGS GS BWGS BJAC AS. For 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 @@ -26,8 +26,8 @@ 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 %%%%%%%%%%%%%%%% -FBGS ! Second (post) smoother, ignored if NONE -1 ! Number of sweeps for (post) smoother +L1-JACOBI ! Second (post) smoother, ignored if NONE +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