Additional fields for rho(BA) estimate.

amgprec/amg_base_prec_type.F90

@@ -326,6 +326,9 @@ module amg_base_prec_type
   integer(psb_ipk_), parameter :: amg_poly_lottes_      = 0
   integer(psb_ipk_), parameter :: amg_poly_lottes_beta_ = 1
   integer(psb_ipk_), parameter :: amg_poly_new_         = 2
+  
+  integer(psb_ipk_), parameter :: amg_poly_rho_est_power_ = 0
+
   !
   ! Legal values for entry: amg_prec_status_
   !
@@ -572,6 +575,8 @@ contains
       val = amg_poly_lottes_beta_
     case('POLY_NEW')
       val = amg_poly_new_
+    case('POLY_RHO_EST_POWER')
+      val =  amg_poly_rho_est_power_
     case('A_NORMI')
       val = amg_max_norm_
     case('USER_CHOICE')

amgprec/amg_d_poly_smoother.f90

@@ -59,9 +59,11 @@ module amg_d_poly_smoother
     !    class(amg_d_base_solver_type), allocatable :: sv
     !
     integer(psb_ipk_)              :: pdegree, variant
+    integer(psb_ipk_)              :: rho_estimate=amg_poly_rho_est_power_
+    integer(psb_ipk_)              :: rho_estimate_iterations=10
     type(psb_dspmat_type), pointer :: pa => null()
     real(psb_dpk_), allocatable    :: poly_beta(:)
-    real(psb_dpk_)                 :: rho_ba
+    real(psb_dpk_)                 :: rho_ba = -done
   contains
     procedure, pass(sm) :: apply_v => amg_d_poly_smoother_apply_vect
 !!$    procedure, pass(sm) :: apply_a => amg_d_poly_smoother_apply
@@ -317,10 +319,11 @@ contains
     !
     ! Default: BJAC with no residual check
     !
-    sm%pdegree = 1
+    sm%pdegree      = 1
-    sm%rho_ba  = dzero
+    sm%rho_ba       = -done
-    sm%variant = amg_poly_lottes_
+    sm%variant      = amg_poly_lottes_
-
+    sm%rho_estimate = amg_poly_rho_est_power_
+    sm%rho_estimate_iterations = 20
     if (allocated(sm%sv)) then
       call sm%sv%default()
     end if

amgprec/impl/smoother/amg_d_poly_smoother_bld.f90

@@ -98,44 +98,50 @@ subroutine amg_d_poly_smoother_bld(a,desc_a,sm,info,amold,vmold,imold)
     goto 9999
   end if
   
-  if (.false.) then 
+!!$  if (.false.) then 
-    select type(ssv => sm%sv)
+!!$    select type(ssv => sm%sv)
-    class is(amg_d_l1_diag_solver_type)
+!!$    class is(amg_d_l1_diag_solver_type)
-      da  = a%arwsum(info)
+!!$      da  = a%arwsum(info)
-      dsv = ssv%dv%get_vect()
+!!$      dsv = ssv%dv%get_vect()
-      sm%rho_ba = maxval(da(1:n_row)*dsv(1:n_row))
+!!$      sm%rho_ba = maxval(da(1:n_row)*dsv(1:n_row))
-    class default
+!!$    class default
-      write(0,*) 'PolySmoother BUILD: only L1-Jacobi/L1-DIAG for now ',ssv%get_fmt()
+!!$      write(0,*) 'PolySmoother BUILD: only L1-Jacobi/L1-DIAG for now ',ssv%get_fmt()
+!!$      sm%rho_ba = done          
+!!$    end select
+!!$  else
+  if (sm%rho_ba <= dzero) then
+    select case(sm%rho_estimate)
+    case(amg_poly_rho_est_power_)
+      block
+        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,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,sm%rho_estimate_iterations
+          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
+      end block
+    case default
+      write(0,*) ' Unknown algorithm for RHO(BA) estimate, defaulting to a value of 1.0 '
       sm%rho_ba = done
     end select
-  else
-    block
-      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,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
   
   if (debug_level >= psb_debug_outer_) &

amgprec/impl/smoother/amg_d_poly_smoother_csetc.f90

@@ -55,6 +55,8 @@ subroutine amg_d_poly_smoother_csetc(sm,what,val,info,idx)
   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_d_base_smoother_type%set(what,val,info,idx=idx)
   end select

amgprec/impl/smoother/amg_d_poly_smoother_cseti.f90

@@ -64,6 +64,21 @@ subroutine amg_d_poly_smoother_cseti(sm,what,val,info,idx)
       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_d_base_smoother_type%set(what,val,info,idx=idx)
   end select

amgprec/impl/smoother/amg_d_poly_smoother_csetr.f90

@@ -54,8 +54,13 @@ subroutine amg_d_poly_smoother_csetr(sm,what,val,info,idx)
   call psb_erractionsave(err_act)
 
   select case(psb_toupper(what))
-  case('RHO_BA')
+  case('POLY_RHO_BA')
-    sm%rho_ba  = val
+    if ((dzero<val).and.(val<=done)) then 
+      sm%rho_ba  = val
+    else
+      write(0,*) 'Invalid choice for POLY_RHO_BA, defaulting to compute estimate'
+      sm%rho_ba = -done
+    end if
   case default
     call sm%amg_d_base_smoother_type%set(what,val,info,idx=idx)
   end select

samples/advanced/pdegen/runs/amg_pde3d.inp

@@ -4,13 +4,19 @@ CSR                         ! Storage format CSR COO JAD
 CONST                       ! PDECOEFF: CONST, EXP, GAUSS Coefficients of the PDE
 CG                          ! Iterative method: BiCGSTAB BiCGSTABL BiCG CG CGS FCG GCR RGMRES
 2                           ! ISTOPC
-00050                       ! ITMAX
+00500                       ! ITMAX
 1                           ! ITRACE
 30                          ! IRST (restart for RGMRES and BiCGSTABL)
 1.d-6                       ! EPS
 %%%%%%%%%%%  Main preconditioner choices %%%%%%%%%%%%%%%%
 ML-VBM-VCYCLE-FBGS-D-BJAC       ! Longer descriptive name for preconditioner (up to 20 chars)
 ML                        ! Preconditioner type: NONE JACOBI GS FBGS BJAC AS ML POLY
+%
+%  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