Merge branch 'dev-openmp' into TestFerdous

(Might have solved stucked runs)

Changelog

@@ -1,5 +1,6 @@
 Changelog. A lot less detailed than usual, at least for past
 	history.
+2022/05/20: Restart ChangeLog. Updated to new name AMG4PSBLAS, now using PSB3.8	
 2018/10/28: Fix interface to MUMPS and configry machinery. Require PSB 3.6.
 2018/10/10: ICTXT argument in prec%init().	
 2018/07/30: Fixes for Intel compilers. BootCMatch interface in examples. 	

LICENSE

@@ -1,10 +1,10 @@
   
    
-                             AMG4PSBLAS version 1.0
+                             AMG4PSBLAS version 1.1
     Algebraic Multigrid Package
-               based on PSBLAS (Parallel Sparse BLAS version 3.7)
+               based on PSBLAS (Parallel Sparse BLAS version 3.8)
     
-    (C) Copyright 2021
+    (C) Copyright 2022
   
         Salvatore Filippone  
         Pasqua D'Ambra   

Make.inc.in

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

Makefile

@@ -1,10 +1,13 @@
 include Make.inc
 
 
-all:  library 
+all:  objs lib
 
-library: libdir amgp cbnd
-#cbnd
+objs: amgp cbnd
+
+lib: libdir objs
+	cd amgprec && $(MAKE) lib
+	cd cbind && $(MAKE)  lib
 
 libdir:
 	(if test ! -d lib ; then mkdir lib; fi)
@@ -14,10 +17,11 @@ libdir:
         
 
 amgp:
-	$(MAKE) -C amgprec all
+	cd amgprec && $(MAKE) objs
 cbnd: amgp
-	$(MAKE) -C cbind all
-install: all
+	cd cbind && $(MAKE)  objs
+
+install: lib
 	mkdir -p $(INSTALL_LIBDIR) &&\
 	   $(INSTALL_DATA) lib/*.a  $(INSTALL_LIBDIR)
 	mkdir -p  $(INSTALL_INCLUDEDIR) &&\
@@ -41,14 +45,14 @@ cleanlib:
 	(cd modules; /bin/rm -f *.a *$(.mod) *$(.fh))
 
 veryclean: cleanlib
-	(cd amgprec; make veryclean)
-	(cd samples/simple/fileread; make clean)
-	(cd samples/simple/pdegen; make clean)
-	(cd samples/advanced/fileread; make clean)
-	(cd samples/advanced/pdegen; make clean)
+	(cd amgprec && $(MAKE) veryclean)
+	(cd samples/simple/fileread && $(MAKE) clean)
+	(cd samples/simple/pdegen && $(MAKE) clean)
+	(cd samples/advanced/fileread && $(MAKE) clean)
+	(cd samples/advanced/pdegen && $(MAKE) clean)
 
 check: all
 	make check -C samples/advanced/pdegen
 
 clean:
-	(cd amgprec; make clean)
+	(cd amgprec && $(MAKE) clean)

README.md

@@ -1,6 +1,5 @@
-
                          AMG4PSBLAS  
- Algebraic Multigrid Package  based on PSBLAS (Parallel Sparse BLAS version 3.7)
+ 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)

amgprec/Makefile

@@ -16,7 +16,8 @@ DMODOBJS=amg_d_prec_type.o \
 	amg_d_dec_aggregator_mod.o amg_d_symdec_aggregator_mod.o \
 	amg_d_ainv_solver.o amg_d_base_ainv_mod.o \
 	amg_d_invk_solver.o amg_d_invt_solver.o amg_d_krm_solver.o \
-	amg_d_matchboxp_mod.o amg_d_parmatch_aggregator_mod.o
+	amg_d_matchboxp_mod.o amg_d_parmatch_aggregator_mod.o \
+        amg_d_newmatch_aggregator_mod.o amg_d_decmatch_mod.o
 
 SMODOBJS=amg_s_prec_type.o amg_s_ilu_fact_mod.o \
 	amg_s_inner_mod.o amg_s_ilu_solver.o amg_s_diag_solver.o amg_s_jac_smoother.o amg_s_as_smoother.o \
@@ -62,17 +63,20 @@ OBJS=$(MODOBJS)
 LOCAL_MODS=$(MODOBJS:.o=$(.mod))
 LIBNAME=libamg_prec.a
 
-all: lib impld
+all: objs impld
 
-impld: $(OBJS)
-	$(MAKE) -C impl
+objs: $(OBJS)
+	/bin/cp -p amg_const.h $(INCDIR)
+	/bin/cp -p *$(.mod) $(MODDIR)
+
+impld: objs
+	cd impl && $(MAKE) 
 
 lib: $(OBJS) impld
+	cd impl && $(MAKE) lib
 	$(AR) $(HERE)/$(LIBNAME) $(OBJS)
 	$(RANLIB) $(HERE)/$(LIBNAME)
 	/bin/cp -p $(HERE)/$(LIBNAME) $(LIBDIR)
-	/bin/cp -p amg_const.h $(INCDIR)
-	/bin/cp -p *$(.mod) $(MODDIR)
 
 
 $(MODOBJS): $(PSBLAS_MODDIR)/$(PSBBASEMODNAME)$(.mod)
@@ -113,7 +117,7 @@ amg_c_prec_type.o: amg_c_onelev_mod.o
 amg_z_prec_type.o: amg_z_onelev_mod.o
 
 amg_s_onelev_mod.o: amg_s_base_smoother_mod.o amg_s_dec_aggregator_mod.o amg_s_parmatch_aggregator_mod.o 
-amg_d_onelev_mod.o: amg_d_base_smoother_mod.o amg_d_dec_aggregator_mod.o amg_d_parmatch_aggregator_mod.o 
+amg_d_onelev_mod.o: amg_d_base_smoother_mod.o amg_d_dec_aggregator_mod.o amg_d_parmatch_aggregator_mod.o amg_d_newmatch_aggregator_mod.o
 amg_c_onelev_mod.o: amg_c_base_smoother_mod.o amg_c_dec_aggregator_mod.o
 amg_z_onelev_mod.o: amg_z_base_smoother_mod.o amg_z_dec_aggregator_mod.o
 
@@ -126,6 +130,8 @@ amg_d_base_aggregator_mod.o: amg_base_prec_type.o
 amg_d_parmatch_aggregator_mod.o amg_d_dec_aggregator_mod.o: amg_d_base_aggregator_mod.o
 amg_d_hybrid_aggregator_mod.o amg_d_symdec_aggregator_mod.o: amg_d_dec_aggregator_mod.o
 amg_d_parmatch_aggregator_mod.o: amg_d_matchboxp_mod.o 
+amg_d_newmatch_aggregator_mod.o: amg_d_base_aggregator_mod.o
+amg_d_newmatch_aggregator_mod.o: amg_d_decmatch_mod.o
 
 amg_c_base_aggregator_mod.o: amg_base_prec_type.o
 amg_c_parmatch_aggregator_mod.o amg_c_dec_aggregator_mod.o: amg_c_base_aggregator_mod.o
@@ -218,4 +224,4 @@ clean: implclean
 	/bin/rm -f $(OBJS) $(LOCAL_MODS) *$(.mod)
 
 implclean:
-	$(MAKE) -C impl clean
+	cd impl && $(MAKE) clean

amgprec/amg_base_prec_type.F90

@@ -64,7 +64,7 @@ module amg_base_prec_type
   !
   use psb_const_mod
   use psb_base_mod, only :&
-       & psb_desc_type, psb_i_vect_type, psb_i_base_vect_type,&
+       & psb_desc_type, psb_ctxt_type,&
        & psb_ipk_, psb_dpk_, psb_spk_, psb_epk_,  &
        & psb_cdfree, psb_halo_, psb_none_, psb_sum_, psb_avg_, &
        & psb_nohalo_, psb_square_root_, psb_toupper, psb_root_,&
@@ -81,10 +81,10 @@ module amg_base_prec_type
   !
   ! Version numbers
   !
-  character(len=*), parameter   :: amg_version_string_ = "1.0.1"
+  character(len=*), parameter   :: amg_version_string_ = "1.1.0"
   integer(psb_ipk_), parameter  :: amg_version_major_  = 1
-  integer(psb_ipk_), parameter  :: amg_version_minor_  = 0
-  integer(psb_ipk_), parameter  :: amg_patchlevel_     = 1
+  integer(psb_ipk_), parameter  :: amg_version_minor_  = 1
+  integer(psb_ipk_), parameter  :: amg_patchlevel_     = 0
 
   type amg_ml_parms
     integer(psb_ipk_) :: sweeps_pre, sweeps_post
@@ -275,7 +275,8 @@ module amg_base_prec_type
   integer(psb_ipk_), parameter :: amg_sym_dec_aggr_  = 1
   integer(psb_ipk_), parameter :: amg_ext_aggr_      = 2
   integer(psb_ipk_), parameter :: amg_coupled_aggr_  = 3
-  integer(psb_ipk_), parameter :: amg_max_par_aggr_alg_  = amg_coupled_aggr_
+  integer(psb_ipk_), parameter :: amg_newmtc_aggr_   = 4
+  integer(psb_ipk_), parameter :: amg_max_par_aggr_alg_  = amg_newmtc_aggr_
   !
   ! Legal values for entry: amg_aggr_type_
   !
@@ -283,6 +284,7 @@ module amg_base_prec_type
   integer(psb_ipk_), parameter :: amg_soc1_        = 1
   integer(psb_ipk_), parameter :: amg_soc2_        = 2
   integer(psb_ipk_), parameter :: amg_matchboxp_   = 3
+  integer(psb_ipk_), parameter :: amg_newmatch_    = 4
   !
   ! Legal values for entry: amg_aggr_prol_
   !
@@ -371,13 +373,14 @@ module amg_base_prec_type
   character(len=15), parameter, private :: &
        &  matrix_names(0:1)=(/'distributed   ','replicated    '/)
   character(len=18), parameter, private :: &
-       &  aggr_type_names(0:3)=(/'None              ',&
+       &  aggr_type_names(0:4)=(/'None              ',&
        &  'SOC measure 1     ', 'SOC Measure 2     ',&
-       &  'Parallel Matching '/)
+       &  'Parallel Matching ','Decoupled Matching'/)
   character(len=18), parameter, private :: &
-       &  par_aggr_alg_names(0:3)=(/&
+       &  par_aggr_alg_names(0:4)=(/&
        & 'decoupled aggr.   ', 'sym. dec. aggr.   ',&
-       & 'user defined aggr.', 'coupled aggr.     '/)
+       & 'user defined aggr.', 'coupled aggr.     ',&
+       & 'new matching aggr.'/)
   character(len=18), parameter, private :: &
        &  ord_names(0:1)=(/'Natural ordering  ','Desc. degree ord. '/)
   character(len=6), parameter, private :: &
@@ -516,6 +519,10 @@ contains
       val = amg_soc2_
     case('SOC1')
       val = amg_soc1_
+    case('NEWMATCH')
+      val = amg_newmatch_
+    case('NEWMTC')
+      val = amg_newmtc_aggr_
     case('MATCHBOXP','PARMATCH')
       val = amg_matchboxp_
     case('COUPLED','COUP')
@@ -649,43 +656,52 @@ contains
     end if
   end subroutine ml_parms_mlcycledsc
 
-  subroutine ml_parms_mldescr(pm,iout,info)
+  subroutine ml_parms_mldescr(pm,iout,info,prefix)
 
     Implicit None
 
     ! Arguments
-    class(amg_ml_parms), intent(in) :: pm
-    integer(psb_ipk_), intent(in)             :: iout
-    integer(psb_ipk_), intent(out)            :: info
+    class(amg_ml_parms), intent(in)        :: pm
+    integer(psb_ipk_), intent(in)          :: iout
+    integer(psb_ipk_), intent(out)         :: info
+    character(len=*), intent(in), optional :: prefix
+
+    character(1024)  :: prefix_ 
     info = psb_success_
+    if (present(prefix)) then
+      prefix_ = prefix
+    else
+      prefix_ = ""
+    end if
+    
     if ((pm%ml_cycle>=amg_no_ml_).and.(pm%ml_cycle<=amg_max_ml_cycle_)) then
 
 
-      write(iout,*) '  Parallel aggregation algorithm: ',&
+      write(iout,*) trim(prefix),'  Parallel aggregation algorithm: ',&
            &   par_aggr_alg_names(pm%par_aggr_alg)
-      if (pm%aggr_type>0) write(iout,*) '  Aggregation type: ',&
+      if (pm%aggr_type>0) write(iout,*) trim(prefix),'  Aggregation type: ',&
            & aggr_type_names(pm%aggr_type)
       !if (pm%par_aggr_alg /= amg_ext_aggr_) then
         if ( pm%aggr_ord /= amg_aggr_ord_nat_) &
-             & write(iout,*) '               with initial ordering: ',&
+             & write(iout,*) trim(prefix),'               with initial ordering: ',&
              &   ord_names(pm%aggr_ord)
-        write(iout,*) '  Aggregation prolongator: ', &
+        write(iout,*) trim(prefix),'  Aggregation prolongator: ', &
              &  aggr_prols(pm%aggr_prol)
         if (pm%aggr_prol /= amg_no_smooth_) then
-          write(iout,*) '              with: ', aggr_filters(pm%aggr_filter)
+          write(iout,*) trim(prefix),'              with: ', aggr_filters(pm%aggr_filter)
           if (pm%aggr_omega_alg == amg_eig_est_) then
-            write(iout,*) '  Damping omega computation: spectral radius estimate'
-            write(iout,*) '  Spectral radius estimate: ', &
+            write(iout,*) trim(prefix),'  Damping omega computation: spectral radius estimate'
+            write(iout,*) trim(prefix),'  Spectral radius estimate: ', &
                  & eigen_estimates(pm%aggr_eig)
           else if (pm%aggr_omega_alg == amg_user_choice_) then
-            write(iout,*) '  Damping omega computation: user defined value.'
+            write(iout,*) trim(prefix),'  Damping omega computation: user defined value.'
           else
-            write(iout,*) '  Damping omega computation: unknown value in iprcparm!!'
+            write(iout,*) trim(prefix),'  Damping omega computation: unknown value in iprcparm!!'
           end if
         end if
       !end if
     else
-      write(iout,*) '  Multilevel type: Unkonwn value. Something is amiss....',&
+      write(iout,*) trim(prefix),'  Multilevel type: Unkonwn value. Something is amiss....',&
            & pm%ml_cycle
     end if
 
@@ -693,15 +709,16 @@ contains
 
   end subroutine ml_parms_mldescr
 
-  subroutine ml_parms_descr(pm,iout,info,coarse)
+  subroutine ml_parms_descr(pm,iout,info,coarse,prefix)
 
     Implicit None
 
     ! Arguments
-    class(amg_ml_parms), intent(in) :: pm
-    integer(psb_ipk_), intent(in)             :: iout
-    integer(psb_ipk_), intent(out)            :: info
-    logical, intent(in), optional   :: coarse
+    class(amg_ml_parms), intent(in)        :: pm
+    integer(psb_ipk_), intent(in)          :: iout
+    integer(psb_ipk_), intent(out)         :: info
+    logical, intent(in), optional          :: coarse
+    character(len=*), intent(in), optional :: prefix
     logical :: coarse_
 
     info = psb_success_
@@ -712,7 +729,7 @@ contains
     end if
 
     if (coarse_) then
-      call pm%coarsedescr(iout,info)
+      call pm%coarsedescr(iout,info,prefix=prefix)
     end if
 
     return
@@ -720,81 +737,126 @@ contains
   end subroutine ml_parms_descr
 
 
-  subroutine ml_parms_coarsedescr(pm,iout,info)
+  subroutine ml_parms_coarsedescr(pm,iout,info,prefix)
 
 
     Implicit None
 
     ! Arguments
-    class(amg_ml_parms), intent(in) :: pm
-    integer(psb_ipk_), intent(in)             :: iout
-    integer(psb_ipk_), intent(out)            :: info
+    class(amg_ml_parms), intent(in)       :: pm
+    integer(psb_ipk_), intent(in)         :: iout
+    integer(psb_ipk_), intent(out)        :: info
+    character(len=*), intent(in), optional :: prefix
+
+    character(1024)  :: prefix_ 
 
     info = psb_success_
-    write(iout,*) '  Coarse matrix: ',&
+    if (present(prefix)) then
+      prefix_ = prefix
+    else
+      prefix_ = ""
+    end if
+
+    write(iout,*) trim(prefix),'  Coarse matrix: ',&
          & matrix_names(pm%coarse_mat)
     select case(pm%coarse_solve)
     case (amg_bjac_,amg_as_)
-      write(iout,*) '  Number of sweeps : ',&
-           & pm%sweeps_pre
-      write(iout,*) '  Coarse solver: ',&
+      write(iout,*) trim(prefix),'  Coarse solver: ',&
            & 'Block Jacobi'
+      write(iout,*) trim(prefix),'  Number of sweeps : ',&
+           & pm%sweeps_pre
     case (amg_l1_bjac_)
-      write(iout,*) '  Number of sweeps : ',&
-           & pm%sweeps_pre
-      write(iout,*) '  Coarse solver: ',&
+      write(iout,*) trim(prefix),'  Coarse solver: ',&
            & 'L1-Block Jacobi'
-    case (amg_jac_)
-      write(iout,*) '  Number of sweeps : ',&
+      write(iout,*) trim(prefix),'  Number of sweeps : ',&
            & pm%sweeps_pre
-      write(iout,*) '  Coarse solver: ',&
+    case (amg_jac_)
+      write(iout,*) trim(prefix),'  Coarse solver: ',&
            & 'Point Jacobi'
+      write(iout,*) trim(prefix),'  Number of sweeps : ',&
+           & pm%sweeps_pre
+    case (amg_l1_jac_)
+      write(iout,*) trim(prefix),'  Coarse solver: ',&
+           & 'L1-Jacobi'
+      write(iout,*) trim(prefix),'  Number of sweeps : ',&
+           & pm%sweeps_pre
+    case (amg_l1_fbgs_)
+      write(iout,*) trim(prefix),'  Coarse solver: ',&
+           & 'L1 Forward-Backward Gauss-Seidel (Hybrid)'
+      write(iout,*) trim(prefix),'  Number of sweeps : ',&
+           & pm%sweeps_pre
+    case (amg_l1_gs_)
+      write(iout,*) trim(prefix),'  Coarse solver: ',&
+           & 'L1 Gauss-Seidel (Hybrid)'
+      write(iout,*) trim(prefix),'  Number of sweeps : ',&
+           & pm%sweeps_pre
+    case (amg_fbgs_)
+      write(iout,*) trim(prefix),'  Coarse solver: ',&
+           & 'Forward-Backward Gauss-Seidel (Hybrid)'
+      write(iout,*) trim(prefix),'  Number of sweeps : ',&
+           & pm%sweeps_pre
     case default
-      write(iout,*) '  Coarse solver: ',&
+      write(iout,*) trim(prefix),'  Coarse solver: ',&
            & amg_fact_names(pm%coarse_solve)
     end select
-
+    
   end subroutine ml_parms_coarsedescr
 
-  subroutine s_ml_parms_descr(pm,iout,info,coarse)
+  subroutine s_ml_parms_descr(pm,iout,info,coarse,prefix)
 
     Implicit None
 
     ! Arguments
-    class(amg_sml_parms), intent(in) :: pm
-    integer(psb_ipk_), intent(in)             :: iout
-    integer(psb_ipk_), intent(out)            :: info
-    logical, intent(in), optional   :: coarse
+    class(amg_sml_parms), intent(in)       :: pm
+    integer(psb_ipk_), intent(in)          :: iout
+    integer(psb_ipk_), intent(out)         :: info
+    logical, intent(in), optional          :: coarse
+    character(len=*), intent(in), optional :: prefix
+
+    character(1024)  :: prefix_ 
 
     info = psb_success_
-
-    call pm%amg_ml_parms%descr(iout,info,coarse)
-    if (pm%aggr_prol /= amg_no_smooth_) then
-      write(iout,*) '  Damping omega value  :',pm%aggr_omega_val
+    if (present(prefix)) then
+      prefix_ = prefix
+    else
+      prefix_ = ""
     end if
-    write(iout,*) '  Aggregation threshold:',pm%aggr_thresh
+
+    call pm%amg_ml_parms%descr(iout,info,coarse,prefix=prefix)
+    if (pm%aggr_prol /= amg_no_smooth_) then
+      write(iout,*) trim(prefix),'  Damping omega value  :',pm%aggr_omega_val
+    end if
+    write(iout,*) trim(prefix),'  Aggregation threshold:',pm%aggr_thresh
 
     return
 
   end subroutine s_ml_parms_descr
 
-  subroutine d_ml_parms_descr(pm,iout,info,coarse)
+  subroutine d_ml_parms_descr(pm,iout,info,coarse,prefix)
 
     Implicit None
 
     ! Arguments
-    class(amg_dml_parms), intent(in) :: pm
-    integer(psb_ipk_), intent(in)             :: iout
-    integer(psb_ipk_), intent(out)            :: info
-    logical, intent(in), optional   :: coarse
+    class(amg_dml_parms), intent(in)       :: pm
+    integer(psb_ipk_), intent(in)          :: iout
+    integer(psb_ipk_), intent(out)         :: info
+    logical, intent(in), optional          :: coarse
+    character(len=*), intent(in), optional :: prefix
+
+    character(1024)  :: prefix_ 
 
     info = psb_success_
-
-    call pm%amg_ml_parms%descr(iout,info,coarse)
-    if (pm%aggr_prol /= amg_no_smooth_) then
-      write(iout,*) '  Damping omega value  :',pm%aggr_omega_val
+    if (present(prefix)) then
+      prefix_ = prefix
+    else
+      prefix_ = ""
     end if
-    write(iout,*) '  Aggregation threshold:',pm%aggr_thresh
+
+    call pm%amg_ml_parms%descr(iout,info,coarse,prefix=prefix)
+    if (pm%aggr_prol /= amg_no_smooth_) then
+      write(iout,*) trim(prefix),'  Damping omega value  :',pm%aggr_omega_val
+    end if
+    write(iout,*) trim(prefix),'  Aggregation threshold:',pm%aggr_thresh
 
     return
 

amgprec/amg_c_ainv_solver.F90

@@ -198,7 +198,7 @@ module amg_c_ainv_solver
 !!$  end interface
 
   interface
-    subroutine amg_c_ainv_solver_descr(sv,info,iout,coarse)
+    subroutine amg_c_ainv_solver_descr(sv,info,iout,coarse,prefix)
       import :: psb_dpk_, amg_c_ainv_solver_type, psb_ipk_
 
       Implicit None
@@ -208,7 +208,7 @@ module amg_c_ainv_solver
       integer(psb_ipk_), intent(out)            :: info
       integer(psb_ipk_), intent(in), optional   :: iout
       logical, intent(in), optional             :: coarse
-
+      character(len=*), intent(in), optional  :: prefix
     end subroutine amg_c_ainv_solver_descr
   end interface
 

amgprec/amg_c_as_smoother.f90

@@ -396,21 +396,23 @@ contains
   end subroutine c_as_smoother_default
 
 
-  subroutine c_as_smoother_descr(sm,info,iout,coarse)
+  subroutine c_as_smoother_descr(sm,info,iout,coarse,prefix)
 
     Implicit None
 
     ! Arguments
     class(amg_c_as_smoother_type), intent(in) :: sm
-    integer(psb_ipk_), intent(out)                      :: info
-    integer(psb_ipk_), intent(in), optional             :: iout
-    logical, intent(in), optional             :: coarse
+    integer(psb_ipk_), intent(out)              :: info
+    integer(psb_ipk_), intent(in), optional     :: iout
+    logical, intent(in), optional              :: coarse
+      character(len=*), intent(in), optional   :: prefix
 
     ! Local variables
     integer(psb_ipk_)      :: err_act
     character(len=20), parameter :: name='amg_c_as_smoother_descr'
     integer(psb_ipk_) :: iout_
     logical      :: coarse_
+    character(1024)    :: prefix_
 
     call psb_erractionsave(err_act)
     info = psb_success_
@@ -424,16 +426,21 @@ contains
     else
       iout_ = psb_out_unit
     endif
+    if (present(prefix)) then
+      prefix_ = prefix
+    else
+      prefix_ = ""
+    end if
 
     if (.not.coarse_) then 
-      write(iout_,*) '  Additive Schwarz with  ',&
+      write(iout_,*) trim(prefix_), '  Additive Schwarz with  ',&
            &  sm%novr, ' overlap layers.'
-      write(iout_,*) '  Restrictor:  ',restrict_names(sm%restr)
-      write(iout_,*) '  Prolongator: ',prolong_names(sm%prol)
-      write(iout_,*) '  Local solver:'
+      write(iout_,*) trim(prefix_), '  Restrictor:  ',restrict_names(sm%restr)
+      write(iout_,*) trim(prefix_), '  Prolongator: ',prolong_names(sm%prol)
+      write(iout_,*) trim(prefix_), '  Local solver:'
     endif
     if (allocated(sm%sv)) then 
-      call sm%sv%descr(info,iout_,coarse=coarse)
+      call sm%sv%descr(info,iout_,coarse=coarse,prefix=prefix)
     end if
 
     call psb_erractionrestore(err_act)

amgprec/amg_c_base_aggregator_mod.f90

@@ -275,15 +275,22 @@ contains
     val = .false.
   end function amg_c_base_aggregator_xt_desc
 
-  subroutine  amg_c_base_aggregator_descr(ag,parms,iout,info)
+  subroutine  amg_c_base_aggregator_descr(ag,parms,iout,info,prefix)
     implicit none 
     class(amg_c_base_aggregator_type), intent(in) :: ag
     type(amg_sml_parms), intent(in)   :: parms
     integer(psb_ipk_), intent(in)  :: iout
     integer(psb_ipk_), intent(out) :: info
+    character(len=*), intent(in), optional  :: prefix
+    character(1024)    :: prefix_
+    if (present(prefix)) then
+      prefix_ = prefix
+    else
+      prefix_ = ""
+    end if
 
-    write(iout,*) 'Aggregator object type: ',ag%fmt()
-    call parms%mldescr(iout,info)
+    write(iout,*) trim(prefix_),' ', 'Aggregator object type: ',ag%fmt()
+    call parms%mldescr(iout,info,prefix=prefix)
     
     return
   end subroutine amg_c_base_aggregator_descr

amgprec/amg_c_base_smoother_mod.f90

@@ -272,7 +272,7 @@ module amg_c_base_smoother_mod
   end interface
   
   interface 
-    subroutine amg_c_base_smoother_descr(sm,info,iout,coarse)
+    subroutine amg_c_base_smoother_descr(sm,info,iout,coarse,prefix)
       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_base_smoother_type, psb_ipk_
@@ -281,6 +281,7 @@ module amg_c_base_smoother_mod
       integer(psb_ipk_), intent(out)                :: info
       integer(psb_ipk_), intent(in), optional       :: iout
       logical, intent(in), optional                 :: coarse
+      character(len=*), intent(in), optional        :: prefix
     end subroutine amg_c_base_smoother_descr
   end interface
   

amgprec/amg_c_base_solver_mod.f90

@@ -270,7 +270,7 @@ module amg_c_base_solver_mod
   end interface
   
   interface
-    subroutine amg_c_base_solver_descr(sv,info,iout,coarse)
+    subroutine amg_c_base_solver_descr(sv,info,iout,coarse,prefix)
       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_base_solver_type, psb_ipk_
@@ -281,7 +281,7 @@ module amg_c_base_solver_mod
       integer(psb_ipk_), intent(out)              :: info
       integer(psb_ipk_), intent(in), optional     :: iout
       logical, intent(in), optional               :: coarse
-
+      character(len=*), intent(in), optional      :: prefix
     end subroutine amg_c_base_solver_descr
   end interface 
   

amgprec/amg_c_dec_aggregator_mod.f90

@@ -184,16 +184,23 @@ contains
     val = "Decoupled aggregation"
   end function amg_c_dec_aggregator_fmt
   
-  subroutine  amg_c_dec_aggregator_descr(ag,parms,iout,info)
+  subroutine  amg_c_dec_aggregator_descr(ag,parms,iout,info,prefix)
     implicit none 
     class(amg_c_dec_aggregator_type), intent(in) :: ag
     type(amg_sml_parms), intent(in)   :: parms
     integer(psb_ipk_), intent(in)  :: iout
     integer(psb_ipk_), intent(out) :: info
+    character(len=*), intent(in), optional  :: prefix
+    character(1024)    :: prefix_
+    if (present(prefix)) then
+      prefix_ = prefix
+    else
+      prefix_ = ""
+    end if
 
-    write(iout,*) 'Decoupled Aggregator'
-    write(iout,*) 'Aggregator object type: ',ag%fmt()
-    call parms%mldescr(iout,info)
+    write(iout,*) trim(prefix_),' ','Decoupled Aggregator'
+    write(iout,*) trim(prefix_),' ','Aggregator object type: ',ag%fmt()
+    call parms%mldescr(iout,info,prefix=prefix)
     
     return
   end subroutine amg_c_dec_aggregator_descr

amgprec/amg_c_diag_solver.f90

@@ -219,7 +219,7 @@ contains
 
   end subroutine c_diag_solver_free
 
-  subroutine c_diag_solver_descr(sv,info,iout,coarse)
+  subroutine c_diag_solver_descr(sv,info,iout,coarse,prefix)
 
     Implicit None
 
@@ -228,11 +228,13 @@ contains
     integer(psb_ipk_), intent(out)              :: info
     integer(psb_ipk_), intent(in), optional     :: iout
     logical, intent(in), optional               :: coarse
+    character(len=*), intent(in), optional      :: prefix
 
     ! Local variables
     integer(psb_ipk_)      :: err_act
     character(len=20), parameter :: name='amg_c_diag_solver_descr'
     integer(psb_ipk_) :: iout_
+    character(1024)    :: prefix_
 
     info = psb_success_
     if (present(iout)) then 
@@ -240,8 +242,13 @@ contains
     else
       iout_ = psb_out_unit
     endif
-    
-    write(iout_,*) '  Diagonal local solver '
+    if (present(prefix)) then
+      prefix_ = prefix
+    else
+      prefix_ = ""
+    end if
+
+    write(iout_,*) trim(prefix_), '  Diagonal local solver '
 
     return
 
@@ -352,7 +359,7 @@ module amg_c_l1_diag_solver
   
 contains
 
-  subroutine c_l1_diag_solver_descr(sv,info,iout,coarse)
+  subroutine c_l1_diag_solver_descr(sv,info,iout,coarse,prefix)
 
     Implicit None
 
@@ -361,11 +368,13 @@ contains
     integer(psb_ipk_), intent(out)              :: info
     integer(psb_ipk_), intent(in), optional     :: iout
     logical, intent(in), optional               :: coarse
+    character(len=*), intent(in), optional  :: prefix
 
     ! Local variables
     integer(psb_ipk_)      :: err_act
     character(len=20), parameter :: name='amg_c_l1_diag_solver_descr'
     integer(psb_ipk_) :: iout_
+    character(1024)    :: prefix_
 
     info = psb_success_
     if (present(iout)) then 
@@ -373,8 +382,13 @@ contains
     else
       iout_ = psb_out_unit
     endif
-    
-    write(iout_,*) '  L1 Diagonal solver '
+    if (present(prefix)) then
+      prefix_ = prefix
+    else
+      prefix_ = ""
+    end if
+
+    write(iout_,*) trim(prefix_), '  L1 Diagonal solver '
 
     return
 

amgprec/amg_c_gs_solver.f90

@@ -433,20 +433,22 @@ contains
     return
   end subroutine c_gs_solver_free
 
-  subroutine c_gs_solver_descr(sv,info,iout,coarse)
+  subroutine c_gs_solver_descr(sv,info,iout,coarse,prefix)
 
     Implicit None
 
     ! Arguments
     class(amg_c_gs_solver_type), intent(in) :: sv
-    integer(psb_ipk_), intent(out)             :: info
-    integer(psb_ipk_), intent(in), optional    :: iout
-    logical, intent(in), optional       :: coarse
+    integer(psb_ipk_), intent(out)            :: info
+    integer(psb_ipk_), intent(in), optional   :: iout
+    logical, intent(in), optional             :: coarse
+    character(len=*), intent(in), optional    :: prefix
 
     ! Local variables
     integer(psb_ipk_)      :: err_act
     character(len=20), parameter :: name='amg_c_gs_solver_descr'
     integer(psb_ipk_) :: iout_
+    character(1024)    :: prefix_
 
     call psb_erractionsave(err_act)
     info = psb_success_
@@ -455,12 +457,17 @@ contains
     else
       iout_ = psb_out_unit
     endif
+    if (present(prefix)) then
+      prefix_ = prefix
+    else
+      prefix_ = ""
+    end if
 
     if (sv%eps<=dzero) then 
-      write(iout_,*) '  Forward Gauss-Seidel iterative solver with  ',&
+      write(iout_,*) trim(prefix_), '  Forward Gauss-Seidel iterative solver with  ',&
            &  sv%sweeps,' sweeps'
     else
-      write(iout_,*) '  Forward Gauss-Seidel iterative solver with  tolerance',&
+      write(iout_,*) trim(prefix_), '  Forward Gauss-Seidel iterative solver with  tolerance',&
            &  sv%eps,' and maxit', sv%sweeps
     end if
     
@@ -526,20 +533,22 @@ contains
     val = .true.
   end function c_gs_solver_is_iterative
     
-  subroutine c_bwgs_solver_descr(sv,info,iout,coarse)
+  subroutine c_bwgs_solver_descr(sv,info,iout,coarse,prefix)
 
     Implicit None
 
     ! Arguments
     class(amg_c_bwgs_solver_type), intent(in) :: sv
-    integer(psb_ipk_), intent(out)             :: info
-    integer(psb_ipk_), intent(in), optional    :: iout
-    logical, intent(in), optional       :: coarse
+    integer(psb_ipk_), intent(out)              :: info
+    integer(psb_ipk_), intent(in), optional     :: iout
+    logical, intent(in), optional               :: coarse
+    character(len=*), intent(in), optional      :: prefix
 
     ! Local variables
     integer(psb_ipk_)      :: err_act
     character(len=20), parameter :: name='amg_c_bwgs_solver_descr'
     integer(psb_ipk_) :: iout_
+    character(1024)    :: prefix_
 
     call psb_erractionsave(err_act)
     info = psb_success_
@@ -548,12 +557,17 @@ contains
     else
       iout_ = psb_out_unit
     endif
+    if (present(prefix)) then
+      prefix_ = prefix
+    else
+      prefix_ = ""
+    end if
 
     if (sv%eps<=dzero) then 
-      write(iout_,*) '  Backward Gauss-Seidel iterative solver with  ',&
+      write(iout_,*) trim(prefix_), '  Backward Gauss-Seidel iterative solver with  ',&
            &  sv%sweeps,' sweeps'
     else
-      write(iout_,*) '  Backward Gauss-Seidel iterative solver with  tolerance',&
+      write(iout_,*) trim(prefix_), '  Backward Gauss-Seidel iterative solver with  tolerance',&
            &  sv%eps,' and maxit', sv%sweeps
     end if
     

amgprec/amg_c_id_solver.f90

@@ -157,7 +157,7 @@ contains
     return
   end subroutine c_id_solver_free
 
-  subroutine c_id_solver_descr(sv,info,iout,coarse)
+  subroutine c_id_solver_descr(sv,info,iout,coarse,prefix)
 
     Implicit None
 
@@ -165,12 +165,14 @@ contains
     class(amg_c_id_solver_type), intent(in) :: sv
     integer(psb_ipk_), intent(out)            :: info
     integer(psb_ipk_), intent(in), optional   :: iout
-    logical, intent(in), optional       :: coarse
+    logical, intent(in), optional             :: coarse
+    character(len=*), intent(in), optional    :: prefix
 
     ! Local variables
     integer(psb_ipk_)      :: err_act
     character(len=20), parameter :: name='amg_c_id_solver_descr'
     integer(psb_ipk_) :: iout_
+    character(1024)    :: prefix_
 
     info = psb_success_
     if (present(iout)) then 
@@ -178,8 +180,13 @@ contains
     else
       iout_ = psb_out_unit
     endif
+    if (present(prefix)) then
+      prefix_ = prefix
+    else
+      prefix_ = ""
+    end if
     
-    write(iout_,*) '  Identity local solver '
+    write(iout_,*) trim(prefix_), '  Identity local solver '
 
     return
 

amgprec/amg_c_ilu_solver.f90

@@ -406,7 +406,7 @@ contains
     return
   end subroutine c_ilu_solver_free
 
-  subroutine c_ilu_solver_descr(sv,info,iout,coarse)
+  subroutine c_ilu_solver_descr(sv,info,iout,coarse,prefix)
 
     Implicit None
 
@@ -414,12 +414,14 @@ contains
     class(amg_c_ilu_solver_type), intent(in) :: sv
     integer(psb_ipk_), intent(out)             :: info
     integer(psb_ipk_), intent(in), optional    :: iout
-    logical, intent(in), optional       :: coarse
-
+    logical, intent(in), optional              :: coarse
+    character(len=*), intent(in), optional  :: prefix
+      
     ! Local variables
     integer(psb_ipk_) :: err_act
     character(len=20), parameter :: name='amg_c_ilu_solver_descr'
     integer(psb_ipk_) :: iout_
+    character(1024)    :: prefix_
 
     call psb_erractionsave(err_act)
     info = psb_success_
@@ -428,15 +430,20 @@ contains
     else
       iout_ = psb_out_unit
     endif
+    if (present(prefix)) then
+      prefix_ = prefix
+    else
+      prefix_ = ""
+    end if
 
-    write(iout_,*) '  Incomplete factorization solver: ',&
+    write(iout_,*) trim(prefix_), '  Incomplete factorization solver: ',&
          &  amg_fact_names(sv%fact_type)
     select case(sv%fact_type)
     case(psb_ilu_n_,psb_milu_n_)      
-      write(iout_,*) '  Fill level:',sv%fill_in
+      write(iout_,*) trim(prefix_), '  Fill level:',sv%fill_in
     case(psb_ilu_t_)         
-      write(iout_,*) '  Fill level:',sv%fill_in
-      write(iout_,*) '  Fill threshold :',sv%thresh
+      write(iout_,*) trim(prefix_), '  Fill level:',sv%fill_in
+      write(iout_,*) trim(prefix_), '  Fill threshold :',sv%thresh
     end select
 
     call psb_erractionrestore(err_act)

amgprec/amg_c_invk_solver.f90

@@ -123,7 +123,7 @@ module amg_c_invk_solver
   end interface
 
   interface
-    subroutine amg_c_invk_solver_descr(sv,info,iout,coarse)
+    subroutine amg_c_invk_solver_descr(sv,info,iout,coarse,prefix)
       import :: psb_spk_, amg_c_invk_solver_type, psb_ipk_
 
       Implicit None
@@ -133,7 +133,7 @@ module amg_c_invk_solver
       integer(psb_ipk_), intent(out)            :: info
       integer(psb_ipk_), intent(in), optional   :: iout
       logical, intent(in), optional             :: coarse
-
+      character(len=*), intent(in), optional  :: prefix
     end subroutine amg_c_invk_solver_descr
   end interface
 

amgprec/amg_c_invt_solver.f90

@@ -134,16 +134,17 @@ module amg_c_invt_solver
   end interface
 
   interface
-    subroutine amg_c_invt_solver_descr(sv,info,iout,coarse)
+    subroutine amg_c_invt_solver_descr(sv,info,iout,coarse,prefix)
       import :: psb_spk_, amg_c_invt_solver_type, psb_ipk_
 
       Implicit None
 
       ! Arguments
       class(amg_c_invt_solver_type), intent(in) :: sv
-      integer(psb_ipk_), intent(out)            :: info
-      integer(psb_ipk_), intent(in), optional   :: iout
-      logical, intent(in), optional             :: coarse
+      integer(psb_ipk_), intent(out)              :: info
+      integer(psb_ipk_), intent(in), optional     :: iout
+      logical, intent(in), optional               :: coarse
+      character(len=*), intent(in), optional      :: prefix
 
     end subroutine amg_c_invt_solver_descr
   end interface

amgprec/amg_c_jac_smoother.f90

@@ -219,12 +219,13 @@ module amg_c_jac_smoother
   end interface
 
   interface
-    subroutine amg_c_jac_smoother_descr(sm,info,iout,coarse)
+    subroutine amg_c_jac_smoother_descr(sm,info,iout,coarse,prefix)
       import :: amg_c_jac_smoother_type, psb_ipk_
       class(amg_c_jac_smoother_type), intent(in) :: sm
       integer(psb_ipk_), intent(out)               :: info
       integer(psb_ipk_), intent(in), optional      :: iout
-      logical, intent(in), optional              :: coarse
+      logical, intent(in), optional                :: coarse
+      character(len=*), intent(in), optional       :: prefix
     end subroutine amg_c_jac_smoother_descr
   end interface
 
@@ -313,12 +314,13 @@ module amg_c_jac_smoother
   end interface
 
   interface
-    subroutine amg_c_l1_jac_smoother_descr(sm,info,iout,coarse)
+    subroutine amg_c_l1_jac_smoother_descr(sm,info,iout,coarse,prefix)
       import :: amg_c_l1_jac_smoother_type, psb_ipk_
       class(amg_c_l1_jac_smoother_type), intent(in) :: sm
-      integer(psb_ipk_), intent(out)               :: info
-      integer(psb_ipk_), intent(in), optional      :: iout
-      logical, intent(in), optional              :: coarse
+      integer(psb_ipk_), intent(out)                  :: info
+      integer(psb_ipk_), intent(in), optional         :: iout
+      logical, intent(in), optional                   :: coarse
+      character(len=*), intent(in), optional          :: prefix
     end subroutine amg_c_l1_jac_smoother_descr
   end interface
 

amgprec/amg_c_krm_solver.f90

@@ -436,7 +436,7 @@ contains
     val = "KRM solver"
   end function c_krm_solver_get_fmt
 
-  subroutine c_krm_solver_descr(sv,info,iout,coarse)
+  subroutine c_krm_solver_descr(sv,info,iout,coarse,prefix)
 
     Implicit None
 
@@ -444,12 +444,14 @@ contains
     class(amg_c_krm_solver_type), intent(in) :: sv
     integer(psb_ipk_), intent(out)             :: info
     integer(psb_ipk_), intent(in), optional    :: iout
-    logical, intent(in), optional       :: coarse
+    logical, intent(in), optional              :: coarse
+    character(len=*), intent(in), optional  :: prefix
 
     ! Local variables
     integer(psb_ipk_)      :: err_act
     character(len=20), parameter :: name='amg_c_krm_solver_descr'
     integer(psb_ipk_) :: iout_
+    character(1024)    :: prefix_
 
     call psb_erractionsave(err_act)
     info = psb_success_
@@ -458,23 +460,22 @@ contains
     else
       iout_ = psb_out_unit
     endif
+    if (present(prefix)) then
+      prefix_ = prefix
+    else
+      prefix_ = ""
+    end if
 
     if (sv%global) then
-      write(iout_,*) '  Krylov solver (global)'
+      write(iout_,*) trim(prefix_), '  Krylov solver (global)'
     else
-      write(iout_,*) '  Krylov solver (local) '
+      write(iout_,*) trim(prefix_), '  Krylov solver (local) '
     end if
-    write(iout_,*) '    method: ',sv%method
-    write(iout_,*) '     kprec: ',sv%kprec
-    if (sv%i_sub_solve > 0) then
-      write(iout_,*) ' sub_solve: ',amg_fact_names(sv%i_sub_solve)
-    else
-      write(iout_,*) ' sub_solve: ',sv%sub_solve
-    end if
-    write(iout_,*) '     itmax: ',sv%itmax
-    write(iout_,*) '       eps: ',sv%eps
-    write(iout_,*) '    fillin: ',sv%fillin
-
+    write(iout_,*) trim(prefix_), '    method: ',sv%method
+    write(iout_,*) trim(prefix_), '     kprec: ',sv%kprec
+    call sv%prec%descr(info,iout_,prefix='KRM : '//prefix_)
+    write(iout_,*) trim(prefix_), '     itmax: ',sv%itmax
+    write(iout_,*) trim(prefix_), '       eps: ',sv%eps
 
     call psb_erractionrestore(err_act)
     return

amgprec/amg_c_mumps_solver.F90

@@ -313,22 +313,24 @@ subroutine c_mumps_solver_finalize(sv)
 
 end subroutine c_mumps_solver_finalize
 
-subroutine c_mumps_solver_descr(sv,info,iout,coarse)
+subroutine c_mumps_solver_descr(sv,info,iout,coarse,prefix)
 
   Implicit None
 
   ! Arguments
   class(amg_c_mumps_solver_type), intent(in) :: sv
-  integer(psb_ipk_), intent(out)                     :: info
-  integer(psb_ipk_), intent(in), optional            :: iout
-  logical, intent(in), optional       :: coarse
-
+  integer(psb_ipk_), intent(out)               :: info
+  integer(psb_ipk_), intent(in), optional      :: iout
+  logical, intent(in), optional                :: coarse
+  character(len=*), intent(in), optional       :: prefix
+      
   ! Local variables
   integer(psb_ipk_)   :: err_act
   type(psb_ctxt_type) :: ctxt
   integer(psb_ipk_)   :: me, np
   character(len=20), parameter :: name='amg_z_mumps_solver_descr'
   integer(psb_ipk_) :: iout_
+  character(1024)    :: prefix_
 
   call psb_erractionsave(err_act)
   info = psb_success_
@@ -337,8 +339,13 @@ subroutine c_mumps_solver_descr(sv,info,iout,coarse)
   else
     iout_ = psb_out_unit
   endif
+  if (present(prefix)) then
+    prefix_ = prefix
+  else
+    prefix_ = ""
+  end if
 
-  write(iout_,*) '  MUMPS  Solver. '
+  write(iout_,*) trim(prefix_), '  MUMPS  Solver. '
 
   call psb_erractionrestore(err_act)
   return

amgprec/amg_c_onelev_mod.f90

@@ -257,7 +257,7 @@ module amg_c_onelev_mod
   end interface
 
   interface
-    subroutine amg_c_base_onelev_descr(lv,il,nl,ilmin,info,iout, verbosity)
+    subroutine amg_c_base_onelev_descr(lv,il,nl,ilmin,info,iout, verbosity,prefix)
       import :: psb_cspmat_type, psb_c_vect_type, psb_c_base_vect_type, &
            & psb_clinmap_type, psb_spk_, amg_c_onelev_type, &
            & psb_ipk_, psb_epk_, psb_desc_type
@@ -268,6 +268,7 @@ module amg_c_onelev_mod
       integer(psb_ipk_), intent(out)          :: info
       integer(psb_ipk_), intent(in), optional :: iout
       integer(psb_ipk_), intent(in), optional :: verbosity
+      character(len=*), intent(in), optional  :: prefix
     end subroutine amg_c_base_onelev_descr
   end interface
 

amgprec/amg_c_prec_type.f90

@@ -155,15 +155,16 @@ module amg_c_prec_type
 
 
   interface amg_precdescr
-    subroutine amg_cfile_prec_descr(prec,info,iout,root,verbosity)
+    subroutine amg_cfile_prec_descr(prec,info,iout,root,verbosity,prefix)
       import :: amg_cprec_type, psb_ipk_
       implicit none
       ! Arguments
-      class(amg_cprec_type), intent(in)      :: prec
-      integer(psb_ipk_), intent(out) :: info
+      class(amg_cprec_type), intent(in)     :: prec
+      integer(psb_ipk_), intent(out)          :: info
       integer(psb_ipk_), intent(in), optional :: iout
       integer(psb_ipk_), intent(in), optional :: root
       integer(psb_ipk_), intent(in), optional :: verbosity
+      character(len=*), intent(in), optional  :: prefix
     end subroutine amg_cfile_prec_descr
   end interface
 

amgprec/amg_c_slu_solver.F90

@@ -385,20 +385,22 @@ contains
   
   end subroutine c_slu_solver_finalize
 
-  subroutine c_slu_solver_descr(sv,info,iout,coarse)
+  subroutine c_slu_solver_descr(sv,info,iout,coarse,prefix)
 
     Implicit None
 
     ! Arguments
     class(amg_c_slu_solver_type), intent(in) :: sv
-    integer, intent(out)                     :: info
-    integer, intent(in), optional            :: iout
-    logical, intent(in), optional       :: coarse
+    integer, intent(out)                       :: info
+    integer, intent(in), optional              :: iout
+    logical, intent(in), optional              :: coarse
+    character(len=*), intent(in), optional     :: prefix
 
     ! Local variables
     integer      :: err_act
     character(len=20), parameter :: name='amg_c_slu_solver_descr'
     integer :: iout_
+    character(1024)    :: prefix_
 
     call psb_erractionsave(err_act)
     info = psb_success_
@@ -407,8 +409,13 @@ contains
     else
       iout_ = psb_out_unit
     endif
+    if (present(prefix)) then
+      prefix_ = prefix
+    else
+      prefix_ = ""
+    end if
     
-    write(iout_,*) '  SuperLU Sparse Factorization Solver. '
+    write(iout_,*) trim(prefix_), '  SuperLU Sparse Factorization Solver. '
 
     call psb_erractionrestore(err_act)
     return

amgprec/amg_c_symdec_aggregator_mod.f90

@@ -88,16 +88,25 @@ contains
     val = "Symmetric Decoupled aggregation"
   end function amg_c_symdec_aggregator_fmt
   
-  subroutine  amg_c_symdec_aggregator_descr(ag,parms,iout,info)
+  subroutine  amg_c_symdec_aggregator_descr(ag,parms,iout,info,prefix)
     implicit none 
     class(amg_c_symdec_aggregator_type), intent(in) :: ag
     type(amg_sml_parms), intent(in)   :: parms
     integer(psb_ipk_), intent(in)  :: iout
     integer(psb_ipk_), intent(out) :: info
+    character(len=*), intent(in), optional  :: prefix
 
-    write(iout,*) 'Decoupled Aggregator locally-symmetrized'
-    write(iout,*) 'Aggregator object type: ',ag%fmt()
-    call parms%mldescr(iout,info)
+    character(1024)    :: prefix_
+
+    if (present(prefix)) then
+      prefix_ = prefix
+    else
+      prefix_ = ""
+    end if
+
+    write(iout,*) trim(prefix_),' ','Decoupled Aggregator locally-symmetrized'
+    write(iout,*) trim(prefix_),' ','Aggregator object type: ',ag%fmt()
+    call parms%mldescr(iout,info,prefix=prefix)
     
     return
   end subroutine amg_c_symdec_aggregator_descr

amgprec/amg_d_ainv_solver.F90

@@ -198,7 +198,7 @@ module amg_d_ainv_solver
 !!$  end interface
 
   interface
-    subroutine amg_d_ainv_solver_descr(sv,info,iout,coarse)
+    subroutine amg_d_ainv_solver_descr(sv,info,iout,coarse,prefix)
       import :: psb_dpk_, amg_d_ainv_solver_type, psb_ipk_
 
       Implicit None
@@ -208,7 +208,7 @@ module amg_d_ainv_solver
       integer(psb_ipk_), intent(out)            :: info
       integer(psb_ipk_), intent(in), optional   :: iout
       logical, intent(in), optional             :: coarse
-
+      character(len=*), intent(in), optional  :: prefix
     end subroutine amg_d_ainv_solver_descr
   end interface
 

amgprec/amg_d_as_smoother.f90

@@ -396,21 +396,23 @@ contains
   end subroutine d_as_smoother_default
 
 
-  subroutine d_as_smoother_descr(sm,info,iout,coarse)
+  subroutine d_as_smoother_descr(sm,info,iout,coarse,prefix)
 
     Implicit None
 
     ! Arguments
     class(amg_d_as_smoother_type), intent(in) :: sm
-    integer(psb_ipk_), intent(out)                      :: info
-    integer(psb_ipk_), intent(in), optional             :: iout
-    logical, intent(in), optional             :: coarse
+    integer(psb_ipk_), intent(out)              :: info
+    integer(psb_ipk_), intent(in), optional     :: iout
+    logical, intent(in), optional              :: coarse
+      character(len=*), intent(in), optional   :: prefix
 
     ! Local variables
     integer(psb_ipk_)      :: err_act
     character(len=20), parameter :: name='amg_d_as_smoother_descr'
     integer(psb_ipk_) :: iout_
     logical      :: coarse_
+    character(1024)    :: prefix_
 
     call psb_erractionsave(err_act)
     info = psb_success_
@@ -424,16 +426,21 @@ contains
     else
       iout_ = psb_out_unit
     endif
+    if (present(prefix)) then
+      prefix_ = prefix
+    else
+      prefix_ = ""
+    end if
 
     if (.not.coarse_) then 
-      write(iout_,*) '  Additive Schwarz with  ',&
+      write(iout_,*) trim(prefix_), '  Additive Schwarz with  ',&
            &  sm%novr, ' overlap layers.'
-      write(iout_,*) '  Restrictor:  ',restrict_names(sm%restr)
-      write(iout_,*) '  Prolongator: ',prolong_names(sm%prol)
-      write(iout_,*) '  Local solver:'
+      write(iout_,*) trim(prefix_), '  Restrictor:  ',restrict_names(sm%restr)
+      write(iout_,*) trim(prefix_), '  Prolongator: ',prolong_names(sm%prol)
+      write(iout_,*) trim(prefix_), '  Local solver:'
     endif
     if (allocated(sm%sv)) then 
-      call sm%sv%descr(info,iout_,coarse=coarse)
+      call sm%sv%descr(info,iout_,coarse=coarse,prefix=prefix)
     end if
 
     call psb_erractionrestore(err_act)

amgprec/amg_d_base_aggregator_mod.f90

@@ -275,15 +275,22 @@ contains
     val = .false.
   end function amg_d_base_aggregator_xt_desc
 
-  subroutine  amg_d_base_aggregator_descr(ag,parms,iout,info)
+  subroutine  amg_d_base_aggregator_descr(ag,parms,iout,info,prefix)
     implicit none 
     class(amg_d_base_aggregator_type), intent(in) :: ag
     type(amg_dml_parms), intent(in)   :: parms
     integer(psb_ipk_), intent(in)  :: iout
     integer(psb_ipk_), intent(out) :: info
+    character(len=*), intent(in), optional  :: prefix
+    character(1024)    :: prefix_
+    if (present(prefix)) then
+      prefix_ = prefix
+    else
+      prefix_ = ""
+    end if
 
-    write(iout,*) 'Aggregator object type: ',ag%fmt()
-    call parms%mldescr(iout,info)
+    write(iout,*) trim(prefix_),' ', 'Aggregator object type: ',ag%fmt()
+    call parms%mldescr(iout,info,prefix=prefix)
     
     return
   end subroutine amg_d_base_aggregator_descr

amgprec/amg_d_base_smoother_mod.f90

@@ -272,7 +272,7 @@ module amg_d_base_smoother_mod
   end interface
   
   interface 
-    subroutine amg_d_base_smoother_descr(sm,info,iout,coarse)
+    subroutine amg_d_base_smoother_descr(sm,info,iout,coarse,prefix)
       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_base_smoother_type, psb_ipk_
@@ -281,6 +281,7 @@ module amg_d_base_smoother_mod
       integer(psb_ipk_), intent(out)                :: info
       integer(psb_ipk_), intent(in), optional       :: iout
       logical, intent(in), optional                 :: coarse
+      character(len=*), intent(in), optional        :: prefix
     end subroutine amg_d_base_smoother_descr
   end interface
   

amgprec/amg_d_base_solver_mod.f90

@@ -270,7 +270,7 @@ module amg_d_base_solver_mod
   end interface
   
   interface
-    subroutine amg_d_base_solver_descr(sv,info,iout,coarse)
+    subroutine amg_d_base_solver_descr(sv,info,iout,coarse,prefix)
       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_base_solver_type, psb_ipk_
@@ -281,7 +281,7 @@ module amg_d_base_solver_mod
       integer(psb_ipk_), intent(out)              :: info
       integer(psb_ipk_), intent(in), optional     :: iout
       logical, intent(in), optional               :: coarse
-
+      character(len=*), intent(in), optional      :: prefix
     end subroutine amg_d_base_solver_descr
   end interface 
   

amgprec/amg_d_dec_aggregator_mod.f90

@@ -184,16 +184,23 @@ contains
     val = "Decoupled aggregation"
   end function amg_d_dec_aggregator_fmt
   
-  subroutine  amg_d_dec_aggregator_descr(ag,parms,iout,info)
+  subroutine  amg_d_dec_aggregator_descr(ag,parms,iout,info,prefix)
     implicit none 
     class(amg_d_dec_aggregator_type), intent(in) :: ag
     type(amg_dml_parms), intent(in)   :: parms
     integer(psb_ipk_), intent(in)  :: iout
     integer(psb_ipk_), intent(out) :: info
+    character(len=*), intent(in), optional  :: prefix
+    character(1024)    :: prefix_
+    if (present(prefix)) then
+      prefix_ = prefix
+    else
+      prefix_ = ""
+    end if
 
-    write(iout,*) 'Decoupled Aggregator'
-    write(iout,*) 'Aggregator object type: ',ag%fmt()
-    call parms%mldescr(iout,info)
+    write(iout,*) trim(prefix_),' ','Decoupled Aggregator'
+    write(iout,*) trim(prefix_),' ','Aggregator object type: ',ag%fmt()
+    call parms%mldescr(iout,info,prefix=prefix)
     
     return
   end subroutine amg_d_dec_aggregator_descr

amgprec/amg_d_decmatch_mod.f90

@@ -0,0 +1,529 @@
+!
+!
+!                             AMG4PSBLAS version 1.0
+!    Algebraic Multigrid Package
+!               based on PSBLAS (Parallel Sparse BLAS version 3.7)
+!
+!    (C) Copyright 2021
+!
+!        Salvatore Filippone
+!        Pasqua D'Ambra
+!        Fabio Durastante
+!
+!    Redistribution and use in source and binary forms, with or without
+!    modification, are permitted provided that the following conditions
+!    are met:
+!      1. Redistributions of source code must retain the above copyright
+!         notice, this list of conditions and the following disclaimer.
+!      2. Redistributions in binary form must reproduce the above copyright
+!         notice, this list of conditions, and the following disclaimer in the
+!         documentation and/or other materials provided with the distribution.
+!      3. The name of the AMG4PSBLAS group or the names of its contributors may
+!         not be used to endorse or promote products derived from this
+!         software without specific written permission.
+!
+!    THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+!    ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
+!    TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
+!    PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE AMG4PSBLAS GROUP OR ITS CONTRIBUTORS
+!    BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
+!    CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
+!    SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
+!    INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
+!    CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
+!    ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+!    POSSIBILITY OF SUCH DAMAGE.
+!
+module amg_d_decmatch_mod
+
+  use iso_c_binding
+  use psb_base_cbind_mod
+
+  interface new_Match_If
+    function dnew_Match_If(ipar,matching,lambda,nr, irp, ja, val, diag, w, mate) &
+         & bind(c,name="dnew_Match_If") result(res)
+      use iso_c_binding
+      import :: psb_c_ipk_, psb_c_lpk_, psb_c_mpk_, psb_c_epk_
+      implicit none
+
+      integer(psb_c_ipk_) :: res
+      integer(psb_c_ipk_), value :: nr,ipar,matching
+      real(c_double), value :: lambda      
+      type(c_ptr), value :: irp, ja, mate
+      type(c_ptr), value :: val, diag, w 
+    end function dnew_Match_If
+  end interface new_Match_If
+
+  interface amg_build_decmatch
+    module procedure amg_dbuild_decmatch
+  end interface amg_build_decmatch
+
+  logical, parameter, private :: print_statistics=.false.
+contains
+
+  subroutine amg_ddecmatch_build_prol(w,a,desc_a,ilaggr,nlaggr,prol,info,&
+       & symmetrize,reproducible,display_inp, display_out, print_out, &
+       & parallel, matching,lambda)
+    use psb_base_mod
+    use psb_util_mod
+    use iso_c_binding
+    implicit none
+    real(psb_dpk_), allocatable, intent(inout)  :: w(:)
+    type(psb_dspmat_type), intent(inout)        :: a
+    type(psb_desc_type)                         :: desc_a
+    integer(psb_lpk_), allocatable, intent(out) :: ilaggr(:)
+    integer(psb_lpk_), allocatable, intent(out) :: nlaggr(:)
+    type(psb_ldspmat_type), intent(out)         :: prol
+    integer(psb_ipk_), intent(out)              :: info
+    logical, optional, intent(in)               :: display_inp, display_out, reproducible
+    logical, optional, intent(in)               :: symmetrize, print_out, parallel
+    integer(psb_ipk_), optional, intent(in)     :: matching
+    real(psb_dpk_), optional, intent(in)        :: lambda
+
+    !
+    !
+    type(psb_ctxt_type) :: ictxt
+    integer(psb_ipk_) :: iam, np, iown
+    integer(psb_ipk_) :: nr, nc, sweep, nzl, ncsave, nct, idx
+    integer(psb_lpk_) :: i, k, kg, idxg, ntaggr, naggrm1, naggrp1, &
+         & ip, nlpairs, nlsingl, nunmatched, lnr
+    real(psb_dpk_)    :: wk, widx, wmax, nrmagg
+    real(psb_dpk_), allocatable    :: wtemp(:)
+    integer(psb_ipk_), allocatable :: mate(:)
+    integer(psb_lpk_), allocatable ::  ilv(:)
+    integer(psb_ipk_), save     :: cnt=1
+    character(len=256)          :: aname
+    type(psb_ld_coo_sparse_mat) :: tmpcoo
+    logical  :: display_out_, print_out_, reproducible_, parallel_
+    integer(psb_ipk_) :: matching_
+    real(psb_dpk_)    :: lambda_    
+    logical, parameter  :: dump=.false., debug=.false., dump_mate=.false., &
+         & debug_ilaggr=.false., debug_sync=.false.
+    integer(psb_ipk_), save :: idx_bldmtc=-1, idx_phase1=-1, idx_phase2=-1, idx_phase3=-1
+    logical, parameter :: do_timings=.true.
+
+    ictxt = desc_a%get_ctxt()
+    call psb_info(ictxt,iam,np)
+
+    if ((do_timings).and.(idx_phase1==-1))       &
+         & idx_phase1 = psb_get_timer_idx("MBP_BLDP: phase1       ")
+    if ((do_timings).and.(idx_bldmtc==-1))       &
+         & idx_bldmtc = psb_get_timer_idx("MBP_BLDP: buil_matching")
+    if ((do_timings).and.(idx_phase2==-1))       &
+         & idx_phase2 = psb_get_timer_idx("MBP_BLDP: phase2       ")
+    if ((do_timings).and.(idx_phase3==-1))       &
+         & idx_phase3 = psb_get_timer_idx("MBP_BLDP: phase3       ")
+
+    if (do_timings) call psb_tic(idx_phase1)
+
+    if (present(display_out)) then
+      display_out_ = display_out
+    else
+      display_out_ = .false.
+    end if
+    if (present(print_out)) then
+      print_out_ = print_out
+    else
+      print_out_ = .false.
+    end if
+    if (present(reproducible)) then
+      reproducible_ = reproducible
+    else
+      reproducible_ = .false.
+    end if
+
+    if (present(parallel)) then
+      parallel_ = parallel
+    else
+      parallel_ = .true.
+    end if
+
+    if (present(matching)) then
+      matching_ = matching
+    else
+      matching_ = 2
+    end if
+
+    if (present(lambda)) then
+      lambda_ = lambda
+    else
+      lambda_ = 2.0
+    end if
+
+    allocate(nlaggr(0:np-1),stat=info)
+    if (info /= 0) then
+      return
+    end if
+
+    nlaggr = 0
+    ilv    = [(i,i=1,desc_a%get_local_cols())]
+    call desc_a%l2gip(ilv,info,owned=.false.)
+
+    call psb_geall(ilaggr,desc_a,info)
+    ilaggr = -1
+    call psb_geasb(ilaggr,desc_a,info)
+    nr = a%get_nrows()
+    nc = a%get_ncols()
+    if (size(w) < nc) then
+      call psb_realloc(nc,w,info)
+    end if
+    call psb_halo(w,desc_a,info)
+
+    if (debug) write(0,*) iam,' buildprol into buildmatching:',&
+         & nr, nc
+    if (debug_sync) then
+      call psb_barrier(ictxt)
+      if (iam == 0) write(0,*)' buildprol into buildmatching:',&
+           & nr, nc
+    end if
+    if (do_timings) call psb_toc(idx_phase1)
+    if (do_timings) call psb_tic(idx_bldmtc)
+    call amg_dbuild_decmatch(parallel_,matching_,lambda_,w,a,desc_a,mate,info)
+    if (do_timings) call psb_toc(idx_bldmtc)
+    if (debug) write(0,*) iam,' buildprol from buildmatching:',&
+         & info
+    if (debug_sync) then
+      call psb_barrier(ictxt)
+      if (iam == 0) write(0,*)' out from buildmatching:', info
+    end if
+
+    if (info == 0) then
+      if (do_timings) call psb_tic(idx_phase2)
+      if (debug_sync) then
+        call psb_barrier(ictxt)
+        if (iam == 0) write(0,*)' Into building the tentative prol:'
+      end if
+
+      call psb_geall(wtemp,desc_a,info)
+      wtemp      = dzero
+      call psb_geasb(wtemp,desc_a,info)
+
+      nlaggr(iam) = 0
+      nlpairs    = 0
+      nlsingl    = 0
+      nunmatched = 0
+      !
+      ! First sweep
+      ! On return from build_matching, mate has been converted to local numbering,
+      ! so assigning to idx is OK.
+      !
+      do k=1, nr
+        idx = mate(k)
+        !
+        ! Figure out an allocation of aggregates to processes
+        !
+        if (idx < 0) then
+          !
+          !  Unmatched vertex, potential  singleton.
+          !
+          nunmatched = nunmatched + 1
+          if (abs(w(k))<epsilon(nrmagg)) then
+            ! Keep it unaggregated
+            wtemp(k) = dzero
+          else
+            ! Create a singleton aggregate
+            nlaggr(iam) = nlaggr(iam) + 1
+            ilaggr(k)   = nlaggr(iam)
+            wtemp(k)    = w(k)/abs(w(k))
+            nlsingl     = nlsingl + 1
+          end if
+!!$          write(0,*) k,mate(k),ilaggr(k),' negative match ',abs(w(k)), epsilon(nrmagg)
+        else  if (idx > nc) then
+          write(0,*) 'Impossible: mate(k) > nc'
+          cycle
+        else
+
+          if (ilaggr(k) == -1) then
+
+            wk   = w(k)
+            widx = w(idx)
+            wmax = max(abs(wk),abs(widx))
+            nrmagg = wmax*sqrt((wk/wmax)**2+(widx/wmax)**2)
+            if (nrmagg > epsilon(nrmagg)) then
+              if (idx <= nr) then
+                if (ilaggr(idx) == -1) then
+                  ! Now, if both vertices are local, the aggregate is local
+                  ! (kinda obvious).
+                  nlaggr(iam) = nlaggr(iam) + 1
+                  ilaggr(k)   = nlaggr(iam)
+                  ilaggr(idx) = nlaggr(iam)
+                  wtemp(k)    = w(k)/nrmagg
+                  wtemp(idx)  = w(idx)/nrmagg
+                end if
+                nlpairs = nlpairs+1
+              else
+                write(0,*) 'Really?  mate(k) > nr? ',mate(k),nr
+              end if
+            else
+              if (abs(w(k))<epsilon(nrmagg)) then
+                ! Keep it unaggregated
+                wtemp(k) = dzero
+              else
+                ! Create a singleton aggregate
+                nlaggr(iam) = nlaggr(iam) + 1
+                ilaggr(k)   = nlaggr(iam)
+                wtemp(k)    = w(k)/abs(w(k))
+                nlsingl     = nlsingl + 1
+              end if
+            end if
+          end if
+        end if
+      end do
+      if (do_timings) call psb_toc(idx_phase2)
+      if (do_timings) call psb_tic(idx_phase3)
+
+      ! Ok, now compute offsets, gather halo and fix non-local
+      ! aggregates  (those where ilaggr == -2)
+      call psb_sum(ictxt,nlaggr)
+      ntaggr  = sum(nlaggr(0:np-1))
+      naggrm1 = sum(nlaggr(0:iam-1))
+      naggrp1 = sum(nlaggr(0:iam))
+      !
+      ! Shift all indices already assigned (i.e. >0)
+      !
+      do k=1,nr
+        if (ilaggr(k) > 0) then
+          ilaggr(k) = ilaggr(k) + naggrm1
+!!$        else
+!!$          write(0,*) 'Leftover ILAGGR',k,ilaggr(k),mate(k),abs(w(k)),epsilon(nrmagg)
+        end if
+      end do
+      call psb_halo(ilaggr,desc_a,info)
+      call psb_halo(wtemp,desc_a,info)
+      ! Cleanup as yet unmarked entries
+      do k=1,nr
+        if (ilaggr(k) == -2) then
+          idx = mate(k)
+          if (idx > nr) then
+            i   = ilaggr(idx)
+            if (i > 0) then
+              ilaggr(k) = i
+            else
+              write(0,*) 'Error : unresolved (paired) index ',k,idx,i,nr,nc, ilv(k),ilv(idx)
+            end if
+          else
+            write(0,*) 'Error : unresolved (paired) index ',k,idx,i,nr,nc, ilv(k),ilv(idx)
+          end if
+        end if
+        if (ilaggr(k) <0) then
+          write(0,*) 'Decmatch: Funny number: ',k,ilv(k),ilaggr(k),wtemp(k)
+        end if
+      end do
+      if (debug_sync) then
+        call psb_barrier(ictxt)
+        if (iam == 0) write(0,*)' Done building the tentative prol:'
+      end if
+
+
+      if (dump_mate) then
+        block
+          integer(psb_lpk_), allocatable :: glaggr(:)
+          write(aname,'(a,i3.3,a,i3.3,a)') 'mateg-',cnt,'-p',iam,'.mtx'
+          open(20,file=aname)
+          write(20,'(a,I3,a)') '% sparse vector on process ',iam,' '
+          do k=1, nr
+            write(20,'(3(I8,1X))') ilv(k),ilv(mate(k))
+          end do
+          close(20)
+          write(aname,'(a,i3.3,a,i3.3,a)') 'nloc-',cnt,'-p',iam,'.mtx'
+          open(20,file=aname)
+          write(20,'(a,I3,a)') '% sparse vector on process ',iam,' '
+          write(20,'(a,I12,a)') 'nlpairs ',nlpairs
+          write(20,'(a,I12,a)') 'nlsingl ',nlsingl
+          write(20,'(a,I12,a)') 'nlaggr(iam) ',nlaggr(iam)
+          close(20)
+          write(aname,'(a,i3.3,a,i3.3,a,i3.3,a)') 'ilaggr-',cnt,'-i',iam,'-p',np,'.mtx'
+          open(20,file=aname)
+          write(20,'(a,I3,a)') '% sparse vector on process ',iam,' '
+          do k=1, nr
+            write(20,'(3(I8,1X))') ilv(k),ilaggr(k)
+          end do
+          close(20)
+
+          write(aname,'(a,i3.3,a,i3.3,a)') 'glaggr-',cnt,'-p',np,'.mtx'
+          call psb_gather(glaggr,ilaggr,desc_a,info,root=izero)
+          if (iam==0) call mm_array_write(glaggr,'Aggregates ',info,filename=aname)
+
+          cnt=cnt+1
+        end block
+      end if
+      block
+        integer(psb_lpk_) :: v(3)
+        v(1) = nunmatched
+        v(2) = nlsingl
+        v(3) = nlpairs
+        call psb_sum(ictxt,v)
+        nunmatched = v(1)
+        nlsingl    = v(2)
+        nlpairs    = v(3)
+
+      end block
+      if (print_statistics) then
+        if (iam == 0) then
+          write(0,*) 'Matching statistics: Unmatched nodes ',&
+               & nunmatched,' Singletons:',nlsingl,' Pairs:',nlpairs
+        end if
+      end if
+      
+      if (display_out_) then
+        block
+          integer(psb_ipk_) :: idx
+          !
+          !  And finally print out
+          !
+          do i=0,np-1
+            call psb_barrier(ictxt)
+            if (iam == i) then
+              write(0,*) 'Process ', iam,' hosts aggregates:  (',naggrm1+1,' : ',naggrp1-1,')'
+              do k=1, nr
+                idx  = mate(k)
+                kg    = ilv(k)
+                if (idx >0) then
+                  idxg  = ilv(idx)
+                else
+                  idxg = -1
+                end if
+                if (idx < 0) then
+                  write(0,*) kg,': singleton   (',kg,' ( Proc',iam,') ) into aggregate => ', ilaggr(k)
+                else if (idx <= nr) then
+                  write(0,*) kg,': paired with (',idxg,' ( Proc',iam,') ) into aggregate => ', ilaggr(k)
+                else
+                  call desc_a%indxmap%qry_halo_owner(idx,iown,info)
+                  write(0,*) kg,': paired with (',idxg,' ( Proc',iown,') ) into aggregate => ', ilaggr(k)
+                end if
+              end do
+              flush(0)
+            end if
+          end do
+        end block
+      end if
+
+      ! Dirty trick: allocate tmpcoo with local
+      ! number of aggregates, then change to ntaggr.
+      ! Just to make sure the allocation is not global
+      lnr = nr
+      call tmpcoo%allocate(lnr,nlaggr(iam),lnr)
+      k = 0
+      do i=1,nr
+        !
+        ! Note: at this point, a value ilaggr(i)<=0
+        ! tags an unaggregated row, and it has to be
+        ! left alone (i.e.: it should stay at fine level only)
+        !
+        if (ilaggr(i)>0) then
+          k = k + 1
+          tmpcoo%val(k) = wtemp(i)
+          tmpcoo%ia(k)  = i
+          tmpcoo%ja(k)  = ilaggr(i)
+        end if
+      end do
+      call tmpcoo%set_nzeros(k)
+      call tmpcoo%set_dupl(psb_dupl_add_)
+      call tmpcoo%set_sorted() ! This is now in row-major
+
+      if (display_out_) then
+        call psb_barrier(ictxt)
+        flush(0)
+
+        if (iam == 0) write(0,*) 'Prolongator: '
+        flush(0)
+        do i=0,np-1
+          call psb_barrier(ictxt)
+          if (iam == i) then
+            do k=1, nr
+              write(0,*) ilv(tmpcoo%ia(k)),tmpcoo%ja(k), tmpcoo%val(k)
+            end do
+            flush(0)
+          end if
+        end do
+      end if
+
+      call prol%mv_from(tmpcoo)
+      if (do_timings) call psb_toc(idx_phase3)
+
+      if (print_out_) then
+        write(aname,'(a,i3.3,a)') 'prol-g-',iam,'.mtx'
+        call prol%print(fname=aname,head='Test ',ivr=ilv)
+        write(aname,'(a,i3.3,a)') 'prol-',iam,'.mtx'
+        call prol%print(fname=aname,head='Test ')
+      end if
+
+    else
+      write(0,*) iam,' :  error from Matching: ',info
+    end if
+
+  end subroutine amg_ddecmatch_build_prol
+
+  subroutine amg_dbuild_decmatch(parallel,matching,lambda,w,a,desc_a,mate,info)
+    use psb_base_mod
+    use psb_util_mod
+    use iso_c_binding
+    implicit none
+    logical, intent(in)               :: parallel
+    integer(psb_ipk_), intent(in)     :: matching
+    real(psb_dpk_), intent(in)        :: lambda
+    real(psb_dpk_), target                      :: w(:)
+    type(psb_dspmat_type), intent(inout)        :: a
+    type(psb_desc_type)                         :: desc_a
+    integer(psb_ipk_), allocatable, intent(out), target :: mate(:)
+    integer(psb_ipk_), intent(out)              :: info
+
+    type(psb_d_csr_sparse_mat), target :: tcsr
+    real(psb_dpk_), allocatable, target :: diag(:)
+    real(psb_dpk_) :: ph0t,ph1t,ph2t
+    type(psb_ctxt_type) :: ictxt
+    integer(psb_ipk_) :: iam, np
+    integer(psb_ipk_)  :: nr, nc, nz, i, nunmatch, ipar 
+    integer(psb_ipk_), save :: cnt=2
+    logical, parameter :: debug=.false., dump_ahat=.false., debug_sync=.false.
+    logical, parameter :: old_style=.false., sort_minp=.true.
+    character(len=40) :: name='build_matching', fname
+    integer(psb_ipk_), save :: idx_cmboxp=-1, idx_bldahat=-1, idx_phase2=-1, idx_phase3=-1
+    logical, parameter :: do_timings=.true.
+
+    ictxt = desc_a%get_ctxt()
+    call psb_info(ictxt,iam,np)
+
+    nr = a%get_nrows()
+    call a%cp_to(tcsr)
+    call psb_realloc(nr,mate,info)
+    diag = a%get_diag(info)
+    if (parallel) then
+      ipar = 2
+    else
+      ipar = 1
+    end if
+    !
+    ! Now call matching!
+    !
+    if (debug) write(0,*) iam,' buildmatching into NewMatch:'
+    if (do_timings) call psb_tic(idx_cmboxp)
+    info = dnew_Match_If(ipar,matching,lambda,nr,c_loc(tcsr%irp),c_loc(tcsr%ja),&
+         & c_loc(tcsr%val),c_loc(diag),c_loc(w),c_loc(mate))
+    if (do_timings) call psb_toc(idx_cmboxp)
+    if (debug) write(0,*) iam,' buildmatching from NewMatch:', info
+    if (debug_sync) then
+      call psb_max(ictxt,info)
+      if (iam == 0) write(0,*)' done NewMatch', info
+    end if
+
+    if (do_timings) call psb_tic(idx_phase3)
+    nunmatch = count(mate(1:nr)<=0)
+        ! call psb_sum(ictxt,nunmatch)
+    !if (nunmatch /= 0) write(0,*) iam,' Unmatched nodes local imbalance ',nunmatch
+    !    if (count(mate(1:nr)<0) /= nunmatch) write(0,*) iam,' Matching results ?',&
+    !        & nunmatch, count(mate(1:nr)<0)
+    if (debug_sync) then
+      call psb_barrier(ictxt)
+      if (iam == 0) write(0,*)' done build_matching '
+    end if
+
+    if (do_timings) call psb_toc(idx_phase3)
+    return
+
+9999 continue
+    call psb_error(ictxt)
+
+  end subroutine amg_dbuild_decmatch
+
+end module amg_d_decmatch_mod

amgprec/amg_d_diag_solver.f90

@@ -219,7 +219,7 @@ contains
 
   end subroutine d_diag_solver_free
 
-  subroutine d_diag_solver_descr(sv,info,iout,coarse)
+  subroutine d_diag_solver_descr(sv,info,iout,coarse,prefix)
 
     Implicit None
 
@@ -228,11 +228,13 @@ contains
     integer(psb_ipk_), intent(out)              :: info
     integer(psb_ipk_), intent(in), optional     :: iout
     logical, intent(in), optional               :: coarse
+    character(len=*), intent(in), optional      :: prefix
 
     ! Local variables
     integer(psb_ipk_)      :: err_act
     character(len=20), parameter :: name='amg_d_diag_solver_descr'
     integer(psb_ipk_) :: iout_
+    character(1024)    :: prefix_
 
     info = psb_success_
     if (present(iout)) then 
@@ -240,8 +242,13 @@ contains
     else
       iout_ = psb_out_unit
     endif
-    
-    write(iout_,*) '  Diagonal local solver '
+    if (present(prefix)) then
+      prefix_ = prefix
+    else
+      prefix_ = ""
+    end if
+
+    write(iout_,*) trim(prefix_), '  Diagonal local solver '
 
     return
 
@@ -352,7 +359,7 @@ module amg_d_l1_diag_solver
   
 contains
 
-  subroutine d_l1_diag_solver_descr(sv,info,iout,coarse)
+  subroutine d_l1_diag_solver_descr(sv,info,iout,coarse,prefix)
 
     Implicit None
 
@@ -361,11 +368,13 @@ contains
     integer(psb_ipk_), intent(out)              :: info
     integer(psb_ipk_), intent(in), optional     :: iout
     logical, intent(in), optional               :: coarse
+    character(len=*), intent(in), optional  :: prefix
 
     ! Local variables
     integer(psb_ipk_)      :: err_act
     character(len=20), parameter :: name='amg_d_l1_diag_solver_descr'
     integer(psb_ipk_) :: iout_
+    character(1024)    :: prefix_
 
     info = psb_success_
     if (present(iout)) then 
@@ -373,8 +382,13 @@ contains
     else
       iout_ = psb_out_unit
     endif
-    
-    write(iout_,*) '  L1 Diagonal solver '
+    if (present(prefix)) then
+      prefix_ = prefix
+    else
+      prefix_ = ""
+    end if
+
+    write(iout_,*) trim(prefix_), '  L1 Diagonal solver '
 
     return
 

amgprec/amg_d_gs_solver.f90

@@ -433,20 +433,22 @@ contains
     return
   end subroutine d_gs_solver_free
 
-  subroutine d_gs_solver_descr(sv,info,iout,coarse)
+  subroutine d_gs_solver_descr(sv,info,iout,coarse,prefix)
 
     Implicit None
 
     ! Arguments
     class(amg_d_gs_solver_type), intent(in) :: sv
-    integer(psb_ipk_), intent(out)             :: info
-    integer(psb_ipk_), intent(in), optional    :: iout
-    logical, intent(in), optional       :: coarse
+    integer(psb_ipk_), intent(out)            :: info
+    integer(psb_ipk_), intent(in), optional   :: iout
+    logical, intent(in), optional             :: coarse
+    character(len=*), intent(in), optional    :: prefix
 
     ! Local variables
     integer(psb_ipk_)      :: err_act
     character(len=20), parameter :: name='amg_d_gs_solver_descr'
     integer(psb_ipk_) :: iout_
+    character(1024)    :: prefix_
 
     call psb_erractionsave(err_act)
     info = psb_success_
@@ -455,12 +457,17 @@ contains
     else
       iout_ = psb_out_unit
     endif
+    if (present(prefix)) then
+      prefix_ = prefix
+    else
+      prefix_ = ""
+    end if
 
     if (sv%eps<=dzero) then 
-      write(iout_,*) '  Forward Gauss-Seidel iterative solver with  ',&
+      write(iout_,*) trim(prefix_), '  Forward Gauss-Seidel iterative solver with  ',&
            &  sv%sweeps,' sweeps'
     else
-      write(iout_,*) '  Forward Gauss-Seidel iterative solver with  tolerance',&
+      write(iout_,*) trim(prefix_), '  Forward Gauss-Seidel iterative solver with  tolerance',&
            &  sv%eps,' and maxit', sv%sweeps
     end if
     
@@ -526,20 +533,22 @@ contains
     val = .true.
   end function d_gs_solver_is_iterative
     
-  subroutine d_bwgs_solver_descr(sv,info,iout,coarse)
+  subroutine d_bwgs_solver_descr(sv,info,iout,coarse,prefix)
 
     Implicit None
 
     ! Arguments
     class(amg_d_bwgs_solver_type), intent(in) :: sv
-    integer(psb_ipk_), intent(out)             :: info
-    integer(psb_ipk_), intent(in), optional    :: iout
-    logical, intent(in), optional       :: coarse
+    integer(psb_ipk_), intent(out)              :: info
+    integer(psb_ipk_), intent(in), optional     :: iout
+    logical, intent(in), optional               :: coarse
+    character(len=*), intent(in), optional      :: prefix
 
     ! Local variables
     integer(psb_ipk_)      :: err_act
     character(len=20), parameter :: name='amg_d_bwgs_solver_descr'
     integer(psb_ipk_) :: iout_
+    character(1024)    :: prefix_
 
     call psb_erractionsave(err_act)
     info = psb_success_
@@ -548,12 +557,17 @@ contains
     else
       iout_ = psb_out_unit
     endif
+    if (present(prefix)) then
+      prefix_ = prefix
+    else
+      prefix_ = ""
+    end if
 
     if (sv%eps<=dzero) then 
-      write(iout_,*) '  Backward Gauss-Seidel iterative solver with  ',&
+      write(iout_,*) trim(prefix_), '  Backward Gauss-Seidel iterative solver with  ',&
            &  sv%sweeps,' sweeps'
     else
-      write(iout_,*) '  Backward Gauss-Seidel iterative solver with  tolerance',&
+      write(iout_,*) trim(prefix_), '  Backward Gauss-Seidel iterative solver with  tolerance',&
            &  sv%eps,' and maxit', sv%sweeps
     end if
     

amgprec/amg_d_id_solver.f90

@@ -157,7 +157,7 @@ contains
     return
   end subroutine d_id_solver_free
 
-  subroutine d_id_solver_descr(sv,info,iout,coarse)
+  subroutine d_id_solver_descr(sv,info,iout,coarse,prefix)
 
     Implicit None
 
@@ -165,12 +165,14 @@ contains
     class(amg_d_id_solver_type), intent(in) :: sv
     integer(psb_ipk_), intent(out)            :: info
     integer(psb_ipk_), intent(in), optional   :: iout
-    logical, intent(in), optional       :: coarse
+    logical, intent(in), optional             :: coarse
+    character(len=*), intent(in), optional    :: prefix
 
     ! Local variables
     integer(psb_ipk_)      :: err_act
     character(len=20), parameter :: name='amg_d_id_solver_descr'
     integer(psb_ipk_) :: iout_
+    character(1024)    :: prefix_
 
     info = psb_success_
     if (present(iout)) then 
@@ -178,8 +180,13 @@ contains
     else
       iout_ = psb_out_unit
     endif
+    if (present(prefix)) then
+      prefix_ = prefix
+    else
+      prefix_ = ""
+    end if
     
-    write(iout_,*) '  Identity local solver '
+    write(iout_,*) trim(prefix_), '  Identity local solver '
 
     return
 

amgprec/amg_d_ilu_solver.f90

@@ -406,7 +406,7 @@ contains
     return
   end subroutine d_ilu_solver_free
 
-  subroutine d_ilu_solver_descr(sv,info,iout,coarse)
+  subroutine d_ilu_solver_descr(sv,info,iout,coarse,prefix)
 
     Implicit None
 
@@ -414,12 +414,14 @@ contains
     class(amg_d_ilu_solver_type), intent(in) :: sv
     integer(psb_ipk_), intent(out)             :: info
     integer(psb_ipk_), intent(in), optional    :: iout
-    logical, intent(in), optional       :: coarse
-
+    logical, intent(in), optional              :: coarse
+    character(len=*), intent(in), optional  :: prefix
+      
     ! Local variables
     integer(psb_ipk_) :: err_act
     character(len=20), parameter :: name='amg_d_ilu_solver_descr'
     integer(psb_ipk_) :: iout_
+    character(1024)    :: prefix_
 
     call psb_erractionsave(err_act)
     info = psb_success_
@@ -428,15 +430,20 @@ contains
     else
       iout_ = psb_out_unit
     endif
+    if (present(prefix)) then
+      prefix_ = prefix
+    else
+      prefix_ = ""
+    end if
 
-    write(iout_,*) '  Incomplete factorization solver: ',&
+    write(iout_,*) trim(prefix_), '  Incomplete factorization solver: ',&
          &  amg_fact_names(sv%fact_type)
     select case(sv%fact_type)
     case(psb_ilu_n_,psb_milu_n_)      
-      write(iout_,*) '  Fill level:',sv%fill_in
+      write(iout_,*) trim(prefix_), '  Fill level:',sv%fill_in
     case(psb_ilu_t_)         
-      write(iout_,*) '  Fill level:',sv%fill_in
-      write(iout_,*) '  Fill threshold :',sv%thresh
+      write(iout_,*) trim(prefix_), '  Fill level:',sv%fill_in
+      write(iout_,*) trim(prefix_), '  Fill threshold :',sv%thresh
     end select
 
     call psb_erractionrestore(err_act)

amgprec/amg_d_invk_solver.f90

@@ -123,7 +123,7 @@ module amg_d_invk_solver
   end interface
 
   interface
-    subroutine amg_d_invk_solver_descr(sv,info,iout,coarse)
+    subroutine amg_d_invk_solver_descr(sv,info,iout,coarse,prefix)
       import :: psb_dpk_, amg_d_invk_solver_type, psb_ipk_
 
       Implicit None
@@ -133,7 +133,7 @@ module amg_d_invk_solver
       integer(psb_ipk_), intent(out)            :: info
       integer(psb_ipk_), intent(in), optional   :: iout
       logical, intent(in), optional             :: coarse
-
+      character(len=*), intent(in), optional  :: prefix
     end subroutine amg_d_invk_solver_descr
   end interface
 

amgprec/amg_d_invt_solver.f90

@@ -134,16 +134,17 @@ module amg_d_invt_solver
   end interface
 
   interface
-    subroutine amg_d_invt_solver_descr(sv,info,iout,coarse)
+    subroutine amg_d_invt_solver_descr(sv,info,iout,coarse,prefix)
       import :: psb_dpk_, amg_d_invt_solver_type, psb_ipk_
 
       Implicit None
 
       ! Arguments
       class(amg_d_invt_solver_type), intent(in) :: sv
-      integer(psb_ipk_), intent(out)            :: info
-      integer(psb_ipk_), intent(in), optional   :: iout
-      logical, intent(in), optional             :: coarse
+      integer(psb_ipk_), intent(out)              :: info
+      integer(psb_ipk_), intent(in), optional     :: iout
+      logical, intent(in), optional               :: coarse
+      character(len=*), intent(in), optional      :: prefix
 
     end subroutine amg_d_invt_solver_descr
   end interface

amgprec/amg_d_jac_smoother.f90

@@ -219,12 +219,13 @@ module amg_d_jac_smoother
   end interface
 
   interface
-    subroutine amg_d_jac_smoother_descr(sm,info,iout,coarse)
+    subroutine amg_d_jac_smoother_descr(sm,info,iout,coarse,prefix)
       import :: amg_d_jac_smoother_type, psb_ipk_
       class(amg_d_jac_smoother_type), intent(in) :: sm
       integer(psb_ipk_), intent(out)               :: info
       integer(psb_ipk_), intent(in), optional      :: iout
-      logical, intent(in), optional              :: coarse
+      logical, intent(in), optional                :: coarse
+      character(len=*), intent(in), optional       :: prefix
     end subroutine amg_d_jac_smoother_descr
   end interface
 
@@ -313,12 +314,13 @@ module amg_d_jac_smoother
   end interface
 
   interface
-    subroutine amg_d_l1_jac_smoother_descr(sm,info,iout,coarse)
+    subroutine amg_d_l1_jac_smoother_descr(sm,info,iout,coarse,prefix)
       import :: amg_d_l1_jac_smoother_type, psb_ipk_
       class(amg_d_l1_jac_smoother_type), intent(in) :: sm
-      integer(psb_ipk_), intent(out)               :: info
-      integer(psb_ipk_), intent(in), optional      :: iout
-      logical, intent(in), optional              :: coarse
+      integer(psb_ipk_), intent(out)                  :: info
+      integer(psb_ipk_), intent(in), optional         :: iout
+      logical, intent(in), optional                   :: coarse
+      character(len=*), intent(in), optional          :: prefix
     end subroutine amg_d_l1_jac_smoother_descr
   end interface
 

amgprec/amg_d_krm_solver.f90

@@ -436,7 +436,7 @@ contains
     val = "KRM solver"
   end function d_krm_solver_get_fmt
 
-  subroutine d_krm_solver_descr(sv,info,iout,coarse)
+  subroutine d_krm_solver_descr(sv,info,iout,coarse,prefix)
 
     Implicit None
 
@@ -444,12 +444,14 @@ contains
     class(amg_d_krm_solver_type), intent(in) :: sv
     integer(psb_ipk_), intent(out)             :: info
     integer(psb_ipk_), intent(in), optional    :: iout
-    logical, intent(in), optional       :: coarse
+    logical, intent(in), optional              :: coarse
+    character(len=*), intent(in), optional  :: prefix
 
     ! Local variables
     integer(psb_ipk_)      :: err_act
     character(len=20), parameter :: name='amg_d_krm_solver_descr'
     integer(psb_ipk_) :: iout_
+    character(1024)    :: prefix_
 
     call psb_erractionsave(err_act)
     info = psb_success_
@@ -458,23 +460,22 @@ contains
     else
       iout_ = psb_out_unit
     endif
+    if (present(prefix)) then
+      prefix_ = prefix
+    else
+      prefix_ = ""
+    end if
 
     if (sv%global) then
-      write(iout_,*) '  Krylov solver (global)'
+      write(iout_,*) trim(prefix_), '  Krylov solver (global)'
     else
-      write(iout_,*) '  Krylov solver (local) '
+      write(iout_,*) trim(prefix_), '  Krylov solver (local) '
     end if
-    write(iout_,*) '    method: ',sv%method
-    write(iout_,*) '     kprec: ',sv%kprec
-    if (sv%i_sub_solve > 0) then
-      write(iout_,*) ' sub_solve: ',amg_fact_names(sv%i_sub_solve)
-    else
-      write(iout_,*) ' sub_solve: ',sv%sub_solve
-    end if
-    write(iout_,*) '     itmax: ',sv%itmax
-    write(iout_,*) '       eps: ',sv%eps
-    write(iout_,*) '    fillin: ',sv%fillin
-
+    write(iout_,*) trim(prefix_), '    method: ',sv%method
+    write(iout_,*) trim(prefix_), '     kprec: ',sv%kprec
+    call sv%prec%descr(info,iout_,prefix='KRM : '//prefix_)
+    write(iout_,*) trim(prefix_), '     itmax: ',sv%itmax
+    write(iout_,*) trim(prefix_), '       eps: ',sv%eps
 
     call psb_erractionrestore(err_act)
     return

amgprec/amg_d_matchboxp_mod.f90

@@ -143,9 +143,10 @@ contains
     type(psb_ld_coo_sparse_mat) :: tmpcoo
     logical  :: display_out_, print_out_, reproducible_
     logical, parameter  :: dump=.false., debug=.false., dump_mate=.false., &
-         & debug_ilaggr=.false., debug_sync=.false.
+         & debug_ilaggr=.false., debug_sync=.false., debug_mate=.false.
     integer(psb_ipk_), save :: idx_bldmtc=-1, idx_phase1=-1, idx_phase2=-1, idx_phase3=-1
     logical, parameter :: do_timings=.true.
+    integer, parameter :: ilaggr_neginit=-1, ilaggr_nonlocal=-2
 
     ictxt = desc_a%get_ctxt()
     call psb_info(ictxt,iam,np)
@@ -187,7 +188,7 @@ contains
     call desc_a%l2gip(ilv,info,owned=.false.)
 
     call psb_geall(ilaggr,desc_a,info)
-    ilaggr = -1
+    ilaggr = ilaggr_neginit
     call psb_geasb(ilaggr,desc_a,info)
     nr = a%get_nrows()
     nc = a%get_ncols()
@@ -213,7 +214,20 @@ contains
       call psb_barrier(ictxt)
       if (iam == 0) write(0,*)' out from buildmatching:', info
     end if
-
+    if (debug_mate) then
+      block
+        integer(psb_lpk_), allocatable :: ckmate(:)
+        allocate(ckmate(nr))
+        ckmate(1:nr) = mate(1:nr)
+        call psb_msort(ckmate(1:nr))
+        do i=1,nr-1
+          if ((ckmate(i)>0) .and. (ckmate(i) == ckmate(i+1))) then
+            write(0,*) iam,' Duplicate mate entry  at',i,' :',ckmate(i)
+          end if
+        end do
+      end block
+    end if
+     
     if (info == 0) then
       if (do_timings) call psb_tic(idx_phase2)
       if (debug_sync) then
@@ -259,7 +273,7 @@ contains
           cycle
         else
 
-          if (ilaggr(k) == -1) then
+          if (ilaggr(k) == ilaggr_neginit) then
 
             wk   = w(k)
             widx = w(idx)
@@ -267,7 +281,7 @@ contains
             nrmagg = wmax*sqrt((wk/wmax)**2+(widx/wmax)**2)
             if (nrmagg > epsilon(nrmagg)) then
               if (idx <= nr) then
-                if (ilaggr(idx) == -1) then
+                if (ilaggr(idx) == ilaggr_neginit) then
                   ! Now, if both vertices are local, the aggregate is local
                   ! (kinda obvious).
                   nlaggr(iam) = nlaggr(iam) + 1
@@ -275,6 +289,9 @@ contains
                   ilaggr(idx) = nlaggr(iam)
                   wtemp(k)    = w(k)/nrmagg
                   wtemp(idx)  = w(idx)/nrmagg
+                else
+                  write(0,*) iam,' Inconsistent mate? ',k,mate(k),idx,&
+                       &mate(idx),ilaggr(idx)
                 end if
                 nlpairs = nlpairs+1
               else if (idx <= nc) then
@@ -294,7 +311,7 @@ contains
                     ilaggr(k)   = nlaggr(iam)
                     nlpairs = nlpairs+1
                   else
-                    ilaggr(k) = -2
+                    ilaggr(k) = ilaggr_nonlocal
                   end if
                 else
                   ! Use a statistically unbiased tie-breaking rule,
@@ -309,7 +326,7 @@ contains
                     ilaggr(k)   = nlaggr(iam)
                     nlpairs = nlpairs+1
                   else
-                    ilaggr(k) = -2
+                    ilaggr(k) = ilaggr_nonlocal
                   end if
                 end if
               end if
@@ -325,6 +342,12 @@ contains
                 nlsingl     = nlsingl + 1
               end if
             end if
+            if (ilaggr(k) == ilaggr_neginit) then
+              write(0,*) iam,' Error: no update to ',k,mate(k),&
+                   & abs(w(k)),nrmagg,epsilon(nrmagg),wtemp(k)
+            end if
+          else
+            if (ilaggr(k)<0) write(0,*) 'Strange? ',k,ilaggr(k)
           end if
         end if
       end do
@@ -332,7 +355,7 @@ contains
       if (do_timings) call psb_tic(idx_phase3)
 
       ! Ok, now compute offsets, gather halo and fix non-local
-      ! aggregates  (those where ilaggr == -2)
+      ! aggregates  (those where ilaggr == ilaggr_nonlocal)
       call psb_sum(ictxt,nlaggr)
       ntaggr  = sum(nlaggr(0:np-1))
       naggrm1 = sum(nlaggr(0:iam-1))
@@ -347,7 +370,7 @@ contains
       call psb_halo(wtemp,desc_a,info)
       ! Cleanup as yet unmarked entries
       do k=1,nr
-        if (ilaggr(k) == -2) then
+        if (ilaggr(k) == ilaggr_nonlocal) then
           idx = mate(k)
           if (idx > nr) then
             i   = ilaggr(idx)
@@ -359,9 +382,14 @@ contains
           else
             write(0,*) 'Error : unresolved (paired) index ',k,idx,i,nr,nc, ilv(k),ilv(idx)
           end if
-        end if
-        if (ilaggr(k) <0) then
-          write(0,*) 'Matchboxp: Funny number: ',k,ilv(k),ilaggr(k),wtemp(k)
+        else if (ilaggr(k) <0) then
+          write(0,*) iam,'Matchboxp: Funny number: ',k,ilv(k),ilaggr(k),wtemp(k)
+          write(0,*) iam,'         :             : ',nr,nc,mate(k)
+          if (mate(k) <= nr) then
+            write(0,*) iam,'         :             : ',ilaggr(mate(k)),mate(mate(k)),&
+                 & ilv(k),ilv(mate(k)), ilv(mate(mate(k))),ilaggr(mate(mate(k)))
+          end if
+          flush(0)
         end if
       end do
       if (debug_sync) then
@@ -414,7 +442,7 @@ contains
 
       end block
       if (iam == 0) then
-        write(0,*) 'Matching statistics: Unmatched nodes ',&
+        write(0,*) iam,'Matching statistics: Unmatched nodes ',&
              & nunmatched,' Singletons:',nlsingl,' Pairs:',nlpairs
       end if
 

amgprec/amg_d_mumps_solver.F90

@@ -313,22 +313,24 @@ subroutine d_mumps_solver_finalize(sv)
 
 end subroutine d_mumps_solver_finalize
 
-subroutine d_mumps_solver_descr(sv,info,iout,coarse)
+subroutine d_mumps_solver_descr(sv,info,iout,coarse,prefix)
 
   Implicit None
 
   ! Arguments
   class(amg_d_mumps_solver_type), intent(in) :: sv
-  integer(psb_ipk_), intent(out)                     :: info
-  integer(psb_ipk_), intent(in), optional            :: iout
-  logical, intent(in), optional       :: coarse
-
+  integer(psb_ipk_), intent(out)               :: info
+  integer(psb_ipk_), intent(in), optional      :: iout
+  logical, intent(in), optional                :: coarse
+  character(len=*), intent(in), optional       :: prefix
+      
   ! Local variables
   integer(psb_ipk_)   :: err_act
   type(psb_ctxt_type) :: ctxt
   integer(psb_ipk_)   :: me, np
   character(len=20), parameter :: name='amg_z_mumps_solver_descr'
   integer(psb_ipk_) :: iout_
+  character(1024)    :: prefix_
 
   call psb_erractionsave(err_act)
   info = psb_success_
@@ -337,8 +339,13 @@ subroutine d_mumps_solver_descr(sv,info,iout,coarse)
   else
     iout_ = psb_out_unit
   endif
+  if (present(prefix)) then
+    prefix_ = prefix
+  else
+    prefix_ = ""
+  end if
 
-  write(iout_,*) '  MUMPS  Solver. '
+  write(iout_,*) trim(prefix_), '  MUMPS  Solver. '
 
   call psb_erractionrestore(err_act)
   return

amgprec/amg_d_newmatch_aggregator_mod.F90

@@ -0,0 +1,585 @@
+!
+!
+!  The aggregator object hosts the aggregation method for building
+!  the multilevel hierarchy. This variant is based on the hybrid method
+!  presented in 
+!
+!
+!   sm           -  class(amg_T_base_smoother_type), allocatable
+!                   The current level preconditioner (aka smoother).
+!   parms        -  type(amg_RTml_parms)
+!                   The parameters defining the multilevel strategy.
+!   ac           -  The local part of the current-level matrix, built by
+!                   coarsening the previous-level matrix.
+!   desc_ac      -  type(psb_desc_type).
+!                   The communication descriptor associated to the matrix
+!                   stored in ac.
+!   base_a       -  type(psb_Tspmat_type), pointer.
+!                   Pointer (really a pointer!) to the local part of the current 
+!                   matrix (so we have a unified treatment of residuals).
+!                   We need this to avoid passing explicitly the current matrix
+!                   to the routine which applies the preconditioner.
+!   base_desc    -  type(psb_desc_type), pointer.
+!                   Pointer to the communication descriptor associated to the
+!                   matrix pointed by base_a.
+!   map          -  Stores the maps (restriction and prolongation) between the
+!                   vector spaces associated to the index spaces of the previous
+!                   and current levels.
+!
+!   Methods:  
+!     Most methods follow the encapsulation hierarchy: they take whatever action
+!     is appropriate for the current object, then call the corresponding method for
+!     the contained object.
+!     As an example: the descr() method prints out a description of the
+!     level. It starts by invoking the descr() method of the parms object,
+!     then calls the descr() method of the smoother object. 
+!
+!    descr      -   Prints a description of the object.
+!    default    -   Set default values
+!    dump       -   Dump to file object contents
+!    set        -   Sets various parameters; when a request is unknown
+!                   it is passed to the smoother object for further processing.
+!    check      -   Sanity checks.
+!    sizeof     -   Total memory occupation in bytes
+!    get_nzeros -   Number of nonzeros 
+!
+!
+
+module amg_d_newmatch_aggregator_mod
+  use amg_d_base_aggregator_mod
+  use iso_c_binding
+  
+  type, bind(c)::  nwm_Vector
+    type(c_ptr) :: data
+    integer(c_int) :: size
+    integer(c_int) :: owns_data
+  end type nwm_Vector
+  
+  type, bind(c)::  nwm_CSRMatrix
+    type(c_ptr) :: i
+    type(c_ptr) :: j
+    integer(c_int) :: num_rows
+    integer(c_int) :: num_cols
+    integer(c_int) :: num_nonzeros
+    integer(c_int) :: owns_data
+    type(c_ptr) :: data
+  end type nwm_CSRMatrix
+  
+  type, extends(amg_d_base_aggregator_type) :: amg_d_newmatch_aggregator_type
+    integer(psb_ipk_) :: matching_alg
+    integer(psb_ipk_) :: n_sweeps
+    !
+    !  Note: the BootCMatch kernel we invoke  overwrites
+    !  the W argument with its update. Hence, copy it in w_nxt
+    !  before passing it to the matching
+    !
+    integer(psb_ipk_) :: orig_aggr_size
+    integer(psb_ipk_) :: jacobi_sweeps
+    real(psb_dpk_), allocatable :: w(:), w_nxt(:)
+    type(psb_dspmat_type), allocatable  :: prol, restr
+    type(psb_dspmat_type), allocatable  :: ac, base_a, rwa
+    type(psb_desc_type), allocatable    :: desc_ac, desc_ax, base_desc, rwdesc
+    type(nwm_Vector)  :: w_c_nxt
+    integer(psb_ipk_) :: max_csize
+    integer(psb_ipk_) :: max_nlevels
+    real(psb_dpk_)    :: lambda
+    logical           :: reproducible_matching = .false.
+    logical           :: need_symmetrize       = .false.
+    logical           :: unsmoothed_hierarchy  = .true.
+    logical           :: parallel_matching     = .true. 
+  contains
+    procedure, pass(ag) :: bld_tprol    => amg_d_newmatch_aggregator_build_tprol
+    procedure, pass(ag) :: csetc        => amg_d_newmatch_aggr_csetc
+    procedure, pass(ag) :: csetr        => amg_d_newmatch_aggr_csetr
+    procedure, pass(ag) :: cseti        => d_newmatch_aggr_cseti
+    procedure, pass(ag) :: default      => d_newmatch_aggr_set_default
+    procedure, pass(ag) :: mat_asb      => amg_d_newmatch_aggregator_mat_asb
+    procedure, pass(ag) :: mat_bld      => amg_d_newmatch_aggregator_mat_bld
+    procedure, pass(ag) :: inner_mat_asb => amg_d_newmatch_aggregator_inner_mat_asb
+    procedure, pass(ag) :: update_next  => d_newmatch_aggregator_update_next
+    procedure, pass(ag) :: bld_wnxt     => d_newmatch_bld_wnxt
+    procedure, pass(ag) :: bld_default_w    => d_bld_default_w
+    procedure, pass(ag) :: set_c_default_w  => d_set_default_nwm_w
+    procedure, pass(ag) :: descr        => d_newmatch_aggregator_descr
+    procedure, pass(ag) :: clone        => d_newmatch_aggregator_clone
+    procedure, pass(ag) :: free         => d_newmatch_aggregator_free
+    procedure, nopass   :: fmt          => d_newmatch_aggregator_fmt
+  end type amg_d_newmatch_aggregator_type
+
+
+  interface
+    subroutine  amg_d_newmatch_aggregator_build_tprol(ag,parms,ag_data,&
+         & a,desc_a,ilaggr,nlaggr,t_prol,info)
+      import :: amg_d_newmatch_aggregator_type, psb_desc_type, &
+           & psb_dspmat_type, psb_ldspmat_type, psb_dpk_,  &
+           & psb_ipk_, psb_lpk_, psb_epk_, amg_dml_parms, amg_daggr_data
+      implicit none
+      class(amg_d_newmatch_aggregator_type), target, intent(inout) :: ag
+      type(amg_dml_parms), intent(inout)   :: parms 
+      type(amg_daggr_data), intent(in)     :: ag_data
+      type(psb_dspmat_type), intent(inout) :: a
+      type(psb_desc_type), intent(inout)   :: desc_a
+      integer(psb_lpk_), allocatable, intent(out) :: ilaggr(:),nlaggr(:)
+      type(psb_ldspmat_type), intent(out)  :: t_prol
+      integer(psb_ipk_), intent(out)      :: info
+    end subroutine amg_d_newmatch_aggregator_build_tprol
+  end interface
+
+  interface
+    subroutine  amg_d_newmatch_aggregator_mat_bld(ag,parms,a,desc_a,ilaggr,nlaggr,&
+         & ac,desc_ac,op_prol,op_restr,t_prol,info)
+      import :: amg_d_newmatch_aggregator_type, psb_desc_type, &
+           & psb_dspmat_type, psb_ldspmat_type, psb_dpk_,  &
+           & psb_ipk_, psb_lpk_, psb_epk_, amg_dml_parms
+      implicit none
+      class(amg_d_newmatch_aggregator_type), target, intent(inout) :: ag
+      type(amg_dml_parms), intent(inout)    :: parms 
+      type(psb_dspmat_type), intent(in)     :: a
+      type(psb_desc_type), intent(inout)    :: desc_a
+      integer(psb_lpk_), intent(inout)      :: ilaggr(:), nlaggr(:)
+      type(psb_dspmat_type), intent(out)    :: op_prol,ac,op_restr
+      type(psb_ldspmat_type), intent(inout) :: t_prol
+      type(psb_desc_type), intent(inout)    :: desc_ac
+      integer(psb_ipk_), intent(out)        :: info
+    end subroutine amg_d_newmatch_aggregator_mat_bld
+  end interface
+
+  interface
+    subroutine  amg_d_newmatch_aggregator_mat_asb(ag,parms,a,desc_a,&
+         & ac,desc_ac, op_prol,op_restr,info)
+      import :: amg_d_newmatch_aggregator_type, psb_desc_type, &
+           & psb_dspmat_type, psb_ldspmat_type, psb_dpk_,  &
+           & psb_ipk_, psb_lpk_, psb_epk_, amg_dml_parms
+      implicit none
+      class(amg_d_newmatch_aggregator_type), target, intent(inout) :: ag
+      type(amg_dml_parms), intent(inout)    :: parms 
+      type(psb_dspmat_type), intent(in)     :: a
+      type(psb_desc_type), intent(inout)    :: desc_a
+      type(psb_dspmat_type), intent(inout) :: op_prol,ac,op_restr
+      type(psb_desc_type), intent(inout)     :: desc_ac
+      integer(psb_ipk_), intent(out)        :: info
+    end subroutine amg_d_newmatch_aggregator_mat_asb
+  end interface  
+  
+
+  interface
+    subroutine amg_d_newmatch_map_to_tprol(desc_a,ilaggr,nlaggr,valaggr, op_prol,info)
+      import :: amg_d_newmatch_aggregator_type, psb_desc_type, &
+           & psb_dspmat_type, psb_ldspmat_type, psb_dpk_,  &
+           & psb_ipk_, psb_lpk_, psb_epk_, amg_dml_parms
+      implicit none
+      type(psb_desc_type), intent(in)    :: desc_a
+      integer(psb_lpk_), allocatable, intent(inout)  :: ilaggr(:),nlaggr(:)
+      real(psb_dpk_), allocatable, intent(inout)  :: valaggr(:)
+      type(psb_ldspmat_type), intent(out)  :: op_prol
+      integer(psb_ipk_), intent(out)               :: info
+    end subroutine amg_d_newmatch_map_to_tprol
+  end interface
+  
+  interface
+    subroutine amg_daggrmat_unsmth_spmm_asb(a,desc_a,ilaggr,nlaggr,parms,&
+         & ac,op_prol,op_restr,info)
+      import :: amg_d_newmatch_aggregator_type, psb_desc_type, &
+           & psb_dspmat_type, psb_ldspmat_type, psb_dpk_,  &
+           & psb_ipk_, psb_lpk_, psb_epk_, amg_dml_parms
+      implicit none
+      type(psb_dspmat_type), intent(in)        :: a
+      type(psb_desc_type), intent(in)            :: desc_a
+      integer(psb_lpk_), intent(inout)           :: ilaggr(:), nlaggr(:)
+      type(amg_dml_parms), intent(inout)      :: parms 
+      type(psb_ldspmat_type), intent(inout)     :: op_prol
+      type(psb_ldspmat_type), intent(out)       :: ac,op_restr
+      integer(psb_ipk_), intent(out)             :: info
+    end subroutine amg_daggrmat_unsmth_spmm_asb
+  end interface
+
+  interface
+    subroutine  amg_d_newmatch_aggregator_inner_mat_asb(ag,parms,a,desc_a,&
+         & ac,desc_ac, op_prol,op_restr,info)
+      import :: amg_d_newmatch_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
+      class(amg_d_newmatch_aggregator_type), target, intent(inout) :: ag
+      type(amg_dml_parms), intent(inout)    :: parms
+      type(psb_dspmat_type), intent(in)     :: a
+      type(psb_desc_type), intent(in)       :: desc_a
+      type(psb_dspmat_type), intent(inout) :: op_prol,op_restr
+      type(psb_dspmat_type), intent(inout)  :: ac
+      type(psb_desc_type), intent(inout)    :: desc_ac
+      integer(psb_ipk_), intent(out)        :: info
+    end subroutine amg_d_newmatch_aggregator_inner_mat_asb
+  end interface
+
+  
+!!$  interface
+!!$    subroutine amg_d_newmatch_unsmth_spmm_bld(a,desc_a,ilaggr,nlaggr,parms,&
+!!$         & ac,desc_ac,op_prol,op_restr,t_prol,info)
+!!$      import :: amg_d_newmatch_aggregator_type, psb_desc_type, &
+!!$           & psb_dspmat_type, psb_ldspmat_type, psb_dpk_,  &
+!!$           & psb_ipk_, psb_lpk_, psb_epk_, amg_dml_parms
+!!$      
+!!$      implicit none
+!!$      
+!!$      ! Arguments
+!!$      type(psb_dspmat_type), intent(in)        :: a
+!!$      type(psb_desc_type), intent(in)            :: desc_a
+!!$      integer(psb_lpk_), intent(inout)           :: ilaggr(:), nlaggr(:)
+!!$      type(amg_dml_parms), intent(inout)      :: parms 
+!!$      type(psb_ldspmat_type), intent(inout)     :: t_prol
+!!$      type(psb_dspmat_type), intent(inout)       :: op_prol,ac,op_restr
+!!$      type(psb_desc_type), intent(inout)    :: desc_ac
+!!$      integer(psb_ipk_), intent(out)             :: info
+!!$    end subroutine amg_d_newmatch_unsmth_spmm_bld
+!!$  end interface
+
+  interface
+    subroutine amg_d_newmatch_spmm_bld_ov(a,desc_a,ilaggr,nlaggr,parms,&
+         & ac,desc_ac,op_prol,op_restr,t_prol,info)
+      import :: amg_d_newmatch_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
+      type(psb_dspmat_type), intent(inout)    :: a
+      type(psb_desc_type), intent(inout)       :: desc_a
+      integer(psb_lpk_), intent(inout)        :: ilaggr(:), nlaggr(:)
+      type(amg_dml_parms), intent(inout)      :: parms
+      type(psb_ldspmat_type), intent(inout)   :: t_prol
+      type(psb_dspmat_type), intent(inout)    :: op_prol,ac, op_restr
+      type(psb_desc_type), intent(out)          :: desc_ac
+      integer(psb_ipk_), intent(out)          :: info
+    end subroutine amg_d_newmatch_spmm_bld_ov
+  end interface
+
+  interface
+    subroutine amg_d_newmatch_spmm_bld_inner(a,desc_a,ilaggr,nlaggr,parms,&
+         & ac,desc_ac,op_prol,op_restr,t_prol,info)
+      import :: amg_d_newmatch_aggregator_type, psb_desc_type, psb_dspmat_type,&
+           & psb_ldspmat_type, psb_dpk_, psb_ipk_, psb_lpk_, amg_dml_parms, amg_daggr_data,&
+           & psb_d_csr_sparse_mat, psb_ld_csr_sparse_mat
+      implicit none
+      type(psb_d_csr_sparse_mat), intent(inout) :: a
+      type(psb_desc_type), intent(inout)          :: desc_a
+      integer(psb_lpk_), intent(inout)          :: ilaggr(:), nlaggr(:)
+      type(amg_dml_parms), intent(inout)        :: parms
+      type(psb_ldspmat_type), intent(inout)     :: t_prol
+      type(psb_dspmat_type), intent(inout)      :: op_prol,ac, op_restr
+      type(psb_desc_type), intent(out)          :: desc_ac
+      integer(psb_ipk_), intent(out)            :: info
+    end subroutine amg_d_newmatch_spmm_bld_inner
+  end interface
+  
+  private :: is_legal_malg, is_legal_csize, is_legal_nsweeps, is_legal_nlevels
+  
+
+contains
+
+  subroutine d_bld_default_w(ag,nr)
+    use psb_realloc_mod
+    implicit none 
+    class(amg_d_newmatch_aggregator_type), target, intent(inout) :: ag
+    integer(psb_ipk_), intent(in) :: nr
+    integer(psb_ipk_) :: info
+    call psb_realloc(nr,ag%w,info)
+    if (info /= psb_success_) return
+    ag%w = done
+    call ag%set_c_default_w()
+  end subroutine d_bld_default_w
+
+  subroutine d_set_default_nwm_w(ag)
+    use psb_realloc_mod
+    use iso_c_binding
+    implicit none 
+    class(amg_d_newmatch_aggregator_type), target, intent(inout) :: ag
+    integer(psb_ipk_) :: info
+    
+    call psb_safe_ab_cpy(ag%w,ag%w_nxt,info)
+    ag%w_c_nxt%size      = psb_size(ag%w_nxt)
+    ag%w_c_nxt%owns_data = 0
+    if (ag%w_c_nxt%size > 0) call set_cloc(ag%w_nxt, ag%w_c_nxt)
+
+  end subroutine d_set_default_nwm_w
+
+  subroutine set_cloc(vect,w_c_nxt)
+    use iso_c_binding
+    real(psb_dpk_), target :: vect(:)
+    type(nwm_Vector) :: w_c_nxt
+    
+    w_c_nxt%data = c_loc(vect)
+  end subroutine set_cloc
+
+
+  subroutine d_newmatch_bld_wnxt(ag,ilaggr,valaggr,nx)
+    use psb_realloc_mod
+    implicit none 
+    class(amg_d_newmatch_aggregator_type), target, intent(inout) :: ag
+    integer(psb_lpk_), intent(in)  :: ilaggr(:)
+    real(psb_dpk_), intent(in)     :: valaggr(:)
+    integer(psb_ipk_), intent(in)  :: nx
+
+    integer(psb_ipk_) :: info,i,j
+
+    ! The vector was already fixed in the call to Newmatch.
+    call psb_realloc(nx,ag%w_nxt,info)
+    
+  end subroutine d_newmatch_bld_wnxt
+  
+  function d_newmatch_aggregator_fmt() result(val)
+    implicit none 
+    character(len=32)  :: val
+
+    val = "new matching aggregation"
+  end function d_newmatch_aggregator_fmt
+
+  subroutine  d_newmatch_aggregator_descr(ag,parms,iout,info,prefix)
+    implicit none 
+    class(amg_d_newmatch_aggregator_type), intent(in) :: ag
+    type(amg_dml_parms), intent(in)   :: parms
+    integer(psb_ipk_), intent(in)  :: iout
+    integer(psb_ipk_), intent(out) :: info
+    character(len=*), intent(in), optional  :: prefix
+    character(1024)    :: prefix_
+    if (present(prefix)) then
+      prefix_ = prefix
+    else
+      prefix_ = ""
+    end if
+
+    write(iout,*) trim(prefix_),' ','NewMatch Aggregator'
+    write(iout,*) trim(prefix_),' ','   Number of Matching   sweeps: ',ag%n_sweeps
+    write(iout,*) trim(prefix_),' ','   Matching algorithm         : ',ag%matching_alg
+    write(iout,*) trim(prefix_),' ','Aggregator object type: ',ag%fmt()
+    call parms%mldescr(iout,info)
+    
+    return
+  end subroutine d_newmatch_aggregator_descr
+
+  function is_legal_malg(alg) result(val)
+    logical :: val
+    integer(psb_ipk_) :: alg
+
+    val = ((0<=alg).and.(alg<=2))
+  end function is_legal_malg
+
+  function is_legal_csize(csize) result(val)
+    logical :: val
+    integer(psb_ipk_) :: csize
+
+    val = ((-1==csize).or.(csize >0))
+  end function is_legal_csize
+
+  function is_legal_nsweeps(nsw) result(val)
+    logical :: val
+    integer(psb_ipk_) :: nsw
+
+    val = (1<=nsw)
+  end function is_legal_nsweeps
+
+  function is_legal_nlevels(nlv) result(val)
+    logical :: val
+    integer(psb_ipk_) :: nlv
+
+    val = (1<=nlv)
+  end function is_legal_nlevels
+  
+  subroutine  d_newmatch_aggregator_update_next(ag,agnext,info)
+    use psb_realloc_mod
+    implicit none 
+    class(amg_d_newmatch_aggregator_type), target, intent(inout) :: ag
+    class(amg_d_base_aggregator_type), target, intent(inout) :: agnext
+    integer(psb_ipk_), intent(out)       :: info
+
+    !
+    !
+    select type(agnext)
+    class is (amg_d_newmatch_aggregator_type)
+      if (.not.is_legal_malg(agnext%matching_alg)) &
+           & agnext%matching_alg = ag%matching_alg
+      if (.not.is_legal_nsweeps(agnext%n_sweeps))&
+           & agnext%n_sweeps     = ag%n_sweeps 
+      if (.not.is_legal_csize(agnext%max_csize))&
+           & agnext%max_csize    = ag%max_csize 
+      if (.not.is_legal_nlevels(agnext%max_nlevels))&
+           & agnext%max_nlevels  = ag%max_nlevels
+      ! Is this going to generate shallow copies/memory leaks/double frees?
+      ! To be investigated further.
+      call psb_safe_ab_cpy(ag%w_nxt,agnext%w,info)
+      call agnext%set_c_default_w()
+    class default
+      ! What should we do here? 
+    end select
+    info = 0 
+  end subroutine d_newmatch_aggregator_update_next
+
+  
+  subroutine amg_d_newmatch_aggr_csetr(ag,what,val,info,idx)
+
+    Implicit None
+
+    ! Arguments
+    class(amg_d_newmatch_aggregator_type), intent(inout) :: ag
+    character(len=*), intent(in)               :: what
+    real(psb_dpk_), intent(in)                 :: val
+    integer(psb_ipk_), intent(out)             :: info
+    integer(psb_ipk_), intent(in), optional    :: idx
+    integer(psb_ipk_)  :: err_act, iwhat
+    character(len=20)  :: name='d_newmatch_aggr_csetr'
+    info = psb_success_
+
+    ! For now we ignore IDX
+
+    select case(psb_toupper(trim(what)))
+    case('NWM_LAMBDA')
+      ag%lambda = val
+    case default
+      ! Do nothing
+    end select
+    return
+  end subroutine amg_d_newmatch_aggr_csetr
+  
+  subroutine amg_d_newmatch_aggr_csetc(ag,what,val,info,idx)
+
+    Implicit None
+
+    ! Arguments
+    class(amg_d_newmatch_aggregator_type), intent(inout) :: ag
+    character(len=*), intent(in)                  :: what
+    character(len=*), intent(in)                 :: val
+    integer(psb_ipk_), intent(out)                :: info
+    integer(psb_ipk_), intent(in), optional       :: idx
+    integer(psb_ipk_)  :: err_act, iwhat
+    character(len=20)  :: name='d_newmatch_aggr_csetc'
+    info = psb_success_
+
+    ! For now we ignore IDX
+
+    select case(psb_toupper(trim(what)))
+    case('NWM_PARALLEL_MATCHING')
+      select case(psb_toupper(trim(val)))
+      case('SEQUENTIAL','F','FALSE')
+        ag%parallel_matching = .false.
+      case('PARALLEL','TRUE','T')
+        ag%parallel_matching =.true.
+      end select
+    case('NWM_REPRODUCIBLE_MATCHING')
+      select case(psb_toupper(trim(val)))
+      case('F','FALSE')
+        ag%reproducible_matching = .false.
+      case('REPRODUCIBLE','TRUE','T')
+        ag%reproducible_matching =.true.
+      end select
+    case('NWM_NEED_SYMMETRIZE')
+      select case(psb_toupper(trim(val)))
+      case('FALSE','F')
+        ag%need_symmetrize = .false.
+      case('SYMMETRIZE','TRUE','T')
+        ag%need_symmetrize =.true.
+      end select
+    case('NWM_UNSMOOTHED_HIERARCHY')
+      select case(psb_toupper(trim(val)))
+      case('F','FALSE')
+        ag%unsmoothed_hierarchy = .false.
+      case('T','TRUE')
+        ag%unsmoothed_hierarchy =.true.
+      end select
+    case default
+      ! Do nothing
+    end select
+    return
+  end subroutine amg_d_newmatch_aggr_csetc
+  
+  subroutine d_newmatch_aggr_cseti(ag,what,val,info,idx)
+
+    Implicit None
+
+    ! Arguments
+    class(amg_d_newmatch_aggregator_type), intent(inout) :: ag
+    character(len=*), intent(in)                  :: what
+    integer(psb_ipk_), intent(in)                 :: val
+    integer(psb_ipk_), intent(out)                :: info
+    integer(psb_ipk_), intent(in), optional       :: idx
+    integer(psb_ipk_)  :: err_act, iwhat
+    character(len=20)  :: name='d_newmatch_aggr_cseti'
+    info = psb_success_
+
+    ! For now we ignore IDX
+    
+    select case(psb_toupper(trim(what)))
+    case('NWM_MATCH_ALG','NWM_MATCHING_ALG')
+      ag%matching_alg = val
+    case('NWM_SWEEPS')
+      ag%n_sweeps=val
+    case('NWM_MAX_CSIZE')
+      ag%max_csize = val
+    case('NWM_MAX_NLEVELS')
+      ag%max_nlevels = val
+    case('NWM_W_SIZE')
+      call ag%bld_default_w(val)
+    case('AGGR_SIZE')
+      ag%orig_aggr_size = val
+      ag%n_sweeps=max(1,ceiling(log(val*1.0)/log(2.0)))
+    case default
+      ! Do nothing
+    end select
+    return
+  end subroutine d_newmatch_aggr_cseti
+
+  subroutine d_newmatch_aggr_set_default(ag)
+
+    Implicit None
+
+    ! Arguments
+    class(amg_d_newmatch_aggregator_type), intent(inout) :: ag
+    character(len=20)  :: name='d_newmatch_aggr_set_default'
+    ag%matching_alg = 1
+    ag%n_sweeps     = 1
+    ag%max_nlevels  = 36
+    ag%max_csize    = -1
+    ag%lambda       = -1
+    !
+    ! Apparently newMatch works better
+    ! by keeping all entries
+    ! 
+    ag%do_clean_zeros = .false.
+
+    return
+
+  end subroutine d_newmatch_aggr_set_default
+
+  subroutine  d_newmatch_aggregator_free(ag,info)
+    use iso_c_binding
+    implicit none 
+    class(amg_d_newmatch_aggregator_type), intent(inout) :: ag
+    integer(psb_ipk_), intent(out)       :: info
+
+    info = 0
+    if (allocated(ag%w)) deallocate(ag%w,stat=info)
+    if (info /= 0) return
+    if (allocated(ag%w_nxt)) deallocate(ag%w_nxt,stat=info)
+    if (info /= 0) return
+    ag%w_c_nxt%size = 0
+    ag%w_c_nxt%data = c_null_ptr
+    ag%w_c_nxt%owns_data = 0
+  end subroutine d_newmatch_aggregator_free
+
+  subroutine  d_newmatch_aggregator_clone(ag,agnext,info)
+    implicit none 
+    class(amg_d_newmatch_aggregator_type), intent(inout) :: ag
+    class(amg_d_base_aggregator_type), allocatable, intent(inout) :: agnext
+    integer(psb_ipk_), intent(out)       :: info
+
+    info = 0 
+    if (allocated(agnext)) then
+      call agnext%free(info)
+      if (info == 0) deallocate(agnext,stat=info)
+    end if
+    if (info /= 0) return
+    allocate(agnext,source=ag,stat=info)
+    select type(agnext)
+    class is (amg_d_newmatch_aggregator_type)
+      call agnext%set_c_default_w()
+    class default
+      ! Should never ever get here
+      info = -1
+    end select
+  end subroutine d_newmatch_aggregator_clone  
+  
+end module amg_d_newmatch_aggregator_mod

amgprec/amg_d_onelev_mod.f90

@@ -57,6 +57,7 @@ module amg_d_onelev_mod
   use amg_d_base_smoother_mod
   use amg_d_dec_aggregator_mod
   use amg_d_parmatch_aggregator_mod
+  use amg_d_newmatch_aggregator_mod
 
   use psb_base_mod, only : psb_dspmat_type, psb_d_vect_type, &
        & psb_d_base_vect_type, psb_ldspmat_type,  psb_dlinmap_type, psb_dpk_, &
@@ -258,7 +259,7 @@ module amg_d_onelev_mod
   end interface
 
   interface
-    subroutine amg_d_base_onelev_descr(lv,il,nl,ilmin,info,iout, verbosity)
+    subroutine amg_d_base_onelev_descr(lv,il,nl,ilmin,info,iout, verbosity,prefix)
       import :: psb_dspmat_type, psb_d_vect_type, psb_d_base_vect_type, &
            & psb_dlinmap_type, psb_dpk_, amg_d_onelev_type, &
            & psb_ipk_, psb_epk_, psb_desc_type
@@ -269,6 +270,7 @@ module amg_d_onelev_mod
       integer(psb_ipk_), intent(out)          :: info
       integer(psb_ipk_), intent(in), optional :: iout
       integer(psb_ipk_), intent(in), optional :: verbosity
+      character(len=*), intent(in), optional  :: prefix
     end subroutine amg_d_base_onelev_descr
   end interface
 

amgprec/amg_d_parmatch_aggregator_mod.F90

@@ -390,18 +390,25 @@ contains
 
   end function amg_d_parmatch_aggregator_sizeof
 
-  subroutine  amg_d_parmatch_aggregator_descr(ag,parms,iout,info)
+  subroutine  amg_d_parmatch_aggregator_descr(ag,parms,iout,info,prefix)
     implicit none
     class(amg_d_parmatch_aggregator_type), intent(in) :: ag
     type(amg_dml_parms), intent(in)   :: parms
     integer(psb_ipk_), intent(in)  :: iout
     integer(psb_ipk_), intent(out) :: info
+    character(len=*), intent(in), optional  :: prefix
+    character(1024)    :: prefix_
+    if (present(prefix)) then
+      prefix_ = prefix
+    else
+      prefix_ = ""
+    end if
 
-    write(iout,*) 'Parallel Matching Aggregator'
-    write(iout,*) '   Number of matching  sweeps: ',ag%n_sweeps
-    write(iout,*) '   Matching algorithm         : MatchBoxP (PREIS)'
-    write(iout,*) 'Aggregator object type: ',ag%fmt()
-    call parms%mldescr(iout,info)
+    write(iout,*) trim(prefix_),' ','Parallel Matching Aggregator'
+    write(iout,*) trim(prefix_),' ','   Number of matching  sweeps: ',ag%n_sweeps
+    write(iout,*) trim(prefix_),' ','   Matching algorithm         : MatchBoxP (PREIS)'
+    write(iout,*) trim(prefix_),' ','Aggregator object type: ',ag%fmt()
+    call parms%mldescr(iout,info,prefix=prefix)
 
     return
   end subroutine amg_d_parmatch_aggregator_descr

amgprec/amg_d_prec_type.f90

@@ -155,15 +155,16 @@ module amg_d_prec_type
 
 
   interface amg_precdescr
-    subroutine amg_dfile_prec_descr(prec,info,iout,root,verbosity)
+    subroutine amg_dfile_prec_descr(prec,info,iout,root,verbosity,prefix)
       import :: amg_dprec_type, psb_ipk_
       implicit none
       ! Arguments
-      class(amg_dprec_type), intent(in)      :: prec
-      integer(psb_ipk_), intent(out) :: info
+      class(amg_dprec_type), intent(in)     :: prec
+      integer(psb_ipk_), intent(out)          :: info
       integer(psb_ipk_), intent(in), optional :: iout
       integer(psb_ipk_), intent(in), optional :: root
       integer(psb_ipk_), intent(in), optional :: verbosity
+      character(len=*), intent(in), optional  :: prefix
     end subroutine amg_dfile_prec_descr
   end interface
 

amgprec/amg_d_slu_solver.F90

@@ -385,20 +385,22 @@ contains
   
   end subroutine d_slu_solver_finalize
 
-  subroutine d_slu_solver_descr(sv,info,iout,coarse)
+  subroutine d_slu_solver_descr(sv,info,iout,coarse,prefix)
 
     Implicit None
 
     ! Arguments
     class(amg_d_slu_solver_type), intent(in) :: sv
-    integer, intent(out)                     :: info
-    integer, intent(in), optional            :: iout
-    logical, intent(in), optional       :: coarse
+    integer, intent(out)                       :: info
+    integer, intent(in), optional              :: iout
+    logical, intent(in), optional              :: coarse
+    character(len=*), intent(in), optional     :: prefix
 
     ! Local variables
     integer      :: err_act
     character(len=20), parameter :: name='amg_d_slu_solver_descr'
     integer :: iout_
+    character(1024)    :: prefix_
 
     call psb_erractionsave(err_act)
     info = psb_success_
@@ -407,8 +409,13 @@ contains
     else
       iout_ = psb_out_unit
     endif
+    if (present(prefix)) then
+      prefix_ = prefix
+    else
+      prefix_ = ""
+    end if
     
-    write(iout_,*) '  SuperLU Sparse Factorization Solver. '
+    write(iout_,*) trim(prefix_), '  SuperLU Sparse Factorization Solver. '
 
     call psb_erractionrestore(err_act)
     return

amgprec/amg_d_sludist_solver.F90

@@ -270,13 +270,13 @@ contains
     ! Local variables
     type(psb_dspmat_type) :: atmp
     type(psb_d_csr_sparse_mat) :: acsr
-    integer(psb_lpk_), allocatable :: gia(:), gja(:)
     type(psb_ctxt_type) :: ctxt
+    integer(psb_lpk_), allocatable :: gia(:), gja(:)
     integer(psb_lpk_) :: lfrst
     integer(psb_ipk_) :: n_row,n_col, nrow_a, nztota, nglob, nzt, npr, npc
     integer(psb_ipk_) :: ifrst, ibcheck
     integer(psb_ipk_) :: np,me,i, err_act, debug_unit, debug_level
-    character(len=20) :: name='d_sludist_solver_bld', ch_err
+    character(len=20)  :: name='d_sludist_solver_bld', ch_err
     
     info=psb_success_
     call psb_erractionsave(err_act)
@@ -325,7 +325,6 @@ contains
     acsr%ja(:)  = acsr%ja(:)  - 1
     acsr%irp(:) = acsr%irp(:) - 1
     ifrst = lfrst - 1
-    
     info = amg_dsludist_fact(nglob,nrow_a,nztota,ifrst,&
          & acsr%val,acsr%irp,acsr%ja,sv%lufactors,&
          & npr,npc)
@@ -422,15 +421,16 @@ contains
   
   end subroutine d_sludist_solver_finalize
 
-  subroutine d_sludist_solver_descr(sv,info,iout,coarse)
+  subroutine d_sludist_solver_descr(sv,info,iout,coarse,prefix)
 
     Implicit None
 
     ! Arguments
     class(amg_d_sludist_solver_type), intent(in) :: sv
-    integer, intent(out)                     :: info
-    integer, intent(in), optional            :: iout
-    logical, intent(in), optional       :: coarse
+    integer, intent(out)                           :: info
+    integer, intent(in), optional                  :: iout
+    logical, intent(in), optional                  :: coarse
+    character(len=*), intent(in), optional         :: prefix
 
     ! Local variables
     integer      :: err_act
@@ -438,6 +438,7 @@ contains
     integer      :: me, np
     character(len=20), parameter :: name='amg_d_sludist_solver_descr'
     integer :: iout_
+    character(1024)    :: prefix_
 
     call psb_erractionsave(err_act)
     info = psb_success_
@@ -446,8 +447,13 @@ contains
     else
       iout_ = psb_out_unit
     endif
+    if (present(prefix)) then
+      prefix_ = prefix
+    else
+      prefix_ = ""
+    end if
     
-    write(iout_,*) '  SuperLU_Dist Sparse Factorization Solver. '
+    write(iout_,*) trim(prefix_), '  SuperLU_Dist Sparse Factorization Solver. '
 
     call psb_erractionrestore(err_act)
     return

amgprec/amg_d_symdec_aggregator_mod.f90

@@ -88,16 +88,25 @@ contains
     val = "Symmetric Decoupled aggregation"
   end function amg_d_symdec_aggregator_fmt
   
-  subroutine  amg_d_symdec_aggregator_descr(ag,parms,iout,info)
+  subroutine  amg_d_symdec_aggregator_descr(ag,parms,iout,info,prefix)
     implicit none 
     class(amg_d_symdec_aggregator_type), intent(in) :: ag
     type(amg_dml_parms), intent(in)   :: parms
     integer(psb_ipk_), intent(in)  :: iout
     integer(psb_ipk_), intent(out) :: info
+    character(len=*), intent(in), optional  :: prefix
 
-    write(iout,*) 'Decoupled Aggregator locally-symmetrized'
-    write(iout,*) 'Aggregator object type: ',ag%fmt()
-    call parms%mldescr(iout,info)
+    character(1024)    :: prefix_
+
+    if (present(prefix)) then
+      prefix_ = prefix
+    else
+      prefix_ = ""
+    end if
+
+    write(iout,*) trim(prefix_),' ','Decoupled Aggregator locally-symmetrized'
+    write(iout,*) trim(prefix_),' ','Aggregator object type: ',ag%fmt()
+    call parms%mldescr(iout,info,prefix=prefix)
     
     return
   end subroutine amg_d_symdec_aggregator_descr

amgprec/amg_d_umf_solver.F90

@@ -390,20 +390,22 @@ contains
   
   end subroutine d_umf_solver_finalize
 
-  subroutine d_umf_solver_descr(sv,info,iout,coarse)
+  subroutine d_umf_solver_descr(sv,info,iout,coarse,prefix)
 
     Implicit None
 
     ! Arguments
     class(amg_d_umf_solver_type), intent(in) :: sv
-    integer, intent(out)                     :: info
-    integer, intent(in), optional            :: iout
-    logical, intent(in), optional       :: coarse
+    integer, intent(out)                       :: info
+    integer, intent(in), optional              :: iout
+    logical, intent(in), optional              :: coarse
+      character(len=*), intent(in), optional   :: prefix
 
     ! Local variables
     integer      :: err_act
     character(len=20), parameter :: name='amg_d_umf_solver_descr'
     integer :: iout_
+    character(1024)    :: prefix_
 
     call psb_erractionsave(err_act)
     info = psb_success_
@@ -412,8 +414,13 @@ contains
     else
       iout_ = psb_out_unit
     endif
+    if (present(prefix)) then
+      prefix_ = prefix
+    else
+      prefix_ = ""
+    end if
     
-    write(iout_,*) '  UMFPACK Sparse Factorization Solver. '
+    write(iout_,*) trim(prefix_), '  UMFPACK Sparse Factorization Solver. '
 
     call psb_erractionrestore(err_act)
     return

amgprec/amg_s_ainv_solver.F90

@@ -198,7 +198,7 @@ module amg_s_ainv_solver
 !!$  end interface
 
   interface
-    subroutine amg_s_ainv_solver_descr(sv,info,iout,coarse)
+    subroutine amg_s_ainv_solver_descr(sv,info,iout,coarse,prefix)
       import :: psb_dpk_, amg_s_ainv_solver_type, psb_ipk_
 
       Implicit None
@@ -208,7 +208,7 @@ module amg_s_ainv_solver
       integer(psb_ipk_), intent(out)            :: info
       integer(psb_ipk_), intent(in), optional   :: iout
       logical, intent(in), optional             :: coarse
-
+      character(len=*), intent(in), optional  :: prefix
     end subroutine amg_s_ainv_solver_descr
   end interface
 

amgprec/amg_s_as_smoother.f90

@@ -396,21 +396,23 @@ contains
   end subroutine s_as_smoother_default
 
 
-  subroutine s_as_smoother_descr(sm,info,iout,coarse)
+  subroutine s_as_smoother_descr(sm,info,iout,coarse,prefix)
 
     Implicit None
 
     ! Arguments
     class(amg_s_as_smoother_type), intent(in) :: sm
-    integer(psb_ipk_), intent(out)                      :: info
-    integer(psb_ipk_), intent(in), optional             :: iout
-    logical, intent(in), optional             :: coarse
+    integer(psb_ipk_), intent(out)              :: info
+    integer(psb_ipk_), intent(in), optional     :: iout
+    logical, intent(in), optional              :: coarse
+      character(len=*), intent(in), optional   :: prefix
 
     ! Local variables
     integer(psb_ipk_)      :: err_act
     character(len=20), parameter :: name='amg_s_as_smoother_descr'
     integer(psb_ipk_) :: iout_
     logical      :: coarse_
+    character(1024)    :: prefix_
 
     call psb_erractionsave(err_act)
     info = psb_success_
@@ -424,16 +426,21 @@ contains
     else
       iout_ = psb_out_unit
     endif
+    if (present(prefix)) then
+      prefix_ = prefix
+    else
+      prefix_ = ""
+    end if
 
     if (.not.coarse_) then 
-      write(iout_,*) '  Additive Schwarz with  ',&
+      write(iout_,*) trim(prefix_), '  Additive Schwarz with  ',&
            &  sm%novr, ' overlap layers.'
-      write(iout_,*) '  Restrictor:  ',restrict_names(sm%restr)
-      write(iout_,*) '  Prolongator: ',prolong_names(sm%prol)
-      write(iout_,*) '  Local solver:'
+      write(iout_,*) trim(prefix_), '  Restrictor:  ',restrict_names(sm%restr)
+      write(iout_,*) trim(prefix_), '  Prolongator: ',prolong_names(sm%prol)
+      write(iout_,*) trim(prefix_), '  Local solver:'
     endif
     if (allocated(sm%sv)) then 
-      call sm%sv%descr(info,iout_,coarse=coarse)
+      call sm%sv%descr(info,iout_,coarse=coarse,prefix=prefix)
     end if
 
     call psb_erractionrestore(err_act)

amgprec/amg_s_base_aggregator_mod.f90

@@ -275,15 +275,22 @@ contains
     val = .false.
   end function amg_s_base_aggregator_xt_desc
 
-  subroutine  amg_s_base_aggregator_descr(ag,parms,iout,info)
+  subroutine  amg_s_base_aggregator_descr(ag,parms,iout,info,prefix)
     implicit none 
     class(amg_s_base_aggregator_type), intent(in) :: ag
     type(amg_sml_parms), intent(in)   :: parms
     integer(psb_ipk_), intent(in)  :: iout
     integer(psb_ipk_), intent(out) :: info
+    character(len=*), intent(in), optional  :: prefix
+    character(1024)    :: prefix_
+    if (present(prefix)) then
+      prefix_ = prefix
+    else
+      prefix_ = ""
+    end if
 
-    write(iout,*) 'Aggregator object type: ',ag%fmt()
-    call parms%mldescr(iout,info)
+    write(iout,*) trim(prefix_),' ', 'Aggregator object type: ',ag%fmt()
+    call parms%mldescr(iout,info,prefix=prefix)
     
     return
   end subroutine amg_s_base_aggregator_descr

amgprec/amg_s_base_smoother_mod.f90

@@ -272,7 +272,7 @@ module amg_s_base_smoother_mod
   end interface
   
   interface 
-    subroutine amg_s_base_smoother_descr(sm,info,iout,coarse)
+    subroutine amg_s_base_smoother_descr(sm,info,iout,coarse,prefix)
       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_base_smoother_type, psb_ipk_
@@ -281,6 +281,7 @@ module amg_s_base_smoother_mod
       integer(psb_ipk_), intent(out)                :: info
       integer(psb_ipk_), intent(in), optional       :: iout
       logical, intent(in), optional                 :: coarse
+      character(len=*), intent(in), optional        :: prefix
     end subroutine amg_s_base_smoother_descr
   end interface
   

amgprec/amg_s_base_solver_mod.f90

@@ -270,7 +270,7 @@ module amg_s_base_solver_mod
   end interface
   
   interface
-    subroutine amg_s_base_solver_descr(sv,info,iout,coarse)
+    subroutine amg_s_base_solver_descr(sv,info,iout,coarse,prefix)
       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_base_solver_type, psb_ipk_
@@ -281,7 +281,7 @@ module amg_s_base_solver_mod
       integer(psb_ipk_), intent(out)              :: info
       integer(psb_ipk_), intent(in), optional     :: iout
       logical, intent(in), optional               :: coarse
-
+      character(len=*), intent(in), optional      :: prefix
     end subroutine amg_s_base_solver_descr
   end interface 
   

amgprec/amg_s_dec_aggregator_mod.f90

@@ -184,16 +184,23 @@ contains
     val = "Decoupled aggregation"
   end function amg_s_dec_aggregator_fmt
   
-  subroutine  amg_s_dec_aggregator_descr(ag,parms,iout,info)
+  subroutine  amg_s_dec_aggregator_descr(ag,parms,iout,info,prefix)
     implicit none 
     class(amg_s_dec_aggregator_type), intent(in) :: ag
     type(amg_sml_parms), intent(in)   :: parms
     integer(psb_ipk_), intent(in)  :: iout
     integer(psb_ipk_), intent(out) :: info
+    character(len=*), intent(in), optional  :: prefix
+    character(1024)    :: prefix_
+    if (present(prefix)) then
+      prefix_ = prefix
+    else
+      prefix_ = ""
+    end if
 
-    write(iout,*) 'Decoupled Aggregator'
-    write(iout,*) 'Aggregator object type: ',ag%fmt()
-    call parms%mldescr(iout,info)
+    write(iout,*) trim(prefix_),' ','Decoupled Aggregator'
+    write(iout,*) trim(prefix_),' ','Aggregator object type: ',ag%fmt()
+    call parms%mldescr(iout,info,prefix=prefix)
     
     return
   end subroutine amg_s_dec_aggregator_descr

amgprec/amg_s_diag_solver.f90

@@ -219,7 +219,7 @@ contains
 
   end subroutine s_diag_solver_free
 
-  subroutine s_diag_solver_descr(sv,info,iout,coarse)
+  subroutine s_diag_solver_descr(sv,info,iout,coarse,prefix)
 
     Implicit None
 
@@ -228,11 +228,13 @@ contains
     integer(psb_ipk_), intent(out)              :: info
     integer(psb_ipk_), intent(in), optional     :: iout
     logical, intent(in), optional               :: coarse
+    character(len=*), intent(in), optional      :: prefix
 
     ! Local variables
     integer(psb_ipk_)      :: err_act
     character(len=20), parameter :: name='amg_s_diag_solver_descr'
     integer(psb_ipk_) :: iout_
+    character(1024)    :: prefix_
 
     info = psb_success_
     if (present(iout)) then 
@@ -240,8 +242,13 @@ contains
     else
       iout_ = psb_out_unit
     endif
-    
-    write(iout_,*) '  Diagonal local solver '
+    if (present(prefix)) then
+      prefix_ = prefix
+    else
+      prefix_ = ""
+    end if
+
+    write(iout_,*) trim(prefix_), '  Diagonal local solver '
 
     return
 
@@ -352,7 +359,7 @@ module amg_s_l1_diag_solver
   
 contains
 
-  subroutine s_l1_diag_solver_descr(sv,info,iout,coarse)
+  subroutine s_l1_diag_solver_descr(sv,info,iout,coarse,prefix)
 
     Implicit None
 
@@ -361,11 +368,13 @@ contains
     integer(psb_ipk_), intent(out)              :: info
     integer(psb_ipk_), intent(in), optional     :: iout
     logical, intent(in), optional               :: coarse
+    character(len=*), intent(in), optional  :: prefix
 
     ! Local variables
     integer(psb_ipk_)      :: err_act
     character(len=20), parameter :: name='amg_s_l1_diag_solver_descr'
     integer(psb_ipk_) :: iout_
+    character(1024)    :: prefix_
 
     info = psb_success_
     if (present(iout)) then 
@@ -373,8 +382,13 @@ contains
     else
       iout_ = psb_out_unit
     endif
-    
-    write(iout_,*) '  L1 Diagonal solver '
+    if (present(prefix)) then
+      prefix_ = prefix
+    else
+      prefix_ = ""
+    end if
+
+    write(iout_,*) trim(prefix_), '  L1 Diagonal solver '
 
     return
 

amgprec/amg_s_gs_solver.f90

@@ -433,20 +433,22 @@ contains
     return
   end subroutine s_gs_solver_free
 
-  subroutine s_gs_solver_descr(sv,info,iout,coarse)
+  subroutine s_gs_solver_descr(sv,info,iout,coarse,prefix)
 
     Implicit None
 
     ! Arguments
     class(amg_s_gs_solver_type), intent(in) :: sv
-    integer(psb_ipk_), intent(out)             :: info
-    integer(psb_ipk_), intent(in), optional    :: iout
-    logical, intent(in), optional       :: coarse
+    integer(psb_ipk_), intent(out)            :: info
+    integer(psb_ipk_), intent(in), optional   :: iout
+    logical, intent(in), optional             :: coarse
+    character(len=*), intent(in), optional    :: prefix
 
     ! Local variables
     integer(psb_ipk_)      :: err_act
     character(len=20), parameter :: name='amg_s_gs_solver_descr'
     integer(psb_ipk_) :: iout_
+    character(1024)    :: prefix_
 
     call psb_erractionsave(err_act)
     info = psb_success_
@@ -455,12 +457,17 @@ contains
     else
       iout_ = psb_out_unit
     endif
+    if (present(prefix)) then
+      prefix_ = prefix
+    else
+      prefix_ = ""
+    end if
 
     if (sv%eps<=dzero) then 
-      write(iout_,*) '  Forward Gauss-Seidel iterative solver with  ',&
+      write(iout_,*) trim(prefix_), '  Forward Gauss-Seidel iterative solver with  ',&
            &  sv%sweeps,' sweeps'
     else
-      write(iout_,*) '  Forward Gauss-Seidel iterative solver with  tolerance',&
+      write(iout_,*) trim(prefix_), '  Forward Gauss-Seidel iterative solver with  tolerance',&
            &  sv%eps,' and maxit', sv%sweeps
     end if
     
@@ -526,20 +533,22 @@ contains
     val = .true.
   end function s_gs_solver_is_iterative
     
-  subroutine s_bwgs_solver_descr(sv,info,iout,coarse)
+  subroutine s_bwgs_solver_descr(sv,info,iout,coarse,prefix)
 
     Implicit None
 
     ! Arguments
     class(amg_s_bwgs_solver_type), intent(in) :: sv
-    integer(psb_ipk_), intent(out)             :: info
-    integer(psb_ipk_), intent(in), optional    :: iout
-    logical, intent(in), optional       :: coarse
+    integer(psb_ipk_), intent(out)              :: info
+    integer(psb_ipk_), intent(in), optional     :: iout
+    logical, intent(in), optional               :: coarse
+    character(len=*), intent(in), optional      :: prefix
 
     ! Local variables
     integer(psb_ipk_)      :: err_act
     character(len=20), parameter :: name='amg_s_bwgs_solver_descr'
     integer(psb_ipk_) :: iout_
+    character(1024)    :: prefix_
 
     call psb_erractionsave(err_act)
     info = psb_success_
@@ -548,12 +557,17 @@ contains
     else
       iout_ = psb_out_unit
     endif
+    if (present(prefix)) then
+      prefix_ = prefix
+    else
+      prefix_ = ""
+    end if
 
     if (sv%eps<=dzero) then 
-      write(iout_,*) '  Backward Gauss-Seidel iterative solver with  ',&
+      write(iout_,*) trim(prefix_), '  Backward Gauss-Seidel iterative solver with  ',&
            &  sv%sweeps,' sweeps'
     else
-      write(iout_,*) '  Backward Gauss-Seidel iterative solver with  tolerance',&
+      write(iout_,*) trim(prefix_), '  Backward Gauss-Seidel iterative solver with  tolerance',&
            &  sv%eps,' and maxit', sv%sweeps
     end if
     

amgprec/amg_s_id_solver.f90

@@ -157,7 +157,7 @@ contains
     return
   end subroutine s_id_solver_free
 
-  subroutine s_id_solver_descr(sv,info,iout,coarse)
+  subroutine s_id_solver_descr(sv,info,iout,coarse,prefix)
 
     Implicit None
 
@@ -165,12 +165,14 @@ contains
     class(amg_s_id_solver_type), intent(in) :: sv
     integer(psb_ipk_), intent(out)            :: info
     integer(psb_ipk_), intent(in), optional   :: iout
-    logical, intent(in), optional       :: coarse
+    logical, intent(in), optional             :: coarse
+    character(len=*), intent(in), optional    :: prefix
 
     ! Local variables
     integer(psb_ipk_)      :: err_act
     character(len=20), parameter :: name='amg_s_id_solver_descr'
     integer(psb_ipk_) :: iout_
+    character(1024)    :: prefix_
 
     info = psb_success_
     if (present(iout)) then 
@@ -178,8 +180,13 @@ contains
     else
       iout_ = psb_out_unit
     endif
+    if (present(prefix)) then
+      prefix_ = prefix
+    else
+      prefix_ = ""
+    end if
     
-    write(iout_,*) '  Identity local solver '
+    write(iout_,*) trim(prefix_), '  Identity local solver '
 
     return
 

amgprec/amg_s_ilu_solver.f90

@@ -406,7 +406,7 @@ contains
     return
   end subroutine s_ilu_solver_free
 
-  subroutine s_ilu_solver_descr(sv,info,iout,coarse)
+  subroutine s_ilu_solver_descr(sv,info,iout,coarse,prefix)
 
     Implicit None
 
@@ -414,12 +414,14 @@ contains
     class(amg_s_ilu_solver_type), intent(in) :: sv
     integer(psb_ipk_), intent(out)             :: info
     integer(psb_ipk_), intent(in), optional    :: iout
-    logical, intent(in), optional       :: coarse
-
+    logical, intent(in), optional              :: coarse
+    character(len=*), intent(in), optional  :: prefix
+      
     ! Local variables
     integer(psb_ipk_) :: err_act
     character(len=20), parameter :: name='amg_s_ilu_solver_descr'
     integer(psb_ipk_) :: iout_
+    character(1024)    :: prefix_
 
     call psb_erractionsave(err_act)
     info = psb_success_
@@ -428,15 +430,20 @@ contains
     else
       iout_ = psb_out_unit
     endif
+    if (present(prefix)) then
+      prefix_ = prefix
+    else
+      prefix_ = ""
+    end if
 
-    write(iout_,*) '  Incomplete factorization solver: ',&
+    write(iout_,*) trim(prefix_), '  Incomplete factorization solver: ',&
          &  amg_fact_names(sv%fact_type)
     select case(sv%fact_type)
     case(psb_ilu_n_,psb_milu_n_)      
-      write(iout_,*) '  Fill level:',sv%fill_in
+      write(iout_,*) trim(prefix_), '  Fill level:',sv%fill_in
     case(psb_ilu_t_)         
-      write(iout_,*) '  Fill level:',sv%fill_in
-      write(iout_,*) '  Fill threshold :',sv%thresh
+      write(iout_,*) trim(prefix_), '  Fill level:',sv%fill_in
+      write(iout_,*) trim(prefix_), '  Fill threshold :',sv%thresh
     end select
 
     call psb_erractionrestore(err_act)

amgprec/amg_s_invk_solver.f90

@@ -123,7 +123,7 @@ module amg_s_invk_solver
   end interface
 
   interface
-    subroutine amg_s_invk_solver_descr(sv,info,iout,coarse)
+    subroutine amg_s_invk_solver_descr(sv,info,iout,coarse,prefix)
       import :: psb_spk_, amg_s_invk_solver_type, psb_ipk_
 
       Implicit None
@@ -133,7 +133,7 @@ module amg_s_invk_solver
       integer(psb_ipk_), intent(out)            :: info
       integer(psb_ipk_), intent(in), optional   :: iout
       logical, intent(in), optional             :: coarse
-
+      character(len=*), intent(in), optional  :: prefix
     end subroutine amg_s_invk_solver_descr
   end interface
 

amgprec/amg_s_invt_solver.f90

@@ -134,16 +134,17 @@ module amg_s_invt_solver
   end interface
 
   interface
-    subroutine amg_s_invt_solver_descr(sv,info,iout,coarse)
+    subroutine amg_s_invt_solver_descr(sv,info,iout,coarse,prefix)
       import :: psb_spk_, amg_s_invt_solver_type, psb_ipk_
 
       Implicit None
 
       ! Arguments
       class(amg_s_invt_solver_type), intent(in) :: sv
-      integer(psb_ipk_), intent(out)            :: info
-      integer(psb_ipk_), intent(in), optional   :: iout
-      logical, intent(in), optional             :: coarse
+      integer(psb_ipk_), intent(out)              :: info
+      integer(psb_ipk_), intent(in), optional     :: iout
+      logical, intent(in), optional               :: coarse
+      character(len=*), intent(in), optional      :: prefix
 
     end subroutine amg_s_invt_solver_descr
   end interface

amgprec/amg_s_jac_smoother.f90

@@ -219,12 +219,13 @@ module amg_s_jac_smoother
   end interface
 
   interface
-    subroutine amg_s_jac_smoother_descr(sm,info,iout,coarse)
+    subroutine amg_s_jac_smoother_descr(sm,info,iout,coarse,prefix)
       import :: amg_s_jac_smoother_type, psb_ipk_
       class(amg_s_jac_smoother_type), intent(in) :: sm
       integer(psb_ipk_), intent(out)               :: info
       integer(psb_ipk_), intent(in), optional      :: iout
-      logical, intent(in), optional              :: coarse
+      logical, intent(in), optional                :: coarse
+      character(len=*), intent(in), optional       :: prefix
     end subroutine amg_s_jac_smoother_descr
   end interface
 
@@ -313,12 +314,13 @@ module amg_s_jac_smoother
   end interface
 
   interface
-    subroutine amg_s_l1_jac_smoother_descr(sm,info,iout,coarse)
+    subroutine amg_s_l1_jac_smoother_descr(sm,info,iout,coarse,prefix)
       import :: amg_s_l1_jac_smoother_type, psb_ipk_
       class(amg_s_l1_jac_smoother_type), intent(in) :: sm
-      integer(psb_ipk_), intent(out)               :: info
-      integer(psb_ipk_), intent(in), optional      :: iout
-      logical, intent(in), optional              :: coarse
+      integer(psb_ipk_), intent(out)                  :: info
+      integer(psb_ipk_), intent(in), optional         :: iout
+      logical, intent(in), optional                   :: coarse
+      character(len=*), intent(in), optional          :: prefix
     end subroutine amg_s_l1_jac_smoother_descr
   end interface
 

amgprec/amg_s_krm_solver.f90

@@ -436,7 +436,7 @@ contains
     val = "KRM solver"
   end function s_krm_solver_get_fmt
 
-  subroutine s_krm_solver_descr(sv,info,iout,coarse)
+  subroutine s_krm_solver_descr(sv,info,iout,coarse,prefix)
 
     Implicit None
 
@@ -444,12 +444,14 @@ contains
     class(amg_s_krm_solver_type), intent(in) :: sv
     integer(psb_ipk_), intent(out)             :: info
     integer(psb_ipk_), intent(in), optional    :: iout
-    logical, intent(in), optional       :: coarse
+    logical, intent(in), optional              :: coarse
+    character(len=*), intent(in), optional  :: prefix
 
     ! Local variables
     integer(psb_ipk_)      :: err_act
     character(len=20), parameter :: name='amg_s_krm_solver_descr'
     integer(psb_ipk_) :: iout_
+    character(1024)    :: prefix_
 
     call psb_erractionsave(err_act)
     info = psb_success_
@@ -458,23 +460,22 @@ contains
     else
       iout_ = psb_out_unit
     endif
+    if (present(prefix)) then
+      prefix_ = prefix
+    else
+      prefix_ = ""
+    end if
 
     if (sv%global) then
-      write(iout_,*) '  Krylov solver (global)'
+      write(iout_,*) trim(prefix_), '  Krylov solver (global)'
     else
-      write(iout_,*) '  Krylov solver (local) '
+      write(iout_,*) trim(prefix_), '  Krylov solver (local) '
     end if
-    write(iout_,*) '    method: ',sv%method
-    write(iout_,*) '     kprec: ',sv%kprec
-    if (sv%i_sub_solve > 0) then
-      write(iout_,*) ' sub_solve: ',amg_fact_names(sv%i_sub_solve)
-    else
-      write(iout_,*) ' sub_solve: ',sv%sub_solve
-    end if
-    write(iout_,*) '     itmax: ',sv%itmax
-    write(iout_,*) '       eps: ',sv%eps
-    write(iout_,*) '    fillin: ',sv%fillin
-
+    write(iout_,*) trim(prefix_), '    method: ',sv%method
+    write(iout_,*) trim(prefix_), '     kprec: ',sv%kprec
+    call sv%prec%descr(info,iout_,prefix='KRM : '//prefix_)
+    write(iout_,*) trim(prefix_), '     itmax: ',sv%itmax
+    write(iout_,*) trim(prefix_), '       eps: ',sv%eps
 
     call psb_erractionrestore(err_act)
     return

amgprec/amg_s_matchboxp_mod.f90

@@ -143,9 +143,10 @@ contains
     type(psb_ls_coo_sparse_mat) :: tmpcoo
     logical  :: display_out_, print_out_, reproducible_
     logical, parameter  :: dump=.false., debug=.false., dump_mate=.false., &
-         & debug_ilaggr=.false., debug_sync=.false.
+         & debug_ilaggr=.false., debug_sync=.false., debug_mate=.false.
     integer(psb_ipk_), save :: idx_bldmtc=-1, idx_phase1=-1, idx_phase2=-1, idx_phase3=-1
     logical, parameter :: do_timings=.true.
+    integer, parameter :: ilaggr_neginit=-1, ilaggr_nonlocal=-2
 
     ictxt = desc_a%get_ctxt()
     call psb_info(ictxt,iam,np)
@@ -187,7 +188,7 @@ contains
     call desc_a%l2gip(ilv,info,owned=.false.)
 
     call psb_geall(ilaggr,desc_a,info)
-    ilaggr = -1
+    ilaggr = ilaggr_neginit
     call psb_geasb(ilaggr,desc_a,info)
     nr = a%get_nrows()
     nc = a%get_ncols()
@@ -213,7 +214,20 @@ contains
       call psb_barrier(ictxt)
       if (iam == 0) write(0,*)' out from buildmatching:', info
     end if
-
+    if (debug_mate) then
+      block
+        integer(psb_lpk_), allocatable :: ckmate(:)
+        allocate(ckmate(nr))
+        ckmate(1:nr) = mate(1:nr)
+        call psb_msort(ckmate(1:nr))
+        do i=1,nr-1
+          if ((ckmate(i)>0) .and. (ckmate(i) == ckmate(i+1))) then
+            write(0,*) iam,' Duplicate mate entry  at',i,' :',ckmate(i)
+          end if
+        end do
+      end block
+    end if
+     
     if (info == 0) then
       if (do_timings) call psb_tic(idx_phase2)
       if (debug_sync) then
@@ -259,7 +273,7 @@ contains
           cycle
         else
 
-          if (ilaggr(k) == -1) then
+          if (ilaggr(k) == ilaggr_neginit) then
 
             wk   = w(k)
             widx = w(idx)
@@ -267,7 +281,7 @@ contains
             nrmagg = wmax*sqrt((wk/wmax)**2+(widx/wmax)**2)
             if (nrmagg > epsilon(nrmagg)) then
               if (idx <= nr) then
-                if (ilaggr(idx) == -1) then
+                if (ilaggr(idx) == ilaggr_neginit) then
                   ! Now, if both vertices are local, the aggregate is local
                   ! (kinda obvious).
                   nlaggr(iam) = nlaggr(iam) + 1
@@ -275,6 +289,9 @@ contains
                   ilaggr(idx) = nlaggr(iam)
                   wtemp(k)    = w(k)/nrmagg
                   wtemp(idx)  = w(idx)/nrmagg
+                else
+                  write(0,*) iam,' Inconsistent mate? ',k,mate(k),idx,&
+                       &mate(idx),ilaggr(idx)
                 end if
                 nlpairs = nlpairs+1
               else if (idx <= nc) then
@@ -294,7 +311,7 @@ contains
                     ilaggr(k)   = nlaggr(iam)
                     nlpairs = nlpairs+1
                   else
-                    ilaggr(k) = -2
+                    ilaggr(k) = ilaggr_nonlocal
                   end if
                 else
                   ! Use a statistically unbiased tie-breaking rule,
@@ -309,7 +326,7 @@ contains
                     ilaggr(k)   = nlaggr(iam)
                     nlpairs = nlpairs+1
                   else
-                    ilaggr(k) = -2
+                    ilaggr(k) = ilaggr_nonlocal
                   end if
                 end if
               end if
@@ -325,6 +342,12 @@ contains
                 nlsingl     = nlsingl + 1
               end if
             end if
+            if (ilaggr(k) == ilaggr_neginit) then
+              write(0,*) iam,' Error: no update to ',k,mate(k),&
+                   & abs(w(k)),nrmagg,epsilon(nrmagg),wtemp(k)
+            end if
+          else
+            if (ilaggr(k)<0) write(0,*) 'Strange? ',k,ilaggr(k)
           end if
         end if
       end do
@@ -332,7 +355,7 @@ contains
       if (do_timings) call psb_tic(idx_phase3)
 
       ! Ok, now compute offsets, gather halo and fix non-local
-      ! aggregates  (those where ilaggr == -2)
+      ! aggregates  (those where ilaggr == ilaggr_nonlocal)
       call psb_sum(ictxt,nlaggr)
       ntaggr  = sum(nlaggr(0:np-1))
       naggrm1 = sum(nlaggr(0:iam-1))
@@ -347,7 +370,7 @@ contains
       call psb_halo(wtemp,desc_a,info)
       ! Cleanup as yet unmarked entries
       do k=1,nr
-        if (ilaggr(k) == -2) then
+        if (ilaggr(k) == ilaggr_nonlocal) then
           idx = mate(k)
           if (idx > nr) then
             i   = ilaggr(idx)
@@ -359,9 +382,14 @@ contains
           else
             write(0,*) 'Error : unresolved (paired) index ',k,idx,i,nr,nc, ilv(k),ilv(idx)
           end if
-        end if
-        if (ilaggr(k) <0) then
-          write(0,*) 'Matchboxp: Funny number: ',k,ilv(k),ilaggr(k),wtemp(k)
+        else if (ilaggr(k) <0) then
+          write(0,*) iam,'Matchboxp: Funny number: ',k,ilv(k),ilaggr(k),wtemp(k)
+          write(0,*) iam,'         :             : ',nr,nc,mate(k)
+          if (mate(k) <= nr) then
+            write(0,*) iam,'         :             : ',ilaggr(mate(k)),mate(mate(k)),&
+                 & ilv(k),ilv(mate(k)), ilv(mate(mate(k))),ilaggr(mate(mate(k)))
+          end if
+          flush(0)
         end if
       end do
       if (debug_sync) then
@@ -414,7 +442,7 @@ contains
 
       end block
       if (iam == 0) then
-        write(0,*) 'Matching statistics: Unmatched nodes ',&
+        write(0,*) iam,'Matching statistics: Unmatched nodes ',&
              & nunmatched,' Singletons:',nlsingl,' Pairs:',nlpairs
       end if
 

amgprec/amg_s_mumps_solver.F90

@@ -313,22 +313,24 @@ subroutine s_mumps_solver_finalize(sv)
 
 end subroutine s_mumps_solver_finalize
 
-subroutine s_mumps_solver_descr(sv,info,iout,coarse)
+subroutine s_mumps_solver_descr(sv,info,iout,coarse,prefix)
 
   Implicit None
 
   ! Arguments
   class(amg_s_mumps_solver_type), intent(in) :: sv
-  integer(psb_ipk_), intent(out)                     :: info
-  integer(psb_ipk_), intent(in), optional            :: iout
-  logical, intent(in), optional       :: coarse
-
+  integer(psb_ipk_), intent(out)               :: info
+  integer(psb_ipk_), intent(in), optional      :: iout
+  logical, intent(in), optional                :: coarse
+  character(len=*), intent(in), optional       :: prefix
+      
   ! Local variables
   integer(psb_ipk_)   :: err_act
   type(psb_ctxt_type) :: ctxt
   integer(psb_ipk_)   :: me, np
   character(len=20), parameter :: name='amg_z_mumps_solver_descr'
   integer(psb_ipk_) :: iout_
+  character(1024)    :: prefix_
 
   call psb_erractionsave(err_act)
   info = psb_success_
@@ -337,8 +339,13 @@ subroutine s_mumps_solver_descr(sv,info,iout,coarse)
   else
     iout_ = psb_out_unit
   endif
+  if (present(prefix)) then
+    prefix_ = prefix
+  else
+    prefix_ = ""
+  end if
 
-  write(iout_,*) '  MUMPS  Solver. '
+  write(iout_,*) trim(prefix_), '  MUMPS  Solver. '
 
   call psb_erractionrestore(err_act)
   return

amgprec/amg_s_onelev_mod.f90

@@ -258,7 +258,7 @@ module amg_s_onelev_mod
   end interface
 
   interface
-    subroutine amg_s_base_onelev_descr(lv,il,nl,ilmin,info,iout, verbosity)
+    subroutine amg_s_base_onelev_descr(lv,il,nl,ilmin,info,iout, verbosity,prefix)
       import :: psb_sspmat_type, psb_s_vect_type, psb_s_base_vect_type, &
            & psb_slinmap_type, psb_spk_, amg_s_onelev_type, &
            & psb_ipk_, psb_epk_, psb_desc_type
@@ -269,6 +269,7 @@ module amg_s_onelev_mod
       integer(psb_ipk_), intent(out)          :: info
       integer(psb_ipk_), intent(in), optional :: iout
       integer(psb_ipk_), intent(in), optional :: verbosity
+      character(len=*), intent(in), optional  :: prefix
     end subroutine amg_s_base_onelev_descr
   end interface
 

amgprec/amg_s_parmatch_aggregator_mod.F90

@@ -390,18 +390,25 @@ contains
 
   end function amg_s_parmatch_aggregator_sizeof
 
-  subroutine  amg_s_parmatch_aggregator_descr(ag,parms,iout,info)
+  subroutine  amg_s_parmatch_aggregator_descr(ag,parms,iout,info,prefix)
     implicit none
     class(amg_s_parmatch_aggregator_type), intent(in) :: ag
     type(amg_sml_parms), intent(in)   :: parms
     integer(psb_ipk_), intent(in)  :: iout
     integer(psb_ipk_), intent(out) :: info
+    character(len=*), intent(in), optional  :: prefix
+    character(1024)    :: prefix_
+    if (present(prefix)) then
+      prefix_ = prefix
+    else
+      prefix_ = ""
+    end if
 
-    write(iout,*) 'Parallel Matching Aggregator'
-    write(iout,*) '   Number of matching  sweeps: ',ag%n_sweeps
-    write(iout,*) '   Matching algorithm         : MatchBoxP (PREIS)'
-    write(iout,*) 'Aggregator object type: ',ag%fmt()
-    call parms%mldescr(iout,info)
+    write(iout,*) trim(prefix_),' ','Parallel Matching Aggregator'
+    write(iout,*) trim(prefix_),' ','   Number of matching  sweeps: ',ag%n_sweeps
+    write(iout,*) trim(prefix_),' ','   Matching algorithm         : MatchBoxP (PREIS)'
+    write(iout,*) trim(prefix_),' ','Aggregator object type: ',ag%fmt()
+    call parms%mldescr(iout,info,prefix=prefix)
 
     return
   end subroutine amg_s_parmatch_aggregator_descr

amgprec/amg_s_prec_type.f90

@@ -155,15 +155,16 @@ module amg_s_prec_type
 
 
   interface amg_precdescr
-    subroutine amg_sfile_prec_descr(prec,info,iout,root,verbosity)
+    subroutine amg_sfile_prec_descr(prec,info,iout,root,verbosity,prefix)
       import :: amg_sprec_type, psb_ipk_
       implicit none
       ! Arguments
-      class(amg_sprec_type), intent(in)      :: prec
-      integer(psb_ipk_), intent(out) :: info
+      class(amg_sprec_type), intent(in)     :: prec
+      integer(psb_ipk_), intent(out)          :: info
       integer(psb_ipk_), intent(in), optional :: iout
       integer(psb_ipk_), intent(in), optional :: root
       integer(psb_ipk_), intent(in), optional :: verbosity
+      character(len=*), intent(in), optional  :: prefix
     end subroutine amg_sfile_prec_descr
   end interface
 

amgprec/amg_s_slu_solver.F90

@@ -385,20 +385,22 @@ contains
   
   end subroutine s_slu_solver_finalize
 
-  subroutine s_slu_solver_descr(sv,info,iout,coarse)
+  subroutine s_slu_solver_descr(sv,info,iout,coarse,prefix)
 
     Implicit None
 
     ! Arguments
     class(amg_s_slu_solver_type), intent(in) :: sv
-    integer, intent(out)                     :: info
-    integer, intent(in), optional            :: iout
-    logical, intent(in), optional       :: coarse
+    integer, intent(out)                       :: info
+    integer, intent(in), optional              :: iout
+    logical, intent(in), optional              :: coarse
+    character(len=*), intent(in), optional     :: prefix
 
     ! Local variables
     integer      :: err_act
     character(len=20), parameter :: name='amg_s_slu_solver_descr'
     integer :: iout_
+    character(1024)    :: prefix_
 
     call psb_erractionsave(err_act)
     info = psb_success_
@@ -407,8 +409,13 @@ contains
     else
       iout_ = psb_out_unit
     endif
+    if (present(prefix)) then
+      prefix_ = prefix
+    else
+      prefix_ = ""
+    end if
     
-    write(iout_,*) '  SuperLU Sparse Factorization Solver. '
+    write(iout_,*) trim(prefix_), '  SuperLU Sparse Factorization Solver. '
 
     call psb_erractionrestore(err_act)
     return

amgprec/amg_s_symdec_aggregator_mod.f90

@@ -88,16 +88,25 @@ contains
     val = "Symmetric Decoupled aggregation"
   end function amg_s_symdec_aggregator_fmt
   
-  subroutine  amg_s_symdec_aggregator_descr(ag,parms,iout,info)
+  subroutine  amg_s_symdec_aggregator_descr(ag,parms,iout,info,prefix)
     implicit none 
     class(amg_s_symdec_aggregator_type), intent(in) :: ag
     type(amg_sml_parms), intent(in)   :: parms
     integer(psb_ipk_), intent(in)  :: iout
     integer(psb_ipk_), intent(out) :: info
+    character(len=*), intent(in), optional  :: prefix
 
-    write(iout,*) 'Decoupled Aggregator locally-symmetrized'
-    write(iout,*) 'Aggregator object type: ',ag%fmt()
-    call parms%mldescr(iout,info)
+    character(1024)    :: prefix_
+
+    if (present(prefix)) then
+      prefix_ = prefix
+    else
+      prefix_ = ""
+    end if
+
+    write(iout,*) trim(prefix_),' ','Decoupled Aggregator locally-symmetrized'
+    write(iout,*) trim(prefix_),' ','Aggregator object type: ',ag%fmt()
+    call parms%mldescr(iout,info,prefix=prefix)
     
     return
   end subroutine amg_s_symdec_aggregator_descr

amgprec/amg_z_ainv_solver.F90

@@ -198,7 +198,7 @@ module amg_z_ainv_solver
 !!$  end interface
 
   interface
-    subroutine amg_z_ainv_solver_descr(sv,info,iout,coarse)
+    subroutine amg_z_ainv_solver_descr(sv,info,iout,coarse,prefix)
       import :: psb_dpk_, amg_z_ainv_solver_type, psb_ipk_
 
       Implicit None
@@ -208,7 +208,7 @@ module amg_z_ainv_solver
       integer(psb_ipk_), intent(out)            :: info
       integer(psb_ipk_), intent(in), optional   :: iout
       logical, intent(in), optional             :: coarse
-
+      character(len=*), intent(in), optional  :: prefix
     end subroutine amg_z_ainv_solver_descr
   end interface
 

amgprec/amg_z_as_smoother.f90

@@ -396,21 +396,23 @@ contains
   end subroutine z_as_smoother_default
 
 
-  subroutine z_as_smoother_descr(sm,info,iout,coarse)
+  subroutine z_as_smoother_descr(sm,info,iout,coarse,prefix)
 
     Implicit None
 
     ! Arguments
     class(amg_z_as_smoother_type), intent(in) :: sm
-    integer(psb_ipk_), intent(out)                      :: info
-    integer(psb_ipk_), intent(in), optional             :: iout
-    logical, intent(in), optional             :: coarse
+    integer(psb_ipk_), intent(out)              :: info
+    integer(psb_ipk_), intent(in), optional     :: iout
+    logical, intent(in), optional              :: coarse
+      character(len=*), intent(in), optional   :: prefix
 
     ! Local variables
     integer(psb_ipk_)      :: err_act
     character(len=20), parameter :: name='amg_z_as_smoother_descr'
     integer(psb_ipk_) :: iout_
     logical      :: coarse_
+    character(1024)    :: prefix_
 
     call psb_erractionsave(err_act)
     info = psb_success_
@@ -424,16 +426,21 @@ contains
     else
       iout_ = psb_out_unit
     endif
+    if (present(prefix)) then
+      prefix_ = prefix
+    else
+      prefix_ = ""
+    end if
 
     if (.not.coarse_) then 
-      write(iout_,*) '  Additive Schwarz with  ',&
+      write(iout_,*) trim(prefix_), '  Additive Schwarz with  ',&
            &  sm%novr, ' overlap layers.'
-      write(iout_,*) '  Restrictor:  ',restrict_names(sm%restr)
-      write(iout_,*) '  Prolongator: ',prolong_names(sm%prol)
-      write(iout_,*) '  Local solver:'
+      write(iout_,*) trim(prefix_), '  Restrictor:  ',restrict_names(sm%restr)
+      write(iout_,*) trim(prefix_), '  Prolongator: ',prolong_names(sm%prol)
+      write(iout_,*) trim(prefix_), '  Local solver:'
     endif
     if (allocated(sm%sv)) then 
-      call sm%sv%descr(info,iout_,coarse=coarse)
+      call sm%sv%descr(info,iout_,coarse=coarse,prefix=prefix)
     end if
 
     call psb_erractionrestore(err_act)

amgprec/amg_z_base_aggregator_mod.f90

@@ -275,15 +275,22 @@ contains
     val = .false.
   end function amg_z_base_aggregator_xt_desc
 
-  subroutine  amg_z_base_aggregator_descr(ag,parms,iout,info)
+  subroutine  amg_z_base_aggregator_descr(ag,parms,iout,info,prefix)
     implicit none 
     class(amg_z_base_aggregator_type), intent(in) :: ag
     type(amg_dml_parms), intent(in)   :: parms
     integer(psb_ipk_), intent(in)  :: iout
     integer(psb_ipk_), intent(out) :: info
+    character(len=*), intent(in), optional  :: prefix
+    character(1024)    :: prefix_
+    if (present(prefix)) then
+      prefix_ = prefix
+    else
+      prefix_ = ""
+    end if
 
-    write(iout,*) 'Aggregator object type: ',ag%fmt()
-    call parms%mldescr(iout,info)
+    write(iout,*) trim(prefix_),' ', 'Aggregator object type: ',ag%fmt()
+    call parms%mldescr(iout,info,prefix=prefix)
     
     return
   end subroutine amg_z_base_aggregator_descr

amgprec/amg_z_base_smoother_mod.f90

@@ -272,7 +272,7 @@ module amg_z_base_smoother_mod
   end interface
   
   interface 
-    subroutine amg_z_base_smoother_descr(sm,info,iout,coarse)
+    subroutine amg_z_base_smoother_descr(sm,info,iout,coarse,prefix)
       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_base_smoother_type, psb_ipk_
@@ -281,6 +281,7 @@ module amg_z_base_smoother_mod
       integer(psb_ipk_), intent(out)                :: info
       integer(psb_ipk_), intent(in), optional       :: iout
       logical, intent(in), optional                 :: coarse
+      character(len=*), intent(in), optional        :: prefix
     end subroutine amg_z_base_smoother_descr
   end interface
   

amgprec/amg_z_base_solver_mod.f90

@@ -270,7 +270,7 @@ module amg_z_base_solver_mod
   end interface
   
   interface
-    subroutine amg_z_base_solver_descr(sv,info,iout,coarse)
+    subroutine amg_z_base_solver_descr(sv,info,iout,coarse,prefix)
       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_base_solver_type, psb_ipk_
@@ -281,7 +281,7 @@ module amg_z_base_solver_mod
       integer(psb_ipk_), intent(out)              :: info
       integer(psb_ipk_), intent(in), optional     :: iout
       logical, intent(in), optional               :: coarse
-
+      character(len=*), intent(in), optional      :: prefix
     end subroutine amg_z_base_solver_descr
   end interface 
   

amgprec/amg_z_dec_aggregator_mod.f90

@@ -184,16 +184,23 @@ contains
     val = "Decoupled aggregation"
   end function amg_z_dec_aggregator_fmt
   
-  subroutine  amg_z_dec_aggregator_descr(ag,parms,iout,info)
+  subroutine  amg_z_dec_aggregator_descr(ag,parms,iout,info,prefix)
     implicit none 
     class(amg_z_dec_aggregator_type), intent(in) :: ag
     type(amg_dml_parms), intent(in)   :: parms
     integer(psb_ipk_), intent(in)  :: iout
     integer(psb_ipk_), intent(out) :: info
+    character(len=*), intent(in), optional  :: prefix
+    character(1024)    :: prefix_
+    if (present(prefix)) then
+      prefix_ = prefix
+    else
+      prefix_ = ""
+    end if
 
-    write(iout,*) 'Decoupled Aggregator'
-    write(iout,*) 'Aggregator object type: ',ag%fmt()
-    call parms%mldescr(iout,info)
+    write(iout,*) trim(prefix_),' ','Decoupled Aggregator'
+    write(iout,*) trim(prefix_),' ','Aggregator object type: ',ag%fmt()
+    call parms%mldescr(iout,info,prefix=prefix)
     
     return
   end subroutine amg_z_dec_aggregator_descr

amgprec/amg_z_diag_solver.f90

@@ -219,7 +219,7 @@ contains
 
   end subroutine z_diag_solver_free
 
-  subroutine z_diag_solver_descr(sv,info,iout,coarse)
+  subroutine z_diag_solver_descr(sv,info,iout,coarse,prefix)
 
     Implicit None
 
@@ -228,11 +228,13 @@ contains
     integer(psb_ipk_), intent(out)              :: info
     integer(psb_ipk_), intent(in), optional     :: iout
     logical, intent(in), optional               :: coarse
+    character(len=*), intent(in), optional      :: prefix
 
     ! Local variables
     integer(psb_ipk_)      :: err_act
     character(len=20), parameter :: name='amg_z_diag_solver_descr'
     integer(psb_ipk_) :: iout_
+    character(1024)    :: prefix_
 
     info = psb_success_
     if (present(iout)) then 
@@ -240,8 +242,13 @@ contains
     else
       iout_ = psb_out_unit
     endif
-    
-    write(iout_,*) '  Diagonal local solver '
+    if (present(prefix)) then
+      prefix_ = prefix
+    else
+      prefix_ = ""
+    end if
+
+    write(iout_,*) trim(prefix_), '  Diagonal local solver '
 
     return
 
@@ -352,7 +359,7 @@ module amg_z_l1_diag_solver
   
 contains
 
-  subroutine z_l1_diag_solver_descr(sv,info,iout,coarse)
+  subroutine z_l1_diag_solver_descr(sv,info,iout,coarse,prefix)
 
     Implicit None
 
@@ -361,11 +368,13 @@ contains
     integer(psb_ipk_), intent(out)              :: info
     integer(psb_ipk_), intent(in), optional     :: iout
     logical, intent(in), optional               :: coarse
+    character(len=*), intent(in), optional  :: prefix
 
     ! Local variables
     integer(psb_ipk_)      :: err_act
     character(len=20), parameter :: name='amg_z_l1_diag_solver_descr'
     integer(psb_ipk_) :: iout_
+    character(1024)    :: prefix_
 
     info = psb_success_
     if (present(iout)) then 
@@ -373,8 +382,13 @@ contains
     else
       iout_ = psb_out_unit
     endif
-    
-    write(iout_,*) '  L1 Diagonal solver '
+    if (present(prefix)) then
+      prefix_ = prefix
+    else
+      prefix_ = ""
+    end if
+
+    write(iout_,*) trim(prefix_), '  L1 Diagonal solver '
 
     return
 

amgprec/amg_z_gs_solver.f90

@@ -433,20 +433,22 @@ contains
     return
   end subroutine z_gs_solver_free
 
-  subroutine z_gs_solver_descr(sv,info,iout,coarse)
+  subroutine z_gs_solver_descr(sv,info,iout,coarse,prefix)
 
     Implicit None
 
     ! Arguments
     class(amg_z_gs_solver_type), intent(in) :: sv
-    integer(psb_ipk_), intent(out)             :: info
-    integer(psb_ipk_), intent(in), optional    :: iout
-    logical, intent(in), optional       :: coarse
+    integer(psb_ipk_), intent(out)            :: info
+    integer(psb_ipk_), intent(in), optional   :: iout
+    logical, intent(in), optional             :: coarse
+    character(len=*), intent(in), optional    :: prefix
 
     ! Local variables
     integer(psb_ipk_)      :: err_act
     character(len=20), parameter :: name='amg_z_gs_solver_descr'
     integer(psb_ipk_) :: iout_
+    character(1024)    :: prefix_
 
     call psb_erractionsave(err_act)
     info = psb_success_
@@ -455,12 +457,17 @@ contains
     else
       iout_ = psb_out_unit
     endif
+    if (present(prefix)) then
+      prefix_ = prefix
+    else
+      prefix_ = ""
+    end if
 
     if (sv%eps<=dzero) then 
-      write(iout_,*) '  Forward Gauss-Seidel iterative solver with  ',&
+      write(iout_,*) trim(prefix_), '  Forward Gauss-Seidel iterative solver with  ',&
            &  sv%sweeps,' sweeps'
     else
-      write(iout_,*) '  Forward Gauss-Seidel iterative solver with  tolerance',&
+      write(iout_,*) trim(prefix_), '  Forward Gauss-Seidel iterative solver with  tolerance',&
            &  sv%eps,' and maxit', sv%sweeps
     end if
     
@@ -526,20 +533,22 @@ contains
     val = .true.
   end function z_gs_solver_is_iterative
     
-  subroutine z_bwgs_solver_descr(sv,info,iout,coarse)
+  subroutine z_bwgs_solver_descr(sv,info,iout,coarse,prefix)
 
     Implicit None
 
     ! Arguments
     class(amg_z_bwgs_solver_type), intent(in) :: sv
-    integer(psb_ipk_), intent(out)             :: info
-    integer(psb_ipk_), intent(in), optional    :: iout
-    logical, intent(in), optional       :: coarse
+    integer(psb_ipk_), intent(out)              :: info
+    integer(psb_ipk_), intent(in), optional     :: iout
+    logical, intent(in), optional               :: coarse
+    character(len=*), intent(in), optional      :: prefix
 
     ! Local variables
     integer(psb_ipk_)      :: err_act
     character(len=20), parameter :: name='amg_z_bwgs_solver_descr'
     integer(psb_ipk_) :: iout_
+    character(1024)    :: prefix_
 
     call psb_erractionsave(err_act)
     info = psb_success_
@@ -548,12 +557,17 @@ contains
     else
       iout_ = psb_out_unit
     endif
+    if (present(prefix)) then
+      prefix_ = prefix
+    else
+      prefix_ = ""
+    end if
 
     if (sv%eps<=dzero) then 
-      write(iout_,*) '  Backward Gauss-Seidel iterative solver with  ',&
+      write(iout_,*) trim(prefix_), '  Backward Gauss-Seidel iterative solver with  ',&
            &  sv%sweeps,' sweeps'
     else
-      write(iout_,*) '  Backward Gauss-Seidel iterative solver with  tolerance',&
+      write(iout_,*) trim(prefix_), '  Backward Gauss-Seidel iterative solver with  tolerance',&
            &  sv%eps,' and maxit', sv%sweeps
     end if
     

amgprec/amg_z_id_solver.f90

@@ -157,7 +157,7 @@ contains
     return
   end subroutine z_id_solver_free
 
-  subroutine z_id_solver_descr(sv,info,iout,coarse)
+  subroutine z_id_solver_descr(sv,info,iout,coarse,prefix)
 
     Implicit None
 
@@ -165,12 +165,14 @@ contains
     class(amg_z_id_solver_type), intent(in) :: sv
     integer(psb_ipk_), intent(out)            :: info
     integer(psb_ipk_), intent(in), optional   :: iout
-    logical, intent(in), optional       :: coarse
+    logical, intent(in), optional             :: coarse
+    character(len=*), intent(in), optional    :: prefix
 
     ! Local variables
     integer(psb_ipk_)      :: err_act
     character(len=20), parameter :: name='amg_z_id_solver_descr'
     integer(psb_ipk_) :: iout_
+    character(1024)    :: prefix_
 
     info = psb_success_
     if (present(iout)) then 
@@ -178,8 +180,13 @@ contains
     else
       iout_ = psb_out_unit
     endif
+    if (present(prefix)) then
+      prefix_ = prefix
+    else
+      prefix_ = ""
+    end if
     
-    write(iout_,*) '  Identity local solver '
+    write(iout_,*) trim(prefix_), '  Identity local solver '
 
     return
 

amgprec/amg_z_ilu_solver.f90

@@ -406,7 +406,7 @@ contains
     return
   end subroutine z_ilu_solver_free
 
-  subroutine z_ilu_solver_descr(sv,info,iout,coarse)
+  subroutine z_ilu_solver_descr(sv,info,iout,coarse,prefix)
 
     Implicit None
 
@@ -414,12 +414,14 @@ contains
     class(amg_z_ilu_solver_type), intent(in) :: sv
     integer(psb_ipk_), intent(out)             :: info
     integer(psb_ipk_), intent(in), optional    :: iout
-    logical, intent(in), optional       :: coarse
-
+    logical, intent(in), optional              :: coarse
+    character(len=*), intent(in), optional  :: prefix
+      
     ! Local variables
     integer(psb_ipk_) :: err_act
     character(len=20), parameter :: name='amg_z_ilu_solver_descr'
     integer(psb_ipk_) :: iout_
+    character(1024)    :: prefix_
 
     call psb_erractionsave(err_act)
     info = psb_success_
@@ -428,15 +430,20 @@ contains
     else
       iout_ = psb_out_unit
     endif
+    if (present(prefix)) then
+      prefix_ = prefix
+    else
+      prefix_ = ""
+    end if
 
-    write(iout_,*) '  Incomplete factorization solver: ',&
+    write(iout_,*) trim(prefix_), '  Incomplete factorization solver: ',&
          &  amg_fact_names(sv%fact_type)
     select case(sv%fact_type)
     case(psb_ilu_n_,psb_milu_n_)      
-      write(iout_,*) '  Fill level:',sv%fill_in
+      write(iout_,*) trim(prefix_), '  Fill level:',sv%fill_in
     case(psb_ilu_t_)         
-      write(iout_,*) '  Fill level:',sv%fill_in
-      write(iout_,*) '  Fill threshold :',sv%thresh
+      write(iout_,*) trim(prefix_), '  Fill level:',sv%fill_in
+      write(iout_,*) trim(prefix_), '  Fill threshold :',sv%thresh
     end select
 
     call psb_erractionrestore(err_act)

amgprec/amg_z_invk_solver.f90

@@ -123,7 +123,7 @@ module amg_z_invk_solver
   end interface
 
   interface
-    subroutine amg_z_invk_solver_descr(sv,info,iout,coarse)
+    subroutine amg_z_invk_solver_descr(sv,info,iout,coarse,prefix)
       import :: psb_dpk_, amg_z_invk_solver_type, psb_ipk_
 
       Implicit None
@@ -133,7 +133,7 @@ module amg_z_invk_solver
       integer(psb_ipk_), intent(out)            :: info
       integer(psb_ipk_), intent(in), optional   :: iout
       logical, intent(in), optional             :: coarse
-
+      character(len=*), intent(in), optional  :: prefix
     end subroutine amg_z_invk_solver_descr
   end interface
 

amgprec/amg_z_invt_solver.f90

@@ -134,16 +134,17 @@ module amg_z_invt_solver
   end interface
 
   interface
-    subroutine amg_z_invt_solver_descr(sv,info,iout,coarse)
+    subroutine amg_z_invt_solver_descr(sv,info,iout,coarse,prefix)
       import :: psb_dpk_, amg_z_invt_solver_type, psb_ipk_
 
       Implicit None
 
       ! Arguments
       class(amg_z_invt_solver_type), intent(in) :: sv
-      integer(psb_ipk_), intent(out)            :: info
-      integer(psb_ipk_), intent(in), optional   :: iout
-      logical, intent(in), optional             :: coarse
+      integer(psb_ipk_), intent(out)              :: info
+      integer(psb_ipk_), intent(in), optional     :: iout
+      logical, intent(in), optional               :: coarse
+      character(len=*), intent(in), optional      :: prefix
 
     end subroutine amg_z_invt_solver_descr
   end interface

amgprec/amg_z_jac_smoother.f90

@@ -219,12 +219,13 @@ module amg_z_jac_smoother
   end interface
 
   interface
-    subroutine amg_z_jac_smoother_descr(sm,info,iout,coarse)
+    subroutine amg_z_jac_smoother_descr(sm,info,iout,coarse,prefix)
       import :: amg_z_jac_smoother_type, psb_ipk_
       class(amg_z_jac_smoother_type), intent(in) :: sm
       integer(psb_ipk_), intent(out)               :: info
       integer(psb_ipk_), intent(in), optional      :: iout
-      logical, intent(in), optional              :: coarse
+      logical, intent(in), optional                :: coarse
+      character(len=*), intent(in), optional       :: prefix
     end subroutine amg_z_jac_smoother_descr
   end interface
 
@@ -313,12 +314,13 @@ module amg_z_jac_smoother
   end interface
 
   interface
-    subroutine amg_z_l1_jac_smoother_descr(sm,info,iout,coarse)
+    subroutine amg_z_l1_jac_smoother_descr(sm,info,iout,coarse,prefix)
       import :: amg_z_l1_jac_smoother_type, psb_ipk_
       class(amg_z_l1_jac_smoother_type), intent(in) :: sm
-      integer(psb_ipk_), intent(out)               :: info
-      integer(psb_ipk_), intent(in), optional      :: iout
-      logical, intent(in), optional              :: coarse
+      integer(psb_ipk_), intent(out)                  :: info
+      integer(psb_ipk_), intent(in), optional         :: iout
+      logical, intent(in), optional                   :: coarse
+      character(len=*), intent(in), optional          :: prefix
     end subroutine amg_z_l1_jac_smoother_descr
   end interface
 

amgprec/amg_z_krm_solver.f90

@@ -436,7 +436,7 @@ contains
     val = "KRM solver"
   end function z_krm_solver_get_fmt
 
-  subroutine z_krm_solver_descr(sv,info,iout,coarse)
+  subroutine z_krm_solver_descr(sv,info,iout,coarse,prefix)
 
     Implicit None
 
@@ -444,12 +444,14 @@ contains
     class(amg_z_krm_solver_type), intent(in) :: sv
     integer(psb_ipk_), intent(out)             :: info
     integer(psb_ipk_), intent(in), optional    :: iout
-    logical, intent(in), optional       :: coarse
+    logical, intent(in), optional              :: coarse
+    character(len=*), intent(in), optional  :: prefix
 
     ! Local variables
     integer(psb_ipk_)      :: err_act
     character(len=20), parameter :: name='amg_z_krm_solver_descr'
     integer(psb_ipk_) :: iout_
+    character(1024)    :: prefix_
 
     call psb_erractionsave(err_act)
     info = psb_success_
@@ -458,23 +460,22 @@ contains
     else
       iout_ = psb_out_unit
     endif
+    if (present(prefix)) then
+      prefix_ = prefix
+    else
+      prefix_ = ""
+    end if
 
     if (sv%global) then
-      write(iout_,*) '  Krylov solver (global)'
+      write(iout_,*) trim(prefix_), '  Krylov solver (global)'
     else
-      write(iout_,*) '  Krylov solver (local) '
+      write(iout_,*) trim(prefix_), '  Krylov solver (local) '
     end if
-    write(iout_,*) '    method: ',sv%method
-    write(iout_,*) '     kprec: ',sv%kprec
-    if (sv%i_sub_solve > 0) then
-      write(iout_,*) ' sub_solve: ',amg_fact_names(sv%i_sub_solve)
-    else
-      write(iout_,*) ' sub_solve: ',sv%sub_solve
-    end if
-    write(iout_,*) '     itmax: ',sv%itmax
-    write(iout_,*) '       eps: ',sv%eps
-    write(iout_,*) '    fillin: ',sv%fillin
-
+    write(iout_,*) trim(prefix_), '    method: ',sv%method
+    write(iout_,*) trim(prefix_), '     kprec: ',sv%kprec
+    call sv%prec%descr(info,iout_,prefix='KRM : '//prefix_)
+    write(iout_,*) trim(prefix_), '     itmax: ',sv%itmax
+    write(iout_,*) trim(prefix_), '       eps: ',sv%eps
 
     call psb_erractionrestore(err_act)
     return

amgprec/amg_z_mumps_solver.F90

@@ -313,22 +313,24 @@ subroutine z_mumps_solver_finalize(sv)
 
 end subroutine z_mumps_solver_finalize
 
-subroutine z_mumps_solver_descr(sv,info,iout,coarse)
+subroutine z_mumps_solver_descr(sv,info,iout,coarse,prefix)
 
   Implicit None
 
   ! Arguments
   class(amg_z_mumps_solver_type), intent(in) :: sv
-  integer(psb_ipk_), intent(out)                     :: info
-  integer(psb_ipk_), intent(in), optional            :: iout
-  logical, intent(in), optional       :: coarse
-
+  integer(psb_ipk_), intent(out)               :: info
+  integer(psb_ipk_), intent(in), optional      :: iout
+  logical, intent(in), optional                :: coarse
+  character(len=*), intent(in), optional       :: prefix
+      
   ! Local variables
   integer(psb_ipk_)   :: err_act
   type(psb_ctxt_type) :: ctxt
   integer(psb_ipk_)   :: me, np
   character(len=20), parameter :: name='amg_z_mumps_solver_descr'
   integer(psb_ipk_) :: iout_
+  character(1024)    :: prefix_
 
   call psb_erractionsave(err_act)
   info = psb_success_
@@ -337,8 +339,13 @@ subroutine z_mumps_solver_descr(sv,info,iout,coarse)
   else
     iout_ = psb_out_unit
   endif
+  if (present(prefix)) then
+    prefix_ = prefix
+  else
+    prefix_ = ""
+  end if
 
-  write(iout_,*) '  MUMPS  Solver. '
+  write(iout_,*) trim(prefix_), '  MUMPS  Solver. '
 
   call psb_erractionrestore(err_act)
   return

amgprec/amg_z_onelev_mod.f90

@@ -257,7 +257,7 @@ module amg_z_onelev_mod
   end interface
 
   interface
-    subroutine amg_z_base_onelev_descr(lv,il,nl,ilmin,info,iout, verbosity)
+    subroutine amg_z_base_onelev_descr(lv,il,nl,ilmin,info,iout, verbosity,prefix)
       import :: psb_zspmat_type, psb_z_vect_type, psb_z_base_vect_type, &
            & psb_zlinmap_type, psb_dpk_, amg_z_onelev_type, &
            & psb_ipk_, psb_epk_, psb_desc_type
@@ -268,6 +268,7 @@ module amg_z_onelev_mod
       integer(psb_ipk_), intent(out)          :: info
       integer(psb_ipk_), intent(in), optional :: iout
       integer(psb_ipk_), intent(in), optional :: verbosity
+      character(len=*), intent(in), optional  :: prefix
     end subroutine amg_z_base_onelev_descr
   end interface
 

amgprec/amg_z_prec_type.f90

@@ -155,15 +155,16 @@ module amg_z_prec_type
 
 
   interface amg_precdescr
-    subroutine amg_zfile_prec_descr(prec,info,iout,root,verbosity)
+    subroutine amg_zfile_prec_descr(prec,info,iout,root,verbosity,prefix)
       import :: amg_zprec_type, psb_ipk_
       implicit none
       ! Arguments
-      class(amg_zprec_type), intent(in)      :: prec
-      integer(psb_ipk_), intent(out) :: info
+      class(amg_zprec_type), intent(in)     :: prec
+      integer(psb_ipk_), intent(out)          :: info
       integer(psb_ipk_), intent(in), optional :: iout
       integer(psb_ipk_), intent(in), optional :: root
       integer(psb_ipk_), intent(in), optional :: verbosity
+      character(len=*), intent(in), optional  :: prefix
     end subroutine amg_zfile_prec_descr
   end interface
 

amgprec/amg_z_slu_solver.F90

@@ -385,20 +385,22 @@ contains
   
   end subroutine z_slu_solver_finalize
 
-  subroutine z_slu_solver_descr(sv,info,iout,coarse)
+  subroutine z_slu_solver_descr(sv,info,iout,coarse,prefix)
 
     Implicit None
 
     ! Arguments
     class(amg_z_slu_solver_type), intent(in) :: sv
-    integer, intent(out)                     :: info
-    integer, intent(in), optional            :: iout
-    logical, intent(in), optional       :: coarse
+    integer, intent(out)                       :: info
+    integer, intent(in), optional              :: iout
+    logical, intent(in), optional              :: coarse
+    character(len=*), intent(in), optional     :: prefix
 
     ! Local variables
     integer      :: err_act
     character(len=20), parameter :: name='amg_z_slu_solver_descr'
     integer :: iout_
+    character(1024)    :: prefix_
 
     call psb_erractionsave(err_act)
     info = psb_success_
@@ -407,8 +409,13 @@ contains
     else
       iout_ = psb_out_unit
     endif
+    if (present(prefix)) then
+      prefix_ = prefix
+    else
+      prefix_ = ""
+    end if
     
-    write(iout_,*) '  SuperLU Sparse Factorization Solver. '
+    write(iout_,*) trim(prefix_), '  SuperLU Sparse Factorization Solver. '
 
     call psb_erractionrestore(err_act)
     return

amgprec/amg_z_sludist_solver.F90

@@ -270,13 +270,13 @@ contains
     ! Local variables
     type(psb_zspmat_type) :: atmp
     type(psb_z_csr_sparse_mat) :: acsr
-    integer(psb_lpk_), allocatable :: gia(:), gja(:)
     type(psb_ctxt_type) :: ctxt
+    integer(psb_lpk_), allocatable :: gia(:), gja(:)
     integer(psb_lpk_) :: lfrst
     integer(psb_ipk_) :: n_row,n_col, nrow_a, nztota, nglob, nzt, npr, npc
     integer(psb_ipk_) :: ifrst, ibcheck
     integer(psb_ipk_) :: np,me,i, err_act, debug_unit, debug_level
-    character(len=20) :: name='z_sludist_solver_bld', ch_err
+    character(len=20)  :: name='z_sludist_solver_bld', ch_err
     
     info=psb_success_
     call psb_erractionsave(err_act)
@@ -324,7 +324,7 @@ contains
     acsr%ja(1:nztota) = gja(1:nztota)
     acsr%ja(:)  = acsr%ja(:)  - 1
     acsr%irp(:) = acsr%irp(:) - 1
-    ifrst = lfrst - 1    
+    ifrst = lfrst - 1
     info = amg_zsludist_fact(nglob,nrow_a,nztota,ifrst,&
          & acsr%val,acsr%irp,acsr%ja,sv%lufactors,&
          & npr,npc)
@@ -421,15 +421,16 @@ contains
   
   end subroutine z_sludist_solver_finalize
 
-  subroutine z_sludist_solver_descr(sv,info,iout,coarse)
+  subroutine z_sludist_solver_descr(sv,info,iout,coarse,prefix)
 
     Implicit None
 
     ! Arguments
     class(amg_z_sludist_solver_type), intent(in) :: sv
-    integer, intent(out)                     :: info
-    integer, intent(in), optional            :: iout
-    logical, intent(in), optional       :: coarse
+    integer, intent(out)                           :: info
+    integer, intent(in), optional                  :: iout
+    logical, intent(in), optional                  :: coarse
+    character(len=*), intent(in), optional         :: prefix
 
     ! Local variables
     integer      :: err_act
@@ -437,6 +438,7 @@ contains
     integer      :: me, np
     character(len=20), parameter :: name='amg_z_sludist_solver_descr'
     integer :: iout_
+    character(1024)    :: prefix_
 
     call psb_erractionsave(err_act)
     info = psb_success_
@@ -445,8 +447,13 @@ contains
     else
       iout_ = psb_out_unit
     endif
+    if (present(prefix)) then
+      prefix_ = prefix
+    else
+      prefix_ = ""
+    end if
     
-    write(iout_,*) '  SuperLU_Dist Sparse Factorization Solver. '
+    write(iout_,*) trim(prefix_), '  SuperLU_Dist Sparse Factorization Solver. '
 
     call psb_erractionrestore(err_act)
     return

amgprec/amg_z_symdec_aggregator_mod.f90

@@ -88,16 +88,25 @@ contains
     val = "Symmetric Decoupled aggregation"
   end function amg_z_symdec_aggregator_fmt
   
-  subroutine  amg_z_symdec_aggregator_descr(ag,parms,iout,info)
+  subroutine  amg_z_symdec_aggregator_descr(ag,parms,iout,info,prefix)
     implicit none 
     class(amg_z_symdec_aggregator_type), intent(in) :: ag
     type(amg_dml_parms), intent(in)   :: parms
     integer(psb_ipk_), intent(in)  :: iout
     integer(psb_ipk_), intent(out) :: info
+    character(len=*), intent(in), optional  :: prefix
 
-    write(iout,*) 'Decoupled Aggregator locally-symmetrized'
-    write(iout,*) 'Aggregator object type: ',ag%fmt()
-    call parms%mldescr(iout,info)
+    character(1024)    :: prefix_
+
+    if (present(prefix)) then
+      prefix_ = prefix
+    else
+      prefix_ = ""
+    end if
+
+    write(iout,*) trim(prefix_),' ','Decoupled Aggregator locally-symmetrized'
+    write(iout,*) trim(prefix_),' ','Aggregator object type: ',ag%fmt()
+    call parms%mldescr(iout,info,prefix=prefix)
     
     return
   end subroutine amg_z_symdec_aggregator_descr

amgprec/amg_z_umf_solver.F90

@@ -390,20 +390,22 @@ contains
   
   end subroutine z_umf_solver_finalize
 
-  subroutine z_umf_solver_descr(sv,info,iout,coarse)
+  subroutine z_umf_solver_descr(sv,info,iout,coarse,prefix)
 
     Implicit None
 
     ! Arguments
     class(amg_z_umf_solver_type), intent(in) :: sv
-    integer, intent(out)                     :: info
-    integer, intent(in), optional            :: iout
-    logical, intent(in), optional       :: coarse
+    integer, intent(out)                       :: info
+    integer, intent(in), optional              :: iout
+    logical, intent(in), optional              :: coarse
+      character(len=*), intent(in), optional   :: prefix
 
     ! Local variables
     integer      :: err_act
     character(len=20), parameter :: name='amg_z_umf_solver_descr'
     integer :: iout_
+    character(1024)    :: prefix_
 
     call psb_erractionsave(err_act)
     info = psb_success_
@@ -412,8 +414,13 @@ contains
     else
       iout_ = psb_out_unit
     endif
+    if (present(prefix)) then
+      prefix_ = prefix
+    else
+      prefix_ = ""
+    end if
     
-    write(iout_,*) '  UMFPACK Sparse Factorization Solver. '
+    write(iout_,*) trim(prefix_), '  UMFPACK Sparse Factorization Solver. '
 
     call psb_erractionrestore(err_act)
     return

amgprec/impl/Makefile

@@ -67,22 +67,25 @@ OBJS=$(F90OBJS) $(COBJS) $(MPCOBJS)
 
 LIBNAME=libamg_prec.a
 
+objs: $(OBJS) aggrd levd smoothd solvd
+
 lib:  $(OBJS) aggrd levd smoothd solvd
+	cd aggregator && $(MAKE) lib
+	cd level && $(MAKE) lib
+	cd smoother && $(MAKE) lib
+	cd solver && $(MAKE) lib
 	$(AR) $(HERE)/$(LIBNAME) $(OBJS)
 	$(RANLIB) $(HERE)/$(LIBNAME)
 
 aggrd:
-	$(MAKE) -C aggregator
+	cd aggregator && $(MAKE) objs
 levd:
-	$(MAKE) -C level
+	cd level && $(MAKE)  objs
 smoothd:
-	$(MAKE) -C smoother
+	cd smoother   && $(MAKE)  objs
 solvd:
-	$(MAKE) -C solver
+	cd solver && $(MAKE)  objs
 
-mpobjs: 
-	(make $(MPFOBJS) FC="$(MPFC)" FCOPT="$(FCOPT)")
-	(make $(MPCOBJS) CC="$(MPCC)" CCOPT="$(CCOPT)")
 
 veryclean: clean
 	/bin/rm -f $(LIBNAME)
@@ -91,10 +94,10 @@ clean: solvclean smoothclean levclean aggrclean
 	/bin/rm -f $(OBJS) $(LOCAL_MODS)
 
 aggrclean:
-		$(MAKE) -C aggregator clean
+	cd aggregator && $(MAKE) clean
 levclean:
-		$(MAKE) -C level clean
+	cd level && $(MAKE) clean
 smoothclean:
-		$(MAKE) -C smoother clean
+	cd smoother && 	$(MAKE)  clean
 solvclean:
-		$(MAKE) -C solver clean
+	cd solver && $(MAKE)  clean

amgprec/impl/aggregator/Makefile

@@ -5,7 +5,7 @@ MODDIR=../../../modules
 HERE=../..
 
 FINCLUDES=$(FMFLAG)$(HERE) $(FMFLAG)$(MODDIR) $(FMFLAG)$(INCDIR) $(PSBLAS_INCLUDES)
-CXXINCLUDES=$(FMFLAG)$(HERE) $(FMFLAG)$(INCDIR) $(FMFLAG)/.
+CXXINCLUDES=$(FIFLAG)$(HERE) $(FIFLAG)$(INCDIR) $(FIFLAG)/. -I../../../../ParallelRomaF-main/include -I$(PSBLAS_INCDIR)
 
 #CINCLUDES= -I${SUPERLU_INCDIR} -I${HSL_INCDIR}  -I${SPRAL_INCDIR} -I/home/users/pasqua/Ambra/BootCMatch/include -lBCM -L/home/users/pasqua/Ambra/BootCMatch/lib  -lm
 
@@ -59,23 +59,45 @@ amg_s_parmatch_spmm_bld.o \
 amg_s_parmatch_spmm_bld_ov.o \
 amg_s_parmatch_unsmth_bld.o \
 amg_s_parmatch_smth_bld.o \
-amg_s_parmatch_spmm_bld_inner.o
+amg_s_parmatch_spmm_bld_inner.o \
+amg_d_newmatch_aggregator_inner_mat_asb.o\
+amg_d_newmatch_aggregator_mat_asb.o \
+amg_d_newmatch_aggregator_mat_bld.o \
+amg_d_newmatch_aggregator_tprol.o  \
+amg_d_newmatch_map_to_tprol.o \
+amg_d_newmatch_spmm_bld_inner.o \
+amg_d_newmatch_spmm_bld_ov.o
 
-MPCOBJS=MatchBoxPC.o \
-algoDistEdgeApproxDomEdgesLinearSearchMesgBndlSmallMateC.o
+MPCXXOBJS=MatchBoxPC.o \
+algoDistEdgeApproxDomEdgesLinearSearchMesgBndlSmallMateC.o \
+newmatch_interface.o \
+sendBundledMessages.o \
+initialize.o \
+extractUChunk.o \
+isAlreadyMatched.o \
+findOwnerOfGhost.o \
+clean.o \
+computeCandidateMate.o \
+parallelComputeCandidateMateB.o \
+processMatchedVertices.o \
+processMatchedVerticesAndSendMessages.o \
+processCrossEdge.o \
+queueTransfer.o \
+processMessages.o \
+processExposedVertex.o \
+algoDistEdgeApproxDomEdgesLinearSearchMesgBndlSmallMateCMP.o \
+MatchingAlgorithms.o
 
-OBJS = $(FOBJS) $(MPCOBJS)
+OBJS = $(FOBJS) $(MPCOBJS) $(MPCXXOBJS) 
 
 LIBNAME=libamg_prec.a
 
-lib:  $(OBJS)
+objs: $(OBJS)
+
+lib:  objs
 	$(AR) $(HERE)/$(LIBNAME) $(OBJS)
 	$(RANLIB) $(HERE)/$(LIBNAME)
 
-mpobjs:
-	(make $(MPFOBJS) F90="$(MPF90)" F90COPT="$(F90COPT)")
-	(make $(MPCOBJS) CC="$(MPCC)" CCOPT="$(CCOPT)")
-
 veryclean: clean
 	/bin/rm -f $(LIBNAME)
 

amgprec/impl/aggregator/MatchBoxPC.cpp

@@ -60,17 +60,43 @@ void dMatchBoxPC(MilanLongInt NLVer, MilanLongInt NLEdge,
 		MilanLongInt* ph1_card, MilanLongInt* ph2_card ) {
 #if !defined(SERIAL_MPI)
   MPI_Comm C_comm=MPI_Comm_f2c(icomm);
+
 #ifdef DEBUG
   fprintf(stderr,"MatchBoxPC: rank %d nlver %ld nledge %ld [ %ld %ld ]\n",
 	  myRank,NLVer, NLEdge,verDistance[0],verDistance[1]);
 #endif
-  dalgoDistEdgeApproxDomEdgesLinearSearchMesgBndlSmallMateC(NLVer, NLEdge,
+
+
+#define TIME_TRACKER
+    #ifdef TIME_TRACKER
+        double tmr = MPI_Wtime();
+    #endif
+
+#define OMP
+#ifdef OMP
+        dalgoDistEdgeApproxDomEdgesLinearSearchMesgBndlSmallMateCMP(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
+        dalgoDistEdgeApproxDomEdgesLinearSearchMesgBndlSmallMateC(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
+
+
+  #ifdef TIME_TRACKER
+    tmr = MPI_Wtime() - tmr;
+    fprintf(stderr, "Elaboration time: %f for %ld nodes\n", tmr, NLVer);
+  #endif
+
 #endif
 }
 

amgprec/impl/aggregator/MatchBoxPC.h

@@ -52,145 +52,412 @@
 
 #ifndef _matchboxpC_H_
 #define _matchboxpC_H_
-//Turn on a lot of debugging information with this switch:
+// Turn on a lot of debugging information with this switch:
 //#define PRINT_DEBUG_INFO_
 #include <stdio.h>
 #include <iostream>
 #include <assert.h>
 #include <map>
 #include <vector>
-// #include "matchboxp.h"
+#include "omp.h"
 #include "primitiveDataTypeDefinitions.h"
 #include "dataStrStaticQueue.h"
 
 using namespace std;
 
+const int NUM_THREAD = 4;
+const int UCHUNK = 10;
+
+const MilanLongInt REQUEST = 1;
+const MilanLongInt SUCCESS = 2;
+const MilanLongInt FAILURE = 3;
+const MilanLongInt SIZEINFO = 4;
+
+const int ComputeTag = 7; // Predefined tag
+const int BundleTag = 9;  // Predefined tag
+
+static vector<MilanLongInt> DEFAULT_VECTOR;
+
+// MPI type map
+template <typename T>
+MPI_Datatype TypeMap();
+template <>
+inline MPI_Datatype TypeMap<int64_t>() { return MPI_LONG_LONG; }
+template <>
+inline MPI_Datatype TypeMap<int>() { return MPI_INT; }
+template <>
+inline MPI_Datatype TypeMap<double>() { return MPI_DOUBLE; }
+template <>
+inline MPI_Datatype TypeMap<float>() { return MPI_FLOAT; }
+
 #ifdef __cplusplus
-extern "C" {
+extern "C"
+{
 #endif
 
 #if !defined(SERIAL_MPI)
-  
-#define MilanMpiLongInt  MPI_LONG_LONG
+
+#define MilanMpiLongInt MPI_LONG_LONG
 
 #ifndef _primitiveDataType_Definition_
 #define _primitiveDataType_Definition_
-    //Regular integer:
-    #ifndef INTEGER_H
-    #define INTEGER_H
-        typedef int32_t MilanInt;
-    #endif
+// Regular integer:
+#ifndef INTEGER_H
+#define INTEGER_H
+    typedef int32_t MilanInt;
+#endif
 
-    //Regular long integer:
-    #ifndef LONG_INT_H
-    #define LONG_INT_H
-        #ifdef BIT64
-            typedef int64_t MilanLongInt;
-            typedef MPI_LONG MilanMpiLongInt;
-        #else
-            typedef int32_t MilanLongInt;
-            typedef MPI_INT MilanMpiLongInt;
-        #endif
-    #endif
+// Regular long integer:
+#ifndef LONG_INT_H
+#define LONG_INT_H
+#ifdef BIT64
+    typedef int64_t MilanLongInt;
+    typedef MPI_LONG MilanMpiLongInt;
+#else
+    typedef int32_t MilanLongInt;
+    typedef MPI_INT MilanMpiLongInt;
+#endif
+#endif
 
-    //Regular boolean
-    #ifndef BOOL_H
-    #define BOOL_H
-        typedef bool MilanBool;
-    #endif
+// Regular boolean
+#ifndef BOOL_H
+#define BOOL_H
+    typedef bool MilanBool;
+#endif
 
-    //Regular double and absolute value computation:
-    #ifndef REAL_H
-    #define REAL_H
-        typedef double MilanReal;
-        typedef MPI_DOUBLE MilanMpiReal;
-        inline MilanReal MilanAbs(MilanReal value)
-        {
-            return fabs(value);
-        }
-    #endif
+// Regular double and absolute value computation:
+#ifndef REAL_H
+#define REAL_H
+    typedef double MilanReal;
+    typedef MPI_DOUBLE MilanMpiReal;
+    inline MilanReal MilanAbs(MilanReal value)
+    {
+        return fabs(value);
+    }
+#endif
 
-    //Regular float and absolute value computation:
-    #ifndef FLOAT_H
-    #define FLOAT_H
-        typedef float MilanFloat;
-        typedef MPI_FLOAT MilanMpiFloat;
-        inline MilanFloat MilanAbsFloat(MilanFloat value)
-        {
-            return fabs(value);
-        }
-    #endif
+// Regular float and absolute value computation:
+#ifndef FLOAT_H
+#define FLOAT_H
+    typedef float MilanFloat;
+    typedef MPI_FLOAT MilanMpiFloat;
+    inline MilanFloat MilanAbsFloat(MilanFloat value)
+    {
+        return fabs(value);
+    }
+#endif
 
-    //// Define the limits:
-    #ifndef LIMITS_H
-    #define LIMITS_H
-    //Integer Maximum and Minimum:
-  //      #define MilanIntMax INT_MAX
-  //    #define MilanIntMin INT_MIN
-        #define MilanIntMax INT32_MAX
-        #define MilanIntMin INT32_MIN
+//// Define the limits:
+#ifndef LIMITS_H
+#define LIMITS_H
+    // Integer Maximum and Minimum:
+    //      #define MilanIntMax INT_MAX
+    //    #define MilanIntMin INT_MIN
+#define MilanIntMax INT32_MAX
+#define MilanIntMin INT32_MIN
 
-        #ifdef BIT64
-            #define MilanLongIntMax INT64_MAX
-            #define MilanLongIntMin -INT64_MAX
-        #else
-            #define MilanLongIntMax INT32_MAX
-            #define MilanLongIntMin -INT32_MAX
-        #endif
+#ifdef BIT64
+#define MilanLongIntMax INT64_MAX
+#define MilanLongIntMin -INT64_MAX
+#else
+#define MilanLongIntMax INT32_MAX
+#define MilanLongIntMin -INT32_MAX
+#endif
 
-    #endif
+#endif
 
     // +INFINITY
     const double PLUS_INFINITY = numeric_limits<int>::infinity();
     const double MINUS_INFINITY = -PLUS_INFINITY;
-    //#define MilanRealMax LDBL_MAX
-    #define MilanRealMax PLUS_INFINITY
-    #define MilanRealMin MINUS_INFINITY
+//#define MilanRealMax LDBL_MAX
+#define MilanRealMax PLUS_INFINITY
+#define MilanRealMin MINUS_INFINITY
 #endif
 
-//Function of find the owner of a ghost vertex using binary search:
-inline MilanInt findOwnerOfGhost(MilanLongInt vtxIndex, MilanLongInt *mVerDistance,
-                                     MilanInt myRank, MilanInt numProcs);
+    // Function of find the owner of a ghost vertex using binary search:
+    MilanInt findOwnerOfGhost(MilanLongInt vtxIndex, MilanLongInt *mVerDistance,
+                              MilanInt myRank, MilanInt numProcs);
 
-  void dalgoDistEdgeApproxDomEdgesLinearSearchMesgBndlSmallMateC
-(
- MilanLongInt NLVer, MilanLongInt NLEdge,
- MilanLongInt* verLocPtr, MilanLongInt* verLocInd, MilanReal* edgeLocWeight,
- MilanLongInt* verDistance,
- MilanLongInt* Mate,
- MilanInt myRank, MilanInt numProcs, MPI_Comm comm,
- MilanLongInt* msgIndSent, MilanLongInt* msgActualSent, MilanReal* msgPercent,
- MilanReal* ph0_time, MilanReal* ph1_time, MilanReal* ph2_time,
- MilanLongInt* ph1_card, MilanLongInt* ph2_card );
+    MilanLongInt firstComputeCandidateMate(MilanLongInt adj1,
+                                           MilanLongInt adj2,
+                                           MilanLongInt *verLocInd,
+                                           MilanReal *edgeLocWeight);
 
- 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 );
+    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);
 
-void dMatchBoxPC(MilanLongInt NLVer, MilanLongInt NLEdge,
-		MilanLongInt* verLocPtr, MilanLongInt* verLocInd, MilanReal* edgeLocWeight,
-		MilanLongInt* verDistance,
-		MilanLongInt* Mate,
-		MilanInt myRank, MilanInt numProcs, MilanInt icomm,
-		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 sMatchBoxPC(MilanLongInt NLVer, MilanLongInt NLEdge,
-		MilanLongInt* verLocPtr, MilanLongInt* verLocInd, MilanFloat* edgeLocWeight,
-		MilanLongInt* verDistance,
-		MilanLongInt* Mate,
-		MilanInt myRank, MilanInt numProcs, MilanInt icomm,
-		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);
+
+    void dalgoDistEdgeApproxDomEdgesLinearSearchMesgBndlSmallMateC(
+        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);
+
+    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);
+
+    void dMatchBoxPC(MilanLongInt NLVer, MilanLongInt NLEdge,
+                     MilanLongInt *verLocPtr, MilanLongInt *verLocInd, MilanReal *edgeLocWeight,
+                     MilanLongInt *verDistance,
+                     MilanLongInt *Mate,
+                     MilanInt myRank, MilanInt numProcs, MilanInt icomm,
+                     MilanLongInt *msgIndSent, MilanLongInt *msgActualSent, MilanReal *msgPercent,
+                     MilanReal *ph0_time, MilanReal *ph1_time, MilanReal *ph2_time,
+                     MilanLongInt *ph1_card, MilanLongInt *ph2_card);
+
+    void sMatchBoxPC(MilanLongInt NLVer, MilanLongInt NLEdge,
+                     MilanLongInt *verLocPtr, MilanLongInt *verLocInd, MilanFloat *edgeLocWeight,
+                     MilanLongInt *verDistance,
+                     MilanLongInt *Mate,
+                     MilanInt myRank, MilanInt numProcs, MilanInt icomm,
+                     MilanLongInt *msgIndSent, MilanLongInt *msgActualSent, MilanReal *msgPercent,
+                     MilanReal *ph0_time, MilanReal *ph1_time, MilanReal *ph2_time,
+                     MilanLongInt *ph1_card, MilanLongInt *ph2_card);
 
 #endif
 #ifdef __cplusplus

amgprec/impl/aggregator/algoDistEdgeApproxDomEdgesLinearSearchMesgBndlSmallMateC.cpp

@@ -72,12 +72,6 @@
 
 #ifdef SERIAL_MPI
 #else
-//MPI type map
-template<typename T> MPI_Datatype TypeMap();
-template<> inline MPI_Datatype TypeMap<int64_t>() { return MPI_LONG_LONG; }
-template<> inline MPI_Datatype TypeMap<int>() { return MPI_INT; }
-template<> inline MPI_Datatype TypeMap<double>() { return MPI_DOUBLE; }
-template<> inline MPI_Datatype TypeMap<float>() { return MPI_FLOAT; }
 
 // DOUBLE PRECISION VERSION
 //WARNING: The vertex block on a given rank is contiguous

amgprec/impl/aggregator/algoDistEdgeApproxDomEdgesLinearSearchMesgBndlSmallMateCMP.cpp

@@ -0,0 +1,554 @@
+#include "MatchBoxPC.h"
+
+// ***********************************************************************
+//
+//        MatchboxP: A C++ library for approximate weighted matching
+//               Mahantesh Halappanavar (hala@pnnl.gov)
+//               Pacific Northwest National Laboratory
+//
+// ***********************************************************************
+//
+//       Copyright (2021) Battelle Memorial Institute
+//                      All rights reserved.
+//
+// Redistribution and use in source and binary forms, with or without
+// modification, are permitted provided that the following conditions
+// are met:
+//
+// 1. Redistributions of source code must retain the above copyright
+// notice, this list of conditions and the following disclaimer.
+//
+// 2. Redistributions in binary form must reproduce the above copyright
+// notice, this list of conditions and the following disclaimer in the
+// documentation and/or other materials provided with the distribution.
+//
+// 3. Neither the name of the copyright holder nor the names of its
+// contributors may be used to endorse or promote products derived from
+// this software without specific prior written permission.
+//
+// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+// "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+// LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
+// FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
+// COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
+// INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
+// BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
+// LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
+// CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
+// LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
+// ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+// POSSIBILITY OF SUCH DAMAGE.
+//
+// ************************************************************************
+//////////////////////////////////////////////////////////////////////////////////////
+/////////////////////////// DOMINATING EDGES MODEL ///////////////////////////////////
+//////////////////////////////////////////////////////////////////////////////////////
+/* Function	: algoDistEdgeApproxDomEdgesLinearSearchMesgBndlSmallMate()
+ *
+ * Date     : New update: Feb 17, 2019, Richland, Washington.
+ * Date		: Original development: May 17, 2009, E&CS Bldg.
+ *
+ * Purpose	: Compute Approximate Maximum Weight Matching in Linear Time
+ *
+ * Args		: inputMatrix - instance of Compressed-Col format of Matrix
+ *                Mate - The Mate array
+ *
+ * Returns	: By Value: (void)
+ *            By Reference: Mate
+ *
+ * Comments	: 1/2 Approx Algorithm. Picks the locally available heaviest edge.
+ *                Assumption: The Mate Array is empty.
+ */
+
+/*
+ NLVer = #of vertices, NLEdge = #of edges
+ CSR/CSC/Compressed format: verLocPtr = Pointer, verLocInd = Index, edgeLocWeight = edge weights (positive real numbers)
+ verDistance = A vector of size |P|+1 containing the cumulative number of vertices per process
+ Mate = A vector of size |V_p| (local subgraph) to store the output (matching)
+ MPI: myRank, numProcs, comm,
+ Statistics: msgIndSent, msgActualSent, msgPercent : Size: |P| number of processes in the comm-world
+ Statistics: ph0_time, ph1_time, ph2_time: Runtimes
+ Statistics: ph1_card, ph2_card : Size: |P| number of processes in the comm-world (number of matched edges in Phase 1 and Phase 2)
+ */
+//#define DEBUG_HANG_
+#ifdef SERIAL_MPI
+#else
+
+// DOUBLE PRECISION VERSION
+// WARNING: The vertex block on a given rank is contiguous
+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)
+{
+
+    /*
+     * 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
+
+    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;
+
+    MilanLongInt msgActual = 0, msgInd = 0;
+    MilanReal 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;
+    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++)
+        PCounter[i] = 0;
+
+    MilanLongInt NumMessagesBundled = 0;
+    // TODO when the last computational section will be refactored this could be eliminated
+    MilanInt ghostOwner = 0; // Changed by Fabio to be an integer, addresses needs to be integers!
+    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<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...";
+    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_B(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_B(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);
+
+    processMatchedVertices(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 //////////////////////////////
+        ///////////////////////////////////////////////////////////////////////////////////
+      
+      processMatchedVerticesAndSendMessages(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 //////////////////////////////////////
+        ///////////////////////////////////////////////////////////////////////////////////
+
+        processMessages(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
+}
+// End of algoDistEdgeApproxDomEdgesLinearSearchMesgBndlSmallMate
+#endif
+
+#endif

amgprec/impl/aggregator/amg_c_dec_aggregator_tprol.f90

@@ -97,6 +97,8 @@ subroutine  amg_c_dec_aggregator_build_tprol(ag,parms,ag_data,&
   integer(psb_lpk_)   :: ntaggr
   integer(psb_ipk_)   :: debug_level, debug_unit
   logical             :: clean_zeros
+  integer(psb_ipk_), save :: idx_map_bld=-1, idx_map_tprol=-1
+  logical, parameter      :: do_timings=.false.
 
   name='amg_c_dec_aggregator_tprol'
   call psb_erractionsave(err_act)
@@ -108,6 +110,10 @@ subroutine  amg_c_dec_aggregator_build_tprol(ag,parms,ag_data,&
   info  = psb_success_
   ctxt = desc_a%get_context()
   call psb_info(ctxt,me,np)
+  if ((do_timings).and.(idx_map_bld==-1))       &
+       & idx_map_bld = psb_get_timer_idx("DEC_TPROL: map_bld")
+  if ((do_timings).and.(idx_map_tprol==-1))     &
+       & idx_map_tprol = psb_get_timer_idx("DEC_TPROL: map_tprol")
 
   call amg_check_def(parms%ml_cycle,'Multilevel cycle',&
        &   amg_mult_ml_,is_legal_ml_cycle)
@@ -121,10 +127,14 @@ subroutine  amg_c_dec_aggregator_build_tprol(ag,parms,ag_data,&
   ! The decoupled aggregator based on SOC measures ignores
   ! ag_data except for clean_zeros; soc_map_bld is a procedure pointer.
   !
+  if (do_timings) call psb_tic(idx_map_bld)
   clean_zeros = ag%do_clean_zeros
   call ag%soc_map_bld(parms%aggr_ord,parms%aggr_thresh,clean_zeros,a,desc_a,nlaggr,ilaggr,info)
+  if (do_timings) call psb_toc(idx_map_bld)
+  if (do_timings) call psb_tic(idx_map_tprol)
 
   if (info==psb_success_) call amg_map_to_tprol(desc_a,ilaggr,nlaggr,t_prol,info)
+  if (do_timings) call psb_toc(idx_map_tprol)
   if (info /= psb_success_) then
     info=psb_err_from_subroutine_
     call psb_errpush(info,name,a_err='soc_map_bld/map_to_tprol')