Adjusted CBIND test program.

This reverts commit a17f503486.

Make.inc.in

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

Make_n.inc.in

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

Makefile

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

amgprec/amg_base_prec_type.F90

@@ -323,9 +323,10 @@ module amg_base_prec_type
   !
   ! Legal values for entry: amg_poly_variant_
   !
-  integer(psb_ipk_), parameter :: amg_poly_lottes_      = 0
-  integer(psb_ipk_), parameter :: amg_poly_lottes_beta_ = 1
-  integer(psb_ipk_), parameter :: amg_poly_new_         = 2
+  integer(psb_ipk_), parameter :: amg_cheb_4_     = 0
+  integer(psb_ipk_), parameter :: amg_cheb_4_opt_ = 1
+  integer(psb_ipk_), parameter :: amg_cheb_1_opt_ = 2
+  integer(psb_ipk_), parameter :: amg_poly_dbg_   = 8
   
   integer(psb_ipk_), parameter :: amg_poly_rho_est_power_ = 0
 
@@ -569,12 +570,14 @@ contains
       val = amg_as_
     case('POLY')
       val = amg_poly_
-    case('POLY_LOTTES')
-      val = amg_poly_lottes_
-    case('POLY_LOTTES_BETA')
-      val = amg_poly_lottes_beta_
-    case('POLY_NEW')
-      val = amg_poly_new_
+    case('CHEB_4')
+      val = amg_cheb_4_
+    case('CHEB_4_OPT')
+      val = amg_cheb_4_opt_
+    case('CHEB_1_OPT')
+      val = amg_cheb_1_opt_
+    case('POLY_DBG')
+      val = amg_poly_dbg_
     case('POLY_RHO_EST_POWER')
       val =  amg_poly_rho_est_power_
     case('A_NORMI')

amgprec/amg_d_poly_smoother.f90

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

amgprec/amg_s_poly_smoother.f90

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

amgprec/impl/aggregator/MatchBoxPC.cpp

@@ -67,14 +67,13 @@ void dMatchBoxPC(MilanLongInt NLVer, MilanLongInt NLEdge,
 #endif
 
 
-#define TIME_TRACKER
-    #ifdef TIME_TRACKER
-        double tmr = MPI_Wtime();
-    #endif
+#undef TIME_TRACKER
+#ifdef TIME_TRACKER
+  double tmr = MPI_Wtime();
+#endif
 
-// Rimosso per tornare al vecchio matching #define OMP
-#ifdef OMP
-				fprintf(stderr,"Warning: using buggy OpenMP matching!\n");
+#if defined(OPENMP)
+	//fprintf(stderr,"Warning: using buggy OpenMP matching!\n");
         dalgoDistEdgeApproxDomEdgesLinearSearchMesgBndlSmallMateCMP(NLVer, NLEdge,
 							   verLocPtr, verLocInd, edgeLocWeight,
 							   verDistance,  Mate,
@@ -93,11 +92,11 @@ void dMatchBoxPC(MilanLongInt NLVer, MilanLongInt NLEdge,
 #endif
 
 
-  #ifdef TIME_TRACKER
-    tmr = MPI_Wtime() - tmr;
-    fprintf(stderr, "Elaboration time: %f for %ld nodes\n", tmr, NLVer);
-  #endif
-
+#ifdef TIME_TRACKER
+	tmr = MPI_Wtime() - tmr;
+	fprintf(stderr, "Elaboration time: %f for %ld nodes\n", tmr, NLVer);
+#endif
+	
 #endif
 }
 
@@ -115,13 +114,24 @@ void sMatchBoxPC(MilanLongInt NLVer, MilanLongInt NLEdge,
   fprintf(stderr,"MatchBoxPC: rank %d nlver %ld nledge %ld [ %ld %ld ]\n",
 	  myRank,NLVer, NLEdge,verDistance[0],verDistance[1]);
 #endif
+#if defined(OPENMP)
+	//fprintf(stderr,"Warning: using buggy OpenMP matching!\n");
+  salgoDistEdgeApproxDomEdgesLinearSearchMesgBndlSmallMateCMP(NLVer, NLEdge,
+							      verLocPtr, verLocInd, edgeLocWeight,
+							      verDistance,  Mate,
+							      myRank, numProcs, C_comm,
+							      msgIndSent, msgActualSent, msgPercent,
+							      ph0_time, ph1_time, ph2_time,
+							      ph1_card, ph2_card );
+#else
   salgoDistEdgeApproxDomEdgesLinearSearchMesgBndlSmallMateC(NLVer, NLEdge,
-							   verLocPtr, verLocInd, edgeLocWeight,
-							   verDistance,  Mate,
-							   myRank, numProcs, C_comm,
-							   msgIndSent, msgActualSent, msgPercent,
-							   ph0_time, ph1_time, ph2_time,
-							   ph1_card, ph2_card );
+							    verLocPtr, verLocInd, edgeLocWeight,
+							    verDistance,  Mate,
+							    myRank, numProcs, C_comm,
+							    msgIndSent, msgActualSent, msgPercent,
+							    ph0_time, ph1_time, ph2_time,
+							    ph1_card, ph2_card );
+#endif
 #endif
 }
 

amgprec/impl/aggregator/MatchBoxPC.h

@@ -59,7 +59,7 @@
 #include <assert.h>
 #include <map>
 #include <vector>
-#ifdef OMP
+#ifdef OPENMP
 // OpenMP is included and used if and only if the OpenMP version of the matching
 // is required
 #include "omp.h"
@@ -82,6 +82,8 @@ const int BundleTag = 9;  // Predefined tag
 
 static vector<MilanLongInt> DEFAULT_VECTOR;
 
+#if !defined(SERIAL_MPI)
+
 // MPI type map
 template <typename T>
 MPI_Datatype TypeMap();
@@ -93,6 +95,7 @@ template <>
 inline MPI_Datatype TypeMap<double>() { return MPI_DOUBLE; }
 template <>
 inline MPI_Datatype TypeMap<float>() { return MPI_FLOAT; }
+#endif
 
 #ifdef __cplusplus
 extern "C"
@@ -178,264 +181,416 @@ extern "C"
 #define MilanRealMin MINUS_INFINITY
 #endif
 
-#ifdef OMP
 /* These functions are only used in the experimental OMP implementation, if that
 is disabled there is no reason to actually compile or reference them. */
 
-    // Function of find the owner of a ghost vertex using binary search:
-    MilanInt findOwnerOfGhost(MilanLongInt vtxIndex, MilanLongInt *mVerDistance,
-                              MilanInt myRank, MilanInt numProcs);
+  // Function of find the owner of a ghost vertex using binary search:
+  MilanInt findOwnerOfGhost(MilanLongInt vtxIndex, MilanLongInt *mVerDistance,
+			    MilanInt myRank, MilanInt numProcs);
+  
+  MilanLongInt firstComputeCandidateMateD(MilanLongInt adj1,
+					 MilanLongInt adj2,
+					 MilanLongInt *verLocInd,
+					 MilanReal *edgeLocWeight);
 
-    MilanLongInt firstComputeCandidateMate(MilanLongInt adj1,
-                                           MilanLongInt adj2,
-                                           MilanLongInt *verLocInd,
-                                           MilanReal *edgeLocWeight);
+  
+  void queuesTransfer(vector<MilanLongInt> &U,
+		      vector<MilanLongInt> &privateU,
+		      vector<MilanLongInt> &QLocalVtx,
+		      vector<MilanLongInt> &QGhostVtx,
+		      vector<MilanLongInt> &QMsgType,
+		      vector<MilanInt> &QOwner,
+		      vector<MilanLongInt> &privateQLocalVtx,
+		      vector<MilanLongInt> &privateQGhostVtx,
+		      vector<MilanLongInt> &privateQMsgType,
+		      vector<MilanInt> &privateQOwner);
+  
+  bool isAlreadyMatched(MilanLongInt node,
+			MilanLongInt StartIndex,
+			MilanLongInt EndIndex,
+			vector<MilanLongInt> &GMate,
+			MilanLongInt *Mate,
+			map<MilanLongInt, MilanLongInt> &Ghost2LocalMap);
+  
+  MilanLongInt computeCandidateMateD(MilanLongInt adj1,
+				    MilanLongInt adj2,
+				    MilanReal *edgeLocWeight,
+				    MilanLongInt k,
+				    MilanLongInt *verLocInd,
+				    MilanLongInt StartIndex,
+				    MilanLongInt EndIndex,
+				    vector<MilanLongInt> &GMate,
+				    MilanLongInt *Mate,
+				    map<MilanLongInt, MilanLongInt> &Ghost2LocalMap);
+  
+  void initialize(MilanLongInt NLVer, MilanLongInt NLEdge,
+		  MilanLongInt StartIndex, MilanLongInt EndIndex,
+		  MilanLongInt *numGhostEdgesPtr,
+		  MilanLongInt *numGhostVerticesPtr,
+		  MilanLongInt *S,
+		  MilanLongInt *verLocInd,
+		  MilanLongInt *verLocPtr,
+		  map<MilanLongInt, MilanLongInt> &Ghost2LocalMap,
+		  vector<MilanLongInt> &Counter,
+		  vector<MilanLongInt> &verGhostPtr,
+		  vector<MilanLongInt> &verGhostInd,
+		  vector<MilanLongInt> &tempCounter,
+		  vector<MilanLongInt> &GMate,
+		  vector<MilanLongInt> &Message,
+		  vector<MilanLongInt> &QLocalVtx,
+		  vector<MilanLongInt> &QGhostVtx,
+		  vector<MilanLongInt> &QMsgType,
+		  vector<MilanInt> &QOwner,
+		  MilanLongInt *&candidateMate,
+		  vector<MilanLongInt> &U,
+		  vector<MilanLongInt> &privateU,
+		  vector<MilanLongInt> &privateQLocalVtx,
+		  vector<MilanLongInt> &privateQGhostVtx,
+		  vector<MilanLongInt> &privateQMsgType,
+		  vector<MilanInt> &privateQOwner);
+  
+  void clean(MilanLongInt NLVer,
+	     MilanInt myRank,
+	     MilanLongInt MessageIndex,
+	     vector<MPI_Request> &SRequest,
+	     vector<MPI_Status> &SStatus,
+	     MilanInt BufferSize,
+	     MilanLongInt *Buffer,
+	     MilanLongInt msgActual,
+	     MilanLongInt *msgActualSent,
+	     MilanLongInt msgInd,
+	     MilanLongInt *msgIndSent,
+	     MilanLongInt NumMessagesBundled,
+	     MilanReal *msgPercent);
+  
+  void PARALLEL_COMPUTE_CANDIDATE_MATE_BD(MilanLongInt NLVer,
+					 MilanLongInt *verLocPtr,
+					 MilanLongInt *verLocInd,
+					 MilanInt myRank,
+					 MilanReal *edgeLocWeight,
+					 MilanLongInt *candidateMate);
+  
+  void PARALLEL_PROCESS_EXPOSED_VERTEX_BD(MilanLongInt NLVer,
+					 MilanLongInt *candidateMate,
+					 MilanLongInt *verLocInd,
+					 MilanLongInt *verLocPtr,
+					 MilanLongInt StartIndex,
+					 MilanLongInt EndIndex,
+					 MilanLongInt *Mate,
+					 vector<MilanLongInt> &GMate,
+					 map<MilanLongInt, MilanLongInt> &Ghost2LocalMap,
+					 MilanReal *edgeLocWeight,
+					 MilanLongInt *myCardPtr,
+					 MilanLongInt *msgIndPtr,
+					 MilanLongInt *NumMessagesBundledPtr,
+					 MilanLongInt *SPtr,
+					 MilanLongInt *verDistance,
+					 MilanLongInt *PCounter,
+					 vector<MilanLongInt> &Counter,
+					 MilanInt myRank,
+					 MilanInt numProcs,
+					 vector<MilanLongInt> &U,
+					 vector<MilanLongInt> &privateU,
+					 vector<MilanLongInt> &QLocalVtx,
+					 vector<MilanLongInt> &QGhostVtx,
+					 vector<MilanLongInt> &QMsgType,
+					 vector<MilanInt> &QOwner,
+					 vector<MilanLongInt> &privateQLocalVtx,
+					 vector<MilanLongInt> &privateQGhostVtx,
+					 vector<MilanLongInt> &privateQMsgType,
+					 vector<MilanInt> &privateQOwner);
+  
+  void PROCESS_CROSS_EDGE(MilanLongInt *edge,
+			  MilanLongInt *SPtr);
+  
+  void processMatchedVerticesD(
+			      MilanLongInt NLVer,
+			      vector<MilanLongInt> &UChunkBeingProcessed,
+			      vector<MilanLongInt> &U,
+			      vector<MilanLongInt> &privateU,
+			      MilanLongInt StartIndex,
+			      MilanLongInt EndIndex,
+			      MilanLongInt *myCardPtr,
+			      MilanLongInt *msgIndPtr,
+			      MilanLongInt *NumMessagesBundledPtr,
+			      MilanLongInt *SPtr,
+			      MilanLongInt *verLocPtr,
+			      MilanLongInt *verLocInd,
+			      MilanLongInt *verDistance,
+			      MilanLongInt *PCounter,
+			      vector<MilanLongInt> &Counter,
+			      MilanInt myRank,
+			      MilanInt numProcs,
+			      MilanLongInt *candidateMate,
+			      vector<MilanLongInt> &GMate,
+			      MilanLongInt *Mate,
+			      map<MilanLongInt, MilanLongInt> &Ghost2LocalMap,
+			      MilanReal *edgeLocWeight,
+			      vector<MilanLongInt> &QLocalVtx,
+			      vector<MilanLongInt> &QGhostVtx,
+			      vector<MilanLongInt> &QMsgType,
+			      vector<MilanInt> &QOwner,
+			      vector<MilanLongInt> &privateQLocalVtx,
+			      vector<MilanLongInt> &privateQGhostVtx,
+			      vector<MilanLongInt> &privateQMsgType,
+			      vector<MilanInt> &privateQOwner);
+  
+  void processMatchedVerticesAndSendMessagesD(
+					     MilanLongInt NLVer,
+					     vector<MilanLongInt> &UChunkBeingProcessed,
+					     vector<MilanLongInt> &U,
+					     vector<MilanLongInt> &privateU,
+					     MilanLongInt StartIndex,
+					     MilanLongInt EndIndex,
+					     MilanLongInt *myCardPtr,
+					     MilanLongInt *msgIndPtr,
+					     MilanLongInt *NumMessagesBundledPtr,
+					     MilanLongInt *SPtr,
+					     MilanLongInt *verLocPtr,
+					     MilanLongInt *verLocInd,
+					     MilanLongInt *verDistance,
+					     MilanLongInt *PCounter,
+					     vector<MilanLongInt> &Counter,
+					     MilanInt myRank,
+					     MilanInt numProcs,
+					     MilanLongInt *candidateMate,
+					     vector<MilanLongInt> &GMate,
+					     MilanLongInt *Mate,
+					     map<MilanLongInt, MilanLongInt> &Ghost2LocalMap,
+					     MilanReal *edgeLocWeight,
+					     vector<MilanLongInt> &QLocalVtx,
+					     vector<MilanLongInt> &QGhostVtx,
+					     vector<MilanLongInt> &QMsgType,
+					     vector<MilanInt> &QOwner,
+					     vector<MilanLongInt> &privateQLocalVtx,
+					     vector<MilanLongInt> &privateQGhostVtx,
+					     vector<MilanLongInt> &privateQMsgType,
+					     vector<MilanInt> &privateQOwner,
+					     MPI_Comm comm,
+					     MilanLongInt *msgActual,
+					     vector<MilanLongInt> &Message);
+  
+  void sendBundledMessages(MilanLongInt *numGhostEdgesPtr,
+			   MilanInt *BufferSizePtr,
+			   MilanLongInt *Buffer,
+			   vector<MilanLongInt> &PCumulative,
+			   vector<MilanLongInt> &PMessageBundle,
+			   vector<MilanLongInt> &PSizeInfoMessages,
+			   MilanLongInt *PCounter,
+			   MilanLongInt NumMessagesBundled,
+			   MilanLongInt *msgActualPtr,
+			   MilanLongInt *MessageIndexPtr,
+			   MilanInt numProcs,
+			   MilanInt myRank,
+			   MPI_Comm comm,
+			   vector<MilanLongInt> &QLocalVtx,
+			   vector<MilanLongInt> &QGhostVtx,
+			   vector<MilanLongInt> &QMsgType,
+			   vector<MilanInt> &QOwner,
+			   vector<MPI_Request> &SRequest,
+			   vector<MPI_Status> &SStatus);
+  
+  void processMessagesD(
+		       MilanLongInt NLVer,
+		       MilanLongInt *Mate,
+		       MilanLongInt *candidateMate,
+		       map<MilanLongInt, MilanLongInt> &Ghost2LocalMap,
+		       vector<MilanLongInt> &GMate,
+		       vector<MilanLongInt> &Counter,
+		       MilanLongInt StartIndex,
+		       MilanLongInt EndIndex,
+		       MilanLongInt *myCardPtr,
+		       MilanLongInt *msgIndPtr,
+		       MilanLongInt *msgActualPtr,
+		       MilanReal *edgeLocWeight,
+		       MilanLongInt *verDistance,
+		       MilanLongInt *verLocPtr,
+		       MilanLongInt k,
+		       MilanLongInt *verLocInd,
+		       MilanInt numProcs,
+		       MilanInt myRank,
+		       MPI_Comm comm,
+		       vector<MilanLongInt> &Message,
+		       MilanLongInt numGhostEdges,
+		       MilanLongInt u,
+		       MilanLongInt v,
+		       MilanLongInt *SPtr,
+		       vector<MilanLongInt> &U);
+  
+  void extractUChunk(
+		     vector<MilanLongInt> &UChunkBeingProcessed,
+		     vector<MilanLongInt> &U,
+		     vector<MilanLongInt> &privateU);
 
-    void queuesTransfer(vector<MilanLongInt> &U,
-                        vector<MilanLongInt> &privateU,
-                        vector<MilanLongInt> &QLocalVtx,
-                        vector<MilanLongInt> &QGhostVtx,
-                        vector<MilanLongInt> &QMsgType,
-                        vector<MilanInt> &QOwner,
-                        vector<MilanLongInt> &privateQLocalVtx,
-                        vector<MilanLongInt> &privateQGhostVtx,
-                        vector<MilanLongInt> &privateQMsgType,
-                        vector<MilanInt> &privateQOwner);
+  MilanLongInt firstComputeCandidateMateS(MilanLongInt adj1,
+					 MilanLongInt adj2,
+					 MilanLongInt *verLocInd,
+					  MilanFloat *edgeLocWeight);
 
-    bool isAlreadyMatched(MilanLongInt node,
-                          MilanLongInt StartIndex,
-                          MilanLongInt EndIndex,
-                          vector<MilanLongInt> &GMate,
-                          MilanLongInt *Mate,
-                          map<MilanLongInt, MilanLongInt> &Ghost2LocalMap);
+  MilanLongInt computeCandidateMateS(MilanLongInt adj1,
+				    MilanLongInt adj2,
+				    MilanFloat *edgeLocWeight,
+				    MilanLongInt k,
+				    MilanLongInt *verLocInd,
+				    MilanLongInt StartIndex,
+				    MilanLongInt EndIndex,
+				    vector<MilanLongInt> &GMate,
+				    MilanLongInt *Mate,
+				    map<MilanLongInt, MilanLongInt> &Ghost2LocalMap);
+  
+  void PARALLEL_COMPUTE_CANDIDATE_MATE_BS(MilanLongInt NLVer,
+					 MilanLongInt *verLocPtr,
+					 MilanLongInt *verLocInd,
+					 MilanInt myRank,
+					 MilanFloat *edgeLocWeight,
+					 MilanLongInt *candidateMate);
+  
+  void PARALLEL_PROCESS_EXPOSED_VERTEX_BS(MilanLongInt NLVer,
+					 MilanLongInt *candidateMate,
+					 MilanLongInt *verLocInd,
+					 MilanLongInt *verLocPtr,
+					 MilanLongInt StartIndex,
+					 MilanLongInt EndIndex,
+					 MilanLongInt *Mate,
+					 vector<MilanLongInt> &GMate,
+					 map<MilanLongInt, MilanLongInt> &Ghost2LocalMap,
+					 MilanFloat *edgeLocWeight,
+					 MilanLongInt *myCardPtr,
+					 MilanLongInt *msgIndPtr,
+					 MilanLongInt *NumMessagesBundledPtr,
+					 MilanLongInt *SPtr,
+					 MilanLongInt *verDistance,
+					 MilanLongInt *PCounter,
+					 vector<MilanLongInt> &Counter,
+					 MilanInt myRank,
+					 MilanInt numProcs,
+					 vector<MilanLongInt> &U,
+					 vector<MilanLongInt> &privateU,
+					 vector<MilanLongInt> &QLocalVtx,
+					 vector<MilanLongInt> &QGhostVtx,
+					 vector<MilanLongInt> &QMsgType,
+					 vector<MilanInt> &QOwner,
+					 vector<MilanLongInt> &privateQLocalVtx,
+					 vector<MilanLongInt> &privateQGhostVtx,
+					 vector<MilanLongInt> &privateQMsgType,
+					 vector<MilanInt> &privateQOwner);
+  void processMatchedVerticesS(
+			      MilanLongInt NLVer,
+			      vector<MilanLongInt> &UChunkBeingProcessed,
+			      vector<MilanLongInt> &U,
+			      vector<MilanLongInt> &privateU,
+			      MilanLongInt StartIndex,
+			      MilanLongInt EndIndex,
+			      MilanLongInt *myCardPtr,
+			      MilanLongInt *msgIndPtr,
+			      MilanLongInt *NumMessagesBundledPtr,
+			      MilanLongInt *SPtr,
+			      MilanLongInt *verLocPtr,
+			      MilanLongInt *verLocInd,
+			      MilanLongInt *verDistance,
+			      MilanLongInt *PCounter,
+			      vector<MilanLongInt> &Counter,
+			      MilanInt myRank,
+			      MilanInt numProcs,
+			      MilanLongInt *candidateMate,
+			      vector<MilanLongInt> &GMate,
+			      MilanLongInt *Mate,
+			      map<MilanLongInt, MilanLongInt> &Ghost2LocalMap,
+			      MilanFloat *edgeLocWeight,
+			      vector<MilanLongInt> &QLocalVtx,
+			      vector<MilanLongInt> &QGhostVtx,
+			      vector<MilanLongInt> &QMsgType,
+			      vector<MilanInt> &QOwner,
+			      vector<MilanLongInt> &privateQLocalVtx,
+			      vector<MilanLongInt> &privateQGhostVtx,
+			      vector<MilanLongInt> &privateQMsgType,
+			      vector<MilanInt> &privateQOwner);
+  
+  void processMatchedVerticesAndSendMessagesS(
+					     MilanLongInt NLVer,
+					     vector<MilanLongInt> &UChunkBeingProcessed,
+					     vector<MilanLongInt> &U,
+					     vector<MilanLongInt> &privateU,
+					     MilanLongInt StartIndex,
+					     MilanLongInt EndIndex,
+					     MilanLongInt *myCardPtr,
+					     MilanLongInt *msgIndPtr,
+					     MilanLongInt *NumMessagesBundledPtr,
+					     MilanLongInt *SPtr,
+					     MilanLongInt *verLocPtr,
+					     MilanLongInt *verLocInd,
+					     MilanLongInt *verDistance,
+					     MilanLongInt *PCounter,
+					     vector<MilanLongInt> &Counter,
+					     MilanInt myRank,
+					     MilanInt numProcs,
+					     MilanLongInt *candidateMate,
+					     vector<MilanLongInt> &GMate,
+					     MilanLongInt *Mate,
+					     map<MilanLongInt, MilanLongInt> &Ghost2LocalMap,
+					     MilanFloat *edgeLocWeight,
+					     vector<MilanLongInt> &QLocalVtx,
+					     vector<MilanLongInt> &QGhostVtx,
+					     vector<MilanLongInt> &QMsgType,
+					     vector<MilanInt> &QOwner,
+					     vector<MilanLongInt> &privateQLocalVtx,
+					     vector<MilanLongInt> &privateQGhostVtx,
+					     vector<MilanLongInt> &privateQMsgType,
+					     vector<MilanInt> &privateQOwner,
+					     MPI_Comm comm,
+					     MilanLongInt *msgActual,
+					     vector<MilanLongInt> &Message);
+  
+  void processMessagesS(
+		       MilanLongInt NLVer,
+		       MilanLongInt *Mate,
+		       MilanLongInt *candidateMate,
+		       map<MilanLongInt, MilanLongInt> &Ghost2LocalMap,
+		       vector<MilanLongInt> &GMate,
+		       vector<MilanLongInt> &Counter,
+		       MilanLongInt StartIndex,
+		       MilanLongInt EndIndex,
+		       MilanLongInt *myCardPtr,
+		       MilanLongInt *msgIndPtr,
+		       MilanLongInt *msgActualPtr,
+		       MilanFloat *edgeLocWeight,
+		       MilanLongInt *verDistance,
+		       MilanLongInt *verLocPtr,
+		       MilanLongInt k,
+		       MilanLongInt *verLocInd,
+		       MilanInt numProcs,
+		       MilanInt myRank,
+		       MPI_Comm comm,
+		       vector<MilanLongInt> &Message,
+		       MilanLongInt numGhostEdges,
+		       MilanLongInt u,
+		       MilanLongInt v,
+		       MilanLongInt *SPtr,
+		       vector<MilanLongInt> &U);
+  
+  
+  void dalgoDistEdgeApproxDomEdgesLinearSearchMesgBndlSmallMateCMP(
+			 MilanLongInt NLVer, MilanLongInt NLEdge,
+			 MilanLongInt *verLocPtr, MilanLongInt *verLocInd, MilanReal *edgeLocWeight,
+			 MilanLongInt *verDistance,
+			 MilanLongInt *Mate,
+			 MilanInt myRank, MilanInt numProcs, MPI_Comm comm,
+			 MilanLongInt *msgIndSent, MilanLongInt *msgActualSent, MilanReal *msgPercent,
+			 MilanReal *ph0_time, MilanReal *ph1_time, MilanReal *ph2_time,
+			 MilanLongInt *ph1_card, MilanLongInt *ph2_card);
 
-    MilanLongInt computeCandidateMate(MilanLongInt adj1,
-                                      MilanLongInt adj2,
-                                      MilanReal *edgeLocWeight,
-                                      MilanLongInt k,
-                                      MilanLongInt *verLocInd,
-                                      MilanLongInt StartIndex,
-                                      MilanLongInt EndIndex,
-                                      vector<MilanLongInt> &GMate,
-                                      MilanLongInt *Mate,
-                                      map<MilanLongInt, MilanLongInt> &Ghost2LocalMap);
-
-    void initialize(MilanLongInt NLVer, MilanLongInt NLEdge,
-                    MilanLongInt StartIndex, MilanLongInt EndIndex,
-                    MilanLongInt *numGhostEdgesPtr,
-                    MilanLongInt *numGhostVerticesPtr,
-                    MilanLongInt *S,
-                    MilanLongInt *verLocInd,
-                    MilanLongInt *verLocPtr,
-                    map<MilanLongInt, MilanLongInt> &Ghost2LocalMap,
-                    vector<MilanLongInt> &Counter,
-                    vector<MilanLongInt> &verGhostPtr,
-                    vector<MilanLongInt> &verGhostInd,
-                    vector<MilanLongInt> &tempCounter,
-                    vector<MilanLongInt> &GMate,
-                    vector<MilanLongInt> &Message,
-                    vector<MilanLongInt> &QLocalVtx,
-                    vector<MilanLongInt> &QGhostVtx,
-                    vector<MilanLongInt> &QMsgType,
-                    vector<MilanInt> &QOwner,
-                    MilanLongInt *&candidateMate,
-                    vector<MilanLongInt> &U,
-                    vector<MilanLongInt> &privateU,
-                    vector<MilanLongInt> &privateQLocalVtx,
-                    vector<MilanLongInt> &privateQGhostVtx,
-                    vector<MilanLongInt> &privateQMsgType,
-                    vector<MilanInt> &privateQOwner);
-
-    void clean(MilanLongInt NLVer,
-               MilanInt myRank,
-               MilanLongInt MessageIndex,
-               vector<MPI_Request> &SRequest,
-               vector<MPI_Status> &SStatus,
-               MilanInt BufferSize,
-               MilanLongInt *Buffer,
-               MilanLongInt msgActual,
-               MilanLongInt *msgActualSent,
-               MilanLongInt msgInd,
-               MilanLongInt *msgIndSent,
-               MilanLongInt NumMessagesBundled,
-               MilanReal *msgPercent);
-
-    void PARALLEL_COMPUTE_CANDIDATE_MATE_B(MilanLongInt NLVer,
-                                           MilanLongInt *verLocPtr,
-                                           MilanLongInt *verLocInd,
-                                           MilanInt myRank,
-                                           MilanReal *edgeLocWeight,
-                                           MilanLongInt *candidateMate);
-
-    void PARALLEL_PROCESS_EXPOSED_VERTEX_B(MilanLongInt NLVer,
-                                           MilanLongInt *candidateMate,
-                                           MilanLongInt *verLocInd,
-                                           MilanLongInt *verLocPtr,
-                                           MilanLongInt StartIndex,
-                                           MilanLongInt EndIndex,
-                                           MilanLongInt *Mate,
-                                           vector<MilanLongInt> &GMate,
-                                           map<MilanLongInt, MilanLongInt> &Ghost2LocalMap,
-                                           MilanReal *edgeLocWeight,
-                                           MilanLongInt *myCardPtr,
-                                           MilanLongInt *msgIndPtr,
-                                           MilanLongInt *NumMessagesBundledPtr,
-                                           MilanLongInt *SPtr,
-                                           MilanLongInt *verDistance,
-                                           MilanLongInt *PCounter,
-                                           vector<MilanLongInt> &Counter,
-                                           MilanInt myRank,
-                                           MilanInt numProcs,
-                                           vector<MilanLongInt> &U,
-                                           vector<MilanLongInt> &privateU,
-                                           vector<MilanLongInt> &QLocalVtx,
-                                           vector<MilanLongInt> &QGhostVtx,
-                                           vector<MilanLongInt> &QMsgType,
-                                           vector<MilanInt> &QOwner,
-                                           vector<MilanLongInt> &privateQLocalVtx,
-                                           vector<MilanLongInt> &privateQGhostVtx,
-                                           vector<MilanLongInt> &privateQMsgType,
-                                           vector<MilanInt> &privateQOwner);
-
-    void PROCESS_CROSS_EDGE(MilanLongInt *edge,
-                            MilanLongInt *SPtr);
-
-    void processMatchedVertices(
-        MilanLongInt NLVer,
-        vector<MilanLongInt> &UChunkBeingProcessed,
-        vector<MilanLongInt> &U,
-        vector<MilanLongInt> &privateU,
-        MilanLongInt StartIndex,
-        MilanLongInt EndIndex,
-        MilanLongInt *myCardPtr,
-        MilanLongInt *msgIndPtr,
-        MilanLongInt *NumMessagesBundledPtr,
-        MilanLongInt *SPtr,
-        MilanLongInt *verLocPtr,
-        MilanLongInt *verLocInd,
-        MilanLongInt *verDistance,
-        MilanLongInt *PCounter,
-        vector<MilanLongInt> &Counter,
-        MilanInt myRank,
-        MilanInt numProcs,
-        MilanLongInt *candidateMate,
-        vector<MilanLongInt> &GMate,
-        MilanLongInt *Mate,
-        map<MilanLongInt, MilanLongInt> &Ghost2LocalMap,
-        MilanReal *edgeLocWeight,
-        vector<MilanLongInt> &QLocalVtx,
-        vector<MilanLongInt> &QGhostVtx,
-        vector<MilanLongInt> &QMsgType,
-        vector<MilanInt> &QOwner,
-        vector<MilanLongInt> &privateQLocalVtx,
-        vector<MilanLongInt> &privateQGhostVtx,
-        vector<MilanLongInt> &privateQMsgType,
-        vector<MilanInt> &privateQOwner);
-
-    void processMatchedVerticesAndSendMessages(
-        MilanLongInt NLVer,
-        vector<MilanLongInt> &UChunkBeingProcessed,
-        vector<MilanLongInt> &U,
-        vector<MilanLongInt> &privateU,
-        MilanLongInt StartIndex,
-        MilanLongInt EndIndex,
-        MilanLongInt *myCardPtr,
-        MilanLongInt *msgIndPtr,
-        MilanLongInt *NumMessagesBundledPtr,
-        MilanLongInt *SPtr,
-        MilanLongInt *verLocPtr,
-        MilanLongInt *verLocInd,
-        MilanLongInt *verDistance,
-        MilanLongInt *PCounter,
-        vector<MilanLongInt> &Counter,
-        MilanInt myRank,
-        MilanInt numProcs,
-        MilanLongInt *candidateMate,
-        vector<MilanLongInt> &GMate,
-        MilanLongInt *Mate,
-        map<MilanLongInt, MilanLongInt> &Ghost2LocalMap,
-        MilanReal *edgeLocWeight,
-        vector<MilanLongInt> &QLocalVtx,
-        vector<MilanLongInt> &QGhostVtx,
-        vector<MilanLongInt> &QMsgType,
-        vector<MilanInt> &QOwner,
-        vector<MilanLongInt> &privateQLocalVtx,
-        vector<MilanLongInt> &privateQGhostVtx,
-        vector<MilanLongInt> &privateQMsgType,
-        vector<MilanInt> &privateQOwner,
-        MPI_Comm comm,
-        MilanLongInt *msgActual,
-        vector<MilanLongInt> &Message);
-
-    void sendBundledMessages(MilanLongInt *numGhostEdgesPtr,
-                             MilanInt *BufferSizePtr,
-                             MilanLongInt *Buffer,
-                             vector<MilanLongInt> &PCumulative,
-                             vector<MilanLongInt> &PMessageBundle,
-                             vector<MilanLongInt> &PSizeInfoMessages,
-                             MilanLongInt *PCounter,
-                             MilanLongInt NumMessagesBundled,
-                             MilanLongInt *msgActualPtr,
-                             MilanLongInt *MessageIndexPtr,
-                             MilanInt numProcs,
-                             MilanInt myRank,
-                             MPI_Comm comm,
-                             vector<MilanLongInt> &QLocalVtx,
-                             vector<MilanLongInt> &QGhostVtx,
-                             vector<MilanLongInt> &QMsgType,
-                             vector<MilanInt> &QOwner,
-                             vector<MPI_Request> &SRequest,
-                             vector<MPI_Status> &SStatus);
-
-    void processMessages(
-        MilanLongInt NLVer,
-        MilanLongInt *Mate,
-        MilanLongInt *candidateMate,
-        map<MilanLongInt, MilanLongInt> &Ghost2LocalMap,
-        vector<MilanLongInt> &GMate,
-        vector<MilanLongInt> &Counter,
-        MilanLongInt StartIndex,
-        MilanLongInt EndIndex,
-        MilanLongInt *myCardPtr,
-        MilanLongInt *msgIndPtr,
-        MilanLongInt *msgActualPtr,
-        MilanReal *edgeLocWeight,
-        MilanLongInt *verDistance,
-        MilanLongInt *verLocPtr,
-        MilanLongInt k,
-        MilanLongInt *verLocInd,
-        MilanInt numProcs,
-        MilanInt myRank,
-        MPI_Comm comm,
-        vector<MilanLongInt> &Message,
-        MilanLongInt numGhostEdges,
-        MilanLongInt u,
-        MilanLongInt v,
-        MilanLongInt *SPtr,
-        vector<MilanLongInt> &U);
-
-    void extractUChunk(
-        vector<MilanLongInt> &UChunkBeingProcessed,
-        vector<MilanLongInt> &U,
-        vector<MilanLongInt> &privateU);
-
-    void dalgoDistEdgeApproxDomEdgesLinearSearchMesgBndlSmallMateCMP(
-        MilanLongInt NLVer, MilanLongInt NLEdge,
-        MilanLongInt *verLocPtr, MilanLongInt *verLocInd, MilanReal *edgeLocWeight,
-        MilanLongInt *verDistance,
-        MilanLongInt *Mate,
-        MilanInt myRank, MilanInt numProcs, MPI_Comm comm,
-        MilanLongInt *msgIndSent, MilanLongInt *msgActualSent, MilanReal *msgPercent,
-        MilanReal *ph0_time, MilanReal *ph1_time, MilanReal *ph2_time,
-        MilanLongInt *ph1_card, MilanLongInt *ph2_card);
-#endif
+  
+  void salgoDistEdgeApproxDomEdgesLinearSearchMesgBndlSmallMateCMP(
+			 MilanLongInt NLVer, MilanLongInt NLEdge,
+			 MilanLongInt *verLocPtr, MilanLongInt *verLocInd, MilanFloat *edgeLocWeight,
+			 MilanLongInt *verDistance,
+			 MilanLongInt *Mate,
+			 MilanInt myRank, MilanInt numProcs, MPI_Comm comm,
+			 MilanLongInt *msgIndSent, MilanLongInt *msgActualSent, MilanReal *msgPercent,
+			 MilanReal *ph0_time, MilanReal *ph1_time, MilanReal *ph2_time,
+			 MilanLongInt *ph1_card, MilanLongInt *ph2_card);
 
 
-#ifndef OMP
-   //Function of find the owner of a ghost vertex using binary search:
-   inline MilanInt findOwnerOfGhost(MilanLongInt vtxIndex, MilanLongInt *mVerDistance,
-                                        MilanInt myRank, MilanInt numProcs);
-#endif
 
     void dalgoDistEdgeApproxDomEdgesLinearSearchMesgBndlSmallMateC(
         MilanLongInt NLVer, MilanLongInt NLEdge,

amgprec/impl/aggregator/algoDistEdgeApproxDomEdgesLinearSearchMesgBndlSmallMateC.cpp

@@ -1303,16 +1303,16 @@ void dalgoDistEdgeApproxDomEdgesLinearSearchMesgBndlSmallMateC(
 // SINGLE PRECISION VERSION
 
 void salgoDistEdgeApproxDomEdgesLinearSearchMesgBndlSmallMateC(
-                                                              MilanLongInt NLVer, MilanLongInt NLEdge,
-                                                              MilanLongInt* verLocPtr, MilanLongInt* verLocInd,
-							      MilanFloat* edgeLocWeight,
-                                                              MilanLongInt* verDistance,
-                                                              MilanLongInt* Mate,
-                                                              MilanInt myRank, MilanInt numProcs, MPI_Comm comm,
-                                                              MilanLongInt* msgIndSent, MilanLongInt* msgActualSent,
-							      MilanReal* msgPercent,
-                                                              MilanReal* ph0_time, MilanReal* ph1_time, MilanReal* ph2_time,
-                                                              MilanLongInt* ph1_card, MilanLongInt* ph2_card ) {
+      MilanLongInt NLVer, MilanLongInt NLEdge,
+      MilanLongInt* verLocPtr, MilanLongInt* verLocInd,
+      MilanFloat* edgeLocWeight,
+      MilanLongInt* verDistance,
+      MilanLongInt* Mate,
+      MilanInt myRank, MilanInt numProcs, MPI_Comm comm,
+      MilanLongInt* msgIndSent, MilanLongInt* msgActualSent,
+      MilanReal* msgPercent,
+      MilanReal* ph0_time, MilanReal* ph1_time, MilanReal* ph2_time,
+      MilanLongInt* ph1_card, MilanLongInt* ph2_card ) {
 #if !defined(SERIAL_MPI)
 #ifdef PRINT_DEBUG_INFO_
     cout<<"\n("<<myRank<<")Within algoEdgeApproxDominatingEdgesLinearSearchMessageBundling()"; fflush(stdout);

amgprec/impl/aggregator/algoDistEdgeApproxDomEdgesLinearSearchMesgBndlSmallMateCMP.cpp

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

amgprec/impl/aggregator/amg_c_soc1_map_bld.F90

@@ -275,7 +275,8 @@ subroutine amg_c_soc1_map_bld(iorder,theta,clean_zeros,a,desc_a,nlaggr,ilaggr,in
             disjoint = all(ilaggr(icol(1:ip)) == -(nr+1)).or.(ip==0)
             if (disjoint) then
               locnaggr(kk)     = locnaggr(kk) + 1
-              itmp = (bnds(kk)-1+locnaggr(kk))*nths+kk
+              itmp = (bnds(kk)-1+locnaggr(kk)) !be careful about overflow
+              itmp = itmp*nths+kk
               if (itmp < (bnds(kk)-1+locnaggr(kk))) then
                 !$omp atomic update
                 info = max(12345678,info)

amgprec/impl/aggregator/amg_c_soc2_map_bld.F90

@@ -309,7 +309,8 @@ subroutine amg_c_soc2_map_bld(iorder,theta,clean_zeros,a,desc_a,nlaggr,ilaggr,in
             !
             if (disjoint) then
               locnaggr(kk)     = locnaggr(kk) + 1
-              itmp = (bnds(kk)-1+locnaggr(kk))*nths+kk
+              itmp = (bnds(kk)-1+locnaggr(kk)) !be careful about overflow
+              itmp = itmp*nths+kk
               if (itmp < (bnds(kk)-1+locnaggr(kk))) then
                 !$omp atomic update
                 info = max(12345678,info)

amgprec/impl/aggregator/amg_caggrmat_minnrg_bld.f90

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

amgprec/impl/aggregator/amg_caggrmat_smth_bld.f90

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

amgprec/impl/aggregator/amg_d_parmatch_smth_bld.F90

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

amgprec/impl/aggregator/amg_d_soc1_map_bld.F90

@@ -275,7 +275,8 @@ subroutine amg_d_soc1_map_bld(iorder,theta,clean_zeros,a,desc_a,nlaggr,ilaggr,in
             disjoint = all(ilaggr(icol(1:ip)) == -(nr+1)).or.(ip==0)
             if (disjoint) then
               locnaggr(kk)     = locnaggr(kk) + 1
-              itmp = (bnds(kk)-1+locnaggr(kk))*nths+kk
+              itmp = (bnds(kk)-1+locnaggr(kk)) !be careful about overflow
+              itmp = itmp*nths+kk
               if (itmp < (bnds(kk)-1+locnaggr(kk))) then
                 !$omp atomic update
                 info = max(12345678,info)

amgprec/impl/aggregator/amg_d_soc2_map_bld.F90

@@ -309,7 +309,8 @@ subroutine amg_d_soc2_map_bld(iorder,theta,clean_zeros,a,desc_a,nlaggr,ilaggr,in
             !
             if (disjoint) then
               locnaggr(kk)     = locnaggr(kk) + 1
-              itmp = (bnds(kk)-1+locnaggr(kk))*nths+kk
+              itmp = (bnds(kk)-1+locnaggr(kk)) !be careful about overflow
+              itmp = itmp*nths+kk
               if (itmp < (bnds(kk)-1+locnaggr(kk))) then
                 !$omp atomic update
                 info = max(12345678,info)

amgprec/impl/aggregator/amg_daggrmat_minnrg_bld.f90

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

amgprec/impl/aggregator/amg_daggrmat_smth_bld.f90

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

amgprec/impl/aggregator/amg_s_parmatch_smth_bld.F90

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

amgprec/impl/aggregator/amg_s_soc1_map_bld.F90

@@ -275,7 +275,8 @@ subroutine amg_s_soc1_map_bld(iorder,theta,clean_zeros,a,desc_a,nlaggr,ilaggr,in
             disjoint = all(ilaggr(icol(1:ip)) == -(nr+1)).or.(ip==0)
             if (disjoint) then
               locnaggr(kk)     = locnaggr(kk) + 1
-              itmp = (bnds(kk)-1+locnaggr(kk))*nths+kk
+              itmp = (bnds(kk)-1+locnaggr(kk)) !be careful about overflow
+              itmp = itmp*nths+kk
               if (itmp < (bnds(kk)-1+locnaggr(kk))) then
                 !$omp atomic update
                 info = max(12345678,info)

amgprec/impl/aggregator/amg_s_soc2_map_bld.F90

@@ -309,7 +309,8 @@ subroutine amg_s_soc2_map_bld(iorder,theta,clean_zeros,a,desc_a,nlaggr,ilaggr,in
             !
             if (disjoint) then
               locnaggr(kk)     = locnaggr(kk) + 1
-              itmp = (bnds(kk)-1+locnaggr(kk))*nths+kk
+              itmp = (bnds(kk)-1+locnaggr(kk)) !be careful about overflow
+              itmp = itmp*nths+kk
               if (itmp < (bnds(kk)-1+locnaggr(kk))) then
                 !$omp atomic update
                 info = max(12345678,info)

amgprec/impl/aggregator/amg_saggrmat_minnrg_bld.f90

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

amgprec/impl/aggregator/amg_saggrmat_smth_bld.f90

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

amgprec/impl/aggregator/amg_z_soc1_map_bld.F90

@@ -275,7 +275,8 @@ subroutine amg_z_soc1_map_bld(iorder,theta,clean_zeros,a,desc_a,nlaggr,ilaggr,in
             disjoint = all(ilaggr(icol(1:ip)) == -(nr+1)).or.(ip==0)
             if (disjoint) then
               locnaggr(kk)     = locnaggr(kk) + 1
-              itmp = (bnds(kk)-1+locnaggr(kk))*nths+kk
+              itmp = (bnds(kk)-1+locnaggr(kk)) !be careful about overflow
+              itmp = itmp*nths+kk
               if (itmp < (bnds(kk)-1+locnaggr(kk))) then
                 !$omp atomic update
                 info = max(12345678,info)

amgprec/impl/aggregator/amg_z_soc2_map_bld.F90

@@ -309,7 +309,8 @@ subroutine amg_z_soc2_map_bld(iorder,theta,clean_zeros,a,desc_a,nlaggr,ilaggr,in
             !
             if (disjoint) then
               locnaggr(kk)     = locnaggr(kk) + 1
-              itmp = (bnds(kk)-1+locnaggr(kk))*nths+kk
+              itmp = (bnds(kk)-1+locnaggr(kk)) !be careful about overflow
+              itmp = itmp*nths+kk
               if (itmp < (bnds(kk)-1+locnaggr(kk))) then
                 !$omp atomic update
                 info = max(12345678,info)

amgprec/impl/aggregator/amg_zaggrmat_minnrg_bld.f90

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

amgprec/impl/aggregator/amg_zaggrmat_smth_bld.f90

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

amgprec/impl/aggregator/clean.cpp

@@ -1,6 +1,6 @@
 #include "MatchBoxPC.h"
-#ifdef OMP
 // TODO comment
+#if !defined(SERIAL_MPI)
 
 void clean(MilanLongInt NLVer,
            MilanInt myRank,

amgprec/impl/aggregator/computeCandidateMate.cpp

@@ -1,5 +1,4 @@
 #include "MatchBoxPC.h"
-#ifdef OMP
 /**
  * Execute the research fr the Candidate Mate without controlling if the vertices are already matched.
  * Returns the vertices with the highest weight
@@ -9,7 +8,8 @@
  * @param edgeLocWeight
  * @return
  */
-MilanLongInt firstComputeCandidateMate(MilanLongInt adj1,
+
+MilanLongInt firstComputeCandidateMateD(MilanLongInt adj1,
                                        MilanLongInt adj2,
                                        MilanLongInt *verLocInd,
                                        MilanReal *edgeLocWeight)
@@ -42,7 +42,7 @@ MilanLongInt firstComputeCandidateMate(MilanLongInt adj1,
  * @param Ghost2LocalMap
  * @return
  */
-MilanLongInt computeCandidateMate(MilanLongInt adj1,
+MilanLongInt computeCandidateMateD(MilanLongInt adj1,
                                   MilanLongInt adj2,
                                   MilanReal *edgeLocWeight,
                                   MilanLongInt k,
@@ -71,4 +71,68 @@ MilanLongInt computeCandidateMate(MilanLongInt adj1,
 
     return w;
 }
-#endif
+
+
+MilanLongInt firstComputeCandidateMateS(MilanLongInt adj1,
+                                       MilanLongInt adj2,
+                                       MilanLongInt *verLocInd,
+                                       MilanFloat *edgeLocWeight)
+{
+    MilanInt w = -1;
+    MilanFloat heaviestEdgeWt = 0.0f; // Assign the smallest 
+    int finalK;
+    for (int k = adj1; k < adj2; k++) {
+      if ((edgeLocWeight[k] > heaviestEdgeWt) ||
+	  ((edgeLocWeight[k] == heaviestEdgeWt) && (w < verLocInd[k])))     {
+	heaviestEdgeWt = edgeLocWeight[k];
+	w = verLocInd[k];
+	finalK = k;
+      }
+    } // End of for loop
+    return finalK;
+}
+
+/**
+ * //TODO documentation
+ * @param adj1
+ * @param adj2
+ * @param edgeLocWeight
+ * @param k
+ * @param verLocInd
+ * @param StartIndex
+ * @param EndIndex
+ * @param GMate
+ * @param Mate
+ * @param Ghost2LocalMap
+ * @return
+ */
+MilanLongInt computeCandidateMateS(MilanLongInt adj1,
+                                  MilanLongInt adj2,
+				  MilanFloat *edgeLocWeight,
+                                  MilanLongInt k,
+                                  MilanLongInt *verLocInd,
+                                  MilanLongInt StartIndex,
+                                  MilanLongInt EndIndex,
+                                  vector<MilanLongInt> &GMate,
+                                  MilanLongInt *Mate,
+                                  map<MilanLongInt, MilanLongInt> &Ghost2LocalMap)
+{
+    // Start: PARALLEL_COMPUTE_CANDIDATE_MATE_B(v)
+
+    MilanInt w = -1;
+    MilanFloat heaviestEdgeWt = 0.0f; // Assign the smallest Value
+    for (k = adj1; k < adj2; k++)   {
+      if (isAlreadyMatched(verLocInd[k], StartIndex, EndIndex, GMate, Mate, Ghost2LocalMap))
+	continue;
+
+      if ((edgeLocWeight[k] > heaviestEdgeWt) ||
+	  ((edgeLocWeight[k] == heaviestEdgeWt) && (w < verLocInd[k])))     {
+	heaviestEdgeWt = edgeLocWeight[k];
+	w = verLocInd[k];
+      }
+    } // End of for loop
+      //  End: PARALLEL_COMPUTE_CANDIDATE_MATE_B(v)
+
+    return w;
+}
+

amgprec/impl/aggregator/extractUChunk.cpp

@@ -1,5 +1,4 @@
 #include "MatchBoxPC.h"
-#ifdef OMP
 void extractUChunk(
     vector<MilanLongInt> &UChunkBeingProcessed,
     vector<MilanLongInt> &U,
@@ -29,4 +28,3 @@ void extractUChunk(
 
     } // End of critical U // End of critical U
 }
-#endif

amgprec/impl/aggregator/findOwnerOfGhost.cpp

@@ -1,5 +1,4 @@
 #include "MatchBoxPC.h"
-#ifdef OMP
 /// Find the owner of a ghost node:
 MilanInt findOwnerOfGhost(MilanLongInt vtxIndex, MilanLongInt *mVerDistance,
                           MilanInt myRank, MilanInt numProcs)
@@ -27,4 +26,3 @@ MilanInt findOwnerOfGhost(MilanLongInt vtxIndex, MilanLongInt *mVerDistance,
 
   return Current;
 } // End of findOwnerOfGhost()
-#endif

amgprec/impl/aggregator/initialize.cpp

@@ -1,5 +1,4 @@
 #include "MatchBoxPC.h"
-#ifdef OMP
 void initialize(MilanLongInt NLVer, MilanLongInt NLEdge,
                 MilanLongInt StartIndex, MilanLongInt EndIndex,
                 MilanLongInt *numGhostEdges,
@@ -302,4 +301,3 @@ void initialize(MilanLongInt NLVer, MilanLongInt NLEdge,
         } // End of single region
     }     // End of parallel region
 }
-#endif

amgprec/impl/aggregator/isAlreadyMatched.cpp

@@ -1,5 +1,4 @@
 #include "MatchBoxPC.h"
-#ifdef OMP
 /**
  * //TODO documentation
  * @param k
@@ -44,4 +43,3 @@ bool isAlreadyMatched(MilanLongInt node,
 
     return val >= 0; // Already matched
 }
-#endif

amgprec/impl/aggregator/parallelComputeCandidateMateB.cpp

@@ -1,6 +1,6 @@
 #include "MatchBoxPC.h"
-#ifdef OMP
-void PARALLEL_COMPUTE_CANDIDATE_MATE_B(MilanLongInt NLVer,
+
+void PARALLEL_COMPUTE_CANDIDATE_MATE_BD(MilanLongInt NLVer,
                                               MilanLongInt *verLocPtr,
                                               MilanLongInt *verLocInd,
                                               MilanInt myRank,
@@ -20,9 +20,37 @@ void PARALLEL_COMPUTE_CANDIDATE_MATE_B(MilanLongInt NLVer,
             fflush(stdout);
 #endif
             // Start: PARALLEL_COMPUTE_CANDIDATE_MATE_B(v)
-            candidateMate[v] = firstComputeCandidateMate(verLocPtr[v], verLocPtr[v + 1], verLocInd, edgeLocWeight);
+            candidateMate[v] = firstComputeCandidateMateD(verLocPtr[v], verLocPtr[v + 1],
+							 verLocInd, edgeLocWeight);
             // End: PARALLEL_COMPUTE_CANDIDATE_MATE_B(v)
         }
     }
 }
+
+void PARALLEL_COMPUTE_CANDIDATE_MATE_BS(MilanLongInt NLVer,
+                                              MilanLongInt *verLocPtr,
+                                              MilanLongInt *verLocInd,
+                                              MilanInt myRank,
+                                              MilanFloat *edgeLocWeight,
+                                              MilanLongInt *candidateMate)
+{
+
+    MilanLongInt v = -1;
+
+#pragma omp parallel private(v) default(shared) num_threads(NUM_THREAD)
+    {
+
+#pragma omp for schedule(static)
+        for (v = 0; v < NLVer; v++) {
+#ifdef PRINT_DEBUG_INFO_
+            cout << "\n(" << myRank << ")Processing: " << v + StartIndex << endl;
+            fflush(stdout);
 #endif
+            // Start: PARALLEL_COMPUTE_CANDIDATE_MATE_B(v)
+            candidateMate[v] = firstComputeCandidateMateS(verLocPtr[v], verLocPtr[v + 1],
+							 verLocInd, edgeLocWeight);
+            // End: PARALLEL_COMPUTE_CANDIDATE_MATE_B(v)
+        }
+    }
+}
+

amgprec/impl/aggregator/processCrossEdge.cpp

@@ -1,5 +1,4 @@
 #include "MatchBoxPC.h"
-#ifdef OMP
 void PROCESS_CROSS_EDGE(MilanLongInt *edge,
                         MilanLongInt *S)
 {
@@ -22,4 +21,3 @@ void PROCESS_CROSS_EDGE(MilanLongInt *edge,
 
     // End: PARALLEL_PROCESS_CROSS_EDGE_B
 }
-#endif

amgprec/impl/aggregator/processExposedVertex.cpp

@@ -1,6 +1,5 @@
-#include "MatchBoxPC.h"
-#ifdef OMP
-void PARALLEL_PROCESS_EXPOSED_VERTEX_B(MilanLongInt NLVer,
+#include "MatchBoxPC.h" 
+void PARALLEL_PROCESS_EXPOSED_VERTEX_BD(MilanLongInt NLVer,
                                        MilanLongInt *candidateMate,
                                        MilanLongInt *verLocInd,
                                        MilanLongInt *verLocPtr,
@@ -31,15 +30,16 @@ void PARALLEL_PROCESS_EXPOSED_VERTEX_B(MilanLongInt NLVer,
                                        vector<MilanInt> &privateQOwner)
 {
 
-    MilanLongInt v = -1, k = -1, w = -1, adj11 = 0, adj12 = 0, k1 = 0;
+  MilanLongInt v = -1, k = -1, w = -1, adj11 = 0, adj12 = 0, k1 = 0;
     MilanInt ghostOwner = 0, option, igw;
 
-#pragma omp parallel private(option, k, w, v, k1, adj11, adj12, ghostOwner)    \
-    firstprivate(privateU, StartIndex, EndIndex, privateQLocalVtx, privateQGhostVtx, privateQMsgType, privateQOwner) \
-     default(shared)   num_threads(NUM_THREAD)
+    //#pragma omp parallel private(option, k, w, v, k1, adj11, adj12, ghostOwner) \
+  firstprivate(privateU, StartIndex, EndIndex, privateQLocalVtx,	\
+	       privateQGhostVtx, privateQMsgType, privateQOwner)	\
+  default(shared)   num_threads(NUM_THREAD)
 
     {
-#pragma omp for reduction(+                             \
+      //#pragma omp for reduction(+			\
                           : PCounter[:numProcs], myCard \
                           [:1], msgInd                  \
                           [:1], NumMessagesBundled      \
@@ -62,102 +62,89 @@ void PARALLEL_PROCESS_EXPOSED_VERTEX_B(MilanLongInt NLVer,
             fflush(stdout);
 #endif
             // If found a dominating edge:
-            if (w >= 0)
-            {
-
-#pragma omp critical(processExposed)
-                {
-                    if (isAlreadyMatched(verLocInd[k], StartIndex, EndIndex, GMate, Mate, Ghost2LocalMap)) {
-		      w = computeCandidateMate(verLocPtr[v],
-					       verLocPtr[v + 1],
-					       edgeLocWeight, 0,
-					       verLocInd,
-					       StartIndex,
-					       EndIndex,
-					       GMate,
-					       Mate,
-					       Ghost2LocalMap);
-		      candidateMate[v] = w;
-                    }
-
-                    if (w >= 0) {
-		      (*myCard)++;
-		      if ((w < StartIndex) || (w > EndIndex)) { // w is a ghost vertex
-			option = 2;
-			if (candidateMate[NLVer + Ghost2LocalMap[w]] == v + StartIndex) {
-			  option = 1;
-			  Mate[v] = w;
-			  GMate[Ghost2LocalMap[w]] = v + StartIndex; // w is a Ghost
-
-			} // End of if CandidateMate[w] = v
-
-		      } // End of if a Ghost Vertex
-		      else  { // w is a local vertex
-
-			if (candidateMate[w - StartIndex] == (v + StartIndex)) {
-			  option = 3;
-			  Mate[v] = w;                           // v is local
-			  Mate[w - StartIndex] = v + StartIndex; // w is local
-
+            if (w >= 0) {
+#pragma omp critical(Matching)
+	      {
+		if (isAlreadyMatched(verLocInd[k], StartIndex, EndIndex, GMate, Mate, Ghost2LocalMap)) {
+		  w = computeCandidateMateD(verLocPtr[v], verLocPtr[v + 1], edgeLocWeight, 0,
+					   verLocInd, StartIndex, EndIndex,
+					   GMate,  Mate,  Ghost2LocalMap);
+		  candidateMate[v] = w;
+		}
+	      }
+	      if (w >= 0) {
+		(*myCard)++;
+		if ((w < StartIndex) || (w > EndIndex)) { // w is a ghost vertex
+		  option = 2;
+		  if (candidateMate[NLVer + Ghost2LocalMap[w]] == v + StartIndex) {
+		    option = 1;
+		    Mate[v] = w;
+		    GMate[Ghost2LocalMap[w]] = v + StartIndex; // w is a Ghost
+		  } // End of if CandidateMate[w] = v		  
+		} // End of if a Ghost Vertex
+		else  { // w is a local vertex
+		  if (candidateMate[w - StartIndex] == (v + StartIndex)) {
+		    option = 3;
+		    Mate[v] = w;                           // v is local
+		    Mate[w - StartIndex] = v + StartIndex; // w is local
 #ifdef PRINT_DEBUG_INFO_
-			  cout << "\n(" << myRank << ")MATCH: (" << v + StartIndex << "," << w << ") ";
-			  fflush(stdout);
+		    cout << "\n(" << myRank << ")MATCH: (" << v + StartIndex << "," << w << ") ";
+		    fflush(stdout);
 #endif
-
-			} // End of if ( candidateMate[w-StartIndex] == (v+StartIndex) )
-		      }     // End of Else
-
-                    } // End of second if
-
-                } // End critical processExposed
-
+		  } // End of if ( candidateMate[w-StartIndex] == (v+StartIndex) )
+		}     // End of Else
+	      } // End of second if
             } // End of if(w >=0)
             else  {
-	      // This piece of code is executed a really small amount of times
-	      adj11 = verLocPtr[v];
-	      adj12 = verLocPtr[v + 1];
-	      for (k1 = adj11; k1 < adj12; k1++) {
-		w = verLocInd[k1];
-		if ((w < StartIndex) || (w > EndIndex)) { // A ghost
-
+	      //#pragma omp critical(adjuse)
+	      {
+		// This piece of code is executed a really small number of times
+		adj11 = verLocPtr[v];
+		adj12 = verLocPtr[v + 1];
+		for (k1 = adj11; k1 < adj12; k1++) {
+		  w = verLocInd[k1];
+		  if ((w < StartIndex) || (w > EndIndex)) { // A ghost
+		    
 #ifdef PRINT_DEBUG_INFO_
-		  cout << "\n(" << myRank << ")Sending a failure message: ";
-		  cout << "\n(" << myRank << ")Ghost is " << w << " Owner is: " << findOwnerOfGhost(w, verDistance, myRank, numProcs);
-		  fflush(stdout);
+		    cout << "\n(" << myRank << ")Sending a failure message: ";
+		    cout << "\n(" << myRank << ")Ghost is " << w << " Owner is: " << findOwnerOfGhost(w, verDistance, myRank, numProcs);
+		    fflush(stdout);
 #endif
-		  (*msgInd)++;
-		  (*NumMessagesBundled)++;
-		  ghostOwner = findOwnerOfGhost(w, verDistance, myRank, numProcs);
-		  // assert(ghostOwner != -1);
-		  // assert(ghostOwner != myRank);
-		  PCounter[ghostOwner]++;
-
-		  privateQLocalVtx.push_back(v + StartIndex);
-		  privateQGhostVtx.push_back(w);
-		  privateQMsgType.push_back(FAILURE);
-		  privateQOwner.push_back(ghostOwner);
-
-		} // End of if(GHOST)
-	      }     // End of for loop
+		    (*msgInd)++;
+		    (*NumMessagesBundled)++;
+		    ghostOwner = findOwnerOfGhost(w, verDistance, myRank, numProcs);
+		    // assert(ghostOwner != -1);
+		    // assert(ghostOwner != myRank);
+#pragma omp atomic
+		    PCounter[ghostOwner]++;
+		    
+		    privateQLocalVtx.push_back(v + StartIndex);
+		    privateQGhostVtx.push_back(w);
+		    privateQMsgType.push_back(FAILURE);
+		    privateQOwner.push_back(ghostOwner);
+		    
+		  } // End of if(GHOST)
+		}     // End of for loop
+	      }
             }
             // End:   PARALLEL_PROCESS_EXPOSED_VERTEX_B(v)
-
+	    
             switch (option)
-            {
-            case -1:
+	      {
+	      case -1:
                 break;
-            case 1:
+	      case 1:
                 privateU.push_back(v + StartIndex);
                 privateU.push_back(w);
-
+		
 #ifdef PRINT_DEBUG_INFO_
                 cout << "\n(" << myRank << ")MATCH: (" << v + StartIndex << "," << w << ")";
                 fflush(stdout);
 #endif
-
+		
                 //  Decrement the counter:
                 PROCESS_CROSS_EDGE(&Counter[Ghost2LocalMap[w]], S);
-            case 2:
+	      case 2:
 #ifdef PRINT_DEBUG_INFO_
                 cout << "\n(" << myRank << ")Sending a request message (291):";
                 cout << "\n(" << myRank << ")Local is: " << v + StartIndex << " Ghost is " << w << " Owner is: " << findOwnerOfGhost(w, verDistance, myRank, numProcs) << endl;
@@ -168,29 +155,215 @@ void PARALLEL_PROCESS_EXPOSED_VERTEX_B(MilanLongInt NLVer,
                 ghostOwner = findOwnerOfGhost(w, verDistance, myRank, numProcs);
                 // assert(ghostOwner != -1);
                 // assert(ghostOwner != myRank);
+#pragma omp atomic
                 PCounter[ghostOwner]++;
-
+		
                 privateQLocalVtx.push_back(v + StartIndex);
                 privateQGhostVtx.push_back(w);
                 privateQMsgType.push_back(REQUEST);
                 privateQOwner.push_back(ghostOwner);
                 break;
-            case 3:
-            default:
+	      case 3:
+	      default:
                 privateU.push_back(v + StartIndex);
                 privateU.push_back(w);
                 break;
-            }
-
+	      }
+	    
         } // End of for ( v=0; v < NLVer; v++ )
-
+	
         queuesTransfer(U, privateU, QLocalVtx,
                        QGhostVtx,
                        QMsgType, QOwner, privateQLocalVtx,
                        privateQGhostVtx,
                        privateQMsgType,
                        privateQOwner);
-
+	
     } // End of parallel region
 }
+
+
+void PARALLEL_PROCESS_EXPOSED_VERTEX_BS(MilanLongInt NLVer,
+                                       MilanLongInt *candidateMate,
+                                       MilanLongInt *verLocInd,
+                                       MilanLongInt *verLocPtr,
+                                       MilanLongInt StartIndex,
+                                       MilanLongInt EndIndex,
+                                       MilanLongInt *Mate,
+                                       vector<MilanLongInt> &GMate,
+                                       map<MilanLongInt, MilanLongInt> &Ghost2LocalMap,
+                                       MilanFloat *edgeLocWeight,
+                                       MilanLongInt *myCard,
+                                       MilanLongInt *msgInd,
+                                       MilanLongInt *NumMessagesBundled,
+                                       MilanLongInt *S,
+                                       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)
+{
+
+    MilanLongInt v = -1, k = -1, w = -1, adj11 = 0, adj12 = 0, k1 = 0;
+    MilanInt ghostOwner = 0, option, igw;
+
+    //#pragma omp parallel private(option, k, w, v, k1, adj11, adj12, ghostOwner) \
+    firstprivate(privateU, StartIndex, EndIndex, privateQLocalVtx, \
+		 privateQGhostVtx, privateQMsgType, privateQOwner) \
+     default(shared)   num_threads(NUM_THREAD)
+
+    {
+      //#pragma omp for reduction(+			\
+                          : PCounter[:numProcs], myCard \
+                          [:1], msgInd                  \
+                          [:1], NumMessagesBundled      \
+                          [:1])				\
+  schedule(static)
+        for (v = 0; v < NLVer; v++) {
+            option = -1;
+            // Start: PARALLEL_PROCESS_EXPOSED_VERTEX_B(v)
+            k = candidateMate[v];
+            candidateMate[v] = verLocInd[k];
+            w = candidateMate[v];
+
+#ifdef PRINT_DEBUG_INFO_
+            cout << "\n(" << myRank << ")Processing: " << v + StartIndex << endl;
+            fflush(stdout);
 #endif
+
+#ifdef PRINT_DEBUG_INFO_
+            cout << "\n(" << myRank << ")" << v + StartIndex << " Points to: " << w;
+            fflush(stdout);
+#endif
+            // If found a dominating edge:
+            if (w >= 0) {
+#pragma omp critical(Matching)
+	      {
+		if (isAlreadyMatched(verLocInd[k], StartIndex, EndIndex, GMate, Mate, Ghost2LocalMap)) {
+		  w = computeCandidateMateS(verLocPtr[v], verLocPtr[v + 1], edgeLocWeight, 0,
+					    verLocInd, StartIndex, EndIndex,
+					    GMate,  Mate,  Ghost2LocalMap);
+		  candidateMate[v] = w;
+		}
+		if (w >= 0) {
+		  (*myCard)++;
+		  if ((w < StartIndex) || (w > EndIndex)) { // w is a ghost vertex
+		    option = 2;
+		    if (candidateMate[NLVer + Ghost2LocalMap[w]] == v + StartIndex) {
+		      option = 1;
+		      Mate[v] = w;
+		      GMate[Ghost2LocalMap[w]] = v + StartIndex; // w is a Ghost
+		    } // End of if CandidateMate[w] = v		  
+		  } // End of if a Ghost Vertex
+		  else  { // w is a local vertex
+		    if (candidateMate[w - StartIndex] == (v + StartIndex)) {
+		      option = 3;
+		      Mate[v] = w;                           // v is local
+		      Mate[w - StartIndex] = v + StartIndex; // w is local
+#ifdef PRINT_DEBUG_INFO_
+		      cout << "\n(" << myRank << ")MATCH: (" << v + StartIndex << "," << w << ") ";
+		      fflush(stdout);
+#endif
+		    } // End of if ( candidateMate[w-StartIndex] == (v+StartIndex) )
+		  }     // End of Else
+		} // End of second if
+	      }
+	      
+            } // End of if(w >=0)
+            else  {
+	      //#pragma omp critical(adjuse)
+	      {
+		// This piece of code is executed a really small number of times
+		adj11 = verLocPtr[v];
+		adj12 = verLocPtr[v + 1];
+		for (k1 = adj11; k1 < adj12; k1++) {
+		  w = verLocInd[k1];
+		  if ((w < StartIndex) || (w > EndIndex)) { // A ghost
+		    
+#ifdef PRINT_DEBUG_INFO_
+		    cout << "\n(" << myRank << ")Sending a failure message: ";
+		    cout << "\n(" << myRank << ")Ghost is " << w << " Owner is: " << findOwnerOfGhost(w, verDistance, myRank, numProcs);
+		    fflush(stdout);
+#endif
+		    (*msgInd)++;
+		    (*NumMessagesBundled)++;
+		    ghostOwner = findOwnerOfGhost(w, verDistance, myRank, numProcs);
+		    // assert(ghostOwner != -1);
+		    // assert(ghostOwner != myRank);
+#pragma omp atomic
+		    PCounter[ghostOwner]++;
+		    
+		    privateQLocalVtx.push_back(v + StartIndex);
+		    privateQGhostVtx.push_back(w);
+		    privateQMsgType.push_back(FAILURE);
+		    privateQOwner.push_back(ghostOwner);
+		    
+		  } // End of if(GHOST)
+		}     // End of for loop
+	      }
+            }
+            // End:   PARALLEL_PROCESS_EXPOSED_VERTEX_B(v)
+	    
+            switch (option)
+	      {
+	      case -1:
+                break;
+	      case 1:
+                privateU.push_back(v + StartIndex);
+                privateU.push_back(w);
+		
+#ifdef PRINT_DEBUG_INFO_
+                cout << "\n(" << myRank << ")MATCH: (" << v + StartIndex << "," << w << ")";
+                fflush(stdout);
+#endif
+		
+                //  Decrement the counter:
+                PROCESS_CROSS_EDGE(&Counter[Ghost2LocalMap[w]], S);
+	      case 2:
+#ifdef PRINT_DEBUG_INFO_
+                cout << "\n(" << myRank << ")Sending a request message (291):";
+                cout << "\n(" << myRank << ")Local is: " << v + StartIndex << " Ghost is " << w << " Owner is: " << findOwnerOfGhost(w, verDistance, myRank, numProcs) << endl;
+                fflush(stdout);
+#endif
+                (*msgInd)++;
+                (*NumMessagesBundled)++;
+                ghostOwner = findOwnerOfGhost(w, verDistance, myRank, numProcs);
+                // assert(ghostOwner != -1);
+                // assert(ghostOwner != myRank);
+#pragma omp atomic
+                PCounter[ghostOwner]++;
+		
+                privateQLocalVtx.push_back(v + StartIndex);
+                privateQGhostVtx.push_back(w);
+                privateQMsgType.push_back(REQUEST);
+                privateQOwner.push_back(ghostOwner);
+                break;
+	      case 3:
+	      default:
+                privateU.push_back(v + StartIndex);
+                privateU.push_back(w);
+                break;
+	      }
+	    
+        } // End of for ( v=0; v < NLVer; v++ )
+	
+        queuesTransfer(U, privateU, QLocalVtx,
+                       QGhostVtx,
+                       QMsgType, QOwner, privateQLocalVtx,
+                       privateQGhostVtx,
+                       privateQMsgType,
+                       privateQOwner);
+	
+    } // End of parallel region
+}
+

amgprec/impl/aggregator/processMatchedVertices.cpp

@@ -1,6 +1,5 @@
 #include "MatchBoxPC.h"
-#ifdef OMP
-void processMatchedVertices(
+void processMatchedVerticesD(
     MilanLongInt NLVer,
     vector<MilanLongInt> &UChunkBeingProcessed,
     vector<MilanLongInt> &U,
@@ -91,57 +90,57 @@ void processMatchedVertices(
 		  if (mateVal < 0)  {
 #pragma omp critical
 		    {
-		      if (candidateMate[v - StartIndex] == u) {
+#pragma omp atomic read
+		      mateVal = Mate[v - StartIndex];
+		      // If the current vertex is pointing to a matched vertex and is not matched
+		      if (mateVal < 0)  {
+			
+			if (candidateMate[v - StartIndex] == u) {
 			  // Start: PARALLEL_PROCESS_EXPOSED_VERTEX_B(v)
-			w = computeCandidateMate(verLocPtr[v - StartIndex],
-						 verLocPtr[v - StartIndex + 1],
-						 edgeLocWeight, 0,
-						 verLocInd,
-						 StartIndex,
-						 EndIndex,
-						 GMate,
-						 Mate,
-						 Ghost2LocalMap);
-
-			candidateMate[v - StartIndex] = w;
-
+			  w = computeCandidateMateD(verLocPtr[v - StartIndex],
+						   verLocPtr[v - StartIndex + 1],
+						   edgeLocWeight, 0,
+						   verLocInd, StartIndex, EndIndex,
+						   GMate, Mate, Ghost2LocalMap);
+			  candidateMate[v - StartIndex] = w;
 #ifdef PRINT_DEBUG_INFO_
-			cout << "\n(" << myRank << ")" << v << " Points to: " << w;
-			fflush(stdout);
+			  cout << "\n(" << myRank << ")" << v << " Points to: " << w;
+			  fflush(stdout);
 #endif
-			// If found a dominating edge:
-			if (w >= 0)   {
-			  if ((w < StartIndex) || (w > EndIndex))  { // A ghost
+			  // If found a dominating edge:
+			  if (w >= 0)   {
+			    if ((w < StartIndex) || (w > EndIndex))  { // A ghost
 #ifdef PRINT_DEBUG_INFO_
-			    cout << "\n(" << myRank << ")Sending a request message:";
-			    cout << "\n(" << myRank << ")Ghost is " << w << " Owner is: " << findOwnerOfGhost(w, verDistance, myRank, numProcs);
+			      cout << "\n(" << myRank << ")Sending a request message:";
+			      cout << "\n(" << myRank << ")Ghost is " << w << " Owner is: " << findOwnerOfGhost(w, verDistance, myRank, numProcs);
 #endif
-			    option = 2;
-
-			    if (candidateMate[NLVer + Ghost2LocalMap[w]] == v)  {
-			      option = 1;
-			      Mate[v - StartIndex] = w;     // v is a local vertex
-			      GMate[Ghost2LocalMap[w]] = v; // w is a ghost vertex
-
-			    } // End of if CandidateMate[w] = v
-			  }     // End of if a Ghost Vertex
-			  else    { // w is a local vertex
-			    if (candidateMate[w - StartIndex] == v)  {
-			      option = 3;
-			      Mate[v - StartIndex] = w; // v is a local vertex
-			      Mate[w - StartIndex] = v; // w is a local vertex
-
+			      option = 2;
+			      
+			      if (candidateMate[NLVer + Ghost2LocalMap[w]] == v)  {
+				option = 1;
+				Mate[v - StartIndex] = w;     // v is a local vertex
+				GMate[Ghost2LocalMap[w]] = v; // w is a ghost vertex
+				
+			      } // End of if CandidateMate[w] = v
+			    }     // End of if a Ghost Vertex
+			    else    { // w is a local vertex
+			      if (candidateMate[w - StartIndex] == v)  {
+				option = 3;
+				Mate[v - StartIndex] = w; // v is a local vertex
+				Mate[w - StartIndex] = v; // w is a local vertex
+				
 #ifdef PRINT_DEBUG_INFO_
-			      cout << "\n(" << myRank << ")MATCH: (" << v << "," << w << ") ";
-			      fflush(stdout);
+				cout << "\n(" << myRank << ")MATCH: (" << v << "," << w << ") ";
+				fflush(stdout);
 #endif
-			    } // End of if(CandidateMate(w) = v
-			  }     // End of Else
-			}         // End of if(w >=0)
-			else
-			  option = 4; // End of Else: w == -1
-			// End:   PARALLEL_PROCESS_EXPOSED_VERTEX_B(v)
-		      } // End of If (candidateMate[v-StartIndex] == u
+			      } // End of if(CandidateMate(w) = v
+			    }     // End of Else
+			  }         // End of if(w >=0)
+			  else
+			    option = 4; // End of Else: w == -1
+			  // End:   PARALLEL_PROCESS_EXPOSED_VERTEX_B(v)
+			} // End of If (candidateMate[v-StartIndex] == u
+		      }
 		    }     // End of task
 		  }         // mateval < 0
 		}             // End of if ( (v >= StartIndex) && (v <= EndIndex) ) //If Local Vertex:
@@ -179,6 +178,7 @@ void processMatchedVertices(
 		    ghostOwner = findOwnerOfGhost(w, verDistance, myRank, numProcs);
 		    // assert(ghostOwner != -1);
 		    // assert(ghostOwner != myRank);
+#pragma omp atomic		    
 		    PCounter[ghostOwner]++;
 		    (*NumMessagesBundled)++;
 		    (*msgInd)++;
@@ -211,7 +211,7 @@ void processMatchedVertices(
 			ghostOwner = findOwnerOfGhost(w, verDistance, myRank, numProcs);
 			// assert(ghostOwner != -1);
 			// assert(ghostOwner != myRank);
-
+#pragma omp atomic		    
 			PCounter[ghostOwner]++;
 			(*NumMessagesBundled)++;
 			(*msgInd)++;
@@ -248,7 +248,6 @@ void processMatchedVertices(
 
 		    break;
 		  } // End of switch
-
 	      } // End of inner for
 	    }
 	  } // End of outer for
@@ -265,17 +264,15 @@ void processMatchedVertices(
 	    U.insert(U.end(), privateU.begin(), privateU.end());
 	  }
 
-	  privateU.clear();
-
 #pragma omp critical(sendMessageTransfer)
 	  {
-
 	    QLocalVtx.insert(QLocalVtx.end(), privateQLocalVtx.begin(), privateQLocalVtx.end());
 	    QGhostVtx.insert(QGhostVtx.end(), privateQGhostVtx.begin(), privateQGhostVtx.end());
 	    QMsgType.insert(QMsgType.end(), privateQMsgType.begin(), privateQMsgType.end());
 	    QOwner.insert(QOwner.end(), privateQOwner.begin(), privateQOwner.end());
 	  }
 
+	  privateU.clear();
 	  privateQLocalVtx.clear();
 	  privateQGhostVtx.clear();
 	  privateQMsgType.clear();
@@ -292,4 +289,297 @@ void processMatchedVertices(
 #endif
     } // End of parallel region
 }
+
+
+
+void processMatchedVerticesS(
+    MilanLongInt NLVer,
+    vector<MilanLongInt> &UChunkBeingProcessed,
+    vector<MilanLongInt> &U,
+    vector<MilanLongInt> &privateU,
+    MilanLongInt StartIndex,
+    MilanLongInt EndIndex,
+    MilanLongInt *myCard,
+    MilanLongInt *msgInd,
+    MilanLongInt *NumMessagesBundled,
+    MilanLongInt *SPtr,
+    MilanLongInt *verLocPtr,
+    MilanLongInt *verLocInd,
+    MilanLongInt *verDistance,
+    MilanLongInt *PCounter,
+    vector<MilanLongInt> &Counter,
+    MilanInt myRank,
+    MilanInt numProcs,
+    MilanLongInt *candidateMate,
+    vector<MilanLongInt> &GMate,
+    MilanLongInt *Mate,
+    map<MilanLongInt, MilanLongInt> &Ghost2LocalMap,
+    MilanFloat *edgeLocWeight,
+    vector<MilanLongInt> &QLocalVtx,
+    vector<MilanLongInt> &QGhostVtx,
+    vector<MilanLongInt> &QMsgType,
+    vector<MilanInt> &QOwner,
+    vector<MilanLongInt> &privateQLocalVtx,
+    vector<MilanLongInt> &privateQGhostVtx,
+    vector<MilanLongInt> &privateQMsgType,
+    vector<MilanInt> &privateQOwner)
+{
+
+    MilanLongInt adj1, adj2, adj11, adj12, k, k1, v = -1, w = -1, ghostOwner;
+    int option;
+    MilanLongInt mateVal;
+
+#ifdef PRINT_DEBUG_INFO_
+    cout << "\n(" << myRank << "=========================************===============================" << endl;
+    fflush(stdout);
+    fflush(stdout);
 #endif
+
+#ifdef COUNT_LOCAL_VERTEX
+    MilanLongInt localVertices = 0;
+#endif
+    //#pragma omp parallel private(k, w, v, k1, adj1, adj2, adj11, adj12, ghostOwner, option) \
+    firstprivate(privateU, StartIndex, EndIndex, privateQLocalVtx, privateQGhostVtx, \
+       privateQMsgType, privateQOwner, UChunkBeingProcessed) \
+       default(shared)  num_threads(NUM_THREAD)   \
+            reduction(+                           \
+                      : msgInd[:1], PCounter      \
+                      [:numProcs], myCard         \
+                      [:1], NumMessagesBundled    \
+                      [:1])
+    {
+
+        while (!U.empty()) {
+
+	  extractUChunk(UChunkBeingProcessed, U, privateU);
+
+	  for (MilanLongInt u : UChunkBeingProcessed)   {
+#ifdef PRINT_DEBUG_INFO_
+	    cout << "\n(" << myRank << ")u: " << u;
+	    fflush(stdout);
+#endif
+	    if ((u >= StartIndex) && (u <= EndIndex))     { // Process Only the Local Vertices
+
+#ifdef COUNT_LOCAL_VERTEX
+	      localVertices++;
+#endif
+
+	      // Get the Adjacency list for u
+	      adj1 = verLocPtr[u - StartIndex]; // Pointer
+	      adj2 = verLocPtr[u - StartIndex + 1];
+	      for (k = adj1; k < adj2; k++) {
+		option = -1;
+		v = verLocInd[k];
+
+		if ((v >= StartIndex) && (v <= EndIndex)) { // If Local Vertex:
+
+#ifdef PRINT_DEBUG_INFO_
+		  cout << "\n(" << myRank << ")v: " << v << " c(v)= " << candidateMate[v - StartIndex] << " Mate[v]: " << Mate[v];
+		  fflush(stdout);
+#endif
+#pragma omp atomic read
+		  mateVal = Mate[v - StartIndex];
+		  // If the current vertex is pointing to a matched vertex and is not matched
+		  if (mateVal < 0)  {
+#pragma omp critical
+		    {
+#pragma omp atomic read
+		      mateVal = Mate[v - StartIndex];
+		      // If the current vertex is pointing to a matched vertex and is not matched
+		      if (mateVal < 0)  {
+			
+			if (candidateMate[v - StartIndex] == u) {
+			  // Start: PARALLEL_PROCESS_EXPOSED_VERTEX_B(v)
+			  w = computeCandidateMateS(verLocPtr[v - StartIndex],
+						   verLocPtr[v - StartIndex + 1],
+						   edgeLocWeight, 0,
+						   verLocInd, StartIndex, EndIndex,
+						   GMate, Mate, Ghost2LocalMap);
+			  candidateMate[v - StartIndex] = w;
+#ifdef PRINT_DEBUG_INFO_
+			  cout << "\n(" << myRank << ")" << v << " Points to: " << w;
+			  fflush(stdout);
+#endif
+			  // If found a dominating edge:
+			  if (w >= 0)   {
+			    if ((w < StartIndex) || (w > EndIndex))  { // A ghost
+#ifdef PRINT_DEBUG_INFO_
+			      cout << "\n(" << myRank << ")Sending a request message:";
+			      cout << "\n(" << myRank << ")Ghost is " << w << " Owner is: " << findOwnerOfGhost(w, verDistance, myRank, numProcs);
+#endif
+			      option = 2;
+			      
+			      if (candidateMate[NLVer + Ghost2LocalMap[w]] == v)  {
+				option = 1;
+				Mate[v - StartIndex] = w;     // v is a local vertex
+				GMate[Ghost2LocalMap[w]] = v; // w is a ghost vertex
+				
+			      } // End of if CandidateMate[w] = v
+			    }     // End of if a Ghost Vertex
+			    else    { // w is a local vertex
+			      if (candidateMate[w - StartIndex] == v)  {
+				option = 3;
+				Mate[v - StartIndex] = w; // v is a local vertex
+				Mate[w - StartIndex] = v; // w is a local vertex
+				
+#ifdef PRINT_DEBUG_INFO_
+				cout << "\n(" << myRank << ")MATCH: (" << v << "," << w << ") ";
+				fflush(stdout);
+#endif
+			      } // End of if(CandidateMate(w) = v
+			    }     // End of Else
+			  }         // End of if(w >=0)
+			  else
+			    option = 4; // End of Else: w == -1
+			  // End:   PARALLEL_PROCESS_EXPOSED_VERTEX_B(v)
+			} // End of If (candidateMate[v-StartIndex] == u
+		      }
+		    }     // End of task
+		  }         // mateval < 0
+		}             // End of if ( (v >= StartIndex) && (v <= EndIndex) ) //If Local Vertex:
+		else  { // Neighbor is a ghost vertex
+
+#pragma omp critical
+		  {
+		    if (candidateMate[NLVer + Ghost2LocalMap[v]] == u)
+		      candidateMate[NLVer + Ghost2LocalMap[v]] = -1;
+		    if (v != Mate[u - StartIndex])
+		      option = 5; // u is local
+		  }                   // End of critical
+		}                       // End of Else //A Ghost Vertex
+
+		switch (option)
+		  {
+		  case -1:
+		    // No things to do
+		    break;
+		  case 1:
+		    // Found a dominating edge, it is a ghost and candidateMate[NLVer + Ghost2LocalMap[w]] == v
+		    privateU.push_back(v);
+		    privateU.push_back(w);
+
+		    (*myCard)++;
+#ifdef PRINT_DEBUG_INFO_
+		    cout << "\n(" << myRank << ")MATCH: (" << v << "," << w << ") ";
+		    fflush(stdout);
+#endif
+		    // Decrement the counter:
+		    PROCESS_CROSS_EDGE(&Counter[Ghost2LocalMap[w]], SPtr);
+		  case 2:
+
+		    // Found a dominating edge, it is a ghost
+		    ghostOwner = findOwnerOfGhost(w, verDistance, myRank, numProcs);
+		    // assert(ghostOwner != -1);
+		    // assert(ghostOwner != myRank);
+#pragma omp atomic		    
+		    PCounter[ghostOwner]++;
+		    (*NumMessagesBundled)++;
+		    (*msgInd)++;
+
+		    privateQLocalVtx.push_back(v);
+		    privateQGhostVtx.push_back(w);
+		    privateQMsgType.push_back(REQUEST);
+		    privateQOwner.push_back(ghostOwner);
+		    break;
+		  case 3:
+		    privateU.push_back(v);
+		    privateU.push_back(w);
+
+		    (*myCard)++;
+		    break;
+		  case 4:
+		    // Could not find a dominating vertex
+		    adj11 = verLocPtr[v - StartIndex];
+		    adj12 = verLocPtr[v - StartIndex + 1];
+		    for (k1 = adj11; k1 < adj12; k1++) {
+		      w = verLocInd[k1];
+		      if ((w < StartIndex) || (w > EndIndex)) { // A ghost
+
+#ifdef PRINT_DEBUG_INFO_
+			cout << "\n(" << myRank << ")Sending a failure message: ";
+			cout << "\n(" << myRank << ")Ghost is " << w << " Owner is: " << findOwnerOfGhost(w, verDistance, myRank, numProcs);
+			fflush(stdout);
+#endif
+
+			ghostOwner = findOwnerOfGhost(w, verDistance, myRank, numProcs);
+			// assert(ghostOwner != -1);
+			// assert(ghostOwner != myRank);
+#pragma omp atomic		    
+			PCounter[ghostOwner]++;
+			(*NumMessagesBundled)++;
+			(*msgInd)++;
+
+			privateQLocalVtx.push_back(v);
+			privateQGhostVtx.push_back(w);
+			privateQMsgType.push_back(FAILURE);
+			privateQOwner.push_back(ghostOwner);
+
+		      } // End of if(GHOST)
+		    }     // End of for loop
+		    break;
+		  case 5:
+		  default:
+
+#ifdef PRINT_DEBUG_INFO_
+		    cout << "\n(" << myRank << ")Sending a success message: ";
+		    cout << "\n(" << myRank << ")Ghost is " << v << " Owner is: " << findOwnerOfGhost(v, verDistance, myRank, numProcs) << "\n";
+		    fflush(stdout);
+#endif
+
+		    ghostOwner = findOwnerOfGhost(v, verDistance, myRank, numProcs);
+		    // assert(ghostOwner != -1);
+		    // assert(ghostOwner != myRank);
+
+		    (*NumMessagesBundled)++;
+		    PCounter[ghostOwner]++;
+		    (*msgInd)++;
+
+		    privateQLocalVtx.push_back(u);
+		    privateQGhostVtx.push_back(v);
+		    privateQMsgType.push_back(SUCCESS);
+		    privateQOwner.push_back(ghostOwner);
+
+		    break;
+		  } // End of switch
+	      } // End of inner for
+	    }
+	  } // End of outer for
+
+	  queuesTransfer(U, privateU, QLocalVtx,
+			 QGhostVtx,
+			 QMsgType, QOwner, privateQLocalVtx,
+			 privateQGhostVtx,
+			 privateQMsgType,
+			 privateQOwner);
+
+#pragma omp critical(U)
+	  {
+	    U.insert(U.end(), privateU.begin(), privateU.end());
+	  }
+
+#pragma omp critical(sendMessageTransfer)
+	  {
+	    QLocalVtx.insert(QLocalVtx.end(), privateQLocalVtx.begin(), privateQLocalVtx.end());
+	    QGhostVtx.insert(QGhostVtx.end(), privateQGhostVtx.begin(), privateQGhostVtx.end());
+	    QMsgType.insert(QMsgType.end(), privateQMsgType.begin(), privateQMsgType.end());
+	    QOwner.insert(QOwner.end(), privateQOwner.begin(), privateQOwner.end());
+	  }
+
+	  privateU.clear();
+	  privateQLocalVtx.clear();
+	  privateQGhostVtx.clear();
+	  privateQMsgType.clear();
+	  privateQOwner.clear();
+
+        } // End of while ( !U.empty() )
+
+#ifdef COUNT_LOCAL_VERTEX
+        printf("Count local vertexes: %ld for thread %d of processor %d\n",
+               localVertices,
+               omp_get_thread_num(),
+               myRank);
+
+#endif
+    } // End of parallel region
+}
+

amgprec/impl/aggregator/processMatchedVerticesAndSendMessages.cpp

@@ -1,7 +1,6 @@
 #include "MatchBoxPC.h"
-#ifdef OMP
 //#define DEBUG_HANG_
-void processMatchedVerticesAndSendMessages(
+void processMatchedVerticesAndSendMessagesD(
 					   MilanLongInt NLVer,
 					   vector<MilanLongInt> &UChunkBeingProcessed,
 					   vector<MilanLongInt> &U,
@@ -27,7 +26,7 @@ void processMatchedVerticesAndSendMessages(
 					   vector<MilanLongInt> &QLocalVtx,
 					   vector<MilanLongInt> &QGhostVtx,
 					   vector<MilanLongInt> &QMsgType,
-					   vector<MilanInt> &QOwner,
+ 					   vector<MilanInt> &QOwner,
 					   vector<MilanLongInt> &privateQLocalVtx,
 					   vector<MilanLongInt> &privateQGhostVtx,
 					   vector<MilanLongInt> &privateQMsgType,
@@ -52,15 +51,16 @@ void processMatchedVerticesAndSendMessages(
     MilanLongInt localVertices = 0;
 #endif
     //#pragma omp parallel private(k, w, v, k1, adj1, adj2, adj11, adj12, ghostOwner, option) \
-    firstprivate(Message, privateU, StartIndex, EndIndex, privateQLocalVtx, privateQGhostVtx,\
-		 privateQMsgType, privateQOwner, UChunkBeingProcessed) default(shared) \
-        num_threads(NUM_THREAD)                                                        \
-            reduction(+                                                                \
-                      : msgInd[:1], PCounter                                           \
-                      [:numProcs], myCard                                              \
-                      [:1], NumMessagesBundled                                         \
-                      [:1], msgActual                                                  \
-                      [:1])
+  firstprivate(Message, privateU, StartIndex, EndIndex, privateQLocalVtx, \
+	       privateQGhostVtx, privateQMsgType, privateQOwner, UChunkBeingProcessed) \
+  default(shared)							\
+      num_threads(NUM_THREAD)						\
+  reduction(+								\
+	    : msgInd[:1], PCounter					\
+	    [:numProcs], myCard						\
+	    [:1], NumMessagesBundled					\
+	    [:1], msgActual						\
+	    [:1])
     {
 
         while (!U.empty()) {
@@ -77,7 +77,6 @@ void processMatchedVerticesAndSendMessages(
 #ifdef COUNT_LOCAL_VERTEX
 	      localVertices++;
 #endif
-
 	      // Get the Adjacency list for u
 	      adj1 = verLocPtr[u - StartIndex]; // Pointer
 	      adj2 = verLocPtr[u - StartIndex + 1];
@@ -97,58 +96,57 @@ void processMatchedVerticesAndSendMessages(
 		  if (mateVal < 0) {
 #pragma omp critical
 		    {
-		      if (candidateMate[v - StartIndex] == u)  {
-			// Start: PARALLEL_PROCESS_EXPOSED_VERTEX_B(v)
-			w = computeCandidateMate(verLocPtr[v - StartIndex],
-						 verLocPtr[v - StartIndex + 1],
-						 edgeLocWeight, 0,
-						 verLocInd,
-						 StartIndex,
-						 EndIndex,
-						 GMate,
-						 Mate,
-						 Ghost2LocalMap);
-
-			candidateMate[v - StartIndex] = w;
-
+#pragma omp atomic read
+		      mateVal = Mate[v - StartIndex];
+		      // If the current vertex is pointing to a matched vertex and is not matched
+		      if (mateVal < 0) {
+			
+			if (candidateMate[v - StartIndex] == u)  {
+			  // Start: PARALLEL_PROCESS_EXPOSED_VERTEX_B(v)
+			  w = computeCandidateMateD(verLocPtr[v - StartIndex],
+						   verLocPtr[v - StartIndex + 1],
+						   edgeLocWeight, 0,
+						   verLocInd, StartIndex, EndIndex,
+						   GMate, Mate, Ghost2LocalMap);
+			  candidateMate[v - StartIndex] = w;
 #ifdef PRINT_DEBUG_INFO_
-			cout << "\n(" << myRank << ")" << v << " Points to: " << w;
-			fflush(stdout);
+			  cout << "\n(" << myRank << ")" << v << " Points to: " << w;
+			  fflush(stdout);
 #endif
-			// If found a dominating edge:
-			if (w >= 0)  {
-
-			  if ((w < StartIndex) || (w > EndIndex))  { // A ghost
+			  // If found a dominating edge:
+			  if (w >= 0)  {
+			    
+			    if ((w < StartIndex) || (w > EndIndex))  { // A ghost
 #ifdef PRINT_DEBUG_INFO_
-			    cout << "\n(" << myRank << ")Sending a request message:";
-			    cout << "\n(" << myRank << ")Ghost is " << w << " Owner is: " << findOwnerOfGhost(w, verDistance, myRank, numProcs);
+			      cout << "\n(" << myRank << ")Sending a request message:";
+			      cout << "\n(" << myRank << ")Ghost is " << w << " Owner is: " << findOwnerOfGhost(w, verDistance, myRank, numProcs);
 #endif
-			    option = 2;
-
-			    if (candidateMate[NLVer + Ghost2LocalMap[w]] == v) {
-			      option = 1;
-			      Mate[v - StartIndex] = w;     // v is a local vertex
-			      GMate[Ghost2LocalMap[w]] = v; // w is a ghost vertex
-
-			    } // End of if CandidateMate[w] = v
-			  }     // End of if a Ghost Vertex
-			  else   { // w is a local vertex
-			    if (candidateMate[w - StartIndex] == v)  {
-			      option = 3;
-			      Mate[v - StartIndex] = w; // v is a local vertex
-			      Mate[w - StartIndex] = v; // w is a local vertex
-
+			      option = 2;
+			      
+			      if (candidateMate[NLVer + Ghost2LocalMap[w]] == v) {
+				option = 1;
+				Mate[v - StartIndex] = w;     // v is a local vertex
+				GMate[Ghost2LocalMap[w]] = v; // w is a ghost vertex
+				
+			      } // End of if CandidateMate[w] = v
+			    }     // End of if a Ghost Vertex
+			    else   { // w is a local vertex
+			      if (candidateMate[w - StartIndex] == v)  {
+				option = 3;
+				Mate[v - StartIndex] = w; // v is a local vertex
+				Mate[w - StartIndex] = v; // w is a local vertex
 #ifdef PRINT_DEBUG_INFO_
-			      cout << "\n(" << myRank << ")MATCH: (" << v << "," << w << ") ";
-			      fflush(stdout);
+				cout << "\n(" << myRank << ")MATCH: (" << v << "," << w << ") ";
+				fflush(stdout);
 #endif
-			    } // End of if(CandidateMate(w) = v
-			  }     // End of Else
-			}         // End of if(w >=0)
-			else
-			  option = 4; // End of Else: w == -1
-			// End:   PARALLEL_PROCESS_EXPOSED_VERTEX_B(v)
-		      } // End of If (candidateMate[v-StartIndex] == u
+			      } // End of if(CandidateMate(w) = v
+			    }     // End of Else
+			  }         // End of if(w >=0)
+			  else
+			    option = 4; // End of Else: w == -1
+			  // End:   PARALLEL_PROCESS_EXPOSED_VERTEX_B(v)
+			} // End of If (candidateMate[v-StartIndex] == u
+		      }
 		    }     // End of task
 		  }         // mateval < 0
 		}             // End of if ( (v >= StartIndex) && (v <= EndIndex) ) //If Local Vertex:
@@ -211,15 +209,12 @@ void processMatchedVerticesAndSendMessages(
 		    for (k1 = adj11; k1 < adj12; k1++) {
 		      w = verLocInd[k1];
 		      if ((w < StartIndex) || (w > EndIndex))  { // A ghost
-
 #ifdef PRINT_DEBUG_INFO_
 			cout << "\n(" << myRank << ")Sending a failure message: ";
 			cout << "\n(" << myRank << ")Ghost is " << w << " Owner is: " << findOwnerOfGhost(w, verDistance, myRank, numProcs);
 			fflush(stdout);
 #endif
-
 			ghostOwner = findOwnerOfGhost(w, verDistance, myRank, numProcs);
-
 			// Build the Message Packet:
 			// Message[0] = v;       // LOCAL
 			// Message[1] = w;       // GHOST
@@ -229,7 +224,6 @@ void processMatchedVerticesAndSendMessages(
 
 			(*msgActual)++;
 			(*msgInd)++;
-
 			privateQLocalVtx.push_back(v);
 			privateQGhostVtx.push_back(w);
 			privateQMsgType.push_back(FAILURE);
@@ -248,17 +242,14 @@ void processMatchedVerticesAndSendMessages(
 #endif
 
 		    ghostOwner = findOwnerOfGhost(v, verDistance, myRank, numProcs);
-
 		    // Build the Message Packet:
 		    // Message[0] = u;       // LOCAL
 		    // Message[1] = v;       // GHOST
 		    // Message[2] = SUCCESS; // TYPE
 		    // Send a Request (Asynchronous)
 		    // MPI_Bsend(&Message[0], 3, TypeMap<MilanLongInt>(), ghostOwner, ComputeTag, comm);
-
 		    (*msgActual)++;
 		    (*msgInd)++;
-
 		    privateQLocalVtx.push_back(u);
 		    privateQGhostVtx.push_back(v);
 		    privateQMsgType.push_back(SUCCESS);
@@ -281,10 +272,7 @@ void processMatchedVerticesAndSendMessages(
 
 #ifdef COUNT_LOCAL_VERTEX
         printf("Count local vertexes: %ld for thread %d of processor %d\n",
-               localVertices,
-               omp_get_thread_num(),
-               myRank);
-
+               localVertices, mp_get_thread_num(), myRank);
 #endif
     } // End of parallel region
 
@@ -293,12 +281,10 @@ void processMatchedVerticesAndSendMessages(
     cout << myRank<<" Sending: "<<QOwner.size()-initialSize<<" messages" <<endl;
 #endif
     for (int i = initialSize; i < QOwner.size(); i++)   {
-
       Message[0] = QLocalVtx[i];
       Message[1] = QGhostVtx[i];
       Message[2] = QMsgType[i];
       ghostOwner = QOwner[i];
-
       //MPI_Bsend(&Message[0], 3, TypeMap<MilanLongInt>(), ghostOwner, ComputeTag, comm);
       //cout << myRank<<" Sending to "<<ghostOwner<<endl;
       MPI_Bsend(&Message[0], 3, TypeMap<MilanLongInt>(), ghostOwner, ComputeTag, comm);
@@ -307,4 +293,299 @@ void processMatchedVerticesAndSendMessages(
     cout << myRank<<" Done sending messages"<<endl;
 #endif
 }
+
+
+
+void processMatchedVerticesAndSendMessagesS(
+					   MilanLongInt NLVer,
+					   vector<MilanLongInt> &UChunkBeingProcessed,
+					   vector<MilanLongInt> &U,
+					   vector<MilanLongInt> &privateU,
+					   MilanLongInt StartIndex,
+					   MilanLongInt EndIndex,
+					   MilanLongInt *myCard,
+					   MilanLongInt *msgInd,
+					   MilanLongInt *NumMessagesBundled,
+					   MilanLongInt *SPtr,
+					   MilanLongInt *verLocPtr,
+					   MilanLongInt *verLocInd,
+					   MilanLongInt *verDistance,
+					   MilanLongInt *PCounter,
+					   vector<MilanLongInt> &Counter,
+					   MilanInt myRank,
+					   MilanInt numProcs,
+					   MilanLongInt *candidateMate,
+					   vector<MilanLongInt> &GMate,
+					   MilanLongInt *Mate,
+					   map<MilanLongInt, MilanLongInt> &Ghost2LocalMap,
+					   MilanFloat *edgeLocWeight,
+					   vector<MilanLongInt> &QLocalVtx,
+					   vector<MilanLongInt> &QGhostVtx,
+					   vector<MilanLongInt> &QMsgType,
+					   vector<MilanInt> &QOwner,
+					   vector<MilanLongInt> &privateQLocalVtx,
+					   vector<MilanLongInt> &privateQGhostVtx,
+					   vector<MilanLongInt> &privateQMsgType,
+					   vector<MilanInt> &privateQOwner,
+					   MPI_Comm comm,
+					   MilanLongInt *msgActual,
+					   vector<MilanLongInt> &Message)
+{
+
+    MilanLongInt initialSize = QLocalVtx.size();
+    MilanLongInt adj1, adj2, adj11, adj12, k, k1, v = -1, w = -1, ghostOwner;
+    int option;
+    MilanLongInt mateVal;
+
+#ifdef PRINT_DEBUG_INFO_
+    cout << "\n(" << myRank << "=========================************===============================" << endl;
+    fflush(stdout);
+    fflush(stdout);
 #endif
+
+#ifdef COUNT_LOCAL_VERTEX
+    MilanLongInt localVertices = 0;
+#endif
+    //#pragma omp parallel private(k, w, v, k1, adj1, adj2, adj11, adj12, ghostOwner, option) \
+    firstprivate(Message, privateU, StartIndex, EndIndex, privateQLocalVtx, privateQGhostVtx,\
+		 privateQMsgType, privateQOwner, UChunkBeingProcessed) default(shared) \
+        num_threads(NUM_THREAD)                                                        \
+            reduction(+                                                                \
+                      : msgInd[:1], PCounter                                           \
+                      [:numProcs], myCard                                              \
+                      [:1], NumMessagesBundled                                         \
+                      [:1], msgActual                                                  \
+                      [:1])
+    {
+
+        while (!U.empty()) {
+
+	  extractUChunk(UChunkBeingProcessed, U, privateU);
+
+	  for (MilanLongInt u : UChunkBeingProcessed) {
+#ifdef PRINT_DEBUG_INFO_
+	    cout << "\n(" << myRank << ")u: " << u;
+	    fflush(stdout);
+#endif
+	    if ((u >= StartIndex) && (u <= EndIndex)) { // Process Only the Local Vertices
+
+#ifdef COUNT_LOCAL_VERTEX
+	      localVertices++;
+#endif
+	      // Get the Adjacency list for u
+	      adj1 = verLocPtr[u - StartIndex]; // Pointer
+	      adj2 = verLocPtr[u - StartIndex + 1];
+	      for (k = adj1; k < adj2; k++) {
+		option = -1;
+		v = verLocInd[k];
+
+		if ((v >= StartIndex) && (v <= EndIndex)) { // If Local Vertex:
+
+#ifdef PRINT_DEBUG_INFO_
+		  cout << "\n(" << myRank << ")v: " << v << " c(v)= " << candidateMate[v - StartIndex] << " Mate[v]: " << Mate[v];
+		  fflush(stdout);
+#endif
+#pragma omp atomic read
+		  mateVal = Mate[v - StartIndex];
+		  // If the current vertex is pointing to a matched vertex and is not matched
+		  if (mateVal < 0) {
+#pragma omp critical
+		    {
+#pragma omp atomic read
+		      mateVal = Mate[v - StartIndex];
+		      // If the current vertex is pointing to a matched vertex and is not matched
+		      if (mateVal < 0) {
+			
+			if (candidateMate[v - StartIndex] == u)  {
+			  // Start: PARALLEL_PROCESS_EXPOSED_VERTEX_B(v)
+			  w = computeCandidateMateS(verLocPtr[v - StartIndex],
+						   verLocPtr[v - StartIndex + 1],
+						   edgeLocWeight, 0,
+						   verLocInd, StartIndex, EndIndex,
+						   GMate, Mate, Ghost2LocalMap);
+			  candidateMate[v - StartIndex] = w;
+#ifdef PRINT_DEBUG_INFO_
+			  cout << "\n(" << myRank << ")" << v << " Points to: " << w;
+			  fflush(stdout);
+#endif
+			  // If found a dominating edge:
+			  if (w >= 0)  {
+			    
+			    if ((w < StartIndex) || (w > EndIndex))  { // A ghost
+#ifdef PRINT_DEBUG_INFO_
+			      cout << "\n(" << myRank << ")Sending a request message:";
+			      cout << "\n(" << myRank << ")Ghost is " << w << " Owner is: " << findOwnerOfGhost(w, verDistance, myRank, numProcs);
+#endif
+			      option = 2;
+			      
+			      if (candidateMate[NLVer + Ghost2LocalMap[w]] == v) {
+				option = 1;
+				Mate[v - StartIndex] = w;     // v is a local vertex
+				GMate[Ghost2LocalMap[w]] = v; // w is a ghost vertex
+				
+			      } // End of if CandidateMate[w] = v
+			    }     // End of if a Ghost Vertex
+			    else   { // w is a local vertex
+			      if (candidateMate[w - StartIndex] == v)  {
+				option = 3;
+				Mate[v - StartIndex] = w; // v is a local vertex
+				Mate[w - StartIndex] = v; // w is a local vertex
+#ifdef PRINT_DEBUG_INFO_
+				cout << "\n(" << myRank << ")MATCH: (" << v << "," << w << ") ";
+				fflush(stdout);
+#endif
+			      } // End of if(CandidateMate(w) = v
+			    }     // End of Else
+			  }         // End of if(w >=0)
+			  else
+			    option = 4; // End of Else: w == -1
+			  // End:   PARALLEL_PROCESS_EXPOSED_VERTEX_B(v)
+			} // End of If (candidateMate[v-StartIndex] == u
+		      }
+		    }     // End of task
+		  }         // mateval < 0
+		}             // End of if ( (v >= StartIndex) && (v <= EndIndex) ) //If Local Vertex:
+		else { // Neighbor is a ghost vertex
+
+#pragma omp critical
+		  {
+		    if (candidateMate[NLVer + Ghost2LocalMap[v]] == u)
+		      candidateMate[NLVer + Ghost2LocalMap[v]] = -1;
+		    if (v != Mate[u - StartIndex])
+		      option = 5; // u is local
+		  }                   // End of critical
+		}                       // End of Else //A Ghost Vertex
+
+		switch (option)
+		  {
+		  case -1:
+		    // No things to do
+		    break;
+		  case 1:
+		    // Found a dominating edge, it is a ghost and candidateMate[NLVer + Ghost2LocalMap[w]] == v
+		    privateU.push_back(v);
+		    privateU.push_back(w);
+		    (*myCard)++;
+#ifdef PRINT_DEBUG_INFO_
+		    cout << "\n(" << myRank << ")MATCH: (" << v << "," << w << ") ";
+		    fflush(stdout);
+#endif
+		    // Decrement the counter:
+		    PROCESS_CROSS_EDGE(&Counter[Ghost2LocalMap[w]], SPtr);
+		  case 2:
+
+		    // Found a dominating edge, it is a ghost
+		    ghostOwner = findOwnerOfGhost(w, verDistance, myRank, numProcs);
+
+		    // Build the Message Packet:
+		    // Message[0] = v;       // LOCAL
+		    // Message[1] = w;       // GHOST
+		    // Message[2] = REQUEST; // TYPE
+		    // Send a Request (Asynchronous)
+		    // MPI_Bsend(&Message[0], 3, TypeMap<MilanLongInt>(), ghostOwner, ComputeTag, comm);
+
+		    (*msgActual)++;
+		    (*msgInd)++;
+
+		    privateQLocalVtx.push_back(v);
+		    privateQGhostVtx.push_back(w);
+		    privateQMsgType.push_back(REQUEST);
+		    privateQOwner.push_back(ghostOwner);
+		    break;
+		  case 3:
+		    privateU.push_back(v);
+		    privateU.push_back(w);
+		    (*myCard)++;
+		    break;
+		  case 4:
+		    // Could not find a dominating vertex
+		    adj11 = verLocPtr[v - StartIndex];
+		    adj12 = verLocPtr[v - StartIndex + 1];
+		    for (k1 = adj11; k1 < adj12; k1++) {
+		      w = verLocInd[k1];
+		      if ((w < StartIndex) || (w > EndIndex))  { // A ghost
+#ifdef PRINT_DEBUG_INFO_
+			cout << "\n(" << myRank << ")Sending a failure message: ";
+			cout << "\n(" << myRank << ")Ghost is " << w << " Owner is: " << findOwnerOfGhost(w, verDistance, myRank, numProcs);
+			fflush(stdout);
+#endif
+			ghostOwner = findOwnerOfGhost(w, verDistance, myRank, numProcs);
+			// Build the Message Packet:
+			// Message[0] = v;       // LOCAL
+			// Message[1] = w;       // GHOST
+			// Message[2] = FAILURE; // TYPE
+			// Send a Request (Asynchronous)
+			// MPI_Bsend(&Message[0], 3, TypeMap<MilanLongInt>(), ghostOwner, ComputeTag, comm);
+
+			(*msgActual)++;
+			(*msgInd)++;
+			privateQLocalVtx.push_back(v);
+			privateQGhostVtx.push_back(w);
+			privateQMsgType.push_back(FAILURE);
+			privateQOwner.push_back(ghostOwner);
+
+		      } // End of if(GHOST)
+		    }     // End of for loop
+		    break;
+		  case 5:
+		  default:
+
+#ifdef PRINT_DEBUG_INFO_
+		    cout << "\n(" << myRank << ")Sending a success message: ";
+		    cout << "\n(" << myRank << ")Ghost is " << v << " Owner is: " << findOwnerOfGhost(v, verDistance, myRank, numProcs) << "\n";
+		    fflush(stdout);
+#endif
+
+		    ghostOwner = findOwnerOfGhost(v, verDistance, myRank, numProcs);
+		    // Build the Message Packet:
+		    // Message[0] = u;       // LOCAL
+		    // Message[1] = v;       // GHOST
+		    // Message[2] = SUCCESS; // TYPE
+		    // Send a Request (Asynchronous)
+		    // MPI_Bsend(&Message[0], 3, TypeMap<MilanLongInt>(), ghostOwner, ComputeTag, comm);
+		    (*msgActual)++;
+		    (*msgInd)++;
+		    privateQLocalVtx.push_back(u);
+		    privateQGhostVtx.push_back(v);
+		    privateQMsgType.push_back(SUCCESS);
+		    privateQOwner.push_back(ghostOwner);
+
+		    break;
+		  } // End of switch
+	      }     // End of inner for
+	    }
+	  } // End of outer for
+
+	  queuesTransfer(U, privateU, QLocalVtx,
+			 QGhostVtx,
+			 QMsgType, QOwner, privateQLocalVtx,
+			 privateQGhostVtx,
+			 privateQMsgType,
+			 privateQOwner);
+
+        } // End of while ( !U.empty() )
+
+#ifdef COUNT_LOCAL_VERTEX
+        printf("Count local vertexes: %ld for thread %d of processor %d\n",
+               localVertices, mp_get_thread_num(), myRank);
+#endif
+    } // End of parallel region
+
+    // Send the messages
+#ifdef DEBUG_HANG_
+    cout << myRank<<" Sending: "<<QOwner.size()-initialSize<<" messages" <<endl;
+#endif
+    for (int i = initialSize; i < QOwner.size(); i++)   {
+      Message[0] = QLocalVtx[i];
+      Message[1] = QGhostVtx[i];
+      Message[2] = QMsgType[i];
+      ghostOwner = QOwner[i];
+      //MPI_Bsend(&Message[0], 3, TypeMap<MilanLongInt>(), ghostOwner, ComputeTag, comm);
+      //cout << myRank<<" Sending to "<<ghostOwner<<endl;
+      MPI_Bsend(&Message[0], 3, TypeMap<MilanLongInt>(), ghostOwner, ComputeTag, comm);
+    }
+#ifdef DEBUG_HANG_
+    cout << myRank<<" Done sending messages"<<endl;
+#endif
+}
+

amgprec/impl/aggregator/processMessages.cpp

@@ -1,8 +1,8 @@
 #include "MatchBoxPC.h"
-#ifdef OMP
 //#define DEBUG_HANG_
+#if !defined(SERIAL_MPI)
 
-void processMessages(
+void processMessagesD(
     MilanLongInt NLVer,
     MilanLongInt *Mate,
     MilanLongInt *candidateMate,
@@ -212,8 +212,324 @@ void processMessages(
 	    // Process only if not already matched ( v is local)
 	    if (candidateMate[v - StartIndex] == u)  {
 	      // Start: PARALLEL_PROCESS_EXPOSED_VERTEX_B(v)
-	      w = computeCandidateMate(verLocPtr[v - StartIndex], verLocPtr[v - StartIndex + 1], edgeLocWeight, k,
-				       verLocInd, StartIndex, EndIndex, GMate, Mate, Ghost2LocalMap);
+	      w = computeCandidateMateD(verLocPtr[v - StartIndex], verLocPtr[v - StartIndex + 1],
+					edgeLocWeight, k,verLocInd, StartIndex, EndIndex,
+					GMate, Mate, Ghost2LocalMap);
+	      candidateMate[v - StartIndex] = w;
+#ifdef PRINT_DEBUG_INFO_
+	      cout << "\n(" << myRank << ")" << v << " Points to: " << w << endl;
+	      fflush(stdout);
+#endif
+	      // If found a dominating edge:
+	      if (w >= 0)  {
+		if ((w < StartIndex) || (w > EndIndex))   {
+		  // w is a ghost
+		  // Build the Message Packet:
+		  Message[0] = v;       // LOCAL
+		  Message[1] = w;       // GHOST
+		  Message[2] = REQUEST; // TYPE
+		  // Send a Request (Asynchronous)
+#ifdef PRINT_DEBUG_INFO_
+		  cout << "\n(" << myRank << ")Sending a request message: ";
+		  cout << "\n(" << myRank << ")Ghost is " << w << " Owner is: " << findOwnerOfGhost(w, verDistance, myRank, numProcs) << endl;
+		  fflush(stdout);
+#endif
+		  ghostOwner = findOwnerOfGhost(w, verDistance, myRank, numProcs);
+		  //assert(ghostOwner != -1);
+		  //assert(ghostOwner != myRank);
+		  //cout << myRank<<" Sending to "<<ghostOwner<<endl;
+		  MPI_Bsend(&Message[0], 3, TypeMap<MilanLongInt>(), ghostOwner, ComputeTag, comm);
+		  (*msgInd)++;
+		  (*msgActual)++;
+		  if (candidateMate[NLVer + Ghost2LocalMap[w]] == v) {
+		    Mate[v - StartIndex] = w;     // v is local
+		    GMate[Ghost2LocalMap[w]] = v; // w is ghost
+		    U.push_back(v);
+		    U.push_back(w);
+		    (*myCard)++;
+#ifdef PRINT_DEBUG_INFO_
+		    cout << "\n(" << myRank << ")MATCH: (" << v << "," << w << ") " << endl;
+		    fflush(stdout);
+#endif
+
+		    PROCESS_CROSS_EDGE(&Counter[Ghost2LocalMap[w]], S);
+		  } // End of if CandidateMate[w] = v
+		}     // End of if a Ghost Vertex
+		else   { // w is a local vertex
+		  if (candidateMate[w - StartIndex] == v) {
+		    Mate[v - StartIndex] = w; // v is local
+		    Mate[w - StartIndex] = v; // w is local
+		    // Q.push_back(u);
+		    U.push_back(v);
+		    U.push_back(w);
+		    (*myCard)++;
+#ifdef PRINT_DEBUG_INFO_
+		    cout << "\n(" << myRank << ")MATCH: (" << v << "," << w << ") " << endl;
+		    fflush(stdout);
+#endif
+		  } // End of if(CandidateMate(w) = v
+		}     // End of Else
+	      }         // End of if(w >=0)
+	      else  { // No dominant edge found
+		adj11 = verLocPtr[v - StartIndex];
+		adj12 = verLocPtr[v - StartIndex + 1];
+		for (k1 = adj11; k1 < adj12; k1++) {
+		  w = verLocInd[k1];
+		  if ((w < StartIndex) || (w > EndIndex))  {
+		    // A ghost
+		    // Build the Message Packet:
+		    Message[0] = v;       // LOCAL
+		    Message[1] = w;       // GHOST
+		    Message[2] = FAILURE; // TYPE
+		    // Send a Request (Asynchronous)
+#ifdef PRINT_DEBUG_INFO_
+		    cout << "\n(" << myRank << ")Sending a failure message: ";
+		    cout << "\n(" << myRank << ")Ghost is " << w << " Owner is: " << findOwnerOfGhost(w, verDistance, myRank, numProcs) << endl;
+		    fflush(stdout);
+#endif
+		    ghostOwner = findOwnerOfGhost(w, verDistance, myRank, numProcs);
+		    //assert(ghostOwner != -1);
+		    //assert(ghostOwner != myRank);
+		    //cout << myRank<<" Sending to "<<ghostOwner<<endl;
+		    MPI_Bsend(&Message[0], 3, TypeMap<MilanLongInt>(), ghostOwner, ComputeTag, comm);
+		    (*msgInd)++;
+		    (*msgActual)++;
+		  } // End of if(GHOST)
+		}     // End of for loop
+	      }         // End of Else: w == -1
+                        // End:   PARALLEL_PROCESS_EXPOSED_VERTEX_B(v)
+	    } // End of if ( candidateMate[v-StartIndex] == u )
+	  }     // End of if ( Mate[v] == -1 )
+	}         // End of if ( message_type == SUCCESS )
+	else  {
+	  // CASE III: FAILURE
+#ifdef PRINT_DEBUG_INFO_
+	  cout << "\n(" << myRank << ")Message type is FAILURE" << endl;
+	  fflush(stdout);
+#endif
+	  GMate[Ghost2LocalMap[u]] = EndIndex + 1;            // Set a Dummy Mate to make sure that we do not (u is a ghost) process this anymore
+	  PROCESS_CROSS_EDGE(&Counter[Ghost2LocalMap[u]], S); // Decrease the counter
+	}                                                       // End of else: CASE III
+      }                                                           // End of else: CASE I
+    }
+
+    return;
+}
+
+
+void processMessagesS(
+    MilanLongInt NLVer,
+    MilanLongInt *Mate,
+    MilanLongInt *candidateMate,
+    map<MilanLongInt, MilanLongInt> &Ghost2LocalMap,
+    vector<MilanLongInt> &GMate,
+    vector<MilanLongInt> &Counter,
+    MilanLongInt StartIndex,
+    MilanLongInt EndIndex,
+    MilanLongInt *myCard,
+    MilanLongInt *msgInd,
+    MilanLongInt *msgActual,
+    MilanFloat *edgeLocWeight,
+    MilanLongInt *verDistance,
+    MilanLongInt *verLocPtr,
+    MilanLongInt k,
+    MilanLongInt *verLocInd,
+    MilanInt numProcs,
+    MilanInt myRank,
+    MPI_Comm comm,
+    vector<MilanLongInt> &Message,
+    MilanLongInt numGhostEdges,
+    MilanLongInt u,
+    MilanLongInt v,
+    MilanLongInt *S,
+    vector<MilanLongInt> &U)
+{
+
+    //#define PRINT_DEBUG_INFO_
+
+    MilanInt Sender;
+    MPI_Status computeStatus;
+    MilanLongInt bundleSize, w;
+    MilanLongInt adj11, adj12, k1;
+    MilanLongInt ghostOwner;
+    int error_codeC;
+    error_codeC = MPI_Comm_set_errhandler(MPI_COMM_WORLD, MPI_ERRORS_RETURN);
+    char error_message[MPI_MAX_ERROR_STRING];
+    int message_length;
+    MilanLongInt message_type = 0;
+
+    // Buffer to receive bundled messages
+    // Maximum messages that can be received from any processor is
+    // twice the edge cut: REQUEST; REQUEST+(FAILURE/SUCCESS)
+    vector<MilanLongInt> ReceiveBuffer;
+    try
+    {
+        ReceiveBuffer.reserve(numGhostEdges * 2 * 3); // Three integers per cross edge
+    }
+    catch (length_error)
+    {
+        cout << "Error in function algoDistEdgeApproxDominatingEdgesMessageBundling: \n";
+        cout << "Not enough memory to allocate the internal variables \n";
+        exit(1);
+    }
+
+#ifdef PRINT_DEBUG_INFO_
+    cout
+        << "\n(" << myRank << "=========================************===============================" << endl;
+    fflush(stdout);
+    fflush(stdout);
+#endif
+#ifdef PRINT_DEBUG_INFO_
+    cout << "\n(" << myRank << ")About to begin Message processing phase ... *S=" << *S << endl;
+    fflush(stdout);
+#endif
+#ifdef PRINT_DEBUG_INFO_
+    cout << "\n(" << myRank << "=========================************===============================" << endl;
+    fflush(stdout);
+    fflush(stdout);
+#endif
+    // BLOCKING RECEIVE:
+#ifdef PRINT_DEBUG_INFO_
+    cout << "\n(" << myRank << " Waiting for blocking receive..." << endl;
+    fflush(stdout);
+    fflush(stdout);
+#endif
+
+    //cout << myRank<<" Receiving ...";
+    error_codeC = MPI_Recv(&Message[0], 3, TypeMap<MilanLongInt>(), MPI_ANY_SOURCE, ComputeTag, comm, &computeStatus);
+    if (error_codeC != MPI_SUCCESS)
+    {
+        MPI_Error_string(error_codeC, error_message, &message_length);
+        cout << "\n*Error in call to MPI_Receive on Slave: " << error_message << "\n";
+        fflush(stdout);
+    }
+    Sender = computeStatus.MPI_SOURCE;
+    //cout << " ...from "<<Sender << endl;
+#ifdef PRINT_DEBUG_INFO_
+    cout << "\n(" << myRank << ")Received message from Process " << Sender << " Type= " << Message[2] << endl;
+    fflush(stdout);
+#endif
+
+    if (Message[2] == SIZEINFO)  {
+
+#ifdef PRINT_DEBUG_INFO_
+      cout << "\n(" << myRank << ")Received bundled message from Process " << Sender << " Size= " << Message[0] << endl;
+      fflush(stdout);
+#endif
+      bundleSize = Message[0]; //#of integers in the message
+      // Build the Message Buffer:
+      if (!ReceiveBuffer.empty())
+	ReceiveBuffer.clear();            // Empty it out first
+      ReceiveBuffer.resize(bundleSize, -1); // Initialize
+#ifdef PRINT_DEBUG_INFO_
+      cout << "\n(" << myRank << ")Message Bundle Before: " << endl;
+      for (int i = 0; i < bundleSize; i++)
+	cout << ReceiveBuffer[i] << ",";
+      cout << endl;
+      fflush(stdout);
+#endif
+      // Receive the message
+      //cout << myRank<<" Receiving from "<<Sender<<endl;
+      error_codeC = MPI_Recv(&ReceiveBuffer[0], bundleSize, TypeMap<MilanLongInt>(), Sender, BundleTag, comm, &computeStatus);
+      if (error_codeC != MPI_SUCCESS) {
+	MPI_Error_string(error_codeC, error_message, &message_length);
+	cout << "\n*Error in call to MPI_Receive on processor " << myRank << " Error: " << error_message << "\n";
+	fflush(stdout);
+      }
+#ifdef PRINT_DEBUG_INFO_
+      cout << "\n(" << myRank << ")Message Bundle After: " << endl;
+      for (int i = 0; i < bundleSize; i++)
+	cout << ReceiveBuffer[i] << ",";
+      cout << endl;
+      fflush(stdout);
+#endif
+    } else   { // Just a single message:
+#ifdef PRINT_DEBUG_INFO_
+      cout << "\n(" << myRank << ")Received regular message from Process " << Sender << " u= " << Message[0] << " v= " << Message[1] << endl;
+      fflush(stdout);
+#endif
+      // Add the current message to Queue:
+      bundleSize = 3; //#of integers in the message
+      // Build the Message Buffer:
+      if (!ReceiveBuffer.empty())
+	ReceiveBuffer.clear();            // Empty it out first
+      ReceiveBuffer.resize(bundleSize, -1); // Initialize
+
+      ReceiveBuffer[0] = Message[0]; // u
+      ReceiveBuffer[1] = Message[1]; // v
+      ReceiveBuffer[2] = Message[2]; // message_type
+    }
+
+#ifdef DEBUG_GHOST_
+    if ((v < StartIndex) || (v > EndIndex)) {
+      cout << "\n(" << myRank << ") From ReceiveBuffer: This should not happen: u= " << u << " v= " << v << " Type= " << message_type << " StartIndex " << StartIndex << " EndIndex " << EndIndex << endl;
+      fflush(stdout);
+    }
+#endif
+#ifdef PRINT_DEBUG_INFO_
+    cout << "\n(" << myRank << ")Processing message: u= " << u << " v= " << v << " Type= " << message_type << endl;
+    fflush(stdout);
+#endif
+
+    // Most of the time bundleSize == 3, thus, it's not worth parallelizing thi loop
+    for (MilanLongInt bundleCounter = 3; bundleCounter < bundleSize + 3; bundleCounter += 3)  {
+      u = ReceiveBuffer[bundleCounter - 3];            // GHOST
+      v = ReceiveBuffer[bundleCounter - 2];            // LOCAL
+      message_type = ReceiveBuffer[bundleCounter - 1]; // TYPE
+
+      // CASE I: REQUEST
+      if (message_type == REQUEST)  {
+#ifdef PRINT_DEBUG_INFO_
+	cout << "\n(" << myRank << ")Message type is REQUEST" << endl;
+	fflush(stdout);
+#endif
+#ifdef DEBUG_GHOST_
+	if ((v < 0) || (v < StartIndex) || ((v - StartIndex) > NLVer)) {
+	  cout << "\n(" << myRank << ") case 1 Bad address " << v << " " << StartIndex << " " << v - StartIndex << " " << NLVer << endl;
+	  fflush(stdout);
+	}
+
+#endif
+
+	if (Mate[v - StartIndex] == -1) {
+	  // Process only if not already matched  (v is local)
+	  candidateMate[NLVer + Ghost2LocalMap[u]] = v; // Set CandidateMate for the ghost
+	  if (candidateMate[v - StartIndex] == u)   {
+	    GMate[Ghost2LocalMap[u]] = v; // u is ghost
+	    Mate[v - StartIndex] = u;     // v is local
+	    U.push_back(v);
+	    U.push_back(u);
+	    (*myCard)++;
+#ifdef PRINT_DEBUG_INFO_
+	    cout << "\n(" << myRank << ")MATCH: (" << v << "," << u << ") " << endl;
+	    fflush(stdout);
+#endif
+
+	    PROCESS_CROSS_EDGE(&Counter[Ghost2LocalMap[u]], S);
+	  } // End of if ( candidateMate[v-StartIndex] == u )e
+	}     // End of if ( Mate[v] == -1 )
+      }         // End of REQUEST
+      else  { // CASE II: SUCCESS
+	if (message_type == SUCCESS) {
+#ifdef PRINT_DEBUG_INFO_
+	  cout << "\n(" << myRank << ")Message type is SUCCESS" << endl;
+	  fflush(stdout);
+#endif
+	  GMate[Ghost2LocalMap[u]] = EndIndex + 1; // Set a Dummy Mate to make sure that we do not (u is a ghost) process it again
+	  PROCESS_CROSS_EDGE(&Counter[Ghost2LocalMap[u]], S);
+#ifdef DEBUG_GHOST_
+	  if ((v < 0) || (v < StartIndex) || ((v - StartIndex) > NLVer)) {
+	    cout << "\n(" << myRank << ") case 2  Bad address " << v << " " << StartIndex << " " << v - StartIndex << " " << NLVer << endl;
+	    fflush(stdout);
+	  }
+#endif
+	  if (Mate[v - StartIndex] == -1) {
+	    // Process only if not already matched ( v is local)
+	    if (candidateMate[v - StartIndex] == u)  {
+	      // Start: PARALLEL_PROCESS_EXPOSED_VERTEX_B(v)
+	      w = computeCandidateMateS(verLocPtr[v - StartIndex], verLocPtr[v - StartIndex + 1],
+					edgeLocWeight, k,verLocInd, StartIndex, EndIndex,
+					GMate, Mate, Ghost2LocalMap);
 	      candidateMate[v - StartIndex] = w;
 #ifdef PRINT_DEBUG_INFO_
 	      cout << "\n(" << myRank << ")" << v << " Points to: " << w << endl;

amgprec/impl/aggregator/queueTransfer.cpp

@@ -1,5 +1,4 @@
 #include "MatchBoxPC.h"
-#ifdef OMP
 void queuesTransfer(vector<MilanLongInt> &U,
                     vector<MilanLongInt> &privateU,
                     vector<MilanLongInt> &QLocalVtx,
@@ -15,23 +14,20 @@ void queuesTransfer(vector<MilanLongInt> &U,
 #pragma omp critical(U)
     {
         U.insert(U.end(), privateU.begin(), privateU.end());
-    }
-
-    privateU.clear();
-
-#pragma omp critical(sendMessageTransfer)
-    {
 
+    //#pragma omp critical(sendMessageTransfer)
+    //{
         QLocalVtx.insert(QLocalVtx.end(), privateQLocalVtx.begin(), privateQLocalVtx.end());
         QGhostVtx.insert(QGhostVtx.end(), privateQGhostVtx.begin(), privateQGhostVtx.end());
         QMsgType.insert(QMsgType.end(), privateQMsgType.begin(), privateQMsgType.end());
         QOwner.insert(QOwner.end(), privateQOwner.begin(), privateQOwner.end());
     }
 
+    privateU.clear();
     privateQLocalVtx.clear();
     privateQGhostVtx.clear();
     privateQMsgType.clear();
     privateQOwner.clear();
 
 }
-#endif
+

amgprec/impl/aggregator/sendBundledMessages.cpp

@@ -1,24 +1,23 @@
 #include "MatchBoxPC.h"
-#ifdef OMP
 void sendBundledMessages(MilanLongInt *numGhostEdges,
-                                MilanInt *BufferSize,
-                                MilanLongInt *Buffer,
-                                vector<MilanLongInt> &PCumulative,
-                                vector<MilanLongInt> &PMessageBundle,
-                                vector<MilanLongInt> &PSizeInfoMessages,
-                                MilanLongInt *PCounter,
-                                MilanLongInt NumMessagesBundled,
-                                MilanLongInt *msgActual,
-                                MilanLongInt *msgInd,
-                                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)
+			 MilanInt *BufferSize,
+			 MilanLongInt *Buffer,
+			 vector<MilanLongInt> &PCumulative,
+			 vector<MilanLongInt> &PMessageBundle,
+			 vector<MilanLongInt> &PSizeInfoMessages,
+			 MilanLongInt *PCounter,
+			 MilanLongInt NumMessagesBundled,
+			 MilanLongInt *msgActual,
+			 MilanLongInt *msgInd,
+			 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)
 {
 
     MilanLongInt myIndex = 0, numMessagesToSend;
@@ -207,4 +206,3 @@ void sendBundledMessages(MilanLongInt *numGhostEdges,
 }
 }
 }
-#endif

amgprec/impl/amg_c_hierarchy_bld.F90

@@ -451,11 +451,16 @@ subroutine amg_c_hierarchy_bld(a,desc_a,prec,info)
     if (.not.associated(prec%precv(1)%base_desc,desc_a)) then
       prec%precv(1)%base_desc => prec%precv(1)%desc_ac
     end if
-    do i=2, iszv 
+    do i=2, iszv
       prec%precv(i)%base_a       => prec%precv(i)%ac
       prec%precv(i)%base_desc    => prec%precv(i)%desc_ac
-      prec%precv(i)%linmap%p_desc_U => prec%precv(i-1)%base_desc
-      prec%precv(i)%linmap%p_desc_V => prec%precv(i)%base_desc
+      ! This is needed when the linmap object has been built
+      ! reusing the base_desc descriptor through a pointer.
+      ! With PSBLAS 4 we will have a better solution
+      if (associated(prec%precv(i)%linmap%p_desc_U)) &
+           & prec%precv(i)%linmap%p_desc_U => prec%precv(i-1)%base_desc
+      if (associated(prec%precv(i)%linmap%p_desc_V))&
+           & prec%precv(i)%linmap%p_desc_V => prec%precv(i)%base_desc
     end do
   end if
 
@@ -507,6 +512,71 @@ subroutine amg_c_hierarchy_bld(a,desc_a,prec,info)
   return
 
 contains
+
+#if ( __GNUC__ == 13 && __GNUC_MINOR__ == 3) 
+  ! gfortran 13.3.0 generates a strange error here with MOLD
+  ! moving to SOURCE but only for this version, since it's heavier
+  subroutine save_smoothers(level,save1, save2,info)
+    type(amg_c_onelev_type), intent(inout) :: level
+    class(amg_c_base_smoother_type), allocatable , intent(inout) :: save1, save2
+    integer(psb_ipk_), intent(out) :: info
+
+    info  = 0 
+    if (allocated(save1)) then
+      call save1%free(info)
+      if (info  == 0) deallocate(save1,stat=info)
+      if (info /= 0) return
+    end if
+    if (allocated(save2)) then
+      call save2%free(info)
+      if (info  == 0) deallocate(save2,stat=info)
+      if (info /= 0) return
+    end if
+    allocate(save1, source=level%sm,stat=info)
+    if (info == 0) call level%sm%clone_settings(save1,info)
+    if ((info == 0).and.allocated(level%sm2a)) then
+      allocate(save2, source=level%sm2a,stat=info) 
+      if (info == 0) call level%sm2a%clone_settings(save2,info)
+    end if
+
+    return
+  end subroutine save_smoothers
+
+  subroutine restore_smoothers(level,save1, save2,info)
+    type(amg_c_onelev_type), intent(inout), target :: level
+    class(amg_c_base_smoother_type), allocatable, intent(inout) :: save1, save2
+    integer(psb_ipk_), intent(out) :: info
+
+    info  = 0
+
+    if (allocated(level%sm)) then
+      if (info  == 0) call level%sm%free(info)
+      if (info  == 0) deallocate(level%sm,stat=info)
+    end if
+    if (allocated(save1)) then
+      if (info  == 0) allocate(level%sm,source=save1,stat=info)
+      if (info == 0) call save1%clone_settings(level%sm,info)
+    end if
+
+    if (info /= 0) return
+
+    if (allocated(level%sm2a)) then
+      if (info  == 0) call level%sm2a%free(info)
+      if (info  == 0) deallocate(level%sm2a,stat=info)
+    end if
+    if (allocated(save2)) then
+      if (info  == 0) allocate(level%sm2a,source=save2,stat=info)
+      if (info == 0) call save2%clone_settings(level%sm2a,info)
+      if (info  == 0) level%sm2 => level%sm2a
+    else
+      if (allocated(level%sm)) level%sm2 => level%sm
+    end if
+
+    return
+  end subroutine restore_smoothers
+
+#else
+
   subroutine save_smoothers(level,save1, save2,info)
     type(amg_c_onelev_type), intent(inout) :: level
     class(amg_c_base_smoother_type), allocatable , intent(inout) :: save1, save2
@@ -565,5 +635,5 @@ contains
 
     return
   end subroutine restore_smoothers
-
+#endif
 end subroutine amg_c_hierarchy_bld

amgprec/impl/amg_d_hierarchy_bld.F90

@@ -451,11 +451,16 @@ subroutine amg_d_hierarchy_bld(a,desc_a,prec,info)
     if (.not.associated(prec%precv(1)%base_desc,desc_a)) then
       prec%precv(1)%base_desc => prec%precv(1)%desc_ac
     end if
-    do i=2, iszv 
+    do i=2, iszv
       prec%precv(i)%base_a       => prec%precv(i)%ac
       prec%precv(i)%base_desc    => prec%precv(i)%desc_ac
-      prec%precv(i)%linmap%p_desc_U => prec%precv(i-1)%base_desc
-      prec%precv(i)%linmap%p_desc_V => prec%precv(i)%base_desc
+      ! This is needed when the linmap object has been built
+      ! reusing the base_desc descriptor through a pointer.
+      ! With PSBLAS 4 we will have a better solution
+      if (associated(prec%precv(i)%linmap%p_desc_U)) &
+           & prec%precv(i)%linmap%p_desc_U => prec%precv(i-1)%base_desc
+      if (associated(prec%precv(i)%linmap%p_desc_V))&
+           & prec%precv(i)%linmap%p_desc_V => prec%precv(i)%base_desc
     end do
   end if
 
@@ -507,6 +512,71 @@ subroutine amg_d_hierarchy_bld(a,desc_a,prec,info)
   return
 
 contains
+
+#if ( __GNUC__ == 13 && __GNUC_MINOR__ == 3) 
+  ! gfortran 13.3.0 generates a strange error here with MOLD
+  ! moving to SOURCE but only for this version, since it's heavier
+  subroutine save_smoothers(level,save1, save2,info)
+    type(amg_d_onelev_type), intent(inout) :: level
+    class(amg_d_base_smoother_type), allocatable , intent(inout) :: save1, save2
+    integer(psb_ipk_), intent(out) :: info
+
+    info  = 0 
+    if (allocated(save1)) then
+      call save1%free(info)
+      if (info  == 0) deallocate(save1,stat=info)
+      if (info /= 0) return
+    end if
+    if (allocated(save2)) then
+      call save2%free(info)
+      if (info  == 0) deallocate(save2,stat=info)
+      if (info /= 0) return
+    end if
+    allocate(save1, source=level%sm,stat=info)
+    if (info == 0) call level%sm%clone_settings(save1,info)
+    if ((info == 0).and.allocated(level%sm2a)) then
+      allocate(save2, source=level%sm2a,stat=info) 
+      if (info == 0) call level%sm2a%clone_settings(save2,info)
+    end if
+
+    return
+  end subroutine save_smoothers
+
+  subroutine restore_smoothers(level,save1, save2,info)
+    type(amg_d_onelev_type), intent(inout), target :: level
+    class(amg_d_base_smoother_type), allocatable, intent(inout) :: save1, save2
+    integer(psb_ipk_), intent(out) :: info
+
+    info  = 0
+
+    if (allocated(level%sm)) then
+      if (info  == 0) call level%sm%free(info)
+      if (info  == 0) deallocate(level%sm,stat=info)
+    end if
+    if (allocated(save1)) then
+      if (info  == 0) allocate(level%sm,source=save1,stat=info)
+      if (info == 0) call save1%clone_settings(level%sm,info)
+    end if
+
+    if (info /= 0) return
+
+    if (allocated(level%sm2a)) then
+      if (info  == 0) call level%sm2a%free(info)
+      if (info  == 0) deallocate(level%sm2a,stat=info)
+    end if
+    if (allocated(save2)) then
+      if (info  == 0) allocate(level%sm2a,source=save2,stat=info)
+      if (info == 0) call save2%clone_settings(level%sm2a,info)
+      if (info  == 0) level%sm2 => level%sm2a
+    else
+      if (allocated(level%sm)) level%sm2 => level%sm
+    end if
+
+    return
+  end subroutine restore_smoothers
+
+#else
+
   subroutine save_smoothers(level,save1, save2,info)
     type(amg_d_onelev_type), intent(inout) :: level
     class(amg_d_base_smoother_type), allocatable , intent(inout) :: save1, save2
@@ -565,5 +635,5 @@ contains
 
     return
   end subroutine restore_smoothers
-
+#endif
 end subroutine amg_d_hierarchy_bld

amgprec/impl/amg_s_hierarchy_bld.F90

@@ -451,11 +451,16 @@ subroutine amg_s_hierarchy_bld(a,desc_a,prec,info)
     if (.not.associated(prec%precv(1)%base_desc,desc_a)) then
       prec%precv(1)%base_desc => prec%precv(1)%desc_ac
     end if
-    do i=2, iszv 
+    do i=2, iszv
       prec%precv(i)%base_a       => prec%precv(i)%ac
       prec%precv(i)%base_desc    => prec%precv(i)%desc_ac
-      prec%precv(i)%linmap%p_desc_U => prec%precv(i-1)%base_desc
-      prec%precv(i)%linmap%p_desc_V => prec%precv(i)%base_desc
+      ! This is needed when the linmap object has been built
+      ! reusing the base_desc descriptor through a pointer.
+      ! With PSBLAS 4 we will have a better solution
+      if (associated(prec%precv(i)%linmap%p_desc_U)) &
+           & prec%precv(i)%linmap%p_desc_U => prec%precv(i-1)%base_desc
+      if (associated(prec%precv(i)%linmap%p_desc_V))&
+           & prec%precv(i)%linmap%p_desc_V => prec%precv(i)%base_desc
     end do
   end if
 
@@ -507,6 +512,71 @@ subroutine amg_s_hierarchy_bld(a,desc_a,prec,info)
   return
 
 contains
+
+#if ( __GNUC__ == 13 && __GNUC_MINOR__ == 3) 
+  ! gfortran 13.3.0 generates a strange error here with MOLD
+  ! moving to SOURCE but only for this version, since it's heavier
+  subroutine save_smoothers(level,save1, save2,info)
+    type(amg_s_onelev_type), intent(inout) :: level
+    class(amg_s_base_smoother_type), allocatable , intent(inout) :: save1, save2
+    integer(psb_ipk_), intent(out) :: info
+
+    info  = 0 
+    if (allocated(save1)) then
+      call save1%free(info)
+      if (info  == 0) deallocate(save1,stat=info)
+      if (info /= 0) return
+    end if
+    if (allocated(save2)) then
+      call save2%free(info)
+      if (info  == 0) deallocate(save2,stat=info)
+      if (info /= 0) return
+    end if
+    allocate(save1, source=level%sm,stat=info)
+    if (info == 0) call level%sm%clone_settings(save1,info)
+    if ((info == 0).and.allocated(level%sm2a)) then
+      allocate(save2, source=level%sm2a,stat=info) 
+      if (info == 0) call level%sm2a%clone_settings(save2,info)
+    end if
+
+    return
+  end subroutine save_smoothers
+
+  subroutine restore_smoothers(level,save1, save2,info)
+    type(amg_s_onelev_type), intent(inout), target :: level
+    class(amg_s_base_smoother_type), allocatable, intent(inout) :: save1, save2
+    integer(psb_ipk_), intent(out) :: info
+
+    info  = 0
+
+    if (allocated(level%sm)) then
+      if (info  == 0) call level%sm%free(info)
+      if (info  == 0) deallocate(level%sm,stat=info)
+    end if
+    if (allocated(save1)) then
+      if (info  == 0) allocate(level%sm,source=save1,stat=info)
+      if (info == 0) call save1%clone_settings(level%sm,info)
+    end if
+
+    if (info /= 0) return
+
+    if (allocated(level%sm2a)) then
+      if (info  == 0) call level%sm2a%free(info)
+      if (info  == 0) deallocate(level%sm2a,stat=info)
+    end if
+    if (allocated(save2)) then
+      if (info  == 0) allocate(level%sm2a,source=save2,stat=info)
+      if (info == 0) call save2%clone_settings(level%sm2a,info)
+      if (info  == 0) level%sm2 => level%sm2a
+    else
+      if (allocated(level%sm)) level%sm2 => level%sm
+    end if
+
+    return
+  end subroutine restore_smoothers
+
+#else
+
   subroutine save_smoothers(level,save1, save2,info)
     type(amg_s_onelev_type), intent(inout) :: level
     class(amg_s_base_smoother_type), allocatable , intent(inout) :: save1, save2
@@ -565,5 +635,5 @@ contains
 
     return
   end subroutine restore_smoothers
-
+#endif
 end subroutine amg_s_hierarchy_bld

amgprec/impl/amg_z_hierarchy_bld.F90

@@ -451,11 +451,16 @@ subroutine amg_z_hierarchy_bld(a,desc_a,prec,info)
     if (.not.associated(prec%precv(1)%base_desc,desc_a)) then
       prec%precv(1)%base_desc => prec%precv(1)%desc_ac
     end if
-    do i=2, iszv 
+    do i=2, iszv
       prec%precv(i)%base_a       => prec%precv(i)%ac
       prec%precv(i)%base_desc    => prec%precv(i)%desc_ac
-      prec%precv(i)%linmap%p_desc_U => prec%precv(i-1)%base_desc
-      prec%precv(i)%linmap%p_desc_V => prec%precv(i)%base_desc
+      ! This is needed when the linmap object has been built
+      ! reusing the base_desc descriptor through a pointer.
+      ! With PSBLAS 4 we will have a better solution
+      if (associated(prec%precv(i)%linmap%p_desc_U)) &
+           & prec%precv(i)%linmap%p_desc_U => prec%precv(i-1)%base_desc
+      if (associated(prec%precv(i)%linmap%p_desc_V))&
+           & prec%precv(i)%linmap%p_desc_V => prec%precv(i)%base_desc
     end do
   end if
 
@@ -507,6 +512,71 @@ subroutine amg_z_hierarchy_bld(a,desc_a,prec,info)
   return
 
 contains
+
+#if ( __GNUC__ == 13 && __GNUC_MINOR__ == 3) 
+  ! gfortran 13.3.0 generates a strange error here with MOLD
+  ! moving to SOURCE but only for this version, since it's heavier
+  subroutine save_smoothers(level,save1, save2,info)
+    type(amg_z_onelev_type), intent(inout) :: level
+    class(amg_z_base_smoother_type), allocatable , intent(inout) :: save1, save2
+    integer(psb_ipk_), intent(out) :: info
+
+    info  = 0 
+    if (allocated(save1)) then
+      call save1%free(info)
+      if (info  == 0) deallocate(save1,stat=info)
+      if (info /= 0) return
+    end if
+    if (allocated(save2)) then
+      call save2%free(info)
+      if (info  == 0) deallocate(save2,stat=info)
+      if (info /= 0) return
+    end if
+    allocate(save1, source=level%sm,stat=info)
+    if (info == 0) call level%sm%clone_settings(save1,info)
+    if ((info == 0).and.allocated(level%sm2a)) then
+      allocate(save2, source=level%sm2a,stat=info) 
+      if (info == 0) call level%sm2a%clone_settings(save2,info)
+    end if
+
+    return
+  end subroutine save_smoothers
+
+  subroutine restore_smoothers(level,save1, save2,info)
+    type(amg_z_onelev_type), intent(inout), target :: level
+    class(amg_z_base_smoother_type), allocatable, intent(inout) :: save1, save2
+    integer(psb_ipk_), intent(out) :: info
+
+    info  = 0
+
+    if (allocated(level%sm)) then
+      if (info  == 0) call level%sm%free(info)
+      if (info  == 0) deallocate(level%sm,stat=info)
+    end if
+    if (allocated(save1)) then
+      if (info  == 0) allocate(level%sm,source=save1,stat=info)
+      if (info == 0) call save1%clone_settings(level%sm,info)
+    end if
+
+    if (info /= 0) return
+
+    if (allocated(level%sm2a)) then
+      if (info  == 0) call level%sm2a%free(info)
+      if (info  == 0) deallocate(level%sm2a,stat=info)
+    end if
+    if (allocated(save2)) then
+      if (info  == 0) allocate(level%sm2a,source=save2,stat=info)
+      if (info == 0) call save2%clone_settings(level%sm2a,info)
+      if (info  == 0) level%sm2 => level%sm2a
+    else
+      if (allocated(level%sm)) level%sm2 => level%sm
+    end if
+
+    return
+  end subroutine restore_smoothers
+
+#else
+
   subroutine save_smoothers(level,save1, save2,info)
     type(amg_z_onelev_type), intent(inout) :: level
     class(amg_z_base_smoother_type), allocatable , intent(inout) :: save1, save2
@@ -565,5 +635,5 @@ contains
 
     return
   end subroutine restore_smoothers
-
+#endif
 end subroutine amg_z_hierarchy_bld

amgprec/impl/level/amg_c_base_onelev_csetc.F90

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

amgprec/impl/level/amg_c_base_onelev_dump.f90

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

amgprec/impl/level/amg_c_base_onelev_memory_use.f90

@@ -89,7 +89,7 @@ subroutine amg_c_base_onelev_memory_use(lv,il,nl,ilmin,info,iout,verbosity,prefi
   if (present(global)) then
     global_ = global
   else
-    global_ = .false.
+    global_ = .true.
   end if
 
   if (present(prefix)) then
@@ -98,8 +98,7 @@ subroutine amg_c_base_onelev_memory_use(lv,il,nl,ilmin,info,iout,verbosity,prefi
     prefix_ = ""
   end if
 
-  write(iout_,*) trim(prefix_)
-
+  if ((me == 0).or.(verbosity_>0)) write(iout_,*) trim(prefix_)
 
   if (global_) then
     allocate(sz(6))

amgprec/impl/level/amg_d_base_onelev_csetc.F90

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

amgprec/impl/level/amg_d_base_onelev_dump.f90

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

amgprec/impl/level/amg_d_base_onelev_memory_use.f90

@@ -89,7 +89,7 @@ subroutine amg_d_base_onelev_memory_use(lv,il,nl,ilmin,info,iout,verbosity,prefi
   if (present(global)) then
     global_ = global
   else
-    global_ = .false.
+    global_ = .true.
   end if
 
   if (present(prefix)) then
@@ -98,8 +98,7 @@ subroutine amg_d_base_onelev_memory_use(lv,il,nl,ilmin,info,iout,verbosity,prefi
     prefix_ = ""
   end if
 
-  write(iout_,*) trim(prefix_)
-
+  if ((me == 0).or.(verbosity_>0)) write(iout_,*) trim(prefix_)
 
   if (global_) then
     allocate(sz(6))

amgprec/impl/level/amg_s_base_onelev_csetc.F90

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

amgprec/impl/level/amg_s_base_onelev_dump.f90

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

amgprec/impl/level/amg_s_base_onelev_memory_use.f90

@@ -89,7 +89,7 @@ subroutine amg_s_base_onelev_memory_use(lv,il,nl,ilmin,info,iout,verbosity,prefi
   if (present(global)) then
     global_ = global
   else
-    global_ = .false.
+    global_ = .true.
   end if
 
   if (present(prefix)) then
@@ -98,8 +98,7 @@ subroutine amg_s_base_onelev_memory_use(lv,il,nl,ilmin,info,iout,verbosity,prefi
     prefix_ = ""
   end if
 
-  write(iout_,*) trim(prefix_)
-
+  if ((me == 0).or.(verbosity_>0)) write(iout_,*) trim(prefix_)
 
   if (global_) then
     allocate(sz(6))

amgprec/impl/level/amg_z_base_onelev_csetc.F90

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

amgprec/impl/level/amg_z_base_onelev_dump.f90

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

amgprec/impl/level/amg_z_base_onelev_memory_use.f90

@@ -89,7 +89,7 @@ subroutine amg_z_base_onelev_memory_use(lv,il,nl,ilmin,info,iout,verbosity,prefi
   if (present(global)) then
     global_ = global
   else
-    global_ = .false.
+    global_ = .true.
   end if
 
   if (present(prefix)) then
@@ -98,8 +98,7 @@ subroutine amg_z_base_onelev_memory_use(lv,il,nl,ilmin,info,iout,verbosity,prefi
     prefix_ = ""
   end if
 
-  write(iout_,*) trim(prefix_)
-
+  if ((me == 0).or.(verbosity_>0)) write(iout_,*) trim(prefix_)
 
   if (global_) then
     allocate(sz(6))

amgprec/impl/smoother/amg_c_jac_smoother_apply_vect.f90

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

amgprec/impl/smoother/amg_d_jac_smoother_apply_vect.f90

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

amgprec/impl/smoother/amg_d_poly_smoother_apply_vect.f90

@@ -36,11 +36,11 @@
 !
 !
 subroutine amg_d_poly_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
-     & sweeps,work,wv,info,init,initu)
+     & sweeps,work,wv,info,init,initu) 
 
   use psb_base_mod
   use amg_d_diag_solver
-  use psb_base_krylov_conv_mod, only : log_conv
+  use psb_base_linsolve_conv_mod, only : log_conv
   use amg_d_poly_smoother, amg_protect_name => amg_d_poly_smoother_apply_vect
   implicit none
   type(psb_desc_type), intent(in)                 :: desc_data
@@ -55,6 +55,10 @@ subroutine amg_d_poly_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
   integer(psb_ipk_), intent(out)                  :: info
   character, intent(in), optional                :: init
   type(psb_d_vect_type),intent(inout), optional   :: initu
+  ! Timers
+  logical, parameter  :: do_timings=.true.
+  integer(psb_ipk_), save  :: poly_1=-1, poly_2=-1, poly_3=-1
+  integer(psb_ipk_), save  :: poly_mv=-1, poly_sv=-1, poly_vect=-1
   !
   integer(psb_ipk_)    :: n_row,n_col
   type(psb_d_vect_type)  :: tx, ty, tz, r
@@ -92,7 +96,19 @@ subroutine amg_d_poly_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
     call psb_errpush(info,name)
     goto 9999
   end if
-
+  
+  if ((do_timings).and.(poly_1==-1))       &
+       & poly_1 = psb_get_timer_idx("POLY: Chebychev4")
+  if ((do_timings).and.(poly_2==-1))       &
+       & poly_2 = psb_get_timer_idx("POLY: OptChebychev4")
+  if ((do_timings).and.(poly_3==-1))       &
+       & poly_3 = psb_get_timer_idx("POLY: OptChebychev1")
+  if ((do_timings).and.(poly_mv==-1))       &
+       & poly_mv = psb_get_timer_idx("POLY: spMV")
+  if ((do_timings).and.(poly_vect==-1))       &
+       & poly_vect = psb_get_timer_idx("POLY: Vectors")
+  if ((do_timings).and.(poly_sv==-1))       &
+       & poly_sv = psb_get_timer_idx("POLY: solver")
   n_row = desc_data%get_local_rows()
   n_col = desc_data%get_local_cols()
 
@@ -124,39 +140,40 @@ subroutine amg_d_poly_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
     call tz%zero()
 
     select case(sm%variant)
-    case(amg_poly_lottes_)
+    case(amg_cheb_4_)
+      if (do_timings) call psb_tic(poly_1)
       block
         real(psb_dpk_)      :: cz, cr
         !  b == x
         !  x == tx
         !
-        do i=1, sm%pdegree
+        do i=1, sm%pdegree-1
           !   B r_{k-1}
-          call sm%sv%apply(done,r,dzero,ty,desc_data,trans_,aux,wv(5:),info,init='Z')
+          if (do_timings) call psb_tic(poly_sv)
+          call sm%sv%apply(done,r,dzero,ty,desc_data,trans_,aux,wv(5:),info,init='Z') ! ty = M^{-1} r
+          if (do_timings) call psb_toc(poly_sv)
           cz = (2*i*done-3)/(2*i*done+done)
           cr = (8*i*done-4)/((2*i*done+done)*sm%rho_ba)
-          if (.false.) then
-            ! z_k =  cz z_{k-1} + cr ty = cz z_{k-1} + cr Br_{k-1}
-            call psb_geaxpby(cr,ty,cz,tz,desc_data,info)
-            ! r_k =    b-Ax_k  = x -A tx
-            call psb_geaxpby(done,tz,done,tx,desc_data,info)
-          else
-            call psb_abgdxyz(cr,cz,done,done,ty,tz,tx,desc_data,info)
-          end if
-          if (.false.) then
-            call psb_geaxpby(done,x,dzero,r,desc_data,info)
-            call psb_spmm(-done,sm%pa,tx,done,r,desc_data,info,work=aux,trans=trans_)
-          else
-            call psb_spmm(-done,sm%pa,tz,done,r,desc_data,info,work=aux,trans=trans_)
-          end if
-!!$        res  = psb_genrm2(r,desc_data,info)
-!!$        write(0,*) 'Polynomial smoother  LOTTES',i,res
-          ! x_k =  x_{k-1} + z_k
+          if (do_timings) call psb_tic(poly_vect)
+          call psb_upd_xyz(cr,cz,done,done,ty,tz,tx,desc_data,info) ! zk = cz * zk-1 + cr * rk-1 
+          if (do_timings) call psb_toc(poly_vect)
+          if (do_timings) call psb_tic(poly_mv)
+          call psb_spmm(-done,sm%pa,tz,done,r,desc_data,info,work=aux,trans=trans_)
+          if (do_timings) call psb_toc(poly_mv)
         end do
+        if (do_timings) call psb_tic(poly_sv)
+        call sm%sv%apply(done,r,dzero,ty,desc_data,trans_,aux,wv(5:),info,init='Z') ! ty = M^{-1} r
+        if (do_timings) call psb_toc(poly_sv)
+        cz = (2*sm%pdegree*done-3)/(2*sm%pdegree*done+done)
+        cr = (8*sm%pdegree*done-4)/((2*sm%pdegree*done+done)*sm%rho_ba)
+        if (do_timings) call psb_tic(poly_vect)
+        call psb_upd_xyz(cr,cz,done,done,ty,tz,tx,desc_data,info)
+        if (do_timings) call psb_toc(poly_vect)
       end block
+      if (do_timings) call psb_toc(poly_1)
 
-    case(amg_poly_lottes_beta_)
-
+    case(amg_cheb_4_opt_)
+      if (do_timings) call psb_tic(poly_2)
       block
         real(psb_dpk_)      :: cz, cr
         !  b == x
@@ -170,32 +187,30 @@ subroutine amg_d_poly_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
           sm%poly_beta(1:sm%pdegree) = amg_d_poly_beta_mat(1:sm%pdegree,sm%pdegree)
         end if
 
-        do i=1, sm%pdegree
+        do i=1, sm%pdegree-1
           !   B r_{k-1}
+          if (do_timings) call psb_tic(poly_sv)
           call sm%sv%apply(done,r,dzero,ty,desc_data,trans_,aux,wv(5:),info,init='Z')
+          if (do_timings) call psb_toc(poly_sv)
           cz = (2*i*done-3)/(2*i*done+done)
           cr = (8*i*done-4)/((2*i*done+done)*sm%rho_ba)
-          if (.false.) then
-            ! z_k =  cz z_{k-1} + cr ty = cz z_{k-1} + cr Br_{k-1}
-            call psb_geaxpby(cr,ty,cz,tz,desc_data,info)
-            ! r_k =    b-Ax_k  = x -A tx
-            call psb_geaxpby(sm%poly_beta(i),tz,done,tx,desc_data,info)
-          else
-            call psb_abgdxyz(cr,cz,sm%poly_beta(i),done,ty,tz,tx,desc_data,info)
-          end if
-          if (.false.) then
-            call psb_geaxpby(done,x,dzero,r,desc_data,info)
-            call psb_spmm(-done,sm%pa,tx,done,r,desc_data,info,work=aux,trans=trans_)
-          else
-            call psb_spmm(-done,sm%pa,tz,done,r,desc_data,info,work=aux,trans=trans_)
-          end if
-!!$        res  = psb_genrm2(r,desc_data,info)
-!!$        write(0,*) 'Polynomial smoother LOTTES_BETA',i,res
-          ! x_k =  x_{k-1} + z_k
+          if (do_timings) call psb_tic(poly_vect)
+          call psb_upd_xyz(cr,cz,sm%poly_beta(i),done,ty,tz,tx,desc_data,info)
+          if (do_timings) call psb_toc(poly_vect)
+          if (do_timings) call psb_tic(poly_mv)
+          call psb_spmm(-done,sm%pa,tz,done,r,desc_data,info,work=aux,trans=trans_)
+          if (do_timings) call psb_toc(poly_mv)
         end do
+        call sm%sv%apply(done,r,dzero,ty,desc_data,trans_,aux,wv(5:),info,init='Z')
+        cz = (2*sm%pdegree*done-3)/(2*sm%pdegree*done+done)
+        cr = (8*sm%pdegree*done-4)/((2*sm%pdegree*done+done)*sm%rho_ba)
+        if (do_timings) call psb_tic(poly_vect)
+        call psb_upd_xyz(cr,cz,sm%poly_beta(sm%pdegree),done,ty,tz,tx,desc_data,info)
+        if (do_timings) call psb_toc(poly_vect)
       end block
-
-    case(amg_poly_new_)
+      if (do_timings) call psb_toc(poly_2)
+    case(amg_cheb_1_opt_)
+      if (do_timings) call psb_tic(poly_3)
       block
         real(psb_dpk_)      :: sigma, theta, delta, rho_old, rho
         !  b == x
@@ -206,40 +221,35 @@ subroutine amg_d_poly_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
         delta = (done-sm%cf_a)/2
         sigma = theta/delta
         rho_old = done/sigma
+        if (do_timings) call psb_tic(poly_sv)
         call sm%sv%apply(done,r,dzero,ty,desc_data,trans_,aux,wv(5:),info,init='Z')
+        if (do_timings) call psb_toc(poly_sv)
         call psb_geaxpby((done/sm%rho_ba),ty,dzero,r,desc_data,info)
-        if (.false.) then
-          call psb_geaxpby((done/theta),r,dzero,tz,desc_data,info)
-          call psb_geaxpby(done,tz,done,tx,desc_data,info)
-        else
-          call psb_abgdxyz((done/theta),dzero,done,done,r,tz,tx,desc_data,info)
-        end if
+        if (do_timings) call psb_tic(poly_vect)
+        call psb_upd_xyz((done/theta),dzero,done,done,r,tz,tx,desc_data,info)
+        if (do_timings) call psb_toc(poly_vect)
 
         ! tz == d
         do i=1, sm%pdegree-1
           !
-           !
+          !
           ! r_{k-1} = r_k - (1/rho(BA)) B A d_k
+          if (do_timings) call psb_tic(poly_mv)
           call psb_spmm(done,sm%pa,tz,dzero,ty,desc_data,info,work=aux,trans=trans_)
+          if (do_timings) call psb_toc(poly_mv)
+          if (do_timings) call psb_tic(poly_sv)
           call sm%sv%apply(-(done/sm%rho_ba),ty,done,r,desc_data,trans_,aux,wv(5:),info,init='Z')
-
+          if (do_timings) call psb_toc(poly_sv)
           !
           ! d_{k+1} = (rho rho_old) d_k + 2(rho/delta) r_{k+1}
           rho = done/(2*sigma - rho_old)
-          if (.false.) then
-            call psb_geaxpby((2*rho/delta),r,(rho*rho_old),tz,desc_data,info)
-            call psb_geaxpby(done,tz,done,tx,desc_data,info)
-          else
-            call psb_abgdxyz((2*rho/delta),(rho*rho_old),done,done,r,tz,tx,desc_data,info)
-          end if
-!!$          res  = psb_genrm2(r,desc_data,info)
-!!$          write(0,*) 'Polynomial smoother NEW ',i,res
-          ! x_k =  x_{k-1} + z_k
+          if (do_timings) call psb_tic(poly_vect)
+          call psb_upd_xyz((2*rho/delta),(rho*rho_old),done,done,r,tz,tx,desc_data,info)
+          if (do_timings) call psb_toc(poly_vect)
           rho_old = rho
         end do
       end block
-
-
+      if (do_timings) call psb_toc(poly_3)
     case default
       info=psb_err_internal_error_
       call psb_errpush(info,name,&

amgprec/impl/smoother/amg_d_poly_smoother_bld.f90

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

amgprec/impl/smoother/amg_d_poly_smoother_cseti.f90

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

amgprec/impl/smoother/amg_d_poly_smoother_descr.f90

@@ -77,17 +77,17 @@ subroutine amg_d_poly_smoother_descr(sm,info,iout,coarse,prefix)
 
   write(iout_,*) trim(prefix_), '  Polynomial smoother   '
   select case(sm%variant)
-  case(amg_poly_lottes_)
-    write(iout_,*) trim(prefix_), '        variant:        ','POLY_LOTTES'
+  case(amg_cheb_4_)
+    write(iout_,*) trim(prefix_), '        variant:        ','CHEB_4'
     write(iout_,*) trim(prefix_), '         Degree:        ',sm%pdegree
     write(iout_,*) trim(prefix_), '         rho_ba:        ',sm%rho_ba
-  case(amg_poly_lottes_beta_)
-    write(iout_,*) trim(prefix_), '        variant:        ','POLY_LOTTES_BETA'
+  case(amg_cheb_4_opt_)
+    write(iout_,*) trim(prefix_), '        variant:        ','CHEB_4_OPT'
     write(iout_,*) trim(prefix_), '         Degree:        ',sm%pdegree
     write(iout_,*) trim(prefix_), '         rho_ba:        ',sm%rho_ba
     if (allocated(sm%poly_beta)) write(iout_,*) trim(prefix_), '   Coefficients:        ',sm%poly_beta(1:sm%pdegree)
-  case(amg_poly_new_)
-    write(iout_,*) trim(prefix_), '        variant:        ','POLY_NEW'
+  case(amg_cheb_1_opt_)
+    write(iout_,*) trim(prefix_), '        variant:        ','CHEB_1_OPT'
     write(iout_,*) trim(prefix_), '         Degree:        ',sm%pdegree
     write(iout_,*) trim(prefix_), '         rho_ba:        ',sm%rho_ba
     write(iout_,*) trim(prefix_), '    Coefficient:        ',sm%cf_a

amgprec/impl/smoother/amg_s_jac_smoother_apply_vect.f90

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

amgprec/impl/smoother/amg_s_poly_smoother_apply_vect.f90

@@ -36,11 +36,11 @@
 !
 !
 subroutine amg_s_poly_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
-     & sweeps,work,wv,info,init,initu)
+     & sweeps,work,wv,info,init,initu) 
 
   use psb_base_mod
   use amg_s_diag_solver
-  use psb_base_krylov_conv_mod, only : log_conv
+  use psb_base_linsolve_conv_mod, only : log_conv
   use amg_s_poly_smoother, amg_protect_name => amg_s_poly_smoother_apply_vect
   implicit none
   type(psb_desc_type), intent(in)                 :: desc_data
@@ -55,6 +55,10 @@ subroutine amg_s_poly_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
   integer(psb_ipk_), intent(out)                  :: info
   character, intent(in), optional                :: init
   type(psb_s_vect_type),intent(inout), optional   :: initu
+  ! Timers
+  logical, parameter  :: do_timings=.true.
+  integer(psb_ipk_), save  :: poly_1=-1, poly_2=-1, poly_3=-1
+  integer(psb_ipk_), save  :: poly_mv=-1, poly_sv=-1, poly_vect=-1
   !
   integer(psb_ipk_)    :: n_row,n_col
   type(psb_s_vect_type)  :: tx, ty, tz, r
@@ -92,7 +96,19 @@ subroutine amg_s_poly_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
     call psb_errpush(info,name)
     goto 9999
   end if
-
+  
+  if ((do_timings).and.(poly_1==-1))       &
+       & poly_1 = psb_get_timer_idx("POLY: Chebychev4")
+  if ((do_timings).and.(poly_2==-1))       &
+       & poly_2 = psb_get_timer_idx("POLY: OptChebychev4")
+  if ((do_timings).and.(poly_3==-1))       &
+       & poly_3 = psb_get_timer_idx("POLY: OptChebychev1")
+  if ((do_timings).and.(poly_mv==-1))       &
+       & poly_mv = psb_get_timer_idx("POLY: spMV")
+  if ((do_timings).and.(poly_vect==-1))       &
+       & poly_vect = psb_get_timer_idx("POLY: Vectors")
+  if ((do_timings).and.(poly_sv==-1))       &
+       & poly_sv = psb_get_timer_idx("POLY: solver")
   n_row = desc_data%get_local_rows()
   n_col = desc_data%get_local_cols()
 
@@ -124,39 +140,40 @@ subroutine amg_s_poly_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
     call tz%zero()
 
     select case(sm%variant)
-    case(amg_poly_lottes_)
+    case(amg_cheb_4_)
+      if (do_timings) call psb_tic(poly_1)
       block
         real(psb_spk_)      :: cz, cr
         !  b == x
         !  x == tx
         !
-        do i=1, sm%pdegree
+        do i=1, sm%pdegree-1
           !   B r_{k-1}
-          call sm%sv%apply(sone,r,szero,ty,desc_data,trans_,aux,wv(5:),info,init='Z')
+          if (do_timings) call psb_tic(poly_sv)
+          call sm%sv%apply(sone,r,szero,ty,desc_data,trans_,aux,wv(5:),info,init='Z') ! ty = M^{-1} r
+          if (do_timings) call psb_toc(poly_sv)
           cz = (2*i*sone-3)/(2*i*sone+sone)
           cr = (8*i*sone-4)/((2*i*sone+sone)*sm%rho_ba)
-          if (.false.) then
-            ! z_k =  cz z_{k-1} + cr ty = cz z_{k-1} + cr Br_{k-1}
-            call psb_geaxpby(cr,ty,cz,tz,desc_data,info)
-            ! r_k =    b-Ax_k  = x -A tx
-            call psb_geaxpby(sone,tz,sone,tx,desc_data,info)
-          else
-            call psb_abgdxyz(cr,cz,sone,sone,ty,tz,tx,desc_data,info)
-          end if
-          if (.false.) then
-            call psb_geaxpby(sone,x,szero,r,desc_data,info)
-            call psb_spmm(-sone,sm%pa,tx,sone,r,desc_data,info,work=aux,trans=trans_)
-          else
-            call psb_spmm(-sone,sm%pa,tz,sone,r,desc_data,info,work=aux,trans=trans_)
-          end if
-!!$        res  = psb_genrm2(r,desc_data,info)
-!!$        write(0,*) 'Polynomial smoother  LOTTES',i,res
-          ! x_k =  x_{k-1} + z_k
+          if (do_timings) call psb_tic(poly_vect)
+          call psb_upd_xyz(cr,cz,sone,sone,ty,tz,tx,desc_data,info) ! zk = cz * zk-1 + cr * rk-1 
+          if (do_timings) call psb_toc(poly_vect)
+          if (do_timings) call psb_tic(poly_mv)
+          call psb_spmm(-sone,sm%pa,tz,sone,r,desc_data,info,work=aux,trans=trans_)
+          if (do_timings) call psb_toc(poly_mv)
         end do
+        if (do_timings) call psb_tic(poly_sv)
+        call sm%sv%apply(sone,r,szero,ty,desc_data,trans_,aux,wv(5:),info,init='Z') ! ty = M^{-1} r
+        if (do_timings) call psb_toc(poly_sv)
+        cz = (2*sm%pdegree*sone-3)/(2*sm%pdegree*sone+sone)
+        cr = (8*sm%pdegree*sone-4)/((2*sm%pdegree*sone+sone)*sm%rho_ba)
+        if (do_timings) call psb_tic(poly_vect)
+        call psb_upd_xyz(cr,cz,sone,sone,ty,tz,tx,desc_data,info)
+        if (do_timings) call psb_toc(poly_vect)
       end block
+      if (do_timings) call psb_toc(poly_1)
 
-    case(amg_poly_lottes_beta_)
-
+    case(amg_cheb_4_opt_)
+      if (do_timings) call psb_tic(poly_2)
       block
         real(psb_spk_)      :: cz, cr
         !  b == x
@@ -170,32 +187,30 @@ subroutine amg_s_poly_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
           sm%poly_beta(1:sm%pdegree) = amg_d_poly_beta_mat(1:sm%pdegree,sm%pdegree)
         end if
 
-        do i=1, sm%pdegree
+        do i=1, sm%pdegree-1
           !   B r_{k-1}
+          if (do_timings) call psb_tic(poly_sv)
           call sm%sv%apply(sone,r,szero,ty,desc_data,trans_,aux,wv(5:),info,init='Z')
+          if (do_timings) call psb_toc(poly_sv)
           cz = (2*i*sone-3)/(2*i*sone+sone)
           cr = (8*i*sone-4)/((2*i*sone+sone)*sm%rho_ba)
-          if (.false.) then
-            ! z_k =  cz z_{k-1} + cr ty = cz z_{k-1} + cr Br_{k-1}
-            call psb_geaxpby(cr,ty,cz,tz,desc_data,info)
-            ! r_k =    b-Ax_k  = x -A tx
-            call psb_geaxpby(sm%poly_beta(i),tz,sone,tx,desc_data,info)
-          else
-            call psb_abgdxyz(cr,cz,sm%poly_beta(i),sone,ty,tz,tx,desc_data,info)
-          end if
-          if (.false.) then
-            call psb_geaxpby(sone,x,szero,r,desc_data,info)
-            call psb_spmm(-sone,sm%pa,tx,sone,r,desc_data,info,work=aux,trans=trans_)
-          else
-            call psb_spmm(-sone,sm%pa,tz,sone,r,desc_data,info,work=aux,trans=trans_)
-          end if
-!!$        res  = psb_genrm2(r,desc_data,info)
-!!$        write(0,*) 'Polynomial smoother LOTTES_BETA',i,res
-          ! x_k =  x_{k-1} + z_k
+          if (do_timings) call psb_tic(poly_vect)
+          call psb_upd_xyz(cr,cz,sm%poly_beta(i),sone,ty,tz,tx,desc_data,info)
+          if (do_timings) call psb_toc(poly_vect)
+          if (do_timings) call psb_tic(poly_mv)
+          call psb_spmm(-sone,sm%pa,tz,sone,r,desc_data,info,work=aux,trans=trans_)
+          if (do_timings) call psb_toc(poly_mv)
         end do
+        call sm%sv%apply(sone,r,szero,ty,desc_data,trans_,aux,wv(5:),info,init='Z')
+        cz = (2*sm%pdegree*sone-3)/(2*sm%pdegree*sone+sone)
+        cr = (8*sm%pdegree*sone-4)/((2*sm%pdegree*sone+sone)*sm%rho_ba)
+        if (do_timings) call psb_tic(poly_vect)
+        call psb_upd_xyz(cr,cz,sm%poly_beta(sm%pdegree),sone,ty,tz,tx,desc_data,info)
+        if (do_timings) call psb_toc(poly_vect)
       end block
-
-    case(amg_poly_new_)
+      if (do_timings) call psb_toc(poly_2)
+    case(amg_cheb_1_opt_)
+      if (do_timings) call psb_tic(poly_3)
       block
         real(psb_spk_)      :: sigma, theta, delta, rho_old, rho
         !  b == x
@@ -206,40 +221,35 @@ subroutine amg_s_poly_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
         delta = (sone-sm%cf_a)/2
         sigma = theta/delta
         rho_old = sone/sigma
+        if (do_timings) call psb_tic(poly_sv)
         call sm%sv%apply(sone,r,szero,ty,desc_data,trans_,aux,wv(5:),info,init='Z')
+        if (do_timings) call psb_toc(poly_sv)
         call psb_geaxpby((sone/sm%rho_ba),ty,szero,r,desc_data,info)
-        if (.false.) then
-          call psb_geaxpby((sone/theta),r,szero,tz,desc_data,info)
-          call psb_geaxpby(sone,tz,sone,tx,desc_data,info)
-        else
-          call psb_abgdxyz((sone/theta),szero,sone,sone,r,tz,tx,desc_data,info)
-        end if
+        if (do_timings) call psb_tic(poly_vect)
+        call psb_upd_xyz((sone/theta),szero,sone,sone,r,tz,tx,desc_data,info)
+        if (do_timings) call psb_toc(poly_vect)
 
         ! tz == d
         do i=1, sm%pdegree-1
           !
-           !
+          !
           ! r_{k-1} = r_k - (1/rho(BA)) B A d_k
+          if (do_timings) call psb_tic(poly_mv)
           call psb_spmm(sone,sm%pa,tz,szero,ty,desc_data,info,work=aux,trans=trans_)
+          if (do_timings) call psb_toc(poly_mv)
+          if (do_timings) call psb_tic(poly_sv)
           call sm%sv%apply(-(sone/sm%rho_ba),ty,sone,r,desc_data,trans_,aux,wv(5:),info,init='Z')
-
+          if (do_timings) call psb_toc(poly_sv)
           !
           ! d_{k+1} = (rho rho_old) d_k + 2(rho/delta) r_{k+1}
           rho = sone/(2*sigma - rho_old)
-          if (.false.) then
-            call psb_geaxpby((2*rho/delta),r,(rho*rho_old),tz,desc_data,info)
-            call psb_geaxpby(sone,tz,sone,tx,desc_data,info)
-          else
-            call psb_abgdxyz((2*rho/delta),(rho*rho_old),sone,sone,r,tz,tx,desc_data,info)
-          end if
-!!$          res  = psb_genrm2(r,desc_data,info)
-!!$          write(0,*) 'Polynomial smoother NEW ',i,res
-          ! x_k =  x_{k-1} + z_k
+          if (do_timings) call psb_tic(poly_vect)
+          call psb_upd_xyz((2*rho/delta),(rho*rho_old),sone,sone,r,tz,tx,desc_data,info)
+          if (do_timings) call psb_toc(poly_vect)
           rho_old = rho
         end do
       end block
-
-
+      if (do_timings) call psb_toc(poly_3)
     case default
       info=psb_err_internal_error_
       call psb_errpush(info,name,&

amgprec/impl/smoother/amg_s_poly_smoother_bld.f90

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

amgprec/impl/smoother/amg_s_poly_smoother_cseti.f90

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

amgprec/impl/smoother/amg_s_poly_smoother_descr.f90

@@ -77,17 +77,17 @@ subroutine amg_s_poly_smoother_descr(sm,info,iout,coarse,prefix)
 
   write(iout_,*) trim(prefix_), '  Polynomial smoother   '
   select case(sm%variant)
-  case(amg_poly_lottes_)
-    write(iout_,*) trim(prefix_), '        variant:        ','POLY_LOTTES'
+  case(amg_cheb_4_)
+    write(iout_,*) trim(prefix_), '        variant:        ','CHEB_4'
     write(iout_,*) trim(prefix_), '         Degree:        ',sm%pdegree
     write(iout_,*) trim(prefix_), '         rho_ba:        ',sm%rho_ba
-  case(amg_poly_lottes_beta_)
-    write(iout_,*) trim(prefix_), '        variant:        ','POLY_LOTTES_BETA'
+  case(amg_cheb_4_opt_)
+    write(iout_,*) trim(prefix_), '        variant:        ','CHEB_4_OPT'
     write(iout_,*) trim(prefix_), '         Degree:        ',sm%pdegree
     write(iout_,*) trim(prefix_), '         rho_ba:        ',sm%rho_ba
     if (allocated(sm%poly_beta)) write(iout_,*) trim(prefix_), '   Coefficients:        ',sm%poly_beta(1:sm%pdegree)
-  case(amg_poly_new_)
-    write(iout_,*) trim(prefix_), '        variant:        ','POLY_NEW'
+  case(amg_cheb_1_opt_)
+    write(iout_,*) trim(prefix_), '        variant:        ','CHEB_1_OPT'
     write(iout_,*) trim(prefix_), '         Degree:        ',sm%pdegree
     write(iout_,*) trim(prefix_), '         rho_ba:        ',sm%rho_ba
     write(iout_,*) trim(prefix_), '    Coefficient:        ',sm%cf_a

amgprec/impl/smoother/amg_z_jac_smoother_apply_vect.f90

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

amgprec/impl/solver/amg_c_ilu_solver_bld.f90

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

amgprec/impl/solver/amg_c_jac_solver_apply.f90

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

amgprec/impl/solver/amg_c_jac_solver_apply_vect.f90

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

amgprec/impl/solver/amg_c_krm_solver_impl.f90

@@ -169,7 +169,7 @@ subroutine amg_c_krm_solver_apply_vect(alpha,sv,x,beta,y,desc_data,&
      & trans,work,wv,info,init,initu)
 
   use psb_base_mod
-  use psb_krylov_mod
+  use psb_linsolve_mod
   use amg_c_krm_solver, amg_protect_name => amg_c_krm_solver_apply_vect
 
   Implicit None

amgprec/impl/solver/amg_d_ilu_solver_bld.f90

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

amgprec/impl/solver/amg_d_jac_solver_apply.f90

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

amgprec/impl/solver/amg_d_jac_solver_apply_vect.f90

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

amgprec/impl/solver/amg_d_krm_solver_impl.f90

@@ -169,7 +169,7 @@ subroutine amg_d_krm_solver_apply_vect(alpha,sv,x,beta,y,desc_data,&
      & trans,work,wv,info,init,initu)
 
   use psb_base_mod
-  use psb_krylov_mod
+  use psb_linsolve_mod
   use amg_d_krm_solver, amg_protect_name => amg_d_krm_solver_apply_vect
 
   Implicit None

amgprec/impl/solver/amg_s_ilu_solver_bld.f90

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

amgprec/impl/solver/amg_s_jac_solver_apply.f90

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

amgprec/impl/solver/amg_s_jac_solver_apply_vect.f90

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

amgprec/impl/solver/amg_s_krm_solver_impl.f90

@@ -169,7 +169,7 @@ subroutine amg_s_krm_solver_apply_vect(alpha,sv,x,beta,y,desc_data,&
      & trans,work,wv,info,init,initu)
 
   use psb_base_mod
-  use psb_krylov_mod
+  use psb_linsolve_mod
   use amg_s_krm_solver, amg_protect_name => amg_s_krm_solver_apply_vect
 
   Implicit None

amgprec/impl/solver/amg_z_ilu_solver_bld.f90

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

amgprec/impl/solver/amg_z_jac_solver_apply.f90

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

amgprec/impl/solver/amg_z_jac_solver_apply_vect.f90

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

amgprec/impl/solver/amg_z_krm_solver_impl.f90

@@ -169,7 +169,7 @@ subroutine amg_z_krm_solver_apply_vect(alpha,sv,x,beta,y,desc_data,&
      & trans,work,wv,info,init,initu)
 
   use psb_base_mod
-  use psb_krylov_mod
+  use psb_linsolve_mod
   use amg_z_krm_solver, amg_protect_name => amg_z_krm_solver_apply_vect
 
   Implicit None

cbind/amgprec/amg_dprec_cbind_mod.F90

@@ -264,7 +264,7 @@ contains
        & ah,ph,bh,xh,cdh,options) bind(c) result(res)
     use psb_base_mod
     use psb_prec_mod
-    use psb_krylov_mod
+    use psb_linsolve_mod
     use psb_prec_cbind_mod
     use psb_dkrylov_cbind_mod
     implicit none
@@ -285,7 +285,7 @@ contains
        & ah,ph,bh,xh,eps,cdh,itmax,iter,err,itrace,irst,istop) bind(c) result(res)
     use psb_base_mod
     use psb_prec_mod
-    use psb_krylov_mod
+    use psb_linsolve_mod
     use psb_objhandle_mod
     use psb_prec_cbind_mod
     use psb_base_string_cbind_mod

cbind/amgprec/amg_zprec_cbind_mod.F90

@@ -264,7 +264,7 @@ contains
        & ah,ph,bh,xh,cdh,options) bind(c) result(res)
     use psb_base_mod
     use psb_prec_mod
-    use psb_krylov_mod
+    use psb_linsolve_mod
     use psb_prec_cbind_mod
     use psb_zkrylov_cbind_mod
     implicit none
@@ -285,7 +285,7 @@ contains
        & ah,ph,bh,xh,eps,cdh,itmax,iter,err,itrace,irst,istop) bind(c) result(res)
     use psb_base_mod
     use psb_prec_mod
-    use psb_krylov_mod
+    use psb_linsolve_mod
     use psb_objhandle_mod
     use psb_prec_cbind_mod
     use psb_base_string_cbind_mod

cbind/test/pargen/Makefile

@@ -9,8 +9,8 @@ HERE=.
 FINCLUDES=$(FMFLAG). $(FMFLAG)$(LIBDIR) $(FMFLAG)$(PSBLAS_INCDIR)
 #PSBLAS_LIBS= -L$(PSBLAS_LIBDIR) -L$(LIBDIR)  $(CPSBLAS_LIB) $(PSBLAS_LIB)
 # -lpsb_krylov_cbind -lpsb_prec_cbind -lpsb_base_cbind 
-PSBC_LIBS= -L$(PSBLAS_LIBDIR) -lpsb_cbind  -lpsb_krylov -lpsb_prec 
-MLDC_LIBS=-L$(LIBDIR) -lmld_cbind  -lmld_prec 
+PSBC_LIBS= -L$(PSBLAS_LIBDIR) -lpsb_cbind  -lpsb_linsolve -lpsb_prec 
+AMGC_LIBS=-L$(LIBDIR) -lamg_cbind  -lamg_prec 
 #
 # Compilers and such
 #
@@ -23,15 +23,15 @@ EXEDIR=./runs
 #UMFLIBS=-lumfpack -lamd -lcholmod -lcolamd -lcamd -lccolamd -L/usr/include/suitesparse
 #UMFFLAGS=-DHave_UMF_ -I/usr/include/suitesparse
 
- all:  mldec
+ all:  amgec
 
-mldec: mldec.o
-	$(MPFC) mldec.o  -o mldec $(MLDC_LIBS)  $(PSBC_LIBS) $(PSBCLDLIBS) $(PSBLAS_LIBS) \
+amgec: amgec.o
+	$(MPFC) amgec.o  -o amgec $(AMGC_LIBS)  $(PSBC_LIBS) $(PSBCLDLIBS) $(PSBLAS_LIBS) \
 	  $(UMFLIBS) $(PSBLDLIBS) $(LDLIBS) -lm -lgfortran
 # \
 #	-lifcore -lifcoremt -lguide -limf -lirc -lintlc -lcxaguard -L/opt/intel/fc/10.0.023/lib/ -lm
 
-	/bin/mv mldec $(EXEDIR)
+	/bin/mv amgec $(EXEDIR)
 
 .f90.o:
 	$(MPFC) $(F90COPT) $(FINCLUDES) $(FDEFINES) -c $<
@@ -40,13 +40,13 @@ mldec: mldec.o
 
 
 clean: 
-	/bin/rm -f mldec.o  $(EXEDIR)/mldec
+	/bin/rm -f amgec.o  $(EXEDIR)/amgec
 verycleanlib: 
 	(cd ../..; make veryclean)
 lib:
 	(cd ../../; make library)
 
 tests: all
-	cd runs ; ./mldec < mlde.inp
+	cd runs ; ./amgec < amge.inp
 
 

cbind/test/pargen/amgec.c

@@ -77,7 +77,7 @@
 #include <math.h>
 
 #include "psb_base_cbind.h"
-#include "mld_cbind.h"
+#include "amg_cbind.h"
 
 
 double  a1(double x, double y, double  z)
@@ -123,7 +123,7 @@ double g(double x, double y, double z)
 
 #define NBMAX       20
 
-psb_i_t matgen(psb_i_t ictxt, psb_i_t nl, psb_i_t idim, psb_l_t vl[],
+psb_i_t matgen(psb_c_ctxt cctxt, psb_i_t nl, psb_i_t idim, psb_l_t vl[],
 	       psb_c_dspmat *ah,psb_c_descriptor *cdh,
 	       psb_c_dvector *xh, psb_c_dvector *bh, psb_c_dvector *rh)
 {
@@ -135,7 +135,7 @@ psb_i_t matgen(psb_i_t ictxt, psb_i_t nl, psb_i_t idim, psb_l_t vl[],
   psb_l_t irow[10*NBMAX], icol[10*NBMAX];
 
   info = 0;
-  psb_c_info(ictxt,&iam,&np);
+  psb_c_info(cctxt,&iam,&np);
   deltah = (double) 1.0/(idim+1);
   sqdeltah = deltah*deltah;
   deltah2  = 2.0* deltah;
@@ -253,11 +253,12 @@ void get_hparm(FILE *fp, char *val)
 
 int main(int argc, char *argv[])
 {
-  psb_i_t ictxt, iam, np;
+  psb_c_ctxt *cctxt;
+  psb_i_t iam, np;
   char methd[40], ptype[40], afmt[8];
   psb_i_t nparms;
   psb_i_t idim,info,istop,itmax,itrace,irst,iter,ret;
-  mld_c_dprec *ph;
+  amg_c_dprec *ph;
   psb_c_dspmat *ah;
   psb_c_dvector *bh, *xh, *rh; 
   psb_i_t nb,nlr, nl;
@@ -269,12 +270,13 @@ int main(int argc, char *argv[])
   psb_c_descriptor *cdh;
   FILE *vectfile;
   
-  ictxt = psb_c_init();
-  psb_c_info(ictxt,&iam,&np);
+  cctxt = psb_c_new_ctxt();
+  psb_c_init(cctxt);
+  psb_c_info(*cctxt,&iam,&np);
   fprintf(stdout,"Initialization: am %d of %d\n",iam,np);
 
   fflush(stdout);
-  psb_c_barrier(ictxt);
+  psb_c_barrier(*cctxt);
   if (iam == 0) {
     get_iparm(stdin,&nparms);
     get_hparm(stdin,methd);
@@ -287,17 +289,17 @@ int main(int argc, char *argv[])
     get_iparm(stdin,&irst);
   } 
   /* Now broadcast the values, and check they're OK */
-  psb_c_ibcast(ictxt,1,&nparms,0);
-  psb_c_hbcast(ictxt,methd,0);
-  psb_c_hbcast(ictxt,ptype,0);
-  psb_c_hbcast(ictxt,afmt,0);
-  psb_c_ibcast(ictxt,1,&idim,0);
-  psb_c_ibcast(ictxt,1,&istop,0);
-  psb_c_ibcast(ictxt,1,&itmax,0);
-  psb_c_ibcast(ictxt,1,&itrace,0);
-  psb_c_ibcast(ictxt,1,&irst,0);
+  psb_c_ibcast(*cctxt,1,&nparms,0);
+  psb_c_hbcast(*cctxt,methd,0);
+  psb_c_hbcast(*cctxt,ptype,0);
+  psb_c_hbcast(*cctxt,afmt,0);
+  psb_c_ibcast(*cctxt,1,&idim,0);
+  psb_c_ibcast(*cctxt,1,&istop,0);
+  psb_c_ibcast(*cctxt,1,&itmax,0);
+  psb_c_ibcast(*cctxt,1,&itrace,0);
+  psb_c_ibcast(*cctxt,1,&irst,0);
   
-  psb_c_barrier(ictxt);
+  psb_c_barrier(*cctxt);
 
   cdh=psb_c_new_descriptor();
   psb_c_set_index_base(0);
@@ -310,15 +312,15 @@ int main(int argc, char *argv[])
   fprintf(stderr,"%d: Input data %d %ld %d %d\n",iam,idim,ng,nb, nl);
   if ((vl=malloc(nb*sizeof(psb_l_t)))==NULL) {
     fprintf(stderr,"On %d: malloc failure\n",iam);
-    psb_c_abort(ictxt);
+    psb_c_abort(*cctxt);
   }
   i = ((psb_l_t)iam) * nb;
   for (k=0; k<nl; k++)
     vl[k] = i+k; 
 
-  if ((info=psb_c_cdall_vl(nl,vl,ictxt,cdh))!=0) {
+  if ((info=psb_c_cdall_vl(nl,vl,*cctxt,cdh))!=0) {
     fprintf(stderr,"From cdall: %d\nBailing out\n",info);
-    psb_c_abort(ictxt);
+    psb_c_abort(*cctxt);
   }
 
   bh  = psb_c_new_dvector();
@@ -337,25 +339,25 @@ int main(int argc, char *argv[])
 
   
   /* Matrix generation  */
-  if (matgen(ictxt,nl,idim,vl,ah,cdh,xh,bh,rh) != 0) {
+  if (matgen(*cctxt,nl,idim,vl,ah,cdh,xh,bh,rh) != 0) {
     fprintf(stderr,"Error during matrix build loop\n");
-    psb_c_abort(ictxt);
+    psb_c_abort(*cctxt);
   }    
-  psb_c_barrier(ictxt);
+  psb_c_barrier(*cctxt);
   /* Set up the preconditioner */ 
-  ph  = mld_c_dprec_new();
-  mld_c_dprecinit(ictxt,ph,ptype);
-  mld_c_dprecseti(ph,"SMOOTHER_SWEEPS",2);
-  mld_c_dprecseti(ph,"SUB_FILLIN",1);
-  mld_c_dprecsetc(ph,"COARSE_SOLVE","BJAC");
-  mld_c_dprecsetc(ph,"COARSE_SUBSOLVE","ILU");
-  mld_c_dprecseti(ph,"COARSE_FILLIN",0);
-  if ((ret=mld_c_dhierarchy_build(ah,cdh,ph))!=0)
+  ph  = amg_c_dprec_new();
+  amg_c_dprecinit(*cctxt,ph,ptype);
+  amg_c_dprecseti(ph,"SMOOTHER_SWEEPS",2);
+  amg_c_dprecseti(ph,"SUB_FILLIN",1);
+  amg_c_dprecsetc(ph,"COARSE_SOLVE","BJAC");
+  amg_c_dprecsetc(ph,"COARSE_SUBSOLVE","ILU");
+  amg_c_dprecseti(ph,"COARSE_FILLIN",0);
+  if ((ret=amg_c_dhierarchy_build(ah,cdh,ph))!=0)
     fprintf(stderr,"From hierarchy_build: %d\n",ret);
-  if ((ret=mld_c_dsmoothers_build(ah,cdh,ph))!=0)
+  if ((ret=amg_c_dsmoothers_build(ah,cdh,ph))!=0)
     fprintf(stderr,"From smoothers_build: %d\n",ret);
 
-  psb_c_barrier(ictxt);
+  psb_c_barrier(*cctxt);
   /* Set up the solver options */ 
   psb_c_DefaultSolverOptions(&options);
   options.eps    = 1.e-6;
@@ -365,7 +367,7 @@ int main(int argc, char *argv[])
   options.istop  = istop;
   psb_c_seterraction_ret();
   t1=psb_c_wtime();
-  ret=mld_c_dkrylov(methd,ah,ph,bh,xh,cdh,&options); 
+  ret=amg_c_dkrylov(methd,ah,ph,bh,xh,cdh,&options); 
   t2=psb_c_wtime();
   iter = options.iter;
   err  = options.err;
@@ -413,20 +415,20 @@ int main(int argc, char *argv[])
   /* Clean up memory */ 
   if ((info=psb_c_dgefree(xh,cdh))!=0) {
     fprintf(stderr,"From dgefree: %d\nBailing out\n",info);
-    psb_c_abort(ictxt);
+    psb_c_abort(*cctxt);
   }
   if ((info=psb_c_dgefree(bh,cdh))!=0) {
     fprintf(stderr,"From dgefree: %d\nBailing out\n",info);
-    psb_c_abort(ictxt);
+    psb_c_abort(*cctxt);
   }
   if ((info=psb_c_dgefree(rh,cdh))!=0) {
     fprintf(stderr,"From dgefree: %d\nBailing out\n",info);
-    psb_c_abort(ictxt);
+    psb_c_abort(*cctxt);
   }
 
   if ((info=psb_c_cdfree(cdh))!=0) {
     fprintf(stderr,"From cdfree: %d\nBailing out\n",info);
-    psb_c_abort(ictxt);
+    psb_c_abort(*cctxt);
   }
   //fprintf(stderr,"pointer  from cdfree: %p\n",cdh->descriptor);
  
@@ -440,6 +442,7 @@ int main(int argc, char *argv[])
 
   if (iam == 0) fprintf(stderr,"program completed successfully\n");
 
-  psb_c_barrier(ictxt);
-  psb_c_exit(ictxt);
+  psb_c_barrier(*cctxt);
+  psb_c_exit(*cctxt);
+  free(cctxt);
 }

config/pac.m4

@@ -409,7 +409,7 @@ save_LDFLAGS=$LDFLAGS;
 ## dnl AC_MSG_NOTICE([psblas dir $pac_cv_psblas_dir])
 ##  PSBLAS_LIBS="-L$pac_cv_psblas_dir/lib"
 ## fi
-PSBLAS_LIBS="-lpsb_krylov -lpsb_prec -lpsb_util -lpsb_base -L$PSBLAS_LIBDIR"
+PSBLAS_LIBS="-lpsb_linsolve -lpsb_prec -lpsb_util -lpsb_base -L$PSBLAS_LIBDIR"
 LDFLAGS=" $PSBLAS_LIBS $save_LDFLAGS"
 
 dnl ac_compile='${MPIFC-$FC} -c -o conftest${ac_objext} $FMFLAG$PSBLAS_DIR/include $FMFLAG$PSBLAS_DIR/lib conftest.$ac_ext  1>&5'
@@ -484,7 +484,7 @@ dnl AC_MSG_NOTICE([psblas dir $pac_cv_psblas_dir])
  PSBLAS_INCLUDES="$FMFLAG$pac_cv_psblas_dir/modules $PSBLAS_INCLUDES"
 fi
 FCFLAGS=" $PSBLAS_INCLUDES $save_FCFLAGS"
-PSBLAS_LIBS="-lpsb_krylov -lpsb_prec -lpsb_util -lpsb_base $PSBLAS_LIBS"
+PSBLAS_LIBS="-lpsb_linsolve -lpsb_prec -lpsb_util -lpsb_base $PSBLAS_LIBS"
 LDFLAGS=" $PSBLAS_LIBS $save_LDFLAGS"
 
 dnl ac_compile='${MPIFC-$FC} -c -o conftest${ac_objext} $FMFLAG$PSBLAS_DIR/include $FMFLAG$PSBLAS_DIR/lib conftest.$ac_ext  1>&5'

configure

@@ -7417,7 +7417,7 @@ then
    FIFLAG="-I"
    BASEMODNAME=PSB_BASE_MOD
    PRECMODNAME=PSB_PREC_MOD
-   METHDMODNAME=PSB_KRYLOV_MOD
+   METHDMODNAME=PSB_LINSOLVE_MOD
    UTILMODNAME=PSB_UTIL_MOD
 
 else
@@ -7550,7 +7550,7 @@ printf "%s\n" "$ax_cv_f90_modflag" >&6; }
    FIFLAG=-I
    BASEMODNAME=psb_base_mod
    PRECMODNAME=psb_prec_mod
-   METHDMODNAME=psb_krylov_mod
+   METHDMODNAME=psb_linsolve_mod
    UTILMODNAME=psb_util_mod
 fi
 
@@ -7683,7 +7683,7 @@ save_LDFLAGS=$LDFLAGS;
 ## dnl AC_MSG_NOTICE([psblas dir $pac_cv_psblas_dir])
 ##  PSBLAS_LIBS="-L$pac_cv_psblas_dir/lib"
 ## fi
-PSBLAS_LIBS="-lpsb_krylov -lpsb_prec -lpsb_util -lpsb_base -L$PSBLAS_LIBDIR"
+PSBLAS_LIBS="-lpsb_linsolve -lpsb_prec -lpsb_util -lpsb_base -L$PSBLAS_LIBDIR"
 LDFLAGS=" $PSBLAS_LIBS $save_LDFLAGS"
 
 ac_link='${MPIFC-$FC} -o conftest${ac_exeext} $FCFLAGS conftest.$ac_ext $LDFLAGS $LIBS  1>&5'
@@ -7787,7 +7787,7 @@ elif test "x$pac_cv_psblas_dir" != "x"; then
  PSBLAS_INCLUDES="$FMFLAG$pac_cv_psblas_dir/modules $PSBLAS_INCLUDES"
 fi
 FCFLAGS=" $PSBLAS_INCLUDES $save_FCFLAGS"
-PSBLAS_LIBS="-lpsb_krylov -lpsb_prec -lpsb_util -lpsb_base $PSBLAS_LIBS"
+PSBLAS_LIBS="-lpsb_linsolve -lpsb_prec -lpsb_util -lpsb_base $PSBLAS_LIBS"
 LDFLAGS=" $PSBLAS_LIBS $save_LDFLAGS"
 
 

configure.ac

@@ -525,7 +525,7 @@ then
    FIFLAG="-I"
    BASEMODNAME=PSB_BASE_MOD
    PRECMODNAME=PSB_PREC_MOD
-   METHDMODNAME=PSB_KRYLOV_MOD
+   METHDMODNAME=PSB_LINSOLVE_MOD
    UTILMODNAME=PSB_UTIL_MOD
 
 else 
@@ -536,7 +536,7 @@ else
    FIFLAG=-I
    BASEMODNAME=psb_base_mod
    PRECMODNAME=psb_prec_mod
-   METHDMODNAME=psb_krylov_mod
+   METHDMODNAME=psb_linsolve_mod
    UTILMODNAME=psb_util_mod
 fi
 

samples/advanced/fileread/Makefile

@@ -3,7 +3,7 @@ AMGINCDIR=$(AMGDIR)/include
 include $(AMGINCDIR)/Make.inc.amg4psblas
 AMGMODDIR=$(AMGDIR)/modules
 AMGLIBDIR=$(AMGDIR)/lib
-AMG_LIBS=-L$(AMGLIBDIR) -lpsb_krylov -lamg_prec -lpsb_prec
+AMG_LIBS=-L$(AMGLIBDIR) -lpsb_linsolve -lamg_prec -lpsb_prec
 FINCLUDES=$(FMFLAG). $(FMFLAG)$(AMGMODDIR) $(FMFLAG)$(AMGINCDIR) $(PSBLAS_INCLUDES) $(FIFLAG).
 
 DFSOBJS=amg_df_sample.o data_input.o

samples/advanced/fileread/amg_cf_sample.f90

@@ -38,7 +38,7 @@
 program amg_cf_sample
   use psb_base_mod
   use amg_prec_mod
-  use psb_krylov_mod
+  use psb_linsolve_mod
   use psb_util_mod
   use data_input
   implicit none

samples/advanced/fileread/amg_df_sample.f90

@@ -38,7 +38,7 @@
 program amg_df_sample
   use psb_base_mod
   use amg_prec_mod
-  use psb_krylov_mod
+  use psb_linsolve_mod
   use psb_util_mod
   use data_input
   implicit none

samples/advanced/fileread/amg_sf_sample.f90

@@ -38,7 +38,7 @@
 program amg_sf_sample
   use psb_base_mod
   use amg_prec_mod
-  use psb_krylov_mod
+  use psb_linsolve_mod
   use psb_util_mod
   use data_input
   implicit none

samples/advanced/fileread/amg_zf_sample.f90

@@ -38,7 +38,7 @@
 program amg_zf_sample
   use psb_base_mod
   use amg_prec_mod
-  use psb_krylov_mod
+  use psb_linsolve_mod
   use psb_util_mod
   use data_input
   implicit none

samples/advanced/pdegen/Makefile

@@ -3,37 +3,42 @@ AMGINCDIR=$(AMGDIR)/include
 include $(AMGINCDIR)/Make.inc.amg4psblas
 AMGMODDIR=$(AMGDIR)/modules
 AMGLIBDIR=$(AMGDIR)/lib
-AMG_LIBS=-L$(AMGLIBDIR) -lpsb_krylov -lamg_prec -lpsb_prec 
+AMG_LIBS=-L$(AMGLIBDIR) -lpsb_linsolve -lamg_prec -lpsb_prec 
 FINCLUDES=$(FMFLAG). $(FMFLAG)$(AMGMODDIR) $(FMFLAG)$(AMGINCDIR) $(PSBLAS_INCLUDES) $(FIFLAG).
 
 LINKOPT=
 EXEDIR=./runs
-DGEN2D=amg_d_pde2d_base_mod.o amg_d_pde2d_exp_mod.o amg_d_pde2d_gauss_mod.o amg_d_pde2d_box_mod.o
-DGEN3D=amg_d_pde3d_base_mod.o amg_d_pde3d_exp_mod.o amg_d_pde3d_gauss_mod.o amg_d_pde3d_box_mod.o 	
-SGEN2D=amg_s_pde2d_base_mod.o amg_s_pde2d_exp_mod.o amg_s_pde2d_gauss_mod.o amg_s_pde2d_box_mod.o
-SGEN3D=amg_s_pde3d_base_mod.o amg_s_pde3d_exp_mod.o amg_s_pde3d_gauss_mod.o amg_s_pde3d_box_mod.o 	
+DGEN2D=amg_d_pde2d_poisson_mod.o amg_d_pde2d_exp_mod.o \
+  amg_d_pde2d_gauss_mod.o amg_d_pde2d_box_mod.o
+DGEN3D=amg_d_pde3d_poisson_mod.o amg_d_pde3d_exp_mod.o \
+  amg_d_pde3d_gauss_mod.o amg_d_pde3d_box_mod.o 	
+SGEN2D=amg_s_pde2d_poisson_mod.o amg_s_pde2d_exp_mod.o \
+  amg_s_pde2d_gauss_mod.o amg_s_pde2d_box_mod.o
+SGEN3D=amg_s_pde3d_poisson_mod.o amg_s_pde3d_exp_mod.o \
+  amg_s_pde3d_gauss_mod.o amg_s_pde3d_box_mod.o 	
 
 all: amg_s_pde3d amg_d_pde3d amg_s_pde2d amg_d_pde2d
 
 amg_d_pde3d: amg_d_pde3d.o amg_d_genpde_mod.o $(DGEN3D) data_input.o
-	$(FLINK) $(LINKOPT) amg_d_pde3d.o  amg_d_genpde_mod.o $(DGEN3D) data_input.o -o amg_d_pde3d $(AMG_LIBS) $(PSBLAS_LIBS) $(LDLIBS)
+	$(FLINK) $(LINKOPT) amg_d_pde3d.o  amg_d_genpde_mod.o $(DGEN3D) data_input.o \
+         -o amg_d_pde3d $(AMG_LIBS) $(PSBLAS_LIBS) $(LDLIBS)
 	/bin/mv amg_d_pde3d $(EXEDIR)
 
 amg_s_pde3d: amg_s_pde3d.o amg_s_genpde_mod.o $(SGEN3D) data_input.o
-	$(FLINK) $(LINKOPT) amg_s_pde3d.o  amg_s_genpde_mod.o  $(SGEN3D) data_input.o -o amg_s_pde3d $(AMG_LIBS) $(PSBLAS_LIBS) $(LDLIBS)
+	$(FLINK) $(LINKOPT) amg_s_pde3d.o  amg_s_genpde_mod.o  $(SGEN3D) data_input.o \
+         -o amg_s_pde3d $(AMG_LIBS) $(PSBLAS_LIBS) $(LDLIBS)
 	/bin/mv amg_s_pde3d $(EXEDIR)
 
 amg_d_pde2d: amg_d_pde2d.o amg_d_genpde_mod.o  $(DGEN2D) data_input.o
-	$(FLINK) $(LINKOPT) amg_d_pde2d.o amg_d_genpde_mod.o $(DGEN2D) data_input.o -o amg_d_pde2d $(AMG_LIBS) $(PSBLAS_LIBS) $(LDLIBS)
+	$(FLINK) $(LINKOPT) amg_d_pde2d.o amg_d_genpde_mod.o $(DGEN2D) data_input.o \
+         -o amg_d_pde2d $(AMG_LIBS) $(PSBLAS_LIBS) $(LDLIBS)
 	/bin/mv amg_d_pde2d $(EXEDIR)
 
 amg_s_pde2d: amg_s_pde2d.o amg_s_genpde_mod.o $(SGEN2D) data_input.o
-	$(FLINK) $(LINKOPT) amg_s_pde2d.o amg_s_genpde_mod.o $(SGEN2D) data_input.o -o amg_s_pde2d $(AMG_LIBS) $(PSBLAS_LIBS) $(LDLIBS)
+	$(FLINK) $(LINKOPT) amg_s_pde2d.o amg_s_genpde_mod.o $(SGEN2D) data_input.o \
+          -o amg_s_pde2d $(AMG_LIBS) $(PSBLAS_LIBS) $(LDLIBS)
 	/bin/mv amg_s_pde2d $(EXEDIR)
 
-amg_d_pde3d_rebld: amg_d_pde3d_rebld.o  data_input.o
-	$(FLINK) $(LINKOPT) amg_d_pde3d_rebld.o  data_input.o -o amg_d_pde3d_rebld $(AMG_LIBS) $(PSBLAS_LIBS) $(LDLIBS)
-	/bin/mv amg_d_pde3d_rebld $(EXEDIR)
 
 amg_d_pde3d.o amg_s_pde3d.o amg_d_pde2d.o amg_s_pde2d.o: data_input.o
 
@@ -42,6 +47,11 @@ amg_s_pde3d.o: amg_s_genpde_mod.o $(SGEN3D)
 amg_d_pde2d.o: amg_d_genpde_mod.o $(DGEN2D)
 amg_s_pde2d.o: amg_s_genpde_mod.o $(SGEN2D)
 
+amg_d_genpde_mod.o: $(DGEN3D)
+amg_s_genpde_mod.o: $(SGEN3D)
+amg_d_genpde_mod.o: $(DGEN2D)
+amg_s_genpde_mod.o: $(SGEN2D)
+
 check: all
 	cd runs && ./amg_d_pde2d <amg_pde2d.inp && ./amg_s_pde2d<amg_pde2d.inp