Fix input file

Makefile

@@ -3,9 +3,9 @@ include Make.inc
 
 all:  objs lib
 
-objs: amgp cbnd
+objs: libdir amgp cbnd
 
-lib: libdir objs
+lib:  objs
 	cd amgprec && $(MAKE) lib
 	cd cbind && $(MAKE)  lib
 

amgprec/amg_base_prec_type.F90

@@ -288,11 +288,12 @@ module amg_base_prec_type
   !
   ! Legal values for entry: amg_aggr_prol_
   !
-  integer(psb_ipk_), parameter :: amg_no_smooth_   = 0
-  integer(psb_ipk_), parameter :: amg_smooth_prol_ = 1
-  integer(psb_ipk_), parameter :: amg_min_energy_  = 2
+  integer(psb_ipk_), parameter :: amg_no_smooth_      = 0
+  integer(psb_ipk_), parameter :: amg_smooth_prol_    = 1
+  integer(psb_ipk_), parameter :: amg_l1_smooth_prol_ = 2
+  integer(psb_ipk_), parameter :: amg_min_energy_     = 3
   ! Disabling min_energy  for the time being.
-  integer(psb_ipk_), parameter :: amg_max_aggr_prol_=amg_smooth_prol_
+  integer(psb_ipk_), parameter :: amg_max_aggr_prol_= amg_l1_smooth_prol_
   !
   ! Legal values for entry: amg_aggr_filter_
   !
@@ -323,10 +324,10 @@ module amg_base_prec_type
   !
   ! Legal values for entry: amg_poly_variant_
   !
-  integer(psb_ipk_), parameter :: amg_poly_lottes_      = 0
-  integer(psb_ipk_), parameter :: amg_poly_lottes_beta_ = 1
-  integer(psb_ipk_), parameter :: amg_poly_new_         = 2
-  integer(psb_ipk_), parameter :: amg_poly_dbg_         = 8
+  integer(psb_ipk_), parameter :: amg_cheb_4_     = 0
+  integer(psb_ipk_), parameter :: amg_cheb_4_opt_ = 1
+  integer(psb_ipk_), parameter :: amg_cheb_1_opt_ = 2
+  integer(psb_ipk_), parameter :: amg_poly_dbg_   = 8
   
   integer(psb_ipk_), parameter :: amg_poly_rho_est_power_ = 0
 
@@ -376,8 +377,8 @@ module amg_base_prec_type
   character(len=19), parameter, private :: &
        &  eigen_estimates(0:0)=(/'infinity norm     '/)
   character(len=15), parameter, private :: &
-       &  aggr_prols(0:3)=(/'unsmoothed    ','smoothed      ',&
-       &           'min energy    ','bizr. smoothed'/)
+       &  aggr_prols(0:4)=(/'unsmoothed    ','smoothed      ',&
+       &  'l1-smoothed   ','min energy    ','bizr. smoothed'/)
   character(len=15), parameter, private :: &
        &  aggr_filters(0:1)=(/'no filtering  ','filtering     '/)
   character(len=15), parameter, private :: &
@@ -548,6 +549,8 @@ contains
       val = amg_no_smooth_
     case('SMOOTHED')
       val = amg_smooth_prol_
+    case('L1-SMOOTHED','L1SMOOTHED')
+      val = amg_l1_smooth_prol_
     case('MINENERGY')
       val = amg_min_energy_
     case('NOPREC')
@@ -570,12 +573,12 @@ contains
       val = amg_as_
     case('POLY')
       val = amg_poly_
-    case('POLY_LOTTES')
-      val = amg_poly_lottes_
-    case('POLY_LOTTES_BETA')
-      val = amg_poly_lottes_beta_
-    case('POLY_NEW')
-      val = amg_poly_new_
+    case('CHEB_4')
+      val = amg_cheb_4_
+    case('CHEB_4_OPT')
+      val = amg_cheb_4_opt_
+    case('CHEB_1_OPT')
+      val = amg_cheb_1_opt_
     case('POLY_DBG')
       val = amg_poly_dbg_
     case('POLY_RHO_EST_POWER')

amgprec/amg_c_inner_mod.f90

@@ -109,11 +109,12 @@ module amg_c_inner_mod
   end interface amg_map_to_tprol
 
   abstract interface
-    subroutine amg_caggrmat_var_bld(a,desc_a,ilaggr,nlaggr,parms,&
+    subroutine amg_caggrmat_var_bld(dol1smoothing,a,desc_a,ilaggr,nlaggr,parms,&
          & ac,desc_ac,op_prol,op_restr,t_prol,info)
       import :: psb_cspmat_type, psb_desc_type, psb_spk_, psb_ipk_, psb_lpk_, psb_lcspmat_type
       import ::  amg_c_onelev_type, amg_sml_parms
       implicit none 
+      integer(psb_ipk_), intent(in)               :: dol1smoothing
       type(psb_cspmat_type), intent(in)         :: a
       type(psb_desc_type), intent(inout)          :: desc_a
       integer(psb_lpk_), intent(inout)            :: ilaggr(:), nlaggr(:)

amgprec/amg_d_inner_mod.f90

@@ -109,11 +109,12 @@ module amg_d_inner_mod
   end interface amg_map_to_tprol
 
   abstract interface
-    subroutine amg_daggrmat_var_bld(a,desc_a,ilaggr,nlaggr,parms,&
+    subroutine amg_daggrmat_var_bld(dol1smoothing,a,desc_a,ilaggr,nlaggr,parms,&
          & ac,desc_ac,op_prol,op_restr,t_prol,info)
       import :: psb_dspmat_type, psb_desc_type, psb_dpk_, psb_ipk_, psb_lpk_, psb_ldspmat_type
       import ::  amg_d_onelev_type, amg_dml_parms
       implicit none 
+      integer(psb_ipk_), intent(in)               :: dol1smoothing
       type(psb_dspmat_type), intent(in)         :: a
       type(psb_desc_type), intent(inout)          :: desc_a
       integer(psb_lpk_), intent(inout)            :: ilaggr(:), nlaggr(:)

amgprec/amg_d_parmatch_aggregator_mod.F90

@@ -244,11 +244,12 @@ module amg_d_parmatch_aggregator_mod
   end interface
 
   interface
-    subroutine amg_d_parmatch_unsmth_bld(ag,a,desc_a,ilaggr,nlaggr,parms,&
+    subroutine amg_d_parmatch_unsmth_bld(dol1smoothing,ag,a,desc_a,ilaggr,nlaggr,parms,&
          & ac,desc_ac,op_prol,op_restr,t_prol,info)
       import :: amg_d_parmatch_aggregator_type, psb_desc_type, psb_dspmat_type,&
            & psb_ldspmat_type, psb_dpk_, psb_ipk_, psb_lpk_, amg_dml_parms, amg_daggr_data
       implicit none
+      integer(psb_ipk_), intent(in) :: dol1smoothing
       class(amg_d_parmatch_aggregator_type), target, intent(inout) :: ag
       type(psb_dspmat_type), intent(in)       :: a
       type(psb_desc_type), intent(inout)      :: desc_a
@@ -262,11 +263,12 @@ module amg_d_parmatch_aggregator_mod
   end interface
 
   interface
-    subroutine amg_d_parmatch_smth_bld(ag,a,desc_a,ilaggr,nlaggr,parms,&
+    subroutine amg_d_parmatch_smth_bld(dol1smoothing,ag,a,desc_a,ilaggr,nlaggr,parms,&
          & ac,desc_ac,op_prol,op_restr,t_prol,info)
       import :: amg_d_parmatch_aggregator_type, psb_desc_type, psb_dspmat_type,&
            & psb_ldspmat_type, psb_dpk_, psb_ipk_, psb_lpk_, amg_dml_parms, amg_daggr_data
       implicit none
+      integer(psb_ipk_), intent(in) :: dol1smoothing
       class(amg_d_parmatch_aggregator_type), target, intent(inout) :: ag
       type(psb_dspmat_type), intent(in)       :: a
       type(psb_desc_type), intent(inout)      :: desc_a

amgprec/amg_d_poly_smoother.f90

@@ -322,7 +322,7 @@ contains
     !
     sm%pdegree      = 1
     sm%rho_ba       = -done
-    sm%variant      = amg_poly_lottes_
+    sm%variant      = amg_cheb_4_
     sm%rho_estimate = amg_poly_rho_est_power_
     sm%rho_estimate_iterations = 20
     if (allocated(sm%sv)) then

amgprec/amg_s_inner_mod.f90

@@ -109,11 +109,12 @@ module amg_s_inner_mod
   end interface amg_map_to_tprol
 
   abstract interface
-    subroutine amg_saggrmat_var_bld(a,desc_a,ilaggr,nlaggr,parms,&
+    subroutine amg_saggrmat_var_bld(dol1smoothing,a,desc_a,ilaggr,nlaggr,parms,&
          & ac,desc_ac,op_prol,op_restr,t_prol,info)
       import :: psb_sspmat_type, psb_desc_type, psb_spk_, psb_ipk_, psb_lpk_, psb_lsspmat_type
       import ::  amg_s_onelev_type, amg_sml_parms
       implicit none 
+      integer(psb_ipk_), intent(in)               :: dol1smoothing
       type(psb_sspmat_type), intent(in)         :: a
       type(psb_desc_type), intent(inout)          :: desc_a
       integer(psb_lpk_), intent(inout)            :: ilaggr(:), nlaggr(:)

amgprec/amg_s_matchboxp_mod.f90

@@ -272,7 +272,9 @@ contains
           write(0,*) 'Impossible: mate(k) > nc'
           cycle
         else
+
           if (ilaggr(k) == ilaggr_neginit) then
+
             wk   = w(k)
             widx = w(idx)
             wmax = max(abs(wk),abs(widx))

amgprec/amg_s_parmatch_aggregator_mod.F90

@@ -244,11 +244,12 @@ module amg_s_parmatch_aggregator_mod
   end interface
 
   interface
-    subroutine amg_s_parmatch_unsmth_bld(ag,a,desc_a,ilaggr,nlaggr,parms,&
+    subroutine amg_s_parmatch_unsmth_bld(dol1smoothing,ag,a,desc_a,ilaggr,nlaggr,parms,&
          & ac,desc_ac,op_prol,op_restr,t_prol,info)
       import :: amg_s_parmatch_aggregator_type, psb_desc_type, psb_sspmat_type,&
            & psb_lsspmat_type, psb_dpk_, psb_ipk_, psb_lpk_, amg_sml_parms, amg_saggr_data
       implicit none
+      integer(psb_ipk_), intent(in) :: dol1smoothing
       class(amg_s_parmatch_aggregator_type), target, intent(inout) :: ag
       type(psb_sspmat_type), intent(in)       :: a
       type(psb_desc_type), intent(inout)      :: desc_a
@@ -262,11 +263,12 @@ module amg_s_parmatch_aggregator_mod
   end interface
 
   interface
-    subroutine amg_s_parmatch_smth_bld(ag,a,desc_a,ilaggr,nlaggr,parms,&
+    subroutine amg_s_parmatch_smth_bld(dol1smoothing,ag,a,desc_a,ilaggr,nlaggr,parms,&
          & ac,desc_ac,op_prol,op_restr,t_prol,info)
       import :: amg_s_parmatch_aggregator_type, psb_desc_type, psb_sspmat_type,&
            & psb_lsspmat_type, psb_dpk_, psb_ipk_, psb_lpk_, amg_sml_parms, amg_saggr_data
       implicit none
+      integer(psb_ipk_), intent(in) :: dol1smoothing
       class(amg_s_parmatch_aggregator_type), target, intent(inout) :: ag
       type(psb_sspmat_type), intent(in)       :: a
       type(psb_desc_type), intent(inout)      :: desc_a

amgprec/amg_s_poly_smoother.f90

@@ -322,7 +322,7 @@ contains
     !
     sm%pdegree      = 1
     sm%rho_ba       = -sone
-    sm%variant      = amg_poly_lottes_
+    sm%variant      = amg_cheb_4_
     sm%rho_estimate = amg_poly_rho_est_power_
     sm%rho_estimate_iterations = 20
     if (allocated(sm%sv)) then

amgprec/amg_z_inner_mod.f90

@@ -109,11 +109,12 @@ module amg_z_inner_mod
   end interface amg_map_to_tprol
 
   abstract interface
-    subroutine amg_zaggrmat_var_bld(a,desc_a,ilaggr,nlaggr,parms,&
+    subroutine amg_zaggrmat_var_bld(dol1smoothing,a,desc_a,ilaggr,nlaggr,parms,&
          & ac,desc_ac,op_prol,op_restr,t_prol,info)
       import :: psb_zspmat_type, psb_desc_type, psb_dpk_, psb_ipk_, psb_lpk_, psb_lzspmat_type
       import ::  amg_z_onelev_type, amg_dml_parms
       implicit none 
+      integer(psb_ipk_), intent(in)               :: dol1smoothing
       type(psb_zspmat_type), intent(in)         :: a
       type(psb_desc_type), intent(inout)          :: desc_a
       integer(psb_lpk_), intent(inout)            :: ilaggr(:), nlaggr(:)

amgprec/impl/aggregator/MatchBoxPC.h

@@ -82,6 +82,8 @@ const int BundleTag = 9;  // Predefined tag
 
 static vector<MilanLongInt> DEFAULT_VECTOR;
 
+#if !defined(SERIAL_MPI)
+
 // MPI type map
 template <typename T>
 MPI_Datatype TypeMap();
@@ -93,6 +95,7 @@ template <>
 inline MPI_Datatype TypeMap<double>() { return MPI_DOUBLE; }
 template <>
 inline MPI_Datatype TypeMap<float>() { return MPI_FLOAT; }
+#endif
 
 #ifdef __cplusplus
 extern "C"

amgprec/impl/aggregator/amg_c_dec_aggregator_mat_bld.f90

@@ -177,23 +177,24 @@ subroutine  amg_c_dec_aggregator_mat_bld(ag,parms,a,desc_a,ilaggr,nlaggr,&
   select case (parms%aggr_prol)
   case (amg_no_smooth_) 
 
-    call amg_caggrmat_nosmth_bld(a,desc_a,ilaggr,nlaggr,&
-         & parms,ac,desc_ac,op_prol,op_restr,t_prol,info)
+    call amg_caggrmat_nosmth_bld(parms%aggr_prol,a,desc_a,ilaggr,&
+            nlaggr,parms,ac,desc_ac,op_prol,op_restr,t_prol,info)
 
-  case(amg_smooth_prol_) 
+  case(amg_smooth_prol_,amg_l1_smooth_prol_) 
 
-    call amg_caggrmat_smth_bld(a,desc_a,ilaggr,nlaggr, &
-         & parms,ac,desc_ac,op_prol,op_restr,t_prol,info)
+    call amg_caggrmat_smth_bld(parms%aggr_prol,a,desc_a,&
+            ilaggr,nlaggr,parms,ac,desc_ac,op_prol,&
+            op_restr,t_prol,info)
 
 !!$  case(amg_biz_prol_) 
 !!$
 !!$    call amg_caggrmat_biz_bld(a,desc_a,ilaggr,nlaggr, &
 !!$         & parms,ac,desc_ac,op_prol,op_restr,t_prol,info)
-
+  
   case(amg_min_energy_) 
 
-    call amg_caggrmat_minnrg_bld(a,desc_a,ilaggr,nlaggr, &
-         & parms,ac,desc_ac,op_prol,op_restr,t_prol,info)
+    call amg_caggrmat_minnrg_bld(parms%aggr_prol,a,desc_a,ilaggr,&
+            nlaggr,parms,ac,desc_ac,op_prol,op_restr,t_prol,info)
 
   case default
     info = psb_err_internal_error_

amgprec/impl/aggregator/amg_c_soc1_map_bld.F90

@@ -250,7 +250,7 @@ subroutine amg_c_soc1_map_bld(iorder,theta,clean_zeros,a,desc_a,nlaggr,ilaggr,in
               ! we will not reset. 
               if (j>nr) cycle step1
               if (ilaggr(j) > 0) cycle step1
-              if (abs(val(k)) > theta*sqrt(abs(diag(i)*diag(j)))) then
+              if ((abs(val(k)) > theta*sqrt(abs(diag(i)*diag(j)))).and.(diag(i).ne.czero)) then
                 ip = ip + 1
                 icol(ip) = icol(k)
               end if
@@ -357,7 +357,7 @@ subroutine amg_c_soc1_map_bld(iorder,theta,clean_zeros,a,desc_a,nlaggr,ilaggr,in
       do k=1, nz
         j   = icol(k)
         if ((1<=j).and.(j<=nr)) then 
-          if (abs(val(k)) > theta*sqrt(abs(diag(i)*diag(j)))) then
+          if ((abs(val(k)) > theta*sqrt(abs(diag(i)*diag(j)))).and.(diag(i).ne.czero)) then
             ip = ip + 1
             icol(ip) = icol(k)
           end if
@@ -545,4 +545,3 @@ subroutine amg_c_soc1_map_bld(iorder,theta,clean_zeros,a,desc_a,nlaggr,ilaggr,in
   return
 
 end subroutine amg_c_soc1_map_bld
-

amgprec/impl/aggregator/amg_caggrmat_minnrg_bld.f90

@@ -69,6 +69,7 @@
 !
 !
 ! Arguments:
+!    dol1smoothing - fictitious integer argument, it is not used inside
 !    a          -  type(psb_cspmat_type), input.     
 !                  The sparse matrix structure containing the local part of
 !                  the fine-level matrix.
@@ -104,8 +105,8 @@
 !                  Error code.
 !
 !
-subroutine amg_caggrmat_minnrg_bld(a,desc_a,ilaggr,nlaggr,parms,&
-     & ac,desc_ac,op_prol,op_restr,t_prol,info)
+subroutine amg_caggrmat_minnrg_bld(dol1smoothing,a,desc_a,ilaggr,nlaggr,&
+                parms,ac,desc_ac,op_prol,op_restr,t_prol,info)
   use psb_base_mod
   use amg_base_prec_type
   use amg_c_inner_mod, amg_protect_name => amg_caggrmat_minnrg_bld
@@ -113,6 +114,7 @@ subroutine amg_caggrmat_minnrg_bld(a,desc_a,ilaggr,nlaggr,parms,&
   implicit none
 
   ! Arguments
+  integer(psb_ipk_), intent(in)                 :: dol1smoothing
   type(psb_cspmat_type), intent(in)           :: a
   type(psb_desc_type), intent(inout)            :: desc_a
   integer(psb_lpk_), intent(inout)              :: ilaggr(:), nlaggr(:)
@@ -171,6 +173,13 @@ subroutine amg_caggrmat_minnrg_bld(a,desc_a,ilaggr,nlaggr,parms,&
 
   filter_mat = (parms%aggr_filter == amg_filter_mat_)
 
+  if (dol1smoothing.ne.amg_no_smooth_) then
+    info=psb_err_fatal_; 
+    call psb_errpush(info,name,a_err='Are you trying to smooth an unsmoothed aggregation?')
+    goto 9999
+  end if
+
+
   !NEEDS TO BE REWORKED !!
   
   ! naggr: number of local aggregates
@@ -300,7 +309,7 @@ subroutine amg_caggrmat_minnrg_bld(a,desc_a,ilaggr,nlaggr,parms,&
 !!$        endif
 !!$      enddo
 !!$      if (jd == -1) then 
-!!$        write(0,*) 'Wrong input: we need the diagonal!!!!', i
+!!$                write(0,*) name,': Warning: there is no diagonal element', i
 !!$      else
 !!$        acsrf%val(jd)=acsrf%val(jd)-tmp
 !!$      end if

amgprec/impl/aggregator/amg_caggrmat_nosmth_bld.f90

@@ -94,10 +94,11 @@
 !               
 !    info       -  integer, output.
 !                  Error code.
+!    dol1smoothing - optional, this is here just for interfacing reasons. It is not used by the
+!                    code
 !
-!
-subroutine amg_caggrmat_nosmth_bld(a,desc_a,ilaggr,nlaggr,parms,&
-     & ac,desc_ac,op_prol,op_restr,t_prol,info)
+subroutine amg_caggrmat_nosmth_bld(dol1smoothing,a,desc_a,ilaggr,nlaggr,&
+                parms,ac,desc_ac,op_prol,op_restr,t_prol,info)
   use psb_base_mod
   use amg_base_prec_type
   use amg_c_inner_mod, amg_protect_name => amg_caggrmat_nosmth_bld
@@ -105,6 +106,7 @@ subroutine amg_caggrmat_nosmth_bld(a,desc_a,ilaggr,nlaggr,parms,&
   implicit none
 
   ! Arguments
+  integer(psb_ipk_), intent(in)              :: dol1smoothing 
   type(psb_cspmat_type), intent(in)        :: a
   type(psb_desc_type), intent(inout)         :: desc_a
   integer(psb_lpk_), intent(inout)           :: ilaggr(:), nlaggr(:)
@@ -137,6 +139,12 @@ subroutine amg_caggrmat_nosmth_bld(a,desc_a,ilaggr,nlaggr,parms,&
 
   ctxt = desc_a%get_context()
   call psb_info(ctxt, me, np)
+  if (dol1smoothing.ne.amg_no_smooth_) then
+    info=psb_err_fatal_;         
+    call psb_errpush(info,name,a_err='Are you trying to smooth an unsmoothed aggregation?')
+    goto 9999
+  end if
+
   nglob = desc_a%get_global_rows()
   nrow  = desc_a%get_local_rows()
   ncol  = desc_a%get_local_cols()

amgprec/impl/aggregator/amg_caggrmat_smth_bld.f90

@@ -69,6 +69,8 @@
 !
 !
 ! Arguments:
+!    dol1smooth -  Integer taking the type of smoother that has to be used
+!                  on the tentative prolongator
 !    a          -  type(psb_cspmat_type), input.     
 !                  The sparse matrix structure containing the local part of
 !                  the fine-level matrix.
@@ -102,8 +104,8 @@
 !    info       -  integer, output.
 !                  Error code.
 !
-subroutine amg_caggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,&
-     & ac,desc_ac,op_prol,op_restr,t_prol,info)
+subroutine amg_caggrmat_smth_bld(dol1smoothing,a,desc_a,ilaggr,nlaggr,&
+                parms,ac,desc_ac,op_prol,op_restr,t_prol,info)
   use psb_base_mod
   use amg_base_prec_type
   use amg_c_inner_mod, amg_protect_name => amg_caggrmat_smth_bld
@@ -112,6 +114,7 @@ subroutine amg_caggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,&
   implicit none
 
   ! Arguments
+  integer(psb_ipk_), intent(in)            :: dol1smoothing
   type(psb_cspmat_type), intent(in)      :: a
   type(psb_desc_type), intent(inout)       :: desc_a
   integer(psb_lpk_), intent(inout)         :: ilaggr(:), nlaggr(:)
@@ -132,7 +135,7 @@ subroutine amg_caggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,&
   type(psb_c_coo_sparse_mat) :: coo_prol, coo_restr
   type(psb_c_csr_sparse_mat) :: acsr1,  acsrf, csr_prol, acsr
   complex(psb_spk_), allocatable :: adiag(:)
-  real(psb_spk_), allocatable :: arwsum(:)
+  real(psb_spk_), allocatable :: arwsum(:),l1rwsum(:)
   integer(psb_ipk_)  :: ierr(5)
   logical            :: filter_mat
   integer(psb_ipk_)            :: debug_level, debug_unit, err_act
@@ -141,6 +144,7 @@ subroutine amg_caggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,&
   logical, parameter :: debug_new=.false.
   character(len=80) :: filename
   logical, parameter :: do_timings=.false.
+  logical :: do_l1correction=.false.
   integer(psb_ipk_), save :: idx_spspmm=-1, idx_phase1=-1, idx_gtrans=-1, idx_phase2=-1, idx_refine=-1
   integer(psb_ipk_), save :: idx_phase3=-1, idx_cdasb=-1, idx_ptap=-1
 
@@ -173,6 +177,9 @@ subroutine amg_caggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,&
   if ((do_timings).and.(idx_ptap==-1)) &
        & idx_ptap = psb_get_timer_idx("DEC_SMTH_BLD: ptap_bld  ")
 
+  ! check if we have to use Jacobi or l1-Jacobi to smooth the tentative prolongator
+  if (dol1smoothing.eq.amg_l1_smooth_prol_) do_l1correction=.true.
+
 
   nglob = desc_a%get_global_rows()
   nrow  = desc_a%get_local_rows()
@@ -200,6 +207,24 @@ subroutine amg_caggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,&
   if (info == psb_success_) &
        & call psb_halo(adiag,desc_a,info)
   if (info == psb_success_) call a%cp_to(acsr)
+  !
+  ! Do the l1-correction on the diagonal if it is requested
+  !
+  if (do_l1correction) then
+    allocate(l1rwsum(nrow))
+    call acsr%arwsum(l1rwsum)
+    if (info == psb_success_) &
+      & call psb_realloc(ncol,l1rwsum,info)
+    if (info == psb_success_) &
+      & call psb_halo(l1rwsum,desc_a,info)
+    ! \tilde{D}_{i,i} = \sum_{j \ne i} |a_{i,j}|
+    !$OMP parallel do private(i) schedule(static)
+    do i=1,size(adiag)
+      adiag(i) = adiag(i) + l1rwsum(i) - abs(adiag(i))
+    end do
+    !$OMP end parallel do
+  end if
+
 
   if(info /= psb_success_) then
     call psb_errpush(psb_err_from_subroutine_,name,a_err='sp_getdiag')
@@ -230,6 +255,7 @@ subroutine amg_caggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,&
 
       enddo
       if (jd == -1) then 
+        ! if (.not.do_l1correction) 
         write(0,*) 'Wrong input: we need the diagonal!!!!', i
       else
         acsrf%val(jd)=acsrf%val(jd)-tmp
@@ -252,6 +278,10 @@ subroutine amg_caggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,&
   !$OMP end parallel do 
   if (parms%aggr_omega_alg == amg_eig_est_) then 
 
+    !if (do_l1correction) then
+    !  ! For l1-Jacobi this can be estimated with 1
+    !  parms%aggr_omega_val = done
+    !
     if (parms%aggr_eig == amg_max_norm_) then 
       allocate(arwsum(nrow))
       call acsr%arwsum(arwsum)      
@@ -259,7 +289,6 @@ subroutine amg_caggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,&
       call psb_amx(ctxt,anorm)
       omega = 4.d0/(3.d0*anorm)
       parms%aggr_omega_val = omega 
-
     else 
       info = psb_err_internal_error_
       call psb_errpush(info,name,a_err='invalid amg_aggr_eig_')
@@ -322,6 +351,7 @@ subroutine amg_caggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,&
   if (debug_level >= psb_debug_outer_) &
        & write(debug_unit,*) me,' ',trim(name),&
        & 'Done smooth_aggregate '
+  if (allocated(l1rwsum)) deallocate(l1rwsum)
   call psb_erractionrestore(err_act)
   return
 

amgprec/impl/aggregator/amg_d_dec_aggregator_mat_bld.f90

@@ -177,23 +177,24 @@ subroutine  amg_d_dec_aggregator_mat_bld(ag,parms,a,desc_a,ilaggr,nlaggr,&
   select case (parms%aggr_prol)
   case (amg_no_smooth_) 
 
-    call amg_daggrmat_nosmth_bld(a,desc_a,ilaggr,nlaggr,&
-         & parms,ac,desc_ac,op_prol,op_restr,t_prol,info)
+    call amg_daggrmat_nosmth_bld(parms%aggr_prol,a,desc_a,ilaggr,&
+            nlaggr,parms,ac,desc_ac,op_prol,op_restr,t_prol,info)
 
-  case(amg_smooth_prol_) 
+  case(amg_smooth_prol_,amg_l1_smooth_prol_) 
 
-    call amg_daggrmat_smth_bld(a,desc_a,ilaggr,nlaggr, &
-         & parms,ac,desc_ac,op_prol,op_restr,t_prol,info)
+    call amg_daggrmat_smth_bld(parms%aggr_prol,a,desc_a,&
+            ilaggr,nlaggr,parms,ac,desc_ac,op_prol,&
+            op_restr,t_prol,info)
 
 !!$  case(amg_biz_prol_) 
 !!$
 !!$    call amg_daggrmat_biz_bld(a,desc_a,ilaggr,nlaggr, &
 !!$         & parms,ac,desc_ac,op_prol,op_restr,t_prol,info)
-
+  
   case(amg_min_energy_) 
 
-    call amg_daggrmat_minnrg_bld(a,desc_a,ilaggr,nlaggr, &
-         & parms,ac,desc_ac,op_prol,op_restr,t_prol,info)
+    call amg_daggrmat_minnrg_bld(parms%aggr_prol,a,desc_a,ilaggr,&
+            nlaggr,parms,ac,desc_ac,op_prol,op_restr,t_prol,info)
 
   case default
     info = psb_err_internal_error_

amgprec/impl/aggregator/amg_d_parmatch_aggregator_mat_bld.F90

@@ -184,20 +184,23 @@ subroutine  amg_d_parmatch_aggregator_mat_bld(ag,parms,a,desc_a,ilaggr,nlaggr,&
   !
   select case (parms%aggr_prol)
   case (amg_no_smooth_)
-    call amg_d_parmatch_unsmth_bld(ag,a,desc_a,ilaggr,nlaggr,parms,&
-         & ac,desc_ac,op_prol,op_restr,t_prol,info)
+    call amg_d_parmatch_unsmth_bld(parms%aggr_prol,ag,a,desc_a,&
+            ilaggr,nlaggr,parms,ac,desc_ac,op_prol,op_restr,&
+            t_prol,info)
 
-  case(amg_smooth_prol_)
-    call amg_d_parmatch_smth_bld(ag,a,desc_a,ilaggr,nlaggr,parms,&
-         & ac,desc_ac,op_prol,op_restr,t_prol,info)
+  case(amg_smooth_prol_,amg_l1_smooth_prol_)
+    call amg_d_parmatch_smth_bld(parms%aggr_prol,ag,a,desc_a,&
+            ilaggr,nlaggr,parms,ac,desc_ac,op_prol,op_restr,&
+            t_prol,info)
 
 !!$  case(amg_biz_prol_)
 !!$    call amg_daggrmat_biz_bld(a,desc_a,ilaggr,nlaggr, &
 !!$         & parms,ac,desc_ac,op_prol,op_restr,info)
 
   case(amg_min_energy_)
-    call amg_daggrmat_minnrg_bld(a,desc_a,ilaggr,nlaggr, &
-         & parms,ac,desc_ac,op_prol,op_restr,t_prol,info)
+    call amg_daggrmat_minnrg_bld(parms%aggr_prol,a,desc_a,&
+            ilaggr,nlaggr,parms,ac,desc_ac,op_prol,op_restr,&
+            t_prol,info)
 
   case default
     info = psb_err_internal_error_

amgprec/impl/aggregator/amg_d_parmatch_smth_bld.F90

@@ -69,6 +69,8 @@
 !
 !
 ! Arguments:
+!    dol1smoothing - Select between l1-Jacobi and Jacobi as smoother for the
+!                  tentative prolongator
 !    a          -  type(psb_dspmat_type), input.
 !                  The sparse matrix structure containing the local part of
 !                  the fine-level matrix.
@@ -102,8 +104,8 @@
 !    info       -  integer, output.
 !                  Error code.
 !
-subroutine amg_d_parmatch_smth_bld(ag,a,desc_a,ilaggr,nlaggr,parms,&
-     & ac,desc_ac,op_prol,op_restr,t_prol,info)
+subroutine amg_d_parmatch_smth_bld(dol1smoothing,ag,a,desc_a,ilaggr,nlaggr,&
+                parms,ac,desc_ac,op_prol,op_restr,t_prol,info)
   use psb_base_mod
   use amg_base_prec_type
   use amg_d_inner_mod
@@ -116,6 +118,7 @@ subroutine amg_d_parmatch_smth_bld(ag,a,desc_a,ilaggr,nlaggr,parms,&
   implicit none
 
   ! Arguments
+  integer(psb_ipk_), intent(in)            :: dol1smoothing
   class(amg_d_parmatch_aggregator_type), target, intent(inout) :: ag
   type(psb_dspmat_type), intent(in)      :: a
   type(psb_desc_type), intent(inout)     :: desc_a
@@ -137,7 +140,7 @@ subroutine amg_d_parmatch_smth_bld(ag,a,desc_a,ilaggr,nlaggr,parms,&
   type(psb_d_coo_sparse_mat)  :: coo_prol, coo_restr
   type(psb_d_csr_sparse_mat)  :: acsrf, csr_prol, acsr, tcsr
   real(psb_dpk_), allocatable :: adiag(:)
-  real(psb_dpk_), allocatable :: arwsum(:)
+  real(psb_dpk_), allocatable :: arwsum(:),l1rwsum(:)
   logical            :: filter_mat
   integer(psb_ipk_)            :: debug_level, debug_unit, err_act
   integer(psb_ipk_), parameter :: ncmax=16
@@ -145,6 +148,7 @@ subroutine amg_d_parmatch_smth_bld(ag,a,desc_a,ilaggr,nlaggr,parms,&
   logical, parameter :: debug_new=.false., dump_r=.false., dump_p=.false., debug=.false.
   character(len=80) :: filename
   logical, parameter :: do_timings=.false.
+  logical :: do_l1correction=.false.
   integer(psb_ipk_), save :: idx_spspmm=-1, idx_phase1=-1, idx_gtrans=-1, idx_phase2=-1, idx_refine=-1, idx_phase3=-1
   integer(psb_ipk_), save :: idx_cdasb=-1, idx_ptap=-1
 
@@ -166,6 +170,10 @@ subroutine amg_d_parmatch_smth_bld(ag,a,desc_a,ilaggr,nlaggr,parms,&
   ncol  = desc_a%get_local_cols()
 
   theta = parms%aggr_thresh
+  ! Check if we have to perform l1-Jacobi or Jacobi as smoother
+  if(dol1smoothing.eq.amg_l1_smooth_prol_) do_l1correction=.true.
+
+
   !write(0,*) me,' ',trim(name),' Start ',idx_spspmm
   if ((do_timings).and.(idx_spspmm==-1)) &
        & idx_spspmm = psb_get_timer_idx("PMC_SMTH_BLD: par_spspmm")
@@ -217,6 +225,19 @@ subroutine amg_d_parmatch_smth_bld(ag,a,desc_a,ilaggr,nlaggr,parms,&
   if (info == psb_success_) &
        & call psb_halo(adiag,desc_a,info)
   if (info == psb_success_) call a%cp_to(acsr)
+  ! Get the l1-diagonal of D
+  if (do_l1correction) then
+    allocate(l1rwsum(nrow))
+    call acsr%arwsum(l1rwsum)
+    if (info == psb_success_) &
+      & call psb_realloc(ncol,l1rwsum,info)
+    if (info == psb_success_) &
+      & call psb_halo(l1rwsum,desc_a,info)
+    ! \tilde{D}_{i,i} = \sum_{j \ne i} |a_{i,j}|
+    do i=1,size(adiag)
+      adiag(i) = adiag(i) + l1rwsum(i) - abs(adiag(i))
+    end do
+  end if
 
   if(info /= psb_success_) then
     call psb_errpush(psb_err_from_subroutine_,name,a_err='sp_getdiag')
@@ -246,7 +267,7 @@ subroutine amg_d_parmatch_smth_bld(ag,a,desc_a,ilaggr,nlaggr,parms,&
 
       enddo
       if (jd == -1) then
-        write(0,*) 'Wrong input: we need the diagonal!!!!', i
+        write(0,*) name,': Warning: there is no diagonal element', i
       else
         acsrf%val(jd)=acsrf%val(jd)-tmp
       end if
@@ -267,7 +288,10 @@ subroutine amg_d_parmatch_smth_bld(ag,a,desc_a,ilaggr,nlaggr,parms,&
 
   if (parms%aggr_omega_alg == amg_eig_est_) then
 
-    if (parms%aggr_eig == amg_max_norm_) then
+    if (do_l1correction) then
+      ! For l1-Jacobi this can be estimated with 1
+      parms%aggr_omega_val = done
+    else if (parms%aggr_eig == amg_max_norm_) then
       allocate(arwsum(nrow))
       call acsr%arwsum(arwsum)
       anorm = maxval(abs(adiag(1:nrow)*arwsum(1:nrow)))
@@ -373,6 +397,7 @@ subroutine amg_d_parmatch_smth_bld(ag,a,desc_a,ilaggr,nlaggr,parms,&
 
     end block
   end if
+  if (allocated(l1rwsum)) deallocate(l1rwsum)
   if (do_timings) call psb_toc(idx_phase2)
 
   if (debug_level >= psb_debug_outer_) &

amgprec/impl/aggregator/amg_d_parmatch_unsmth_bld.F90

@@ -68,6 +68,8 @@
 !
 !
 ! Arguments:
+!    dol1smoothing - this not actually used inside unsmoothed aggregation, it
+!                  is used just to perform a check
 !    a          -  type(psb_dspmat_type), input.
 !                  The sparse matrix structure containing the local part of
 !                  the fine-level matrix.
@@ -101,8 +103,8 @@
 !    info       -  integer, output.
 !                  Error code.
 !
-subroutine amg_d_parmatch_unsmth_bld(ag,a,desc_a,ilaggr,nlaggr,parms,&
-     & ac,desc_ac,op_prol,op_restr,t_prol,info)
+subroutine amg_d_parmatch_unsmth_bld(dol1smoothing,ag,a,desc_a,ilaggr,nlaggr,&
+                parms,ac,desc_ac,op_prol,op_restr,t_prol,info)
   use psb_base_mod
   use amg_base_prec_type
   use amg_d_inner_mod
@@ -115,6 +117,7 @@ subroutine amg_d_parmatch_unsmth_bld(ag,a,desc_a,ilaggr,nlaggr,parms,&
   implicit none
 
   ! Arguments
+  integer(psb_ipk_), intent(in) :: dol1smoothing
   class(amg_d_parmatch_aggregator_type), target, intent(inout) :: ag
   type(psb_dspmat_type), intent(in)      :: a
   type(psb_desc_type), intent(inout)     :: desc_a
@@ -159,6 +162,11 @@ subroutine amg_d_parmatch_unsmth_bld(ag,a,desc_a,ilaggr,nlaggr,parms,&
   ictxt = desc_a%get_context()
 
   call psb_info(ictxt, me, np)
+  if (dol1smoothing.ne.amg_no_smooth_) then
+    info=psb_err_fatal_;         
+    call psb_errpush(info,name,a_err='Are you trying to smooth an unsmoothed aggregation?')
+    goto 9999
+  end if
 
 #if !defined(SERIAL_MPI)
   nglob = desc_a%get_global_rows()

amgprec/impl/aggregator/amg_d_soc1_map_bld.F90

@@ -250,7 +250,7 @@ subroutine amg_d_soc1_map_bld(iorder,theta,clean_zeros,a,desc_a,nlaggr,ilaggr,in
               ! we will not reset. 
               if (j>nr) cycle step1
               if (ilaggr(j) > 0) cycle step1
-              if (abs(val(k)) > theta*sqrt(abs(diag(i)*diag(j)))) then
+              if ((abs(val(k)) > theta*sqrt(abs(diag(i)*diag(j)))).and.(diag(i).ne.dzero)) then
                 ip = ip + 1
                 icol(ip) = icol(k)
               end if
@@ -357,7 +357,7 @@ subroutine amg_d_soc1_map_bld(iorder,theta,clean_zeros,a,desc_a,nlaggr,ilaggr,in
       do k=1, nz
         j   = icol(k)
         if ((1<=j).and.(j<=nr)) then 
-          if (abs(val(k)) > theta*sqrt(abs(diag(i)*diag(j)))) then
+          if ((abs(val(k)) > theta*sqrt(abs(diag(i)*diag(j)))).and.(diag(i).ne.dzero)) then
             ip = ip + 1
             icol(ip) = icol(k)
           end if
@@ -545,4 +545,3 @@ subroutine amg_d_soc1_map_bld(iorder,theta,clean_zeros,a,desc_a,nlaggr,ilaggr,in
   return
 
 end subroutine amg_d_soc1_map_bld
-

amgprec/impl/aggregator/amg_daggrmat_minnrg_bld.f90

@@ -69,6 +69,7 @@
 !
 !
 ! Arguments:
+!    dol1smoothing - fictitious integer argument, it is not used inside
 !    a          -  type(psb_dspmat_type), input.     
 !                  The sparse matrix structure containing the local part of
 !                  the fine-level matrix.
@@ -104,8 +105,8 @@
 !                  Error code.
 !
 !
-subroutine amg_daggrmat_minnrg_bld(a,desc_a,ilaggr,nlaggr,parms,&
-     & ac,desc_ac,op_prol,op_restr,t_prol,info)
+subroutine amg_daggrmat_minnrg_bld(dol1smoothing,a,desc_a,ilaggr,nlaggr,&
+                parms,ac,desc_ac,op_prol,op_restr,t_prol,info)
   use psb_base_mod
   use amg_base_prec_type
   use amg_d_inner_mod, amg_protect_name => amg_daggrmat_minnrg_bld
@@ -113,6 +114,7 @@ subroutine amg_daggrmat_minnrg_bld(a,desc_a,ilaggr,nlaggr,parms,&
   implicit none
 
   ! Arguments
+  integer(psb_ipk_), intent(in)                 :: dol1smoothing
   type(psb_dspmat_type), intent(in)           :: a
   type(psb_desc_type), intent(inout)            :: desc_a
   integer(psb_lpk_), intent(inout)              :: ilaggr(:), nlaggr(:)
@@ -171,6 +173,13 @@ subroutine amg_daggrmat_minnrg_bld(a,desc_a,ilaggr,nlaggr,parms,&
 
   filter_mat = (parms%aggr_filter == amg_filter_mat_)
 
+  if (dol1smoothing.ne.amg_no_smooth_) then
+    info=psb_err_fatal_; 
+    call psb_errpush(info,name,a_err='Are you trying to smooth an unsmoothed aggregation?')
+    goto 9999
+  end if
+
+
   !NEEDS TO BE REWORKED !!
   
   ! naggr: number of local aggregates
@@ -300,7 +309,7 @@ subroutine amg_daggrmat_minnrg_bld(a,desc_a,ilaggr,nlaggr,parms,&
 !!$        endif
 !!$      enddo
 !!$      if (jd == -1) then 
-!!$        write(0,*) 'Wrong input: we need the diagonal!!!!', i
+!!$                write(0,*) name,': Warning: there is no diagonal element', i
 !!$      else
 !!$        acsrf%val(jd)=acsrf%val(jd)-tmp
 !!$      end if

amgprec/impl/aggregator/amg_daggrmat_nosmth_bld.f90

@@ -94,10 +94,11 @@
 !               
 !    info       -  integer, output.
 !                  Error code.
+!    dol1smoothing - optional, this is here just for interfacing reasons. It is not used by the
+!                    code
 !
-!
-subroutine amg_daggrmat_nosmth_bld(a,desc_a,ilaggr,nlaggr,parms,&
-     & ac,desc_ac,op_prol,op_restr,t_prol,info)
+subroutine amg_daggrmat_nosmth_bld(dol1smoothing,a,desc_a,ilaggr,nlaggr,&
+                parms,ac,desc_ac,op_prol,op_restr,t_prol,info)
   use psb_base_mod
   use amg_base_prec_type
   use amg_d_inner_mod, amg_protect_name => amg_daggrmat_nosmth_bld
@@ -105,6 +106,7 @@ subroutine amg_daggrmat_nosmth_bld(a,desc_a,ilaggr,nlaggr,parms,&
   implicit none
 
   ! Arguments
+  integer(psb_ipk_), intent(in)              :: dol1smoothing 
   type(psb_dspmat_type), intent(in)        :: a
   type(psb_desc_type), intent(inout)         :: desc_a
   integer(psb_lpk_), intent(inout)           :: ilaggr(:), nlaggr(:)
@@ -137,6 +139,12 @@ subroutine amg_daggrmat_nosmth_bld(a,desc_a,ilaggr,nlaggr,parms,&
 
   ctxt = desc_a%get_context()
   call psb_info(ctxt, me, np)
+  if (dol1smoothing.ne.amg_no_smooth_) then
+    info=psb_err_fatal_;         
+    call psb_errpush(info,name,a_err='Are you trying to smooth an unsmoothed aggregation?')
+    goto 9999
+  end if
+
   nglob = desc_a%get_global_rows()
   nrow  = desc_a%get_local_rows()
   ncol  = desc_a%get_local_cols()

amgprec/impl/aggregator/amg_daggrmat_smth_bld.f90

@@ -69,6 +69,8 @@
 !
 !
 ! Arguments:
+!    dol1smooth -  Integer taking the type of smoother that has to be used
+!                  on the tentative prolongator
 !    a          -  type(psb_dspmat_type), input.     
 !                  The sparse matrix structure containing the local part of
 !                  the fine-level matrix.
@@ -102,8 +104,8 @@
 !    info       -  integer, output.
 !                  Error code.
 !
-subroutine amg_daggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,&
-     & ac,desc_ac,op_prol,op_restr,t_prol,info)
+subroutine amg_daggrmat_smth_bld(dol1smoothing,a,desc_a,ilaggr,nlaggr,&
+                parms,ac,desc_ac,op_prol,op_restr,t_prol,info)
   use psb_base_mod
   use amg_base_prec_type
   use amg_d_inner_mod, amg_protect_name => amg_daggrmat_smth_bld
@@ -112,6 +114,7 @@ subroutine amg_daggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,&
   implicit none
 
   ! Arguments
+  integer(psb_ipk_), intent(in)            :: dol1smoothing
   type(psb_dspmat_type), intent(in)      :: a
   type(psb_desc_type), intent(inout)       :: desc_a
   integer(psb_lpk_), intent(inout)         :: ilaggr(:), nlaggr(:)
@@ -132,7 +135,7 @@ subroutine amg_daggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,&
   type(psb_d_coo_sparse_mat) :: coo_prol, coo_restr
   type(psb_d_csr_sparse_mat) :: acsr1,  acsrf, csr_prol, acsr
   real(psb_dpk_), allocatable :: adiag(:)
-  real(psb_dpk_), allocatable :: arwsum(:)
+  real(psb_dpk_), allocatable :: arwsum(:),l1rwsum(:)
   integer(psb_ipk_)  :: ierr(5)
   logical            :: filter_mat
   integer(psb_ipk_)            :: debug_level, debug_unit, err_act
@@ -141,6 +144,7 @@ subroutine amg_daggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,&
   logical, parameter :: debug_new=.false.
   character(len=80) :: filename
   logical, parameter :: do_timings=.false.
+  logical :: do_l1correction=.false.
   integer(psb_ipk_), save :: idx_spspmm=-1, idx_phase1=-1, idx_gtrans=-1, idx_phase2=-1, idx_refine=-1
   integer(psb_ipk_), save :: idx_phase3=-1, idx_cdasb=-1, idx_ptap=-1
 
@@ -173,6 +177,9 @@ subroutine amg_daggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,&
   if ((do_timings).and.(idx_ptap==-1)) &
        & idx_ptap = psb_get_timer_idx("DEC_SMTH_BLD: ptap_bld  ")
 
+  ! check if we have to use Jacobi or l1-Jacobi to smooth the tentative prolongator
+  if (dol1smoothing.eq.amg_l1_smooth_prol_) do_l1correction=.true.
+
 
   nglob = desc_a%get_global_rows()
   nrow  = desc_a%get_local_rows()
@@ -200,6 +207,24 @@ subroutine amg_daggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,&
   if (info == psb_success_) &
        & call psb_halo(adiag,desc_a,info)
   if (info == psb_success_) call a%cp_to(acsr)
+  !
+  ! Do the l1-correction on the diagonal if it is requested
+  !
+  if (do_l1correction) then
+    allocate(l1rwsum(nrow))
+    call acsr%arwsum(l1rwsum)
+    if (info == psb_success_) &
+      & call psb_realloc(ncol,l1rwsum,info)
+    if (info == psb_success_) &
+      & call psb_halo(l1rwsum,desc_a,info)
+    ! \tilde{D}_{i,i} = \sum_{j \ne i} |a_{i,j}|
+    !$OMP parallel do private(i) schedule(static)
+    do i=1,size(adiag)
+      adiag(i) = adiag(i) + l1rwsum(i) - abs(adiag(i))
+    end do
+    !$OMP end parallel do
+  end if
+
 
   if(info /= psb_success_) then
     call psb_errpush(psb_err_from_subroutine_,name,a_err='sp_getdiag')
@@ -230,6 +255,7 @@ subroutine amg_daggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,&
 
       enddo
       if (jd == -1) then 
+        ! if (.not.do_l1correction) 
         write(0,*) 'Wrong input: we need the diagonal!!!!', i
       else
         acsrf%val(jd)=acsrf%val(jd)-tmp
@@ -252,6 +278,10 @@ subroutine amg_daggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,&
   !$OMP end parallel do 
   if (parms%aggr_omega_alg == amg_eig_est_) then 
 
+    !if (do_l1correction) then
+    !  ! For l1-Jacobi this can be estimated with 1
+    !  parms%aggr_omega_val = done
+    !
     if (parms%aggr_eig == amg_max_norm_) then 
       allocate(arwsum(nrow))
       call acsr%arwsum(arwsum)      
@@ -259,7 +289,6 @@ subroutine amg_daggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,&
       call psb_amx(ctxt,anorm)
       omega = 4.d0/(3.d0*anorm)
       parms%aggr_omega_val = omega 
-
     else 
       info = psb_err_internal_error_
       call psb_errpush(info,name,a_err='invalid amg_aggr_eig_')
@@ -322,6 +351,7 @@ subroutine amg_daggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,&
   if (debug_level >= psb_debug_outer_) &
        & write(debug_unit,*) me,' ',trim(name),&
        & 'Done smooth_aggregate '
+  if (allocated(l1rwsum)) deallocate(l1rwsum)
   call psb_erractionrestore(err_act)
   return
 

amgprec/impl/aggregator/amg_s_dec_aggregator_mat_bld.f90

@@ -177,23 +177,24 @@ subroutine  amg_s_dec_aggregator_mat_bld(ag,parms,a,desc_a,ilaggr,nlaggr,&
   select case (parms%aggr_prol)
   case (amg_no_smooth_) 
 
-    call amg_saggrmat_nosmth_bld(a,desc_a,ilaggr,nlaggr,&
-         & parms,ac,desc_ac,op_prol,op_restr,t_prol,info)
+    call amg_saggrmat_nosmth_bld(parms%aggr_prol,a,desc_a,ilaggr,&
+            nlaggr,parms,ac,desc_ac,op_prol,op_restr,t_prol,info)
 
-  case(amg_smooth_prol_) 
+  case(amg_smooth_prol_,amg_l1_smooth_prol_) 
 
-    call amg_saggrmat_smth_bld(a,desc_a,ilaggr,nlaggr, &
-         & parms,ac,desc_ac,op_prol,op_restr,t_prol,info)
+    call amg_saggrmat_smth_bld(parms%aggr_prol,a,desc_a,&
+            ilaggr,nlaggr,parms,ac,desc_ac,op_prol,&
+            op_restr,t_prol,info)
 
 !!$  case(amg_biz_prol_) 
 !!$
 !!$    call amg_saggrmat_biz_bld(a,desc_a,ilaggr,nlaggr, &
 !!$         & parms,ac,desc_ac,op_prol,op_restr,t_prol,info)
-
+  
   case(amg_min_energy_) 
 
-    call amg_saggrmat_minnrg_bld(a,desc_a,ilaggr,nlaggr, &
-         & parms,ac,desc_ac,op_prol,op_restr,t_prol,info)
+    call amg_saggrmat_minnrg_bld(parms%aggr_prol,a,desc_a,ilaggr,&
+            nlaggr,parms,ac,desc_ac,op_prol,op_restr,t_prol,info)
 
   case default
     info = psb_err_internal_error_

amgprec/impl/aggregator/amg_s_parmatch_aggregator_mat_bld.F90

@@ -184,20 +184,23 @@ subroutine  amg_s_parmatch_aggregator_mat_bld(ag,parms,a,desc_a,ilaggr,nlaggr,&
   !
   select case (parms%aggr_prol)
   case (amg_no_smooth_)
-    call amg_s_parmatch_unsmth_bld(ag,a,desc_a,ilaggr,nlaggr,parms,&
-         & ac,desc_ac,op_prol,op_restr,t_prol,info)
+    call amg_s_parmatch_unsmth_bld(parms%aggr_prol,ag,a,desc_a,&
+            ilaggr,nlaggr,parms,ac,desc_ac,op_prol,op_restr,&
+            t_prol,info)
 
-  case(amg_smooth_prol_)
-    call amg_s_parmatch_smth_bld(ag,a,desc_a,ilaggr,nlaggr,parms,&
-         & ac,desc_ac,op_prol,op_restr,t_prol,info)
+  case(amg_smooth_prol_,amg_l1_smooth_prol_)
+    call amg_s_parmatch_smth_bld(parms%aggr_prol,ag,a,desc_a,&
+            ilaggr,nlaggr,parms,ac,desc_ac,op_prol,op_restr,&
+            t_prol,info)
 
 !!$  case(amg_biz_prol_)
 !!$    call amg_saggrmat_biz_bld(a,desc_a,ilaggr,nlaggr, &
 !!$         & parms,ac,desc_ac,op_prol,op_restr,info)
 
   case(amg_min_energy_)
-    call amg_saggrmat_minnrg_bld(a,desc_a,ilaggr,nlaggr, &
-         & parms,ac,desc_ac,op_prol,op_restr,t_prol,info)
+    call amg_saggrmat_minnrg_bld(parms%aggr_prol,a,desc_a,&
+            ilaggr,nlaggr,parms,ac,desc_ac,op_prol,op_restr,&
+            t_prol,info)
 
   case default
     info = psb_err_internal_error_

amgprec/impl/aggregator/amg_s_parmatch_smth_bld.F90

@@ -69,6 +69,8 @@
 !
 !
 ! Arguments:
+!    dol1smoothing - Select between l1-Jacobi and Jacobi as smoother for the
+!                  tentative prolongator
 !    a          -  type(psb_sspmat_type), input.
 !                  The sparse matrix structure containing the local part of
 !                  the fine-level matrix.
@@ -102,8 +104,8 @@
 !    info       -  integer, output.
 !                  Error code.
 !
-subroutine amg_s_parmatch_smth_bld(ag,a,desc_a,ilaggr,nlaggr,parms,&
-     & ac,desc_ac,op_prol,op_restr,t_prol,info)
+subroutine amg_s_parmatch_smth_bld(dol1smoothing,ag,a,desc_a,ilaggr,nlaggr,&
+                parms,ac,desc_ac,op_prol,op_restr,t_prol,info)
   use psb_base_mod
   use amg_base_prec_type
   use amg_s_inner_mod
@@ -116,6 +118,7 @@ subroutine amg_s_parmatch_smth_bld(ag,a,desc_a,ilaggr,nlaggr,parms,&
   implicit none
 
   ! Arguments
+  integer(psb_ipk_), intent(in)            :: dol1smoothing
   class(amg_s_parmatch_aggregator_type), target, intent(inout) :: ag
   type(psb_sspmat_type), intent(in)      :: a
   type(psb_desc_type), intent(inout)     :: desc_a
@@ -137,7 +140,7 @@ subroutine amg_s_parmatch_smth_bld(ag,a,desc_a,ilaggr,nlaggr,parms,&
   type(psb_s_coo_sparse_mat)  :: coo_prol, coo_restr
   type(psb_s_csr_sparse_mat)  :: acsrf, csr_prol, acsr, tcsr
   real(psb_spk_), allocatable :: adiag(:)
-  real(psb_spk_), allocatable :: arwsum(:)
+  real(psb_spk_), allocatable :: arwsum(:),l1rwsum(:)
   logical            :: filter_mat
   integer(psb_ipk_)            :: debug_level, debug_unit, err_act
   integer(psb_ipk_), parameter :: ncmax=16
@@ -145,6 +148,7 @@ subroutine amg_s_parmatch_smth_bld(ag,a,desc_a,ilaggr,nlaggr,parms,&
   logical, parameter :: debug_new=.false., dump_r=.false., dump_p=.false., debug=.false.
   character(len=80) :: filename
   logical, parameter :: do_timings=.false.
+  logical :: do_l1correction=.false.
   integer(psb_ipk_), save :: idx_spspmm=-1, idx_phase1=-1, idx_gtrans=-1, idx_phase2=-1, idx_refine=-1, idx_phase3=-1
   integer(psb_ipk_), save :: idx_cdasb=-1, idx_ptap=-1
 
@@ -166,6 +170,10 @@ subroutine amg_s_parmatch_smth_bld(ag,a,desc_a,ilaggr,nlaggr,parms,&
   ncol  = desc_a%get_local_cols()
 
   theta = parms%aggr_thresh
+  ! Check if we have to perform l1-Jacobi or Jacobi as smoother
+  if(dol1smoothing.eq.amg_l1_smooth_prol_) do_l1correction=.true.
+
+
   !write(0,*) me,' ',trim(name),' Start ',idx_spspmm
   if ((do_timings).and.(idx_spspmm==-1)) &
        & idx_spspmm = psb_get_timer_idx("PMC_SMTH_BLD: par_spspmm")
@@ -217,6 +225,19 @@ subroutine amg_s_parmatch_smth_bld(ag,a,desc_a,ilaggr,nlaggr,parms,&
   if (info == psb_success_) &
        & call psb_halo(adiag,desc_a,info)
   if (info == psb_success_) call a%cp_to(acsr)
+  ! Get the l1-diagonal of D
+  if (do_l1correction) then
+    allocate(l1rwsum(nrow))
+    call acsr%arwsum(l1rwsum)
+    if (info == psb_success_) &
+      & call psb_realloc(ncol,l1rwsum,info)
+    if (info == psb_success_) &
+      & call psb_halo(l1rwsum,desc_a,info)
+    ! \tilde{D}_{i,i} = \sum_{j \ne i} |a_{i,j}|
+    do i=1,size(adiag)
+      adiag(i) = adiag(i) + l1rwsum(i) - abs(adiag(i))
+    end do
+  end if
 
   if(info /= psb_success_) then
     call psb_errpush(psb_err_from_subroutine_,name,a_err='sp_getdiag')
@@ -246,7 +267,7 @@ subroutine amg_s_parmatch_smth_bld(ag,a,desc_a,ilaggr,nlaggr,parms,&
 
       enddo
       if (jd == -1) then
-        write(0,*) 'Wrong input: we need the diagonal!!!!', i
+        write(0,*) name,': Warning: there is no diagonal element', i
       else
         acsrf%val(jd)=acsrf%val(jd)-tmp
       end if
@@ -267,7 +288,10 @@ subroutine amg_s_parmatch_smth_bld(ag,a,desc_a,ilaggr,nlaggr,parms,&
 
   if (parms%aggr_omega_alg == amg_eig_est_) then
 
-    if (parms%aggr_eig == amg_max_norm_) then
+    if (do_l1correction) then
+      ! For l1-Jacobi this can be estimated with 1
+      parms%aggr_omega_val = done
+    else if (parms%aggr_eig == amg_max_norm_) then
       allocate(arwsum(nrow))
       call acsr%arwsum(arwsum)
       anorm = maxval(abs(adiag(1:nrow)*arwsum(1:nrow)))
@@ -373,6 +397,7 @@ subroutine amg_s_parmatch_smth_bld(ag,a,desc_a,ilaggr,nlaggr,parms,&
 
     end block
   end if
+  if (allocated(l1rwsum)) deallocate(l1rwsum)
   if (do_timings) call psb_toc(idx_phase2)
 
   if (debug_level >= psb_debug_outer_) &

amgprec/impl/aggregator/amg_s_parmatch_unsmth_bld.F90

@@ -68,6 +68,8 @@
 !
 !
 ! Arguments:
+!    dol1smoothing - this not actually used inside unsmoothed aggregation, it
+!                  is used just to perform a check
 !    a          -  type(psb_sspmat_type), input.
 !                  The sparse matrix structure containing the local part of
 !                  the fine-level matrix.
@@ -101,8 +103,8 @@
 !    info       -  integer, output.
 !                  Error code.
 !
-subroutine amg_s_parmatch_unsmth_bld(ag,a,desc_a,ilaggr,nlaggr,parms,&
-     & ac,desc_ac,op_prol,op_restr,t_prol,info)
+subroutine amg_s_parmatch_unsmth_bld(dol1smoothing,ag,a,desc_a,ilaggr,nlaggr,&
+                parms,ac,desc_ac,op_prol,op_restr,t_prol,info)
   use psb_base_mod
   use amg_base_prec_type
   use amg_s_inner_mod
@@ -115,6 +117,7 @@ subroutine amg_s_parmatch_unsmth_bld(ag,a,desc_a,ilaggr,nlaggr,parms,&
   implicit none
 
   ! Arguments
+  integer(psb_ipk_), intent(in) :: dol1smoothing
   class(amg_s_parmatch_aggregator_type), target, intent(inout) :: ag
   type(psb_sspmat_type), intent(in)      :: a
   type(psb_desc_type), intent(inout)     :: desc_a
@@ -159,6 +162,11 @@ subroutine amg_s_parmatch_unsmth_bld(ag,a,desc_a,ilaggr,nlaggr,parms,&
   ictxt = desc_a%get_context()
 
   call psb_info(ictxt, me, np)
+  if (dol1smoothing.ne.amg_no_smooth_) then
+    info=psb_err_fatal_;         
+    call psb_errpush(info,name,a_err='Are you trying to smooth an unsmoothed aggregation?')
+    goto 9999
+  end if
 
 #if !defined(SERIAL_MPI)
   nglob = desc_a%get_global_rows()

amgprec/impl/aggregator/amg_s_soc1_map_bld.F90

@@ -250,7 +250,7 @@ subroutine amg_s_soc1_map_bld(iorder,theta,clean_zeros,a,desc_a,nlaggr,ilaggr,in
               ! we will not reset. 
               if (j>nr) cycle step1
               if (ilaggr(j) > 0) cycle step1
-              if (abs(val(k)) > theta*sqrt(abs(diag(i)*diag(j)))) then
+              if ((abs(val(k)) > theta*sqrt(abs(diag(i)*diag(j)))).and.(diag(i).ne.szero)) then
                 ip = ip + 1
                 icol(ip) = icol(k)
               end if
@@ -357,7 +357,7 @@ subroutine amg_s_soc1_map_bld(iorder,theta,clean_zeros,a,desc_a,nlaggr,ilaggr,in
       do k=1, nz
         j   = icol(k)
         if ((1<=j).and.(j<=nr)) then 
-          if (abs(val(k)) > theta*sqrt(abs(diag(i)*diag(j)))) then
+          if ((abs(val(k)) > theta*sqrt(abs(diag(i)*diag(j)))).and.(diag(i).ne.szero)) then
             ip = ip + 1
             icol(ip) = icol(k)
           end if
@@ -545,4 +545,3 @@ subroutine amg_s_soc1_map_bld(iorder,theta,clean_zeros,a,desc_a,nlaggr,ilaggr,in
   return
 
 end subroutine amg_s_soc1_map_bld
-

amgprec/impl/aggregator/amg_saggrmat_minnrg_bld.f90

@@ -69,6 +69,7 @@
 !
 !
 ! Arguments:
+!    dol1smoothing - fictitious integer argument, it is not used inside
 !    a          -  type(psb_sspmat_type), input.     
 !                  The sparse matrix structure containing the local part of
 !                  the fine-level matrix.
@@ -104,8 +105,8 @@
 !                  Error code.
 !
 !
-subroutine amg_saggrmat_minnrg_bld(a,desc_a,ilaggr,nlaggr,parms,&
-     & ac,desc_ac,op_prol,op_restr,t_prol,info)
+subroutine amg_saggrmat_minnrg_bld(dol1smoothing,a,desc_a,ilaggr,nlaggr,&
+                parms,ac,desc_ac,op_prol,op_restr,t_prol,info)
   use psb_base_mod
   use amg_base_prec_type
   use amg_s_inner_mod, amg_protect_name => amg_saggrmat_minnrg_bld
@@ -113,6 +114,7 @@ subroutine amg_saggrmat_minnrg_bld(a,desc_a,ilaggr,nlaggr,parms,&
   implicit none
 
   ! Arguments
+  integer(psb_ipk_), intent(in)                 :: dol1smoothing
   type(psb_sspmat_type), intent(in)           :: a
   type(psb_desc_type), intent(inout)            :: desc_a
   integer(psb_lpk_), intent(inout)              :: ilaggr(:), nlaggr(:)
@@ -171,6 +173,13 @@ subroutine amg_saggrmat_minnrg_bld(a,desc_a,ilaggr,nlaggr,parms,&
 
   filter_mat = (parms%aggr_filter == amg_filter_mat_)
 
+  if (dol1smoothing.ne.amg_no_smooth_) then
+    info=psb_err_fatal_; 
+    call psb_errpush(info,name,a_err='Are you trying to smooth an unsmoothed aggregation?')
+    goto 9999
+  end if
+
+
   !NEEDS TO BE REWORKED !!
   
   ! naggr: number of local aggregates
@@ -300,7 +309,7 @@ subroutine amg_saggrmat_minnrg_bld(a,desc_a,ilaggr,nlaggr,parms,&
 !!$        endif
 !!$      enddo
 !!$      if (jd == -1) then 
-!!$        write(0,*) 'Wrong input: we need the diagonal!!!!', i
+!!$                write(0,*) name,': Warning: there is no diagonal element', i
 !!$      else
 !!$        acsrf%val(jd)=acsrf%val(jd)-tmp
 !!$      end if

amgprec/impl/aggregator/amg_saggrmat_nosmth_bld.f90

@@ -94,10 +94,11 @@
 !               
 !    info       -  integer, output.
 !                  Error code.
+!    dol1smoothing - optional, this is here just for interfacing reasons. It is not used by the
+!                    code
 !
-!
-subroutine amg_saggrmat_nosmth_bld(a,desc_a,ilaggr,nlaggr,parms,&
-     & ac,desc_ac,op_prol,op_restr,t_prol,info)
+subroutine amg_saggrmat_nosmth_bld(dol1smoothing,a,desc_a,ilaggr,nlaggr,&
+                parms,ac,desc_ac,op_prol,op_restr,t_prol,info)
   use psb_base_mod
   use amg_base_prec_type
   use amg_s_inner_mod, amg_protect_name => amg_saggrmat_nosmth_bld
@@ -105,6 +106,7 @@ subroutine amg_saggrmat_nosmth_bld(a,desc_a,ilaggr,nlaggr,parms,&
   implicit none
 
   ! Arguments
+  integer(psb_ipk_), intent(in)              :: dol1smoothing 
   type(psb_sspmat_type), intent(in)        :: a
   type(psb_desc_type), intent(inout)         :: desc_a
   integer(psb_lpk_), intent(inout)           :: ilaggr(:), nlaggr(:)
@@ -137,6 +139,12 @@ subroutine amg_saggrmat_nosmth_bld(a,desc_a,ilaggr,nlaggr,parms,&
 
   ctxt = desc_a%get_context()
   call psb_info(ctxt, me, np)
+  if (dol1smoothing.ne.amg_no_smooth_) then
+    info=psb_err_fatal_;         
+    call psb_errpush(info,name,a_err='Are you trying to smooth an unsmoothed aggregation?')
+    goto 9999
+  end if
+
   nglob = desc_a%get_global_rows()
   nrow  = desc_a%get_local_rows()
   ncol  = desc_a%get_local_cols()

amgprec/impl/aggregator/amg_saggrmat_smth_bld.f90

@@ -69,6 +69,8 @@
 !
 !
 ! Arguments:
+!    dol1smooth -  Integer taking the type of smoother that has to be used
+!                  on the tentative prolongator
 !    a          -  type(psb_sspmat_type), input.     
 !                  The sparse matrix structure containing the local part of
 !                  the fine-level matrix.
@@ -102,8 +104,8 @@
 !    info       -  integer, output.
 !                  Error code.
 !
-subroutine amg_saggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,&
-     & ac,desc_ac,op_prol,op_restr,t_prol,info)
+subroutine amg_saggrmat_smth_bld(dol1smoothing,a,desc_a,ilaggr,nlaggr,&
+                parms,ac,desc_ac,op_prol,op_restr,t_prol,info)
   use psb_base_mod
   use amg_base_prec_type
   use amg_s_inner_mod, amg_protect_name => amg_saggrmat_smth_bld
@@ -112,6 +114,7 @@ subroutine amg_saggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,&
   implicit none
 
   ! Arguments
+  integer(psb_ipk_), intent(in)            :: dol1smoothing
   type(psb_sspmat_type), intent(in)      :: a
   type(psb_desc_type), intent(inout)       :: desc_a
   integer(psb_lpk_), intent(inout)         :: ilaggr(:), nlaggr(:)
@@ -132,7 +135,7 @@ subroutine amg_saggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,&
   type(psb_s_coo_sparse_mat) :: coo_prol, coo_restr
   type(psb_s_csr_sparse_mat) :: acsr1,  acsrf, csr_prol, acsr
   real(psb_spk_), allocatable :: adiag(:)
-  real(psb_spk_), allocatable :: arwsum(:)
+  real(psb_spk_), allocatable :: arwsum(:),l1rwsum(:)
   integer(psb_ipk_)  :: ierr(5)
   logical            :: filter_mat
   integer(psb_ipk_)            :: debug_level, debug_unit, err_act
@@ -141,6 +144,7 @@ subroutine amg_saggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,&
   logical, parameter :: debug_new=.false.
   character(len=80) :: filename
   logical, parameter :: do_timings=.false.
+  logical :: do_l1correction=.false.
   integer(psb_ipk_), save :: idx_spspmm=-1, idx_phase1=-1, idx_gtrans=-1, idx_phase2=-1, idx_refine=-1
   integer(psb_ipk_), save :: idx_phase3=-1, idx_cdasb=-1, idx_ptap=-1
 
@@ -173,6 +177,9 @@ subroutine amg_saggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,&
   if ((do_timings).and.(idx_ptap==-1)) &
        & idx_ptap = psb_get_timer_idx("DEC_SMTH_BLD: ptap_bld  ")
 
+  ! check if we have to use Jacobi or l1-Jacobi to smooth the tentative prolongator
+  if (dol1smoothing.eq.amg_l1_smooth_prol_) do_l1correction=.true.
+
 
   nglob = desc_a%get_global_rows()
   nrow  = desc_a%get_local_rows()
@@ -200,6 +207,24 @@ subroutine amg_saggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,&
   if (info == psb_success_) &
        & call psb_halo(adiag,desc_a,info)
   if (info == psb_success_) call a%cp_to(acsr)
+  !
+  ! Do the l1-correction on the diagonal if it is requested
+  !
+  if (do_l1correction) then
+    allocate(l1rwsum(nrow))
+    call acsr%arwsum(l1rwsum)
+    if (info == psb_success_) &
+      & call psb_realloc(ncol,l1rwsum,info)
+    if (info == psb_success_) &
+      & call psb_halo(l1rwsum,desc_a,info)
+    ! \tilde{D}_{i,i} = \sum_{j \ne i} |a_{i,j}|
+    !$OMP parallel do private(i) schedule(static)
+    do i=1,size(adiag)
+      adiag(i) = adiag(i) + l1rwsum(i) - abs(adiag(i))
+    end do
+    !$OMP end parallel do
+  end if
+
 
   if(info /= psb_success_) then
     call psb_errpush(psb_err_from_subroutine_,name,a_err='sp_getdiag')
@@ -230,6 +255,7 @@ subroutine amg_saggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,&
 
       enddo
       if (jd == -1) then 
+        ! if (.not.do_l1correction) 
         write(0,*) 'Wrong input: we need the diagonal!!!!', i
       else
         acsrf%val(jd)=acsrf%val(jd)-tmp
@@ -252,6 +278,10 @@ subroutine amg_saggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,&
   !$OMP end parallel do 
   if (parms%aggr_omega_alg == amg_eig_est_) then 
 
+    !if (do_l1correction) then
+    !  ! For l1-Jacobi this can be estimated with 1
+    !  parms%aggr_omega_val = done
+    !
     if (parms%aggr_eig == amg_max_norm_) then 
       allocate(arwsum(nrow))
       call acsr%arwsum(arwsum)      
@@ -259,7 +289,6 @@ subroutine amg_saggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,&
       call psb_amx(ctxt,anorm)
       omega = 4.d0/(3.d0*anorm)
       parms%aggr_omega_val = omega 
-
     else 
       info = psb_err_internal_error_
       call psb_errpush(info,name,a_err='invalid amg_aggr_eig_')
@@ -322,6 +351,7 @@ subroutine amg_saggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,&
   if (debug_level >= psb_debug_outer_) &
        & write(debug_unit,*) me,' ',trim(name),&
        & 'Done smooth_aggregate '
+  if (allocated(l1rwsum)) deallocate(l1rwsum)
   call psb_erractionrestore(err_act)
   return
 

amgprec/impl/aggregator/amg_z_dec_aggregator_mat_bld.f90

@@ -177,23 +177,24 @@ subroutine  amg_z_dec_aggregator_mat_bld(ag,parms,a,desc_a,ilaggr,nlaggr,&
   select case (parms%aggr_prol)
   case (amg_no_smooth_) 
 
-    call amg_zaggrmat_nosmth_bld(a,desc_a,ilaggr,nlaggr,&
-         & parms,ac,desc_ac,op_prol,op_restr,t_prol,info)
+    call amg_zaggrmat_nosmth_bld(parms%aggr_prol,a,desc_a,ilaggr,&
+            nlaggr,parms,ac,desc_ac,op_prol,op_restr,t_prol,info)
 
-  case(amg_smooth_prol_) 
+  case(amg_smooth_prol_,amg_l1_smooth_prol_) 
 
-    call amg_zaggrmat_smth_bld(a,desc_a,ilaggr,nlaggr, &
-         & parms,ac,desc_ac,op_prol,op_restr,t_prol,info)
+    call amg_zaggrmat_smth_bld(parms%aggr_prol,a,desc_a,&
+            ilaggr,nlaggr,parms,ac,desc_ac,op_prol,&
+            op_restr,t_prol,info)
 
 !!$  case(amg_biz_prol_) 
 !!$
 !!$    call amg_zaggrmat_biz_bld(a,desc_a,ilaggr,nlaggr, &
 !!$         & parms,ac,desc_ac,op_prol,op_restr,t_prol,info)
-
+  
   case(amg_min_energy_) 
 
-    call amg_zaggrmat_minnrg_bld(a,desc_a,ilaggr,nlaggr, &
-         & parms,ac,desc_ac,op_prol,op_restr,t_prol,info)
+    call amg_zaggrmat_minnrg_bld(parms%aggr_prol,a,desc_a,ilaggr,&
+            nlaggr,parms,ac,desc_ac,op_prol,op_restr,t_prol,info)
 
   case default
     info = psb_err_internal_error_

amgprec/impl/aggregator/amg_z_soc1_map_bld.F90

@@ -250,7 +250,7 @@ subroutine amg_z_soc1_map_bld(iorder,theta,clean_zeros,a,desc_a,nlaggr,ilaggr,in
               ! we will not reset. 
               if (j>nr) cycle step1
               if (ilaggr(j) > 0) cycle step1
-              if (abs(val(k)) > theta*sqrt(abs(diag(i)*diag(j)))) then
+              if ((abs(val(k)) > theta*sqrt(abs(diag(i)*diag(j)))).and.(diag(i).ne.zzero)) then
                 ip = ip + 1
                 icol(ip) = icol(k)
               end if
@@ -357,7 +357,7 @@ subroutine amg_z_soc1_map_bld(iorder,theta,clean_zeros,a,desc_a,nlaggr,ilaggr,in
       do k=1, nz
         j   = icol(k)
         if ((1<=j).and.(j<=nr)) then 
-          if (abs(val(k)) > theta*sqrt(abs(diag(i)*diag(j)))) then
+          if ((abs(val(k)) > theta*sqrt(abs(diag(i)*diag(j)))).and.(diag(i).ne.zzero)) then
             ip = ip + 1
             icol(ip) = icol(k)
           end if
@@ -545,4 +545,3 @@ subroutine amg_z_soc1_map_bld(iorder,theta,clean_zeros,a,desc_a,nlaggr,ilaggr,in
   return
 
 end subroutine amg_z_soc1_map_bld
-

amgprec/impl/aggregator/amg_zaggrmat_minnrg_bld.f90

@@ -69,6 +69,7 @@
 !
 !
 ! Arguments:
+!    dol1smoothing - fictitious integer argument, it is not used inside
 !    a          -  type(psb_zspmat_type), input.     
 !                  The sparse matrix structure containing the local part of
 !                  the fine-level matrix.
@@ -104,8 +105,8 @@
 !                  Error code.
 !
 !
-subroutine amg_zaggrmat_minnrg_bld(a,desc_a,ilaggr,nlaggr,parms,&
-     & ac,desc_ac,op_prol,op_restr,t_prol,info)
+subroutine amg_zaggrmat_minnrg_bld(dol1smoothing,a,desc_a,ilaggr,nlaggr,&
+                parms,ac,desc_ac,op_prol,op_restr,t_prol,info)
   use psb_base_mod
   use amg_base_prec_type
   use amg_z_inner_mod, amg_protect_name => amg_zaggrmat_minnrg_bld
@@ -113,6 +114,7 @@ subroutine amg_zaggrmat_minnrg_bld(a,desc_a,ilaggr,nlaggr,parms,&
   implicit none
 
   ! Arguments
+  integer(psb_ipk_), intent(in)                 :: dol1smoothing
   type(psb_zspmat_type), intent(in)           :: a
   type(psb_desc_type), intent(inout)            :: desc_a
   integer(psb_lpk_), intent(inout)              :: ilaggr(:), nlaggr(:)
@@ -171,6 +173,13 @@ subroutine amg_zaggrmat_minnrg_bld(a,desc_a,ilaggr,nlaggr,parms,&
 
   filter_mat = (parms%aggr_filter == amg_filter_mat_)
 
+  if (dol1smoothing.ne.amg_no_smooth_) then
+    info=psb_err_fatal_; 
+    call psb_errpush(info,name,a_err='Are you trying to smooth an unsmoothed aggregation?')
+    goto 9999
+  end if
+
+
   !NEEDS TO BE REWORKED !!
   
   ! naggr: number of local aggregates
@@ -300,7 +309,7 @@ subroutine amg_zaggrmat_minnrg_bld(a,desc_a,ilaggr,nlaggr,parms,&
 !!$        endif
 !!$      enddo
 !!$      if (jd == -1) then 
-!!$        write(0,*) 'Wrong input: we need the diagonal!!!!', i
+!!$                write(0,*) name,': Warning: there is no diagonal element', i
 !!$      else
 !!$        acsrf%val(jd)=acsrf%val(jd)-tmp
 !!$      end if

amgprec/impl/aggregator/amg_zaggrmat_nosmth_bld.f90

@@ -94,10 +94,11 @@
 !               
 !    info       -  integer, output.
 !                  Error code.
+!    dol1smoothing - optional, this is here just for interfacing reasons. It is not used by the
+!                    code
 !
-!
-subroutine amg_zaggrmat_nosmth_bld(a,desc_a,ilaggr,nlaggr,parms,&
-     & ac,desc_ac,op_prol,op_restr,t_prol,info)
+subroutine amg_zaggrmat_nosmth_bld(dol1smoothing,a,desc_a,ilaggr,nlaggr,&
+                parms,ac,desc_ac,op_prol,op_restr,t_prol,info)
   use psb_base_mod
   use amg_base_prec_type
   use amg_z_inner_mod, amg_protect_name => amg_zaggrmat_nosmth_bld
@@ -105,6 +106,7 @@ subroutine amg_zaggrmat_nosmth_bld(a,desc_a,ilaggr,nlaggr,parms,&
   implicit none
 
   ! Arguments
+  integer(psb_ipk_), intent(in)              :: dol1smoothing 
   type(psb_zspmat_type), intent(in)        :: a
   type(psb_desc_type), intent(inout)         :: desc_a
   integer(psb_lpk_), intent(inout)           :: ilaggr(:), nlaggr(:)
@@ -137,6 +139,12 @@ subroutine amg_zaggrmat_nosmth_bld(a,desc_a,ilaggr,nlaggr,parms,&
 
   ctxt = desc_a%get_context()
   call psb_info(ctxt, me, np)
+  if (dol1smoothing.ne.amg_no_smooth_) then
+    info=psb_err_fatal_;         
+    call psb_errpush(info,name,a_err='Are you trying to smooth an unsmoothed aggregation?')
+    goto 9999
+  end if
+
   nglob = desc_a%get_global_rows()
   nrow  = desc_a%get_local_rows()
   ncol  = desc_a%get_local_cols()

amgprec/impl/aggregator/amg_zaggrmat_smth_bld.f90

@@ -69,6 +69,8 @@
 !
 !
 ! Arguments:
+!    dol1smooth -  Integer taking the type of smoother that has to be used
+!                  on the tentative prolongator
 !    a          -  type(psb_zspmat_type), input.     
 !                  The sparse matrix structure containing the local part of
 !                  the fine-level matrix.
@@ -102,8 +104,8 @@
 !    info       -  integer, output.
 !                  Error code.
 !
-subroutine amg_zaggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,&
-     & ac,desc_ac,op_prol,op_restr,t_prol,info)
+subroutine amg_zaggrmat_smth_bld(dol1smoothing,a,desc_a,ilaggr,nlaggr,&
+                parms,ac,desc_ac,op_prol,op_restr,t_prol,info)
   use psb_base_mod
   use amg_base_prec_type
   use amg_z_inner_mod, amg_protect_name => amg_zaggrmat_smth_bld
@@ -112,6 +114,7 @@ subroutine amg_zaggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,&
   implicit none
 
   ! Arguments
+  integer(psb_ipk_), intent(in)            :: dol1smoothing
   type(psb_zspmat_type), intent(in)      :: a
   type(psb_desc_type), intent(inout)       :: desc_a
   integer(psb_lpk_), intent(inout)         :: ilaggr(:), nlaggr(:)
@@ -132,7 +135,7 @@ subroutine amg_zaggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,&
   type(psb_z_coo_sparse_mat) :: coo_prol, coo_restr
   type(psb_z_csr_sparse_mat) :: acsr1,  acsrf, csr_prol, acsr
   complex(psb_dpk_), allocatable :: adiag(:)
-  real(psb_dpk_), allocatable :: arwsum(:)
+  real(psb_dpk_), allocatable :: arwsum(:),l1rwsum(:)
   integer(psb_ipk_)  :: ierr(5)
   logical            :: filter_mat
   integer(psb_ipk_)            :: debug_level, debug_unit, err_act
@@ -141,6 +144,7 @@ subroutine amg_zaggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,&
   logical, parameter :: debug_new=.false.
   character(len=80) :: filename
   logical, parameter :: do_timings=.false.
+  logical :: do_l1correction=.false.
   integer(psb_ipk_), save :: idx_spspmm=-1, idx_phase1=-1, idx_gtrans=-1, idx_phase2=-1, idx_refine=-1
   integer(psb_ipk_), save :: idx_phase3=-1, idx_cdasb=-1, idx_ptap=-1
 
@@ -173,6 +177,9 @@ subroutine amg_zaggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,&
   if ((do_timings).and.(idx_ptap==-1)) &
        & idx_ptap = psb_get_timer_idx("DEC_SMTH_BLD: ptap_bld  ")
 
+  ! check if we have to use Jacobi or l1-Jacobi to smooth the tentative prolongator
+  if (dol1smoothing.eq.amg_l1_smooth_prol_) do_l1correction=.true.
+
 
   nglob = desc_a%get_global_rows()
   nrow  = desc_a%get_local_rows()
@@ -200,6 +207,24 @@ subroutine amg_zaggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,&
   if (info == psb_success_) &
        & call psb_halo(adiag,desc_a,info)
   if (info == psb_success_) call a%cp_to(acsr)
+  !
+  ! Do the l1-correction on the diagonal if it is requested
+  !
+  if (do_l1correction) then
+    allocate(l1rwsum(nrow))
+    call acsr%arwsum(l1rwsum)
+    if (info == psb_success_) &
+      & call psb_realloc(ncol,l1rwsum,info)
+    if (info == psb_success_) &
+      & call psb_halo(l1rwsum,desc_a,info)
+    ! \tilde{D}_{i,i} = \sum_{j \ne i} |a_{i,j}|
+    !$OMP parallel do private(i) schedule(static)
+    do i=1,size(adiag)
+      adiag(i) = adiag(i) + l1rwsum(i) - abs(adiag(i))
+    end do
+    !$OMP end parallel do
+  end if
+
 
   if(info /= psb_success_) then
     call psb_errpush(psb_err_from_subroutine_,name,a_err='sp_getdiag')
@@ -230,6 +255,7 @@ subroutine amg_zaggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,&
 
       enddo
       if (jd == -1) then 
+        ! if (.not.do_l1correction) 
         write(0,*) 'Wrong input: we need the diagonal!!!!', i
       else
         acsrf%val(jd)=acsrf%val(jd)-tmp
@@ -252,6 +278,10 @@ subroutine amg_zaggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,&
   !$OMP end parallel do 
   if (parms%aggr_omega_alg == amg_eig_est_) then 
 
+    !if (do_l1correction) then
+    !  ! For l1-Jacobi this can be estimated with 1
+    !  parms%aggr_omega_val = done
+    !
     if (parms%aggr_eig == amg_max_norm_) then 
       allocate(arwsum(nrow))
       call acsr%arwsum(arwsum)      
@@ -259,7 +289,6 @@ subroutine amg_zaggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,&
       call psb_amx(ctxt,anorm)
       omega = 4.d0/(3.d0*anorm)
       parms%aggr_omega_val = omega 
-
     else 
       info = psb_err_internal_error_
       call psb_errpush(info,name,a_err='invalid amg_aggr_eig_')
@@ -322,6 +351,7 @@ subroutine amg_zaggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,&
   if (debug_level >= psb_debug_outer_) &
        & write(debug_unit,*) me,' ',trim(name),&
        & 'Done smooth_aggregate '
+  if (allocated(l1rwsum)) deallocate(l1rwsum)
   call psb_erractionrestore(err_act)
   return
 

amgprec/impl/aggregator/clean.cpp

@@ -1,5 +1,6 @@
 #include "MatchBoxPC.h"
 // TODO comment
+#if !defined(SERIAL_MPI)
 
 void clean(MilanLongInt NLVer,
            MilanInt myRank,
@@ -88,3 +89,4 @@ void clean(MilanLongInt NLVer,
         }
     }
 }
+#endif

amgprec/impl/aggregator/computeCandidateMate.cpp

@@ -8,6 +8,7 @@
  * @param edgeLocWeight
  * @return
  */
+
 MilanLongInt firstComputeCandidateMateD(MilanLongInt adj1,
                                        MilanLongInt adj2,
                                        MilanLongInt *verLocInd,
@@ -134,3 +135,4 @@ MilanLongInt computeCandidateMateS(MilanLongInt adj1,
 
     return w;
 }
+

amgprec/impl/aggregator/parallelComputeCandidateMateB.cpp

@@ -1,4 +1,5 @@
 #include "MatchBoxPC.h"
+
 void PARALLEL_COMPUTE_CANDIDATE_MATE_BD(MilanLongInt NLVer,
                                               MilanLongInt *verLocPtr,
                                               MilanLongInt *verLocInd,
@@ -52,3 +53,4 @@ void PARALLEL_COMPUTE_CANDIDATE_MATE_BS(MilanLongInt NLVer,
         }
     }
 }
+

amgprec/impl/aggregator/processExposedVertex.cpp

@@ -366,3 +366,4 @@ void PARALLEL_PROCESS_EXPOSED_VERTEX_BS(MilanLongInt NLVer,
 	
     } // End of parallel region
 }
+

amgprec/impl/aggregator/processMatchedVertices.cpp

@@ -582,3 +582,4 @@ void processMatchedVerticesS(
 #endif
     } // End of parallel region
 }
+

amgprec/impl/aggregator/processMatchedVerticesAndSendMessages.cpp

@@ -588,3 +588,4 @@ void processMatchedVerticesAndSendMessagesS(
     cout << myRank<<" Done sending messages"<<endl;
 #endif
 }
+

amgprec/impl/aggregator/processMessages.cpp

@@ -1,5 +1,6 @@
 #include "MatchBoxPC.h"
 //#define DEBUG_HANG_
+#if !defined(SERIAL_MPI)
 
 void processMessagesD(
     MilanLongInt NLVer,
@@ -629,3 +630,4 @@ void processMessagesS(
 
     return;
 }
+#endif

amgprec/impl/aggregator/queueTransfer.cpp

@@ -30,3 +30,4 @@ void queuesTransfer(vector<MilanLongInt> &U,
     privateQOwner.clear();
 
 }
+

amgprec/impl/amg_c_hierarchy_bld.F90

@@ -451,11 +451,16 @@ subroutine amg_c_hierarchy_bld(a,desc_a,prec,info)
     if (.not.associated(prec%precv(1)%base_desc,desc_a)) then
       prec%precv(1)%base_desc => prec%precv(1)%desc_ac
     end if
-    do i=2, iszv 
+    do i=2, iszv
       prec%precv(i)%base_a       => prec%precv(i)%ac
       prec%precv(i)%base_desc    => prec%precv(i)%desc_ac
-      prec%precv(i)%linmap%p_desc_U => prec%precv(i-1)%base_desc
-      prec%precv(i)%linmap%p_desc_V => prec%precv(i)%base_desc
+      ! This is needed when the linmap object has been built
+      ! reusing the base_desc descriptor through a pointer.
+      ! With PSBLAS 4 we will have a better solution
+      if (associated(prec%precv(i)%linmap%p_desc_U)) &
+           & prec%precv(i)%linmap%p_desc_U => prec%precv(i-1)%base_desc
+      if (associated(prec%precv(i)%linmap%p_desc_V))&
+           & prec%precv(i)%linmap%p_desc_V => prec%precv(i)%base_desc
     end do
   end if
 

amgprec/impl/amg_d_hierarchy_bld.F90

@@ -451,11 +451,16 @@ subroutine amg_d_hierarchy_bld(a,desc_a,prec,info)
     if (.not.associated(prec%precv(1)%base_desc,desc_a)) then
       prec%precv(1)%base_desc => prec%precv(1)%desc_ac
     end if
-    do i=2, iszv 
+    do i=2, iszv
       prec%precv(i)%base_a       => prec%precv(i)%ac
       prec%precv(i)%base_desc    => prec%precv(i)%desc_ac
-      prec%precv(i)%linmap%p_desc_U => prec%precv(i-1)%base_desc
-      prec%precv(i)%linmap%p_desc_V => prec%precv(i)%base_desc
+      ! This is needed when the linmap object has been built
+      ! reusing the base_desc descriptor through a pointer.
+      ! With PSBLAS 4 we will have a better solution
+      if (associated(prec%precv(i)%linmap%p_desc_U)) &
+           & prec%precv(i)%linmap%p_desc_U => prec%precv(i-1)%base_desc
+      if (associated(prec%precv(i)%linmap%p_desc_V))&
+           & prec%precv(i)%linmap%p_desc_V => prec%precv(i)%base_desc
     end do
   end if
 

amgprec/impl/amg_s_hierarchy_bld.F90

@@ -451,11 +451,16 @@ subroutine amg_s_hierarchy_bld(a,desc_a,prec,info)
     if (.not.associated(prec%precv(1)%base_desc,desc_a)) then
       prec%precv(1)%base_desc => prec%precv(1)%desc_ac
     end if
-    do i=2, iszv 
+    do i=2, iszv
       prec%precv(i)%base_a       => prec%precv(i)%ac
       prec%precv(i)%base_desc    => prec%precv(i)%desc_ac
-      prec%precv(i)%linmap%p_desc_U => prec%precv(i-1)%base_desc
-      prec%precv(i)%linmap%p_desc_V => prec%precv(i)%base_desc
+      ! This is needed when the linmap object has been built
+      ! reusing the base_desc descriptor through a pointer.
+      ! With PSBLAS 4 we will have a better solution
+      if (associated(prec%precv(i)%linmap%p_desc_U)) &
+           & prec%precv(i)%linmap%p_desc_U => prec%precv(i-1)%base_desc
+      if (associated(prec%precv(i)%linmap%p_desc_V))&
+           & prec%precv(i)%linmap%p_desc_V => prec%precv(i)%base_desc
     end do
   end if
 

amgprec/impl/amg_z_hierarchy_bld.F90

@@ -451,11 +451,16 @@ subroutine amg_z_hierarchy_bld(a,desc_a,prec,info)
     if (.not.associated(prec%precv(1)%base_desc,desc_a)) then
       prec%precv(1)%base_desc => prec%precv(1)%desc_ac
     end if
-    do i=2, iszv 
+    do i=2, iszv
       prec%precv(i)%base_a       => prec%precv(i)%ac
       prec%precv(i)%base_desc    => prec%precv(i)%desc_ac
-      prec%precv(i)%linmap%p_desc_U => prec%precv(i-1)%base_desc
-      prec%precv(i)%linmap%p_desc_V => prec%precv(i)%base_desc
+      ! This is needed when the linmap object has been built
+      ! reusing the base_desc descriptor through a pointer.
+      ! With PSBLAS 4 we will have a better solution
+      if (associated(prec%precv(i)%linmap%p_desc_U)) &
+           & prec%precv(i)%linmap%p_desc_U => prec%precv(i-1)%base_desc
+      if (associated(prec%precv(i)%linmap%p_desc_V))&
+           & prec%precv(i)%linmap%p_desc_V => prec%precv(i)%base_desc
     end do
   end if
 

amgprec/impl/level/amg_c_base_onelev_csetc.F90

@@ -241,6 +241,8 @@ subroutine amg_c_base_onelev_csetc(lv,what,val,info,pos,idx)
       if (info == 0) deallocate(lv%aggr,stat=info)
       if (info /= 0) then
         info = psb_err_internal_error_
+        call psb_errpush(info,name,a_err='aggregator deallocation?')
+        goto 9999
         return
       end if
     end if
@@ -250,8 +252,16 @@ subroutine amg_c_base_onelev_csetc(lv,what,val,info,pos,idx)
       allocate(amg_c_dec_aggregator_type :: lv%aggr, stat=info)
     case('SYMDEC')
       allocate(amg_c_symdec_aggregator_type :: lv%aggr, stat=info)
-    case default
+#if !defined(SERIAL_MPI)
+#endif
+  case default
       info =  psb_err_internal_error_
+#if !defined(SERIAL_MPI)
+      call psb_errpush(info,name,a_err='Unsupported PAR_AGGR_ALG')
+#else
+      call psb_errpush(info,name,a_err='PAR_AGGR_ALG unsupported (SERIAL_MPI on)')
+#endif
+      goto 9999
     end select
     if (info == psb_success_) call lv%aggr%default()
 

amgprec/impl/level/amg_c_base_onelev_dump.f90

@@ -127,8 +127,7 @@ subroutine amg_c_base_onelev_dump(lv,level,info,prefix,head,ac,rp,&
         ivr = lv%linmap%p_desc_U%get_global_indices(owned=.false.)
         write(fname(lname+1:),'(a,i3.3,a)')'_l',level,'_tprol.mtx'
         !
-        !  This is not implemented yet.
-        !call lv%tprol%print(fname,head=head,ivr=ivr)
+        call lv%tprol%print(fname,head=head,ivr=ivr)
       end if
     end if
   else
@@ -151,8 +150,7 @@ subroutine amg_c_base_onelev_dump(lv,level,info,prefix,head,ac,rp,&
       if (tprol_) then 
         write(fname(lname+1:),'(a,i3.3,a)')'_l',level,'_tprol.mtx'
         !
-        !  This is not implemented yet.
-        !call lv%tprol%print(fname,head=head)
+        call lv%tprol%print(fname,head=head)
       end if
     end if
   end if

amgprec/impl/level/amg_d_base_onelev_csetc.F90

@@ -42,7 +42,9 @@ subroutine amg_d_base_onelev_csetc(lv,what,val,info,pos,idx)
   use amg_d_base_aggregator_mod
   use amg_d_dec_aggregator_mod
   use amg_d_symdec_aggregator_mod
+#if !defined(SERIAL_MPI)
   use amg_d_parmatch_aggregator_mod
+#endif
   use amg_d_poly_smoother
   use amg_d_jac_smoother
   use amg_d_as_smoother
@@ -267,6 +269,8 @@ subroutine amg_d_base_onelev_csetc(lv,what,val,info,pos,idx)
       if (info == 0) deallocate(lv%aggr,stat=info)
       if (info /= 0) then
         info = psb_err_internal_error_
+        call psb_errpush(info,name,a_err='aggregator deallocation?')
+        goto 9999
         return
       end if
     end if
@@ -276,10 +280,18 @@ subroutine amg_d_base_onelev_csetc(lv,what,val,info,pos,idx)
       allocate(amg_d_dec_aggregator_type :: lv%aggr, stat=info)
     case('SYMDEC')
       allocate(amg_d_symdec_aggregator_type :: lv%aggr, stat=info)
+#if !defined(SERIAL_MPI)
     case('COUP','COUPLED')
       allocate(amg_d_parmatch_aggregator_type :: lv%aggr, stat=info)
-    case default
+#endif
+  case default
       info =  psb_err_internal_error_
+#if !defined(SERIAL_MPI)
+      call psb_errpush(info,name,a_err='Unsupported PAR_AGGR_ALG')
+#else
+      call psb_errpush(info,name,a_err='PAR_AGGR_ALG unsupported (SERIAL_MPI on)')
+#endif
+      goto 9999
     end select
     if (info == psb_success_) call lv%aggr%default()
 

amgprec/impl/level/amg_d_base_onelev_dump.f90

@@ -127,8 +127,7 @@ subroutine amg_d_base_onelev_dump(lv,level,info,prefix,head,ac,rp,&
         ivr = lv%linmap%p_desc_U%get_global_indices(owned=.false.)
         write(fname(lname+1:),'(a,i3.3,a)')'_l',level,'_tprol.mtx'
         !
-        !  This is not implemented yet.
-        !call lv%tprol%print(fname,head=head,ivr=ivr)
+        call lv%tprol%print(fname,head=head,ivr=ivr)
       end if
     end if
   else
@@ -151,8 +150,7 @@ subroutine amg_d_base_onelev_dump(lv,level,info,prefix,head,ac,rp,&
       if (tprol_) then 
         write(fname(lname+1:),'(a,i3.3,a)')'_l',level,'_tprol.mtx'
         !
-        !  This is not implemented yet.
-        !call lv%tprol%print(fname,head=head)
+        call lv%tprol%print(fname,head=head)
       end if
     end if
   end if

amgprec/impl/level/amg_s_base_onelev_csetc.F90

@@ -42,7 +42,9 @@ subroutine amg_s_base_onelev_csetc(lv,what,val,info,pos,idx)
   use amg_s_base_aggregator_mod
   use amg_s_dec_aggregator_mod
   use amg_s_symdec_aggregator_mod
+#if !defined(SERIAL_MPI)
   use amg_s_parmatch_aggregator_mod
+#endif
   use amg_s_poly_smoother
   use amg_s_jac_smoother
   use amg_s_as_smoother
@@ -247,6 +249,8 @@ subroutine amg_s_base_onelev_csetc(lv,what,val,info,pos,idx)
       if (info == 0) deallocate(lv%aggr,stat=info)
       if (info /= 0) then
         info = psb_err_internal_error_
+        call psb_errpush(info,name,a_err='aggregator deallocation?')
+        goto 9999
         return
       end if
     end if
@@ -256,10 +260,18 @@ subroutine amg_s_base_onelev_csetc(lv,what,val,info,pos,idx)
       allocate(amg_s_dec_aggregator_type :: lv%aggr, stat=info)
     case('SYMDEC')
       allocate(amg_s_symdec_aggregator_type :: lv%aggr, stat=info)
+#if !defined(SERIAL_MPI)
     case('COUP','COUPLED')
       allocate(amg_s_parmatch_aggregator_type :: lv%aggr, stat=info)
-    case default
+#endif
+  case default
       info =  psb_err_internal_error_
+#if !defined(SERIAL_MPI)
+      call psb_errpush(info,name,a_err='Unsupported PAR_AGGR_ALG')
+#else
+      call psb_errpush(info,name,a_err='PAR_AGGR_ALG unsupported (SERIAL_MPI on)')
+#endif
+      goto 9999
     end select
     if (info == psb_success_) call lv%aggr%default()
 

amgprec/impl/level/amg_s_base_onelev_dump.f90

@@ -127,8 +127,7 @@ subroutine amg_s_base_onelev_dump(lv,level,info,prefix,head,ac,rp,&
         ivr = lv%linmap%p_desc_U%get_global_indices(owned=.false.)
         write(fname(lname+1:),'(a,i3.3,a)')'_l',level,'_tprol.mtx'
         !
-        !  This is not implemented yet.
-        !call lv%tprol%print(fname,head=head,ivr=ivr)
+        call lv%tprol%print(fname,head=head,ivr=ivr)
       end if
     end if
   else
@@ -151,8 +150,7 @@ subroutine amg_s_base_onelev_dump(lv,level,info,prefix,head,ac,rp,&
       if (tprol_) then 
         write(fname(lname+1:),'(a,i3.3,a)')'_l',level,'_tprol.mtx'
         !
-        !  This is not implemented yet.
-        !call lv%tprol%print(fname,head=head)
+        call lv%tprol%print(fname,head=head)
       end if
     end if
   end if

amgprec/impl/level/amg_z_base_onelev_csetc.F90

@@ -261,6 +261,8 @@ subroutine amg_z_base_onelev_csetc(lv,what,val,info,pos,idx)
       if (info == 0) deallocate(lv%aggr,stat=info)
       if (info /= 0) then
         info = psb_err_internal_error_
+        call psb_errpush(info,name,a_err='aggregator deallocation?')
+        goto 9999
         return
       end if
     end if
@@ -270,8 +272,16 @@ subroutine amg_z_base_onelev_csetc(lv,what,val,info,pos,idx)
       allocate(amg_z_dec_aggregator_type :: lv%aggr, stat=info)
     case('SYMDEC')
       allocate(amg_z_symdec_aggregator_type :: lv%aggr, stat=info)
-    case default
+#if !defined(SERIAL_MPI)
+#endif
+  case default
       info =  psb_err_internal_error_
+#if !defined(SERIAL_MPI)
+      call psb_errpush(info,name,a_err='Unsupported PAR_AGGR_ALG')
+#else
+      call psb_errpush(info,name,a_err='PAR_AGGR_ALG unsupported (SERIAL_MPI on)')
+#endif
+      goto 9999
     end select
     if (info == psb_success_) call lv%aggr%default()
 

amgprec/impl/level/amg_z_base_onelev_dump.f90

@@ -127,8 +127,7 @@ subroutine amg_z_base_onelev_dump(lv,level,info,prefix,head,ac,rp,&
         ivr = lv%linmap%p_desc_U%get_global_indices(owned=.false.)
         write(fname(lname+1:),'(a,i3.3,a)')'_l',level,'_tprol.mtx'
         !
-        !  This is not implemented yet.
-        !call lv%tprol%print(fname,head=head,ivr=ivr)
+        call lv%tprol%print(fname,head=head,ivr=ivr)
       end if
     end if
   else
@@ -151,8 +150,7 @@ subroutine amg_z_base_onelev_dump(lv,level,info,prefix,head,ac,rp,&
       if (tprol_) then 
         write(fname(lname+1:),'(a,i3.3,a)')'_l',level,'_tprol.mtx'
         !
-        !  This is not implemented yet.
-        !call lv%tprol%print(fname,head=head)
+        call lv%tprol%print(fname,head=head)
       end if
     end if
   end if

amgprec/impl/smoother/amg_d_poly_smoother_apply_vect.f90

@@ -140,7 +140,7 @@ subroutine amg_d_poly_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
     call tz%zero()
 
     select case(sm%variant)
-    case(amg_poly_lottes_)
+    case(amg_cheb_4_)
       if (do_timings) call psb_tic(poly_1)
       block
         real(psb_dpk_)      :: cz, cr
@@ -155,7 +155,7 @@ subroutine amg_d_poly_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
           cz = (2*i*done-3)/(2*i*done+done)
           cr = (8*i*done-4)/((2*i*done+done)*sm%rho_ba)
           if (do_timings) call psb_tic(poly_vect)
-          call psb_abgdxyz(cr,cz,done,done,ty,tz,tx,desc_data,info) ! zk = cz * zk-1 + cr * rk-1 
+          call psb_upd_xyz(cr,cz,done,done,ty,tz,tx,desc_data,info) ! zk = cz * zk-1 + cr * rk-1 
           if (do_timings) call psb_toc(poly_vect)
           if (do_timings) call psb_tic(poly_mv)
           call psb_spmm(-done,sm%pa,tz,done,r,desc_data,info,work=aux,trans=trans_)
@@ -167,12 +167,12 @@ subroutine amg_d_poly_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
         cz = (2*sm%pdegree*done-3)/(2*sm%pdegree*done+done)
         cr = (8*sm%pdegree*done-4)/((2*sm%pdegree*done+done)*sm%rho_ba)
         if (do_timings) call psb_tic(poly_vect)
-        call psb_abgdxyz(cr,cz,done,done,ty,tz,tx,desc_data,info)
+        call psb_upd_xyz(cr,cz,done,done,ty,tz,tx,desc_data,info)
         if (do_timings) call psb_toc(poly_vect)
       end block
       if (do_timings) call psb_toc(poly_1)
 
-    case(amg_poly_lottes_beta_)
+    case(amg_cheb_4_opt_)
       if (do_timings) call psb_tic(poly_2)
       block
         real(psb_dpk_)      :: cz, cr
@@ -195,7 +195,7 @@ subroutine amg_d_poly_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
           cz = (2*i*done-3)/(2*i*done+done)
           cr = (8*i*done-4)/((2*i*done+done)*sm%rho_ba)
           if (do_timings) call psb_tic(poly_vect)
-          call psb_abgdxyz(cr,cz,sm%poly_beta(i),done,ty,tz,tx,desc_data,info)
+          call psb_upd_xyz(cr,cz,sm%poly_beta(i),done,ty,tz,tx,desc_data,info)
           if (do_timings) call psb_toc(poly_vect)
           if (do_timings) call psb_tic(poly_mv)
           call psb_spmm(-done,sm%pa,tz,done,r,desc_data,info,work=aux,trans=trans_)
@@ -205,11 +205,11 @@ subroutine amg_d_poly_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
         cz = (2*sm%pdegree*done-3)/(2*sm%pdegree*done+done)
         cr = (8*sm%pdegree*done-4)/((2*sm%pdegree*done+done)*sm%rho_ba)
         if (do_timings) call psb_tic(poly_vect)
-        call psb_abgdxyz(cr,cz,sm%poly_beta(sm%pdegree),done,ty,tz,tx,desc_data,info)
+        call psb_upd_xyz(cr,cz,sm%poly_beta(sm%pdegree),done,ty,tz,tx,desc_data,info)
         if (do_timings) call psb_toc(poly_vect)
       end block
       if (do_timings) call psb_toc(poly_2)
-    case(amg_poly_new_)
+    case(amg_cheb_1_opt_)
       if (do_timings) call psb_tic(poly_3)
       block
         real(psb_dpk_)      :: sigma, theta, delta, rho_old, rho
@@ -226,7 +226,7 @@ subroutine amg_d_poly_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
         if (do_timings) call psb_toc(poly_sv)
         call psb_geaxpby((done/sm%rho_ba),ty,dzero,r,desc_data,info)
         if (do_timings) call psb_tic(poly_vect)
-        call psb_abgdxyz((done/theta),dzero,done,done,r,tz,tx,desc_data,info)
+        call psb_upd_xyz((done/theta),dzero,done,done,r,tz,tx,desc_data,info)
         if (do_timings) call psb_toc(poly_vect)
 
         ! tz == d
@@ -244,7 +244,7 @@ subroutine amg_d_poly_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
           ! d_{k+1} = (rho rho_old) d_k + 2(rho/delta) r_{k+1}
           rho = done/(2*sigma - rho_old)
           if (do_timings) call psb_tic(poly_vect)
-          call psb_abgdxyz((2*rho/delta),(rho*rho_old),done,done,r,tz,tx,desc_data,info)
+          call psb_upd_xyz((2*rho/delta),(rho*rho_old),done,done,r,tz,tx,desc_data,info)
           if (do_timings) call psb_toc(poly_vect)
           rho_old = rho
         end do

amgprec/impl/smoother/amg_d_poly_smoother_bld.f90

@@ -75,9 +75,9 @@ subroutine amg_d_poly_smoother_bld(a,desc_a,sm,info,amold,vmold,imold)
   nrow_a = a%get_nrows()
   nztota = a%get_nzeros()
   select case(sm%variant)
-  case(amg_poly_lottes_)
+  case(amg_cheb_4_)
     ! do nothing
-  case(amg_poly_lottes_beta_)
+  case(amg_cheb_4_opt_)
     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)
@@ -87,7 +87,7 @@ subroutine amg_d_poly_smoother_bld(a,desc_a,sm,info,amold,vmold,imold)
            & a_err='invalid sm%degree for poly_beta')
       goto 9999
     end if
-  case(amg_poly_new_)
+  case(amg_cheb_1_opt_)
 
     if ((1<=sm%pdegree).and.(sm%pdegree<=30)) then
       !Ok

amgprec/impl/smoother/amg_d_poly_smoother_cseti.f90

@@ -58,11 +58,11 @@ subroutine amg_d_poly_smoother_cseti(sm,what,val,info,idx)
     sm%pdegree = val
   case('POLY_VARIANT')
     select case(val)
-    case(amg_poly_lottes_,amg_poly_lottes_beta_,amg_poly_new_)
+    case(amg_cheb_4_,amg_cheb_4_opt_,amg_cheb_1_opt_)
       sm%variant = val
     case default
-      write(0,*) 'Invalid choice for POLY_VARIANT, defaulting to amg_poly_lottes_',val
-      sm%variant = amg_poly_lottes_
+      write(0,*) 'Invalid choice for POLY_VARIANT, defaulting to amg_cheb_4_',val
+      sm%variant = amg_cheb_4_
     end select
   case('POLY_RHO_ESTIMATE')
     select case(val)

amgprec/impl/smoother/amg_d_poly_smoother_descr.f90

@@ -77,17 +77,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'
+  case(amg_cheb_4_)
+    write(iout_,*) trim(prefix_), '        variant:        ','CHEB_4'
     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'
+  case(amg_cheb_4_opt_)
+    write(iout_,*) trim(prefix_), '        variant:        ','CHEB_4_OPT'
     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'
+  case(amg_cheb_1_opt_)
+    write(iout_,*) trim(prefix_), '        variant:        ','CHEB_1_OPT'
     write(iout_,*) trim(prefix_), '         Degree:        ',sm%pdegree
     write(iout_,*) trim(prefix_), '         rho_ba:        ',sm%rho_ba
     write(iout_,*) trim(prefix_), '    Coefficient:        ',sm%cf_a

amgprec/impl/smoother/amg_d_poly_smoother_dmp.f90

@@ -77,11 +77,9 @@ subroutine amg_d_poly_smoother_dmp(sm,desc,level,info,prefix,head,smoother,solve
   end if
   lname = len_trim(prefix_)
   fname = trim(prefix_)
-  write(fname(lname+1:lname+5),'(a,i3.3)') '_poly',iam
+  write(fname(lname+1:lname+8),'(a,i3.3)') '_poly',iam
   lname = lname + 8
   ! to be completed
-  
-
 
   ! At base level do nothing for the smoother
   if (allocated(sm%sv)) &

amgprec/impl/smoother/amg_s_poly_smoother_apply_vect.f90

@@ -140,7 +140,7 @@ subroutine amg_s_poly_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
     call tz%zero()
 
     select case(sm%variant)
-    case(amg_poly_lottes_)
+    case(amg_cheb_4_)
       if (do_timings) call psb_tic(poly_1)
       block
         real(psb_spk_)      :: cz, cr
@@ -155,7 +155,7 @@ subroutine amg_s_poly_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
           cz = (2*i*sone-3)/(2*i*sone+sone)
           cr = (8*i*sone-4)/((2*i*sone+sone)*sm%rho_ba)
           if (do_timings) call psb_tic(poly_vect)
-          call psb_abgdxyz(cr,cz,sone,sone,ty,tz,tx,desc_data,info) ! zk = cz * zk-1 + cr * rk-1 
+          call psb_upd_xyz(cr,cz,sone,sone,ty,tz,tx,desc_data,info) ! zk = cz * zk-1 + cr * rk-1 
           if (do_timings) call psb_toc(poly_vect)
           if (do_timings) call psb_tic(poly_mv)
           call psb_spmm(-sone,sm%pa,tz,sone,r,desc_data,info,work=aux,trans=trans_)
@@ -167,12 +167,12 @@ subroutine amg_s_poly_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
         cz = (2*sm%pdegree*sone-3)/(2*sm%pdegree*sone+sone)
         cr = (8*sm%pdegree*sone-4)/((2*sm%pdegree*sone+sone)*sm%rho_ba)
         if (do_timings) call psb_tic(poly_vect)
-        call psb_abgdxyz(cr,cz,sone,sone,ty,tz,tx,desc_data,info)
+        call psb_upd_xyz(cr,cz,sone,sone,ty,tz,tx,desc_data,info)
         if (do_timings) call psb_toc(poly_vect)
       end block
       if (do_timings) call psb_toc(poly_1)
 
-    case(amg_poly_lottes_beta_)
+    case(amg_cheb_4_opt_)
       if (do_timings) call psb_tic(poly_2)
       block
         real(psb_spk_)      :: cz, cr
@@ -195,7 +195,7 @@ subroutine amg_s_poly_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
           cz = (2*i*sone-3)/(2*i*sone+sone)
           cr = (8*i*sone-4)/((2*i*sone+sone)*sm%rho_ba)
           if (do_timings) call psb_tic(poly_vect)
-          call psb_abgdxyz(cr,cz,sm%poly_beta(i),sone,ty,tz,tx,desc_data,info)
+          call psb_upd_xyz(cr,cz,sm%poly_beta(i),sone,ty,tz,tx,desc_data,info)
           if (do_timings) call psb_toc(poly_vect)
           if (do_timings) call psb_tic(poly_mv)
           call psb_spmm(-sone,sm%pa,tz,sone,r,desc_data,info,work=aux,trans=trans_)
@@ -205,11 +205,11 @@ subroutine amg_s_poly_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
         cz = (2*sm%pdegree*sone-3)/(2*sm%pdegree*sone+sone)
         cr = (8*sm%pdegree*sone-4)/((2*sm%pdegree*sone+sone)*sm%rho_ba)
         if (do_timings) call psb_tic(poly_vect)
-        call psb_abgdxyz(cr,cz,sm%poly_beta(sm%pdegree),sone,ty,tz,tx,desc_data,info)
+        call psb_upd_xyz(cr,cz,sm%poly_beta(sm%pdegree),sone,ty,tz,tx,desc_data,info)
         if (do_timings) call psb_toc(poly_vect)
       end block
       if (do_timings) call psb_toc(poly_2)
-    case(amg_poly_new_)
+    case(amg_cheb_1_opt_)
       if (do_timings) call psb_tic(poly_3)
       block
         real(psb_spk_)      :: sigma, theta, delta, rho_old, rho
@@ -226,7 +226,7 @@ subroutine amg_s_poly_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
         if (do_timings) call psb_toc(poly_sv)
         call psb_geaxpby((sone/sm%rho_ba),ty,szero,r,desc_data,info)
         if (do_timings) call psb_tic(poly_vect)
-        call psb_abgdxyz((sone/theta),szero,sone,sone,r,tz,tx,desc_data,info)
+        call psb_upd_xyz((sone/theta),szero,sone,sone,r,tz,tx,desc_data,info)
         if (do_timings) call psb_toc(poly_vect)
 
         ! tz == d
@@ -244,7 +244,7 @@ subroutine amg_s_poly_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
           ! d_{k+1} = (rho rho_old) d_k + 2(rho/delta) r_{k+1}
           rho = sone/(2*sigma - rho_old)
           if (do_timings) call psb_tic(poly_vect)
-          call psb_abgdxyz((2*rho/delta),(rho*rho_old),sone,sone,r,tz,tx,desc_data,info)
+          call psb_upd_xyz((2*rho/delta),(rho*rho_old),sone,sone,r,tz,tx,desc_data,info)
           if (do_timings) call psb_toc(poly_vect)
           rho_old = rho
         end do

amgprec/impl/smoother/amg_s_poly_smoother_bld.f90

@@ -75,9 +75,9 @@ subroutine amg_s_poly_smoother_bld(a,desc_a,sm,info,amold,vmold,imold)
   nrow_a = a%get_nrows()
   nztota = a%get_nzeros()
   select case(sm%variant)
-  case(amg_poly_lottes_)
+  case(amg_cheb_4_)
     ! do nothing
-  case(amg_poly_lottes_beta_)
+  case(amg_cheb_4_opt_)
     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)
@@ -87,7 +87,7 @@ subroutine amg_s_poly_smoother_bld(a,desc_a,sm,info,amold,vmold,imold)
            & a_err='invalid sm%degree for poly_beta')
       goto 9999
     end if
-  case(amg_poly_new_)
+  case(amg_cheb_1_opt_)
 
     if ((1<=sm%pdegree).and.(sm%pdegree<=30)) then
       !Ok

amgprec/impl/smoother/amg_s_poly_smoother_cseti.f90

@@ -58,11 +58,11 @@ subroutine amg_s_poly_smoother_cseti(sm,what,val,info,idx)
     sm%pdegree = val
   case('POLY_VARIANT')
     select case(val)
-    case(amg_poly_lottes_,amg_poly_lottes_beta_,amg_poly_new_)
+    case(amg_cheb_4_,amg_cheb_4_opt_,amg_cheb_1_opt_)
       sm%variant = val
     case default
-      write(0,*) 'Invalid choice for POLY_VARIANT, defaulting to amg_poly_lottes_',val
-      sm%variant = amg_poly_lottes_
+      write(0,*) 'Invalid choice for POLY_VARIANT, defaulting to amg_cheb_4_',val
+      sm%variant = amg_cheb_4_
     end select
   case('POLY_RHO_ESTIMATE')
     select case(val)

amgprec/impl/smoother/amg_s_poly_smoother_descr.f90

@@ -77,17 +77,17 @@ subroutine amg_s_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'
+  case(amg_cheb_4_)
+    write(iout_,*) trim(prefix_), '        variant:        ','CHEB_4'
     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'
+  case(amg_cheb_4_opt_)
+    write(iout_,*) trim(prefix_), '        variant:        ','CHEB_4_OPT'
     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'
+  case(amg_cheb_1_opt_)
+    write(iout_,*) trim(prefix_), '        variant:        ','CHEB_1_OPT'
     write(iout_,*) trim(prefix_), '         Degree:        ',sm%pdegree
     write(iout_,*) trim(prefix_), '         rho_ba:        ',sm%rho_ba
     write(iout_,*) trim(prefix_), '    Coefficient:        ',sm%cf_a

amgprec/impl/smoother/amg_s_poly_smoother_dmp.f90

@@ -77,11 +77,9 @@ subroutine amg_s_poly_smoother_dmp(sm,desc,level,info,prefix,head,smoother,solve
   end if
   lname = len_trim(prefix_)
   fname = trim(prefix_)
-  write(fname(lname+1:lname+5),'(a,i3.3)') '_poly',iam
+  write(fname(lname+1:lname+8),'(a,i3.3)') '_poly',iam
   lname = lname + 8
   ! to be completed
-  
-
 
   ! At base level do nothing for the smoother
   if (allocated(sm%sv)) &

amgprec/impl/solver/amg_c_ilu_solver_bld.f90

@@ -53,7 +53,6 @@ subroutine amg_c_ilu_solver_bld(a,desc_a,sv,info,b,amold,vmold,imold)
   class(psb_i_base_vect_type), intent(in), optional   :: imold
   ! Local variables
   integer(psb_ipk_) :: n_row,n_col, nrow_a, nztota, psb_fctype
-!!$    complex(psb_spk_), pointer :: ww(:), aux(:), tx(:),ty(:)
   type(psb_ctxt_type) :: ctxt
   integer(psb_ipk_)   :: np, me, i, err_act, debug_unit, debug_level
   character(len=20)   :: name='c_ilu_solver_bld', ch_err

amgprec/impl/solver/amg_d_ilu_solver_bld.f90

@@ -53,7 +53,6 @@ subroutine amg_d_ilu_solver_bld(a,desc_a,sv,info,b,amold,vmold,imold)
   class(psb_i_base_vect_type), intent(in), optional   :: imold
   ! Local variables
   integer(psb_ipk_) :: n_row,n_col, nrow_a, nztota, psb_fctype
-!!$    real(psb_dpk_), pointer :: ww(:), aux(:), tx(:),ty(:)
   type(psb_ctxt_type) :: ctxt
   integer(psb_ipk_)   :: np, me, i, err_act, debug_unit, debug_level
   character(len=20)   :: name='d_ilu_solver_bld', ch_err

amgprec/impl/solver/amg_s_ilu_solver_bld.f90

@@ -53,7 +53,6 @@ subroutine amg_s_ilu_solver_bld(a,desc_a,sv,info,b,amold,vmold,imold)
   class(psb_i_base_vect_type), intent(in), optional   :: imold
   ! Local variables
   integer(psb_ipk_) :: n_row,n_col, nrow_a, nztota, psb_fctype
-!!$    real(psb_spk_), pointer :: ww(:), aux(:), tx(:),ty(:)
   type(psb_ctxt_type) :: ctxt
   integer(psb_ipk_)   :: np, me, i, err_act, debug_unit, debug_level
   character(len=20)   :: name='s_ilu_solver_bld', ch_err

amgprec/impl/solver/amg_z_ilu_solver_bld.f90

@@ -53,7 +53,6 @@ subroutine amg_z_ilu_solver_bld(a,desc_a,sv,info,b,amold,vmold,imold)
   class(psb_i_base_vect_type), intent(in), optional   :: imold
   ! Local variables
   integer(psb_ipk_) :: n_row,n_col, nrow_a, nztota, psb_fctype
-!!$    complex(psb_dpk_), pointer :: ww(:), aux(:), tx(:),ty(:)
   type(psb_ctxt_type) :: ctxt
   integer(psb_ipk_)   :: np, me, i, err_act, debug_unit, debug_level
   character(len=20)   :: name='z_ilu_solver_bld', ch_err

samples/advanced/pdegen/Makefile

@@ -8,32 +8,37 @@ 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 	
+DGEN2D=amg_d_pde2d_poisson_mod.o amg_d_pde2d_exp_mod.o \
+  amg_d_pde2d_gauss_mod.o amg_d_pde2d_box_mod.o
+DGEN3D=amg_d_pde3d_poisson_mod.o amg_d_pde3d_exp_mod.o \
+  amg_d_pde3d_gauss_mod.o amg_d_pde3d_box_mod.o 	
+SGEN2D=amg_s_pde2d_poisson_mod.o amg_s_pde2d_exp_mod.o \
+  amg_s_pde2d_gauss_mod.o amg_s_pde2d_box_mod.o
+SGEN3D=amg_s_pde3d_poisson_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 $(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)
+	$(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 $(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)
+	$(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  $(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)
+	$(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 $(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)
+	$(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
-	$(FLINK) $(LINKOPT) amg_d_pde3d_rebld.o  data_input.o -o amg_d_pde3d_rebld $(AMG_LIBS) $(PSBLAS_LIBS) $(LDLIBS)
-	/bin/mv amg_d_pde3d_rebld $(EXEDIR)
 
 amg_d_pde3d.o amg_s_pde3d.o amg_d_pde2d.o amg_s_pde2d.o: data_input.o
 
@@ -42,6 +47,11 @@ 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)
 
+amg_d_genpde_mod.o: $(DGEN3D)
+amg_s_genpde_mod.o: $(SGEN3D)
+amg_d_genpde_mod.o: $(DGEN2D)
+amg_s_genpde_mod.o: $(SGEN2D)
+
 check: all
 	cd runs && ./amg_d_pde2d <amg_pde2d.inp && ./amg_s_pde2d<amg_pde2d.inp
 

samples/advanced/pdegen/amg_d_pde2d.F90

@@ -69,7 +69,7 @@ program amg_d_pde2d
   use psb_krylov_mod
   use psb_util_mod
   use data_input
-  use amg_d_pde2d_base_mod
+  use amg_d_pde2d_poisson_mod
   use amg_d_pde2d_exp_mod
   use amg_d_pde2d_box_mod
   use amg_d_pde2d_gauss_mod
@@ -81,7 +81,8 @@ program amg_d_pde2d
 
   ! input parameters
   character(len=20) :: kmethd, ptype
-  character(len=5)  :: afmt, pdecoeff
+  character(len=5)  :: afmt
+  character(len=32)  :: pdecoeff
   integer(psb_ipk_) :: idim
   integer(psb_epk_) :: system_size
 
@@ -146,6 +147,8 @@ program amg_d_pde2d
     integer(psb_ipk_)  :: jsweeps      ! (pre-)smoother / 1-lev prec. sweeps
     integer(psb_ipk_)  :: degree       ! degree for polynomial smoother
     character(len=32)  :: pvariant     ! polynomial  variant
+    character(len=32)  :: prhovariant  ! how to estimate rho(M^{-1}A)
+    real(psb_dpk_)     :: prhovalue    ! if previous is set value, we set it from this one
     integer(psb_ipk_)  :: novr         ! number of overlap layers
     character(len=32)  :: restr        ! restriction over application of AS
     character(len=32)  :: prol         ! prolongation over application of AS
@@ -162,6 +165,8 @@ program amg_d_pde2d
     integer(psb_ipk_)  :: jsweeps2     ! post-smoother sweeps
     integer(psb_ipk_)  :: degree2      ! degree for polynomial smoother
     character(len=32)  :: pvariant2    ! polynomial  variant
+    character(len=32)  :: prhovariant2 ! how to estimate rho(M^{-1}A)
+    real(psb_dpk_)     :: prhovalue2   ! if previous is set value, we set it from this one
     integer(psb_ipk_)  :: novr2        ! number of overlap layers
     character(len=32)  :: restr2       ! restriction  over application of AS
     character(len=32)  :: prol2        ! prolongation over application of AS
@@ -244,18 +249,22 @@ program amg_d_pde2d
   call psb_barrier(ctxt)
   t1 = psb_wtime()
   select case(psb_toupper(trim(pdecoeff)))
-  case("CONST")
+  case("POISSON")
     call amg_gen_pde2d(ctxt,idim,a,b,x,desc_a,afmt,&
-        & a1_base,a2_base,b1_base,b2_base,c_base,g_base,info)
+         & a1_poisson,a2_poisson,&
+         & b1_poisson,b2_poisson,c_poisson,g_poisson,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)
+         & 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)
+         & 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)
+         & 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'
@@ -344,7 +353,12 @@ program amg_d_pde2d
     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)
-    
+    if (p_choice%prhovalue > dzero ) then
+            call prec%set('poly_rho_ba', p_choice%prhovalue, info)
+    else
+            call prec%set('poly_rho_estimate', p_choice%prhovariant, info)
+    end if
+
     select case (psb_toupper(p_choice%smther))
     case ('GS','BWGS','FBGS','JACOBI','L1-JACOBI','L1-FBGS')
       ! do nothing
@@ -376,6 +390,12 @@ program amg_d_pde2d
       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')
+      if (p_choice%prhovalue > dzero ) then
+         call prec%set('poly_rho_ba', p_choice%prhovalue2, info,pos='post')
+      else
+         call prec%set('poly_rho_estimate', p_choice%prhovariant2, info,pos='post')
+      end if
+
       select case (psb_toupper(p_choice%smther2))
       case ('GS','BWGS','FBGS','JACOBI','L1-JACOBI','L1-FBGS')
         ! do nothing
@@ -592,7 +612,9 @@ 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)     ! (pre-)smoother / 1-lev prec sweeps
-      call read_data(prec%pvariant,inp_unit)   ! 
+      call read_data(prec%pvariant,inp_unit)   !
+      call read_data(prec%prhovariant,inp_unit)! how to estimate rho(M^{-1}A)
+      call read_data(prec%prhovalue,inp_unit)  ! if previous is set value, we set it from this one
       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
@@ -607,6 +629,8 @@ contains
       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%prhovariant2,inp_unit)! how to estimate rho(M^{-1}A)
+      call read_data(prec%prhovalue2,inp_unit)  ! if previous is set value, we set it from this one
       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
@@ -679,6 +703,8 @@ contains
     call psb_bcast(ctxt,prec%jsweeps)
     call psb_bcast(ctxt,prec%degree)
     call psb_bcast(ctxt,prec%pvariant)
+    call psb_bcast(ctxt,prec%prhovariant)
+    call psb_bcast(ctxt,prec%prhovalue)  
     call psb_bcast(ctxt,prec%novr)
     call psb_bcast(ctxt,prec%restr)
     call psb_bcast(ctxt,prec%prol)
@@ -693,6 +719,8 @@ contains
     call psb_bcast(ctxt,prec%jsweeps2)
     call psb_bcast(ctxt,prec%degree2)
     call psb_bcast(ctxt,prec%pvariant2)
+    call psb_bcast(ctxt,prec%prhovariant2)
+    call psb_bcast(ctxt,prec%prhovalue2)  
     call psb_bcast(ctxt,prec%novr2)
     call psb_bcast(ctxt,prec%restr2)
     call psb_bcast(ctxt,prec%prol2)

samples/advanced/pdegen/amg_d_pde2d_poisson_mod.f90

@@ -34,56 +34,56 @@
 !    ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
 !    POSSIBILITY OF SUCH DAMAGE.
 !
-module amg_d_pde2d_base_mod
+module amg_d_pde2d_poisson_mod
   use psb_base_mod, only : psb_dpk_, dzero, done
   real(psb_dpk_), save, private :: epsilon=done/80
 contains
-  subroutine pde_set_parm2d_base(dat)
+  subroutine pde_set_parm2d_poisson(dat)
     real(psb_dpk_), intent(in) :: dat
     epsilon = dat
-  end subroutine pde_set_parm2d_base
+  end subroutine pde_set_parm2d_poisson
   !
   ! functions parametrizing the differential equation
   !
-  function b1_base(x,y)
+  function b1_poisson(x,y)
     implicit none 
-    real(psb_dpk_) :: b1_base
+    real(psb_dpk_) :: b1_poisson
     real(psb_dpk_), intent(in) :: x,y
-    b1_base = dzero/1.414_psb_dpk_
-  end function b1_base
-  function b2_base(x,y)
+    b1_poisson = dzero
+  end function b1_poisson
+  function b2_poisson(x,y)
     implicit none 
-    real(psb_dpk_) ::  b2_base
+    real(psb_dpk_) ::  b2_poisson
     real(psb_dpk_), intent(in) :: x,y
-    b2_base = dzero/1.414_psb_dpk_
-  end function b2_base
-  function c_base(x,y)
+    b2_poisson = dzero
+  end function b2_poisson
+  function c_poisson(x,y)
     implicit none 
-    real(psb_dpk_) ::  c_base
+    real(psb_dpk_) ::  c_poisson
     real(psb_dpk_), intent(in) :: x,y
-    c_base = dzero
-  end function c_base
-  function a1_base(x,y)
+    c_poisson = dzero
+  end function c_poisson
+  function a1_poisson(x,y)
     implicit none 
-    real(psb_dpk_) ::  a1_base
+    real(psb_dpk_) ::  a1_poisson
     real(psb_dpk_), intent(in) :: x,y
-    a1_base=done*epsilon
-  end function a1_base
-  function a2_base(x,y)
+    a1_poisson=done*epsilon
+  end function a1_poisson
+  function a2_poisson(x,y)
     implicit none 
-    real(psb_dpk_) ::  a2_base
+    real(psb_dpk_) ::  a2_poisson
     real(psb_dpk_), intent(in) :: x,y
-    a2_base=done*epsilon
-  end function a2_base
-  function g_base(x,y)
+    a2_poisson=done*epsilon
+  end function a2_poisson
+  function g_poisson(x,y)
     implicit none 
-    real(psb_dpk_) ::  g_base
+    real(psb_dpk_) ::  g_poisson
     real(psb_dpk_), intent(in) :: x,y
-    g_base = dzero
+    g_poisson = dzero
     if (x == done) then
-      g_base = done
+      g_poisson = done
     else if (x == dzero) then
-      g_base = done
+      g_poisson = done
     end if
-  end function g_base
-end module amg_d_pde2d_base_mod
+  end function g_poisson
+end module amg_d_pde2d_poisson_mod

samples/advanced/pdegen/amg_d_pde3d.F90

@@ -70,7 +70,7 @@ program amg_d_pde3d
   use psb_krylov_mod
   use psb_util_mod
   use data_input
-  use amg_d_pde3d_base_mod
+  use amg_d_pde3d_poisson_mod
   use amg_d_pde3d_exp_mod
   use amg_d_pde3d_box_mod
   use amg_d_pde3d_gauss_mod
@@ -82,7 +82,8 @@ program amg_d_pde3d
 
   ! input parameters
   character(len=20) :: kmethd, ptype
-  character(len=5)  :: afmt, pdecoeff
+  character(len=5)  :: afmt
+  character(len=32)  :: pdecoeff
   integer(psb_ipk_) :: idim
   integer(psb_epk_) :: system_size
 
@@ -147,6 +148,8 @@ program amg_d_pde3d
     integer(psb_ipk_)  :: jsweeps      ! (pre-)smoother / 1-lev prec. sweeps
     integer(psb_ipk_)  :: degree       ! degree for polynomial smoother
     character(len=32)  :: pvariant     ! polynomial  variant
+    character(len=32)  :: prhovariant  ! how to estimate rho(M^{-1}A)
+    real(psb_dpk_)     :: prhovalue    ! if previous is set value, we set it from this one
     integer(psb_ipk_)  :: novr         ! number of overlap layers
     character(len=32)  :: restr        ! restriction over application of AS
     character(len=32)  :: prol         ! prolongation over application of AS
@@ -163,6 +166,8 @@ program amg_d_pde3d
     integer(psb_ipk_)  :: jsweeps2     ! post-smoother sweeps
     integer(psb_ipk_)  :: degree2      ! degree for polynomial smoother
     character(len=32)  :: pvariant2    ! polynomial  variant
+    character(len=32)  :: prhovariant2 ! how to estimate rho(M^{-1}A)
+    real(psb_dpk_)     :: prhovalue2   ! if previous is set value, we set it from this one
     integer(psb_ipk_)  :: novr2        ! number of overlap layers
     character(len=32)  :: restr2       ! restriction  over application of AS
     character(len=32)  :: prol2        ! prolongation over application of AS
@@ -246,18 +251,22 @@ program amg_d_pde3d
   call psb_barrier(ctxt)
   t1 = psb_wtime()
   select case(psb_toupper(trim(pdecoeff)))
-  case("CONST")
+  case("POISSON")
     call amg_gen_pde3d(ctxt,idim,a,b,x,desc_a,afmt,&
-        & a1_base,a2_base,a3_base,b1_base,b2_base,b3_base,c_base,g_base,info)
+         & a1_poisson,a2_poisson,a3_poisson,&
+         & b1_poisson,b2_poisson,b3_poisson,c_poisson,g_poisson,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)
+         & 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)
+         & 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)
+         & a1_gauss,a2_gauss,a3_gauss,&
+         & b1_gauss,b2_gauss,b3_gauss,c_gauss,g_gauss,info)
   case default
     info=psb_err_from_subroutine_
     ch_err='amg_gen_pdecoeff'
@@ -348,7 +357,12 @@ program amg_d_pde3d
     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)
-    
+    if (p_choice%prhovalue > dzero ) then
+            call prec%set('poly_rho_ba', p_choice%prhovalue, info)
+    else
+            call prec%set('poly_rho_estimate', p_choice%prhovariant, info)
+    end if
+
     select case (psb_toupper(p_choice%smther))
     case ('GS','BWGS','FBGS','JACOBI','L1-JACOBI','L1-FBGS')
       ! do nothing
@@ -380,6 +394,12 @@ program amg_d_pde3d
       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')
+      if (p_choice%prhovalue > dzero ) then
+         call prec%set('poly_rho_ba', p_choice%prhovalue2, info,pos='post')
+      else
+         call prec%set('poly_rho_estimate', p_choice%prhovariant2, info,pos='post')
+      end if
+
       select case (psb_toupper(p_choice%smther2))
       case ('GS','BWGS','FBGS','JACOBI','L1-JACOBI','L1-FBGS')
         ! do nothing
@@ -596,7 +616,9 @@ 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)     ! (pre-)smoother / 1-lev prec sweeps
-      call read_data(prec%pvariant,inp_unit)   ! 
+      call read_data(prec%pvariant,inp_unit)   !
+      call read_data(prec%prhovariant,inp_unit)! how to estimate rho(M^{-1}A)
+      call read_data(prec%prhovalue,inp_unit)  ! if previous is set value, we set it from this one
       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
@@ -611,6 +633,8 @@ contains
       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%prhovariant2,inp_unit)! how to estimate rho(M^{-1}A)
+      call read_data(prec%prhovalue2,inp_unit)  ! if previous is set value, we set it from this one
       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
@@ -683,6 +707,8 @@ contains
     call psb_bcast(ctxt,prec%jsweeps)
     call psb_bcast(ctxt,prec%degree)
     call psb_bcast(ctxt,prec%pvariant)
+    call psb_bcast(ctxt,prec%prhovariant)
+    call psb_bcast(ctxt,prec%prhovalue)  
     call psb_bcast(ctxt,prec%novr)
     call psb_bcast(ctxt,prec%restr)
     call psb_bcast(ctxt,prec%prol)
@@ -697,6 +723,8 @@ contains
     call psb_bcast(ctxt,prec%jsweeps2)
     call psb_bcast(ctxt,prec%degree2)
     call psb_bcast(ctxt,prec%pvariant2)
+    call psb_bcast(ctxt,prec%prhovariant2)
+    call psb_bcast(ctxt,prec%prhovalue2)  
     call psb_bcast(ctxt,prec%novr2)
     call psb_bcast(ctxt,prec%restr2)
     call psb_bcast(ctxt,prec%prol2)

samples/advanced/pdegen/amg_d_pde3d_poisson_mod.f90

@@ -34,68 +34,68 @@
 !    ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
 !    POSSIBILITY OF SUCH DAMAGE.
 !
-module amg_d_pde3d_base_mod
+module amg_d_pde3d_poisson_mod
   use psb_base_mod, only : psb_dpk_, done, dzero
   real(psb_dpk_), save, private :: epsilon=done/80
 contains
-  subroutine pde_set_parm3d_base(dat)
+  subroutine pde_set_parm3d_poisson(dat)
     real(psb_dpk_), intent(in) :: dat
     epsilon = dat
-  end subroutine pde_set_parm3d_base
+  end subroutine pde_set_parm3d_poisson
   !
   ! functions parametrizing the differential equation
   !
-  function b1_base(x,y,z)
+  function b1_poisson(x,y,z)
     implicit none 
-    real(psb_dpk_) :: b1_base
+    real(psb_dpk_) :: b1_poisson
     real(psb_dpk_), intent(in) :: x,y,z
-    b1_base=dzero/sqrt(3.0_psb_dpk_)
-  end function b1_base
-  function b2_base(x,y,z)
+    b1_poisson=dzero
+  end function b1_poisson
+  function b2_poisson(x,y,z)
     implicit none 
-    real(psb_dpk_) ::  b2_base
+    real(psb_dpk_) ::  b2_poisson
     real(psb_dpk_), intent(in) :: x,y,z
-    b2_base=dzero/sqrt(3.0_psb_dpk_)
-  end function b2_base
-  function b3_base(x,y,z)
+    b2_poisson=dzero
+  end function b2_poisson
+  function b3_poisson(x,y,z)
     implicit none 
-    real(psb_dpk_) ::  b3_base
+    real(psb_dpk_) ::  b3_poisson
     real(psb_dpk_), intent(in) :: x,y,z
-    b3_base=dzero/sqrt(3.0_psb_dpk_)
-  end function b3_base
-  function c_base(x,y,z)
+    b3_poisson=dzero
+  end function b3_poisson
+  function c_poisson(x,y,z)
     implicit none 
-    real(psb_dpk_) ::  c_base
+    real(psb_dpk_) ::  c_poisson
     real(psb_dpk_), intent(in) :: x,y,z
-    c_base=dzero
-  end function c_base
-  function a1_base(x,y,z)
+    c_poisson=dzero
+  end function c_poisson
+  function a1_poisson(x,y,z)
     implicit none 
-    real(psb_dpk_) ::  a1_base
+    real(psb_dpk_) ::  a1_poisson
     real(psb_dpk_), intent(in) :: x,y,z
-    a1_base=epsilon
-  end function a1_base
-  function a2_base(x,y,z)
+    a1_poisson=epsilon
+  end function a1_poisson
+  function a2_poisson(x,y,z)
     implicit none 
-    real(psb_dpk_) ::  a2_base
+    real(psb_dpk_) ::  a2_poisson
     real(psb_dpk_), intent(in) :: x,y,z
-    a2_base=epsilon
-  end function a2_base
-  function a3_base(x,y,z)
+    a2_poisson=epsilon
+  end function a2_poisson
+  function a3_poisson(x,y,z)
     implicit none 
-    real(psb_dpk_) ::  a3_base
+    real(psb_dpk_) ::  a3_poisson
     real(psb_dpk_), intent(in) :: x,y,z
-    a3_base=epsilon
-  end function a3_base
-  function g_base(x,y,z)
+    a3_poisson=epsilon
+  end function a3_poisson
+  function g_poisson(x,y,z)
     implicit none 
-    real(psb_dpk_) ::  g_base
+    real(psb_dpk_) ::  g_poisson
     real(psb_dpk_), intent(in) :: x,y,z
-    g_base = dzero
+    g_poisson = dzero
     if (x == done) then
-      g_base = done
+      g_poisson = done
     else if (x == dzero) then
-      g_base = done
+      g_poisson = done
     end if
-  end function g_base
-end module amg_d_pde3d_base_mod
+  end function g_poisson
+end module amg_d_pde3d_poisson_mod

samples/advanced/pdegen/amg_s_pde2d.F90

@@ -69,7 +69,7 @@ program amg_s_pde2d
   use psb_krylov_mod
   use psb_util_mod
   use data_input
-  use amg_s_pde2d_base_mod
+  use amg_s_pde2d_poisson_mod
   use amg_s_pde2d_exp_mod
   use amg_s_pde2d_box_mod
   use amg_s_pde2d_gauss_mod
@@ -81,7 +81,8 @@ program amg_s_pde2d
 
   ! input parameters
   character(len=20) :: kmethd, ptype
-  character(len=5)  :: afmt, pdecoeff
+  character(len=5)  :: afmt
+  character(len=32)  :: pdecoeff
   integer(psb_ipk_) :: idim
   integer(psb_epk_) :: system_size
 
@@ -146,6 +147,8 @@ program amg_s_pde2d
     integer(psb_ipk_)  :: jsweeps      ! (pre-)smoother / 1-lev prec. sweeps
     integer(psb_ipk_)  :: degree       ! degree for polynomial smoother
     character(len=32)  :: pvariant     ! polynomial  variant
+    character(len=32)  :: prhovariant  ! how to estimate rho(M^{-1}A)
+    real(psb_spk_)     :: prhovalue    ! if previous is set value, we set it from this one
     integer(psb_ipk_)  :: novr         ! number of overlap layers
     character(len=32)  :: restr        ! restriction over application of AS
     character(len=32)  :: prol         ! prolongation over application of AS
@@ -162,6 +165,8 @@ program amg_s_pde2d
     integer(psb_ipk_)  :: jsweeps2     ! post-smoother sweeps
     integer(psb_ipk_)  :: degree2      ! degree for polynomial smoother
     character(len=32)  :: pvariant2    ! polynomial  variant
+    character(len=32)  :: prhovariant2 ! how to estimate rho(M^{-1}A)
+    real(psb_spk_)     :: prhovalue2   ! if previous is set value, we set it from this one
     integer(psb_ipk_)  :: novr2        ! number of overlap layers
     character(len=32)  :: restr2       ! restriction  over application of AS
     character(len=32)  :: prol2        ! prolongation over application of AS
@@ -244,18 +249,22 @@ program amg_s_pde2d
   call psb_barrier(ctxt)
   t1 = psb_wtime()
   select case(psb_toupper(trim(pdecoeff)))
-  case("CONST")
+  case("POISSON")
     call amg_gen_pde2d(ctxt,idim,a,b,x,desc_a,afmt,&
-        & a1_base,a2_base,b1_base,b2_base,c_base,g_base,info)
+         & a1_poisson,a2_poisson,&
+         & b1_poisson,b2_poisson,c_poisson,g_poisson,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)
+         & 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)
+         & 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)
+         & 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'
@@ -344,7 +353,12 @@ program amg_s_pde2d
     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)
-    
+    if (p_choice%prhovalue > szero ) then
+            call prec%set('poly_rho_ba', p_choice%prhovalue, info)
+    else
+            call prec%set('poly_rho_estimate', p_choice%prhovariant, info)
+    end if
+
     select case (psb_toupper(p_choice%smther))
     case ('GS','BWGS','FBGS','JACOBI','L1-JACOBI','L1-FBGS')
       ! do nothing
@@ -376,6 +390,12 @@ program amg_s_pde2d
       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')
+      if (p_choice%prhovalue > szero ) then
+         call prec%set('poly_rho_ba', p_choice%prhovalue2, info,pos='post')
+      else
+         call prec%set('poly_rho_estimate', p_choice%prhovariant2, info,pos='post')
+      end if
+
       select case (psb_toupper(p_choice%smther2))
       case ('GS','BWGS','FBGS','JACOBI','L1-JACOBI','L1-FBGS')
         ! do nothing
@@ -592,7 +612,9 @@ 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)     ! (pre-)smoother / 1-lev prec sweeps
-      call read_data(prec%pvariant,inp_unit)   ! 
+      call read_data(prec%pvariant,inp_unit)   !
+      call read_data(prec%prhovariant,inp_unit)! how to estimate rho(M^{-1}A)
+      call read_data(prec%prhovalue,inp_unit)  ! if previous is set value, we set it from this one
       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
@@ -607,6 +629,8 @@ contains
       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%prhovariant2,inp_unit)! how to estimate rho(M^{-1}A)
+      call read_data(prec%prhovalue2,inp_unit)  ! if previous is set value, we set it from this one
       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
@@ -679,6 +703,8 @@ contains
     call psb_bcast(ctxt,prec%jsweeps)
     call psb_bcast(ctxt,prec%degree)
     call psb_bcast(ctxt,prec%pvariant)
+    call psb_bcast(ctxt,prec%prhovariant)
+    call psb_bcast(ctxt,prec%prhovalue)  
     call psb_bcast(ctxt,prec%novr)
     call psb_bcast(ctxt,prec%restr)
     call psb_bcast(ctxt,prec%prol)
@@ -693,6 +719,8 @@ contains
     call psb_bcast(ctxt,prec%jsweeps2)
     call psb_bcast(ctxt,prec%degree2)
     call psb_bcast(ctxt,prec%pvariant2)
+    call psb_bcast(ctxt,prec%prhovariant2)
+    call psb_bcast(ctxt,prec%prhovalue2)  
     call psb_bcast(ctxt,prec%novr2)
     call psb_bcast(ctxt,prec%restr2)
     call psb_bcast(ctxt,prec%prol2)

samples/advanced/pdegen/amg_s_pde2d_poisson_mod.f90

@@ -34,56 +34,56 @@
 !    ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
 !    POSSIBILITY OF SUCH DAMAGE.
 !
-module amg_s_pde2d_base_mod
+module amg_s_pde2d_poisson_mod
   use psb_base_mod, only : psb_spk_, szero, sone
   real(psb_spk_), save, private :: epsilon=sone/80
 contains
-  subroutine pde_set_parm2d_base(dat)
+  subroutine pde_set_parm2d_poisson(dat)
     real(psb_spk_), intent(in) :: dat
     epsilon = dat
-  end subroutine pde_set_parm2d_base
+  end subroutine pde_set_parm2d_poisson
   !
   ! functions parametrizing the differential equation
   !
-  function b1_base(x,y)
+  function b1_poisson(x,y)
     implicit none 
-    real(psb_spk_) :: b1_base
+    real(psb_spk_) :: b1_poisson
     real(psb_spk_), intent(in) :: x,y
-    b1_base = szero/1.414_psb_spk_
-  end function b1_base
-  function b2_base(x,y)
+    b1_poisson = szero
+  end function b1_poisson
+  function b2_poisson(x,y)
     implicit none 
-    real(psb_spk_) ::  b2_base
+    real(psb_spk_) ::  b2_poisson
     real(psb_spk_), intent(in) :: x,y
-    b2_base = szero/1.414_psb_spk_
-  end function b2_base
-  function c_base(x,y)
+    b2_poisson = szero
+  end function b2_poisson
+  function c_poisson(x,y)
     implicit none 
-    real(psb_spk_) ::  c_base
+    real(psb_spk_) ::  c_poisson
     real(psb_spk_), intent(in) :: x,y
-    c_base = szero
-  end function c_base
-  function a1_base(x,y)
+    c_poisson = szero
+  end function c_poisson
+  function a1_poisson(x,y)
     implicit none 
-    real(psb_spk_) ::  a1_base
+    real(psb_spk_) ::  a1_poisson
     real(psb_spk_), intent(in) :: x,y
-    a1_base=sone*epsilon
-  end function a1_base
-  function a2_base(x,y)
+    a1_poisson=sone*epsilon
+  end function a1_poisson
+  function a2_poisson(x,y)
     implicit none 
-    real(psb_spk_) ::  a2_base
+    real(psb_spk_) ::  a2_poisson
     real(psb_spk_), intent(in) :: x,y
-    a2_base=sone*epsilon
-  end function a2_base
-  function g_base(x,y)
+    a2_poisson=sone*epsilon
+  end function a2_poisson
+  function g_poisson(x,y)
     implicit none 
-    real(psb_spk_) ::  g_base
+    real(psb_spk_) ::  g_poisson
     real(psb_spk_), intent(in) :: x,y
-    g_base = szero
+    g_poisson = szero
     if (x == sone) then
-      g_base = sone
+      g_poisson = sone
     else if (x == szero) then
-      g_base = sone
+      g_poisson = sone
     end if
-  end function g_base
-end module amg_s_pde2d_base_mod
+  end function g_poisson
+end module amg_s_pde2d_poisson_mod

samples/advanced/pdegen/amg_s_pde3d.F90

@@ -70,7 +70,7 @@ program amg_s_pde3d
   use psb_krylov_mod
   use psb_util_mod
   use data_input
-  use amg_s_pde3d_base_mod
+  use amg_s_pde3d_poisson_mod
   use amg_s_pde3d_exp_mod
   use amg_s_pde3d_box_mod
   use amg_s_pde3d_gauss_mod
@@ -82,7 +82,8 @@ program amg_s_pde3d
 
   ! input parameters
   character(len=20) :: kmethd, ptype
-  character(len=5)  :: afmt, pdecoeff
+  character(len=5)  :: afmt
+  character(len=32)  :: pdecoeff
   integer(psb_ipk_) :: idim
   integer(psb_epk_) :: system_size
 
@@ -147,6 +148,8 @@ program amg_s_pde3d
     integer(psb_ipk_)  :: jsweeps      ! (pre-)smoother / 1-lev prec. sweeps
     integer(psb_ipk_)  :: degree       ! degree for polynomial smoother
     character(len=32)  :: pvariant     ! polynomial  variant
+    character(len=32)  :: prhovariant  ! how to estimate rho(M^{-1}A)
+    real(psb_spk_)     :: prhovalue    ! if previous is set value, we set it from this one
     integer(psb_ipk_)  :: novr         ! number of overlap layers
     character(len=32)  :: restr        ! restriction over application of AS
     character(len=32)  :: prol         ! prolongation over application of AS
@@ -163,6 +166,8 @@ program amg_s_pde3d
     integer(psb_ipk_)  :: jsweeps2     ! post-smoother sweeps
     integer(psb_ipk_)  :: degree2      ! degree for polynomial smoother
     character(len=32)  :: pvariant2    ! polynomial  variant
+    character(len=32)  :: prhovariant2 ! how to estimate rho(M^{-1}A)
+    real(psb_spk_)     :: prhovalue2   ! if previous is set value, we set it from this one
     integer(psb_ipk_)  :: novr2        ! number of overlap layers
     character(len=32)  :: restr2       ! restriction  over application of AS
     character(len=32)  :: prol2        ! prolongation over application of AS
@@ -246,18 +251,22 @@ program amg_s_pde3d
   call psb_barrier(ctxt)
   t1 = psb_wtime()
   select case(psb_toupper(trim(pdecoeff)))
-  case("CONST")
+  case("POISSON")
     call amg_gen_pde3d(ctxt,idim,a,b,x,desc_a,afmt,&
-        & a1_base,a2_base,a3_base,b1_base,b2_base,b3_base,c_base,g_base,info)
+         & a1_poisson,a2_poisson,a3_poisson,&
+         & b1_poisson,b2_poisson,b3_poisson,c_poisson,g_poisson,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)
+         & 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)
+         & 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)
+         & a1_gauss,a2_gauss,a3_gauss,&
+         & b1_gauss,b2_gauss,b3_gauss,c_gauss,g_gauss,info)
   case default
     info=psb_err_from_subroutine_
     ch_err='amg_gen_pdecoeff'
@@ -348,7 +357,12 @@ program amg_s_pde3d
     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)
-    
+    if (p_choice%prhovalue > szero ) then
+            call prec%set('poly_rho_ba', p_choice%prhovalue, info)
+    else
+            call prec%set('poly_rho_estimate', p_choice%prhovariant, info)
+    end if
+
     select case (psb_toupper(p_choice%smther))
     case ('GS','BWGS','FBGS','JACOBI','L1-JACOBI','L1-FBGS')
       ! do nothing
@@ -380,6 +394,12 @@ program amg_s_pde3d
       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')
+      if (p_choice%prhovalue > szero ) then
+         call prec%set('poly_rho_ba', p_choice%prhovalue2, info,pos='post')
+      else
+         call prec%set('poly_rho_estimate', p_choice%prhovariant2, info,pos='post')
+      end if
+
       select case (psb_toupper(p_choice%smther2))
       case ('GS','BWGS','FBGS','JACOBI','L1-JACOBI','L1-FBGS')
         ! do nothing
@@ -596,7 +616,9 @@ 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)     ! (pre-)smoother / 1-lev prec sweeps
-      call read_data(prec%pvariant,inp_unit)   ! 
+      call read_data(prec%pvariant,inp_unit)   !
+      call read_data(prec%prhovariant,inp_unit)! how to estimate rho(M^{-1}A)
+      call read_data(prec%prhovalue,inp_unit)  ! if previous is set value, we set it from this one
       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
@@ -611,6 +633,8 @@ contains
       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%prhovariant2,inp_unit)! how to estimate rho(M^{-1}A)
+      call read_data(prec%prhovalue2,inp_unit)  ! if previous is set value, we set it from this one
       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
@@ -683,6 +707,8 @@ contains
     call psb_bcast(ctxt,prec%jsweeps)
     call psb_bcast(ctxt,prec%degree)
     call psb_bcast(ctxt,prec%pvariant)
+    call psb_bcast(ctxt,prec%prhovariant)
+    call psb_bcast(ctxt,prec%prhovalue)  
     call psb_bcast(ctxt,prec%novr)
     call psb_bcast(ctxt,prec%restr)
     call psb_bcast(ctxt,prec%prol)
@@ -697,6 +723,8 @@ contains
     call psb_bcast(ctxt,prec%jsweeps2)
     call psb_bcast(ctxt,prec%degree2)
     call psb_bcast(ctxt,prec%pvariant2)
+    call psb_bcast(ctxt,prec%prhovariant2)
+    call psb_bcast(ctxt,prec%prhovalue2)  
     call psb_bcast(ctxt,prec%novr2)
     call psb_bcast(ctxt,prec%restr2)
     call psb_bcast(ctxt,prec%prol2)

samples/advanced/pdegen/amg_s_pde3d_poisson_mod.f90

@@ -34,68 +34,68 @@
 !    ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
 !    POSSIBILITY OF SUCH DAMAGE.
 !
-module amg_s_pde3d_base_mod
+module amg_s_pde3d_poisson_mod
   use psb_base_mod, only : psb_spk_, sone, szero
   real(psb_spk_), save, private :: epsilon=sone/80
 contains
-  subroutine pde_set_parm3d_base(dat)
+  subroutine pde_set_parm3d_poisson(dat)
     real(psb_spk_), intent(in) :: dat
     epsilon = dat
-  end subroutine pde_set_parm3d_base
+  end subroutine pde_set_parm3d_poisson
   !
   ! functions parametrizing the differential equation
   !
-  function b1_base(x,y,z)
+  function b1_poisson(x,y,z)
     implicit none 
-    real(psb_spk_) :: b1_base
+    real(psb_spk_) :: b1_poisson
     real(psb_spk_), intent(in) :: x,y,z
-    b1_base=szero/sqrt(3.0_psb_spk_)
-  end function b1_base
-  function b2_base(x,y,z)
+    b1_poisson=szero
+  end function b1_poisson
+  function b2_poisson(x,y,z)
     implicit none 
-    real(psb_spk_) ::  b2_base
+    real(psb_spk_) ::  b2_poisson
     real(psb_spk_), intent(in) :: x,y,z
-    b2_base=szero/sqrt(3.0_psb_spk_)
-  end function b2_base
-  function b3_base(x,y,z)
+    b2_poisson=szero
+  end function b2_poisson
+  function b3_poisson(x,y,z)
     implicit none 
-    real(psb_spk_) ::  b3_base
+    real(psb_spk_) ::  b3_poisson
     real(psb_spk_), intent(in) :: x,y,z
-    b3_base=szero/sqrt(3.0_psb_spk_)
-  end function b3_base
-  function c_base(x,y,z)
+    b3_poisson=szero
+  end function b3_poisson
+  function c_poisson(x,y,z)
     implicit none 
-    real(psb_spk_) ::  c_base
+    real(psb_spk_) ::  c_poisson
     real(psb_spk_), intent(in) :: x,y,z
-    c_base=szero
-  end function c_base
-  function a1_base(x,y,z)
+    c_poisson=szero
+  end function c_poisson
+  function a1_poisson(x,y,z)
     implicit none 
-    real(psb_spk_) ::  a1_base
+    real(psb_spk_) ::  a1_poisson
     real(psb_spk_), intent(in) :: x,y,z
-    a1_base=epsilon
-  end function a1_base
-  function a2_base(x,y,z)
+    a1_poisson=epsilon
+  end function a1_poisson
+  function a2_poisson(x,y,z)
     implicit none 
-    real(psb_spk_) ::  a2_base
+    real(psb_spk_) ::  a2_poisson
     real(psb_spk_), intent(in) :: x,y,z
-    a2_base=epsilon
-  end function a2_base
-  function a3_base(x,y,z)
+    a2_poisson=epsilon
+  end function a2_poisson
+  function a3_poisson(x,y,z)
     implicit none 
-    real(psb_spk_) ::  a3_base
+    real(psb_spk_) ::  a3_poisson
     real(psb_spk_), intent(in) :: x,y,z
-    a3_base=epsilon
-  end function a3_base
-  function g_base(x,y,z)
+    a3_poisson=epsilon
+  end function a3_poisson
+  function g_poisson(x,y,z)
     implicit none 
-    real(psb_spk_) ::  g_base
+    real(psb_spk_) ::  g_poisson
     real(psb_spk_), intent(in) :: x,y,z
-    g_base = szero
+    g_poisson = szero
     if (x == sone) then
-      g_base = sone
+      g_poisson = sone
     else if (x == szero) then
-      g_base = sone
+      g_poisson = sone
     end if
-  end function g_base
-end module amg_s_pde3d_base_mod
+  end function g_poisson
+end module amg_s_pde3d_poisson_mod

samples/advanced/pdegen/runs/amg_pde2d.inp

@@ -1,8 +1,10 @@
 %%%%%%%%%%%  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, GAUSS Coefficients of the PDE
-CG                          ! Iterative method: BiCGSTAB BiCGSTABL BiCG CG CGS FCG GCR RGMRES
+0150                        ! IDIM; domain size. Linear system size is IDIM**2
+POISSON                     ! PDECOEFF: POISSON, EXP, BOX, GAUSS
+%                           ! Coefficients of the PDE
+CG                          ! Iterative method:
+%                           !    BiCGSTAB BiCGSTABL BiCG CG CGS FCG GCR RGMRES
 2                           ! ISTOPC
 00500                       ! ITMAX
 1                           ! ITRACE
@@ -16,11 +18,12 @@ ML                        ! Preconditioner type: NONE JACOBI GS FBGS BJAC AS ML
 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
+CHEB_4_OPT            ! Polynomial variant
 %  Fields to be added for POLY
-%    POLY_RHO_ESTIMATE             Currently only POLY_RHO_EST_POWER
+POLY_RHO_EST_POWER     !    POLY_RHO_ESTIMATE             Currently only POLY_RHO_EST_POWER
 %    POLY_RHO_ESTIMATE_ITERATIONS  default = 20
-%    POLY_RHO_BA                   set to value
+%    POLY_RHO_BA            set to value
+1.0
 %
 0                           ! Number of overlap layers for AS preconditioner
 HALO                        ! AS restriction operator: NONE HALO
@@ -35,7 +38,11 @@ LLK                         ! AINV variant
 NONE                        ! Second (post) smoother, ignored if NONE
 6                           ! Number of sweeps for (post) smoother
 1                           ! degree for polynomial smoother
-POLY_LOTTES_BETA                 ! Polynomial variant
+CHEB_4_OPT                 ! Polynomial variant
+POLY_RHO_EST_POWER     %    POLY_RHO_ESTIMATE             Currently only POLY_RHO_EST_POWER
+%    POLY_RHO_ESTIMATE_ITERATIONS  default = 20
+%    POLY_RHO_BA            set to value
+1.0
 0                           ! Number of overlap layers for AS preconditioner
 HALO                        ! AS restriction operator: NONE HALO
 NONE                        ! AS prolongation operator: NONE SUM AVG

samples/advanced/pdegen/runs/amg_pde3d.inp

@@ -1,8 +1,10 @@
 %%%%%%%%%%%  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, BOX, GAUSS Coefficients of the PDE
-CG                          ! Iterative method: BiCGSTAB BiCGSTABL BiCG CG CGS FCG GCR RGMRES
+0050                        ! IDIM; domain size. Linear system size is IDIM**3
+POISSON                     ! PDECOEFF: POISSON, EXP, BOX, GAUSS
+%                           ! Coefficients of the PDE
+CG                          ! Iterative method:
+%                           !    BiCGSTAB BiCGSTABL BiCG CG CGS FCG GCR RGMRES
 2                           ! ISTOPC
 00500                       ! ITMAX
 1                           ! ITRACE
@@ -13,14 +15,15 @@ ML-VBM-VCYCLE-FBGS-D-BJAC       ! Longer descriptive name for preconditioner (up
 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.
+POLY                        ! 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
+6                           ! degree for polynomial smoother
+CHEB_4_OPT            ! Polynomial variant
 %  Fields to be added for POLY
-%    POLY_RHO_ESTIMATE             Currently only POLY_RHO_EST_POWER
+POLY_RHO_EST_POWER     !    POLY_RHO_ESTIMATE             Currently only POLY_RHO_EST_POWER
 %    POLY_RHO_ESTIMATE_ITERATIONS  default = 20
-%    POLY_RHO_BA                   set to value
+%    POLY_RHO_BA            set to value
+1.0
 %
 0                           ! Number of overlap layers for AS preconditioner
 HALO                        ! AS restriction operator: NONE HALO
@@ -35,7 +38,11 @@ LLK                         ! AINV variant
 NONE                        ! Second (post) smoother, ignored if NONE
 6                           ! Number of sweeps for (post) smoother
 1                           ! degree for polynomial smoother
-POLY_LOTTES_BETA                 ! Polynomial variant
+CHEB_4_OPT                 ! Polynomial variant
+POLY_RHO_EST_POWER     %    POLY_RHO_ESTIMATE             Currently only POLY_RHO_EST_POWER
+%    POLY_RHO_ESTIMATE_ITERATIONS  default = 20
+%    POLY_RHO_BA            set to value
+1.0
 0                           ! Number of overlap layers for AS preconditioner
 HALO                        ! AS restriction operator: NONE HALO
 NONE                        ! AS prolongation operator: NONE SUM AVG