Fix intent for ovrl

.gitignore

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

Make.inc.in

@@ -75,7 +75,7 @@ CDEFINES=$(AMGCDEFINES)
 AMGFDEFINES=@AMGFDEFINES@ $(PSBFDEFINES)
 FDEFINES=$(AMGFDEFINES)
 
-CXXDEFINES=@AMGCXXDEFINES@ $(PSBCXXDEFINES)
+CXXDEFINES=@AMGCXXDEFINES@
 
 @COMPILERULES@
 

Make_n.inc.in

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

Makefile

@@ -3,9 +3,9 @@ include Make.inc
 
 all:  objs lib
 
-objs: libdir amgp cbnd
+objs: amgp cbnd
 
-lib:  objs
+lib: libdir objs
 	cd amgprec && $(MAKE) lib
 	cd cbind && $(MAKE)  lib
 
@@ -46,7 +46,6 @@ cleanlib:
 
 veryclean: cleanlib
 	(cd amgprec && $(MAKE) veryclean)
-	(cd cbind && $(MAKE) veryclean)
 	(cd samples/simple/fileread && $(MAKE) clean)
 	(cd samples/simple/pdegen && $(MAKE) clean)
 	(cd samples/advanced/fileread && $(MAKE) clean)
@@ -57,4 +56,3 @@ check: all
 
 clean:
 	(cd amgprec && $(MAKE) clean)
-	(cd cbind && $(MAKE) clean)

README.md

@@ -1,76 +1,54 @@
-# AMG4PSBLAS  v1.2
-Algebraic Multigrid Package  based on [PSBLAS](https://github.com/sfilippone/psblas3) (Parallel Sparse BLAS version 3.9)
+                         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 is a package of parallel algebraic multilevel preconditioners included in the PSCToolkit (Parallel Sparse Computation Toolkit) software framework.
+---------------------------------------------------------------------
 
-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.
+AMG4PSBLAS is a package of Algebraic MultiGrid (AMG)
+preconditioners for  the iterative solution of large and sparse linear systems.
 
-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), but it has been
+thoroughly reworked, and it is sufficiently different to warrant a new
+project name.
 
-It is an evolution of MLD2P4 (see [LICENSE.MLD2P4](LICENSE.MLD2P4)), but due to the significant number of changes and the increase in scope, we decided to rename the package as AMG4PSBLAS.
 
-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.
+MAIN REFERENCES:
 
-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.
+     
 
-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.
+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.
 
-## Main Refrerences:
 
-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.
+TO COMPILE
 
 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 [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`; 
+   See MLD2P4 User's and Reference Guide (Section 3) for details.
+2. Tweak Make.inc if you are not satisfied.
+3. 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.
 
-### CUDA, OpeMP, OpenACC
+NOTES
 
-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.
+- 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.
 
-### EoCoE - Software as service portal
-
-In the European project “Energy oriented Center of Excellence: toward exascale for energy” we made available a software as service portal: [https://eocoe.psnc.pl/](https://eocoe.psnc.pl/). This permits to test several cutting-edge computational methods for accelerating the transition to the production, storage and management of clean, decarbonized energy. Among them you have the possibility of running PSBLAS+AMG4PSBLAS on some test problems to become familiar with using the software.
-
-## TODO and bugs
- 
-- [X] Fix all reamining bugs. Bugs? We dont' have any ! 🤓
-
-> [!NOTE]
-> To report bugs 🐛 or issues ❓ please use the [GitHub issue system](https://github.com/sfilippone/amg4psblas/issues).
-
-## The AMG4PSBLAS team. 
-
-- Pasqua D'Ambra         (IAC-CNR, Naples, IT)
-- Fabio Durastante       (University of Pisa and IAC-CNR, IT)
-- Salvatore Filippone    (University of Rome Tor Vergata and IAC-CNR, IT)
 
+The AMG4PSBLAS team. 
+---------------
+Salvatore Filippone
+Pasqua     D'Ambra
+Fabio     Durastante

amgprec/amg_base_prec_type.F90

@@ -288,19 +288,17 @@ 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_l1_smooth_prol_ = 2
-  integer(psb_ipk_), parameter :: amg_min_energy_     = 3
+  integer(psb_ipk_), parameter :: amg_no_smooth_   = 0
+  integer(psb_ipk_), parameter :: amg_smooth_prol_ = 1
+  integer(psb_ipk_), parameter :: amg_min_energy_  = 2
   ! Disabling min_energy  for the time being.
-  integer(psb_ipk_), parameter :: amg_max_aggr_prol_= amg_l1_smooth_prol_
+  integer(psb_ipk_), parameter :: amg_max_aggr_prol_=amg_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_filter_prow_mat_ = 2
-  integer(psb_ipk_), parameter :: amg_max_filter_mat_  = amg_filter_prow_mat_
+  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_
   !
   ! Legal values for entry: amg_aggr_ord_
   !
@@ -325,10 +323,9 @@ module amg_base_prec_type
   !
   ! Legal values for entry: amg_poly_variant_
   !
-  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_lottes_      = 0
+  integer(psb_ipk_), parameter :: amg_poly_lottes_beta_ = 1
+  integer(psb_ipk_), parameter :: amg_poly_new_         = 2
   
   integer(psb_ipk_), parameter :: amg_poly_rho_est_power_ = 0
 
@@ -378,11 +375,10 @@ module amg_base_prec_type
   character(len=19), parameter, private :: &
        &  eigen_estimates(0:0)=(/'infinity norm     '/)
   character(len=15), parameter, private :: &
-       &  aggr_prols(0:4)=(/'unsmoothed    ','smoothed      ',&
-       &  'l1-smoothed   ','min energy    ','bizr. smoothed'/)
+       &  aggr_prols(0:3)=(/'unsmoothed    ','smoothed      ',&
+       &           'min energy    ','bizr. smoothed'/)
   character(len=15), parameter, private :: &
-       &  aggr_filters(0:2)=(/'no filtering  ','filtering     ',&
-       &  'filtering rsum'/)
+       &  aggr_filters(0:1)=(/'no filtering  ','filtering     '/)
   character(len=15), parameter, private :: &
        &  matrix_names(0:1)=(/'distributed   ','replicated    '/)
   character(len=18), parameter, private :: &
@@ -551,8 +547,6 @@ 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')
@@ -575,14 +569,12 @@ contains
       val = amg_as_
     case('POLY')
       val = amg_poly_
-    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_LOTTES')
+      val = amg_poly_lottes_
+    case('POLY_LOTTES_BETA')
+      val = amg_poly_lottes_beta_
+    case('POLY_NEW')
+      val = amg_poly_new_
     case('POLY_RHO_EST_POWER')
       val =  amg_poly_rho_est_power_
     case('A_NORMI')
@@ -593,8 +585,6 @@ contains
       val = amg_eig_est_
     case('FILTER')
       val = amg_filter_mat_
-    case('FILTERROWSUM')
-      val = amg_filter_prow_mat_
     case('NOFILTER','NO_FILTER')
       val = amg_no_filter_mat_
     case('OUTER_SWEEPS')

amgprec/amg_c_ainv_solver.F90

@@ -91,9 +91,9 @@ module amg_c_ainv_solver
       import :: psb_desc_type, psb_cspmat_type,  psb_c_base_sparse_mat, &
        & amg_c_base_solver_type, psb_dpk_, amg_c_ainv_solver_type, psb_ipk_
       Implicit None
-      class(amg_c_ainv_solver_type), intent(inout) :: sv
-      class(amg_c_base_solver_type), intent(inout) :: svout
-      integer(psb_ipk_), intent(out)                 :: info
+      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
     end subroutine amg_c_ainv_solver_clone_settings
   end interface
 
@@ -106,7 +106,7 @@ module amg_c_ainv_solver
       Implicit None
 
       ! Arguments
-      type(psb_cspmat_type), intent(in), target           :: a
+      type(psb_cspmat_type), intent(inout), 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

amgprec/amg_c_as_smoother.f90

@@ -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(in), target        :: a
+      type(psb_cspmat_type), intent(inout), 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

amgprec/amg_c_base_smoother_mod.f90

@@ -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(in), target     :: a
+      type(psb_cspmat_type), intent(inout), 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

amgprec/amg_c_base_solver_mod.f90

@@ -170,7 +170,7 @@ module amg_c_base_solver_mod
       Implicit None
       
       ! Arguments
-      type(psb_cspmat_type), intent(in), target             :: a
+      type(psb_cspmat_type), intent(inout), 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

amgprec/amg_c_diag_solver.f90

@@ -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(in), target           :: a
+      type(psb_cspmat_type), intent(inout), 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(in), target           :: a
+      type(psb_cspmat_type), intent(inout), 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

amgprec/amg_c_gs_solver.f90

@@ -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(in), target           :: a
+      type(psb_cspmat_type), intent(inout), 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(in), target           :: a
+      type(psb_cspmat_type), intent(inout), 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

amgprec/amg_c_id_solver.f90

@@ -123,7 +123,7 @@ contains
     Implicit None
 
     ! Arguments
-    type(psb_cspmat_type), intent(in), target           :: a
+    type(psb_cspmat_type), intent(inout), 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

amgprec/amg_c_ilu_solver.f90

@@ -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(in), target           :: a
+      type(psb_cspmat_type), intent(inout), 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

amgprec/amg_c_inner_mod.f90

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

amgprec/amg_c_invk_solver.f90

@@ -79,9 +79,9 @@ module amg_c_invk_solver
       import :: psb_desc_type, psb_cspmat_type,  psb_c_base_sparse_mat, &
        & amg_c_base_solver_type, psb_spk_, amg_c_invk_solver_type, psb_ipk_
       Implicit None
-      class(amg_c_invk_solver_type), intent(inout) :: sv
-      class(amg_c_base_solver_type), intent(inout) :: svout
-      integer(psb_ipk_), intent(out)                 :: info
+      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
     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(in), target           :: a
+      type(psb_cspmat_type), intent(inout), 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

amgprec/amg_c_invt_solver.f90

@@ -79,9 +79,9 @@ module amg_c_invt_solver
       import :: psb_desc_type, psb_cspmat_type,  psb_c_base_sparse_mat, &
        & amg_c_base_solver_type, psb_spk_, amg_c_invt_solver_type, psb_ipk_
       Implicit None
-      class(amg_c_invt_solver_type), intent(inout) :: sv
-      class(amg_c_base_solver_type), intent(inout) :: svout
-      integer(psb_ipk_), intent(out)                 :: info
+      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
     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(in), target           :: a
+      type(psb_cspmat_type), intent(inout), 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

amgprec/amg_c_jac_smoother.f90

@@ -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(in), target           :: a
+      type(psb_cspmat_type), intent(inout), 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), intent(inout) :: smout
+      class(amg_c_jac_smoother_type), intent(inout)               :: sm
+      class(amg_c_base_smoother_type), allocatable, 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(in), target           :: a
+      type(psb_cspmat_type), intent(inout), 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

amgprec/amg_c_jac_solver.f90

@@ -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(in), target           :: a
+      type(psb_cspmat_type), intent(inout), 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(in), target           :: a
+      type(psb_cspmat_type), intent(inout), 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

amgprec/amg_c_krm_solver.f90

@@ -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(in), target           :: a
+      type(psb_cspmat_type), intent(inout), 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

amgprec/amg_c_mumps_solver.F90

@@ -163,7 +163,7 @@ module amg_c_mumps_solver
       Implicit None
 
       ! Arguments
-      type(psb_cspmat_type), intent(in), target           :: a
+      type(psb_cspmat_type), intent(inout), 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

amgprec/amg_c_slu_solver.F90

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

amgprec/amg_d_ainv_solver.F90

@@ -91,9 +91,9 @@ module amg_d_ainv_solver
       import :: psb_desc_type, psb_dspmat_type,  psb_d_base_sparse_mat, &
        & amg_d_base_solver_type, psb_dpk_, amg_d_ainv_solver_type, psb_ipk_
       Implicit None
-      class(amg_d_ainv_solver_type), intent(inout) :: sv
-      class(amg_d_base_solver_type), intent(inout) :: svout
-      integer(psb_ipk_), intent(out)                 :: info
+      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
     end subroutine amg_d_ainv_solver_clone_settings
   end interface
 
@@ -106,7 +106,7 @@ module amg_d_ainv_solver
       Implicit None
 
       ! Arguments
-      type(psb_dspmat_type), intent(in), target           :: a
+      type(psb_dspmat_type), intent(inout), 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

amgprec/amg_d_as_smoother.f90

@@ -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(in), target        :: a
+      type(psb_dspmat_type), intent(inout), 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

amgprec/amg_d_base_smoother_mod.f90

@@ -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(in), target     :: a
+      type(psb_dspmat_type), intent(inout), 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

amgprec/amg_d_base_solver_mod.f90

@@ -170,7 +170,7 @@ module amg_d_base_solver_mod
       Implicit None
       
       ! Arguments
-      type(psb_dspmat_type), intent(in), target             :: a
+      type(psb_dspmat_type), intent(inout), 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

amgprec/amg_d_diag_solver.f90

@@ -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(in), target           :: a
+      type(psb_dspmat_type), intent(inout), 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(in), target           :: a
+      type(psb_dspmat_type), intent(inout), 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

amgprec/amg_d_gs_solver.f90

@@ -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(in), target           :: a
+      type(psb_dspmat_type), intent(inout), 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(in), target           :: a
+      type(psb_dspmat_type), intent(inout), 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

amgprec/amg_d_id_solver.f90

@@ -123,7 +123,7 @@ contains
     Implicit None
 
     ! Arguments
-    type(psb_dspmat_type), intent(in), target           :: a
+    type(psb_dspmat_type), intent(inout), 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

amgprec/amg_d_ilu_solver.f90

@@ -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(in), target           :: a
+      type(psb_dspmat_type), intent(inout), 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

amgprec/amg_d_inner_mod.f90

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

amgprec/amg_d_invk_solver.f90

@@ -79,9 +79,9 @@ module amg_d_invk_solver
       import :: psb_desc_type, psb_dspmat_type,  psb_d_base_sparse_mat, &
        & amg_d_base_solver_type, psb_dpk_, amg_d_invk_solver_type, psb_ipk_
       Implicit None
-      class(amg_d_invk_solver_type), intent(inout) :: sv
-      class(amg_d_base_solver_type), intent(inout) :: svout
-      integer(psb_ipk_), intent(out)                 :: info
+      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
     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(in), target           :: a
+      type(psb_dspmat_type), intent(inout), 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

amgprec/amg_d_invt_solver.f90

@@ -79,9 +79,9 @@ module amg_d_invt_solver
       import :: psb_desc_type, psb_dspmat_type,  psb_d_base_sparse_mat, &
        & amg_d_base_solver_type, psb_dpk_, amg_d_invt_solver_type, psb_ipk_
       Implicit None
-      class(amg_d_invt_solver_type), intent(inout) :: sv
-      class(amg_d_base_solver_type), intent(inout) :: svout
-      integer(psb_ipk_), intent(out)                 :: info
+      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
     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(in), target           :: a
+      type(psb_dspmat_type), intent(inout), 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

amgprec/amg_d_jac_smoother.f90

@@ -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(in), target           :: a
+      type(psb_dspmat_type), intent(inout), 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), intent(inout) :: smout
+      class(amg_d_jac_smoother_type), intent(inout)               :: sm
+      class(amg_d_base_smoother_type), allocatable, 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(in), target           :: a
+      type(psb_dspmat_type), intent(inout), 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

amgprec/amg_d_jac_solver.f90

@@ -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(in), target           :: a
+      type(psb_dspmat_type), intent(inout), 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(in), target           :: a
+      type(psb_dspmat_type), intent(inout), 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

amgprec/amg_d_krm_solver.f90

@@ -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(in), target           :: a
+      type(psb_dspmat_type), intent(inout), 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

amgprec/amg_d_mumps_solver.F90

@@ -163,7 +163,7 @@ module amg_d_mumps_solver
       Implicit None
 
       ! Arguments
-      type(psb_dspmat_type), intent(in), target           :: a
+      type(psb_dspmat_type), intent(inout), 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

amgprec/amg_d_parmatch_aggregator_mod.F90

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

amgprec/amg_d_poly_smoother.f90

@@ -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(in), target           :: a
+      type(psb_dspmat_type), intent(inout), 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), intent(inout) :: smout
+      class(amg_d_poly_smoother_type), intent(inout)               :: sm
+      class(amg_d_base_smoother_type), allocatable, 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_cheb_4_
+    sm%variant      = amg_poly_lottes_
     sm%rho_estimate = amg_poly_rho_est_power_
     sm%rho_estimate_iterations = 20
     if (allocated(sm%sv)) then

amgprec/amg_d_slu_solver.F90

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

amgprec/amg_d_sludist_solver.F90

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

amgprec/amg_d_umf_solver.F90

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

amgprec/amg_s_ainv_solver.F90

@@ -91,9 +91,9 @@ module amg_s_ainv_solver
       import :: psb_desc_type, psb_sspmat_type,  psb_s_base_sparse_mat, &
        & amg_s_base_solver_type, psb_dpk_, amg_s_ainv_solver_type, psb_ipk_
       Implicit None
-      class(amg_s_ainv_solver_type), intent(inout) :: sv
-      class(amg_s_base_solver_type), intent(inout) :: svout
-      integer(psb_ipk_), intent(out)                 :: info
+      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
     end subroutine amg_s_ainv_solver_clone_settings
   end interface
 
@@ -106,7 +106,7 @@ module amg_s_ainv_solver
       Implicit None
 
       ! Arguments
-      type(psb_sspmat_type), intent(in), target           :: a
+      type(psb_sspmat_type), intent(inout), 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

amgprec/amg_s_as_smoother.f90

@@ -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(in), target        :: a
+      type(psb_sspmat_type), intent(inout), 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

amgprec/amg_s_base_smoother_mod.f90

@@ -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(in), target     :: a
+      type(psb_sspmat_type), intent(inout), 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

amgprec/amg_s_base_solver_mod.f90

@@ -170,7 +170,7 @@ module amg_s_base_solver_mod
       Implicit None
       
       ! Arguments
-      type(psb_sspmat_type), intent(in), target             :: a
+      type(psb_sspmat_type), intent(inout), 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

amgprec/amg_s_diag_solver.f90

@@ -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(in), target           :: a
+      type(psb_sspmat_type), intent(inout), 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(in), target           :: a
+      type(psb_sspmat_type), intent(inout), 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

amgprec/amg_s_gs_solver.f90

@@ -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(in), target           :: a
+      type(psb_sspmat_type), intent(inout), 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(in), target           :: a
+      type(psb_sspmat_type), intent(inout), 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

amgprec/amg_s_id_solver.f90

@@ -123,7 +123,7 @@ contains
     Implicit None
 
     ! Arguments
-    type(psb_sspmat_type), intent(in), target           :: a
+    type(psb_sspmat_type), intent(inout), 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

amgprec/amg_s_ilu_solver.f90

@@ -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(in), target           :: a
+      type(psb_sspmat_type), intent(inout), 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

amgprec/amg_s_inner_mod.f90

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

amgprec/amg_s_invk_solver.f90

@@ -79,9 +79,9 @@ module amg_s_invk_solver
       import :: psb_desc_type, psb_sspmat_type,  psb_s_base_sparse_mat, &
        & amg_s_base_solver_type, psb_spk_, amg_s_invk_solver_type, psb_ipk_
       Implicit None
-      class(amg_s_invk_solver_type), intent(inout) :: sv
-      class(amg_s_base_solver_type), intent(inout) :: svout
-      integer(psb_ipk_), intent(out)                 :: info
+      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
     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(in), target           :: a
+      type(psb_sspmat_type), intent(inout), 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

amgprec/amg_s_invt_solver.f90

@@ -79,9 +79,9 @@ module amg_s_invt_solver
       import :: psb_desc_type, psb_sspmat_type,  psb_s_base_sparse_mat, &
        & amg_s_base_solver_type, psb_spk_, amg_s_invt_solver_type, psb_ipk_
       Implicit None
-      class(amg_s_invt_solver_type), intent(inout) :: sv
-      class(amg_s_base_solver_type), intent(inout) :: svout
-      integer(psb_ipk_), intent(out)                 :: info
+      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
     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(in), target           :: a
+      type(psb_sspmat_type), intent(inout), 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

amgprec/amg_s_jac_smoother.f90

@@ -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(in), target           :: a
+      type(psb_sspmat_type), intent(inout), 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), intent(inout) :: smout
+      class(amg_s_jac_smoother_type), intent(inout)               :: sm
+      class(amg_s_base_smoother_type), allocatable, 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(in), target           :: a
+      type(psb_sspmat_type), intent(inout), 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

amgprec/amg_s_jac_solver.f90

@@ -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(in), target           :: a
+      type(psb_sspmat_type), intent(inout), 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(in), target           :: a
+      type(psb_sspmat_type), intent(inout), 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

amgprec/amg_s_krm_solver.f90

@@ -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(in), target           :: a
+      type(psb_sspmat_type), intent(inout), 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

amgprec/amg_s_mumps_solver.F90

@@ -163,7 +163,7 @@ module amg_s_mumps_solver
       Implicit None
 
       ! Arguments
-      type(psb_sspmat_type), intent(in), target           :: a
+      type(psb_sspmat_type), intent(inout), 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

amgprec/amg_s_parmatch_aggregator_mod.F90

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

amgprec/amg_s_poly_smoother.f90

@@ -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(in), target           :: a
+      type(psb_sspmat_type), intent(inout), 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), intent(inout) :: smout
+      class(amg_s_poly_smoother_type), intent(inout)               :: sm
+      class(amg_s_base_smoother_type), allocatable, 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_cheb_4_
+    sm%variant      = amg_poly_lottes_
     sm%rho_estimate = amg_poly_rho_est_power_
     sm%rho_estimate_iterations = 20
     if (allocated(sm%sv)) then

amgprec/amg_s_slu_solver.F90

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

amgprec/amg_z_ainv_solver.F90

@@ -91,9 +91,9 @@ module amg_z_ainv_solver
       import :: psb_desc_type, psb_zspmat_type,  psb_z_base_sparse_mat, &
        & amg_z_base_solver_type, psb_dpk_, amg_z_ainv_solver_type, psb_ipk_
       Implicit None
-      class(amg_z_ainv_solver_type), intent(inout) :: sv
-      class(amg_z_base_solver_type), intent(inout) :: svout
-      integer(psb_ipk_), intent(out)                 :: info
+      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
     end subroutine amg_z_ainv_solver_clone_settings
   end interface
 
@@ -106,7 +106,7 @@ module amg_z_ainv_solver
       Implicit None
 
       ! Arguments
-      type(psb_zspmat_type), intent(in), target           :: a
+      type(psb_zspmat_type), intent(inout), 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

amgprec/amg_z_as_smoother.f90

@@ -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(in), target        :: a
+      type(psb_zspmat_type), intent(inout), 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

amgprec/amg_z_base_smoother_mod.f90

@@ -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(in), target     :: a
+      type(psb_zspmat_type), intent(inout), 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

amgprec/amg_z_base_solver_mod.f90

@@ -170,7 +170,7 @@ module amg_z_base_solver_mod
       Implicit None
       
       ! Arguments
-      type(psb_zspmat_type), intent(in), target             :: a
+      type(psb_zspmat_type), intent(inout), 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

amgprec/amg_z_diag_solver.f90

@@ -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(in), target           :: a
+      type(psb_zspmat_type), intent(inout), 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(in), target           :: a
+      type(psb_zspmat_type), intent(inout), 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

amgprec/amg_z_gs_solver.f90

@@ -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(in), target           :: a
+      type(psb_zspmat_type), intent(inout), 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(in), target           :: a
+      type(psb_zspmat_type), intent(inout), 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

amgprec/amg_z_id_solver.f90

@@ -123,7 +123,7 @@ contains
     Implicit None
 
     ! Arguments
-    type(psb_zspmat_type), intent(in), target           :: a
+    type(psb_zspmat_type), intent(inout), 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

amgprec/amg_z_ilu_solver.f90

@@ -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(in), target           :: a
+      type(psb_zspmat_type), intent(inout), 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

amgprec/amg_z_inner_mod.f90

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

amgprec/amg_z_invk_solver.f90

@@ -79,9 +79,9 @@ module amg_z_invk_solver
       import :: psb_desc_type, psb_zspmat_type,  psb_z_base_sparse_mat, &
        & amg_z_base_solver_type, psb_dpk_, amg_z_invk_solver_type, psb_ipk_
       Implicit None
-      class(amg_z_invk_solver_type), intent(inout) :: sv
-      class(amg_z_base_solver_type), intent(inout) :: svout
-      integer(psb_ipk_), intent(out)                 :: info
+      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
     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(in), target           :: a
+      type(psb_zspmat_type), intent(inout), 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

amgprec/amg_z_invt_solver.f90

@@ -79,9 +79,9 @@ module amg_z_invt_solver
       import :: psb_desc_type, psb_zspmat_type,  psb_z_base_sparse_mat, &
        & amg_z_base_solver_type, psb_dpk_, amg_z_invt_solver_type, psb_ipk_
       Implicit None
-      class(amg_z_invt_solver_type), intent(inout) :: sv
-      class(amg_z_base_solver_type), intent(inout) :: svout
-      integer(psb_ipk_), intent(out)                 :: info
+      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
     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(in), target           :: a
+      type(psb_zspmat_type), intent(inout), 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

amgprec/amg_z_jac_smoother.f90

@@ -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(in), target           :: a
+      type(psb_zspmat_type), intent(inout), 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), intent(inout) :: smout
+      class(amg_z_jac_smoother_type), intent(inout)               :: sm
+      class(amg_z_base_smoother_type), allocatable, 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(in), target           :: a
+      type(psb_zspmat_type), intent(inout), 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

amgprec/amg_z_jac_solver.f90

@@ -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(in), target           :: a
+      type(psb_zspmat_type), intent(inout), 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(in), target           :: a
+      type(psb_zspmat_type), intent(inout), 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

amgprec/amg_z_krm_solver.f90

@@ -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(in), target           :: a
+      type(psb_zspmat_type), intent(inout), 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

amgprec/amg_z_mumps_solver.F90

@@ -163,7 +163,7 @@ module amg_z_mumps_solver
       Implicit None
 
       ! Arguments
-      type(psb_zspmat_type), intent(in), target           :: a
+      type(psb_zspmat_type), intent(inout), 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

amgprec/amg_z_slu_solver.F90

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

amgprec/amg_z_sludist_solver.F90

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

amgprec/amg_z_umf_solver.F90

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

amgprec/impl/aggregator/MatchBoxPC.cpp

@@ -67,13 +67,14 @@ void dMatchBoxPC(MilanLongInt NLVer, MilanLongInt NLEdge,
 #endif
 
 
-#undef TIME_TRACKER
-#ifdef TIME_TRACKER
-  double tmr = MPI_Wtime();
-#endif
+#define TIME_TRACKER
+    #ifdef TIME_TRACKER
+        double tmr = MPI_Wtime();
+    #endif
 
-#if defined(OPENMP)
-	//fprintf(stderr,"Warning: using buggy OpenMP matching!\n");
+// Rimosso per tornare al vecchio matching #define OMP
+#ifdef OMP
+				fprintf(stderr,"Warning: using buggy OpenMP matching!\n");
         dalgoDistEdgeApproxDomEdgesLinearSearchMesgBndlSmallMateCMP(NLVer, NLEdge,
 							   verLocPtr, verLocInd, edgeLocWeight,
 							   verDistance,  Mate,
@@ -92,11 +93,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
 }
 
@@ -114,24 +115,13 @@ 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 );
-#endif
+							   verLocPtr, verLocInd, edgeLocWeight,
+							   verDistance,  Mate,
+							   myRank, numProcs, C_comm,
+							   msgIndSent, msgActualSent, msgPercent,
+							   ph0_time, ph1_time, ph2_time,
+							   ph1_card, ph2_card );
 #endif
 }
 

amgprec/impl/aggregator/MatchBoxPC.h

@@ -59,7 +59,7 @@
 #include <assert.h>
 #include <map>
 #include <vector>
-#ifdef OPENMP
+#ifdef OMP
 // OpenMP is included and used if and only if the OpenMP version of the matching
 // is required
 #include "omp.h"
@@ -82,8 +82,6 @@ const int BundleTag = 9;  // Predefined tag
 
 static vector<MilanLongInt> DEFAULT_VECTOR;
 
-#if !defined(SERIAL_MPI)
-
 // MPI type map
 template <typename T>
 MPI_Datatype TypeMap();
@@ -95,7 +93,6 @@ template <>
 inline MPI_Datatype TypeMap<double>() { return MPI_DOUBLE; }
 template <>
 inline MPI_Datatype TypeMap<float>() { return MPI_FLOAT; }
-#endif
 
 #ifdef __cplusplus
 extern "C"
@@ -181,416 +178,264 @@ 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);
-  
-  MilanLongInt firstComputeCandidateMateD(MilanLongInt adj1,
-					 MilanLongInt adj2,
-					 MilanLongInt *verLocInd,
-					 MilanReal *edgeLocWeight);
+    // Function of find the owner of a ghost vertex using binary search:
+    MilanInt findOwnerOfGhost(MilanLongInt vtxIndex, MilanLongInt *mVerDistance,
+                              MilanInt myRank, MilanInt numProcs);
 
-  
-  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);
+    MilanLongInt firstComputeCandidateMate(MilanLongInt adj1,
+                                           MilanLongInt adj2,
+                                           MilanLongInt *verLocInd,
+                                           MilanReal *edgeLocWeight);
 
-  MilanLongInt firstComputeCandidateMateS(MilanLongInt adj1,
-					 MilanLongInt adj2,
-					 MilanLongInt *verLocInd,
-					  MilanFloat *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);
 
-  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);
+    bool isAlreadyMatched(MilanLongInt node,
+                          MilanLongInt StartIndex,
+                          MilanLongInt EndIndex,
+                          vector<MilanLongInt> &GMate,
+                          MilanLongInt *Mate,
+                          map<MilanLongInt, MilanLongInt> &Ghost2LocalMap);
 
-  
-  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);
+    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
 
 
+#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,

amgprec/impl/aggregator/algoDistEdgeApproxDomEdgesLinearSearchMesgBndlSmallMateC.cpp

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

amgprec/impl/aggregator/algoDistEdgeApproxDomEdgesLinearSearchMesgBndlSmallMateCMP.cpp

@@ -1,4 +1,5 @@
 #include "MatchBoxPC.h"
+#ifdef OMP
 // ***********************************************************************
 //
 //        MatchboxP: A C++ library for approximate weighted matching
@@ -125,10 +126,8 @@ void dalgoDistEdgeApproxDomEdgesLinearSearchMesgBndlSmallMateCMP(
     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;
+    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
 
     MPI_Status computeStatus;
 
@@ -146,8 +145,7 @@ void dalgoDistEdgeApproxDomEdgesLinearSearchMesgBndlSmallMateCMP(
     // 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; 
+    vector<MilanInt> QOwner; // Changed by Fabio to be an integer, addresses needs to be integers!
 
     MilanLongInt *PCounter = new MilanLongInt[numProcs];
     for (int i = 0; i < numProcs; i++)
@@ -155,8 +153,7 @@ void dalgoDistEdgeApproxDomEdgesLinearSearchMesgBndlSmallMateCMP(
 
     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; 
+    MilanInt ghostOwner = 0; // Changed by Fabio to be an integer, addresses needs to be integers!
     MilanLongInt *candidateMate = nullptr;
 #ifdef PRINT_DEBUG_INFO_
     cout << "\n(" << myRank << ")NV: " << NLVer << "  Edges: " << NLEdge;
@@ -171,12 +168,9 @@ void dalgoDistEdgeApproxDomEdgesLinearSearchMesgBndlSmallMateCMP(
     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; 
+    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
 
 #ifdef PRINT_DEBUG_INFO_
     cout << "\n(" << myRank << ")About to compute Ghost Vertices...";
@@ -243,7 +237,7 @@ void dalgoDistEdgeApproxDomEdgesLinearSearchMesgBndlSmallMateCMP(
      * PARALLEL_COMPUTE_CANDIDATE_MATE_B is now totally parallel.
      */
 
-    PARALLEL_COMPUTE_CANDIDATE_MATE_BD(NLVer,
+    PARALLEL_COMPUTE_CANDIDATE_MATE_B(NLVer,
                                       verLocPtr,
                                       verLocInd,
                                       myRank,
@@ -268,7 +262,7 @@ void dalgoDistEdgeApproxDomEdgesLinearSearchMesgBndlSmallMateCMP(
      * TODO: Test when it's actually more efficient to execute this code
      *       in parallel.
      */
-    PARALLEL_PROCESS_EXPOSED_VERTEX_BD(NLVer,
+    PARALLEL_PROCESS_EXPOSED_VERTEX_B(NLVer,
                                       candidateMate,
                                       verLocInd,
                                       verLocPtr,
@@ -320,7 +314,7 @@ void dalgoDistEdgeApproxDomEdgesLinearSearchMesgBndlSmallMateCMP(
     vector<MilanLongInt> UChunkBeingProcessed;
     UChunkBeingProcessed.reserve(UCHUNK);
 
-    processMatchedVerticesD(NLVer,
+    processMatchedVertices(NLVer,
                            UChunkBeingProcessed,
                            U,
                            privateU,
@@ -429,7 +423,7 @@ void dalgoDistEdgeApproxDomEdgesLinearSearchMesgBndlSmallMateCMP(
         /////////////////////////// PROCESS MATCHED VERTICES //////////////////////////////
         ///////////////////////////////////////////////////////////////////////////////////
 
-      processMatchedVerticesAndSendMessagesD(NLVer,
+      processMatchedVerticesAndSendMessages(NLVer,
 					    UChunkBeingProcessed,
 					    U,
 					    privateU,
@@ -489,8 +483,8 @@ void dalgoDistEdgeApproxDomEdgesLinearSearchMesgBndlSmallMateCMP(
         ///////////////////////////////////////////////////////////////////////////////////
         /////////////////////////// PROCESS MESSAGES //////////////////////////////////////
         ///////////////////////////////////////////////////////////////////////////////////
-	//startTime = MPI_Wtime();
-        processMessagesD(NLVer,
+
+        processMessages(NLVer,
                         Mate,
                         candidateMate,
                         Ghost2LocalMap,
@@ -555,488 +549,6 @@ 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

amgprec/impl/aggregator/amg_c_dec_aggregator_mat_bld.f90

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

amgprec/impl/aggregator/amg_c_ptap_bld.f90

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

amgprec/impl/aggregator/amg_c_rap.f90

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

amgprec/impl/aggregator/amg_c_soc1_map_bld.F90

@@ -250,7 +250,7 @@ subroutine amg_c_soc1_map_bld(iorder,theta,clean_zeros,a,desc_a,nlaggr,ilaggr,in
               ! we will not reset. 
               if (j>nr) cycle step1
               if (ilaggr(j) > 0) cycle step1
-              if ((abs(val(k)) > theta*sqrt(abs(diag(i)*diag(j)))).and.(diag(i).ne.czero)) then
+              if (abs(val(k)) > theta*sqrt(abs(diag(i)*diag(j)))) then
                 ip = ip + 1
                 icol(ip) = icol(k)
               end if
@@ -275,8 +275,7 @@ 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)) !be careful about overflow
-              itmp = itmp*nths+kk
+              itmp = (bnds(kk)-1+locnaggr(kk))*nths+kk
               if (itmp < (bnds(kk)-1+locnaggr(kk))) then
                 !$omp atomic update
                 info = max(12345678,info)
@@ -357,7 +356,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)))).and.(diag(i).ne.czero)) then
+          if (abs(val(k)) > theta*sqrt(abs(diag(i)*diag(j)))) then
             ip = ip + 1
             icol(ip) = icol(k)
           end if
@@ -545,3 +544,4 @@ subroutine amg_c_soc1_map_bld(iorder,theta,clean_zeros,a,desc_a,nlaggr,ilaggr,in
   return
 
 end subroutine amg_c_soc1_map_bld
+

amgprec/impl/aggregator/amg_c_soc2_map_bld.F90

@@ -309,8 +309,7 @@ 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)) !be careful about overflow
-              itmp = itmp*nths+kk
+              itmp = (bnds(kk)-1+locnaggr(kk))*nths+kk
               if (itmp < (bnds(kk)-1+locnaggr(kk))) then
                 !$omp atomic update
                 info = max(12345678,info)

amgprec/impl/aggregator/amg_caggrmat_minnrg_bld.f90

@@ -69,7 +69,6 @@
 !
 !
 ! 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.
@@ -105,8 +104,8 @@
 !                  Error code.
 !
 !
-subroutine amg_caggrmat_minnrg_bld(dol1smoothing,a,desc_a,ilaggr,nlaggr,&
-                parms,ac,desc_ac,op_prol,op_restr,t_prol,info)
+subroutine amg_caggrmat_minnrg_bld(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
@@ -114,7 +113,6 @@ subroutine amg_caggrmat_minnrg_bld(dol1smoothing,a,desc_a,ilaggr,nlaggr,&
   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(:)
@@ -173,13 +171,6 @@ subroutine amg_caggrmat_minnrg_bld(dol1smoothing,a,desc_a,ilaggr,nlaggr,&
 
   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
@@ -309,7 +300,7 @@ subroutine amg_caggrmat_minnrg_bld(dol1smoothing,a,desc_a,ilaggr,nlaggr,&
 !!$        endif
 !!$      enddo
 !!$      if (jd == -1) then 
-!!$                write(0,*) name,': Warning: there is no diagonal element', i
+!!$        write(0,*) 'Wrong input: we need the diagonal!!!!', i
 !!$      else
 !!$        acsrf%val(jd)=acsrf%val(jd)-tmp
 !!$      end if

amgprec/impl/aggregator/amg_caggrmat_nosmth_bld.f90

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

amgprec/impl/aggregator/amg_caggrmat_smth_bld.f90

@@ -69,8 +69,6 @@
 !
 !
 ! 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.
@@ -104,18 +102,16 @@
 !    info       -  integer, output.
 !                  Error code.
 !
-subroutine amg_caggrmat_smth_bld(dol1smoothing,a,desc_a,ilaggr,nlaggr,&
-                parms,ac,desc_ac,op_prol,op_restr,t_prol,info)
+subroutine amg_caggrmat_smth_bld(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(:)
@@ -136,7 +132,7 @@ subroutine amg_caggrmat_smth_bld(dol1smoothing,a,desc_a,ilaggr,nlaggr,&
   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(:),l1rwsum(:)
+  real(psb_spk_), allocatable :: arwsum(:)
   integer(psb_ipk_)  :: ierr(5)
   logical            :: filter_mat
   integer(psb_ipk_)            :: debug_level, debug_unit, err_act
@@ -145,7 +141,6 @@ subroutine amg_caggrmat_smth_bld(dol1smoothing,a,desc_a,ilaggr,nlaggr,&
   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
 
@@ -178,9 +173,6 @@ subroutine amg_caggrmat_smth_bld(dol1smoothing,a,desc_a,ilaggr,nlaggr,&
   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()
@@ -193,7 +185,7 @@ subroutine amg_caggrmat_smth_bld(dol1smoothing,a,desc_a,ilaggr,nlaggr,&
 
   naggrm1 = sum(nlaggr(1:me))
   naggrp1 = sum(nlaggr(1:me+1))
-  filter_mat = (parms%aggr_filter == amg_filter_mat_).or.(parms%aggr_filter == amg_filter_prow_mat_)
+  filter_mat = (parms%aggr_filter == amg_filter_mat_)
 
   !
   ! naggr: number of local aggregates
@@ -208,24 +200,6 @@ subroutine amg_caggrmat_smth_bld(dol1smoothing,a,desc_a,ilaggr,nlaggr,&
   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')
@@ -256,15 +230,9 @@ subroutine amg_caggrmat_smth_bld(dol1smoothing,a,desc_a,ilaggr,nlaggr,&
 
       enddo
       if (jd == -1) then 
-        ! if (.not.do_l1correction) 
         write(0,*) 'Wrong input: we need the diagonal!!!!', i
-      else if (parms%aggr_filter == amg_filter_mat_) then
-        ! We perform filtering in the standard way assuming that A is an M-matrix
+      else
         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 
@@ -272,6 +240,7 @@ subroutine amg_caggrmat_smth_bld(dol1smoothing,a,desc_a,ilaggr,nlaggr,&
     call acsrf%clean_zeros(info)
   end if
 
+
   !$OMP parallel do private(i) schedule(static)
   do i=1,size(adiag)
     if (adiag(i) /= czero) then
@@ -283,17 +252,14 @@ subroutine amg_caggrmat_smth_bld(dol1smoothing,a,desc_a,ilaggr,nlaggr,&
   !$OMP end parallel do 
   if (parms%aggr_omega_alg == amg_eig_est_) 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 
+    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_')
@@ -356,7 +322,6 @@ subroutine amg_caggrmat_smth_bld(dol1smoothing,a,desc_a,ilaggr,nlaggr,&
   if (debug_level >= psb_debug_outer_) &
        & write(debug_unit,*) me,' ',trim(name),&
        & 'Done smooth_aggregate '
-  if (allocated(l1rwsum)) deallocate(l1rwsum)
   call psb_erractionrestore(err_act)
   return
 

amgprec/impl/aggregator/amg_d_dec_aggregator_mat_bld.f90

@@ -177,24 +177,23 @@ 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(parms%aggr_prol,a,desc_a,ilaggr,&
-            nlaggr,parms,ac,desc_ac,op_prol,op_restr,t_prol,info)
+    call amg_daggrmat_nosmth_bld(a,desc_a,ilaggr,nlaggr,&
+         & parms,ac,desc_ac,op_prol,op_restr,t_prol,info)
 
-  case(amg_smooth_prol_,amg_l1_smooth_prol_) 
+  case(amg_smooth_prol_) 
 
-    call amg_daggrmat_smth_bld(parms%aggr_prol,a,desc_a,&
-            ilaggr,nlaggr,parms,ac,desc_ac,op_prol,&
-            op_restr,t_prol,info)
+    call amg_daggrmat_smth_bld(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(parms%aggr_prol,a,desc_a,ilaggr,&
-            nlaggr,parms,ac,desc_ac,op_prol,op_restr,t_prol,info)
+    call amg_daggrmat_minnrg_bld(a,desc_a,ilaggr,nlaggr, &
+         & parms,ac,desc_ac,op_prol,op_restr,t_prol,info)
 
   case default
     info = psb_err_internal_error_

amgprec/impl/aggregator/amg_d_parmatch_aggregator_mat_bld.F90

@@ -184,23 +184,20 @@ 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(parms%aggr_prol,ag,a,desc_a,&
-            ilaggr,nlaggr,parms,ac,desc_ac,op_prol,op_restr,&
-            t_prol,info)
+    call amg_d_parmatch_unsmth_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_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_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(parms%aggr_prol,a,desc_a,&
-            ilaggr,nlaggr,parms,ac,desc_ac,op_prol,op_restr,&
-            t_prol,info)
+    call amg_daggrmat_minnrg_bld(a,desc_a,ilaggr,nlaggr, &
+         & parms,ac,desc_ac,op_prol,op_restr,t_prol,info)
 
   case default
     info = psb_err_internal_error_

amgprec/impl/aggregator/amg_d_parmatch_smth_bld.F90

@@ -69,8 +69,6 @@
 !
 !
 ! 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.
@@ -104,8 +102,8 @@
 !    info       -  integer, output.
 !                  Error code.
 !
-subroutine amg_d_parmatch_smth_bld(dol1smoothing,ag,a,desc_a,ilaggr,nlaggr,&
-                parms,ac,desc_ac,op_prol,op_restr,t_prol,info)
+subroutine amg_d_parmatch_smth_bld(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
@@ -118,7 +116,6 @@ subroutine amg_d_parmatch_smth_bld(dol1smoothing,ag,a,desc_a,ilaggr,nlaggr,&
   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
@@ -140,7 +137,7 @@ subroutine amg_d_parmatch_smth_bld(dol1smoothing,ag,a,desc_a,ilaggr,nlaggr,&
   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(:),l1rwsum(:)
+  real(psb_dpk_), allocatable :: arwsum(:)
   logical            :: filter_mat
   integer(psb_ipk_)            :: debug_level, debug_unit, err_act
   integer(psb_ipk_), parameter :: ncmax=16
@@ -148,7 +145,6 @@ subroutine amg_d_parmatch_smth_bld(dol1smoothing,ag,a,desc_a,ilaggr,nlaggr,&
   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
 
@@ -170,10 +166,6 @@ subroutine amg_d_parmatch_smth_bld(dol1smoothing,ag,a,desc_a,ilaggr,nlaggr,&
   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")
@@ -225,19 +217,6 @@ subroutine amg_d_parmatch_smth_bld(dol1smoothing,ag,a,desc_a,ilaggr,nlaggr,&
   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')
@@ -267,7 +246,7 @@ subroutine amg_d_parmatch_smth_bld(dol1smoothing,ag,a,desc_a,ilaggr,nlaggr,&
 
       enddo
       if (jd == -1) then
-        write(0,*) name,': Warning: there is no diagonal element', i
+        write(0,*) 'Wrong input: we need the diagonal!!!!', i
       else
         acsrf%val(jd)=acsrf%val(jd)-tmp
       end if
@@ -288,10 +267,7 @@ subroutine amg_d_parmatch_smth_bld(dol1smoothing,ag,a,desc_a,ilaggr,nlaggr,&
 
   if (parms%aggr_omega_alg == amg_eig_est_) then
 
-    if (do_l1correction) then
-      ! For l1-Jacobi this can be estimated with 1
-      parms%aggr_omega_val = done
-    else if (parms%aggr_eig == amg_max_norm_) then
+    if (parms%aggr_eig == amg_max_norm_) then
       allocate(arwsum(nrow))
       call acsr%arwsum(arwsum)
       anorm = maxval(abs(adiag(1:nrow)*arwsum(1:nrow)))
@@ -397,7 +373,6 @@ subroutine amg_d_parmatch_smth_bld(dol1smoothing,ag,a,desc_a,ilaggr,nlaggr,&
 
     end block
   end if
-  if (allocated(l1rwsum)) deallocate(l1rwsum)
   if (do_timings) call psb_toc(idx_phase2)
 
   if (debug_level >= psb_debug_outer_) &

amgprec/impl/aggregator/amg_d_parmatch_unsmth_bld.F90

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

amgprec/impl/aggregator/amg_d_ptap_bld.f90

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

amgprec/impl/aggregator/amg_d_rap.f90

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

amgprec/impl/aggregator/amg_d_soc1_map_bld.F90

@@ -250,7 +250,7 @@ subroutine amg_d_soc1_map_bld(iorder,theta,clean_zeros,a,desc_a,nlaggr,ilaggr,in
               ! we will not reset. 
               if (j>nr) cycle step1
               if (ilaggr(j) > 0) cycle step1
-              if ((abs(val(k)) > theta*sqrt(abs(diag(i)*diag(j)))).and.(diag(i).ne.dzero)) then
+              if (abs(val(k)) > theta*sqrt(abs(diag(i)*diag(j)))) then
                 ip = ip + 1
                 icol(ip) = icol(k)
               end if
@@ -275,8 +275,7 @@ subroutine amg_d_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)) !be careful about overflow
-              itmp = itmp*nths+kk
+              itmp = (bnds(kk)-1+locnaggr(kk))*nths+kk
               if (itmp < (bnds(kk)-1+locnaggr(kk))) then
                 !$omp atomic update
                 info = max(12345678,info)
@@ -357,7 +356,7 @@ subroutine amg_d_soc1_map_bld(iorder,theta,clean_zeros,a,desc_a,nlaggr,ilaggr,in
       do k=1, nz
         j   = icol(k)
         if ((1<=j).and.(j<=nr)) then 
-          if ((abs(val(k)) > theta*sqrt(abs(diag(i)*diag(j)))).and.(diag(i).ne.dzero)) then
+          if (abs(val(k)) > theta*sqrt(abs(diag(i)*diag(j)))) then
             ip = ip + 1
             icol(ip) = icol(k)
           end if
@@ -545,3 +544,4 @@ subroutine amg_d_soc1_map_bld(iorder,theta,clean_zeros,a,desc_a,nlaggr,ilaggr,in
   return
 
 end subroutine amg_d_soc1_map_bld
+

amgprec/impl/aggregator/amg_d_soc2_map_bld.F90

@@ -309,8 +309,7 @@ subroutine amg_d_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)) !be careful about overflow
-              itmp = itmp*nths+kk
+              itmp = (bnds(kk)-1+locnaggr(kk))*nths+kk
               if (itmp < (bnds(kk)-1+locnaggr(kk))) then
                 !$omp atomic update
                 info = max(12345678,info)

amgprec/impl/aggregator/amg_daggrmat_minnrg_bld.f90

@@ -69,7 +69,6 @@
 !
 !
 ! Arguments:
-!    dol1smoothing - fictitious integer argument, it is not used inside
 !    a          -  type(psb_dspmat_type), input.     
 !                  The sparse matrix structure containing the local part of
 !                  the fine-level matrix.
@@ -105,8 +104,8 @@
 !                  Error code.
 !
 !
-subroutine amg_daggrmat_minnrg_bld(dol1smoothing,a,desc_a,ilaggr,nlaggr,&
-                parms,ac,desc_ac,op_prol,op_restr,t_prol,info)
+subroutine amg_daggrmat_minnrg_bld(a,desc_a,ilaggr,nlaggr,parms,&
+     & ac,desc_ac,op_prol,op_restr,t_prol,info)
   use psb_base_mod
   use amg_base_prec_type
   use amg_d_inner_mod, amg_protect_name => amg_daggrmat_minnrg_bld
@@ -114,7 +113,6 @@ subroutine amg_daggrmat_minnrg_bld(dol1smoothing,a,desc_a,ilaggr,nlaggr,&
   implicit none
 
   ! Arguments
-  integer(psb_ipk_), intent(in)                 :: dol1smoothing
   type(psb_dspmat_type), intent(in)           :: a
   type(psb_desc_type), intent(inout)            :: desc_a
   integer(psb_lpk_), intent(inout)              :: ilaggr(:), nlaggr(:)
@@ -173,13 +171,6 @@ subroutine amg_daggrmat_minnrg_bld(dol1smoothing,a,desc_a,ilaggr,nlaggr,&
 
   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
@@ -309,7 +300,7 @@ subroutine amg_daggrmat_minnrg_bld(dol1smoothing,a,desc_a,ilaggr,nlaggr,&
 !!$        endif
 !!$      enddo
 !!$      if (jd == -1) then 
-!!$                write(0,*) name,': Warning: there is no diagonal element', i
+!!$        write(0,*) 'Wrong input: we need the diagonal!!!!', i
 !!$      else
 !!$        acsrf%val(jd)=acsrf%val(jd)-tmp
 !!$      end if

amgprec/impl/aggregator/amg_daggrmat_nosmth_bld.f90

@@ -94,11 +94,10 @@
 !               
 !    info       -  integer, output.
 !                  Error code.
-!    dol1smoothing - optional, this is here just for interfacing reasons. It is not used by the
-!                    code
 !
-subroutine amg_daggrmat_nosmth_bld(dol1smoothing,a,desc_a,ilaggr,nlaggr,&
-                parms,ac,desc_ac,op_prol,op_restr,t_prol,info)
+!
+subroutine amg_daggrmat_nosmth_bld(a,desc_a,ilaggr,nlaggr,parms,&
+     & ac,desc_ac,op_prol,op_restr,t_prol,info)
   use psb_base_mod
   use amg_base_prec_type
   use amg_d_inner_mod, amg_protect_name => amg_daggrmat_nosmth_bld
@@ -106,7 +105,6 @@ subroutine amg_daggrmat_nosmth_bld(dol1smoothing,a,desc_a,ilaggr,nlaggr,&
   implicit none
 
   ! Arguments
-  integer(psb_ipk_), intent(in)              :: dol1smoothing 
   type(psb_dspmat_type), intent(in)        :: a
   type(psb_desc_type), intent(inout)         :: desc_a
   integer(psb_lpk_), intent(inout)           :: ilaggr(:), nlaggr(:)
@@ -139,12 +137,6 @@ subroutine amg_daggrmat_nosmth_bld(dol1smoothing,a,desc_a,ilaggr,nlaggr,&
 
   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()