Merge branch 'repackage' of github.com:sfilippone/amg4psblas into repackage

.gitignore

@@ -20,3 +20,5 @@ autom4te.cache
 # the executable from tests
 runs
 
+# Documentation temporary files
+docs/src/userguide.pdf

Make_n.inc.in

@@ -49,7 +49,7 @@ PSBLAS_INCLUDES=@PSBLAS_INCLUDES@
 PSBLAS_LIBS=@PSBLAS_LIBS@	
 PSBBASEMODNAME=psb_base_mod
 PSBPRECMODNAME=psb_prec_mod
-PSBMETHDMODNAME=psb_krylov_mod
+PSBMETHDMODNAME=psb_linsolve_mod
 PSBUTILMODNAME=psb_util_mod
 
 

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
 
@@ -46,6 +46,7 @@ cleanlib:
 
 veryclean: cleanlib
 	(cd amgprec && $(MAKE) veryclean)
+	(cd cbind && $(MAKE) veryclean)
 	(cd samples/simple/fileread && $(MAKE) clean)
 	(cd samples/simple/pdegen && $(MAKE) clean)
 	(cd samples/advanced/fileread && $(MAKE) clean)
@@ -56,3 +57,4 @@ check: all
 
 clean:
 	(cd amgprec && $(MAKE) clean)
+	(cd cbind && $(MAKE) clean)

README.md

@@ -1,54 +1,76 @@
-                         AMG4PSBLAS  
+# AMG4PSBLAS  v1.2
- Algebraic Multigrid Package  based on PSBLAS (Parallel Sparse BLAS version 3.8)
+Algebraic Multigrid Package  based on [PSBLAS](https://github.com/sfilippone/psblas3) (Parallel Sparse BLAS version 3.9)
-    
-Salvatore Filippone    (University of Rome Tor Vergata and IAC-CNR)
-Pasqua D'Ambra         (IAC-CNR, Naples, IT)
-Fabio Durastante       (IAC-CNR, Naples, IT)
 
----------------------------------------------------------------------
+AMG4PSBLAS is a package of parallel algebraic multilevel preconditioners included in the PSCToolkit (Parallel Sparse Computation Toolkit) software framework.
 
-AMG4PSBLAS is a package of Algebraic MultiGrid (AMG)
+It is a progress of a software development project started in 2007, named MLD2P4, which originally implemented a multilevel version of some domain decomposition preconditioners of additive-Schwarz type and was based on a parallel decoupled version of the well known smoothed aggregation method to generate the multilevel hierarchy of coarser matrices.
-preconditioners for  the iterative solution of large and sparse linear systems.
 
-It is an evolution of MLD2P4 (see LICENSE.MLD2P4), but it has been
+In the last years the package was extended for including new algorithms and functionalities for the setup and application new AMG preconditioners with the final aims of improving efficiency and scalability when tens of thousands cores are used and of boosting reliability in dealing with general symmetric positive definite linear systems.
-thoroughly reworked, and it is sufficiently different to warrant a new
-project name.
 
+It is an evolution of MLD2P4 (see [LICENSE.MLD2P4](LICENSE.MLD2P4)), but due to the significant number of changes and the increase in scope, we decided to rename the package as AMG4PSBLAS.
 
-MAIN REFERENCES:
+AMG4PSBLAS has been designed to provide scalable and easy-to-use preconditioners in the context of the PSBLAS (Parallel Sparse Basic Linear Algebra Subprograms) computational framework and can be used in conjuction with the Krylov solvers available in this framework. Our package is based on a completely algebraic approach; therefore users level interfaces assume that the system matrix and preconditioners are represented as PSBLAS distributed sparse matrices.
 
-     
+AMG4PSBLAS enables the user to easily specify different features of an algebraic multilevel preconditioner, thus allowing to experiment with different preconditioners for the problem and parallel computers at hand.
 
-P. D'Ambra, D. di Serafino, S. Filippone,
+The package employs object-oriented design techniques in Fortran 2008, with interfaces to additional third party libraries such as MUMPS, UMFPACK, SuperLU, and SuperLU_Dist, which can be exploited in building multilevel preconditioners. The parallel implementation is based on a Single Program Multiple Data (SPMD) paradigm; the inter-process communication is based on MPI and is managed mainly through PSBLAS.
-MLD2P4: a Package of Parallel Algebraic Multilevel Domain Decomposition
-Preconditioners in Fortran 95,
-ACM Transactions on Mathematical Software, 37 (3), 2010, art. 30,
-doi: 10.1145/1824801.1824808.
 
+## Main Refrerences:
 
-TO COMPILE
+The main reference for this project is
+> D'Ambra, P., Durastante, F., & Filippone, S. (2021). AMG preconditioners for linear solvers towards extreme scale. SIAM Journal on Scientific Computing, 43(5), S679-S703.
+
+AMG4PSBLAS is the suite of preconditioners for the Parallel Sparse Computation Toolkit ([PSCToolkit](https://psctoolkit.github.io/)) suite of libraries. See the paper:
+> D’Ambra, P., Durastante, F., & Filippone, S. (2023). Parallel Sparse Computation Toolkit. Software Impacts, 15, 100463.
+
+The main reference for features inherited from MLD2P4 is
+> P. D'Ambra, D. di Serafino, S. Filippone,
+> MLD2P4: a Package of Parallel Algebraic Multilevel Domain Decomposition
+> Preconditioners in Fortran 95,
+> ACM Transactions on Mathematical Software, 37 (3), 2010, art. 30,
+> doi: 10.1145/1824801.1824808.
+
+## Installing
+
+Installation requires having a working version of the [PSBLAS](https://github.com/sfilippone/psblas3) library installed.
+AMG4PSBLAS has several interfaces to third-party libraries that can be used in the construction and application phases of preconditioners. 
+In particular, it is possible to link AMG4PSBLAS with the libraries: MUMPS, SuperLU, SuperLU_Dist, UMFPACK. This is _not mandatory_ and the library can run 
+in isolation and without these features.
 
 0. Unpack the tar file in a directory of your choice (preferrably
    outside the main PSBLAS directory).
-1. run configure --with-psblas=<ABSOLUTE path of the PSBLAS install directory>
+1. run configure `--with-psblas=<ABSOLUTE path of the PSBLAS install directory>`
    adding the options for MUMPS, SuperLU, SuperLU_Dist, UMFPACK as desired.
-   See MLD2P4 User's and Reference Guide (Section 3) for details.
+   See [AMG4PSBLAS User's and Reference Guide](docs/amg4psblas_1.0-guide.pdf) (Section 3) for details.
-2. Tweak Make.inc if you are not satisfied.
+2. Tweak `Make.inc` if you are not satisfied.
-3. make; 
+3. run `make`; 
 4. Go into the test subdirectory and build the examples of your choice.
-5. (if desired): make install 
+5. (if desired): `make install` 
 
+>[!CAUTION]
+>The single precision version is supported only by MUMPS and SuperLU;
+>thus, even if you specify at configure time to use UMFPACK or SuperLU_Dist,
+>the corresponding preconditioner options will be available only from
+>the double precision version.
 
-NOTES
+### CUDA, OpeMP, OpenACC
 
-- The single precision version is supported only by MUMPS and SuperLU;
+CUDA, OpenMP and OpenACC features are transparently inherited by PSBLAS installation. If PSBLAS has been configured (and installed) with these supports then AMG4PSBLAS will transparently inherit them. It will then be possible to move the computation to GPU accelerator simply by selecting the appropriate variable types. If these have not been activated or installed for PSBLAS then they will not be available for AMG4PSBLAS either and the operation will be purely on CPU/MPI.
-  thus, even if you specify at configure time to use UMFPACK or SuperLU_Dist, 
-  the corresponding preconditioner options will be available only from
-  the double precision version.
 
+### EoCoE - Software as service portal
+
+In the European project “Energy oriented Center of Excellence: toward exascale for energy” we made available a software as service portal: [https://eocoe.psnc.pl/](https://eocoe.psnc.pl/). This permits to test several cutting-edge computational methods for accelerating the transition to the production, storage and management of clean, decarbonized energy. Among them you have the possibility of running PSBLAS+AMG4PSBLAS on some test problems to become familiar with using the software.
+
+## TODO and bugs
+ 
+- [X] Fix all reamining bugs. Bugs? We dont' have any ! 🤓
+
+> [!NOTE]
+> To report bugs 🐛 or issues ❓ please use the [GitHub issue system](https://github.com/sfilippone/amg4psblas/issues).
+
+## The AMG4PSBLAS team. 
+
+- Pasqua D'Ambra         (IAC-CNR, Naples, IT)
+- Fabio Durastante       (University of Pisa and IAC-CNR, IT)
+- Salvatore Filippone    (University of Rome Tor Vergata and IAC-CNR, IT)
 
-The AMG4PSBLAS team. 
----------------
-Salvatore Filippone
-Pasqua     D'Ambra
-Fabio     Durastante

amgprec/amg_base_prec_type.F90

@@ -288,17 +288,19 @@ 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_no_smooth_      = 0
-  integer(psb_ipk_), parameter :: amg_smooth_prol_ = 1
+  integer(psb_ipk_), parameter :: amg_smooth_prol_    = 1
-  integer(psb_ipk_), parameter :: amg_min_energy_  = 2
+  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_
   !
-  integer(psb_ipk_), parameter :: amg_no_filter_mat_  = 0
+  integer(psb_ipk_), parameter :: amg_no_filter_mat_   = 0
-  integer(psb_ipk_), parameter :: amg_filter_mat_     = 1
+  integer(psb_ipk_), parameter :: amg_filter_mat_      = 1
-  integer(psb_ipk_), parameter :: amg_max_filter_mat_ = amg_filter_mat_
+  integer(psb_ipk_), parameter :: amg_filter_prow_mat_ = 2
+  integer(psb_ipk_), parameter :: amg_max_filter_mat_  = amg_filter_prow_mat_
   !
   ! Legal values for entry: amg_aggr_ord_
   !
@@ -323,10 +325,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_cheb_4_     = 0
-  integer(psb_ipk_), parameter :: amg_poly_lottes_beta_ = 1
+  integer(psb_ipk_), parameter :: amg_cheb_4_opt_ = 1
-  integer(psb_ipk_), parameter :: amg_poly_new_         = 2
+  integer(psb_ipk_), parameter :: amg_cheb_1_opt_ = 2
-  integer(psb_ipk_), parameter :: amg_poly_dbg_         = 8
+  integer(psb_ipk_), parameter :: amg_poly_dbg_   = 8
   
   integer(psb_ipk_), parameter :: amg_poly_rho_est_power_ = 0
 
@@ -376,10 +378,11 @@ 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      ',&
+       &  aggr_prols(0:4)=(/'unsmoothed    ','smoothed      ',&
-       &           'min energy    ','bizr. smoothed'/)
+       &  'l1-smoothed   ','min energy    ','bizr. smoothed'/)
   character(len=15), parameter, private :: &
-       &  aggr_filters(0:1)=(/'no filtering  ','filtering     '/)
+       &  aggr_filters(0:2)=(/'no filtering  ','filtering     ',&
+       &  'filtering rsum'/)
   character(len=15), parameter, private :: &
        &  matrix_names(0:1)=(/'distributed   ','replicated    '/)
   character(len=18), parameter, private :: &
@@ -548,6 +551,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 +575,12 @@ contains
       val = amg_as_
     case('POLY')
       val = amg_poly_
-    case('POLY_LOTTES')
+    case('CHEB_4')
-      val = amg_poly_lottes_
+      val = amg_cheb_4_
-    case('POLY_LOTTES_BETA')
+    case('CHEB_4_OPT')
-      val = amg_poly_lottes_beta_
+      val = amg_cheb_4_opt_
-    case('POLY_NEW')
+    case('CHEB_1_OPT')
-      val = amg_poly_new_
+      val = amg_cheb_1_opt_
     case('POLY_DBG')
       val = amg_poly_dbg_
     case('POLY_RHO_EST_POWER')
@@ -588,6 +593,8 @@ contains
       val = amg_eig_est_
     case('FILTER')
       val = amg_filter_mat_
+    case('FILTERROWSUM')
+      val = amg_filter_prow_mat_
     case('NOFILTER','NO_FILTER')
       val = amg_no_filter_mat_
     case('OUTER_SWEEPS')

amgprec/amg_c_ainv_solver.F90

@@ -91,9 +91,9 @@ module amg_c_ainv_solver
       import :: psb_desc_type, psb_cspmat_type,  psb_c_base_sparse_mat, &
        & amg_c_base_solver_type, psb_dpk_, amg_c_ainv_solver_type, psb_ipk_
       Implicit None
-      class(amg_c_ainv_solver_type), intent(inout)              :: sv
+      class(amg_c_ainv_solver_type), intent(inout) :: sv
-      class(amg_c_base_solver_type), allocatable, intent(inout) :: svout
+      class(amg_c_base_solver_type), intent(inout) :: svout
-      integer(psb_ipk_), intent(out)                            :: info
+      integer(psb_ipk_), intent(out)                 :: info
     end subroutine amg_c_ainv_solver_clone_settings
   end interface
 

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_c_invk_solver.f90

@@ -79,9 +79,9 @@ module amg_c_invk_solver
       import :: psb_desc_type, psb_cspmat_type,  psb_c_base_sparse_mat, &
        & amg_c_base_solver_type, psb_spk_, amg_c_invk_solver_type, psb_ipk_
       Implicit None
-      class(amg_c_invk_solver_type), intent(inout)              :: sv
+      class(amg_c_invk_solver_type), intent(inout) :: sv
-      class(amg_c_base_solver_type), allocatable, intent(inout) :: svout
+      class(amg_c_base_solver_type), intent(inout) :: svout
-      integer(psb_ipk_), intent(out)                            :: info
+      integer(psb_ipk_), intent(out)                 :: info
     end subroutine amg_c_invk_solver_clone_settings
   end interface
 

amgprec/amg_c_invt_solver.f90

@@ -79,9 +79,9 @@ module amg_c_invt_solver
       import :: psb_desc_type, psb_cspmat_type,  psb_c_base_sparse_mat, &
        & amg_c_base_solver_type, psb_spk_, amg_c_invt_solver_type, psb_ipk_
       Implicit None
-      class(amg_c_invt_solver_type), intent(inout)              :: sv
+      class(amg_c_invt_solver_type), intent(inout) :: sv
-      class(amg_c_base_solver_type), allocatable, intent(inout) :: svout
+      class(amg_c_base_solver_type), intent(inout) :: svout
-      integer(psb_ipk_), intent(out)                            :: info
+      integer(psb_ipk_), intent(out)                 :: info
     end subroutine amg_c_invt_solver_clone_settings
   end interface
 

amgprec/amg_c_jac_smoother.f90

@@ -203,8 +203,8 @@ module amg_c_jac_smoother
     subroutine amg_c_jac_smoother_clone_settings(sm,smout,info)
       import :: amg_c_jac_smoother_type, psb_spk_, &
            & amg_c_base_smoother_type, psb_ipk_
-      class(amg_c_jac_smoother_type), intent(inout)               :: sm
+      class(amg_c_jac_smoother_type), intent(inout)  :: sm
-      class(amg_c_base_smoother_type), allocatable, intent(inout) :: smout
+      class(amg_c_base_smoother_type), intent(inout) :: smout
       integer(psb_ipk_), intent(out)                :: info
     end subroutine amg_c_jac_smoother_clone_settings
   end interface

amgprec/amg_d_ainv_solver.F90

@@ -91,9 +91,9 @@ module amg_d_ainv_solver
       import :: psb_desc_type, psb_dspmat_type,  psb_d_base_sparse_mat, &
        & amg_d_base_solver_type, psb_dpk_, amg_d_ainv_solver_type, psb_ipk_
       Implicit None
-      class(amg_d_ainv_solver_type), intent(inout)              :: sv
+      class(amg_d_ainv_solver_type), intent(inout) :: sv
-      class(amg_d_base_solver_type), allocatable, intent(inout) :: svout
+      class(amg_d_base_solver_type), intent(inout) :: svout
-      integer(psb_ipk_), intent(out)                            :: info
+      integer(psb_ipk_), intent(out)                 :: info
     end subroutine amg_d_ainv_solver_clone_settings
   end interface
 

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_invk_solver.f90

@@ -79,9 +79,9 @@ module amg_d_invk_solver
       import :: psb_desc_type, psb_dspmat_type,  psb_d_base_sparse_mat, &
        & amg_d_base_solver_type, psb_dpk_, amg_d_invk_solver_type, psb_ipk_
       Implicit None
-      class(amg_d_invk_solver_type), intent(inout)              :: sv
+      class(amg_d_invk_solver_type), intent(inout) :: sv
-      class(amg_d_base_solver_type), allocatable, intent(inout) :: svout
+      class(amg_d_base_solver_type), intent(inout) :: svout
-      integer(psb_ipk_), intent(out)                            :: info
+      integer(psb_ipk_), intent(out)                 :: info
     end subroutine amg_d_invk_solver_clone_settings
   end interface
 

amgprec/amg_d_invt_solver.f90

@@ -79,9 +79,9 @@ module amg_d_invt_solver
       import :: psb_desc_type, psb_dspmat_type,  psb_d_base_sparse_mat, &
        & amg_d_base_solver_type, psb_dpk_, amg_d_invt_solver_type, psb_ipk_
       Implicit None
-      class(amg_d_invt_solver_type), intent(inout)              :: sv
+      class(amg_d_invt_solver_type), intent(inout) :: sv
-      class(amg_d_base_solver_type), allocatable, intent(inout) :: svout
+      class(amg_d_base_solver_type), intent(inout) :: svout
-      integer(psb_ipk_), intent(out)                            :: info
+      integer(psb_ipk_), intent(out)                 :: info
     end subroutine amg_d_invt_solver_clone_settings
   end interface
 

amgprec/amg_d_jac_smoother.f90

@@ -203,8 +203,8 @@ module amg_d_jac_smoother
     subroutine amg_d_jac_smoother_clone_settings(sm,smout,info)
       import :: amg_d_jac_smoother_type, psb_dpk_, &
            & amg_d_base_smoother_type, psb_ipk_
-      class(amg_d_jac_smoother_type), intent(inout)               :: sm
+      class(amg_d_jac_smoother_type), intent(inout)  :: sm
-      class(amg_d_base_smoother_type), allocatable, intent(inout) :: smout
+      class(amg_d_base_smoother_type), intent(inout) :: smout
       integer(psb_ipk_), intent(out)                :: info
     end subroutine amg_d_jac_smoother_clone_settings
   end interface

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

@@ -192,8 +192,8 @@ module amg_d_poly_smoother
     subroutine amg_d_poly_smoother_clone_settings(sm,smout,info)
       import :: amg_d_poly_smoother_type, psb_dpk_, &
            & amg_d_base_smoother_type, psb_ipk_
-      class(amg_d_poly_smoother_type), intent(inout)               :: sm
+      class(amg_d_poly_smoother_type), intent(inout) :: sm
-      class(amg_d_base_smoother_type), allocatable, intent(inout) :: smout
+      class(amg_d_base_smoother_type), intent(inout) :: smout
       integer(psb_ipk_), intent(out)                :: info
     end subroutine amg_d_poly_smoother_clone_settings
   end interface
@@ -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_ainv_solver.F90

@@ -91,9 +91,9 @@ module amg_s_ainv_solver
       import :: psb_desc_type, psb_sspmat_type,  psb_s_base_sparse_mat, &
        & amg_s_base_solver_type, psb_dpk_, amg_s_ainv_solver_type, psb_ipk_
       Implicit None
-      class(amg_s_ainv_solver_type), intent(inout)              :: sv
+      class(amg_s_ainv_solver_type), intent(inout) :: sv
-      class(amg_s_base_solver_type), allocatable, intent(inout) :: svout
+      class(amg_s_base_solver_type), intent(inout) :: svout
-      integer(psb_ipk_), intent(out)                            :: info
+      integer(psb_ipk_), intent(out)                 :: info
     end subroutine amg_s_ainv_solver_clone_settings
   end interface
 

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_invk_solver.f90

@@ -79,9 +79,9 @@ module amg_s_invk_solver
       import :: psb_desc_type, psb_sspmat_type,  psb_s_base_sparse_mat, &
        & amg_s_base_solver_type, psb_spk_, amg_s_invk_solver_type, psb_ipk_
       Implicit None
-      class(amg_s_invk_solver_type), intent(inout)              :: sv
+      class(amg_s_invk_solver_type), intent(inout) :: sv
-      class(amg_s_base_solver_type), allocatable, intent(inout) :: svout
+      class(amg_s_base_solver_type), intent(inout) :: svout
-      integer(psb_ipk_), intent(out)                            :: info
+      integer(psb_ipk_), intent(out)                 :: info
     end subroutine amg_s_invk_solver_clone_settings
   end interface
 

amgprec/amg_s_invt_solver.f90

@@ -79,9 +79,9 @@ module amg_s_invt_solver
       import :: psb_desc_type, psb_sspmat_type,  psb_s_base_sparse_mat, &
        & amg_s_base_solver_type, psb_spk_, amg_s_invt_solver_type, psb_ipk_
       Implicit None
-      class(amg_s_invt_solver_type), intent(inout)              :: sv
+      class(amg_s_invt_solver_type), intent(inout) :: sv
-      class(amg_s_base_solver_type), allocatable, intent(inout) :: svout
+      class(amg_s_base_solver_type), intent(inout) :: svout
-      integer(psb_ipk_), intent(out)                            :: info
+      integer(psb_ipk_), intent(out)                 :: info
     end subroutine amg_s_invt_solver_clone_settings
   end interface
 

amgprec/amg_s_jac_smoother.f90

@@ -203,8 +203,8 @@ module amg_s_jac_smoother
     subroutine amg_s_jac_smoother_clone_settings(sm,smout,info)
       import :: amg_s_jac_smoother_type, psb_spk_, &
            & amg_s_base_smoother_type, psb_ipk_
-      class(amg_s_jac_smoother_type), intent(inout)               :: sm
+      class(amg_s_jac_smoother_type), intent(inout)  :: sm
-      class(amg_s_base_smoother_type), allocatable, intent(inout) :: smout
+      class(amg_s_base_smoother_type), intent(inout) :: smout
       integer(psb_ipk_), intent(out)                :: info
     end subroutine amg_s_jac_smoother_clone_settings
   end interface

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

@@ -192,8 +192,8 @@ module amg_s_poly_smoother
     subroutine amg_s_poly_smoother_clone_settings(sm,smout,info)
       import :: amg_s_poly_smoother_type, psb_spk_, &
            & amg_s_base_smoother_type, psb_ipk_
-      class(amg_s_poly_smoother_type), intent(inout)               :: sm
+      class(amg_s_poly_smoother_type), intent(inout) :: sm
-      class(amg_s_base_smoother_type), allocatable, intent(inout) :: smout
+      class(amg_s_base_smoother_type), intent(inout) :: smout
       integer(psb_ipk_), intent(out)                :: info
     end subroutine amg_s_poly_smoother_clone_settings
   end interface
@@ -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_ainv_solver.F90

@@ -91,9 +91,9 @@ module amg_z_ainv_solver
       import :: psb_desc_type, psb_zspmat_type,  psb_z_base_sparse_mat, &
        & amg_z_base_solver_type, psb_dpk_, amg_z_ainv_solver_type, psb_ipk_
       Implicit None
-      class(amg_z_ainv_solver_type), intent(inout)              :: sv
+      class(amg_z_ainv_solver_type), intent(inout) :: sv
-      class(amg_z_base_solver_type), allocatable, intent(inout) :: svout
+      class(amg_z_base_solver_type), intent(inout) :: svout
-      integer(psb_ipk_), intent(out)                            :: info
+      integer(psb_ipk_), intent(out)                 :: info
     end subroutine amg_z_ainv_solver_clone_settings
   end interface
 

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/amg_z_invk_solver.f90

@@ -79,9 +79,9 @@ module amg_z_invk_solver
       import :: psb_desc_type, psb_zspmat_type,  psb_z_base_sparse_mat, &
        & amg_z_base_solver_type, psb_dpk_, amg_z_invk_solver_type, psb_ipk_
       Implicit None
-      class(amg_z_invk_solver_type), intent(inout)              :: sv
+      class(amg_z_invk_solver_type), intent(inout) :: sv
-      class(amg_z_base_solver_type), allocatable, intent(inout) :: svout
+      class(amg_z_base_solver_type), intent(inout) :: svout
-      integer(psb_ipk_), intent(out)                            :: info
+      integer(psb_ipk_), intent(out)                 :: info
     end subroutine amg_z_invk_solver_clone_settings
   end interface
 

amgprec/amg_z_invt_solver.f90

@@ -79,9 +79,9 @@ module amg_z_invt_solver
       import :: psb_desc_type, psb_zspmat_type,  psb_z_base_sparse_mat, &
        & amg_z_base_solver_type, psb_dpk_, amg_z_invt_solver_type, psb_ipk_
       Implicit None
-      class(amg_z_invt_solver_type), intent(inout)              :: sv
+      class(amg_z_invt_solver_type), intent(inout) :: sv
-      class(amg_z_base_solver_type), allocatable, intent(inout) :: svout
+      class(amg_z_base_solver_type), intent(inout) :: svout
-      integer(psb_ipk_), intent(out)                            :: info
+      integer(psb_ipk_), intent(out)                 :: info
     end subroutine amg_z_invt_solver_clone_settings
   end interface
 

amgprec/amg_z_jac_smoother.f90

@@ -203,8 +203,8 @@ module amg_z_jac_smoother
     subroutine amg_z_jac_smoother_clone_settings(sm,smout,info)
       import :: amg_z_jac_smoother_type, psb_dpk_, &
            & amg_z_base_smoother_type, psb_ipk_
-      class(amg_z_jac_smoother_type), intent(inout)               :: sm
+      class(amg_z_jac_smoother_type), intent(inout)  :: sm
-      class(amg_z_base_smoother_type), allocatable, intent(inout) :: smout
+      class(amg_z_base_smoother_type), intent(inout) :: smout
       integer(psb_ipk_), intent(out)                :: info
     end subroutine amg_z_jac_smoother_clone_settings
   end interface

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,&
+    call amg_caggrmat_nosmth_bld(parms%aggr_prol,a,desc_a,ilaggr,&
-         & parms,ac,desc_ac,op_prol,op_restr,t_prol,info)
+            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, &
+    call amg_caggrmat_smth_bld(parms%aggr_prol,a,desc_a,&
-         & parms,ac,desc_ac,op_prol,op_restr,t_prol,info)
+            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, &
+    call amg_caggrmat_minnrg_bld(parms%aggr_prol,a,desc_a,ilaggr,&
-         & parms,ac,desc_ac,op_prol,op_restr,t_prol,info)
+            nlaggr,parms,ac,desc_ac,op_prol,op_restr,t_prol,info)
 
   case default
     info = psb_err_internal_error_

amgprec/impl/aggregator/amg_c_ptap_bld.f90

@@ -216,6 +216,7 @@ subroutine amg_c_ptap_bld(a_csr,desc_a,nlaggr,parms,ac,&
 
   call ac_csr%set_nrows(desc_ac%get_local_rows())
   call ac_csr%set_ncols(desc_ac%get_local_cols())
+  call ac_csr%clean_zeros(info)
   call ac%mv_from(ac_csr)
   call ac%set_asb()
 
@@ -420,6 +421,7 @@ subroutine amg_c_lc_ptap_bld(a_csr,desc_a,nlaggr,parms,ac,&
 
   call ac_csr%set_nrows(desc_ac%get_local_rows())
   call ac_csr%set_ncols(desc_ac%get_local_cols())
+  call ac_csr%clean_zeros(info)
   call ac%mv_from(ac_csr)
   call ac%set_asb()
 
@@ -629,6 +631,7 @@ subroutine amg_lc_ptap_bld(a_csr,desc_a,nlaggr,parms,ac,&
 
   call ac_csr%set_nrows(desc_ac%get_local_rows())
   call ac_csr%set_ncols(desc_ac%get_local_cols())
+  call ac_csr%clean_zeros(info)
   call ac%mv_from(ac_csr)
   call ac%set_asb()
 

amgprec/impl/aggregator/amg_c_rap.f90

@@ -142,6 +142,7 @@ subroutine amg_c_rap(a_csr,desc_a,nlaggr,parms,ac,&
 
   call ac_csr%set_nrows(desc_ac%get_local_rows())
   call ac_csr%set_ncols(desc_ac%get_local_cols())
+  call ac_csr%clean_zeros(info)
   call ac%mv_from(ac_csr)
   call ac%set_asb()
 

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,&
+subroutine amg_caggrmat_minnrg_bld(dol1smoothing,a,desc_a,ilaggr,nlaggr,&
-     & ac,desc_ac,op_prol,op_restr,t_prol,info)
+                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(dol1smoothing,a,desc_a,ilaggr,nlaggr,&
-subroutine amg_caggrmat_nosmth_bld(a,desc_a,ilaggr,nlaggr,parms,&
+                parms,ac,desc_ac,op_prol,op_restr,t_prol,info)
-     & 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,16 +104,18 @@
 !    info       -  integer, output.
 !                  Error code.
 !
-subroutine amg_caggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,&
+subroutine amg_caggrmat_smth_bld(dol1smoothing,a,desc_a,ilaggr,nlaggr,&
-     & ac,desc_ac,op_prol,op_restr,t_prol,info)
+                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
   use amg_c_base_aggregator_mod
+!  use, intrinsic :: ieee_arithmetic
 
   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 +136,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 +145,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 +178,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()
@@ -185,7 +193,7 @@ subroutine amg_caggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,&
 
   naggrm1 = sum(nlaggr(1:me))
   naggrp1 = sum(nlaggr(1:me+1))
-  filter_mat = (parms%aggr_filter == amg_filter_mat_)
+  filter_mat = (parms%aggr_filter == amg_filter_mat_).or.(parms%aggr_filter == amg_filter_prow_mat_)
 
   !
   ! naggr: number of local aggregates
@@ -200,6 +208,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,9 +256,15 @@ 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
+      else if (parms%aggr_filter == amg_filter_mat_) then
+        ! We perform filtering in the standard way assuming that A is an M-matrix
         acsrf%val(jd)=acsrf%val(jd)-tmp
+      else if (parms%aggr_filter == amg_filter_prow_mat_) then
+        ! We are probably doing l1-correction, hence we want to preserve the
+        ! row sum of the matrix: note the change in sign
+        acsrf%val(jd)=acsrf%val(jd)+tmp
       end if
     enddo
     !$OMP end parallel do 
@@ -240,7 +272,6 @@ subroutine amg_caggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,&
     call acsrf%clean_zeros(info)
   end if
 
-
   !$OMP parallel do private(i) schedule(static)
   do i=1,size(adiag)
     if (adiag(i) /= czero) then
@@ -252,14 +283,17 @@ 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 (parms%aggr_eig == amg_max_norm_) then 
+    if ( (parms%aggr_filter == amg_filter_prow_mat_).and.(do_l1correction) ) then
+      ! For l1-Jacobi this can be estimated with 1:
+      ! this makes sense only if we are preserving the row-sum!
+      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)))
       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 +356,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,&
+    call amg_daggrmat_nosmth_bld(parms%aggr_prol,a,desc_a,ilaggr,&
-         & parms,ac,desc_ac,op_prol,op_restr,t_prol,info)
+            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, &
+    call amg_daggrmat_smth_bld(parms%aggr_prol,a,desc_a,&
-         & parms,ac,desc_ac,op_prol,op_restr,t_prol,info)
+            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, &
+    call amg_daggrmat_minnrg_bld(parms%aggr_prol,a,desc_a,ilaggr,&
-         & parms,ac,desc_ac,op_prol,op_restr,t_prol,info)
+            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,&
+    call amg_d_parmatch_unsmth_bld(parms%aggr_prol,ag,a,desc_a,&
-         & ac,desc_ac,op_prol,op_restr,t_prol,info)
+            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_d_parmatch_smth_bld(ag,a,desc_a,ilaggr,nlaggr,parms,&
+    call amg_d_parmatch_smth_bld(parms%aggr_prol,ag,a,desc_a,&
-         & ac,desc_ac,op_prol,op_restr,t_prol,info)
+            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, &
+    call amg_daggrmat_minnrg_bld(parms%aggr_prol,a,desc_a,&
-         & parms,ac,desc_ac,op_prol,op_restr,t_prol,info)
+            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,&
+subroutine amg_d_parmatch_smth_bld(dol1smoothing,ag,a,desc_a,ilaggr,nlaggr,&
-     & ac,desc_ac,op_prol,op_restr,t_prol,info)
+                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,&
+subroutine amg_d_parmatch_unsmth_bld(dol1smoothing,ag,a,desc_a,ilaggr,nlaggr,&
-     & ac,desc_ac,op_prol,op_restr,t_prol,info)
+                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_ptap_bld.f90

@@ -216,6 +216,7 @@ subroutine amg_d_ptap_bld(a_csr,desc_a,nlaggr,parms,ac,&
 
   call ac_csr%set_nrows(desc_ac%get_local_rows())
   call ac_csr%set_ncols(desc_ac%get_local_cols())
+  call ac_csr%clean_zeros(info)
   call ac%mv_from(ac_csr)
   call ac%set_asb()
 
@@ -420,6 +421,7 @@ subroutine amg_d_ld_ptap_bld(a_csr,desc_a,nlaggr,parms,ac,&
 
   call ac_csr%set_nrows(desc_ac%get_local_rows())
   call ac_csr%set_ncols(desc_ac%get_local_cols())
+  call ac_csr%clean_zeros(info)
   call ac%mv_from(ac_csr)
   call ac%set_asb()
 
@@ -629,6 +631,7 @@ subroutine amg_ld_ptap_bld(a_csr,desc_a,nlaggr,parms,ac,&
 
   call ac_csr%set_nrows(desc_ac%get_local_rows())
   call ac_csr%set_ncols(desc_ac%get_local_cols())
+  call ac_csr%clean_zeros(info)
   call ac%mv_from(ac_csr)
   call ac%set_asb()
 

amgprec/impl/aggregator/amg_d_rap.f90

@@ -142,6 +142,7 @@ subroutine amg_d_rap(a_csr,desc_a,nlaggr,parms,ac,&
 
   call ac_csr%set_nrows(desc_ac%get_local_rows())
   call ac_csr%set_ncols(desc_ac%get_local_cols())
+  call ac_csr%clean_zeros(info)
   call ac%mv_from(ac_csr)
   call ac%set_asb()
 

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,&
+subroutine amg_daggrmat_minnrg_bld(dol1smoothing,a,desc_a,ilaggr,nlaggr,&
-     & ac,desc_ac,op_prol,op_restr,t_prol,info)
+                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(dol1smoothing,a,desc_a,ilaggr,nlaggr,&
-subroutine amg_daggrmat_nosmth_bld(a,desc_a,ilaggr,nlaggr,parms,&
+                parms,ac,desc_ac,op_prol,op_restr,t_prol,info)
-     & 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,16 +104,18 @@
 !    info       -  integer, output.
 !                  Error code.
 !
-subroutine amg_daggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,&
+subroutine amg_daggrmat_smth_bld(dol1smoothing,a,desc_a,ilaggr,nlaggr,&
-     & ac,desc_ac,op_prol,op_restr,t_prol,info)
+                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
   use amg_d_base_aggregator_mod
+!  use, intrinsic :: ieee_arithmetic
 
   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 +136,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 +145,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 +178,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()
@@ -185,7 +193,7 @@ subroutine amg_daggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,&
 
   naggrm1 = sum(nlaggr(1:me))
   naggrp1 = sum(nlaggr(1:me+1))
-  filter_mat = (parms%aggr_filter == amg_filter_mat_)
+  filter_mat = (parms%aggr_filter == amg_filter_mat_).or.(parms%aggr_filter == amg_filter_prow_mat_)
 
   !
   ! naggr: number of local aggregates
@@ -200,6 +208,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,9 +256,15 @@ 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
+      else if (parms%aggr_filter == amg_filter_mat_) then
+        ! We perform filtering in the standard way assuming that A is an M-matrix
         acsrf%val(jd)=acsrf%val(jd)-tmp
+      else if (parms%aggr_filter == amg_filter_prow_mat_) then
+        ! We are probably doing l1-correction, hence we want to preserve the
+        ! row sum of the matrix: note the change in sign
+        acsrf%val(jd)=acsrf%val(jd)+tmp
       end if
     enddo
     !$OMP end parallel do 
@@ -240,7 +272,6 @@ subroutine amg_daggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,&
     call acsrf%clean_zeros(info)
   end if
 
-
   !$OMP parallel do private(i) schedule(static)
   do i=1,size(adiag)
     if (adiag(i) /= dzero) then
@@ -252,14 +283,17 @@ 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 (parms%aggr_eig == amg_max_norm_) then 
+    if ( (parms%aggr_filter == amg_filter_prow_mat_).and.(do_l1correction) ) then
+      ! For l1-Jacobi this can be estimated with 1:
+      ! this makes sense only if we are preserving the row-sum!
+      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)))
       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 +356,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,&
+    call amg_saggrmat_nosmth_bld(parms%aggr_prol,a,desc_a,ilaggr,&
-         & parms,ac,desc_ac,op_prol,op_restr,t_prol,info)
+            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, &
+    call amg_saggrmat_smth_bld(parms%aggr_prol,a,desc_a,&
-         & parms,ac,desc_ac,op_prol,op_restr,t_prol,info)
+            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, &
+    call amg_saggrmat_minnrg_bld(parms%aggr_prol,a,desc_a,ilaggr,&
-         & parms,ac,desc_ac,op_prol,op_restr,t_prol,info)
+            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,&
+    call amg_s_parmatch_unsmth_bld(parms%aggr_prol,ag,a,desc_a,&
-         & ac,desc_ac,op_prol,op_restr,t_prol,info)
+            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_s_parmatch_smth_bld(ag,a,desc_a,ilaggr,nlaggr,parms,&
+    call amg_s_parmatch_smth_bld(parms%aggr_prol,ag,a,desc_a,&
-         & ac,desc_ac,op_prol,op_restr,t_prol,info)
+            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, &
+    call amg_saggrmat_minnrg_bld(parms%aggr_prol,a,desc_a,&
-         & parms,ac,desc_ac,op_prol,op_restr,t_prol,info)
+            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,&
+subroutine amg_s_parmatch_smth_bld(dol1smoothing,ag,a,desc_a,ilaggr,nlaggr,&
-     & ac,desc_ac,op_prol,op_restr,t_prol,info)
+                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,&
+subroutine amg_s_parmatch_unsmth_bld(dol1smoothing,ag,a,desc_a,ilaggr,nlaggr,&
-     & ac,desc_ac,op_prol,op_restr,t_prol,info)
+                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_ptap_bld.f90

@@ -216,6 +216,7 @@ subroutine amg_s_ptap_bld(a_csr,desc_a,nlaggr,parms,ac,&
 
   call ac_csr%set_nrows(desc_ac%get_local_rows())
   call ac_csr%set_ncols(desc_ac%get_local_cols())
+  call ac_csr%clean_zeros(info)
   call ac%mv_from(ac_csr)
   call ac%set_asb()
 
@@ -420,6 +421,7 @@ subroutine amg_s_ls_ptap_bld(a_csr,desc_a,nlaggr,parms,ac,&
 
   call ac_csr%set_nrows(desc_ac%get_local_rows())
   call ac_csr%set_ncols(desc_ac%get_local_cols())
+  call ac_csr%clean_zeros(info)
   call ac%mv_from(ac_csr)
   call ac%set_asb()
 
@@ -629,6 +631,7 @@ subroutine amg_ls_ptap_bld(a_csr,desc_a,nlaggr,parms,ac,&
 
   call ac_csr%set_nrows(desc_ac%get_local_rows())
   call ac_csr%set_ncols(desc_ac%get_local_cols())
+  call ac_csr%clean_zeros(info)
   call ac%mv_from(ac_csr)
   call ac%set_asb()
 

amgprec/impl/aggregator/amg_s_rap.f90

@@ -142,6 +142,7 @@ subroutine amg_s_rap(a_csr,desc_a,nlaggr,parms,ac,&
 
   call ac_csr%set_nrows(desc_ac%get_local_rows())
   call ac_csr%set_ncols(desc_ac%get_local_cols())
+  call ac_csr%clean_zeros(info)
   call ac%mv_from(ac_csr)
   call ac%set_asb()
 

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,&
+subroutine amg_saggrmat_minnrg_bld(dol1smoothing,a,desc_a,ilaggr,nlaggr,&
-     & ac,desc_ac,op_prol,op_restr,t_prol,info)
+                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(dol1smoothing,a,desc_a,ilaggr,nlaggr,&
-subroutine amg_saggrmat_nosmth_bld(a,desc_a,ilaggr,nlaggr,parms,&
+                parms,ac,desc_ac,op_prol,op_restr,t_prol,info)
-     & 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,16 +104,18 @@
 !    info       -  integer, output.
 !                  Error code.
 !
-subroutine amg_saggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,&
+subroutine amg_saggrmat_smth_bld(dol1smoothing,a,desc_a,ilaggr,nlaggr,&
-     & ac,desc_ac,op_prol,op_restr,t_prol,info)
+                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
   use amg_s_base_aggregator_mod
+!  use, intrinsic :: ieee_arithmetic
 
   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 +136,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 +145,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 +178,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()
@@ -185,7 +193,7 @@ subroutine amg_saggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,&
 
   naggrm1 = sum(nlaggr(1:me))
   naggrp1 = sum(nlaggr(1:me+1))
-  filter_mat = (parms%aggr_filter == amg_filter_mat_)
+  filter_mat = (parms%aggr_filter == amg_filter_mat_).or.(parms%aggr_filter == amg_filter_prow_mat_)
 
   !
   ! naggr: number of local aggregates
@@ -200,6 +208,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,9 +256,15 @@ 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
+      else if (parms%aggr_filter == amg_filter_mat_) then
+        ! We perform filtering in the standard way assuming that A is an M-matrix
         acsrf%val(jd)=acsrf%val(jd)-tmp
+      else if (parms%aggr_filter == amg_filter_prow_mat_) then
+        ! We are probably doing l1-correction, hence we want to preserve the
+        ! row sum of the matrix: note the change in sign
+        acsrf%val(jd)=acsrf%val(jd)+tmp
       end if
     enddo
     !$OMP end parallel do 
@@ -240,7 +272,6 @@ subroutine amg_saggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,&
     call acsrf%clean_zeros(info)
   end if
 
-
   !$OMP parallel do private(i) schedule(static)
   do i=1,size(adiag)
     if (adiag(i) /= szero) then
@@ -252,14 +283,17 @@ 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 (parms%aggr_eig == amg_max_norm_) then 
+    if ( (parms%aggr_filter == amg_filter_prow_mat_).and.(do_l1correction) ) then
+      ! For l1-Jacobi this can be estimated with 1:
+      ! this makes sense only if we are preserving the row-sum!
+      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)))
       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 +356,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,&
+    call amg_zaggrmat_nosmth_bld(parms%aggr_prol,a,desc_a,ilaggr,&
-         & parms,ac,desc_ac,op_prol,op_restr,t_prol,info)
+            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, &
+    call amg_zaggrmat_smth_bld(parms%aggr_prol,a,desc_a,&
-         & parms,ac,desc_ac,op_prol,op_restr,t_prol,info)
+            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, &
+    call amg_zaggrmat_minnrg_bld(parms%aggr_prol,a,desc_a,ilaggr,&
-         & parms,ac,desc_ac,op_prol,op_restr,t_prol,info)
+            nlaggr,parms,ac,desc_ac,op_prol,op_restr,t_prol,info)
 
   case default
     info = psb_err_internal_error_

amgprec/impl/aggregator/amg_z_ptap_bld.f90

@@ -216,6 +216,7 @@ subroutine amg_z_ptap_bld(a_csr,desc_a,nlaggr,parms,ac,&
 
   call ac_csr%set_nrows(desc_ac%get_local_rows())
   call ac_csr%set_ncols(desc_ac%get_local_cols())
+  call ac_csr%clean_zeros(info)
   call ac%mv_from(ac_csr)
   call ac%set_asb()
 
@@ -420,6 +421,7 @@ subroutine amg_z_lz_ptap_bld(a_csr,desc_a,nlaggr,parms,ac,&
 
   call ac_csr%set_nrows(desc_ac%get_local_rows())
   call ac_csr%set_ncols(desc_ac%get_local_cols())
+  call ac_csr%clean_zeros(info)
   call ac%mv_from(ac_csr)
   call ac%set_asb()
 
@@ -629,6 +631,7 @@ subroutine amg_lz_ptap_bld(a_csr,desc_a,nlaggr,parms,ac,&
 
   call ac_csr%set_nrows(desc_ac%get_local_rows())
   call ac_csr%set_ncols(desc_ac%get_local_cols())
+  call ac_csr%clean_zeros(info)
   call ac%mv_from(ac_csr)
   call ac%set_asb()
 

amgprec/impl/aggregator/amg_z_rap.f90

@@ -142,6 +142,7 @@ subroutine amg_z_rap(a_csr,desc_a,nlaggr,parms,ac,&
 
   call ac_csr%set_nrows(desc_ac%get_local_rows())
   call ac_csr%set_ncols(desc_ac%get_local_cols())
+  call ac_csr%clean_zeros(info)
   call ac%mv_from(ac_csr)
   call ac%set_asb()
 

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,&
+subroutine amg_zaggrmat_minnrg_bld(dol1smoothing,a,desc_a,ilaggr,nlaggr,&
-     & ac,desc_ac,op_prol,op_restr,t_prol,info)
+                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(dol1smoothing,a,desc_a,ilaggr,nlaggr,&
-subroutine amg_zaggrmat_nosmth_bld(a,desc_a,ilaggr,nlaggr,parms,&
+                parms,ac,desc_ac,op_prol,op_restr,t_prol,info)
-     & 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,16 +104,18 @@
 !    info       -  integer, output.
 !                  Error code.
 !
-subroutine amg_zaggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,&
+subroutine amg_zaggrmat_smth_bld(dol1smoothing,a,desc_a,ilaggr,nlaggr,&
-     & ac,desc_ac,op_prol,op_restr,t_prol,info)
+                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
   use amg_z_base_aggregator_mod
+!  use, intrinsic :: ieee_arithmetic
 
   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 +136,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 +145,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 +178,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()
@@ -185,7 +193,7 @@ subroutine amg_zaggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,&
 
   naggrm1 = sum(nlaggr(1:me))
   naggrp1 = sum(nlaggr(1:me+1))
-  filter_mat = (parms%aggr_filter == amg_filter_mat_)
+  filter_mat = (parms%aggr_filter == amg_filter_mat_).or.(parms%aggr_filter == amg_filter_prow_mat_)
 
   !
   ! naggr: number of local aggregates
@@ -200,6 +208,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,9 +256,15 @@ 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
+      else if (parms%aggr_filter == amg_filter_mat_) then
+        ! We perform filtering in the standard way assuming that A is an M-matrix
         acsrf%val(jd)=acsrf%val(jd)-tmp
+      else if (parms%aggr_filter == amg_filter_prow_mat_) then
+        ! We are probably doing l1-correction, hence we want to preserve the
+        ! row sum of the matrix: note the change in sign
+        acsrf%val(jd)=acsrf%val(jd)+tmp
       end if
     enddo
     !$OMP end parallel do 
@@ -240,7 +272,6 @@ subroutine amg_zaggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,&
     call acsrf%clean_zeros(info)
   end if
 
-
   !$OMP parallel do private(i) schedule(static)
   do i=1,size(adiag)
     if (adiag(i) /= zzero) then
@@ -252,14 +283,17 @@ 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 (parms%aggr_eig == amg_max_norm_) then 
+    if ( (parms%aggr_filter == amg_filter_prow_mat_).and.(do_l1correction) ) then
+      ! For l1-Jacobi this can be estimated with 1:
+      ! this makes sense only if we are preserving the row-sum!
+      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)))
       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 +356,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
+      ! This is needed when the linmap object has been built
-      prec%precv(i)%linmap%p_desc_V => prec%precv(i)%base_desc
+      ! 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
+      ! This is needed when the linmap object has been built
-      prec%precv(i)%linmap%p_desc_V => prec%precv(i)%base_desc
+      ! 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
+      ! This is needed when the linmap object has been built
-      prec%precv(i)%linmap%p_desc_V => prec%precv(i)%base_desc
+      ! 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
+      ! This is needed when the linmap object has been built
-      prec%precv(i)%linmap%p_desc_V => prec%precv(i)%base_desc
+      ! 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_descr.f90

@@ -129,7 +129,7 @@ subroutine amg_c_base_onelev_descr(lv,il,nl,ilmin,info,iout,verbosity,prefix)
                & '                          avg:', &
                & lv%linmap%nagavg
         end if
-        write(iout_,'(a,1xa,1x,f14.2)') trim(prefix_),&
+        write(iout_,'(a,1x,a,1x,f14.2)') trim(prefix_),&
              & '          Aggregation   ratio: ', &
              &  lv%szratio
       end if

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_descr.f90

@@ -129,7 +129,7 @@ subroutine amg_d_base_onelev_descr(lv,il,nl,ilmin,info,iout,verbosity,prefix)
                & '                          avg:', &
                & lv%linmap%nagavg
         end if
-        write(iout_,'(a,1xa,1x,f14.2)') trim(prefix_),&
+        write(iout_,'(a,1x,a,1x,f14.2)') trim(prefix_),&
              & '          Aggregation   ratio: ', &
              &  lv%szratio
       end if

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_descr.f90

@@ -129,7 +129,7 @@ subroutine amg_s_base_onelev_descr(lv,il,nl,ilmin,info,iout,verbosity,prefix)
                & '                          avg:', &
                & lv%linmap%nagavg
         end if
-        write(iout_,'(a,1xa,1x,f14.2)') trim(prefix_),&
+        write(iout_,'(a,1x,a,1x,f14.2)') trim(prefix_),&
              & '          Aggregation   ratio: ', &
              &  lv%szratio
       end if

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_descr.f90

@@ -129,7 +129,7 @@ subroutine amg_z_base_onelev_descr(lv,il,nl,ilmin,info,iout,verbosity,prefix)
                & '                          avg:', &
                & lv%linmap%nagavg
         end if
-        write(iout_,'(a,1xa,1x,f14.2)') trim(prefix_),&
+        write(iout_,'(a,1x,a,1x,f14.2)') trim(prefix_),&
              & '          Aggregation   ratio: ', &
              &  lv%szratio
       end if

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_c_jac_smoother_apply_vect.f90

@@ -40,7 +40,7 @@ subroutine amg_c_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
 
   use psb_base_mod
   use amg_c_diag_solver
-  use psb_base_krylov_conv_mod, only : log_conv
+  use psb_base_linsolve_conv_mod, only : log_conv
   use amg_c_jac_smoother, amg_protect_name => amg_c_jac_smoother_apply_vect
   implicit none
   type(psb_desc_type), intent(in)                 :: desc_data

amgprec/impl/smoother/amg_c_jac_smoother_clone_settings.f90

@@ -41,7 +41,7 @@ subroutine amg_c_jac_smoother_clone_settings(sm,smout,info)
   use amg_c_jac_smoother, amg_protect_name =>  amg_c_jac_smoother_clone_settings
   Implicit None
   ! Arguments
-  class(amg_c_jac_smoother_type), intent(inout)              :: sm
+  class(amg_c_jac_smoother_type), intent(inout)  :: sm
   class(amg_c_base_smoother_type), intent(inout) :: smout
   integer(psb_ipk_), intent(out)                 :: info
   integer(psb_ipk_)  :: err_act

amgprec/impl/smoother/amg_d_jac_smoother_apply_vect.f90

@@ -40,7 +40,7 @@ subroutine amg_d_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
 
   use psb_base_mod
   use amg_d_diag_solver
-  use psb_base_krylov_conv_mod, only : log_conv
+  use psb_base_linsolve_conv_mod, only : log_conv
   use amg_d_jac_smoother, amg_protect_name => amg_d_jac_smoother_apply_vect
   implicit none
   type(psb_desc_type), intent(in)                 :: desc_data

amgprec/impl/smoother/amg_d_jac_smoother_clone_settings.f90

@@ -41,7 +41,7 @@ subroutine amg_d_jac_smoother_clone_settings(sm,smout,info)
   use amg_d_jac_smoother, amg_protect_name =>  amg_d_jac_smoother_clone_settings
   Implicit None
   ! Arguments
-  class(amg_d_jac_smoother_type), intent(inout)              :: sm
+  class(amg_d_jac_smoother_type), intent(inout)  :: sm
   class(amg_d_base_smoother_type), intent(inout) :: smout
   integer(psb_ipk_), intent(out)                 :: info
   integer(psb_ipk_)  :: err_act

amgprec/impl/smoother/amg_d_poly_smoother_apply_vect.f90

@@ -40,7 +40,7 @@ subroutine amg_d_poly_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
 
   use psb_base_mod
   use amg_d_diag_solver
-  use psb_base_krylov_conv_mod, only : log_conv
+  use psb_base_linsolve_conv_mod, only : log_conv
   use amg_d_poly_smoother, amg_protect_name => amg_d_poly_smoother_apply_vect
   implicit none
   type(psb_desc_type), intent(in)                 :: desc_data
@@ -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_clone_settings.f90

@@ -41,7 +41,7 @@ subroutine amg_d_poly_smoother_clone_settings(sm,smout,info)
   use amg_d_poly_smoother, amg_protect_name =>  amg_d_poly_smoother_clone_settings
   Implicit None
   ! Arguments
-  class(amg_d_poly_smoother_type), intent(inout)              :: sm
+  class(amg_d_poly_smoother_type), intent(inout) :: sm
   class(amg_d_base_smoother_type), intent(inout) :: smout
   integer(psb_ipk_), intent(out)                 :: info
   integer(psb_ipk_)  :: err_act

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
+      write(0,*) 'Invalid choice for POLY_VARIANT, defaulting to amg_cheb_4_',val
-      sm%variant = amg_poly_lottes_
+      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_)
+  case(amg_cheb_4_)
-    write(iout_,*) trim(prefix_), '        variant:        ','POLY_LOTTES'
+    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_)
+  case(amg_cheb_4_opt_)
-    write(iout_,*) trim(prefix_), '        variant:        ','POLY_LOTTES_BETA'
+    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_)
+  case(amg_cheb_1_opt_)
-    write(iout_,*) trim(prefix_), '        variant:        ','POLY_NEW'
+    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_jac_smoother_apply_vect.f90

@@ -40,7 +40,7 @@ subroutine amg_s_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
 
   use psb_base_mod
   use amg_s_diag_solver
-  use psb_base_krylov_conv_mod, only : log_conv
+  use psb_base_linsolve_conv_mod, only : log_conv
   use amg_s_jac_smoother, amg_protect_name => amg_s_jac_smoother_apply_vect
   implicit none
   type(psb_desc_type), intent(in)                 :: desc_data