Fix CUDA defines for compilation

Update copyright
Merge branch 'development' of github.com:sfilippone/amg4psblas into development
2026-10-07 07:04:59 +00:00 · 2025-02-21 18:07:58 +01:00 · 2025-02-21 16:41:39 +01:00 · 2025-02-21 16:36:10 +01:00 · 2025-02-21 16:35:45 +01:00 · 2025-02-21 16:33:43 +01:00
377 changed files with 20264 additions and 11068 deletions
@@ -20,3 +20,5 @@ autom4te.cache
 # the executable from tests
 runs

+# Documentation temporary files
+docs/src/userguide.pdf
@@ -1,165 +0,0 @@
-Changelog. A lot less detailed than usual, at least for past
-	history.
-2022/05/20: Restart ChangeLog. Updated to new name AMG4PSBLAS, now using PSB3.8	
-2018/10/28: Fix interface to MUMPS and configry machinery. Require PSB 3.6.
-2018/10/10: ICTXT argument in prec%init().	
-2018/07/30: Fixes for Intel compilers. BootCMatch interface in examples. 	
-2018/05/14: Interface for extension of aggregation methods. 	
-2018/02/28: New cartesian distribution for sample programs.	
-2017/12/15: New WRK component of preconditioner levels, preallocation.	
-2017/10/25: New example input file formats. Added sample matrices.	
-2017/10/02: New CBIND in PSBLAS 3.5.0	
-2017/07/30: Refactored examples. Change default thresholds. 
-2017/05/31: New internal description of ML. 
-2017/05/16: Improve build process. 
-2017/04/20: Force %set interface. Update docs.	
-2017/04/03: Remove obsolete stuff. 	
-2017/03/17: Fixed level%cnv;  add coarse _solver tracker.	
-2017/02/18: Take out clean_zeros; changed NOFILTER; defined FBGS; take out
-	n_prec_levs.
-2017/02/12: Updated mat_dist usage, dubious SP shell for UMFPACK, fixes
-	for RPM packaging.      
-2017/02/02: Fix superlu configury	
-2016/11/12: Fix hierarchy/smoothers build to handle 1 level.	
-2016/10/03: Merged changes to hierearchy building.	
-2016/08/20: Reimplemented decoupled aggregation	
-2016/07/20: Refactored application of multilevel. Defined V,W and
-	K-cycles. 	
-2016/05/18: Reworked internals of PRECSET. Defined  Forward-Backward
-	Gauss-Seidel solver. Now available separate PRE and POST smoother
-	objects. 
-2016/03/30: MUMPS interface.	
-2016/02/28: Hybrid Gauss-Seidel method. 	
-2016/02/03: unify integer argument checks.	
-2015/12/15: defaults single vs. double precision. Use clean_zeros.	
-2015/12/08: new matdist interface	
-2015/10/17: configry fixes	
-2015/10/13: Fixes for SLUDIST versions 3 and 4
-2015/05/03: New heap interface	
-2015/04/21: INTENT fixes
-	
-2014/12/21: New error handling
-	
-2014/10/27: Added versioncheck to configure.
-	
-2014/03/31: New get_diag. 
-	 
-2013/11/07: Merged changes from experimental branch. Fix INCDIR in
-	    makefiles. 
-       
-2013/07/15: Fixes for UMFPACK 5.4, SuperLU 4.3, SuperLU_Dist 3.3
-	
-2013/04/05: CLONE method.	
-	
-2013/03/08: Reworked SET routines. 	
-	
-2012/12/10: Enable long_integers.
-	
-2012/12/05: Split smoother/solver objects.	
-	
-2012/04/30: New scheme to find dynamically the number of level based on
-	    the size of the coarse matrix
-	
-2012/01/10: Done split interface/implementation, plus subdir restructure.
-	
-2011/12/13: Start split interface/implementation to improve build time. 	
-	
-2011/11/25: Now works with _vect methods from PSBLAS. 	
-	
-2011/10/24: New test generation methods.
-	
-2011/06/15: Dump prolongator/restrictor
-	
-2011/04/14: Added MOLD argument(s) to precbld.       
-	
-2011/03/30: Fixed: descriptive methods, example programs.
-	
-2011/03/08: Re-factored modules for ILU methods.	
-	
-2011/03/04: Make X intent(inout) in APPLY to allow for preconditioners
-            using SPMM. 	
-	
-2011/03/02: New set methods.
-	
-2011/01/07: Fixed UMF interfacing for Z data.
-
-2011/01/04: Added UMF inteface for D data.
-	
-2011/01/02: Fix usage of DESC_DATA. Switched all names to F90 ending. 
-
-2010/12/16: Fix usage of replicated space descriptor.
-
-2010/11/16: Fix Jacobi smoother in case of empty off-diagonal. 	
-
-2010/11/04: Defined and tested single real and complex.
-
-2010/11/02: Aligned usage of sparse data type with psblas3.
-
-2009/12/22: Aligned constants with mld2p4 v1.2
-
-2009/12/11: First working version of double multilevel.
-
-2009/12/05: Inttroduction of Smoother/Solver object hierarchy.
-
-2009/09/23: Initial F2003 version. 	
-      
-	
-2009/01/28: Changed names from XbaseprcY to XbaseprecY.
-2009/01/27: Changed names from mld_transfer to mld_move_alloc.
-	
-2009/01/13: Repackaged the one-level preconditioners. Reorganized the
-            build routines, taking out mlprec_bld, and switching the
-	    number of levels when needed. 
-2008/10/27: Changed the definition of prec_type: repackaged with a
-	    onelev-prec-type, containing a baseprec and maps between
-	    index spaces. No performance impact; no changes to
-	    user-level interfaces.
-	
-2008/09/18: Changed mld_sizeof to integer(8); updated samples.	
-	
-2008/08/26: Fixed matrix generation in sample programs. 
-	
-2008/07/25: missing implicit none in mld_prec_type.	
-	
-2008/07/23: added HTML documentation  
-	
-2008/06/13: Fixed aggregation for replicated index spaces.
-
-2008/06/02: Threshold into decoupled aggregation algorithm.
-	
-2008/05/27: Single precision version.
-
-2008/03/09: Introduced configure script. 
-	
-2008/02/08: Merged changes from intermesh branch: we now have an
-	    inter_desc_type object. Cleaned up data allocation and
-	    variable initialization in multilevel prec application.
-
-2008/01/10: Merged various fixes for: prologues, unused variables,
-	    interface details.	
-2007/12/21: Merge version with prologues and internal docs.
-	
-2007/11/15: Created pargen example.	
-
-2007/11/14: Fix INTENT(IN) on X vector in preconditioner routines.
-
-2007/10/19: Merged in ILU(P,T). To be tested extensively.
-
-2007/10/17: Merged ILU(K) into trunk.
-
-2007/10/16: Fixed ILU(K), it now performs satisfactorily. Also updated
-	    ILU(0) to be more legible. 
-	
-2007/10/11: First working version of ILU(K). Still slow, there should
-	    be room for improvement. 
-	
-2007/10/09: Added benchmark code.
-
-2007/10/09: Added MILU_N_. Beware: values for UMF_ etc. have been
-	    shifted. 
-
-2007/10/02: To do: decide whether to name MLD_KRYLOV_MOD or
-	    PSB_KRYLOV_MOD. 
-
-2007/10/01: Start of this changelog. MLD2P4 now has a different
-	    structure, to enable a build not embedded in PSBLAS.
@@ -1,10 +1,10 @@
  
   
-                             AMG4PSBLAS version 1.1
+                             AMG4PSBLAS version 1.2
    Algebraic Multigrid Package
-               based on PSBLAS (Parallel Sparse BLAS version 3.8)
+               based on PSBLAS (Parallel Sparse BLAS version 3.9)
    
-    (C) Copyright 2022
+    (C) Copyright 2025
  
        Salvatore Filippone  
        Pasqua D'Ambra   
@@ -75,7 +75,7 @@ CDEFINES=$(AMGCDEFINES)
 AMGFDEFINES=@AMGFDEFINES@ $(PSBFDEFINES)
 FDEFINES=$(AMGFDEFINES)

-CXXDEFINES=@AMGCXXDEFINES@
+CXXDEFINES=@AMGCXXDEFINES@ $(PSBCXXDEFINES)

@COMPILERULES@

@@ -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


@@ -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

@@ -44,8 +44,10 @@ cleanlib:
 	(cd include; /bin/rm -f *.a *$(.mod) *$(.fh))
 	(cd modules; /bin/rm -f *.a *$(.mod) *$(.fh))

-veryclean: cleanlib
-	(cd amgprec && $(MAKE) veryclean)
+distclean: clean samplesclean
+	/bin/rm -fr Make.inc
+
+samplesclean: clean
 	(cd samples/simple/fileread && $(MAKE) clean)
 	(cd samples/simple/pdegen && $(MAKE) clean)
 	(cd samples/advanced/fileread && $(MAKE) clean)
@@ -54,5 +56,6 @@ veryclean: cleanlib
 check: all
 	make check -C samples/advanced/pdegen

-clean:
+clean:  cleanlib
 	(cd amgprec && $(MAKE) clean)
+	(cd cbind && $(MAKE) clean)
@@ -1,54 +1,76 @@
-                         AMG4PSBLAS  
- Algebraic Multigrid Package  based on PSBLAS (Parallel Sparse BLAS version 3.8)
-    
-Salvatore Filippone    (University of Rome Tor Vergata and IAC-CNR)
-Pasqua D'Ambra         (IAC-CNR, Naples, IT)
-Fabio Durastante       (IAC-CNR, Naples, IT)
+# AMG4PSBLAS  v1.2
+Algebraic Multigrid Package  based on [PSBLAS](https://github.com/sfilippone/psblas3) (Parallel Sparse BLAS version 3.9)

---------------------------------------------------------------------
+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)
-preconditioners for  the iterative solution of large and sparse linear systems.
+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.

-It is an evolution of MLD2P4 (see LICENSE.MLD2P4), but it has been
-thoroughly reworked, and it is sufficiently different to warrant a new
-project name.
+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.

+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,
-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.
+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.

+## 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.
-2. Tweak Make.inc if you are not satisfied.
-3. make; 
+   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.
+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;
-  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.
+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. See also the samples/cuda folder.

+### 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
@@ -1,5 +1,19 @@
 WHAT'S NEW

+AMG4PSBLAS
+Version 1.2
+   1. New polynomial smoothers.
+   2. Introduced L1-variants
+   3. Reorganization of sample programs.
+   
+Version 1.1
+   1. Reworked approximate inverse solvers.
+   
+
+Version 1.0
+   1. Transitioned from MLD2P4
+
+MLD2P4
 Version 2.1
   1. The multigrid preconditioner now include fully general V- and
      W-cycles. We also support K-cycles, both for symmetric and 
@@ -81,9 +81,9 @@ module amg_base_prec_type
  !
  ! Version numbers
  !
-  character(len=*), parameter   :: amg_version_string_ = "1.1.0"
+  character(len=*), parameter   :: amg_version_string_ = "1.2.0"
  integer(psb_ipk_), parameter  :: amg_version_major_  = 1
-  integer(psb_ipk_), parameter  :: amg_version_minor_  = 1
+  integer(psb_ipk_), parameter  :: amg_version_minor_  = 2
  integer(psb_ipk_), parameter  :: amg_patchlevel_     = 0

  type amg_ml_parms
@@ -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_smooth_prol_ = 1
-  integer(psb_ipk_), parameter :: amg_min_energy_  = 2
+  integer(psb_ipk_), parameter :: amg_no_smooth_      = 0
+  integer(psb_ipk_), parameter :: amg_smooth_prol_    = 1
+  integer(psb_ipk_), parameter :: amg_l1_smooth_prol_ = 2
+  integer(psb_ipk_), parameter :: amg_min_energy_     = 3
  ! Disabling min_energy  for the time being.
-  integer(psb_ipk_), parameter :: amg_max_aggr_prol_=amg_smooth_prol_
+  integer(psb_ipk_), parameter :: amg_max_aggr_prol_= amg_l1_smooth_prol_
  !
  ! Legal values for entry: amg_aggr_filter_
  !
-  integer(psb_ipk_), parameter :: amg_no_filter_mat_  = 0
-  integer(psb_ipk_), parameter :: amg_filter_mat_     = 1
-  integer(psb_ipk_), parameter :: amg_max_filter_mat_ = amg_filter_mat_
+  integer(psb_ipk_), parameter :: amg_no_filter_mat_   = 0
+  integer(psb_ipk_), parameter :: amg_filter_mat_      = 1
+  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,9 +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_poly_lottes_beta_ = 1
-  integer(psb_ipk_), parameter :: amg_poly_new_         = 2
+  integer(psb_ipk_), parameter :: amg_cheb_4_     = 0
+  integer(psb_ipk_), parameter :: amg_cheb_4_opt_ = 1
+  integer(psb_ipk_), parameter :: amg_cheb_1_opt_ = 2
+  integer(psb_ipk_), parameter :: amg_poly_dbg_   = 8
  
  integer(psb_ipk_), parameter :: amg_poly_rho_est_power_ = 0

@@ -375,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      ',&
-       &           'min energy    ','bizr. smoothed'/)
+       &  aggr_prols(0:4)=(/'unsmoothed    ','smoothed      ',&
+       &  'l1-smoothed   ','min energy    ','bizr. smoothed'/)
  character(len=15), parameter, private :: &
-       &  aggr_filters(0:1)=(/'no filtering  ','filtering     '/)
+       &  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 :: &
@@ -547,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')
@@ -569,12 +575,14 @@ contains
      val = amg_as_
    case('POLY')
      val = amg_poly_
-    case('POLY_LOTTES')
-      val = amg_poly_lottes_
-    case('POLY_LOTTES_BETA')
-      val = amg_poly_lottes_beta_
-    case('POLY_NEW')
-      val = amg_poly_new_
+    case('CHEB_4')
+      val = amg_cheb_4_
+    case('CHEB_4_OPT')
+      val = amg_cheb_4_opt_
+    case('CHEB_1_OPT')
+      val = amg_cheb_1_opt_
+    case('POLY_DBG')
+      val = amg_poly_dbg_
    case('POLY_RHO_EST_POWER')
      val =  amg_poly_rho_est_power_
    case('A_NORMI')
@@ -585,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')
@@ -62,9 +62,6 @@ module amg_c_ainv_solver
    procedure, pass(sv) :: cseti   => amg_c_ainv_solver_cseti
    procedure, pass(sv) :: csetc   => amg_c_ainv_solver_csetc
    procedure, pass(sv) :: csetr   => amg_c_ainv_solver_csetr
-!!$    procedure, pass(sv) :: seti    => amg_c_ainv_solver_seti
-!!$    procedure, pass(sv) :: setc    => amg_c_ainv_solver_setc
-!!$    procedure, pass(sv) :: setr    => amg_c_ainv_solver_setr
    procedure, pass(sv) :: descr   => amg_c_ainv_solver_descr
    procedure, pass(sv) :: default => c_ainv_solver_default
    procedure, nopass   :: stringval  => c_ainv_stringval
@@ -91,9 +88,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_base_solver_type), allocatable, intent(inout) :: svout
-      integer(psb_ipk_), intent(out)                            :: info
+      class(amg_c_ainv_solver_type), intent(inout) :: sv
+      class(amg_c_base_solver_type), intent(inout) :: svout
+      integer(psb_ipk_), intent(out)                 :: info
    end subroutine amg_c_ainv_solver_clone_settings
  end interface

@@ -106,7 +103,7 @@ module amg_c_ainv_solver
      Implicit None

      ! Arguments
-      type(psb_cspmat_type), intent(inout), target        :: a
+      type(psb_cspmat_type), intent(in), target           :: a
      Type(psb_desc_type), Intent(inout)                  :: desc_a
      class(amg_c_ainv_solver_type), intent(inout)        :: sv
      integer(psb_ipk_), intent(out)                      :: info
@@ -230,7 +230,7 @@ module amg_c_as_smoother
           & psb_desc_type, psb_c_base_sparse_mat, psb_ipk_,&
           & psb_i_base_vect_type
      implicit none 
-      type(psb_cspmat_type), intent(inout), target     :: a
+      type(psb_cspmat_type), intent(in), target        :: a
      Type(psb_desc_type), Intent(inout)                 :: desc_a 
      class(amg_c_as_smoother_type), intent(inout)       :: sm
      integer(psb_ipk_), intent(out)                     :: info
@@ -237,7 +237,7 @@ module amg_c_base_smoother_mod
           & psb_c_vect_type, psb_c_base_vect_type, psb_spk_, &
           & amg_c_base_smoother_type, psb_ipk_, psb_i_base_vect_type
      ! Arguments
-      type(psb_cspmat_type), intent(inout), target  :: a
+      type(psb_cspmat_type), intent(in), target     :: a
      Type(psb_desc_type), Intent(inout)              :: desc_a 
      class(amg_c_base_smoother_type), intent(inout) :: sm 
      integer(psb_ipk_), intent(out)                  :: info
@@ -170,7 +170,7 @@ module amg_c_base_solver_mod
      Implicit None
      
      ! Arguments
-      type(psb_cspmat_type), intent(inout), target          :: a
+      type(psb_cspmat_type), intent(in), target             :: a
      Type(psb_desc_type), Intent(inout)                    :: desc_a 
      class(amg_c_base_solver_type), intent(inout)          :: sv
      integer(psb_ipk_), intent(out)                        :: info
@@ -119,7 +119,7 @@ module amg_c_diag_solver
      import :: psb_desc_type, psb_cspmat_type,  psb_c_base_sparse_mat, &
           & psb_c_vect_type, psb_c_base_vect_type, psb_spk_, &
           & amg_c_diag_solver_type, psb_ipk_, psb_i_base_vect_type      
-      type(psb_cspmat_type), intent(inout), target        :: a
+      type(psb_cspmat_type), intent(in), target           :: a
      Type(psb_desc_type), Intent(inout)                    :: desc_a 
      class(amg_c_diag_solver_type), intent(inout)        :: sv
      integer(psb_ipk_), intent(out)                        :: info
@@ -331,7 +331,7 @@ module amg_c_l1_diag_solver
      import :: psb_desc_type, psb_cspmat_type,  psb_c_base_sparse_mat, &
           & psb_c_vect_type, psb_c_base_vect_type, psb_spk_, &
           & amg_c_l1_diag_solver_type, psb_ipk_, psb_i_base_vect_type      
-      type(psb_cspmat_type), intent(inout), target        :: a
+      type(psb_cspmat_type), intent(in), target           :: a
      Type(psb_desc_type), Intent(inout)                    :: desc_a 
      class(amg_c_l1_diag_solver_type), intent(inout)        :: sv
      integer(psb_ipk_), intent(out)                        :: info
@@ -181,7 +181,7 @@ module amg_c_gs_solver
           & psb_cspmat_type, psb_c_base_sparse_mat, psb_c_base_vect_type,&
           & psb_ipk_, psb_i_base_vect_type
      implicit none 
-      type(psb_cspmat_type), intent(inout), target        :: a
+      type(psb_cspmat_type), intent(in), target           :: a
      Type(psb_desc_type), Intent(inout)                  :: desc_a 
      class(amg_c_gs_solver_type), intent(inout)         :: sv
      integer(psb_ipk_), intent(out)                      :: info
@@ -195,7 +195,7 @@ module amg_c_gs_solver
           & psb_cspmat_type, psb_c_base_sparse_mat, psb_c_base_vect_type,&
           & psb_ipk_, psb_i_base_vect_type
      implicit none 
-      type(psb_cspmat_type), intent(inout), target        :: a
+      type(psb_cspmat_type), intent(in), target           :: a
      Type(psb_desc_type), Intent(inout)                  :: desc_a 
      class(amg_c_bwgs_solver_type), intent(inout)         :: sv
      integer(psb_ipk_), intent(out)                      :: info
@@ -1,125 +0,0 @@
-!  
-!   
-!                             AMG4PSBLAS version 1.0
-!    Algebraic Multigrid Package
-!               based on PSBLAS (Parallel Sparse BLAS version 3.7)
-!    
-!    (C) Copyright 2020 
-!  
-!        Salvatore Filippone  
-!        Pasqua D'Ambra   
-!        Fabio Durastante        
-!   
-!    Redistribution and use in source and binary forms, with or without
-!    modification, are permitted provided that the following conditions
-!    are met:
-!      1. Redistributions of source code must retain the above copyright
-!         notice, this list of conditions and the following disclaimer.
-!      2. Redistributions in binary form must reproduce the above copyright
-!         notice, this list of conditions, and the following disclaimer in the
-!         documentation and/or other materials provided with the distribution.
-!      3. The name of the AMG4PSBLAS group or the names of its contributors may
-!         not be used to endorse or promote products derived from this
-!         software without specific written permission.
-!   
-!    THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
-!    ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
-!    TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
-!    PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE AMG4PSBLAS GROUP OR ITS CONTRIBUTORS
-!    BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
-!    CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
-!    SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
-!    INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
-!    CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
-!    ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
-!    POSSIBILITY OF SUCH DAMAGE.
-!   
-!  
-!
-!
-!  The aggregator object hosts the aggregation method for building
-!  the multilevel hierarchy. This variant is based on the hybrid method
-!  presented in 
-!
-!    S. Gratton, P. Henon, P. Jiranek and X. Vasseur:
-!    Reducing complexity of algebraic multigrid by aggregation
-!    Numerical Lin. Algebra with Applications, 2016, 23:501-518
-!    
-module amg_c_hybrid_aggregator_mod
-
-  use amg_c_dec_aggregator_mod
-  !
-  !   sm           -  class(amg_T_base_smoother_type), allocatable
-  !                   The current level preconditioner (aka smoother).
-  !   parms        -  type(amg_RTml_parms)
-  !                   The parameters defining the multilevel strategy.
-  !   ac           -  The local part of the current-level matrix, built by
-  !                   coarsening the previous-level matrix.
-  !   desc_ac      -  type(psb_desc_type).
-  !                   The communication descriptor associated to the matrix
-  !                   stored in ac.
-  !   base_a       -  type(psb_Tspmat_type), pointer.
-  !                   Pointer (really a pointer!) to the local part of the current 
-  !                   matrix (so we have a unified treatment of residuals).
-  !                   We need this to avoid passing explicitly the current matrix
-  !                   to the routine which applies the preconditioner.
-  !   base_desc    -  type(psb_desc_type), pointer.
-  !                   Pointer to the communication descriptor associated to the
-  !                   matrix pointed by base_a.
-  !   map          -  Stores the maps (restriction and prolongation) between the
-  !                   vector spaces associated to the index spaces of the previous
-  !                   and current levels.
-  !
-  !   Methods:  
-  !     Most methods follow the encapsulation hierarchy: they take whatever action
-  !     is appropriate for the current object, then call the corresponding method for
-  !     the contained object.
-  !     As an example: the descr() method prints out a description of the
-  !     level. It starts by invoking the descr() method of the parms object,
-  !     then calls the descr() method of the smoother object. 
-  !
-  !    descr      -   Prints a description of the object.
-  !    default    -   Set default values
-  !    dump       -   Dump to file object contents
-  !    set        -   Sets various parameters; when a request is unknown
-  !                   it is passed to the smoother object for further processing.
-  !    check      -   Sanity checks.
-  !    sizeof     -   Total memory occupation in bytes
-  !    get_nzeros -   Number of nonzeros 
-  !
-  !
-  type, extends(amg_c_dec_aggregator_type) :: amg_c_hybrid_aggregator_type
-    
-  contains
-    procedure, pass(ag) :: bld_tprol    => amg_c_hybrid_aggregator_build_tprol
-    procedure, nopass   :: fmt          => amg_c_hybrid_aggregator_fmt
-  end type amg_c_hybrid_aggregator_type
-
-
-  interface
-    subroutine  amg_c_hybrid_aggregator_build_tprol(ag,parms,a,desc_a,ilaggr,nlaggr,op_prol,info)
-      import :: amg_c_hybrid_aggregator_type, psb_desc_type, psb_cspmat_type, psb_spk_,  &
-           & psb_ipk_, psb_long_int_k_, amg_sml_parms
-      implicit none
-      class(amg_c_hybrid_aggregator_type), target, intent(inout) :: ag
-      type(amg_sml_parms), intent(inout)  :: parms 
-      type(psb_cspmat_type), intent(in)   :: a
-      type(psb_desc_type), intent(in)     :: desc_a
-      integer(psb_ipk_), allocatable, intent(out) :: ilaggr(:),nlaggr(:)
-      type(psb_cspmat_type), intent(out)  :: op_prol
-      integer(psb_ipk_), intent(out)      :: info
-    end subroutine amg_c_hybrid_aggregator_build_tprol
-  end interface
-
-contains
-
-
-  function amg_c_hybrid_aggregator_fmt() result(val)
-    implicit none 
-    character(len=32)  :: val
-
-    val = "Hybrid Decoupled aggregation"
-  end function amg_c_hybrid_aggregator_fmt
-
-  
-end module amg_c_hybrid_aggregator_mod
@@ -123,7 +123,7 @@ contains
    Implicit None

    ! Arguments
-    type(psb_cspmat_type), intent(inout), target        :: a
+    type(psb_cspmat_type), intent(in), target           :: a
    Type(psb_desc_type), Intent(inout)                  :: desc_a 
    class(amg_c_id_solver_type), intent(inout)          :: sv
    integer(psb_ipk_), intent(out)                      :: info
@@ -144,7 +144,7 @@ module amg_c_ilu_solver
           & psb_cspmat_type, psb_c_base_sparse_mat, psb_c_base_vect_type,&
           & psb_ipk_, psb_i_base_vect_type
      implicit none 
-      type(psb_cspmat_type), intent(inout), target       :: a
+      type(psb_cspmat_type), intent(in), target           :: a
      Type(psb_desc_type), Intent(inout)                  :: desc_a 
      class(amg_c_ilu_solver_type), intent(inout)         :: sv
      integer(psb_ipk_), intent(out)                      :: info
@@ -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(:)
@@ -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_base_solver_type), allocatable, intent(inout) :: svout
-      integer(psb_ipk_), intent(out)                            :: info
+      class(amg_c_invk_solver_type), intent(inout) :: sv
+      class(amg_c_base_solver_type), intent(inout) :: svout
+      integer(psb_ipk_), intent(out)                 :: info
    end subroutine amg_c_invk_solver_clone_settings
  end interface

@@ -94,7 +94,7 @@ module amg_c_invk_solver
      Implicit None

      ! Arguments
-      type(psb_cspmat_type), intent(inout), target        :: a
+      type(psb_cspmat_type), intent(in), target           :: a
      Type(psb_desc_type), Intent(inout)                  :: desc_a
      class(amg_c_invk_solver_type), intent(inout)        :: sv
      integer(psb_ipk_), intent(out)                      :: info
@@ -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_base_solver_type), allocatable, intent(inout) :: svout
-      integer(psb_ipk_), intent(out)                            :: info
+      class(amg_c_invt_solver_type), intent(inout) :: sv
+      class(amg_c_base_solver_type), intent(inout) :: svout
+      integer(psb_ipk_), intent(out)                 :: info
    end subroutine amg_c_invt_solver_clone_settings
  end interface

@@ -94,7 +94,7 @@ module amg_c_invt_solver
      Implicit None

      ! Arguments
-      type(psb_cspmat_type), intent(inout), target        :: a
+      type(psb_cspmat_type), intent(in), target           :: a
      Type(psb_desc_type), Intent(inout)                  :: desc_a
      class(amg_c_invt_solver_type), intent(inout)        :: sv
      integer(psb_ipk_), intent(out)                      :: info
@@ -151,7 +151,7 @@ module amg_c_jac_smoother
      import :: psb_desc_type, amg_c_jac_smoother_type, psb_c_vect_type, psb_spk_, &
           & psb_cspmat_type, psb_c_base_sparse_mat, psb_c_base_vect_type,&
           & psb_ipk_, psb_i_base_vect_type
-      type(psb_cspmat_type), intent(inout), target        :: a
+      type(psb_cspmat_type), intent(in), target           :: a
      Type(psb_desc_type), Intent(inout)                  :: desc_a
      class(amg_c_jac_smoother_type), intent(inout)       :: sm
      integer(psb_ipk_), intent(out)                      :: info
@@ -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_base_smoother_type), allocatable, intent(inout) :: smout
+      class(amg_c_jac_smoother_type), intent(inout)  :: sm
+      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
@@ -274,7 +274,7 @@ module amg_c_jac_smoother
      import :: psb_desc_type, amg_c_l1_jac_smoother_type, psb_c_vect_type, &
           & psb_cspmat_type, psb_c_base_sparse_mat, psb_c_base_vect_type,&
           & psb_ipk_, psb_i_base_vect_type
-      type(psb_cspmat_type), intent(inout), target        :: a
+      type(psb_cspmat_type), intent(in), target           :: a
      Type(psb_desc_type), Intent(inout)                  :: desc_a
      class(amg_c_l1_jac_smoother_type), intent(inout)       :: sm
      integer(psb_ipk_), intent(out)                      :: info
@@ -143,7 +143,7 @@ module amg_c_jac_solver
           & psb_cspmat_type, psb_c_base_sparse_mat, psb_c_base_vect_type,&
           & psb_ipk_, psb_i_base_vect_type
      implicit none 
-      type(psb_cspmat_type), intent(inout), target        :: a
+      type(psb_cspmat_type), intent(in), target           :: a
      Type(psb_desc_type), Intent(inout)                  :: desc_a 
      class(amg_c_jac_solver_type), intent(inout)         :: sv
      integer(psb_ipk_), intent(out)                      :: info
@@ -160,7 +160,7 @@ module amg_c_jac_solver
           & psb_cspmat_type, psb_c_base_sparse_mat, psb_c_base_vect_type,&
           & psb_ipk_, psb_i_base_vect_type
      implicit none 
-      type(psb_cspmat_type), intent(inout), target        :: a
+      type(psb_cspmat_type), intent(in), target           :: a
      Type(psb_desc_type), Intent(inout)                  :: desc_a 
      class(amg_c_l1_jac_solver_type), intent(inout)         :: sv
      integer(psb_ipk_), intent(out)                      :: info
@@ -174,7 +174,7 @@ module amg_c_krm_solver
           & psb_cspmat_type, psb_c_base_sparse_mat, psb_c_base_vect_type,&
           & psb_ipk_, psb_i_base_vect_type
      implicit none
-      type(psb_cspmat_type), intent(inout), target        :: a
+      type(psb_cspmat_type), intent(in), target           :: a
      Type(psb_desc_type), Intent(inout)                  :: desc_a
      class(amg_c_krm_solver_type), intent(inout)         :: sv
      integer(psb_ipk_), intent(out)                      :: info
@@ -163,7 +163,7 @@ module amg_c_mumps_solver
      Implicit None

      ! Arguments
-      type(psb_cspmat_type), intent(inout), target        :: a
+      type(psb_cspmat_type), intent(in), target           :: a
      Type(psb_desc_type), Intent(inout)                  :: desc_a 
      class(amg_c_mumps_solver_type), intent(inout)       :: sv
      integer(psb_ipk_), intent(out)                      :: info
@@ -256,7 +256,7 @@ contains
    Implicit None

    ! Arguments
-    type(psb_cspmat_type), intent(inout), target        :: a
+    type(psb_cspmat_type), intent(in), target           :: a
    Type(psb_desc_type), Intent(inout)                  :: desc_a 
    class(amg_c_slu_solver_type), intent(inout)         :: sv
    integer, intent(out)                                :: info
@@ -62,9 +62,6 @@ module amg_d_ainv_solver
    procedure, pass(sv) :: cseti   => amg_d_ainv_solver_cseti
    procedure, pass(sv) :: csetc   => amg_d_ainv_solver_csetc
    procedure, pass(sv) :: csetr   => amg_d_ainv_solver_csetr
-!!$    procedure, pass(sv) :: seti    => amg_d_ainv_solver_seti
-!!$    procedure, pass(sv) :: setc    => amg_d_ainv_solver_setc
-!!$    procedure, pass(sv) :: setr    => amg_d_ainv_solver_setr
    procedure, pass(sv) :: descr   => amg_d_ainv_solver_descr
    procedure, pass(sv) :: default => d_ainv_solver_default
    procedure, nopass   :: stringval  => d_ainv_stringval
@@ -91,9 +88,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_base_solver_type), allocatable, intent(inout) :: svout
-      integer(psb_ipk_), intent(out)                            :: info
+      class(amg_d_ainv_solver_type), intent(inout) :: sv
+      class(amg_d_base_solver_type), intent(inout) :: svout
+      integer(psb_ipk_), intent(out)                 :: info
    end subroutine amg_d_ainv_solver_clone_settings
  end interface

@@ -106,7 +103,7 @@ module amg_d_ainv_solver
      Implicit None

      ! Arguments
-      type(psb_dspmat_type), intent(inout), target        :: a
+      type(psb_dspmat_type), intent(in), target           :: a
      Type(psb_desc_type), Intent(inout)                  :: desc_a
      class(amg_d_ainv_solver_type), intent(inout)        :: sv
      integer(psb_ipk_), intent(out)                      :: info
@@ -230,7 +230,7 @@ module amg_d_as_smoother
           & psb_desc_type, psb_d_base_sparse_mat, psb_ipk_,&
           & psb_i_base_vect_type
      implicit none 
-      type(psb_dspmat_type), intent(inout), target     :: a
+      type(psb_dspmat_type), intent(in), target        :: a
      Type(psb_desc_type), Intent(inout)                 :: desc_a 
      class(amg_d_as_smoother_type), intent(inout)       :: sm
      integer(psb_ipk_), intent(out)                     :: info
@@ -237,7 +237,7 @@ module amg_d_base_smoother_mod
           & psb_d_vect_type, psb_d_base_vect_type, psb_dpk_, &
           & amg_d_base_smoother_type, psb_ipk_, psb_i_base_vect_type
      ! Arguments
-      type(psb_dspmat_type), intent(inout), target  :: a
+      type(psb_dspmat_type), intent(in), target     :: a
      Type(psb_desc_type), Intent(inout)              :: desc_a 
      class(amg_d_base_smoother_type), intent(inout) :: sm 
      integer(psb_ipk_), intent(out)                  :: info
@@ -170,7 +170,7 @@ module amg_d_base_solver_mod
      Implicit None
      
      ! Arguments
-      type(psb_dspmat_type), intent(inout), target          :: a
+      type(psb_dspmat_type), intent(in), target             :: a
      Type(psb_desc_type), Intent(inout)                    :: desc_a 
      class(amg_d_base_solver_type), intent(inout)          :: sv
      integer(psb_ipk_), intent(out)                        :: info
@@ -119,7 +119,7 @@ module amg_d_diag_solver
      import :: psb_desc_type, psb_dspmat_type,  psb_d_base_sparse_mat, &
           & psb_d_vect_type, psb_d_base_vect_type, psb_dpk_, &
           & amg_d_diag_solver_type, psb_ipk_, psb_i_base_vect_type      
-      type(psb_dspmat_type), intent(inout), target        :: a
+      type(psb_dspmat_type), intent(in), target           :: a
      Type(psb_desc_type), Intent(inout)                    :: desc_a 
      class(amg_d_diag_solver_type), intent(inout)        :: sv
      integer(psb_ipk_), intent(out)                        :: info
@@ -331,7 +331,7 @@ module amg_d_l1_diag_solver
      import :: psb_desc_type, psb_dspmat_type,  psb_d_base_sparse_mat, &
           & psb_d_vect_type, psb_d_base_vect_type, psb_dpk_, &
           & amg_d_l1_diag_solver_type, psb_ipk_, psb_i_base_vect_type      
-      type(psb_dspmat_type), intent(inout), target        :: a
+      type(psb_dspmat_type), intent(in), target           :: a
      Type(psb_desc_type), Intent(inout)                    :: desc_a 
      class(amg_d_l1_diag_solver_type), intent(inout)        :: sv
      integer(psb_ipk_), intent(out)                        :: info
@@ -181,7 +181,7 @@ module amg_d_gs_solver
           & psb_dspmat_type, psb_d_base_sparse_mat, psb_d_base_vect_type,&
           & psb_ipk_, psb_i_base_vect_type
      implicit none 
-      type(psb_dspmat_type), intent(inout), target        :: a
+      type(psb_dspmat_type), intent(in), target           :: a
      Type(psb_desc_type), Intent(inout)                  :: desc_a 
      class(amg_d_gs_solver_type), intent(inout)         :: sv
      integer(psb_ipk_), intent(out)                      :: info
@@ -195,7 +195,7 @@ module amg_d_gs_solver
           & psb_dspmat_type, psb_d_base_sparse_mat, psb_d_base_vect_type,&
           & psb_ipk_, psb_i_base_vect_type
      implicit none 
-      type(psb_dspmat_type), intent(inout), target        :: a
+      type(psb_dspmat_type), intent(in), target           :: a
      Type(psb_desc_type), Intent(inout)                  :: desc_a 
      class(amg_d_bwgs_solver_type), intent(inout)         :: sv
      integer(psb_ipk_), intent(out)                      :: info
@@ -1,125 +0,0 @@
-!  
-!   
-!                             AMG4PSBLAS version 1.0
-!    Algebraic Multigrid Package
-!               based on PSBLAS (Parallel Sparse BLAS version 3.7)
-!    
-!    (C) Copyright 2020 
-!  
-!        Salvatore Filippone  
-!        Pasqua D'Ambra   
-!        Fabio Durastante        
-!   
-!    Redistribution and use in source and binary forms, with or without
-!    modification, are permitted provided that the following conditions
-!    are met:
-!      1. Redistributions of source code must retain the above copyright
-!         notice, this list of conditions and the following disclaimer.
-!      2. Redistributions in binary form must reproduce the above copyright
-!         notice, this list of conditions, and the following disclaimer in the
-!         documentation and/or other materials provided with the distribution.
-!      3. The name of the AMG4PSBLAS group or the names of its contributors may
-!         not be used to endorse or promote products derived from this
-!         software without specific written permission.
-!   
-!    THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
-!    ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
-!    TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
-!    PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE AMG4PSBLAS GROUP OR ITS CONTRIBUTORS
-!    BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
-!    CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
-!    SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
-!    INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
-!    CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
-!    ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
-!    POSSIBILITY OF SUCH DAMAGE.
-!   
-!  
-!
-!
-!  The aggregator object hosts the aggregation method for building
-!  the multilevel hierarchy. This variant is based on the hybrid method
-!  presented in 
-!
-!    S. Gratton, P. Henon, P. Jiranek and X. Vasseur:
-!    Reducing complexity of algebraic multigrid by aggregation
-!    Numerical Lin. Algebra with Applications, 2016, 23:501-518
-!    
-module amg_d_hybrid_aggregator_mod
-
-  use amg_d_dec_aggregator_mod
-  !
-  !   sm           -  class(amg_T_base_smoother_type), allocatable
-  !                   The current level preconditioner (aka smoother).
-  !   parms        -  type(amg_RTml_parms)
-  !                   The parameters defining the multilevel strategy.
-  !   ac           -  The local part of the current-level matrix, built by
-  !                   coarsening the previous-level matrix.
-  !   desc_ac      -  type(psb_desc_type).
-  !                   The communication descriptor associated to the matrix
-  !                   stored in ac.
-  !   base_a       -  type(psb_Tspmat_type), pointer.
-  !                   Pointer (really a pointer!) to the local part of the current 
-  !                   matrix (so we have a unified treatment of residuals).
-  !                   We need this to avoid passing explicitly the current matrix
-  !                   to the routine which applies the preconditioner.
-  !   base_desc    -  type(psb_desc_type), pointer.
-  !                   Pointer to the communication descriptor associated to the
-  !                   matrix pointed by base_a.
-  !   map          -  Stores the maps (restriction and prolongation) between the
-  !                   vector spaces associated to the index spaces of the previous
-  !                   and current levels.
-  !
-  !   Methods:  
-  !     Most methods follow the encapsulation hierarchy: they take whatever action
-  !     is appropriate for the current object, then call the corresponding method for
-  !     the contained object.
-  !     As an example: the descr() method prints out a description of the
-  !     level. It starts by invoking the descr() method of the parms object,
-  !     then calls the descr() method of the smoother object. 
-  !
-  !    descr      -   Prints a description of the object.
-  !    default    -   Set default values
-  !    dump       -   Dump to file object contents
-  !    set        -   Sets various parameters; when a request is unknown
-  !                   it is passed to the smoother object for further processing.
-  !    check      -   Sanity checks.
-  !    sizeof     -   Total memory occupation in bytes
-  !    get_nzeros -   Number of nonzeros 
-  !
-  !
-  type, extends(amg_d_dec_aggregator_type) :: amg_d_hybrid_aggregator_type
-    
-  contains
-    procedure, pass(ag) :: bld_tprol    => amg_d_hybrid_aggregator_build_tprol
-    procedure, nopass   :: fmt          => amg_d_hybrid_aggregator_fmt
-  end type amg_d_hybrid_aggregator_type
-
-
-  interface
-    subroutine  amg_d_hybrid_aggregator_build_tprol(ag,parms,a,desc_a,ilaggr,nlaggr,op_prol,info)
-      import :: amg_d_hybrid_aggregator_type, psb_desc_type, psb_dspmat_type, psb_dpk_,  &
-           & psb_ipk_, psb_long_int_k_, amg_dml_parms
-      implicit none
-      class(amg_d_hybrid_aggregator_type), target, intent(inout) :: ag
-      type(amg_dml_parms), intent(inout)  :: parms 
-      type(psb_dspmat_type), intent(in)   :: a
-      type(psb_desc_type), intent(in)     :: desc_a
-      integer(psb_ipk_), allocatable, intent(out) :: ilaggr(:),nlaggr(:)
-      type(psb_dspmat_type), intent(out)  :: op_prol
-      integer(psb_ipk_), intent(out)      :: info
-    end subroutine amg_d_hybrid_aggregator_build_tprol
-  end interface
-
-contains
-
-
-  function amg_d_hybrid_aggregator_fmt() result(val)
-    implicit none 
-    character(len=32)  :: val
-
-    val = "Hybrid Decoupled aggregation"
-  end function amg_d_hybrid_aggregator_fmt
-
-  
-end module amg_d_hybrid_aggregator_mod
@@ -123,7 +123,7 @@ contains
    Implicit None

    ! Arguments
-    type(psb_dspmat_type), intent(inout), target        :: a
+    type(psb_dspmat_type), intent(in), target           :: a
    Type(psb_desc_type), Intent(inout)                  :: desc_a 
    class(amg_d_id_solver_type), intent(inout)          :: sv
    integer(psb_ipk_), intent(out)                      :: info
@@ -144,7 +144,7 @@ module amg_d_ilu_solver
           & psb_dspmat_type, psb_d_base_sparse_mat, psb_d_base_vect_type,&
           & psb_ipk_, psb_i_base_vect_type
      implicit none 
-      type(psb_dspmat_type), intent(inout), target       :: a
+      type(psb_dspmat_type), intent(in), target           :: a
      Type(psb_desc_type), Intent(inout)                  :: desc_a 
      class(amg_d_ilu_solver_type), intent(inout)         :: sv
      integer(psb_ipk_), intent(out)                      :: info
@@ -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(:)
@@ -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_base_solver_type), allocatable, intent(inout) :: svout
-      integer(psb_ipk_), intent(out)                            :: info
+      class(amg_d_invk_solver_type), intent(inout) :: sv
+      class(amg_d_base_solver_type), intent(inout) :: svout
+      integer(psb_ipk_), intent(out)                 :: info
    end subroutine amg_d_invk_solver_clone_settings
  end interface

@@ -94,7 +94,7 @@ module amg_d_invk_solver
      Implicit None

      ! Arguments
-      type(psb_dspmat_type), intent(inout), target        :: a
+      type(psb_dspmat_type), intent(in), target           :: a
      Type(psb_desc_type), Intent(inout)                  :: desc_a
      class(amg_d_invk_solver_type), intent(inout)        :: sv
      integer(psb_ipk_), intent(out)                      :: info
@@ -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_base_solver_type), allocatable, intent(inout) :: svout
-      integer(psb_ipk_), intent(out)                            :: info
+      class(amg_d_invt_solver_type), intent(inout) :: sv
+      class(amg_d_base_solver_type), intent(inout) :: svout
+      integer(psb_ipk_), intent(out)                 :: info
    end subroutine amg_d_invt_solver_clone_settings
  end interface

@@ -94,7 +94,7 @@ module amg_d_invt_solver
      Implicit None

      ! Arguments
-      type(psb_dspmat_type), intent(inout), target        :: a
+      type(psb_dspmat_type), intent(in), target           :: a
      Type(psb_desc_type), Intent(inout)                  :: desc_a
      class(amg_d_invt_solver_type), intent(inout)        :: sv
      integer(psb_ipk_), intent(out)                      :: info
@@ -151,7 +151,7 @@ module amg_d_jac_smoother
      import :: psb_desc_type, amg_d_jac_smoother_type, psb_d_vect_type, psb_dpk_, &
           & psb_dspmat_type, psb_d_base_sparse_mat, psb_d_base_vect_type,&
           & psb_ipk_, psb_i_base_vect_type
-      type(psb_dspmat_type), intent(inout), target        :: a
+      type(psb_dspmat_type), intent(in), target           :: a
      Type(psb_desc_type), Intent(inout)                  :: desc_a
      class(amg_d_jac_smoother_type), intent(inout)       :: sm
      integer(psb_ipk_), intent(out)                      :: info
@@ -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_base_smoother_type), allocatable, intent(inout) :: smout
+      class(amg_d_jac_smoother_type), intent(inout)  :: sm
+      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
@@ -274,7 +274,7 @@ module amg_d_jac_smoother
      import :: psb_desc_type, amg_d_l1_jac_smoother_type, psb_d_vect_type, &
           & psb_dspmat_type, psb_d_base_sparse_mat, psb_d_base_vect_type,&
           & psb_ipk_, psb_i_base_vect_type
-      type(psb_dspmat_type), intent(inout), target        :: a
+      type(psb_dspmat_type), intent(in), target           :: a
      Type(psb_desc_type), Intent(inout)                  :: desc_a
      class(amg_d_l1_jac_smoother_type), intent(inout)       :: sm
      integer(psb_ipk_), intent(out)                      :: info
@@ -143,7 +143,7 @@ module amg_d_jac_solver
           & psb_dspmat_type, psb_d_base_sparse_mat, psb_d_base_vect_type,&
           & psb_ipk_, psb_i_base_vect_type
      implicit none 
-      type(psb_dspmat_type), intent(inout), target        :: a
+      type(psb_dspmat_type), intent(in), target           :: a
      Type(psb_desc_type), Intent(inout)                  :: desc_a 
      class(amg_d_jac_solver_type), intent(inout)         :: sv
      integer(psb_ipk_), intent(out)                      :: info
@@ -160,7 +160,7 @@ module amg_d_jac_solver
           & psb_dspmat_type, psb_d_base_sparse_mat, psb_d_base_vect_type,&
           & psb_ipk_, psb_i_base_vect_type
      implicit none 
-      type(psb_dspmat_type), intent(inout), target        :: a
+      type(psb_dspmat_type), intent(in), target           :: a
      Type(psb_desc_type), Intent(inout)                  :: desc_a 
      class(amg_d_l1_jac_solver_type), intent(inout)         :: sv
      integer(psb_ipk_), intent(out)                      :: info
@@ -174,7 +174,7 @@ module amg_d_krm_solver
           & psb_dspmat_type, psb_d_base_sparse_mat, psb_d_base_vect_type,&
           & psb_ipk_, psb_i_base_vect_type
      implicit none
-      type(psb_dspmat_type), intent(inout), target        :: a
+      type(psb_dspmat_type), intent(in), target           :: a
      Type(psb_desc_type), Intent(inout)                  :: desc_a
      class(amg_d_krm_solver_type), intent(inout)         :: sv
      integer(psb_ipk_), intent(out)                      :: info
@@ -163,7 +163,7 @@ module amg_d_mumps_solver
      Implicit None

      ! Arguments
-      type(psb_dspmat_type), intent(inout), target        :: a
+      type(psb_dspmat_type), intent(in), target           :: a
      Type(psb_desc_type), Intent(inout)                  :: desc_a 
      class(amg_d_mumps_solver_type), intent(inout)       :: sv
      integer(psb_ipk_), intent(out)                      :: info
@@ -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
@@ -140,7 +140,7 @@ module amg_d_poly_smoother
      import :: psb_desc_type, amg_d_poly_smoother_type, psb_d_vect_type, psb_dpk_, &
           & psb_dspmat_type, psb_d_base_sparse_mat, psb_d_base_vect_type,&
           & psb_ipk_, psb_i_base_vect_type
-      type(psb_dspmat_type), intent(inout), target        :: a
+      type(psb_dspmat_type), intent(in), target           :: a
      Type(psb_desc_type), Intent(inout)                  :: desc_a
      class(amg_d_poly_smoother_type), intent(inout)       :: sm
      integer(psb_ipk_), intent(out)                      :: info
@@ -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_base_smoother_type), allocatable, intent(inout) :: smout
+      class(amg_d_poly_smoother_type), intent(inout) :: sm
+      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
@@ -256,7 +256,7 @@ contains
    Implicit None

    ! Arguments
-    type(psb_dspmat_type), intent(inout), target        :: a
+    type(psb_dspmat_type), intent(in), target           :: a
    Type(psb_desc_type), Intent(inout)                  :: desc_a 
    class(amg_d_slu_solver_type), intent(inout)         :: sv
    integer, intent(out)                                :: info
@@ -259,7 +259,7 @@ contains
    Implicit None

    ! Arguments
-    type(psb_dspmat_type), intent(inout), target        :: a
+    type(psb_dspmat_type), intent(in), target           :: a
    Type(psb_desc_type), Intent(inout)                  :: desc_a 
    class(amg_d_sludist_solver_type), intent(inout)     :: sv
    integer, intent(out)                                :: info
@@ -259,7 +259,7 @@ contains
    Implicit None

    ! Arguments
-    type(psb_dspmat_type), intent(inout), target        :: a
+    type(psb_dspmat_type), intent(in), target           :: a
    Type(psb_desc_type), Intent(inout)                  :: desc_a 
    class(amg_d_umf_solver_type), intent(inout)         :: sv
    integer, intent(out)                                :: info
@@ -62,9 +62,6 @@ module amg_s_ainv_solver
    procedure, pass(sv) :: cseti   => amg_s_ainv_solver_cseti
    procedure, pass(sv) :: csetc   => amg_s_ainv_solver_csetc
    procedure, pass(sv) :: csetr   => amg_s_ainv_solver_csetr
-!!$    procedure, pass(sv) :: seti    => amg_s_ainv_solver_seti
-!!$    procedure, pass(sv) :: setc    => amg_s_ainv_solver_setc
-!!$    procedure, pass(sv) :: setr    => amg_s_ainv_solver_setr
    procedure, pass(sv) :: descr   => amg_s_ainv_solver_descr
    procedure, pass(sv) :: default => s_ainv_solver_default
    procedure, nopass   :: stringval  => s_ainv_stringval
@@ -91,9 +88,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_base_solver_type), allocatable, intent(inout) :: svout
-      integer(psb_ipk_), intent(out)                            :: info
+      class(amg_s_ainv_solver_type), intent(inout) :: sv
+      class(amg_s_base_solver_type), intent(inout) :: svout
+      integer(psb_ipk_), intent(out)                 :: info
    end subroutine amg_s_ainv_solver_clone_settings
  end interface

@@ -106,7 +103,7 @@ module amg_s_ainv_solver
      Implicit None

      ! Arguments
-      type(psb_sspmat_type), intent(inout), target        :: a
+      type(psb_sspmat_type), intent(in), target           :: a
      Type(psb_desc_type), Intent(inout)                  :: desc_a
      class(amg_s_ainv_solver_type), intent(inout)        :: sv
      integer(psb_ipk_), intent(out)                      :: info
@@ -230,7 +230,7 @@ module amg_s_as_smoother
           & psb_desc_type, psb_s_base_sparse_mat, psb_ipk_,&
           & psb_i_base_vect_type
      implicit none 
-      type(psb_sspmat_type), intent(inout), target     :: a
+      type(psb_sspmat_type), intent(in), target        :: a
      Type(psb_desc_type), Intent(inout)                 :: desc_a 
      class(amg_s_as_smoother_type), intent(inout)       :: sm
      integer(psb_ipk_), intent(out)                     :: info
@@ -237,7 +237,7 @@ module amg_s_base_smoother_mod
           & psb_s_vect_type, psb_s_base_vect_type, psb_spk_, &
           & amg_s_base_smoother_type, psb_ipk_, psb_i_base_vect_type
      ! Arguments
-      type(psb_sspmat_type), intent(inout), target  :: a
+      type(psb_sspmat_type), intent(in), target     :: a
      Type(psb_desc_type), Intent(inout)              :: desc_a 
      class(amg_s_base_smoother_type), intent(inout) :: sm 
      integer(psb_ipk_), intent(out)                  :: info
@@ -170,7 +170,7 @@ module amg_s_base_solver_mod
      Implicit None
      
      ! Arguments
-      type(psb_sspmat_type), intent(inout), target          :: a
+      type(psb_sspmat_type), intent(in), target             :: a
      Type(psb_desc_type), Intent(inout)                    :: desc_a 
      class(amg_s_base_solver_type), intent(inout)          :: sv
      integer(psb_ipk_), intent(out)                        :: info
@@ -119,7 +119,7 @@ module amg_s_diag_solver
      import :: psb_desc_type, psb_sspmat_type,  psb_s_base_sparse_mat, &
           & psb_s_vect_type, psb_s_base_vect_type, psb_spk_, &
           & amg_s_diag_solver_type, psb_ipk_, psb_i_base_vect_type      
-      type(psb_sspmat_type), intent(inout), target        :: a
+      type(psb_sspmat_type), intent(in), target           :: a
      Type(psb_desc_type), Intent(inout)                    :: desc_a 
      class(amg_s_diag_solver_type), intent(inout)        :: sv
      integer(psb_ipk_), intent(out)                        :: info
@@ -331,7 +331,7 @@ module amg_s_l1_diag_solver
      import :: psb_desc_type, psb_sspmat_type,  psb_s_base_sparse_mat, &
           & psb_s_vect_type, psb_s_base_vect_type, psb_spk_, &
           & amg_s_l1_diag_solver_type, psb_ipk_, psb_i_base_vect_type      
-      type(psb_sspmat_type), intent(inout), target        :: a
+      type(psb_sspmat_type), intent(in), target           :: a
      Type(psb_desc_type), Intent(inout)                    :: desc_a 
      class(amg_s_l1_diag_solver_type), intent(inout)        :: sv
      integer(psb_ipk_), intent(out)                        :: info
@@ -181,7 +181,7 @@ module amg_s_gs_solver
           & psb_sspmat_type, psb_s_base_sparse_mat, psb_s_base_vect_type,&
           & psb_ipk_, psb_i_base_vect_type
      implicit none 
-      type(psb_sspmat_type), intent(inout), target        :: a
+      type(psb_sspmat_type), intent(in), target           :: a
      Type(psb_desc_type), Intent(inout)                  :: desc_a 
      class(amg_s_gs_solver_type), intent(inout)         :: sv
      integer(psb_ipk_), intent(out)                      :: info
@@ -195,7 +195,7 @@ module amg_s_gs_solver
           & psb_sspmat_type, psb_s_base_sparse_mat, psb_s_base_vect_type,&
           & psb_ipk_, psb_i_base_vect_type
      implicit none 
-      type(psb_sspmat_type), intent(inout), target        :: a
+      type(psb_sspmat_type), intent(in), target           :: a
      Type(psb_desc_type), Intent(inout)                  :: desc_a 
      class(amg_s_bwgs_solver_type), intent(inout)         :: sv
      integer(psb_ipk_), intent(out)                      :: info
@@ -1,125 +0,0 @@
-!  
-!   
-!                             AMG4PSBLAS version 1.0
-!    Algebraic Multigrid Package
-!               based on PSBLAS (Parallel Sparse BLAS version 3.7)
-!    
-!    (C) Copyright 2020 
-!  
-!        Salvatore Filippone  
-!        Pasqua D'Ambra   
-!        Fabio Durastante        
-!   
-!    Redistribution and use in source and binary forms, with or without
-!    modification, are permitted provided that the following conditions
-!    are met:
-!      1. Redistributions of source code must retain the above copyright
-!         notice, this list of conditions and the following disclaimer.
-!      2. Redistributions in binary form must reproduce the above copyright
-!         notice, this list of conditions, and the following disclaimer in the
-!         documentation and/or other materials provided with the distribution.
-!      3. The name of the AMG4PSBLAS group or the names of its contributors may
-!         not be used to endorse or promote products derived from this
-!         software without specific written permission.
-!   
-!    THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
-!    ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
-!    TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
-!    PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE AMG4PSBLAS GROUP OR ITS CONTRIBUTORS
-!    BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
-!    CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
-!    SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
-!    INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
-!    CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
-!    ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
-!    POSSIBILITY OF SUCH DAMAGE.
-!   
-!  
-!
-!
-!  The aggregator object hosts the aggregation method for building
-!  the multilevel hierarchy. This variant is based on the hybrid method
-!  presented in 
-!
-!    S. Gratton, P. Henon, P. Jiranek and X. Vasseur:
-!    Reducing complexity of algebraic multigrid by aggregation
-!    Numerical Lin. Algebra with Applications, 2016, 23:501-518
-!    
-module amg_s_hybrid_aggregator_mod
-
-  use amg_s_dec_aggregator_mod
-  !
-  !   sm           -  class(amg_T_base_smoother_type), allocatable
-  !                   The current level preconditioner (aka smoother).
-  !   parms        -  type(amg_RTml_parms)
-  !                   The parameters defining the multilevel strategy.
-  !   ac           -  The local part of the current-level matrix, built by
-  !                   coarsening the previous-level matrix.
-  !   desc_ac      -  type(psb_desc_type).
-  !                   The communication descriptor associated to the matrix
-  !                   stored in ac.
-  !   base_a       -  type(psb_Tspmat_type), pointer.
-  !                   Pointer (really a pointer!) to the local part of the current 
-  !                   matrix (so we have a unified treatment of residuals).
-  !                   We need this to avoid passing explicitly the current matrix
-  !                   to the routine which applies the preconditioner.
-  !   base_desc    -  type(psb_desc_type), pointer.
-  !                   Pointer to the communication descriptor associated to the
-  !                   matrix pointed by base_a.
-  !   map          -  Stores the maps (restriction and prolongation) between the
-  !                   vector spaces associated to the index spaces of the previous
-  !                   and current levels.
-  !
-  !   Methods:  
-  !     Most methods follow the encapsulation hierarchy: they take whatever action
-  !     is appropriate for the current object, then call the corresponding method for
-  !     the contained object.
-  !     As an example: the descr() method prints out a description of the
-  !     level. It starts by invoking the descr() method of the parms object,
-  !     then calls the descr() method of the smoother object. 
-  !
-  !    descr      -   Prints a description of the object.
-  !    default    -   Set default values
-  !    dump       -   Dump to file object contents
-  !    set        -   Sets various parameters; when a request is unknown
-  !                   it is passed to the smoother object for further processing.
-  !    check      -   Sanity checks.
-  !    sizeof     -   Total memory occupation in bytes
-  !    get_nzeros -   Number of nonzeros 
-  !
-  !
-  type, extends(amg_s_dec_aggregator_type) :: amg_s_hybrid_aggregator_type
-    
-  contains
-    procedure, pass(ag) :: bld_tprol    => amg_s_hybrid_aggregator_build_tprol
-    procedure, nopass   :: fmt          => amg_s_hybrid_aggregator_fmt
-  end type amg_s_hybrid_aggregator_type
-
-
-  interface
-    subroutine  amg_s_hybrid_aggregator_build_tprol(ag,parms,a,desc_a,ilaggr,nlaggr,op_prol,info)
-      import :: amg_s_hybrid_aggregator_type, psb_desc_type, psb_sspmat_type, psb_spk_,  &
-           & psb_ipk_, psb_long_int_k_, amg_sml_parms
-      implicit none
-      class(amg_s_hybrid_aggregator_type), target, intent(inout) :: ag
-      type(amg_sml_parms), intent(inout)  :: parms 
-      type(psb_sspmat_type), intent(in)   :: a
-      type(psb_desc_type), intent(in)     :: desc_a
-      integer(psb_ipk_), allocatable, intent(out) :: ilaggr(:),nlaggr(:)
-      type(psb_sspmat_type), intent(out)  :: op_prol
-      integer(psb_ipk_), intent(out)      :: info
-    end subroutine amg_s_hybrid_aggregator_build_tprol
-  end interface
-
-contains
-
-
-  function amg_s_hybrid_aggregator_fmt() result(val)
-    implicit none 
-    character(len=32)  :: val
-
-    val = "Hybrid Decoupled aggregation"
-  end function amg_s_hybrid_aggregator_fmt
-
-  
-end module amg_s_hybrid_aggregator_mod
@@ -123,7 +123,7 @@ contains
    Implicit None

    ! Arguments
-    type(psb_sspmat_type), intent(inout), target        :: a
+    type(psb_sspmat_type), intent(in), target           :: a
    Type(psb_desc_type), Intent(inout)                  :: desc_a 
    class(amg_s_id_solver_type), intent(inout)          :: sv
    integer(psb_ipk_), intent(out)                      :: info
@@ -144,7 +144,7 @@ module amg_s_ilu_solver
           & psb_sspmat_type, psb_s_base_sparse_mat, psb_s_base_vect_type,&
           & psb_ipk_, psb_i_base_vect_type
      implicit none 
-      type(psb_sspmat_type), intent(inout), target       :: a
+      type(psb_sspmat_type), intent(in), target           :: a
      Type(psb_desc_type), Intent(inout)                  :: desc_a 
      class(amg_s_ilu_solver_type), intent(inout)         :: sv
      integer(psb_ipk_), intent(out)                      :: info
@@ -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(:)
@@ -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_base_solver_type), allocatable, intent(inout) :: svout
-      integer(psb_ipk_), intent(out)                            :: info
+      class(amg_s_invk_solver_type), intent(inout) :: sv
+      class(amg_s_base_solver_type), intent(inout) :: svout
+      integer(psb_ipk_), intent(out)                 :: info
    end subroutine amg_s_invk_solver_clone_settings
  end interface

@@ -94,7 +94,7 @@ module amg_s_invk_solver
      Implicit None

      ! Arguments
-      type(psb_sspmat_type), intent(inout), target        :: a
+      type(psb_sspmat_type), intent(in), target           :: a
      Type(psb_desc_type), Intent(inout)                  :: desc_a
      class(amg_s_invk_solver_type), intent(inout)        :: sv
      integer(psb_ipk_), intent(out)                      :: info
@@ -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_base_solver_type), allocatable, intent(inout) :: svout
-      integer(psb_ipk_), intent(out)                            :: info
+      class(amg_s_invt_solver_type), intent(inout) :: sv
+      class(amg_s_base_solver_type), intent(inout) :: svout
+      integer(psb_ipk_), intent(out)                 :: info
    end subroutine amg_s_invt_solver_clone_settings
  end interface

@@ -94,7 +94,7 @@ module amg_s_invt_solver
      Implicit None

      ! Arguments
-      type(psb_sspmat_type), intent(inout), target        :: a
+      type(psb_sspmat_type), intent(in), target           :: a
      Type(psb_desc_type), Intent(inout)                  :: desc_a
      class(amg_s_invt_solver_type), intent(inout)        :: sv
      integer(psb_ipk_), intent(out)                      :: info
@@ -151,7 +151,7 @@ module amg_s_jac_smoother
      import :: psb_desc_type, amg_s_jac_smoother_type, psb_s_vect_type, psb_spk_, &
           & psb_sspmat_type, psb_s_base_sparse_mat, psb_s_base_vect_type,&
           & psb_ipk_, psb_i_base_vect_type
-      type(psb_sspmat_type), intent(inout), target        :: a
+      type(psb_sspmat_type), intent(in), target           :: a
      Type(psb_desc_type), Intent(inout)                  :: desc_a
      class(amg_s_jac_smoother_type), intent(inout)       :: sm
      integer(psb_ipk_), intent(out)                      :: info
@@ -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_base_smoother_type), allocatable, intent(inout) :: smout
+      class(amg_s_jac_smoother_type), intent(inout)  :: sm
+      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
@@ -274,7 +274,7 @@ module amg_s_jac_smoother
      import :: psb_desc_type, amg_s_l1_jac_smoother_type, psb_s_vect_type, &
           & psb_sspmat_type, psb_s_base_sparse_mat, psb_s_base_vect_type,&
           & psb_ipk_, psb_i_base_vect_type
-      type(psb_sspmat_type), intent(inout), target        :: a
+      type(psb_sspmat_type), intent(in), target           :: a
      Type(psb_desc_type), Intent(inout)                  :: desc_a
      class(amg_s_l1_jac_smoother_type), intent(inout)       :: sm
      integer(psb_ipk_), intent(out)                      :: info
@@ -143,7 +143,7 @@ module amg_s_jac_solver
           & psb_sspmat_type, psb_s_base_sparse_mat, psb_s_base_vect_type,&
           & psb_ipk_, psb_i_base_vect_type
      implicit none 
-      type(psb_sspmat_type), intent(inout), target        :: a
+      type(psb_sspmat_type), intent(in), target           :: a
      Type(psb_desc_type), Intent(inout)                  :: desc_a 
      class(amg_s_jac_solver_type), intent(inout)         :: sv
      integer(psb_ipk_), intent(out)                      :: info
@@ -160,7 +160,7 @@ module amg_s_jac_solver
           & psb_sspmat_type, psb_s_base_sparse_mat, psb_s_base_vect_type,&
           & psb_ipk_, psb_i_base_vect_type
      implicit none 
-      type(psb_sspmat_type), intent(inout), target        :: a
+      type(psb_sspmat_type), intent(in), target           :: a
      Type(psb_desc_type), Intent(inout)                  :: desc_a 
      class(amg_s_l1_jac_solver_type), intent(inout)         :: sv
      integer(psb_ipk_), intent(out)                      :: info
@@ -174,7 +174,7 @@ module amg_s_krm_solver
           & psb_sspmat_type, psb_s_base_sparse_mat, psb_s_base_vect_type,&
           & psb_ipk_, psb_i_base_vect_type
      implicit none
-      type(psb_sspmat_type), intent(inout), target        :: a
+      type(psb_sspmat_type), intent(in), target           :: a
      Type(psb_desc_type), Intent(inout)                  :: desc_a
      class(amg_s_krm_solver_type), intent(inout)         :: sv
      integer(psb_ipk_), intent(out)                      :: info
@@ -163,7 +163,7 @@ module amg_s_mumps_solver
      Implicit None

      ! Arguments
-      type(psb_sspmat_type), intent(inout), target        :: a
+      type(psb_sspmat_type), intent(in), target           :: a
      Type(psb_desc_type), Intent(inout)                  :: desc_a 
      class(amg_s_mumps_solver_type), intent(inout)       :: sv
      integer(psb_ipk_), intent(out)                      :: info
@@ -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
@@ -140,7 +140,7 @@ module amg_s_poly_smoother
      import :: psb_desc_type, amg_s_poly_smoother_type, psb_s_vect_type, psb_spk_, &
           & psb_sspmat_type, psb_s_base_sparse_mat, psb_s_base_vect_type,&
           & psb_ipk_, psb_i_base_vect_type
-      type(psb_sspmat_type), intent(inout), target        :: a
+      type(psb_sspmat_type), intent(in), target           :: a
      Type(psb_desc_type), Intent(inout)                  :: desc_a
      class(amg_s_poly_smoother_type), intent(inout)       :: sm
      integer(psb_ipk_), intent(out)                      :: info
@@ -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_base_smoother_type), allocatable, intent(inout) :: smout
+      class(amg_s_poly_smoother_type), intent(inout) :: sm
+      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
@@ -256,7 +256,7 @@ contains
    Implicit None

    ! Arguments
-    type(psb_sspmat_type), intent(inout), target        :: a
+    type(psb_sspmat_type), intent(in), target           :: a
    Type(psb_desc_type), Intent(inout)                  :: desc_a 
    class(amg_s_slu_solver_type), intent(inout)         :: sv
    integer, intent(out)                                :: info
@@ -62,9 +62,6 @@ module amg_z_ainv_solver
    procedure, pass(sv) :: cseti   => amg_z_ainv_solver_cseti
    procedure, pass(sv) :: csetc   => amg_z_ainv_solver_csetc
    procedure, pass(sv) :: csetr   => amg_z_ainv_solver_csetr
-!!$    procedure, pass(sv) :: seti    => amg_z_ainv_solver_seti
-!!$    procedure, pass(sv) :: setc    => amg_z_ainv_solver_setc
-!!$    procedure, pass(sv) :: setr    => amg_z_ainv_solver_setr
    procedure, pass(sv) :: descr   => amg_z_ainv_solver_descr
    procedure, pass(sv) :: default => z_ainv_solver_default
    procedure, nopass   :: stringval  => z_ainv_stringval
@@ -91,9 +88,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_base_solver_type), allocatable, intent(inout) :: svout
-      integer(psb_ipk_), intent(out)                            :: info
+      class(amg_z_ainv_solver_type), intent(inout) :: sv
+      class(amg_z_base_solver_type), intent(inout) :: svout
+      integer(psb_ipk_), intent(out)                 :: info
    end subroutine amg_z_ainv_solver_clone_settings
  end interface

@@ -106,7 +103,7 @@ module amg_z_ainv_solver
      Implicit None

      ! Arguments
-      type(psb_zspmat_type), intent(inout), target        :: a
+      type(psb_zspmat_type), intent(in), target           :: a
      Type(psb_desc_type), Intent(inout)                  :: desc_a
      class(amg_z_ainv_solver_type), intent(inout)        :: sv
      integer(psb_ipk_), intent(out)                      :: info
@@ -230,7 +230,7 @@ module amg_z_as_smoother
           & psb_desc_type, psb_z_base_sparse_mat, psb_ipk_,&
           & psb_i_base_vect_type
      implicit none 
-      type(psb_zspmat_type), intent(inout), target     :: a
+      type(psb_zspmat_type), intent(in), target        :: a
      Type(psb_desc_type), Intent(inout)                 :: desc_a 
      class(amg_z_as_smoother_type), intent(inout)       :: sm
      integer(psb_ipk_), intent(out)                     :: info
@@ -237,7 +237,7 @@ module amg_z_base_smoother_mod
           & psb_z_vect_type, psb_z_base_vect_type, psb_dpk_, &
           & amg_z_base_smoother_type, psb_ipk_, psb_i_base_vect_type
      ! Arguments
-      type(psb_zspmat_type), intent(inout), target  :: a
+      type(psb_zspmat_type), intent(in), target     :: a
      Type(psb_desc_type), Intent(inout)              :: desc_a 
      class(amg_z_base_smoother_type), intent(inout) :: sm 
      integer(psb_ipk_), intent(out)                  :: info
@@ -170,7 +170,7 @@ module amg_z_base_solver_mod
      Implicit None
      
      ! Arguments
-      type(psb_zspmat_type), intent(inout), target          :: a
+      type(psb_zspmat_type), intent(in), target             :: a
      Type(psb_desc_type), Intent(inout)                    :: desc_a 
      class(amg_z_base_solver_type), intent(inout)          :: sv
      integer(psb_ipk_), intent(out)                        :: info
@@ -119,7 +119,7 @@ module amg_z_diag_solver
      import :: psb_desc_type, psb_zspmat_type,  psb_z_base_sparse_mat, &
           & psb_z_vect_type, psb_z_base_vect_type, psb_dpk_, &
           & amg_z_diag_solver_type, psb_ipk_, psb_i_base_vect_type      
-      type(psb_zspmat_type), intent(inout), target        :: a
+      type(psb_zspmat_type), intent(in), target           :: a
      Type(psb_desc_type), Intent(inout)                    :: desc_a 
      class(amg_z_diag_solver_type), intent(inout)        :: sv
      integer(psb_ipk_), intent(out)                        :: info
@@ -331,7 +331,7 @@ module amg_z_l1_diag_solver
      import :: psb_desc_type, psb_zspmat_type,  psb_z_base_sparse_mat, &
           & psb_z_vect_type, psb_z_base_vect_type, psb_dpk_, &
           & amg_z_l1_diag_solver_type, psb_ipk_, psb_i_base_vect_type      
-      type(psb_zspmat_type), intent(inout), target        :: a
+      type(psb_zspmat_type), intent(in), target           :: a
      Type(psb_desc_type), Intent(inout)                    :: desc_a 
      class(amg_z_l1_diag_solver_type), intent(inout)        :: sv
      integer(psb_ipk_), intent(out)                        :: info
@@ -181,7 +181,7 @@ module amg_z_gs_solver
           & psb_zspmat_type, psb_z_base_sparse_mat, psb_z_base_vect_type,&
           & psb_ipk_, psb_i_base_vect_type
      implicit none 
-      type(psb_zspmat_type), intent(inout), target        :: a
+      type(psb_zspmat_type), intent(in), target           :: a
      Type(psb_desc_type), Intent(inout)                  :: desc_a 
      class(amg_z_gs_solver_type), intent(inout)         :: sv
      integer(psb_ipk_), intent(out)                      :: info
@@ -195,7 +195,7 @@ module amg_z_gs_solver
           & psb_zspmat_type, psb_z_base_sparse_mat, psb_z_base_vect_type,&
           & psb_ipk_, psb_i_base_vect_type
      implicit none 
-      type(psb_zspmat_type), intent(inout), target        :: a
+      type(psb_zspmat_type), intent(in), target           :: a
      Type(psb_desc_type), Intent(inout)                  :: desc_a 
      class(amg_z_bwgs_solver_type), intent(inout)         :: sv
      integer(psb_ipk_), intent(out)                      :: info
@@ -1,125 +0,0 @@
-!  
-!   
-!                             AMG4PSBLAS version 1.0
-!    Algebraic Multigrid Package
-!               based on PSBLAS (Parallel Sparse BLAS version 3.7)
-!    
-!    (C) Copyright 2020 
-!  
-!        Salvatore Filippone  
-!        Pasqua D'Ambra   
-!        Fabio Durastante        
-!   
-!    Redistribution and use in source and binary forms, with or without
-!    modification, are permitted provided that the following conditions
-!    are met:
-!      1. Redistributions of source code must retain the above copyright
-!         notice, this list of conditions and the following disclaimer.
-!      2. Redistributions in binary form must reproduce the above copyright
-!         notice, this list of conditions, and the following disclaimer in the
-!         documentation and/or other materials provided with the distribution.
-!      3. The name of the AMG4PSBLAS group or the names of its contributors may
-!         not be used to endorse or promote products derived from this
-!         software without specific written permission.
-!   
-!    THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
-!    ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
-!    TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
-!    PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE AMG4PSBLAS GROUP OR ITS CONTRIBUTORS
-!    BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
-!    CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
-!    SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
-!    INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
-!    CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
-!    ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
-!    POSSIBILITY OF SUCH DAMAGE.
-!   
-!  
-!
-!
-!  The aggregator object hosts the aggregation method for building
-!  the multilevel hierarchy. This variant is based on the hybrid method
-!  presented in 
-!
-!    S. Gratton, P. Henon, P. Jiranek and X. Vasseur:
-!    Reducing complexity of algebraic multigrid by aggregation
-!    Numerical Lin. Algebra with Applications, 2016, 23:501-518
-!    
-module amg_z_hybrid_aggregator_mod
-
-  use amg_z_dec_aggregator_mod
-  !
-  !   sm           -  class(amg_T_base_smoother_type), allocatable
-  !                   The current level preconditioner (aka smoother).
-  !   parms        -  type(amg_RTml_parms)
-  !                   The parameters defining the multilevel strategy.
-  !   ac           -  The local part of the current-level matrix, built by
-  !                   coarsening the previous-level matrix.
-  !   desc_ac      -  type(psb_desc_type).
-  !                   The communication descriptor associated to the matrix
-  !                   stored in ac.
-  !   base_a       -  type(psb_Tspmat_type), pointer.
-  !                   Pointer (really a pointer!) to the local part of the current 
-  !                   matrix (so we have a unified treatment of residuals).
-  !                   We need this to avoid passing explicitly the current matrix
-  !                   to the routine which applies the preconditioner.
-  !   base_desc    -  type(psb_desc_type), pointer.
-  !                   Pointer to the communication descriptor associated to the
-  !                   matrix pointed by base_a.
-  !   map          -  Stores the maps (restriction and prolongation) between the
-  !                   vector spaces associated to the index spaces of the previous
-  !                   and current levels.
-  !
-  !   Methods:  
-  !     Most methods follow the encapsulation hierarchy: they take whatever action
-  !     is appropriate for the current object, then call the corresponding method for
-  !     the contained object.
-  !     As an example: the descr() method prints out a description of the
-  !     level. It starts by invoking the descr() method of the parms object,
-  !     then calls the descr() method of the smoother object. 
-  !
-  !    descr      -   Prints a description of the object.
-  !    default    -   Set default values
-  !    dump       -   Dump to file object contents
-  !    set        -   Sets various parameters; when a request is unknown
-  !                   it is passed to the smoother object for further processing.
-  !    check      -   Sanity checks.
-  !    sizeof     -   Total memory occupation in bytes
-  !    get_nzeros -   Number of nonzeros 
-  !
-  !
-  type, extends(amg_z_dec_aggregator_type) :: amg_z_hybrid_aggregator_type
-    
-  contains
-    procedure, pass(ag) :: bld_tprol    => amg_z_hybrid_aggregator_build_tprol
-    procedure, nopass   :: fmt          => amg_z_hybrid_aggregator_fmt
-  end type amg_z_hybrid_aggregator_type
-
-
-  interface
-    subroutine  amg_z_hybrid_aggregator_build_tprol(ag,parms,a,desc_a,ilaggr,nlaggr,op_prol,info)
-      import :: amg_z_hybrid_aggregator_type, psb_desc_type, psb_zspmat_type, psb_dpk_,  &
-           & psb_ipk_, psb_long_int_k_, amg_dml_parms
-      implicit none
-      class(amg_z_hybrid_aggregator_type), target, intent(inout) :: ag
-      type(amg_dml_parms), intent(inout)  :: parms 
-      type(psb_zspmat_type), intent(in)   :: a
-      type(psb_desc_type), intent(in)     :: desc_a
-      integer(psb_ipk_), allocatable, intent(out) :: ilaggr(:),nlaggr(:)
-      type(psb_zspmat_type), intent(out)  :: op_prol
-      integer(psb_ipk_), intent(out)      :: info
-    end subroutine amg_z_hybrid_aggregator_build_tprol
-  end interface
-
-contains
-
-
-  function amg_z_hybrid_aggregator_fmt() result(val)
-    implicit none 
-    character(len=32)  :: val
-
-    val = "Hybrid Decoupled aggregation"
-  end function amg_z_hybrid_aggregator_fmt
-
-  
-end module amg_z_hybrid_aggregator_mod
@@ -123,7 +123,7 @@ contains
    Implicit None

    ! Arguments
-    type(psb_zspmat_type), intent(inout), target        :: a
+    type(psb_zspmat_type), intent(in), target           :: a
    Type(psb_desc_type), Intent(inout)                  :: desc_a 
    class(amg_z_id_solver_type), intent(inout)          :: sv
    integer(psb_ipk_), intent(out)                      :: info
@@ -144,7 +144,7 @@ module amg_z_ilu_solver
           & psb_zspmat_type, psb_z_base_sparse_mat, psb_z_base_vect_type,&
           & psb_ipk_, psb_i_base_vect_type
      implicit none 
-      type(psb_zspmat_type), intent(inout), target       :: a
+      type(psb_zspmat_type), intent(in), target           :: a
      Type(psb_desc_type), Intent(inout)                  :: desc_a 
      class(amg_z_ilu_solver_type), intent(inout)         :: sv
      integer(psb_ipk_), intent(out)                      :: info
@@ -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(:)
@@ -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_base_solver_type), allocatable, intent(inout) :: svout
-      integer(psb_ipk_), intent(out)                            :: info
+      class(amg_z_invk_solver_type), intent(inout) :: sv
+      class(amg_z_base_solver_type), intent(inout) :: svout
+      integer(psb_ipk_), intent(out)                 :: info
    end subroutine amg_z_invk_solver_clone_settings
  end interface

@@ -94,7 +94,7 @@ module amg_z_invk_solver
      Implicit None

      ! Arguments
-      type(psb_zspmat_type), intent(inout), target        :: a
+      type(psb_zspmat_type), intent(in), target           :: a
      Type(psb_desc_type), Intent(inout)                  :: desc_a
      class(amg_z_invk_solver_type), intent(inout)        :: sv
      integer(psb_ipk_), intent(out)                      :: info
@@ -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_base_solver_type), allocatable, intent(inout) :: svout
-      integer(psb_ipk_), intent(out)                            :: info
+      class(amg_z_invt_solver_type), intent(inout) :: sv
+      class(amg_z_base_solver_type), intent(inout) :: svout
+      integer(psb_ipk_), intent(out)                 :: info
    end subroutine amg_z_invt_solver_clone_settings
  end interface

@@ -94,7 +94,7 @@ module amg_z_invt_solver
      Implicit None

      ! Arguments
-      type(psb_zspmat_type), intent(inout), target        :: a
+      type(psb_zspmat_type), intent(in), target           :: a
      Type(psb_desc_type), Intent(inout)                  :: desc_a
      class(amg_z_invt_solver_type), intent(inout)        :: sv
      integer(psb_ipk_), intent(out)                      :: info
@@ -151,7 +151,7 @@ module amg_z_jac_smoother
      import :: psb_desc_type, amg_z_jac_smoother_type, psb_z_vect_type, psb_dpk_, &
           & psb_zspmat_type, psb_z_base_sparse_mat, psb_z_base_vect_type,&
           & psb_ipk_, psb_i_base_vect_type
-      type(psb_zspmat_type), intent(inout), target        :: a
+      type(psb_zspmat_type), intent(in), target           :: a
      Type(psb_desc_type), Intent(inout)                  :: desc_a
      class(amg_z_jac_smoother_type), intent(inout)       :: sm
      integer(psb_ipk_), intent(out)                      :: info
@@ -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_base_smoother_type), allocatable, intent(inout) :: smout
+      class(amg_z_jac_smoother_type), intent(inout)  :: sm
+      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
@@ -274,7 +274,7 @@ module amg_z_jac_smoother
      import :: psb_desc_type, amg_z_l1_jac_smoother_type, psb_z_vect_type, &
           & psb_zspmat_type, psb_z_base_sparse_mat, psb_z_base_vect_type,&
           & psb_ipk_, psb_i_base_vect_type
-      type(psb_zspmat_type), intent(inout), target        :: a
+      type(psb_zspmat_type), intent(in), target           :: a
      Type(psb_desc_type), Intent(inout)                  :: desc_a
      class(amg_z_l1_jac_smoother_type), intent(inout)       :: sm
      integer(psb_ipk_), intent(out)                      :: info
@@ -143,7 +143,7 @@ module amg_z_jac_solver
           & psb_zspmat_type, psb_z_base_sparse_mat, psb_z_base_vect_type,&
           & psb_ipk_, psb_i_base_vect_type
      implicit none 
-      type(psb_zspmat_type), intent(inout), target        :: a
+      type(psb_zspmat_type), intent(in), target           :: a
      Type(psb_desc_type), Intent(inout)                  :: desc_a 
      class(amg_z_jac_solver_type), intent(inout)         :: sv
      integer(psb_ipk_), intent(out)                      :: info
@@ -160,7 +160,7 @@ module amg_z_jac_solver
           & psb_zspmat_type, psb_z_base_sparse_mat, psb_z_base_vect_type,&
           & psb_ipk_, psb_i_base_vect_type
      implicit none 
-      type(psb_zspmat_type), intent(inout), target        :: a
+      type(psb_zspmat_type), intent(in), target           :: a
      Type(psb_desc_type), Intent(inout)                  :: desc_a 
      class(amg_z_l1_jac_solver_type), intent(inout)         :: sv
      integer(psb_ipk_), intent(out)                      :: info
@@ -174,7 +174,7 @@ module amg_z_krm_solver
           & psb_zspmat_type, psb_z_base_sparse_mat, psb_z_base_vect_type,&
           & psb_ipk_, psb_i_base_vect_type
      implicit none
-      type(psb_zspmat_type), intent(inout), target        :: a
+      type(psb_zspmat_type), intent(in), target           :: a
      Type(psb_desc_type), Intent(inout)                  :: desc_a
      class(amg_z_krm_solver_type), intent(inout)         :: sv
      integer(psb_ipk_), intent(out)                      :: info
@@ -163,7 +163,7 @@ module amg_z_mumps_solver
      Implicit None

      ! Arguments
-      type(psb_zspmat_type), intent(inout), target        :: a
+      type(psb_zspmat_type), intent(in), target           :: a
      Type(psb_desc_type), Intent(inout)                  :: desc_a 
      class(amg_z_mumps_solver_type), intent(inout)       :: sv
      integer(psb_ipk_), intent(out)                      :: info
@@ -256,7 +256,7 @@ contains
    Implicit None

    ! Arguments
-    type(psb_zspmat_type), intent(inout), target        :: a
+    type(psb_zspmat_type), intent(in), target           :: a
    Type(psb_desc_type), Intent(inout)                  :: desc_a 
    class(amg_z_slu_solver_type), intent(inout)         :: sv
    integer, intent(out)                                :: info
@@ -259,7 +259,7 @@ contains
    Implicit None

    ! Arguments
-    type(psb_zspmat_type), intent(inout), target        :: a
+    type(psb_zspmat_type), intent(in), target           :: a
    Type(psb_desc_type), Intent(inout)                  :: desc_a 
    class(amg_z_sludist_solver_type), intent(inout)     :: sv
    integer, intent(out)                                :: info
@@ -259,7 +259,7 @@ contains
    Implicit None

    ! Arguments
-    type(psb_zspmat_type), intent(inout), target        :: a
+    type(psb_zspmat_type), intent(in), target           :: a
    Type(psb_desc_type), Intent(inout)                  :: desc_a 
    class(amg_z_umf_solver_type), intent(inout)         :: sv
    integer, intent(out)                                :: info
@@ -67,14 +67,13 @@ void dMatchBoxPC(MilanLongInt NLVer, MilanLongInt NLEdge,
 #endif


-#define TIME_TRACKER
-    #ifdef TIME_TRACKER
-        double tmr = MPI_Wtime();
-    #endif
+#undef TIME_TRACKER
+#ifdef TIME_TRACKER
+  double tmr = MPI_Wtime();
+#endif

-// Rimosso per tornare al vecchio matching #define OMP
-#ifdef OMP
-				fprintf(stderr,"Warning: using buggy OpenMP matching!\n");
+#if defined(OPENMP)
+	//fprintf(stderr,"Warning: using buggy OpenMP matching!\n");
        dalgoDistEdgeApproxDomEdgesLinearSearchMesgBndlSmallMateCMP(NLVer, NLEdge,
 							   verLocPtr, verLocInd, edgeLocWeight,
 							   verDistance,  Mate,
@@ -93,11 +92,11 @@ void dMatchBoxPC(MilanLongInt NLVer, MilanLongInt NLEdge,
 #endif


-  #ifdef TIME_TRACKER
-    tmr = MPI_Wtime() - tmr;
-    fprintf(stderr, "Elaboration time: %f for %ld nodes\n", tmr, NLVer);
-  #endif
-
+#ifdef TIME_TRACKER
+	tmr = MPI_Wtime() - tmr;
+	fprintf(stderr, "Elaboration time: %f for %ld nodes\n", tmr, NLVer);
+#endif
+	
 #endif
 }

@@ -115,13 +114,24 @@ void sMatchBoxPC(MilanLongInt NLVer, MilanLongInt NLEdge,
  fprintf(stderr,"MatchBoxPC: rank %d nlver %ld nledge %ld [ %ld %ld ]\n",
 	  myRank,NLVer, NLEdge,verDistance[0],verDistance[1]);
 #endif
+#if defined(OPENMP)
+	//fprintf(stderr,"Warning: using buggy OpenMP matching!\n");
+  salgoDistEdgeApproxDomEdgesLinearSearchMesgBndlSmallMateCMP(NLVer, NLEdge,
+							      verLocPtr, verLocInd, edgeLocWeight,
+							      verDistance,  Mate,
+							      myRank, numProcs, C_comm,
+							      msgIndSent, msgActualSent, msgPercent,
+							      ph0_time, ph1_time, ph2_time,
+							      ph1_card, ph2_card );
+#else
  salgoDistEdgeApproxDomEdgesLinearSearchMesgBndlSmallMateC(NLVer, NLEdge,
-							   verLocPtr, verLocInd, edgeLocWeight,
-							   verDistance,  Mate,
-							   myRank, numProcs, C_comm,
-							   msgIndSent, msgActualSent, msgPercent,
-							   ph0_time, ph1_time, ph2_time,
-							   ph1_card, ph2_card );
+							    verLocPtr, verLocInd, edgeLocWeight,
+							    verDistance,  Mate,
+							    myRank, numProcs, C_comm,
+							    msgIndSent, msgActualSent, msgPercent,
+							    ph0_time, ph1_time, ph2_time,
+							    ph1_card, ph2_card );
+#endif
 #endif
 }

@@ -59,7 +59,7 @@
 #include <assert.h>
 #include <map>
 #include <vector>
-#ifdef OMP
+#ifdef OPENMP
 // OpenMP is included and used if and only if the OpenMP version of the matching
 // is required
 #include "omp.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"
@@ -178,264 +181,416 @@ extern "C"
 #define MilanRealMin MINUS_INFINITY
 #endif

-#ifdef OMP
 /* These functions are only used in the experimental OMP implementation, if that
 is disabled there is no reason to actually compile or reference them. */

-    // Function of find the owner of a ghost vertex using binary search:
-    MilanInt findOwnerOfGhost(MilanLongInt vtxIndex, MilanLongInt *mVerDistance,
-                              MilanInt myRank, MilanInt numProcs);
+  // Function of find the owner of a ghost vertex using binary search:
+  MilanInt findOwnerOfGhost(MilanLongInt vtxIndex, MilanLongInt *mVerDistance,
+			    MilanInt myRank, MilanInt numProcs);
+  
+  MilanLongInt firstComputeCandidateMateD(MilanLongInt adj1,
+					 MilanLongInt adj2,
+					 MilanLongInt *verLocInd,
+					 MilanReal *edgeLocWeight);

-    MilanLongInt firstComputeCandidateMate(MilanLongInt adj1,
-                                           MilanLongInt adj2,
-                                           MilanLongInt *verLocInd,
-                                           MilanReal *edgeLocWeight);
+  
+  void queuesTransfer(vector<MilanLongInt> &U,
+		      vector<MilanLongInt> &privateU,
+		      vector<MilanLongInt> &QLocalVtx,
+		      vector<MilanLongInt> &QGhostVtx,
+		      vector<MilanLongInt> &QMsgType,
+		      vector<MilanInt> &QOwner,
+		      vector<MilanLongInt> &privateQLocalVtx,
+		      vector<MilanLongInt> &privateQGhostVtx,
+		      vector<MilanLongInt> &privateQMsgType,
+		      vector<MilanInt> &privateQOwner);
+  
+  bool isAlreadyMatched(MilanLongInt node,
+			MilanLongInt StartIndex,
+			MilanLongInt EndIndex,
+			vector<MilanLongInt> &GMate,
+			MilanLongInt *Mate,
+			map<MilanLongInt, MilanLongInt> &Ghost2LocalMap);
+  
+  MilanLongInt computeCandidateMateD(MilanLongInt adj1,
+				    MilanLongInt adj2,
+				    MilanReal *edgeLocWeight,
+				    MilanLongInt k,
+				    MilanLongInt *verLocInd,
+				    MilanLongInt StartIndex,
+				    MilanLongInt EndIndex,
+				    vector<MilanLongInt> &GMate,
+				    MilanLongInt *Mate,
+				    map<MilanLongInt, MilanLongInt> &Ghost2LocalMap);
+  
+  void initialize(MilanLongInt NLVer, MilanLongInt NLEdge,
+		  MilanLongInt StartIndex, MilanLongInt EndIndex,
+		  MilanLongInt *numGhostEdgesPtr,
+		  MilanLongInt *numGhostVerticesPtr,
+		  MilanLongInt *S,
+		  MilanLongInt *verLocInd,
+		  MilanLongInt *verLocPtr,
+		  map<MilanLongInt, MilanLongInt> &Ghost2LocalMap,
+		  vector<MilanLongInt> &Counter,
+		  vector<MilanLongInt> &verGhostPtr,
+		  vector<MilanLongInt> &verGhostInd,
+		  vector<MilanLongInt> &tempCounter,
+		  vector<MilanLongInt> &GMate,
+		  vector<MilanLongInt> &Message,
+		  vector<MilanLongInt> &QLocalVtx,
+		  vector<MilanLongInt> &QGhostVtx,
+		  vector<MilanLongInt> &QMsgType,
+		  vector<MilanInt> &QOwner,
+		  MilanLongInt *&candidateMate,
+		  vector<MilanLongInt> &U,
+		  vector<MilanLongInt> &privateU,
+		  vector<MilanLongInt> &privateQLocalVtx,
+		  vector<MilanLongInt> &privateQGhostVtx,
+		  vector<MilanLongInt> &privateQMsgType,
+		  vector<MilanInt> &privateQOwner);
+  
+  void clean(MilanLongInt NLVer,
+	     MilanInt myRank,
+	     MilanLongInt MessageIndex,
+	     vector<MPI_Request> &SRequest,
+	     vector<MPI_Status> &SStatus,
+	     MilanInt BufferSize,
+	     MilanLongInt *Buffer,
+	     MilanLongInt msgActual,
+	     MilanLongInt *msgActualSent,
+	     MilanLongInt msgInd,
+	     MilanLongInt *msgIndSent,
+	     MilanLongInt NumMessagesBundled,
+	     MilanReal *msgPercent);
+  
+  void PARALLEL_COMPUTE_CANDIDATE_MATE_BD(MilanLongInt NLVer,
+					 MilanLongInt *verLocPtr,
+					 MilanLongInt *verLocInd,
+					 MilanInt myRank,
+					 MilanReal *edgeLocWeight,
+					 MilanLongInt *candidateMate);
+  
+  void PARALLEL_PROCESS_EXPOSED_VERTEX_BD(MilanLongInt NLVer,
+					 MilanLongInt *candidateMate,
+					 MilanLongInt *verLocInd,
+					 MilanLongInt *verLocPtr,
+					 MilanLongInt StartIndex,
+					 MilanLongInt EndIndex,
+					 MilanLongInt *Mate,
+					 vector<MilanLongInt> &GMate,
+					 map<MilanLongInt, MilanLongInt> &Ghost2LocalMap,
+					 MilanReal *edgeLocWeight,
+					 MilanLongInt *myCardPtr,
+					 MilanLongInt *msgIndPtr,
+					 MilanLongInt *NumMessagesBundledPtr,
+					 MilanLongInt *SPtr,
+					 MilanLongInt *verDistance,
+					 MilanLongInt *PCounter,
+					 vector<MilanLongInt> &Counter,
+					 MilanInt myRank,
+					 MilanInt numProcs,
+					 vector<MilanLongInt> &U,
+					 vector<MilanLongInt> &privateU,
+					 vector<MilanLongInt> &QLocalVtx,
+					 vector<MilanLongInt> &QGhostVtx,
+					 vector<MilanLongInt> &QMsgType,
+					 vector<MilanInt> &QOwner,
+					 vector<MilanLongInt> &privateQLocalVtx,
+					 vector<MilanLongInt> &privateQGhostVtx,
+					 vector<MilanLongInt> &privateQMsgType,
+					 vector<MilanInt> &privateQOwner);
+  
+  void PROCESS_CROSS_EDGE(MilanLongInt *edge,
+			  MilanLongInt *SPtr);
+  
+  void processMatchedVerticesD(
+			      MilanLongInt NLVer,
+			      vector<MilanLongInt> &UChunkBeingProcessed,
+			      vector<MilanLongInt> &U,
+			      vector<MilanLongInt> &privateU,
+			      MilanLongInt StartIndex,
+			      MilanLongInt EndIndex,
+			      MilanLongInt *myCardPtr,
+			      MilanLongInt *msgIndPtr,
+			      MilanLongInt *NumMessagesBundledPtr,
+			      MilanLongInt *SPtr,
+			      MilanLongInt *verLocPtr,
+			      MilanLongInt *verLocInd,
+			      MilanLongInt *verDistance,
+			      MilanLongInt *PCounter,
+			      vector<MilanLongInt> &Counter,
+			      MilanInt myRank,
+			      MilanInt numProcs,
+			      MilanLongInt *candidateMate,
+			      vector<MilanLongInt> &GMate,
+			      MilanLongInt *Mate,
+			      map<MilanLongInt, MilanLongInt> &Ghost2LocalMap,
+			      MilanReal *edgeLocWeight,
+			      vector<MilanLongInt> &QLocalVtx,
+			      vector<MilanLongInt> &QGhostVtx,
+			      vector<MilanLongInt> &QMsgType,
+			      vector<MilanInt> &QOwner,
+			      vector<MilanLongInt> &privateQLocalVtx,
+			      vector<MilanLongInt> &privateQGhostVtx,
+			      vector<MilanLongInt> &privateQMsgType,
+			      vector<MilanInt> &privateQOwner);
+  
+  void processMatchedVerticesAndSendMessagesD(
+					     MilanLongInt NLVer,
+					     vector<MilanLongInt> &UChunkBeingProcessed,
+					     vector<MilanLongInt> &U,
+					     vector<MilanLongInt> &privateU,
+					     MilanLongInt StartIndex,
+					     MilanLongInt EndIndex,
+					     MilanLongInt *myCardPtr,
+					     MilanLongInt *msgIndPtr,
+					     MilanLongInt *NumMessagesBundledPtr,
+					     MilanLongInt *SPtr,
+					     MilanLongInt *verLocPtr,
+					     MilanLongInt *verLocInd,
+					     MilanLongInt *verDistance,
+					     MilanLongInt *PCounter,
+					     vector<MilanLongInt> &Counter,
+					     MilanInt myRank,
+					     MilanInt numProcs,
+					     MilanLongInt *candidateMate,
+					     vector<MilanLongInt> &GMate,
+					     MilanLongInt *Mate,
+					     map<MilanLongInt, MilanLongInt> &Ghost2LocalMap,
+					     MilanReal *edgeLocWeight,
+					     vector<MilanLongInt> &QLocalVtx,
+					     vector<MilanLongInt> &QGhostVtx,
+					     vector<MilanLongInt> &QMsgType,
+					     vector<MilanInt> &QOwner,
+					     vector<MilanLongInt> &privateQLocalVtx,
+					     vector<MilanLongInt> &privateQGhostVtx,
+					     vector<MilanLongInt> &privateQMsgType,
+					     vector<MilanInt> &privateQOwner,
+					     MPI_Comm comm,
+					     MilanLongInt *msgActual,
+					     vector<MilanLongInt> &Message);
+  
+  void sendBundledMessages(MilanLongInt *numGhostEdgesPtr,
+			   MilanInt *BufferSizePtr,
+			   MilanLongInt *Buffer,
+			   vector<MilanLongInt> &PCumulative,
+			   vector<MilanLongInt> &PMessageBundle,
+			   vector<MilanLongInt> &PSizeInfoMessages,
+			   MilanLongInt *PCounter,
+			   MilanLongInt NumMessagesBundled,
+			   MilanLongInt *msgActualPtr,
+			   MilanLongInt *MessageIndexPtr,
+			   MilanInt numProcs,
+			   MilanInt myRank,
+			   MPI_Comm comm,
+			   vector<MilanLongInt> &QLocalVtx,
+			   vector<MilanLongInt> &QGhostVtx,
+			   vector<MilanLongInt> &QMsgType,
+			   vector<MilanInt> &QOwner,
+			   vector<MPI_Request> &SRequest,
+			   vector<MPI_Status> &SStatus);
+  
+  void processMessagesD(
+		       MilanLongInt NLVer,
+		       MilanLongInt *Mate,
+		       MilanLongInt *candidateMate,
+		       map<MilanLongInt, MilanLongInt> &Ghost2LocalMap,
+		       vector<MilanLongInt> &GMate,
+		       vector<MilanLongInt> &Counter,
+		       MilanLongInt StartIndex,
+		       MilanLongInt EndIndex,
+		       MilanLongInt *myCardPtr,
+		       MilanLongInt *msgIndPtr,
+		       MilanLongInt *msgActualPtr,
+		       MilanReal *edgeLocWeight,
+		       MilanLongInt *verDistance,
+		       MilanLongInt *verLocPtr,
+		       MilanLongInt k,
+		       MilanLongInt *verLocInd,
+		       MilanInt numProcs,
+		       MilanInt myRank,
+		       MPI_Comm comm,
+		       vector<MilanLongInt> &Message,
+		       MilanLongInt numGhostEdges,
+		       MilanLongInt u,
+		       MilanLongInt v,
+		       MilanLongInt *SPtr,
+		       vector<MilanLongInt> &U);
+  
+  void extractUChunk(
+		     vector<MilanLongInt> &UChunkBeingProcessed,
+		     vector<MilanLongInt> &U,
+		     vector<MilanLongInt> &privateU);

-    void queuesTransfer(vector<MilanLongInt> &U,
-                        vector<MilanLongInt> &privateU,
-                        vector<MilanLongInt> &QLocalVtx,
-                        vector<MilanLongInt> &QGhostVtx,
-                        vector<MilanLongInt> &QMsgType,
-                        vector<MilanInt> &QOwner,
-                        vector<MilanLongInt> &privateQLocalVtx,
-                        vector<MilanLongInt> &privateQGhostVtx,
-                        vector<MilanLongInt> &privateQMsgType,
-                        vector<MilanInt> &privateQOwner);
+  MilanLongInt firstComputeCandidateMateS(MilanLongInt adj1,
+					 MilanLongInt adj2,
+					 MilanLongInt *verLocInd,
+					  MilanFloat *edgeLocWeight);

-    bool isAlreadyMatched(MilanLongInt node,
-                          MilanLongInt StartIndex,
-                          MilanLongInt EndIndex,
-                          vector<MilanLongInt> &GMate,
-                          MilanLongInt *Mate,
-                          map<MilanLongInt, MilanLongInt> &Ghost2LocalMap);
+  MilanLongInt computeCandidateMateS(MilanLongInt adj1,
+				    MilanLongInt adj2,
+				    MilanFloat *edgeLocWeight,
+				    MilanLongInt k,
+				    MilanLongInt *verLocInd,
+				    MilanLongInt StartIndex,
+				    MilanLongInt EndIndex,
+				    vector<MilanLongInt> &GMate,
+				    MilanLongInt *Mate,
+				    map<MilanLongInt, MilanLongInt> &Ghost2LocalMap);
+  
+  void PARALLEL_COMPUTE_CANDIDATE_MATE_BS(MilanLongInt NLVer,
+					 MilanLongInt *verLocPtr,
+					 MilanLongInt *verLocInd,
+					 MilanInt myRank,
+					 MilanFloat *edgeLocWeight,
+					 MilanLongInt *candidateMate);
+  
+  void PARALLEL_PROCESS_EXPOSED_VERTEX_BS(MilanLongInt NLVer,
+					 MilanLongInt *candidateMate,
+					 MilanLongInt *verLocInd,
+					 MilanLongInt *verLocPtr,
+					 MilanLongInt StartIndex,
+					 MilanLongInt EndIndex,
+					 MilanLongInt *Mate,
+					 vector<MilanLongInt> &GMate,
+					 map<MilanLongInt, MilanLongInt> &Ghost2LocalMap,
+					 MilanFloat *edgeLocWeight,
+					 MilanLongInt *myCardPtr,
+					 MilanLongInt *msgIndPtr,
+					 MilanLongInt *NumMessagesBundledPtr,
+					 MilanLongInt *SPtr,
+					 MilanLongInt *verDistance,
+					 MilanLongInt *PCounter,
+					 vector<MilanLongInt> &Counter,
+					 MilanInt myRank,
+					 MilanInt numProcs,
+					 vector<MilanLongInt> &U,
+					 vector<MilanLongInt> &privateU,
+					 vector<MilanLongInt> &QLocalVtx,
+					 vector<MilanLongInt> &QGhostVtx,
+					 vector<MilanLongInt> &QMsgType,
+					 vector<MilanInt> &QOwner,
+					 vector<MilanLongInt> &privateQLocalVtx,
+					 vector<MilanLongInt> &privateQGhostVtx,
+					 vector<MilanLongInt> &privateQMsgType,
+					 vector<MilanInt> &privateQOwner);
+  void processMatchedVerticesS(
+			      MilanLongInt NLVer,
+			      vector<MilanLongInt> &UChunkBeingProcessed,
+			      vector<MilanLongInt> &U,
+			      vector<MilanLongInt> &privateU,
+			      MilanLongInt StartIndex,
+			      MilanLongInt EndIndex,
+			      MilanLongInt *myCardPtr,
+			      MilanLongInt *msgIndPtr,
+			      MilanLongInt *NumMessagesBundledPtr,
+			      MilanLongInt *SPtr,
+			      MilanLongInt *verLocPtr,
+			      MilanLongInt *verLocInd,
+			      MilanLongInt *verDistance,
+			      MilanLongInt *PCounter,
+			      vector<MilanLongInt> &Counter,
+			      MilanInt myRank,
+			      MilanInt numProcs,
+			      MilanLongInt *candidateMate,
+			      vector<MilanLongInt> &GMate,
+			      MilanLongInt *Mate,
+			      map<MilanLongInt, MilanLongInt> &Ghost2LocalMap,
+			      MilanFloat *edgeLocWeight,
+			      vector<MilanLongInt> &QLocalVtx,
+			      vector<MilanLongInt> &QGhostVtx,
+			      vector<MilanLongInt> &QMsgType,
+			      vector<MilanInt> &QOwner,
+			      vector<MilanLongInt> &privateQLocalVtx,
+			      vector<MilanLongInt> &privateQGhostVtx,
+			      vector<MilanLongInt> &privateQMsgType,
+			      vector<MilanInt> &privateQOwner);
+  
+  void processMatchedVerticesAndSendMessagesS(
+					     MilanLongInt NLVer,
+					     vector<MilanLongInt> &UChunkBeingProcessed,
+					     vector<MilanLongInt> &U,
+					     vector<MilanLongInt> &privateU,
+					     MilanLongInt StartIndex,
+					     MilanLongInt EndIndex,
+					     MilanLongInt *myCardPtr,
+					     MilanLongInt *msgIndPtr,
+					     MilanLongInt *NumMessagesBundledPtr,
+					     MilanLongInt *SPtr,
+					     MilanLongInt *verLocPtr,
+					     MilanLongInt *verLocInd,
+					     MilanLongInt *verDistance,
+					     MilanLongInt *PCounter,
+					     vector<MilanLongInt> &Counter,
+					     MilanInt myRank,
+					     MilanInt numProcs,
+					     MilanLongInt *candidateMate,
+					     vector<MilanLongInt> &GMate,
+					     MilanLongInt *Mate,
+					     map<MilanLongInt, MilanLongInt> &Ghost2LocalMap,
+					     MilanFloat *edgeLocWeight,
+					     vector<MilanLongInt> &QLocalVtx,
+					     vector<MilanLongInt> &QGhostVtx,
+					     vector<MilanLongInt> &QMsgType,
+					     vector<MilanInt> &QOwner,
+					     vector<MilanLongInt> &privateQLocalVtx,
+					     vector<MilanLongInt> &privateQGhostVtx,
+					     vector<MilanLongInt> &privateQMsgType,
+					     vector<MilanInt> &privateQOwner,
+					     MPI_Comm comm,
+					     MilanLongInt *msgActual,
+					     vector<MilanLongInt> &Message);
+  
+  void processMessagesS(
+		       MilanLongInt NLVer,
+		       MilanLongInt *Mate,
+		       MilanLongInt *candidateMate,
+		       map<MilanLongInt, MilanLongInt> &Ghost2LocalMap,
+		       vector<MilanLongInt> &GMate,
+		       vector<MilanLongInt> &Counter,
+		       MilanLongInt StartIndex,
+		       MilanLongInt EndIndex,
+		       MilanLongInt *myCardPtr,
+		       MilanLongInt *msgIndPtr,
+		       MilanLongInt *msgActualPtr,
+		       MilanFloat *edgeLocWeight,
+		       MilanLongInt *verDistance,
+		       MilanLongInt *verLocPtr,
+		       MilanLongInt k,
+		       MilanLongInt *verLocInd,
+		       MilanInt numProcs,
+		       MilanInt myRank,
+		       MPI_Comm comm,
+		       vector<MilanLongInt> &Message,
+		       MilanLongInt numGhostEdges,
+		       MilanLongInt u,
+		       MilanLongInt v,
+		       MilanLongInt *SPtr,
+		       vector<MilanLongInt> &U);
+  
+  
+  void dalgoDistEdgeApproxDomEdgesLinearSearchMesgBndlSmallMateCMP(
+			 MilanLongInt NLVer, MilanLongInt NLEdge,
+			 MilanLongInt *verLocPtr, MilanLongInt *verLocInd, MilanReal *edgeLocWeight,
+			 MilanLongInt *verDistance,
+			 MilanLongInt *Mate,
+			 MilanInt myRank, MilanInt numProcs, MPI_Comm comm,
+			 MilanLongInt *msgIndSent, MilanLongInt *msgActualSent, MilanReal *msgPercent,
+			 MilanReal *ph0_time, MilanReal *ph1_time, MilanReal *ph2_time,
+			 MilanLongInt *ph1_card, MilanLongInt *ph2_card);

-    MilanLongInt computeCandidateMate(MilanLongInt adj1,
-                                      MilanLongInt adj2,
-                                      MilanReal *edgeLocWeight,
-                                      MilanLongInt k,
-                                      MilanLongInt *verLocInd,
-                                      MilanLongInt StartIndex,
-                                      MilanLongInt EndIndex,
-                                      vector<MilanLongInt> &GMate,
-                                      MilanLongInt *Mate,
-                                      map<MilanLongInt, MilanLongInt> &Ghost2LocalMap);
-
-    void initialize(MilanLongInt NLVer, MilanLongInt NLEdge,
-                    MilanLongInt StartIndex, MilanLongInt EndIndex,
-                    MilanLongInt *numGhostEdgesPtr,
-                    MilanLongInt *numGhostVerticesPtr,
-                    MilanLongInt *S,
-                    MilanLongInt *verLocInd,
-                    MilanLongInt *verLocPtr,
-                    map<MilanLongInt, MilanLongInt> &Ghost2LocalMap,
-                    vector<MilanLongInt> &Counter,
-                    vector<MilanLongInt> &verGhostPtr,
-                    vector<MilanLongInt> &verGhostInd,
-                    vector<MilanLongInt> &tempCounter,
-                    vector<MilanLongInt> &GMate,
-                    vector<MilanLongInt> &Message,
-                    vector<MilanLongInt> &QLocalVtx,
-                    vector<MilanLongInt> &QGhostVtx,
-                    vector<MilanLongInt> &QMsgType,
-                    vector<MilanInt> &QOwner,
-                    MilanLongInt *&candidateMate,
-                    vector<MilanLongInt> &U,
-                    vector<MilanLongInt> &privateU,
-                    vector<MilanLongInt> &privateQLocalVtx,
-                    vector<MilanLongInt> &privateQGhostVtx,
-                    vector<MilanLongInt> &privateQMsgType,
-                    vector<MilanInt> &privateQOwner);
-
-    void clean(MilanLongInt NLVer,
-               MilanInt myRank,
-               MilanLongInt MessageIndex,
-               vector<MPI_Request> &SRequest,
-               vector<MPI_Status> &SStatus,
-               MilanInt BufferSize,
-               MilanLongInt *Buffer,
-               MilanLongInt msgActual,
-               MilanLongInt *msgActualSent,
-               MilanLongInt msgInd,
-               MilanLongInt *msgIndSent,
-               MilanLongInt NumMessagesBundled,
-               MilanReal *msgPercent);
-
-    void PARALLEL_COMPUTE_CANDIDATE_MATE_B(MilanLongInt NLVer,
-                                           MilanLongInt *verLocPtr,
-                                           MilanLongInt *verLocInd,
-                                           MilanInt myRank,
-                                           MilanReal *edgeLocWeight,
-                                           MilanLongInt *candidateMate);
-
-    void PARALLEL_PROCESS_EXPOSED_VERTEX_B(MilanLongInt NLVer,
-                                           MilanLongInt *candidateMate,
-                                           MilanLongInt *verLocInd,
-                                           MilanLongInt *verLocPtr,
-                                           MilanLongInt StartIndex,
-                                           MilanLongInt EndIndex,
-                                           MilanLongInt *Mate,
-                                           vector<MilanLongInt> &GMate,
-                                           map<MilanLongInt, MilanLongInt> &Ghost2LocalMap,
-                                           MilanReal *edgeLocWeight,
-                                           MilanLongInt *myCardPtr,
-                                           MilanLongInt *msgIndPtr,
-                                           MilanLongInt *NumMessagesBundledPtr,
-                                           MilanLongInt *SPtr,
-                                           MilanLongInt *verDistance,
-                                           MilanLongInt *PCounter,
-                                           vector<MilanLongInt> &Counter,
-                                           MilanInt myRank,
-                                           MilanInt numProcs,
-                                           vector<MilanLongInt> &U,
-                                           vector<MilanLongInt> &privateU,
-                                           vector<MilanLongInt> &QLocalVtx,
-                                           vector<MilanLongInt> &QGhostVtx,
-                                           vector<MilanLongInt> &QMsgType,
-                                           vector<MilanInt> &QOwner,
-                                           vector<MilanLongInt> &privateQLocalVtx,
-                                           vector<MilanLongInt> &privateQGhostVtx,
-                                           vector<MilanLongInt> &privateQMsgType,
-                                           vector<MilanInt> &privateQOwner);
-
-    void PROCESS_CROSS_EDGE(MilanLongInt *edge,
-                            MilanLongInt *SPtr);
-
-    void processMatchedVertices(
-        MilanLongInt NLVer,
-        vector<MilanLongInt> &UChunkBeingProcessed,
-        vector<MilanLongInt> &U,
-        vector<MilanLongInt> &privateU,
-        MilanLongInt StartIndex,
-        MilanLongInt EndIndex,
-        MilanLongInt *myCardPtr,
-        MilanLongInt *msgIndPtr,
-        MilanLongInt *NumMessagesBundledPtr,
-        MilanLongInt *SPtr,
-        MilanLongInt *verLocPtr,
-        MilanLongInt *verLocInd,
-        MilanLongInt *verDistance,
-        MilanLongInt *PCounter,
-        vector<MilanLongInt> &Counter,
-        MilanInt myRank,
-        MilanInt numProcs,
-        MilanLongInt *candidateMate,
-        vector<MilanLongInt> &GMate,
-        MilanLongInt *Mate,
-        map<MilanLongInt, MilanLongInt> &Ghost2LocalMap,
-        MilanReal *edgeLocWeight,
-        vector<MilanLongInt> &QLocalVtx,
-        vector<MilanLongInt> &QGhostVtx,
-        vector<MilanLongInt> &QMsgType,
-        vector<MilanInt> &QOwner,
-        vector<MilanLongInt> &privateQLocalVtx,
-        vector<MilanLongInt> &privateQGhostVtx,
-        vector<MilanLongInt> &privateQMsgType,
-        vector<MilanInt> &privateQOwner);
-
-    void processMatchedVerticesAndSendMessages(
-        MilanLongInt NLVer,
-        vector<MilanLongInt> &UChunkBeingProcessed,
-        vector<MilanLongInt> &U,
-        vector<MilanLongInt> &privateU,
-        MilanLongInt StartIndex,
-        MilanLongInt EndIndex,
-        MilanLongInt *myCardPtr,
-        MilanLongInt *msgIndPtr,
-        MilanLongInt *NumMessagesBundledPtr,
-        MilanLongInt *SPtr,
-        MilanLongInt *verLocPtr,
-        MilanLongInt *verLocInd,
-        MilanLongInt *verDistance,
-        MilanLongInt *PCounter,
-        vector<MilanLongInt> &Counter,
-        MilanInt myRank,
-        MilanInt numProcs,
-        MilanLongInt *candidateMate,
-        vector<MilanLongInt> &GMate,
-        MilanLongInt *Mate,
-        map<MilanLongInt, MilanLongInt> &Ghost2LocalMap,
-        MilanReal *edgeLocWeight,
-        vector<MilanLongInt> &QLocalVtx,
-        vector<MilanLongInt> &QGhostVtx,
-        vector<MilanLongInt> &QMsgType,
-        vector<MilanInt> &QOwner,
-        vector<MilanLongInt> &privateQLocalVtx,
-        vector<MilanLongInt> &privateQGhostVtx,
-        vector<MilanLongInt> &privateQMsgType,
-        vector<MilanInt> &privateQOwner,
-        MPI_Comm comm,
-        MilanLongInt *msgActual,
-        vector<MilanLongInt> &Message);
-
-    void sendBundledMessages(MilanLongInt *numGhostEdgesPtr,
-                             MilanInt *BufferSizePtr,
-                             MilanLongInt *Buffer,
-                             vector<MilanLongInt> &PCumulative,
-                             vector<MilanLongInt> &PMessageBundle,
-                             vector<MilanLongInt> &PSizeInfoMessages,
-                             MilanLongInt *PCounter,
-                             MilanLongInt NumMessagesBundled,
-                             MilanLongInt *msgActualPtr,
-                             MilanLongInt *MessageIndexPtr,
-                             MilanInt numProcs,
-                             MilanInt myRank,
-                             MPI_Comm comm,
-                             vector<MilanLongInt> &QLocalVtx,
-                             vector<MilanLongInt> &QGhostVtx,
-                             vector<MilanLongInt> &QMsgType,
-                             vector<MilanInt> &QOwner,
-                             vector<MPI_Request> &SRequest,
-                             vector<MPI_Status> &SStatus);
-
-    void processMessages(
-        MilanLongInt NLVer,
-        MilanLongInt *Mate,
-        MilanLongInt *candidateMate,
-        map<MilanLongInt, MilanLongInt> &Ghost2LocalMap,
-        vector<MilanLongInt> &GMate,
-        vector<MilanLongInt> &Counter,
-        MilanLongInt StartIndex,
-        MilanLongInt EndIndex,
-        MilanLongInt *myCardPtr,
-        MilanLongInt *msgIndPtr,
-        MilanLongInt *msgActualPtr,
-        MilanReal *edgeLocWeight,
-        MilanLongInt *verDistance,
-        MilanLongInt *verLocPtr,
-        MilanLongInt k,
-        MilanLongInt *verLocInd,
-        MilanInt numProcs,
-        MilanInt myRank,
-        MPI_Comm comm,
-        vector<MilanLongInt> &Message,
-        MilanLongInt numGhostEdges,
-        MilanLongInt u,
-        MilanLongInt v,
-        MilanLongInt *SPtr,
-        vector<MilanLongInt> &U);
-
-    void extractUChunk(
-        vector<MilanLongInt> &UChunkBeingProcessed,
-        vector<MilanLongInt> &U,
-        vector<MilanLongInt> &privateU);
-
-    void dalgoDistEdgeApproxDomEdgesLinearSearchMesgBndlSmallMateCMP(
-        MilanLongInt NLVer, MilanLongInt NLEdge,
-        MilanLongInt *verLocPtr, MilanLongInt *verLocInd, MilanReal *edgeLocWeight,
-        MilanLongInt *verDistance,
-        MilanLongInt *Mate,
-        MilanInt myRank, MilanInt numProcs, MPI_Comm comm,
-        MilanLongInt *msgIndSent, MilanLongInt *msgActualSent, MilanReal *msgPercent,
-        MilanReal *ph0_time, MilanReal *ph1_time, MilanReal *ph2_time,
-        MilanLongInt *ph1_card, MilanLongInt *ph2_card);
-#endif
+  
+  void salgoDistEdgeApproxDomEdgesLinearSearchMesgBndlSmallMateCMP(
+			 MilanLongInt NLVer, MilanLongInt NLEdge,
+			 MilanLongInt *verLocPtr, MilanLongInt *verLocInd, MilanFloat *edgeLocWeight,
+			 MilanLongInt *verDistance,
+			 MilanLongInt *Mate,
+			 MilanInt myRank, MilanInt numProcs, MPI_Comm comm,
+			 MilanLongInt *msgIndSent, MilanLongInt *msgActualSent, MilanReal *msgPercent,
+			 MilanReal *ph0_time, MilanReal *ph1_time, MilanReal *ph2_time,
+			 MilanLongInt *ph1_card, MilanLongInt *ph2_card);


-#ifndef OMP
-   //Function of find the owner of a ghost vertex using binary search:
-   inline MilanInt findOwnerOfGhost(MilanLongInt vtxIndex, MilanLongInt *mVerDistance,
-                                        MilanInt myRank, MilanInt numProcs);
-#endif

    void dalgoDistEdgeApproxDomEdgesLinearSearchMesgBndlSmallMateC(
        MilanLongInt NLVer, MilanLongInt NLEdge,
@@ -1303,16 +1303,16 @@ void dalgoDistEdgeApproxDomEdgesLinearSearchMesgBndlSmallMateC(
 // SINGLE PRECISION VERSION

 void salgoDistEdgeApproxDomEdgesLinearSearchMesgBndlSmallMateC(
-                                                              MilanLongInt NLVer, MilanLongInt NLEdge,
-                                                              MilanLongInt* verLocPtr, MilanLongInt* verLocInd,
-							      MilanFloat* edgeLocWeight,
-                                                              MilanLongInt* verDistance,
-                                                              MilanLongInt* Mate,
-                                                              MilanInt myRank, MilanInt numProcs, MPI_Comm comm,
-                                                              MilanLongInt* msgIndSent, MilanLongInt* msgActualSent,
-							      MilanReal* msgPercent,
-                                                              MilanReal* ph0_time, MilanReal* ph1_time, MilanReal* ph2_time,
-                                                              MilanLongInt* ph1_card, MilanLongInt* ph2_card ) {
+      MilanLongInt NLVer, MilanLongInt NLEdge,
+      MilanLongInt* verLocPtr, MilanLongInt* verLocInd,
+      MilanFloat* edgeLocWeight,
+      MilanLongInt* verDistance,
+      MilanLongInt* Mate,
+      MilanInt myRank, MilanInt numProcs, MPI_Comm comm,
+      MilanLongInt* msgIndSent, MilanLongInt* msgActualSent,
+      MilanReal* msgPercent,
+      MilanReal* ph0_time, MilanReal* ph1_time, MilanReal* ph2_time,
+      MilanLongInt* ph1_card, MilanLongInt* ph2_card ) {
 #if !defined(SERIAL_MPI)
 #ifdef PRINT_DEBUG_INFO_
    cout<<"\n("<<myRank<<")Within algoEdgeApproxDominatingEdgesLinearSearchMessageBundling()"; fflush(stdout);
@@ -1,5 +1,4 @@
 #include "MatchBoxPC.h"
-#ifdef OMP
 // ***********************************************************************
 //
 //        MatchboxP: A C++ library for approximate weighted matching
@@ -126,8 +125,10 @@ void dalgoDistEdgeApproxDomEdgesLinearSearchMesgBndlSmallMateCMP(
    fflush(stdout);
 #endif

-    MilanLongInt StartIndex = verDistance[myRank];       // The starting vertex owned by the current rank
-    MilanLongInt EndIndex   = verDistance[myRank + 1] - 1; // The ending vertex owned by the current rank
+    // The starting vertex owned by the current rank
+    MilanLongInt StartIndex = verDistance[myRank];
+    // The ending vertex owned by the current rank
+    MilanLongInt EndIndex   = verDistance[myRank + 1] - 1;

    MPI_Status computeStatus;

@@ -145,7 +146,8 @@ void dalgoDistEdgeApproxDomEdgesLinearSearchMesgBndlSmallMateCMP(
    // only one message will be sent in the initialization phase -
    // one of: REQUEST/FAILURE/SUCCESS
    vector<MilanLongInt> QLocalVtx, QGhostVtx, QMsgType;
-    vector<MilanInt> QOwner; // Changed by Fabio to be an integer, addresses needs to be integers!
+    // Changed by Fabio to be an integer, addresses needs to be integers!
+    vector<MilanInt> QOwner; 

    MilanLongInt *PCounter = new MilanLongInt[numProcs];
    for (int i = 0; i < numProcs; i++)
@@ -153,7 +155,8 @@ void dalgoDistEdgeApproxDomEdgesLinearSearchMesgBndlSmallMateCMP(

    MilanLongInt NumMessagesBundled = 0;
    // TODO when the last computational section will be refactored this could be eliminated
-    MilanInt ghostOwner = 0; // Changed by Fabio to be an integer, addresses needs to be integers!
+    // Changed by Fabio to be an integer, addresses needs to be integers!
+    MilanInt ghostOwner = 0; 
    MilanLongInt *candidateMate = nullptr;
 #ifdef PRINT_DEBUG_INFO_
    cout << "\n(" << myRank << ")NV: " << NLVer << "  Edges: " << NLEdge;
@@ -168,9 +171,12 @@ void dalgoDistEdgeApproxDomEdgesLinearSearchMesgBndlSmallMateCMP(
    MilanLongInt myCard = 0;

    // Build the Ghost Vertex Set: Vg
-    map<MilanLongInt, MilanLongInt> Ghost2LocalMap;       // Map each ghost vertex to a local vertex
-    vector<MilanLongInt> Counter;                         // Store the edge count for each ghost vertex
-    MilanLongInt numGhostVertices = 0, numGhostEdges = 0; // Number of Ghost vertices
+    // Map each ghost vertex to a local vertex
+    map<MilanLongInt, MilanLongInt> Ghost2LocalMap;
+    // Store the edge count for each ghost vertex
+    vector<MilanLongInt> Counter;
+    // Number of Ghost vertices
+    MilanLongInt numGhostVertices = 0, numGhostEdges = 0; 

 #ifdef PRINT_DEBUG_INFO_
    cout << "\n(" << myRank << ")About to compute Ghost Vertices...";
@@ -237,7 +243,7 @@ void dalgoDistEdgeApproxDomEdgesLinearSearchMesgBndlSmallMateCMP(
     * PARALLEL_COMPUTE_CANDIDATE_MATE_B is now totally parallel.
     */

-    PARALLEL_COMPUTE_CANDIDATE_MATE_B(NLVer,
+    PARALLEL_COMPUTE_CANDIDATE_MATE_BD(NLVer,
                                      verLocPtr,
                                      verLocInd,
                                      myRank,
@@ -262,7 +268,7 @@ void dalgoDistEdgeApproxDomEdgesLinearSearchMesgBndlSmallMateCMP(
     * TODO: Test when it's actually more efficient to execute this code
     *       in parallel.
     */
-    PARALLEL_PROCESS_EXPOSED_VERTEX_B(NLVer,
+    PARALLEL_PROCESS_EXPOSED_VERTEX_BD(NLVer,
                                      candidateMate,
                                      verLocInd,
                                      verLocPtr,
@@ -314,7 +320,7 @@ void dalgoDistEdgeApproxDomEdgesLinearSearchMesgBndlSmallMateCMP(
    vector<MilanLongInt> UChunkBeingProcessed;
    UChunkBeingProcessed.reserve(UCHUNK);

-    processMatchedVertices(NLVer,
+    processMatchedVerticesD(NLVer,
                           UChunkBeingProcessed,
                           U,
                           privateU,
@@ -423,7 +429,7 @@ void dalgoDistEdgeApproxDomEdgesLinearSearchMesgBndlSmallMateCMP(
        /////////////////////////// PROCESS MATCHED VERTICES //////////////////////////////
        ///////////////////////////////////////////////////////////////////////////////////

-      processMatchedVerticesAndSendMessages(NLVer,
+      processMatchedVerticesAndSendMessagesD(NLVer,
 					    UChunkBeingProcessed,
 					    U,
 					    privateU,
@@ -483,8 +489,8 @@ void dalgoDistEdgeApproxDomEdgesLinearSearchMesgBndlSmallMateCMP(
        ///////////////////////////////////////////////////////////////////////////////////
        /////////////////////////// PROCESS MESSAGES //////////////////////////////////////
        ///////////////////////////////////////////////////////////////////////////////////
-
-        processMessages(NLVer,
+	//startTime = MPI_Wtime();
+        processMessagesD(NLVer,
                        Mate,
                        candidateMate,
                        Ghost2LocalMap,
@@ -549,6 +555,488 @@ void dalgoDistEdgeApproxDomEdgesLinearSearchMesgBndlSmallMateCMP(
    *ph2_card = myCard;                 // Cardinality at the end of Phase-2
 }
 // End of algoDistEdgeApproxDomEdgesLinearSearchMesgBndlSmallMate
+
+void salgoDistEdgeApproxDomEdgesLinearSearchMesgBndlSmallMateCMP(
+    MilanLongInt NLVer, MilanLongInt NLEdge,
+    MilanLongInt *verLocPtr, MilanLongInt *verLocInd,
+    MilanFloat *edgeLocWeight,
+    MilanLongInt *verDistance,
+    MilanLongInt *Mate,
+    MilanInt myRank, MilanInt numProcs, MPI_Comm comm,
+    MilanLongInt *msgIndSent, MilanLongInt *msgActualSent,
+    MilanReal *msgPercent,
+    MilanReal *ph0_time, MilanReal *ph1_time, MilanReal *ph2_time,
+    MilanLongInt *ph1_card, MilanLongInt *ph2_card)
+{
+
+    /*
+     * verDistance: it's a vector long as the number of processors.
+     *              verDistance[i] contains the first node index of the i-th processor
+     *              verDistance[i + 1] contains the last node index of the i-th processor
+     * NLVer: number of elements in the LocPtr
+     * NLEdge: number of edges assigned to the current processor
+     *
+     * Contains the portion of matrix assigned to the processor in
+     * Yale notation
+     * verLocInd: contains the positions on row of the matrix
+     * verLocPtr: i-th value is the position of the first element on the i-th row and
+     *            i+1-th value is the position of the first element on the i+1-th row
+     */
+
+#if !defined(SERIAL_MPI)
+#ifdef PRINT_DEBUG_INFO_
+    cout << "\n(" << myRank << ")Within algoEdgeApproxDominatingEdgesLinearSearchMessageBundling()";
+    fflush(stdout);
+#endif
+
+#ifdef PRINT_DEBUG_INFO_
+    cout << "\n(" << myRank << ") verDistance [" ;
+    for (int i = 0; i < numProcs; i++)
+      cout << verDistance[i] << "," << verDistance[i+1];
+    cout  << "]\n";
+    fflush(stdout);
+#endif
+#ifdef DEBUG_HANG_
+    if (myRank == 0) {
+      cout << "\n(" << myRank << ") verDistance [" ;
+      for (int i = 0; i < numProcs; i++)
+	cout << verDistance[i] << "," ;
+      cout  << verDistance[numProcs]<< "]\n";
+    }
+    fflush(stdout);
+#endif
+
+    // The starting vertex owned by the current rank
+    MilanLongInt StartIndex = verDistance[myRank];
+    // The ending vertex owned by the current rank
+    MilanLongInt EndIndex   = verDistance[myRank + 1] - 1;
+
+    MPI_Status computeStatus;
+
+    MilanLongInt msgActual = 0, msgInd = 0;
+    MilanFloat heaviestEdgeWt = 0.0f; // Assumes positive weight
+    MilanReal startTime, finishTime;
+
+    startTime = MPI_Wtime();
+
+    // Data structures for sending and receiving messages:
+    vector<MilanLongInt> Message; // [ u, v, message_type ]
+    Message.resize(3, -1);
+    // Data structures for Message Bundling:
+    // Although up to two messages can be sent along any cross edge,
+    // only one message will be sent in the initialization phase -
+    // one of: REQUEST/FAILURE/SUCCESS
+    vector<MilanLongInt> QLocalVtx, QGhostVtx, QMsgType;
+    // Changed by Fabio to be an integer, addresses needs to be integers!
+    vector<MilanInt> QOwner; 
+
+    MilanLongInt *PCounter = new MilanLongInt[numProcs];
+    for (int i = 0; i < numProcs; i++)
+        PCounter[i] = 0;
+
+    MilanLongInt NumMessagesBundled = 0;
+    // TODO when the last computational section will be refactored this could be eliminated
+    // Changed by Fabio to be an integer, addresses needs to be integers!
+    MilanInt ghostOwner = 0; 
+    MilanLongInt *candidateMate = nullptr;
+#ifdef PRINT_DEBUG_INFO_
+    cout << "\n(" << myRank << ")NV: " << NLVer << "  Edges: " << NLEdge;
+    fflush(stdout);
+    cout << "\n(" << myRank << ")StartIndex: " << StartIndex << "  EndIndex: " << EndIndex;
+    fflush(stdout);
+#endif
+    // Other Variables:
+    MilanLongInt u = -1, v = -1, w = -1, i = 0;
+    MilanLongInt k = -1, adj1 = -1, adj2 = -1;
+    MilanLongInt k1 = -1, adj11 = -1, adj12 = -1;
+    MilanLongInt myCard = 0;
+
+    // Build the Ghost Vertex Set: Vg
+    // Map each ghost vertex to a local vertex
+    map<MilanLongInt, MilanLongInt> Ghost2LocalMap;
+    // Store the edge count for each ghost vertex
+    vector<MilanLongInt> Counter;
+    // Number of Ghost vertices
+    MilanLongInt numGhostVertices = 0, numGhostEdges = 0; 
+
+#ifdef PRINT_DEBUG_INFO_
+    cout << "\n(" << myRank << ")About to compute Ghost Vertices...";
+    fflush(stdout);
+#endif
+#ifdef DEBUG_HANG_
+    if (myRank == 0)
+        cout << "\n(" << myRank << ")About to compute Ghost Vertices...";
+    fflush(stdout);
+#endif
+
+    // Define Adjacency Lists for Ghost Vertices:
+    // cout<<"Building Ghost data structures ... \n\n";
+    vector<MilanLongInt> verGhostPtr, verGhostInd, tempCounter;
+    // Mate array for ghost vertices:
+    vector<MilanLongInt> GMate; // Proportional to the number of ghost vertices
+    MilanLongInt S;
+    MilanLongInt privateMyCard = 0;
+    vector<MilanLongInt> PCumulative, PMessageBundle, PSizeInfoMessages;
+    vector<MPI_Request> SRequest;  // Requests that are used for each send message
+    vector<MPI_Status> SStatus;    // Status of sent messages, used in MPI_Wait
+    MilanLongInt MessageIndex = 0; // Pointer for current message
+    MilanInt BufferSize;
+    MilanLongInt *Buffer;
+
+    vector<MilanLongInt> privateQLocalVtx, privateQGhostVtx, privateQMsgType;
+    vector<MilanInt> privateQOwner;
+    vector<MilanLongInt> U, privateU;
+
+    
+    initialize(NLVer, NLEdge, StartIndex,
+               EndIndex, &numGhostEdges,
+               &numGhostVertices, &S,
+               verLocInd, verLocPtr,
+               Ghost2LocalMap, Counter,
+               verGhostPtr, verGhostInd,
+               tempCounter, GMate,
+               Message, QLocalVtx,
+               QGhostVtx, QMsgType, QOwner,
+               candidateMate, U,
+               privateU,
+               privateQLocalVtx,
+               privateQGhostVtx,
+               privateQMsgType,
+               privateQOwner);
+
+    finishTime = MPI_Wtime();
+    *ph0_time = finishTime - startTime; // Time taken for Phase-0: Initialization
+#ifdef DEBUG_HANG_
+    cout << myRank << " Finished initialization" << endl;
+    fflush(stdout);
+#endif
+
+    startTime = MPI_Wtime();
+
+    /////////////////////////////////////////////////////////////////////////////////////////
+    //////////////////////////////////// INITIALIZATION /////////////////////////////////////
+    /////////////////////////////////////////////////////////////////////////////////////////
+    // Compute the Initial Matching Set:
+
+    /*
+     * OMP PARALLEL_COMPUTE_CANDIDATE_MATE_B has been splitted from
+     * PARALLEL_PROCESS_EXPOSED_VERTEX_B in order to better parallelize
+     * the two.
+     * PARALLEL_COMPUTE_CANDIDATE_MATE_B is now totally parallel.
+     */
+
+    PARALLEL_COMPUTE_CANDIDATE_MATE_BS(NLVer,
+                                      verLocPtr,
+                                      verLocInd,
+                                      myRank,
+                                      edgeLocWeight,
+                                      candidateMate);
+
+#ifdef DEBUG_HANG_
+    cout << myRank << " Finished Exposed Vertex" << endl;
+    fflush(stdout);
+#if 0
+    cout << myRank << " candidateMate after parallelCompute " <<endl;
+    for (int i=0; i<NLVer; i++) {
+      cout << candidateMate[i] << " " ;
+    }
+    cout << endl;
+#endif
+#endif
+    /*
+     * PARALLEL_PROCESS_EXPOSED_VERTEX_B
+     * TODO: write comment
+     *
+     * TODO: Test when it's actually more efficient to execute this code
+     *       in parallel.
+     */
+    PARALLEL_PROCESS_EXPOSED_VERTEX_BS(NLVer,
+                                      candidateMate,
+                                      verLocInd,
+                                      verLocPtr,
+                                      StartIndex,
+                                      EndIndex,
+                                      Mate,
+                                      GMate,
+                                      Ghost2LocalMap,
+                                      edgeLocWeight,
+                                      &myCard,
+                                      &msgInd,
+                                      &NumMessagesBundled,
+                                      &S,
+                                      verDistance,
+                                      PCounter,
+                                      Counter,
+                                      myRank,
+                                      numProcs,
+                                      U,
+                                      privateU,
+                                      QLocalVtx,
+                                      QGhostVtx,
+                                      QMsgType,
+                                      QOwner,
+                                      privateQLocalVtx,
+                                      privateQGhostVtx,
+                                      privateQMsgType,
+                                      privateQOwner);
+
+    tempCounter.clear(); // Do not need this any more
+
+#ifdef DEBUG_HANG_
+    cout << myRank << " Finished Exposed Vertex" << endl;
+    fflush(stdout);
+#if 0
+    cout << myRank << " Mate after Exposed Vertices " <<endl;
+    for (int i=0; i<NLVer; i++) {
+      cout << Mate[i] << " " ;
+    }
+    cout << endl;
+#endif
+#endif
+
+    ///////////////////////////////////////////////////////////////////////////////////
+    /////////////////////////// PROCESS MATCHED VERTICES //////////////////////////////
+    ///////////////////////////////////////////////////////////////////////////////////
+
+    // TODO what would be the optimal UCHUNK
+    vector<MilanLongInt> UChunkBeingProcessed;
+    UChunkBeingProcessed.reserve(UCHUNK);
+
+    processMatchedVerticesS(NLVer,
+                           UChunkBeingProcessed,
+                           U,
+                           privateU,
+                           StartIndex,
+                           EndIndex,
+                           &myCard,
+                           &msgInd,
+                           &NumMessagesBundled,
+                           &S,
+                           verLocPtr,
+                           verLocInd,
+                           verDistance,
+                           PCounter,
+                           Counter,
+                           myRank,
+                           numProcs,
+                           candidateMate,
+                           GMate,
+                           Mate,
+                           Ghost2LocalMap,
+                           edgeLocWeight,
+                           QLocalVtx,
+                           QGhostVtx,
+                           QMsgType,
+                           QOwner,
+                           privateQLocalVtx,
+                           privateQGhostVtx,
+                           privateQMsgType,
+                           privateQOwner);
+
+
+#ifdef DEBUG_HANG_
+    cout << myRank << " Finished Process Vertices" << endl;
+    fflush(stdout);
+#if 0
+    cout << myRank << " Mate after Matched Vertices " <<endl;
+    for (int i=0; i<NLVer; i++) {
+      cout << Mate[i] << " " ;
+    }
+    cout << endl;
+#endif
+#endif
+
+    /////////////////////////////////////////////////////////////////////////////////////////
+    ///////////////////////////// SEND BUNDLED MESSAGES /////////////////////////////////////
+    /////////////////////////////////////////////////////////////////////////////////////////
+
+    sendBundledMessages(&numGhostEdges,
+                        &BufferSize,
+                        Buffer,
+                        PCumulative,
+                        PMessageBundle,
+                        PSizeInfoMessages,
+                        PCounter,
+                        NumMessagesBundled,
+                        &msgActual,
+                        &MessageIndex,
+                        numProcs,
+                        myRank,
+                        comm,
+                        QLocalVtx,
+                        QGhostVtx,
+                        QMsgType,
+                        QOwner,
+                        SRequest,
+                        SStatus);
+
+    ///////////////////////// END OF SEND BUNDLED MESSAGES //////////////////////////////////
+
+    finishTime = MPI_Wtime();
+    *ph1_time = finishTime - startTime; // Time taken for Phase-1
+
+#ifdef DEBUG_HANG_
+    cout << myRank << " Finished sendBundles" << endl;
+    fflush(stdout);
+#endif
+
+    *ph1_card = myCard;                 // Cardinality at the end of Phase-1
+    startTime = MPI_Wtime();
+    /////////////////////////////////////////////////////////////////////////////////////////
+    //////////////////////////////////////// MAIN LOOP //////////////////////////////////////
+    /////////////////////////////////////////////////////////////////////////////////////////
+    // Main While Loop:
+#ifdef PRINT_DEBUG_INFO_
+    cout << "\n(" << myRank << "=========================************===============================" << endl;
+    fflush(stdout);
+    fflush(stdout);
+#endif
+#ifdef PRINT_DEBUG_INFO_
+    cout << "\n(" << myRank << ")Entering While(true) loop..";
+    fflush(stdout);
+#endif
+#ifdef PRINT_DEBUG_INFO_
+    cout << "\n(" << myRank << "=========================************===============================" << endl;
+    fflush(stdout);
+    fflush(stdout);
+#endif
+
+    while (true)  {
+#ifdef DEBUG_HANG_
+      //if (myRank == 0)
+      cout << "\n(" << myRank << ") Main loop" << endl;
+      fflush(stdout);
+#endif
+        ///////////////////////////////////////////////////////////////////////////////////
+        /////////////////////////// PROCESS MATCHED VERTICES //////////////////////////////
+        ///////////////////////////////////////////////////////////////////////////////////
+
+      processMatchedVerticesAndSendMessagesS(NLVer,
+					    UChunkBeingProcessed,
+					    U,
+					    privateU,
+					    StartIndex,
+					    EndIndex,
+					    &myCard,
+					    &msgInd,
+					    &NumMessagesBundled,
+					    &S,
+					    verLocPtr,
+					    verLocInd,
+					    verDistance,
+					    PCounter,
+					    Counter,
+					    myRank,
+					    numProcs,
+					    candidateMate,
+					    GMate,
+					    Mate,
+					    Ghost2LocalMap,
+					    edgeLocWeight,
+					    QLocalVtx,
+					    QGhostVtx,
+					    QMsgType,
+					    QOwner,
+					    privateQLocalVtx,
+					    privateQGhostVtx,
+					    privateQMsgType,
+					    privateQOwner,
+					    comm,
+					    &msgActual,
+					    Message);
+
+      ///////////////////////// END OF PROCESS MATCHED VERTICES /////////////////////////
+
+        //// BREAK IF NO MESSAGES EXPECTED /////////
+#ifdef DEBUG_HANG_
+#if 0
+      cout << myRank << " Mate after ProcessMatchedAndSend phase "<<S <<endl;
+      for (int i=0; i<NLVer; i++) {
+	cout << Mate[i] << " " ;
+      }
+      cout << endl;
+#endif
+#endif
+#ifdef PRINT_DEBUG_INFO_
+        cout << "\n(" << myRank << ")Deciding whether to break: S= " << S << endl;
+#endif
+
+        if (S == 0) {
+#ifdef DEBUG_HANG_
+	  cout << "\n(" << myRank << ") Breaking out" << endl;
+	  fflush(stdout);
+#endif
+	  break;
+        }
+        ///////////////////////////////////////////////////////////////////////////////////
+        /////////////////////////// PROCESS MESSAGES //////////////////////////////////////
+        ///////////////////////////////////////////////////////////////////////////////////
+
+        processMessagesS(NLVer,
+                        Mate,
+                        candidateMate,
+                        Ghost2LocalMap,
+                        GMate,
+                        Counter,
+                        StartIndex,
+                        EndIndex,
+                        &myCard,
+                        &msgInd,
+                        &msgActual,
+                        edgeLocWeight,
+                        verDistance,
+                        verLocPtr,
+                        k,
+                        verLocInd,
+                        numProcs,
+                        myRank,
+                        comm,
+                        Message,
+                        numGhostEdges,
+                        u,
+                        v,
+                        &S,
+                        U);
+
+        ///////////////////////// END OF PROCESS MESSAGES /////////////////////////////////
+#ifdef DEBUG_HANG_
+#if 0
+      cout << myRank << " Mate after ProcessMessages phase "<<S <<endl;
+      for (int i=0; i<NLVer; i++) {
+	cout << Mate[i] << " " ;
+      }
+      cout << endl;
+#endif
+#endif
+#ifdef PRINT_DEBUG_INFO_
+        cout << "\n(" << myRank << ")Finished Message processing phase: S= " << S;
+        fflush(stdout);
+        cout << "\n(" << myRank << ")** SENT     : ACTUAL= " << msgActual;
+        fflush(stdout);
+        cout << "\n(" << myRank << ")** SENT     : INDIVIDUAL= " << msgInd << endl;
+        fflush(stdout);
+#endif
+    } // End of while (true)
+
+    clean(NLVer,
+          myRank,
+          MessageIndex,
+          SRequest,
+          SStatus,
+          BufferSize,
+          Buffer,
+          msgActual,
+          msgActualSent,
+          msgInd,
+          msgIndSent,
+          NumMessagesBundled,
+          msgPercent);
+
+    finishTime = MPI_Wtime();
+    *ph2_time = finishTime - startTime; // Time taken for Phase-2
+    *ph2_card = myCard;                 // Cardinality at the end of Phase-2
+}
+
 #endif

 #endif
@@ -177,23 +177,24 @@ subroutine  amg_c_dec_aggregator_mat_bld(ag,parms,a,desc_a,ilaggr,nlaggr,&
  select case (parms%aggr_prol)
  case (amg_no_smooth_) 

-    call amg_caggrmat_nosmth_bld(a,desc_a,ilaggr,nlaggr,&
-         & parms,ac,desc_ac,op_prol,op_restr,t_prol,info)
+    call amg_caggrmat_nosmth_bld(parms%aggr_prol,a,desc_a,ilaggr,&
+            nlaggr,parms,ac,desc_ac,op_prol,op_restr,t_prol,info)

-  case(amg_smooth_prol_) 
+  case(amg_smooth_prol_,amg_l1_smooth_prol_) 

-    call amg_caggrmat_smth_bld(a,desc_a,ilaggr,nlaggr, &
-         & parms,ac,desc_ac,op_prol,op_restr,t_prol,info)
+    call amg_caggrmat_smth_bld(parms%aggr_prol,a,desc_a,&
+            ilaggr,nlaggr,parms,ac,desc_ac,op_prol,&
+            op_restr,t_prol,info)

 !!$  case(amg_biz_prol_) 
 !!$
 !!$    call amg_caggrmat_biz_bld(a,desc_a,ilaggr,nlaggr, &
 !!$         & parms,ac,desc_ac,op_prol,op_restr,t_prol,info)
-
+  
  case(amg_min_energy_) 

-    call amg_caggrmat_minnrg_bld(a,desc_a,ilaggr,nlaggr, &
-         & parms,ac,desc_ac,op_prol,op_restr,t_prol,info)
+    call amg_caggrmat_minnrg_bld(parms%aggr_prol,a,desc_a,ilaggr,&
+            nlaggr,parms,ac,desc_ac,op_prol,op_restr,t_prol,info)

  case default
    info = psb_err_internal_error_
@@ -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()

@@ -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()

@@ -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
@@ -275,7 +275,8 @@ subroutine amg_c_soc1_map_bld(iorder,theta,clean_zeros,a,desc_a,nlaggr,ilaggr,in
            disjoint = all(ilaggr(icol(1:ip)) == -(nr+1)).or.(ip==0)
            if (disjoint) then
              locnaggr(kk)     = locnaggr(kk) + 1
-              itmp = (bnds(kk)-1+locnaggr(kk))*nths+kk
+              itmp = (bnds(kk)-1+locnaggr(kk)) !be careful about overflow
+              itmp = itmp*nths+kk
              if (itmp < (bnds(kk)-1+locnaggr(kk))) then
                !$omp atomic update
                info = max(12345678,info)
@@ -356,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
@@ -544,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
-
@@ -309,7 +309,8 @@ subroutine amg_c_soc2_map_bld(iorder,theta,clean_zeros,a,desc_a,nlaggr,ilaggr,in
            !
            if (disjoint) then
              locnaggr(kk)     = locnaggr(kk) + 1
-              itmp = (bnds(kk)-1+locnaggr(kk))*nths+kk
+              itmp = (bnds(kk)-1+locnaggr(kk)) !be careful about overflow
+              itmp = itmp*nths+kk
              if (itmp < (bnds(kk)-1+locnaggr(kk))) then
                !$omp atomic update
                info = max(12345678,info)
@@ -69,6 +69,7 @@
 !
 !
 ! Arguments:
+!    dol1smoothing - fictitious integer argument, it is not used inside
 !    a          -  type(psb_cspmat_type), input.     
 !                  The sparse matrix structure containing the local part of
 !                  the fine-level matrix.
@@ -104,8 +105,8 @@
 !                  Error code.
 !
 !
-subroutine amg_caggrmat_minnrg_bld(a,desc_a,ilaggr,nlaggr,parms,&
-     & ac,desc_ac,op_prol,op_restr,t_prol,info)
+subroutine amg_caggrmat_minnrg_bld(dol1smoothing,a,desc_a,ilaggr,nlaggr,&
+                parms,ac,desc_ac,op_prol,op_restr,t_prol,info)
  use psb_base_mod
  use amg_base_prec_type
  use amg_c_inner_mod, amg_protect_name => amg_caggrmat_minnrg_bld
@@ -113,6 +114,7 @@ subroutine amg_caggrmat_minnrg_bld(a,desc_a,ilaggr,nlaggr,parms,&
  implicit none

  ! Arguments
+  integer(psb_ipk_), intent(in)                 :: dol1smoothing
  type(psb_cspmat_type), intent(in)           :: a
  type(psb_desc_type), intent(inout)            :: desc_a
  integer(psb_lpk_), intent(inout)              :: ilaggr(:), nlaggr(:)
@@ -171,6 +173,13 @@ subroutine amg_caggrmat_minnrg_bld(a,desc_a,ilaggr,nlaggr,parms,&

  filter_mat = (parms%aggr_filter == amg_filter_mat_)

+  if (dol1smoothing.ne.amg_no_smooth_) then
+    info=psb_err_fatal_; 
+    call psb_errpush(info,name,a_err='Are you trying to smooth an unsmoothed aggregation?')
+    goto 9999
+  end if
+
+
  !NEEDS TO BE REWORKED !!
  
  ! naggr: number of local aggregates
@@ -300,7 +309,7 @@ subroutine amg_caggrmat_minnrg_bld(a,desc_a,ilaggr,nlaggr,parms,&
 !!$        endif
 !!$      enddo
 !!$      if (jd == -1) then 
-!!$        write(0,*) 'Wrong input: we need the diagonal!!!!', i
+!!$                write(0,*) name,': Warning: there is no diagonal element', i
 !!$      else
 !!$        acsrf%val(jd)=acsrf%val(jd)-tmp
 !!$      end if
@@ -94,10 +94,11 @@
 !               
 !    info       -  integer, output.
 !                  Error code.
+!    dol1smoothing - optional, this is here just for interfacing reasons. It is not used by the
+!                    code
 !
-!
-subroutine amg_caggrmat_nosmth_bld(a,desc_a,ilaggr,nlaggr,parms,&
-     & ac,desc_ac,op_prol,op_restr,t_prol,info)
+subroutine amg_caggrmat_nosmth_bld(dol1smoothing,a,desc_a,ilaggr,nlaggr,&
+                parms,ac,desc_ac,op_prol,op_restr,t_prol,info)
  use psb_base_mod
  use amg_base_prec_type
  use amg_c_inner_mod, amg_protect_name => amg_caggrmat_nosmth_bld
@@ -105,6 +106,7 @@ subroutine amg_caggrmat_nosmth_bld(a,desc_a,ilaggr,nlaggr,parms,&
  implicit none

  ! Arguments
+  integer(psb_ipk_), intent(in)              :: dol1smoothing 
  type(psb_cspmat_type), intent(in)        :: a
  type(psb_desc_type), intent(inout)         :: desc_a
  integer(psb_lpk_), intent(inout)           :: ilaggr(:), nlaggr(:)
@@ -137,6 +139,12 @@ subroutine amg_caggrmat_nosmth_bld(a,desc_a,ilaggr,nlaggr,parms,&

  ctxt = desc_a%get_context()
  call psb_info(ctxt, me, np)
+  if (dol1smoothing.ne.amg_no_smooth_) then
+    info=psb_err_fatal_;         
+    call psb_errpush(info,name,a_err='Are you trying to smooth an unsmoothed aggregation?')
+    goto 9999
+  end if
+
  nglob = desc_a%get_global_rows()
  nrow  = desc_a%get_local_rows()
  ncol  = desc_a%get_local_cols()
@@ -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,&
-     & ac,desc_ac,op_prol,op_restr,t_prol,info)
+subroutine amg_caggrmat_smth_bld(dol1smoothing,a,desc_a,ilaggr,nlaggr,&
+                parms,ac,desc_ac,op_prol,op_restr,t_prol,info)
  use psb_base_mod
  use amg_base_prec_type
  use amg_c_inner_mod, amg_protect_name => amg_caggrmat_smth_bld
  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

@@ -177,23 +177,24 @@ subroutine  amg_d_dec_aggregator_mat_bld(ag,parms,a,desc_a,ilaggr,nlaggr,&
  select case (parms%aggr_prol)
  case (amg_no_smooth_) 

-    call amg_daggrmat_nosmth_bld(a,desc_a,ilaggr,nlaggr,&
-         & parms,ac,desc_ac,op_prol,op_restr,t_prol,info)
+    call amg_daggrmat_nosmth_bld(parms%aggr_prol,a,desc_a,ilaggr,&
+            nlaggr,parms,ac,desc_ac,op_prol,op_restr,t_prol,info)

-  case(amg_smooth_prol_) 
+  case(amg_smooth_prol_,amg_l1_smooth_prol_) 

-    call amg_daggrmat_smth_bld(a,desc_a,ilaggr,nlaggr, &
-         & parms,ac,desc_ac,op_prol,op_restr,t_prol,info)
+    call amg_daggrmat_smth_bld(parms%aggr_prol,a,desc_a,&
+            ilaggr,nlaggr,parms,ac,desc_ac,op_prol,&
+            op_restr,t_prol,info)

 !!$  case(amg_biz_prol_) 
 !!$
 !!$    call amg_daggrmat_biz_bld(a,desc_a,ilaggr,nlaggr, &
 !!$         & parms,ac,desc_ac,op_prol,op_restr,t_prol,info)
-
+  
  case(amg_min_energy_) 

-    call amg_daggrmat_minnrg_bld(a,desc_a,ilaggr,nlaggr, &
-         & parms,ac,desc_ac,op_prol,op_restr,t_prol,info)
+    call amg_daggrmat_minnrg_bld(parms%aggr_prol,a,desc_a,ilaggr,&
+            nlaggr,parms,ac,desc_ac,op_prol,op_restr,t_prol,info)

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

-  case(amg_smooth_prol_)
-    call amg_d_parmatch_smth_bld(ag,a,desc_a,ilaggr,nlaggr,parms,&
-         & ac,desc_ac,op_prol,op_restr,t_prol,info)
+  case(amg_smooth_prol_,amg_l1_smooth_prol_)
+    call amg_d_parmatch_smth_bld(parms%aggr_prol,ag,a,desc_a,&
+            ilaggr,nlaggr,parms,ac,desc_ac,op_prol,op_restr,&
+            t_prol,info)

 !!$  case(amg_biz_prol_)
 !!$    call amg_daggrmat_biz_bld(a,desc_a,ilaggr,nlaggr, &
 !!$         & parms,ac,desc_ac,op_prol,op_restr,info)

  case(amg_min_energy_)
-    call amg_daggrmat_minnrg_bld(a,desc_a,ilaggr,nlaggr, &
-         & parms,ac,desc_ac,op_prol,op_restr,t_prol,info)
+    call amg_daggrmat_minnrg_bld(parms%aggr_prol,a,desc_a,&
+            ilaggr,nlaggr,parms,ac,desc_ac,op_prol,op_restr,&
+            t_prol,info)

  case default
    info = psb_err_internal_error_
@@ -69,6 +69,8 @@
 !
 !
 ! Arguments:
+!    dol1smoothing - Select between l1-Jacobi and Jacobi as smoother for the
+!                  tentative prolongator
 !    a          -  type(psb_dspmat_type), input.
 !                  The sparse matrix structure containing the local part of
 !                  the fine-level matrix.
@@ -102,8 +104,8 @@
 !    info       -  integer, output.
 !                  Error code.
 !
-subroutine amg_d_parmatch_smth_bld(ag,a,desc_a,ilaggr,nlaggr,parms,&
-     & ac,desc_ac,op_prol,op_restr,t_prol,info)
+subroutine amg_d_parmatch_smth_bld(dol1smoothing,ag,a,desc_a,ilaggr,nlaggr,&
+                parms,ac,desc_ac,op_prol,op_restr,t_prol,info)
  use psb_base_mod
  use amg_base_prec_type
  use amg_d_inner_mod
@@ -116,6 +118,7 @@ subroutine amg_d_parmatch_smth_bld(ag,a,desc_a,ilaggr,nlaggr,parms,&
  implicit none

  ! Arguments
+  integer(psb_ipk_), intent(in)            :: dol1smoothing
  class(amg_d_parmatch_aggregator_type), target, intent(inout) :: ag
  type(psb_dspmat_type), intent(in)      :: a
  type(psb_desc_type), intent(inout)     :: desc_a
@@ -137,7 +140,7 @@ subroutine amg_d_parmatch_smth_bld(ag,a,desc_a,ilaggr,nlaggr,parms,&
  type(psb_d_coo_sparse_mat)  :: coo_prol, coo_restr
  type(psb_d_csr_sparse_mat)  :: acsrf, csr_prol, acsr, tcsr
  real(psb_dpk_), allocatable :: adiag(:)
-  real(psb_dpk_), allocatable :: arwsum(:)
+  real(psb_dpk_), allocatable :: arwsum(:),l1rwsum(:)
  logical            :: filter_mat
  integer(psb_ipk_)            :: debug_level, debug_unit, err_act
  integer(psb_ipk_), parameter :: ncmax=16
@@ -145,6 +148,7 @@ subroutine amg_d_parmatch_smth_bld(ag,a,desc_a,ilaggr,nlaggr,parms,&
  logical, parameter :: debug_new=.false., dump_r=.false., dump_p=.false., debug=.false.
  character(len=80) :: filename
  logical, parameter :: do_timings=.false.
+  logical :: do_l1correction=.false.
  integer(psb_ipk_), save :: idx_spspmm=-1, idx_phase1=-1, idx_gtrans=-1, idx_phase2=-1, idx_refine=-1, idx_phase3=-1
  integer(psb_ipk_), save :: idx_cdasb=-1, idx_ptap=-1

@@ -166,6 +170,10 @@ subroutine amg_d_parmatch_smth_bld(ag,a,desc_a,ilaggr,nlaggr,parms,&
  ncol  = desc_a%get_local_cols()

  theta = parms%aggr_thresh
+  ! Check if we have to perform l1-Jacobi or Jacobi as smoother
+  if(dol1smoothing.eq.amg_l1_smooth_prol_) do_l1correction=.true.
+
+
  !write(0,*) me,' ',trim(name),' Start ',idx_spspmm
  if ((do_timings).and.(idx_spspmm==-1)) &
       & idx_spspmm = psb_get_timer_idx("PMC_SMTH_BLD: par_spspmm")
@@ -217,6 +225,19 @@ subroutine amg_d_parmatch_smth_bld(ag,a,desc_a,ilaggr,nlaggr,parms,&
  if (info == psb_success_) &
       & call psb_halo(adiag,desc_a,info)
  if (info == psb_success_) call a%cp_to(acsr)
+  ! Get the l1-diagonal of D
+  if (do_l1correction) then
+    allocate(l1rwsum(nrow))
+    call acsr%arwsum(l1rwsum)
+    if (info == psb_success_) &
+      & call psb_realloc(ncol,l1rwsum,info)
+    if (info == psb_success_) &
+      & call psb_halo(l1rwsum,desc_a,info)
+    ! \tilde{D}_{i,i} = \sum_{j \ne i} |a_{i,j}|
+    do i=1,size(adiag)
+      adiag(i) = adiag(i) + l1rwsum(i) - abs(adiag(i))
+    end do
+  end if

  if(info /= psb_success_) then
    call psb_errpush(psb_err_from_subroutine_,name,a_err='sp_getdiag')
@@ -246,7 +267,7 @@ subroutine amg_d_parmatch_smth_bld(ag,a,desc_a,ilaggr,nlaggr,parms,&

      enddo
      if (jd == -1) then
-        write(0,*) 'Wrong input: we need the diagonal!!!!', i
+        write(0,*) name,': Warning: there is no diagonal element', i
      else
        acsrf%val(jd)=acsrf%val(jd)-tmp
      end if
@@ -267,7 +288,10 @@ subroutine amg_d_parmatch_smth_bld(ag,a,desc_a,ilaggr,nlaggr,parms,&

  if (parms%aggr_omega_alg == amg_eig_est_) then

-    if (parms%aggr_eig == amg_max_norm_) then
+    if (do_l1correction) then
+      ! For l1-Jacobi this can be estimated with 1
+      parms%aggr_omega_val = done
+    else if (parms%aggr_eig == amg_max_norm_) then
      allocate(arwsum(nrow))
      call acsr%arwsum(arwsum)
      anorm = maxval(abs(adiag(1:nrow)*arwsum(1:nrow)))
@@ -373,6 +397,7 @@ subroutine amg_d_parmatch_smth_bld(ag,a,desc_a,ilaggr,nlaggr,parms,&

    end block
  end if
+  if (allocated(l1rwsum)) deallocate(l1rwsum)
  if (do_timings) call psb_toc(idx_phase2)

  if (debug_level >= psb_debug_outer_) &
--- a/Show More
+++ b/Show More
Author	SHA1	Message	Date
sfilippone	2e79105695	Fix CUDA defines for compilation	2025-02-21 18:07:58 +01:00
sfilippone	921535c3c9	Update copyright	2025-02-21 16:41:39 +01:00
sfilippone	557809755e	Merge branch 'development' of github.com:sfilippone/amg4psblas into development	2025-02-21 16:36:10 +01:00
sfilippone	64105a17c5	Delete changelog	2025-02-21 16:35:45 +01:00
sfilippone	67c32da222	Update docs	2025-02-21 16:33:43 +01:00
Salvatore Filippone	f317a571f4	Update README.md	2025-02-21 16:12:17 +01:00
sfilippone	b075182ce6	Obsolete AINV_SOLVER_SETirc	2025-02-21 16:06:21 +01:00
sfilippone	be469f8844	Added comments on allocat_wrk	2025-02-21 15:39:43 +01:00
Salvatore Filippone	04bcf04a9c	Update README.md	2025-02-21 15:29:59 +01:00
sfilippone	725e39586d	Fix TLU solver_bld and CUDA test program to print format.	2025-02-21 15:22:51 +01:00
sfilippone	e68304e84f	Merge branch 'development' of github.com:sfilippone/amg4psblas into development	2025-02-21 15:06:09 +01:00
sfilippone	3972f27eb5	Rework cuda example	2025-02-21 15:05:35 +01:00
fdurastante	b21e5aebab	Update README.md	2025-02-21 14:33:39 +01:00
Cirdans-Home	6e2d61b51e	Initial README for examples	2025-02-21 13:03:43 +00:00
sfilippone	4ee170f0ec	Update version 1.2.0	2025-02-21 12:32:36 +01:00
sfilippone	4a9f8e9da0	Move newslv to samples and fix it	2025-02-21 12:32:04 +01:00
sfilippone	94caab8aa9	Removed make veryclean, introduced make distclean	2025-02-20 17:00:32 +01:00
sfilippone	749ab2f6ae	Fix extprol. Remove obsolete hybrid aggregator.	2025-02-20 11:53:58 +01:00
sfilippone	bf96dd554c	Merge branch 'repackage' into development	2025-02-19 17:35:47 +01:00
sfilippone	1aef82023c	Merge branch 'repackage' of github.com:sfilippone/amg4psblas into repackage	2025-02-09 16:32:40 +01:00
sfilippone	161b93da64	Fix interface issues of clone_settings	2025-02-09 16:02:00 +01:00
Salvatore Filippone	8d4af1ba9f	Fix output formatting bug	2025-01-19 11:15:44 +01:00
sfilippone	cd87bdb0c1	Add call to clean_zeros()	2025-01-08 17:19:58 +01:00
Cirdans-Home	53afa87814	Added row-preserving filtering	2024-12-13 17:19:04 +00:00
sfilippone	13c99a0c3f	Cleanup use of RICHARDSON	2024-11-18 16:20:35 +01:00
sfilippone	1aa7c8db59	Fix release version of UG	2024-11-18 15:53:52 +01:00
sfilippone	36b57eec24	Define 1.2 as name of UG	2024-11-18 15:51:16 +01:00
fdurastante	417f8beaf9	Added remark on polynomial smoothers	2024-11-18 14:42:31 +01:00
fdurastante	d66bf1e2f8	Added remark on polynomial smoothers	2024-11-18 14:41:59 +01:00
sfilippone	c590a4088a	Merge branch 'repackage' of github.com:sfilippone/amg4psblas into repackage	2024-11-18 12:52:39 +01:00
sfilippone	80dfd1ad3b	Modified README	2024-11-18 12:47:39 +01:00
Pasqua D'Ambra	857844474a	Update README.md	2024-11-18 10:56:46 +01:00
fdurastante	a85997926c	Update README.md	2024-11-18 10:08:29 +01:00
fdurastante	70fb39ae55	Added Polynomial Smoother infos	2024-11-18 10:03:21 +01:00
sfilippone	82de529c54	Makefile in amg4psblas/docs	2024-11-18 09:17:11 +01:00
sfilippone	db9757a45e	Update README Fortran 2003 to 2008	2024-11-15 17:28:14 +01:00
fdurastante	5b654cc221	Update README.md	2024-11-15 17:13:17 +01:00
sfilippone	d518f5eac7	Merge branch 'l1-and-0-aggr' into repackage	2024-11-12 10:07:44 +01:00
sfilippone	053d5c4bc0	Changed sample input file	2024-11-11 18:10:34 +01:00
sfilippone	032543d625	Adjusted CBIND test program.	2024-11-11 13:06:48 +01:00
sfilippone	07149a02ad	Fixes for krylov->linsolve. Test program to be completed.	2024-11-11 09:11:07 +01:00
sfilippone	08a0c744b1	Switch from KRYLOV to LINSOLVE	2024-11-10 14:00:51 +01:00
sfilippone	60324084d8	New Richardson solver.	2024-11-08 18:25:13 +01:00
sfilippone	83ba79d7ae	Fix input file	2024-11-02 10:33:23 +01:00
sfilippone	b704d50df1	Fix dump of TPROL and POLY smoother	2024-10-30 10:49:05 +01:00
Cirdans-Home	68a9cceaa0	Fixed l1-Jacobi scaling	2024-10-29 15:51:37 +00:00
Salvatore Filippone	aec5a52c7f	Missing CBIND clean target	2024-10-16 10:56:53 +02:00
Salvatore Filippone	b5f5d356fd	Missing CBIND clean target	2024-10-16 10:55:59 +02:00
Cirdans-Home	8966ecb4a6	First implementation of l1-aggregation	2024-10-09 09:43:43 +00:00
Cirdans-Home	3e9a5c0c5b	Remove diagonal zeros from aggregation	2024-10-04 12:28:36 +00:00
sfilippone	dfc261cf34	Modify error message in smth_bld	2024-09-12 18:45:20 +02:00
sfilippone	cab98295e2	Improve handling of pointers in hierarchy_bld	2024-08-06 10:28:02 +02:00
sfilippone	5a83c63810	Unused lines in ilu_solver_bld	2024-08-02 15:00:44 +02:00
sfilippone	2e43f55455	Update sample/advanced/pdegeb	2024-08-02 15:00:28 +02:00
sfilippone	244fcda207	Update test program input	2024-08-02 13:48:42 +02:00
sfilippone	ca6fce0765	Fix for test program	2024-08-02 13:46:41 +02:00
sfilippone	b6f92354d3	Reworked sample data generation	2024-08-02 12:39:09 +02:00
sfilippone	1b7fe6a9a7	Fix input data in samples pdegen	2024-07-30 09:49:46 +02:00
sfilippone	14ea4d9c15	Fix sample programs for polynomial smoothers	2024-07-29 17:35:01 +02:00
sfilippone	3ee333baac	Change name abgdxyz into upd_xyz	2024-07-29 16:59:59 +02:00
sfilippone	ecb41dfbbf	Update pde3d.inp	2024-07-16 13:41:51 +02:00
sfilippone	33ac3f786b	Merge latest changes from polysmooth	2024-07-16 13:12:10 +02:00
sfilippone	474c6a3634	Merge branch 'PolySmooth' into development	2024-07-05 16:50:26 +02:00
sfilippone	c1e8bc0c57	Do not use OpenMP code for serial version	2024-07-05 14:23:44 +02:00
sfilippone	2f5072166d	Switch off OpenMP in certain sections of MatchBOXP	2024-07-05 13:56:52 +02:00
sfilippone	89e2d53e8b	Silence debug print.	2024-07-05 13:56:38 +02:00
sfilippone	bfe0a32e09	Merge branch 'mboxomp' into PolySmooth	2024-07-05 13:08:56 +02:00
sfilippone	e88d176fed	Switch off OpenMP in processExposedVertex	2024-07-05 12:10:26 +02:00
sfilippone	c96727a97c	Merge branch 'PolySmooth' into mboxomp	2024-07-03 11:01:08 +02:00
sfilippone	6362db0cc5	Try improve OpenMP version of matchbox	2024-07-03 11:00:13 +02:00
sfilippone	9239b16175	Merge branch 'PolySmooth' of github.com:sfilippone/amg4psblas into PolySmooth	2024-07-03 10:58:33 +02:00
sfilippone	96a700cb9d	Workaround for save_smoothers bug with gcc-13.3.0	2024-07-03 10:58:00 +02:00
sfilippone	41d91120d4	Merge branch 'PolySmooth' into mboxomp	2024-06-11 09:57:53 +02:00
sfilippone	5d20407b15	Merge branch 'PolySmooth' of github.com:sfilippone/amg4psblas into PolySmooth	2024-06-11 09:56:39 +02:00
sfilippone	322e3f65d1	Merge branch 'PolySmooth' of github.com:sfilippone/amg4psblas into PolySmooth	2024-06-11 09:55:07 +02:00
sfilippone	3ff1ad9372	Fix noise in %memory_use()	2024-06-11 09:54:55 +02:00
sfilippone	818ead5878	Try changes for matching	2024-06-11 09:52:11 +02:00
sfilippone	803d311d1c	S versions. Take out parallel in a few places	2024-06-06 17:12:31 +02:00
sfilippone	677e4fe6bc	Modify MatchBox names with D in preparation for S version	2024-06-05 15:13:18 +02:00
sfilippone	02a83575a2	Reorganize MatchBox (prepare for S OpenMP)	2024-06-05 13:11:41 +02:00
sfilippone	cfbec1f6ea	Cosmetic changes to MatchBoxPC	2024-06-05 10:38:30 +02:00
sfilippone	e11a134a1f	Additional tinkering with OpenMP matchbox	2024-06-04 17:39:31 +02:00
sfilippone	6d05120930	Merge branch 'PolySmooth' of github.com:sfilippone/amg4psblas into PolySmooth	2024-06-04 16:56:07 +02:00
sfilippone	bd2d1e3b26	Additional OpenMP tinkering in matchboxp	2024-06-04 16:53:52 +02:00
sfilippone	c9605d1b29	Mods for OpenMP	2024-06-04 13:25:12 +02:00
sfilippone	67594f8b07	Fixes for OpenMP	2024-05-30 17:35:15 +02:00
sfilippone	301fb57bb1	Merge branch 'PolySmooth' of github.com:sfilippone/amg4psblas into PolySmooth	2024-05-30 17:26:49 +02:00
sfilippone	13eee99ea3	Use ifdef OPENMP	2024-05-30 17:26:34 +02:00
sfilippone	fb802c62cd	Merge PSBCXXDEFINES into AMGCXXDEFINES	2024-05-30 17:25:25 +02:00
sfilippone	767b606bb2	Do away with -DOMP	2024-05-30 16:14:32 +02:00
sfilippone	8492c07521	Merge PSBCXXDEFINES into CXXDEFINES	2024-05-30 16:14:12 +02:00
sfilippone	17698c2725	Changes to OpenMP MathcBox version needs -DOMP	2024-05-30 14:43:26 +02:00
Salvatore Filippone	897c5229a6	Improve behaviour of OpenMP matching	2024-05-30 08:13:08 -04:00
Salvatore Filippone	ab5eaac5ed	Cosmetic changes	2024-05-30 08:12:48 -04:00
sfilippone	234071869d	Merge branch 'PolySmooth' of github.com:sfilippone/amg4psblas into PolySmooth	2024-05-10 15:03:38 +02:00
sfilippone	3e3b343131	Fix potential overflow issue in SOC_MAP_BLD	2024-05-10 15:02:36 +02:00
Cirdans-Home	5790aa0cbd	Revert "First hardcoded implementation of l1 smooth aggregation" This reverts commit `a17f503486`.	2024-04-26 10:55:47 +00:00
Cirdans-Home	a17f503486	First hardcoded implementation of l1 smooth aggregation	2024-04-26 10:54:43 +00:00
Cirdans-Home	74dccb6c44	Added timers and removed unuseful spmm	2024-04-22 07:52:14 +00:00
sfilippone	e83bde6896	New timings	2024-04-17 14:45:20 +02:00
Salvatore Filippone	83d435b49e	Default GLOBAL=.true. for MEMORY_USE	2024-03-16 15:22:24 +01:00
Salvatore Filippone	2ef4459b18	Added "COARSE_INVFILL"	2024-02-21 11:51:18 +01:00
sfilippone	c2fd0ac66d	Disable MATCHBOX with SERIAL_MPI and add error message	2023-11-26 11:43:10 +01:00
sfilippone	5387e206b1	Fixed sample programs	2023-11-24 16:07:22 +01:00
sfilippone	fc34385341	Fix matrix generation for samples	2023-11-16 17:55:42 +01:00
sfilippone	5fbdfb1436	Fix free for jac_solver	2023-11-16 17:55:24 +01:00