OpenMP loop in samples data generation

Added timers to build phases
Debug for matchboxp
2026-10-07 15:15:07 +00:00 · 2023-06-05 11:46:39 +02:00 · 2023-06-02 11:37:58 +02:00 · 2023-02-10 07:53:04 -05:00 · 2023-02-10 07:50:58 -05:00 · 2022-12-13 06:44:12 -05:00
40 changed files with 3733 additions and 573 deletions
@@ -1,4 +1,3 @@
-
                         AMG4PSBLAS  
 Algebraic Multigrid Package  based on PSBLAS (Parallel Sparse BLAS version 3.8)
    
@@ -143,9 +143,10 @@ contains
    type(psb_ld_coo_sparse_mat) :: tmpcoo
    logical  :: display_out_, print_out_, reproducible_
    logical, parameter  :: dump=.false., debug=.false., dump_mate=.false., &
-         & debug_ilaggr=.false., debug_sync=.false.
+         & debug_ilaggr=.false., debug_sync=.false., debug_mate=.false.
    integer(psb_ipk_), save :: idx_bldmtc=-1, idx_phase1=-1, idx_phase2=-1, idx_phase3=-1
    logical, parameter :: do_timings=.true.
+    integer, parameter :: ilaggr_neginit=-1, ilaggr_nonlocal=-2

    ictxt = desc_a%get_ctxt()
    call psb_info(ictxt,iam,np)
@@ -187,7 +188,7 @@ contains
    call desc_a%l2gip(ilv,info,owned=.false.)

    call psb_geall(ilaggr,desc_a,info)
-    ilaggr = -1
+    ilaggr = ilaggr_neginit
    call psb_geasb(ilaggr,desc_a,info)
    nr = a%get_nrows()
    nc = a%get_ncols()
@@ -213,7 +214,20 @@ contains
      call psb_barrier(ictxt)
      if (iam == 0) write(0,*)' out from buildmatching:', info
    end if
-
+    if (debug_mate) then
+      block
+        integer(psb_lpk_), allocatable :: ckmate(:)
+        allocate(ckmate(nr))
+        ckmate(1:nr) = mate(1:nr)
+        call psb_msort(ckmate(1:nr))
+        do i=1,nr-1
+          if ((ckmate(i)>0) .and. (ckmate(i) == ckmate(i+1))) then
+            write(0,*) iam,' Duplicate mate entry  at',i,' :',ckmate(i)
+          end if
+        end do
+      end block
+    end if
+     
    if (info == 0) then
      if (do_timings) call psb_tic(idx_phase2)
      if (debug_sync) then
@@ -259,7 +273,7 @@ contains
          cycle
        else

-          if (ilaggr(k) == -1) then
+          if (ilaggr(k) == ilaggr_neginit) then

            wk   = w(k)
            widx = w(idx)
@@ -267,7 +281,7 @@ contains
            nrmagg = wmax*sqrt((wk/wmax)**2+(widx/wmax)**2)
            if (nrmagg > epsilon(nrmagg)) then
              if (idx <= nr) then
-                if (ilaggr(idx) == -1) then
+                if (ilaggr(idx) == ilaggr_neginit) then
                  ! Now, if both vertices are local, the aggregate is local
                  ! (kinda obvious).
                  nlaggr(iam) = nlaggr(iam) + 1
@@ -275,6 +289,9 @@ contains
                  ilaggr(idx) = nlaggr(iam)
                  wtemp(k)    = w(k)/nrmagg
                  wtemp(idx)  = w(idx)/nrmagg
+                else
+                  write(0,*) iam,' Inconsistent mate? ',k,mate(k),idx,&
+                       &mate(idx),ilaggr(idx)
                end if
                nlpairs = nlpairs+1
              else if (idx <= nc) then
@@ -294,7 +311,7 @@ contains
                    ilaggr(k)   = nlaggr(iam)
                    nlpairs = nlpairs+1
                  else
-                    ilaggr(k) = -2
+                    ilaggr(k) = ilaggr_nonlocal
                  end if
                else
                  ! Use a statistically unbiased tie-breaking rule,
@@ -309,7 +326,7 @@ contains
                    ilaggr(k)   = nlaggr(iam)
                    nlpairs = nlpairs+1
                  else
-                    ilaggr(k) = -2
+                    ilaggr(k) = ilaggr_nonlocal
                  end if
                end if
              end if
@@ -325,6 +342,12 @@ contains
                nlsingl     = nlsingl + 1
              end if
            end if
+            if (ilaggr(k) == ilaggr_neginit) then
+              write(0,*) iam,' Error: no update to ',k,mate(k),&
+                   & abs(w(k)),nrmagg,epsilon(nrmagg),wtemp(k)
+            end if
+          else
+            if (ilaggr(k)<0) write(0,*) 'Strange? ',k,ilaggr(k)
          end if
        end if
      end do
@@ -332,7 +355,7 @@ contains
      if (do_timings) call psb_tic(idx_phase3)

      ! Ok, now compute offsets, gather halo and fix non-local
-      ! aggregates  (those where ilaggr == -2)
+      ! aggregates  (those where ilaggr == ilaggr_nonlocal)
      call psb_sum(ictxt,nlaggr)
      ntaggr  = sum(nlaggr(0:np-1))
      naggrm1 = sum(nlaggr(0:iam-1))
@@ -347,7 +370,7 @@ contains
      call psb_halo(wtemp,desc_a,info)
      ! Cleanup as yet unmarked entries
      do k=1,nr
-        if (ilaggr(k) == -2) then
+        if (ilaggr(k) == ilaggr_nonlocal) then
          idx = mate(k)
          if (idx > nr) then
            i   = ilaggr(idx)
@@ -359,9 +382,14 @@ contains
          else
            write(0,*) 'Error : unresolved (paired) index ',k,idx,i,nr,nc, ilv(k),ilv(idx)
          end if
-        end if
-        if (ilaggr(k) <0) then
-          write(0,*) 'Matchboxp: Funny number: ',k,ilv(k),ilaggr(k),wtemp(k)
+        else if (ilaggr(k) <0) then
+          write(0,*) iam,'Matchboxp: Funny number: ',k,ilv(k),ilaggr(k),wtemp(k)
+          write(0,*) iam,'         :             : ',nr,nc,mate(k)
+          if (mate(k) <= nr) then
+            write(0,*) iam,'         :             : ',ilaggr(mate(k)),mate(mate(k)),&
+                 & ilv(k),ilv(mate(k)), ilv(mate(mate(k))),ilaggr(mate(mate(k)))
+          end if
+          flush(0)
        end if
      end do
      if (debug_sync) then
@@ -414,7 +442,7 @@ contains

      end block
      if (iam == 0) then
-        write(0,*) 'Matching statistics: Unmatched nodes ',&
+        write(0,*) iam,'Matching statistics: Unmatched nodes ',&
             & nunmatched,' Singletons:',nlsingl,' Pairs:',nlpairs
      end if

@@ -143,9 +143,10 @@ contains
    type(psb_ls_coo_sparse_mat) :: tmpcoo
    logical  :: display_out_, print_out_, reproducible_
    logical, parameter  :: dump=.false., debug=.false., dump_mate=.false., &
-         & debug_ilaggr=.false., debug_sync=.false.
+         & debug_ilaggr=.false., debug_sync=.false., debug_mate=.false.
    integer(psb_ipk_), save :: idx_bldmtc=-1, idx_phase1=-1, idx_phase2=-1, idx_phase3=-1
    logical, parameter :: do_timings=.true.
+    integer, parameter :: ilaggr_neginit=-1, ilaggr_nonlocal=-2

    ictxt = desc_a%get_ctxt()
    call psb_info(ictxt,iam,np)
@@ -187,7 +188,7 @@ contains
    call desc_a%l2gip(ilv,info,owned=.false.)

    call psb_geall(ilaggr,desc_a,info)
-    ilaggr = -1
+    ilaggr = ilaggr_neginit
    call psb_geasb(ilaggr,desc_a,info)
    nr = a%get_nrows()
    nc = a%get_ncols()
@@ -213,7 +214,20 @@ contains
      call psb_barrier(ictxt)
      if (iam == 0) write(0,*)' out from buildmatching:', info
    end if
-
+    if (debug_mate) then
+      block
+        integer(psb_lpk_), allocatable :: ckmate(:)
+        allocate(ckmate(nr))
+        ckmate(1:nr) = mate(1:nr)
+        call psb_msort(ckmate(1:nr))
+        do i=1,nr-1
+          if ((ckmate(i)>0) .and. (ckmate(i) == ckmate(i+1))) then
+            write(0,*) iam,' Duplicate mate entry  at',i,' :',ckmate(i)
+          end if
+        end do
+      end block
+    end if
+     
    if (info == 0) then
      if (do_timings) call psb_tic(idx_phase2)
      if (debug_sync) then
@@ -259,7 +273,7 @@ contains
          cycle
        else

-          if (ilaggr(k) == -1) then
+          if (ilaggr(k) == ilaggr_neginit) then

            wk   = w(k)
            widx = w(idx)
@@ -267,7 +281,7 @@ contains
            nrmagg = wmax*sqrt((wk/wmax)**2+(widx/wmax)**2)
            if (nrmagg > epsilon(nrmagg)) then
              if (idx <= nr) then
-                if (ilaggr(idx) == -1) then
+                if (ilaggr(idx) == ilaggr_neginit) then
                  ! Now, if both vertices are local, the aggregate is local
                  ! (kinda obvious).
                  nlaggr(iam) = nlaggr(iam) + 1
@@ -275,6 +289,9 @@ contains
                  ilaggr(idx) = nlaggr(iam)
                  wtemp(k)    = w(k)/nrmagg
                  wtemp(idx)  = w(idx)/nrmagg
+                else
+                  write(0,*) iam,' Inconsistent mate? ',k,mate(k),idx,&
+                       &mate(idx),ilaggr(idx)
                end if
                nlpairs = nlpairs+1
              else if (idx <= nc) then
@@ -294,7 +311,7 @@ contains
                    ilaggr(k)   = nlaggr(iam)
                    nlpairs = nlpairs+1
                  else
-                    ilaggr(k) = -2
+                    ilaggr(k) = ilaggr_nonlocal
                  end if
                else
                  ! Use a statistically unbiased tie-breaking rule,
@@ -309,7 +326,7 @@ contains
                    ilaggr(k)   = nlaggr(iam)
                    nlpairs = nlpairs+1
                  else
-                    ilaggr(k) = -2
+                    ilaggr(k) = ilaggr_nonlocal
                  end if
                end if
              end if
@@ -325,6 +342,12 @@ contains
                nlsingl     = nlsingl + 1
              end if
            end if
+            if (ilaggr(k) == ilaggr_neginit) then
+              write(0,*) iam,' Error: no update to ',k,mate(k),&
+                   & abs(w(k)),nrmagg,epsilon(nrmagg),wtemp(k)
+            end if
+          else
+            if (ilaggr(k)<0) write(0,*) 'Strange? ',k,ilaggr(k)
          end if
        end if
      end do
@@ -332,7 +355,7 @@ contains
      if (do_timings) call psb_tic(idx_phase3)

      ! Ok, now compute offsets, gather halo and fix non-local
-      ! aggregates  (those where ilaggr == -2)
+      ! aggregates  (those where ilaggr == ilaggr_nonlocal)
      call psb_sum(ictxt,nlaggr)
      ntaggr  = sum(nlaggr(0:np-1))
      naggrm1 = sum(nlaggr(0:iam-1))
@@ -347,7 +370,7 @@ contains
      call psb_halo(wtemp,desc_a,info)
      ! Cleanup as yet unmarked entries
      do k=1,nr
-        if (ilaggr(k) == -2) then
+        if (ilaggr(k) == ilaggr_nonlocal) then
          idx = mate(k)
          if (idx > nr) then
            i   = ilaggr(idx)
@@ -359,9 +382,14 @@ contains
          else
            write(0,*) 'Error : unresolved (paired) index ',k,idx,i,nr,nc, ilv(k),ilv(idx)
          end if
-        end if
-        if (ilaggr(k) <0) then
-          write(0,*) 'Matchboxp: Funny number: ',k,ilv(k),ilaggr(k),wtemp(k)
+        else if (ilaggr(k) <0) then
+          write(0,*) iam,'Matchboxp: Funny number: ',k,ilv(k),ilaggr(k),wtemp(k)
+          write(0,*) iam,'         :             : ',nr,nc,mate(k)
+          if (mate(k) <= nr) then
+            write(0,*) iam,'         :             : ',ilaggr(mate(k)),mate(mate(k)),&
+                 & ilv(k),ilv(mate(k)), ilv(mate(mate(k))),ilaggr(mate(mate(k)))
+          end if
+          flush(0)
        end if
      end do
      if (debug_sync) then
@@ -414,7 +442,7 @@ contains

      end block
      if (iam == 0) then
-        write(0,*) 'Matching statistics: Unmatched nodes ',&
+        write(0,*) iam,'Matching statistics: Unmatched nodes ',&
             & nunmatched,' Singletons:',nlsingl,' Pairs:',nlpairs
      end if

@@ -62,7 +62,22 @@ amg_s_parmatch_smth_bld.o \
 amg_s_parmatch_spmm_bld_inner.o

 MPCOBJS=MatchBoxPC.o \
-algoDistEdgeApproxDomEdgesLinearSearchMesgBndlSmallMateC.o
+sendBundledMessages.o \
+initialize.o \
+extractUChunk.o \
+isAlreadyMatched.o \
+findOwnerOfGhost.o \
+clean.o \
+computeCandidateMate.o \
+parallelComputeCandidateMateB.o \
+processMatchedVertices.o \
+processMatchedVerticesAndSendMessages.o \
+processCrossEdge.o \
+queueTransfer.o \
+processMessages.o \
+processExposedVertex.o \
+algoDistEdgeApproxDomEdgesLinearSearchMesgBndlSmallMateC.o \
+algoDistEdgeApproxDomEdgesLinearSearchMesgBndlSmallMateCMP.o

 OBJS = $(FOBJS) $(MPCOBJS)

@@ -60,17 +60,43 @@ void dMatchBoxPC(MilanLongInt NLVer, MilanLongInt NLEdge,
 		MilanLongInt* ph1_card, MilanLongInt* ph2_card ) {
 #if !defined(SERIAL_MPI)
  MPI_Comm C_comm=MPI_Comm_f2c(icomm);
+
 #ifdef DEBUG
  fprintf(stderr,"MatchBoxPC: rank %d nlver %ld nledge %ld [ %ld %ld ]\n",
 	  myRank,NLVer, NLEdge,verDistance[0],verDistance[1]);
 #endif
-  dalgoDistEdgeApproxDomEdgesLinearSearchMesgBndlSmallMateC(NLVer, NLEdge,
+
+
+#define TIME_TRACKER
+    #ifdef TIME_TRACKER
+        double tmr = MPI_Wtime();
+    #endif
+
+#define OMP
+#ifdef OMP
+        dalgoDistEdgeApproxDomEdgesLinearSearchMesgBndlSmallMateCMP(NLVer, NLEdge,
 							   verLocPtr, verLocInd, edgeLocWeight,
 							   verDistance,  Mate,
 							   myRank, numProcs, C_comm,
 							   msgIndSent, msgActualSent, msgPercent,
 							   ph0_time, ph1_time, ph2_time,
 							   ph1_card, ph2_card );
+#else
+        dalgoDistEdgeApproxDomEdgesLinearSearchMesgBndlSmallMateC(NLVer, NLEdge,
+							   verLocPtr, verLocInd, edgeLocWeight,
+							   verDistance,  Mate,
+							   myRank, numProcs, C_comm,
+							   msgIndSent, msgActualSent, msgPercent,
+							   ph0_time, ph1_time, ph2_time,
+							   ph1_card, ph2_card );
+#endif
+
+
+  #ifdef TIME_TRACKER
+    tmr = MPI_Wtime() - tmr;
+    fprintf(stderr, "Elaboration time: %f for %ld nodes\n", tmr, NLVer);
+  #endif
+
 #endif
 }

@@ -52,145 +52,412 @@

 #ifndef _matchboxpC_H_
 #define _matchboxpC_H_
-//Turn on a lot of debugging information with this switch:
+// Turn on a lot of debugging information with this switch:
 //#define PRINT_DEBUG_INFO_
 #include <stdio.h>
 #include <iostream>
 #include <assert.h>
 #include <map>
 #include <vector>
-// #include "matchboxp.h"
+#include "omp.h"
 #include "primitiveDataTypeDefinitions.h"
 #include "dataStrStaticQueue.h"

 using namespace std;

+const int NUM_THREAD = 4;
+const int UCHUNK = 10;
+
+const MilanLongInt REQUEST = 1;
+const MilanLongInt SUCCESS = 2;
+const MilanLongInt FAILURE = 3;
+const MilanLongInt SIZEINFO = 4;
+
+const int ComputeTag = 7; // Predefined tag
+const int BundleTag = 9;  // Predefined tag
+
+static vector<MilanLongInt> DEFAULT_VECTOR;
+
+// MPI type map
+template <typename T>
+MPI_Datatype TypeMap();
+template <>
+inline MPI_Datatype TypeMap<int64_t>() { return MPI_LONG_LONG; }
+template <>
+inline MPI_Datatype TypeMap<int>() { return MPI_INT; }
+template <>
+inline MPI_Datatype TypeMap<double>() { return MPI_DOUBLE; }
+template <>
+inline MPI_Datatype TypeMap<float>() { return MPI_FLOAT; }
+
 #ifdef __cplusplus
-extern "C" {
+extern "C"
+{
 #endif

 #if !defined(SERIAL_MPI)
-  
-#define MilanMpiLongInt  MPI_LONG_LONG
+
+#define MilanMpiLongInt MPI_LONG_LONG

 #ifndef _primitiveDataType_Definition_
 #define _primitiveDataType_Definition_
-    //Regular integer:
-    #ifndef INTEGER_H
-    #define INTEGER_H
-        typedef int32_t MilanInt;
-    #endif
+// Regular integer:
+#ifndef INTEGER_H
+#define INTEGER_H
+    typedef int32_t MilanInt;
+#endif

-    //Regular long integer:
-    #ifndef LONG_INT_H
-    #define LONG_INT_H
-        #ifdef BIT64
-            typedef int64_t MilanLongInt;
-            typedef MPI_LONG MilanMpiLongInt;
-        #else
-            typedef int32_t MilanLongInt;
-            typedef MPI_INT MilanMpiLongInt;
-        #endif
-    #endif
+// Regular long integer:
+#ifndef LONG_INT_H
+#define LONG_INT_H
+#ifdef BIT64
+    typedef int64_t MilanLongInt;
+    typedef MPI_LONG MilanMpiLongInt;
+#else
+    typedef int32_t MilanLongInt;
+    typedef MPI_INT MilanMpiLongInt;
+#endif
+#endif

-    //Regular boolean
-    #ifndef BOOL_H
-    #define BOOL_H
-        typedef bool MilanBool;
-    #endif
+// Regular boolean
+#ifndef BOOL_H
+#define BOOL_H
+    typedef bool MilanBool;
+#endif

-    //Regular double and absolute value computation:
-    #ifndef REAL_H
-    #define REAL_H
-        typedef double MilanReal;
-        typedef MPI_DOUBLE MilanMpiReal;
-        inline MilanReal MilanAbs(MilanReal value)
-        {
-            return fabs(value);
-        }
-    #endif
+// Regular double and absolute value computation:
+#ifndef REAL_H
+#define REAL_H
+    typedef double MilanReal;
+    typedef MPI_DOUBLE MilanMpiReal;
+    inline MilanReal MilanAbs(MilanReal value)
+    {
+        return fabs(value);
+    }
+#endif

-    //Regular float and absolute value computation:
-    #ifndef FLOAT_H
-    #define FLOAT_H
-        typedef float MilanFloat;
-        typedef MPI_FLOAT MilanMpiFloat;
-        inline MilanFloat MilanAbsFloat(MilanFloat value)
-        {
-            return fabs(value);
-        }
-    #endif
+// Regular float and absolute value computation:
+#ifndef FLOAT_H
+#define FLOAT_H
+    typedef float MilanFloat;
+    typedef MPI_FLOAT MilanMpiFloat;
+    inline MilanFloat MilanAbsFloat(MilanFloat value)
+    {
+        return fabs(value);
+    }
+#endif

-    //// Define the limits:
-    #ifndef LIMITS_H
-    #define LIMITS_H
-    //Integer Maximum and Minimum:
-  //      #define MilanIntMax INT_MAX
-  //    #define MilanIntMin INT_MIN
-        #define MilanIntMax INT32_MAX
-        #define MilanIntMin INT32_MIN
+//// Define the limits:
+#ifndef LIMITS_H
+#define LIMITS_H
+    // Integer Maximum and Minimum:
+    //      #define MilanIntMax INT_MAX
+    //    #define MilanIntMin INT_MIN
+#define MilanIntMax INT32_MAX
+#define MilanIntMin INT32_MIN

-        #ifdef BIT64
-            #define MilanLongIntMax INT64_MAX
-            #define MilanLongIntMin -INT64_MAX
-        #else
-            #define MilanLongIntMax INT32_MAX
-            #define MilanLongIntMin -INT32_MAX
-        #endif
+#ifdef BIT64
+#define MilanLongIntMax INT64_MAX
+#define MilanLongIntMin -INT64_MAX
+#else
+#define MilanLongIntMax INT32_MAX
+#define MilanLongIntMin -INT32_MAX
+#endif

-    #endif
+#endif

    // +INFINITY
    const double PLUS_INFINITY = numeric_limits<int>::infinity();
    const double MINUS_INFINITY = -PLUS_INFINITY;
-    //#define MilanRealMax LDBL_MAX
-    #define MilanRealMax PLUS_INFINITY
-    #define MilanRealMin MINUS_INFINITY
+//#define MilanRealMax LDBL_MAX
+#define MilanRealMax PLUS_INFINITY
+#define MilanRealMin MINUS_INFINITY
 #endif

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

-  void dalgoDistEdgeApproxDomEdgesLinearSearchMesgBndlSmallMateC
-(
- MilanLongInt NLVer, MilanLongInt NLEdge,
- MilanLongInt* verLocPtr, MilanLongInt* verLocInd, MilanReal* edgeLocWeight,
- MilanLongInt* verDistance,
- MilanLongInt* Mate,
- MilanInt myRank, MilanInt numProcs, MPI_Comm comm,
- MilanLongInt* msgIndSent, MilanLongInt* msgActualSent, MilanReal* msgPercent,
- MilanReal* ph0_time, MilanReal* ph1_time, MilanReal* ph2_time,
- MilanLongInt* ph1_card, MilanLongInt* ph2_card );
+    MilanLongInt firstComputeCandidateMate(MilanLongInt adj1,
+                                           MilanLongInt adj2,
+                                           MilanLongInt *verLocInd,
+                                           MilanReal *edgeLocWeight);

- void salgoDistEdgeApproxDomEdgesLinearSearchMesgBndlSmallMateC
-(
-MilanLongInt NLVer, MilanLongInt NLEdge,
-MilanLongInt* verLocPtr, MilanLongInt* verLocInd, MilanFloat* edgeLocWeight,
-MilanLongInt* verDistance,
-MilanLongInt* Mate,
-MilanInt myRank, MilanInt numProcs, MPI_Comm comm,
-MilanLongInt* msgIndSent, MilanLongInt* msgActualSent, MilanReal* msgPercent,
-MilanReal* ph0_time, MilanReal* ph1_time, MilanReal* ph2_time,
-MilanLongInt* ph1_card, MilanLongInt* ph2_card );
+    void queuesTransfer(vector<MilanLongInt> &U,
+                        vector<MilanLongInt> &privateU,
+                        vector<MilanLongInt> &QLocalVtx,
+                        vector<MilanLongInt> &QGhostVtx,
+                        vector<MilanLongInt> &QMsgType,
+                        vector<MilanInt> &QOwner,
+                        vector<MilanLongInt> &privateQLocalVtx,
+                        vector<MilanLongInt> &privateQGhostVtx,
+                        vector<MilanLongInt> &privateQMsgType,
+                        vector<MilanInt> &privateQOwner);

-void dMatchBoxPC(MilanLongInt NLVer, MilanLongInt NLEdge,
-		MilanLongInt* verLocPtr, MilanLongInt* verLocInd, MilanReal* edgeLocWeight,
-		MilanLongInt* verDistance,
-		MilanLongInt* Mate,
-		MilanInt myRank, MilanInt numProcs, MilanInt icomm,
-		MilanLongInt* msgIndSent, MilanLongInt* msgActualSent, MilanReal* msgPercent,
-		MilanReal* ph0_time, MilanReal* ph1_time, MilanReal* ph2_time,
-		MilanLongInt* ph1_card, MilanLongInt* ph2_card );
+    bool isAlreadyMatched(MilanLongInt node,
+                          MilanLongInt StartIndex,
+                          MilanLongInt EndIndex,
+                          vector<MilanLongInt> &GMate,
+                          MilanLongInt *Mate,
+                          map<MilanLongInt, MilanLongInt> &Ghost2LocalMap);

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

 #endif
 #ifdef __cplusplus
@@ -72,12 +72,6 @@

 #ifdef SERIAL_MPI
 #else
-//MPI type map
-template<typename T> MPI_Datatype TypeMap();
-template<> inline MPI_Datatype TypeMap<int64_t>() { return MPI_LONG_LONG; }
-template<> inline MPI_Datatype TypeMap<int>() { return MPI_INT; }
-template<> inline MPI_Datatype TypeMap<double>() { return MPI_DOUBLE; }
-template<> inline MPI_Datatype TypeMap<float>() { return MPI_FLOAT; }

 // DOUBLE PRECISION VERSION
 //WARNING: The vertex block on a given rank is contiguous
@@ -0,0 +1,554 @@
+#include "MatchBoxPC.h"
+
+// ***********************************************************************
+//
+//        MatchboxP: A C++ library for approximate weighted matching
+//               Mahantesh Halappanavar (hala@pnnl.gov)
+//               Pacific Northwest National Laboratory
+//
+// ***********************************************************************
+//
+//       Copyright (2021) Battelle Memorial Institute
+//                      All rights reserved.
+//
+// Redistribution and use in source and binary forms, with or without
+// modification, are permitted provided that the following conditions
+// are met:
+//
+// 1. Redistributions of source code must retain the above copyright
+// notice, this list of conditions and the following disclaimer.
+//
+// 2. Redistributions in binary form must reproduce the above copyright
+// notice, this list of conditions and the following disclaimer in the
+// documentation and/or other materials provided with the distribution.
+//
+// 3. Neither the name of the copyright holder nor the names of its
+// contributors may be used to endorse or promote products derived from
+// this software without specific prior written permission.
+//
+// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+// "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+// LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
+// FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
+// COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
+// INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
+// BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
+// LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
+// CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
+// LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
+// ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+// POSSIBILITY OF SUCH DAMAGE.
+//
+// ************************************************************************
+//////////////////////////////////////////////////////////////////////////////////////
+/////////////////////////// DOMINATING EDGES MODEL ///////////////////////////////////
+//////////////////////////////////////////////////////////////////////////////////////
+/* Function	: algoDistEdgeApproxDomEdgesLinearSearchMesgBndlSmallMate()
+ *
+ * Date     : New update: Feb 17, 2019, Richland, Washington.
+ * Date		: Original development: May 17, 2009, E&CS Bldg.
+ *
+ * Purpose	: Compute Approximate Maximum Weight Matching in Linear Time
+ *
+ * Args		: inputMatrix - instance of Compressed-Col format of Matrix
+ *                Mate - The Mate array
+ *
+ * Returns	: By Value: (void)
+ *            By Reference: Mate
+ *
+ * Comments	: 1/2 Approx Algorithm. Picks the locally available heaviest edge.
+ *                Assumption: The Mate Array is empty.
+ */
+
+/*
+ NLVer = #of vertices, NLEdge = #of edges
+ CSR/CSC/Compressed format: verLocPtr = Pointer, verLocInd = Index, edgeLocWeight = edge weights (positive real numbers)
+ verDistance = A vector of size |P|+1 containing the cumulative number of vertices per process
+ Mate = A vector of size |V_p| (local subgraph) to store the output (matching)
+ MPI: myRank, numProcs, comm,
+ Statistics: msgIndSent, msgActualSent, msgPercent : Size: |P| number of processes in the comm-world
+ Statistics: ph0_time, ph1_time, ph2_time: Runtimes
+ Statistics: ph1_card, ph2_card : Size: |P| number of processes in the comm-world (number of matched edges in Phase 1 and Phase 2)
+ */
+//#define DEBUG_HANG_
+#ifdef SERIAL_MPI
+#else
+
+// DOUBLE PRECISION VERSION
+// WARNING: The vertex block on a given rank is contiguous
+void dalgoDistEdgeApproxDomEdgesLinearSearchMesgBndlSmallMateCMP(
+    MilanLongInt NLVer, MilanLongInt NLEdge,
+    MilanLongInt *verLocPtr, MilanLongInt *verLocInd,
+    MilanReal *edgeLocWeight,
+    MilanLongInt *verDistance,
+    MilanLongInt *Mate,
+    MilanInt myRank, MilanInt numProcs, MPI_Comm comm,
+    MilanLongInt *msgIndSent, MilanLongInt *msgActualSent,
+    MilanReal *msgPercent,
+    MilanReal *ph0_time, MilanReal *ph1_time, MilanReal *ph2_time,
+    MilanLongInt *ph1_card, MilanLongInt *ph2_card)
+{
+
+    /*
+     * verDistance: it's a vector long as the number of processors.
+     *              verDistance[i] contains the first node index of the i-th processor
+     *              verDistance[i + 1] contains the last node index of the i-th processor
+     * NLVer: number of elements in the LocPtr
+     * NLEdge: number of edges assigned to the current processor
+     *
+     * Contains the portion of matrix assigned to the processor in
+     * Yale notation
+     * verLocInd: contains the positions on row of the matrix
+     * verLocPtr: i-th value is the position of the first element on the i-th row and
+     *            i+1-th value is the position of the first element on the i+1-th row
+     */
+
+#if !defined(SERIAL_MPI)
+#ifdef PRINT_DEBUG_INFO_
+    cout << "\n(" << myRank << ")Within algoEdgeApproxDominatingEdgesLinearSearchMessageBundling()";
+    fflush(stdout);
+#endif
+
+#ifdef PRINT_DEBUG_INFO_
+    cout << "\n(" << myRank << ") verDistance [" ;
+    for (int i = 0; i < numProcs; i++)
+      cout << verDistance[i] << "," << verDistance[i+1];
+    cout  << "]\n";
+    fflush(stdout);
+#endif
+#ifdef DEBUG_HANG_
+    if (myRank == 0) {
+      cout << "\n(" << myRank << ") verDistance [" ;
+      for (int i = 0; i < numProcs; i++)
+	cout << verDistance[i] << "," ;
+      cout  << verDistance[numProcs]<< "]\n";
+    }
+    fflush(stdout);
+#endif
+
+    MilanLongInt StartIndex = verDistance[myRank];       // The starting vertex owned by the current rank
+    MilanLongInt EndIndex   = verDistance[myRank + 1] - 1; // The ending vertex owned by the current rank
+
+    MPI_Status computeStatus;
+
+    MilanLongInt msgActual = 0, msgInd = 0;
+    MilanReal heaviestEdgeWt = 0.0f; // Assumes positive weight
+    MilanReal startTime, finishTime;
+
+    startTime = MPI_Wtime();
+
+    // Data structures for sending and receiving messages:
+    vector<MilanLongInt> Message; // [ u, v, message_type ]
+    Message.resize(3, -1);
+    // Data structures for Message Bundling:
+    // Although up to two messages can be sent along any cross edge,
+    // only one message will be sent in the initialization phase -
+    // one of: REQUEST/FAILURE/SUCCESS
+    vector<MilanLongInt> QLocalVtx, QGhostVtx, QMsgType;
+    vector<MilanInt> QOwner; // Changed by Fabio to be an integer, addresses needs to be integers!
+
+    MilanLongInt *PCounter = new MilanLongInt[numProcs];
+    for (int i = 0; i < numProcs; i++)
+        PCounter[i] = 0;
+
+    MilanLongInt NumMessagesBundled = 0;
+    // TODO when the last computational section will be refactored this could be eliminated
+    MilanInt ghostOwner = 0; // Changed by Fabio to be an integer, addresses needs to be integers!
+    MilanLongInt *candidateMate = nullptr;
+#ifdef PRINT_DEBUG_INFO_
+    cout << "\n(" << myRank << ")NV: " << NLVer << "  Edges: " << NLEdge;
+    fflush(stdout);
+    cout << "\n(" << myRank << ")StartIndex: " << StartIndex << "  EndIndex: " << EndIndex;
+    fflush(stdout);
+#endif
+    // Other Variables:
+    MilanLongInt u = -1, v = -1, w = -1, i = 0;
+    MilanLongInt k = -1, adj1 = -1, adj2 = -1;
+    MilanLongInt k1 = -1, adj11 = -1, adj12 = -1;
+    MilanLongInt myCard = 0;
+
+    // Build the Ghost Vertex Set: Vg
+    map<MilanLongInt, MilanLongInt> Ghost2LocalMap;       // Map each ghost vertex to a local vertex
+    vector<MilanLongInt> Counter;                         // Store the edge count for each ghost vertex
+    MilanLongInt numGhostVertices = 0, numGhostEdges = 0; // Number of Ghost vertices
+
+#ifdef PRINT_DEBUG_INFO_
+    cout << "\n(" << myRank << ")About to compute Ghost Vertices...";
+    fflush(stdout);
+#endif
+#ifdef DEBUG_HANG_
+    if (myRank == 0)
+        cout << "\n(" << myRank << ")About to compute Ghost Vertices...";
+    fflush(stdout);
+#endif
+
+    // Define Adjacency Lists for Ghost Vertices:
+    // cout<<"Building Ghost data structures ... \n\n";
+    vector<MilanLongInt> verGhostPtr, verGhostInd, tempCounter;
+    // Mate array for ghost vertices:
+    vector<MilanLongInt> GMate; // Proportional to the number of ghost vertices
+    MilanLongInt S;
+    MilanLongInt privateMyCard = 0;
+    vector<MilanLongInt> PCumulative, PMessageBundle, PSizeInfoMessages;
+    vector<MPI_Request> SRequest;  // Requests that are used for each send message
+    vector<MPI_Status> SStatus;    // Status of sent messages, used in MPI_Wait
+    MilanLongInt MessageIndex = 0; // Pointer for current message
+    MilanInt BufferSize;
+    MilanLongInt *Buffer;
+
+    vector<MilanLongInt> privateQLocalVtx, privateQGhostVtx, privateQMsgType;
+    vector<MilanInt> privateQOwner;
+    vector<MilanLongInt> U, privateU;
+
+    initialize(NLVer, NLEdge, StartIndex,
+               EndIndex, &numGhostEdges,
+               &numGhostVertices, &S,
+               verLocInd, verLocPtr,
+               Ghost2LocalMap, Counter,
+               verGhostPtr, verGhostInd,
+               tempCounter, GMate,
+               Message, QLocalVtx,
+               QGhostVtx, QMsgType, QOwner,
+               candidateMate, U,
+               privateU,
+               privateQLocalVtx,
+               privateQGhostVtx,
+               privateQMsgType,
+               privateQOwner);
+
+    finishTime = MPI_Wtime();
+    *ph0_time = finishTime - startTime; // Time taken for Phase-0: Initialization
+#ifdef DEBUG_HANG_
+    cout << myRank << " Finished initialization" << endl;
+    fflush(stdout);
+#endif
+   
+    startTime = MPI_Wtime();
+
+    /////////////////////////////////////////////////////////////////////////////////////////
+    //////////////////////////////////// INITIALIZATION /////////////////////////////////////
+    /////////////////////////////////////////////////////////////////////////////////////////
+    // Compute the Initial Matching Set:
+
+    /*
+     * OMP PARALLEL_COMPUTE_CANDIDATE_MATE_B has been splitted from
+     * PARALLEL_PROCESS_EXPOSED_VERTEX_B in order to better parallelize
+     * the two.
+     * PARALLEL_COMPUTE_CANDIDATE_MATE_B is now totally parallel.
+     */
+
+    PARALLEL_COMPUTE_CANDIDATE_MATE_B(NLVer,
+                                      verLocPtr,
+                                      verLocInd,
+                                      myRank,
+                                      edgeLocWeight,
+                                      candidateMate);
+
+#ifdef DEBUG_HANG_
+    cout << myRank << " Finished Exposed Vertex" << endl;
+    fflush(stdout);
+#if 0
+    cout << myRank << " candidateMate after parallelCompute " <<endl;
+    for (int i=0; i<NLVer; i++) {
+      cout << candidateMate[i] << " " ;
+    }
+    cout << endl;
+#endif
+#endif
+    /*
+     * PARALLEL_PROCESS_EXPOSED_VERTEX_B
+     * TODO: write comment
+     *
+     * TODO: Test when it's actually more efficient to execute this code
+     *       in parallel.
+     */
+    PARALLEL_PROCESS_EXPOSED_VERTEX_B(NLVer,
+                                      candidateMate,
+                                      verLocInd,
+                                      verLocPtr,
+                                      StartIndex,
+                                      EndIndex,
+                                      Mate,
+                                      GMate,
+                                      Ghost2LocalMap,
+                                      edgeLocWeight,
+                                      &myCard,
+                                      &msgInd,
+                                      &NumMessagesBundled,
+                                      &S,
+                                      verDistance,
+                                      PCounter,
+                                      Counter,
+                                      myRank,
+                                      numProcs,
+                                      U,
+                                      privateU,
+                                      QLocalVtx,
+                                      QGhostVtx,
+                                      QMsgType,
+                                      QOwner,
+                                      privateQLocalVtx,
+                                      privateQGhostVtx,
+                                      privateQMsgType,
+                                      privateQOwner);
+
+    tempCounter.clear(); // Do not need this any more
+
+#ifdef DEBUG_HANG_
+    cout << myRank << " Finished Exposed Vertex" << endl;
+    fflush(stdout);
+#if 0
+    cout << myRank << " Mate after Exposed Vertices " <<endl;
+    for (int i=0; i<NLVer; i++) {
+      cout << Mate[i] << " " ;
+    }
+    cout << endl;
+#endif
+#endif
+
+    ///////////////////////////////////////////////////////////////////////////////////
+    /////////////////////////// PROCESS MATCHED VERTICES //////////////////////////////
+    ///////////////////////////////////////////////////////////////////////////////////
+
+    // TODO what would be the optimal UCHUNK
+    vector<MilanLongInt> UChunkBeingProcessed;
+    UChunkBeingProcessed.reserve(UCHUNK);
+
+    processMatchedVertices(NLVer,
+                           UChunkBeingProcessed,
+                           U,
+                           privateU,
+                           StartIndex,
+                           EndIndex,
+                           &myCard,
+                           &msgInd,
+                           &NumMessagesBundled,
+                           &S,
+                           verLocPtr,
+                           verLocInd,
+                           verDistance,
+                           PCounter,
+                           Counter,
+                           myRank,
+                           numProcs,
+                           candidateMate,
+                           GMate,
+                           Mate,
+                           Ghost2LocalMap,
+                           edgeLocWeight,
+                           QLocalVtx,
+                           QGhostVtx,
+                           QMsgType,
+                           QOwner,
+                           privateQLocalVtx,
+                           privateQGhostVtx,
+                           privateQMsgType,
+                           privateQOwner);
+
+
+#ifdef DEBUG_HANG_
+    cout << myRank << " Finished Process Vertices" << endl;
+    fflush(stdout);
+#if 0
+    cout << myRank << " Mate after Matched Vertices " <<endl;
+    for (int i=0; i<NLVer; i++) {
+      cout << Mate[i] << " " ;
+    }
+    cout << endl;
+#endif
+#endif
+
+    /////////////////////////////////////////////////////////////////////////////////////////
+    ///////////////////////////// SEND BUNDLED MESSAGES /////////////////////////////////////
+    /////////////////////////////////////////////////////////////////////////////////////////
+
+    sendBundledMessages(&numGhostEdges,
+                        &BufferSize,
+                        Buffer,
+                        PCumulative,
+                        PMessageBundle,
+                        PSizeInfoMessages,
+                        PCounter,
+                        NumMessagesBundled,
+                        &msgActual,
+                        &MessageIndex,
+                        numProcs,
+                        myRank,
+                        comm,
+                        QLocalVtx,
+                        QGhostVtx,
+                        QMsgType,
+                        QOwner,
+                        SRequest,
+                        SStatus);
+
+    ///////////////////////// END OF SEND BUNDLED MESSAGES //////////////////////////////////
+
+    finishTime = MPI_Wtime();
+    *ph1_time = finishTime - startTime; // Time taken for Phase-1
+
+#ifdef DEBUG_HANG_
+    cout << myRank << " Finished sendBundles" << endl;
+    fflush(stdout);
+#endif
+ 
+    *ph1_card = myCard;                 // Cardinality at the end of Phase-1
+    startTime = MPI_Wtime();
+    /////////////////////////////////////////////////////////////////////////////////////////
+    //////////////////////////////////////// MAIN LOOP //////////////////////////////////////
+    /////////////////////////////////////////////////////////////////////////////////////////
+    // Main While Loop:
+#ifdef PRINT_DEBUG_INFO_
+    cout << "\n(" << myRank << "=========================************===============================" << endl;
+    fflush(stdout);
+    fflush(stdout);
+#endif
+#ifdef PRINT_DEBUG_INFO_
+    cout << "\n(" << myRank << ")Entering While(true) loop..";
+    fflush(stdout);
+#endif
+#ifdef PRINT_DEBUG_INFO_
+    cout << "\n(" << myRank << "=========================************===============================" << endl;
+    fflush(stdout);
+    fflush(stdout);
+#endif
+
+    while (true)  {
+#ifdef DEBUG_HANG_
+      //if (myRank == 0)
+      cout << "\n(" << myRank << ") Main loop" << endl;
+      fflush(stdout);
+#endif
+        ///////////////////////////////////////////////////////////////////////////////////
+        /////////////////////////// PROCESS MATCHED VERTICES //////////////////////////////
+        ///////////////////////////////////////////////////////////////////////////////////
+      
+      processMatchedVerticesAndSendMessages(NLVer,
+					    UChunkBeingProcessed,
+					    U,
+					    privateU,
+					    StartIndex,
+					    EndIndex,
+					    &myCard,
+					    &msgInd,
+					    &NumMessagesBundled,
+					    &S,
+					    verLocPtr,
+					    verLocInd,
+					    verDistance,
+					    PCounter,
+					    Counter,
+					    myRank,
+					    numProcs,
+					    candidateMate,
+					    GMate,
+					    Mate,
+					    Ghost2LocalMap,
+					    edgeLocWeight,
+					    QLocalVtx,
+					    QGhostVtx,
+					    QMsgType,
+					    QOwner,
+					    privateQLocalVtx,
+					    privateQGhostVtx,
+					    privateQMsgType,
+					    privateQOwner,
+					    comm,
+					    &msgActual,
+					    Message);
+      
+      ///////////////////////// END OF PROCESS MATCHED VERTICES /////////////////////////
+
+        //// BREAK IF NO MESSAGES EXPECTED /////////
+#ifdef DEBUG_HANG_
+#if 0
+      cout << myRank << " Mate after ProcessMatchedAndSend phase "<<S <<endl;
+      for (int i=0; i<NLVer; i++) {
+	cout << Mate[i] << " " ;
+      }
+      cout << endl;
+#endif
+#endif
+#ifdef PRINT_DEBUG_INFO_
+        cout << "\n(" << myRank << ")Deciding whether to break: S= " << S << endl;
+#endif
+
+        if (S == 0) {
+#ifdef DEBUG_HANG_
+	  cout << "\n(" << myRank << ") Breaking out" << endl;
+	  fflush(stdout);
+#endif
+	  break;
+        }
+        ///////////////////////////////////////////////////////////////////////////////////
+        /////////////////////////// PROCESS MESSAGES //////////////////////////////////////
+        ///////////////////////////////////////////////////////////////////////////////////
+
+        processMessages(NLVer,
+                        Mate,
+                        candidateMate,
+                        Ghost2LocalMap,
+                        GMate,
+                        Counter,
+                        StartIndex,
+                        EndIndex,
+                        &myCard,
+                        &msgInd,
+                        &msgActual,
+                        edgeLocWeight,
+                        verDistance,
+                        verLocPtr,
+                        k,
+                        verLocInd,
+                        numProcs,
+                        myRank,
+                        comm,
+                        Message,
+                        numGhostEdges,
+                        u,
+                        v,
+                        &S,
+                        U);
+
+        ///////////////////////// END OF PROCESS MESSAGES /////////////////////////////////
+#ifdef DEBUG_HANG_
+#if 0
+      cout << myRank << " Mate after ProcessMessages phase "<<S <<endl;
+      for (int i=0; i<NLVer; i++) {
+	cout << Mate[i] << " " ;
+      }
+      cout << endl;
+#endif
+#endif
+#ifdef PRINT_DEBUG_INFO_
+        cout << "\n(" << myRank << ")Finished Message processing phase: S= " << S;
+        fflush(stdout);
+        cout << "\n(" << myRank << ")** SENT     : ACTUAL= " << msgActual;
+        fflush(stdout);
+        cout << "\n(" << myRank << ")** SENT     : INDIVIDUAL= " << msgInd << endl;
+        fflush(stdout);
+#endif
+    } // End of while (true)
+
+    clean(NLVer,
+          myRank,
+          MessageIndex,
+          SRequest,
+          SStatus,
+          BufferSize,
+          Buffer,
+          msgActual,
+          msgActualSent,
+          msgInd,
+          msgIndSent,
+          NumMessagesBundled,
+          msgPercent);
+
+    finishTime = MPI_Wtime();
+    *ph2_time = finishTime - startTime; // Time taken for Phase-2
+    *ph2_card = myCard;                 // Cardinality at the end of Phase-2
+}
+// End of algoDistEdgeApproxDomEdgesLinearSearchMesgBndlSmallMate
+#endif
+
+#endif
@@ -97,6 +97,8 @@ subroutine  amg_c_dec_aggregator_build_tprol(ag,parms,ag_data,&
  integer(psb_lpk_)   :: ntaggr
  integer(psb_ipk_)   :: debug_level, debug_unit
  logical             :: clean_zeros
+  integer(psb_ipk_), save :: idx_map_bld=-1, idx_map_tprol=-1
+  logical, parameter      :: do_timings=.false.

  name='amg_c_dec_aggregator_tprol'
  call psb_erractionsave(err_act)
@@ -108,6 +110,10 @@ subroutine  amg_c_dec_aggregator_build_tprol(ag,parms,ag_data,&
  info  = psb_success_
  ctxt = desc_a%get_context()
  call psb_info(ctxt,me,np)
+  if ((do_timings).and.(idx_map_bld==-1))       &
+       & idx_map_bld = psb_get_timer_idx("DEC_TPROL: map_bld")
+  if ((do_timings).and.(idx_map_tprol==-1))     &
+       & idx_map_tprol = psb_get_timer_idx("DEC_TPROL: map_tprol")

  call amg_check_def(parms%ml_cycle,'Multilevel cycle',&
       &   amg_mult_ml_,is_legal_ml_cycle)
@@ -121,10 +127,14 @@ subroutine  amg_c_dec_aggregator_build_tprol(ag,parms,ag_data,&
  ! The decoupled aggregator based on SOC measures ignores
  ! ag_data except for clean_zeros; soc_map_bld is a procedure pointer.
  !
+  if (do_timings) call psb_tic(idx_map_bld)
  clean_zeros = ag%do_clean_zeros
  call ag%soc_map_bld(parms%aggr_ord,parms%aggr_thresh,clean_zeros,a,desc_a,nlaggr,ilaggr,info)
+  if (do_timings) call psb_toc(idx_map_bld)
+  if (do_timings) call psb_tic(idx_map_tprol)

  if (info==psb_success_) call amg_map_to_tprol(desc_a,ilaggr,nlaggr,t_prol,info)
+  if (do_timings) call psb_toc(idx_map_tprol)
  if (info /= psb_success_) then
    info=psb_err_from_subroutine_
    call psb_errpush(info,name,a_err='soc_map_bld/map_to_tprol')
@@ -140,6 +140,9 @@ subroutine amg_caggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,&
  real(psb_spk_)     :: anorm, omega, tmp, dg, theta
  logical, parameter :: debug_new=.false.
  character(len=80) :: filename
+  logical, parameter :: do_timings=.false.
+  integer(psb_ipk_), save :: idx_spspmm=-1, idx_phase1=-1, idx_gtrans=-1, idx_phase2=-1, idx_refine=-1
+  integer(psb_ipk_), save :: idx_phase3=-1, idx_cdasb=-1, idx_ptap=-1

  name='amg_aggrmat_smth_bld'
  info=psb_success_
@@ -153,6 +156,23 @@ subroutine amg_caggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,&
  ctxt = desc_a%get_context()

  call psb_info(ctxt, me, np)
+  if ((do_timings).and.(idx_spspmm==-1)) &
+       & idx_spspmm = psb_get_timer_idx("DEC_SMTH_BLD: par_spspmm")
+  if ((do_timings).and.(idx_phase1==-1)) &
+       & idx_phase1 = psb_get_timer_idx("DEC_SMTH_BLD: phase1    ")
+  if ((do_timings).and.(idx_phase2==-1)) &
+       & idx_phase2 = psb_get_timer_idx("DEC_SMTH_BLD: phase2    ")
+  if ((do_timings).and.(idx_phase3==-1)) &
+       & idx_phase3 = psb_get_timer_idx("DEC_SMTH_BLD: phase3    ")
+  if ((do_timings).and.(idx_gtrans==-1)) &
+       & idx_gtrans = psb_get_timer_idx("DEC_SMTH_BLD: gtrans    ")
+  if ((do_timings).and.(idx_refine==-1)) &
+       & idx_refine = psb_get_timer_idx("DEC_SMTH_BLD: refine    ")
+  if ((do_timings).and.(idx_cdasb==-1)) &
+       & idx_cdasb = psb_get_timer_idx("DEC_SMTH_BLD: cdasb     ")
+  if ((do_timings).and.(idx_ptap==-1)) &
+       & idx_ptap = psb_get_timer_idx("DEC_SMTH_BLD: ptap_bld  ")
+

  nglob = desc_a%get_global_rows()
  nrow  = desc_a%get_local_rows()
@@ -171,6 +191,7 @@ subroutine amg_caggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,&
  ! naggr: number of local aggregates
  ! nrow: local rows. 
  ! 
+  if (do_timings) call psb_tic(idx_phase1)

  ! Get the diagonal D
  adiag = a%get_diag(info)
@@ -196,7 +217,7 @@ subroutine amg_caggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,&
    !
    ! Build the filtered matrix Af from A
    ! 
-
+    !$OMP parallel do private(i,j,tmp,jd) schedule(static)
    do i=1, nrow
      tmp = czero
      jd  = -1 
@@ -214,11 +235,13 @@ subroutine amg_caggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,&
        acsrf%val(jd)=acsrf%val(jd)-tmp
      end if
    enddo
+    !$OMP end parallel do 
    ! Take out zeroed terms 
    call acsrf%clean_zeros(info)
  end if


+  !$OMP parallel do private(i) schedule(static)
  do i=1,size(adiag)
    if (adiag(i) /= czero) then
      adiag(i) = cone / adiag(i)
@@ -226,7 +249,7 @@ subroutine amg_caggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,&
      adiag(i) = cone
    end if
  end do
-
+  !$OMP end parallel do 
  if (parms%aggr_omega_alg == amg_eig_est_) then 

    if (parms%aggr_eig == amg_max_norm_) then 
@@ -252,8 +275,9 @@ subroutine amg_caggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,&
    call psb_errpush(info,name,a_err='invalid amg_aggr_omega_alg_')
    goto 9999
  end if
+  if (do_timings) call psb_toc(idx_phase1)

-  
+  if (do_timings) call psb_tic(idx_phase2)
  call acsrf%scal(adiag,info)
  if (info /= psb_success_) goto 9999

@@ -267,6 +291,8 @@ subroutine amg_caggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,&

  call psb_cdasb(desc_ac,info)
  call psb_cd_reinit(desc_ac,info)
+  if (do_timings) call psb_toc(idx_phase2)
+  if (do_timings) call psb_tic(idx_phase3)
  !
  ! Build the smoothed prolongator using either A or Af
  !    acsr1 = (I-w*D*A) Prol      acsr1 = (I-w*D*Af) Prol 
@@ -279,8 +305,8 @@ subroutine amg_caggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,&
    call psb_errpush(psb_err_from_subroutine_,name,a_err='spspmm 1')
    goto 9999
  end if
-
-
+  if (do_timings) call psb_toc(idx_phase3)
+  if (do_timings) call psb_tic(idx_ptap)
  if (debug_level >= psb_debug_outer_) &
       & write(debug_unit,*) me,' ',trim(name),&
       & 'Done SPSPMM 1'
@@ -292,7 +318,7 @@ subroutine amg_caggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,&

  call op_prol%mv_from(coo_prol)
  call op_restr%mv_from(coo_restr)
-
+  if (do_timings) call psb_toc(idx_ptap)
  if (debug_level >= psb_debug_outer_) &
       & write(debug_unit,*) me,' ',trim(name),&
       & 'Done smooth_aggregate '
@@ -97,6 +97,8 @@ subroutine  amg_d_dec_aggregator_build_tprol(ag,parms,ag_data,&
  integer(psb_lpk_)   :: ntaggr
  integer(psb_ipk_)   :: debug_level, debug_unit
  logical             :: clean_zeros
+  integer(psb_ipk_), save :: idx_map_bld=-1, idx_map_tprol=-1
+  logical, parameter      :: do_timings=.false.

  name='amg_d_dec_aggregator_tprol'
  call psb_erractionsave(err_act)
@@ -108,6 +110,10 @@ subroutine  amg_d_dec_aggregator_build_tprol(ag,parms,ag_data,&
  info  = psb_success_
  ctxt = desc_a%get_context()
  call psb_info(ctxt,me,np)
+  if ((do_timings).and.(idx_map_bld==-1))       &
+       & idx_map_bld = psb_get_timer_idx("DEC_TPROL: map_bld")
+  if ((do_timings).and.(idx_map_tprol==-1))     &
+       & idx_map_tprol = psb_get_timer_idx("DEC_TPROL: map_tprol")

  call amg_check_def(parms%ml_cycle,'Multilevel cycle',&
       &   amg_mult_ml_,is_legal_ml_cycle)
@@ -121,10 +127,14 @@ subroutine  amg_d_dec_aggregator_build_tprol(ag,parms,ag_data,&
  ! The decoupled aggregator based on SOC measures ignores
  ! ag_data except for clean_zeros; soc_map_bld is a procedure pointer.
  !
+  if (do_timings) call psb_tic(idx_map_bld)
  clean_zeros = ag%do_clean_zeros
  call ag%soc_map_bld(parms%aggr_ord,parms%aggr_thresh,clean_zeros,a,desc_a,nlaggr,ilaggr,info)
+  if (do_timings) call psb_toc(idx_map_bld)
+  if (do_timings) call psb_tic(idx_map_tprol)

  if (info==psb_success_) call amg_map_to_tprol(desc_a,ilaggr,nlaggr,t_prol,info)
+  if (do_timings) call psb_toc(idx_map_tprol)
  if (info /= psb_success_) then
    info=psb_err_from_subroutine_
    call psb_errpush(info,name,a_err='soc_map_bld/map_to_tprol')
@@ -140,6 +140,9 @@ subroutine amg_daggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,&
  real(psb_dpk_)     :: anorm, omega, tmp, dg, theta
  logical, parameter :: debug_new=.false.
  character(len=80) :: filename
+  logical, parameter :: do_timings=.false.
+  integer(psb_ipk_), save :: idx_spspmm=-1, idx_phase1=-1, idx_gtrans=-1, idx_phase2=-1, idx_refine=-1
+  integer(psb_ipk_), save :: idx_phase3=-1, idx_cdasb=-1, idx_ptap=-1

  name='amg_aggrmat_smth_bld'
  info=psb_success_
@@ -153,6 +156,23 @@ subroutine amg_daggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,&
  ctxt = desc_a%get_context()

  call psb_info(ctxt, me, np)
+  if ((do_timings).and.(idx_spspmm==-1)) &
+       & idx_spspmm = psb_get_timer_idx("DEC_SMTH_BLD: par_spspmm")
+  if ((do_timings).and.(idx_phase1==-1)) &
+       & idx_phase1 = psb_get_timer_idx("DEC_SMTH_BLD: phase1    ")
+  if ((do_timings).and.(idx_phase2==-1)) &
+       & idx_phase2 = psb_get_timer_idx("DEC_SMTH_BLD: phase2    ")
+  if ((do_timings).and.(idx_phase3==-1)) &
+       & idx_phase3 = psb_get_timer_idx("DEC_SMTH_BLD: phase3    ")
+  if ((do_timings).and.(idx_gtrans==-1)) &
+       & idx_gtrans = psb_get_timer_idx("DEC_SMTH_BLD: gtrans    ")
+  if ((do_timings).and.(idx_refine==-1)) &
+       & idx_refine = psb_get_timer_idx("DEC_SMTH_BLD: refine    ")
+  if ((do_timings).and.(idx_cdasb==-1)) &
+       & idx_cdasb = psb_get_timer_idx("DEC_SMTH_BLD: cdasb     ")
+  if ((do_timings).and.(idx_ptap==-1)) &
+       & idx_ptap = psb_get_timer_idx("DEC_SMTH_BLD: ptap_bld  ")
+

  nglob = desc_a%get_global_rows()
  nrow  = desc_a%get_local_rows()
@@ -171,6 +191,7 @@ subroutine amg_daggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,&
  ! naggr: number of local aggregates
  ! nrow: local rows. 
  ! 
+  if (do_timings) call psb_tic(idx_phase1)

  ! Get the diagonal D
  adiag = a%get_diag(info)
@@ -196,7 +217,7 @@ subroutine amg_daggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,&
    !
    ! Build the filtered matrix Af from A
    ! 
-
+    !$OMP parallel do private(i,j,tmp,jd) schedule(static)
    do i=1, nrow
      tmp = dzero
      jd  = -1 
@@ -214,11 +235,13 @@ subroutine amg_daggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,&
        acsrf%val(jd)=acsrf%val(jd)-tmp
      end if
    enddo
+    !$OMP end parallel do 
    ! Take out zeroed terms 
    call acsrf%clean_zeros(info)
  end if


+  !$OMP parallel do private(i) schedule(static)
  do i=1,size(adiag)
    if (adiag(i) /= dzero) then
      adiag(i) = done / adiag(i)
@@ -226,7 +249,7 @@ subroutine amg_daggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,&
      adiag(i) = done
    end if
  end do
-
+  !$OMP end parallel do 
  if (parms%aggr_omega_alg == amg_eig_est_) then 

    if (parms%aggr_eig == amg_max_norm_) then 
@@ -252,8 +275,9 @@ subroutine amg_daggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,&
    call psb_errpush(info,name,a_err='invalid amg_aggr_omega_alg_')
    goto 9999
  end if
+  if (do_timings) call psb_toc(idx_phase1)

-  
+  if (do_timings) call psb_tic(idx_phase2)
  call acsrf%scal(adiag,info)
  if (info /= psb_success_) goto 9999

@@ -267,6 +291,8 @@ subroutine amg_daggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,&

  call psb_cdasb(desc_ac,info)
  call psb_cd_reinit(desc_ac,info)
+  if (do_timings) call psb_toc(idx_phase2)
+  if (do_timings) call psb_tic(idx_phase3)
  !
  ! Build the smoothed prolongator using either A or Af
  !    acsr1 = (I-w*D*A) Prol      acsr1 = (I-w*D*Af) Prol 
@@ -279,8 +305,8 @@ subroutine amg_daggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,&
    call psb_errpush(psb_err_from_subroutine_,name,a_err='spspmm 1')
    goto 9999
  end if
-
-
+  if (do_timings) call psb_toc(idx_phase3)
+  if (do_timings) call psb_tic(idx_ptap)
  if (debug_level >= psb_debug_outer_) &
       & write(debug_unit,*) me,' ',trim(name),&
       & 'Done SPSPMM 1'
@@ -292,7 +318,7 @@ subroutine amg_daggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,&

  call op_prol%mv_from(coo_prol)
  call op_restr%mv_from(coo_restr)
-
+  if (do_timings) call psb_toc(idx_ptap)
  if (debug_level >= psb_debug_outer_) &
       & write(debug_unit,*) me,' ',trim(name),&
       & 'Done smooth_aggregate '
@@ -97,6 +97,8 @@ subroutine  amg_s_dec_aggregator_build_tprol(ag,parms,ag_data,&
  integer(psb_lpk_)   :: ntaggr
  integer(psb_ipk_)   :: debug_level, debug_unit
  logical             :: clean_zeros
+  integer(psb_ipk_), save :: idx_map_bld=-1, idx_map_tprol=-1
+  logical, parameter      :: do_timings=.false.

  name='amg_s_dec_aggregator_tprol'
  call psb_erractionsave(err_act)
@@ -108,6 +110,10 @@ subroutine  amg_s_dec_aggregator_build_tprol(ag,parms,ag_data,&
  info  = psb_success_
  ctxt = desc_a%get_context()
  call psb_info(ctxt,me,np)
+  if ((do_timings).and.(idx_map_bld==-1))       &
+       & idx_map_bld = psb_get_timer_idx("DEC_TPROL: map_bld")
+  if ((do_timings).and.(idx_map_tprol==-1))     &
+       & idx_map_tprol = psb_get_timer_idx("DEC_TPROL: map_tprol")

  call amg_check_def(parms%ml_cycle,'Multilevel cycle',&
       &   amg_mult_ml_,is_legal_ml_cycle)
@@ -121,10 +127,14 @@ subroutine  amg_s_dec_aggregator_build_tprol(ag,parms,ag_data,&
  ! The decoupled aggregator based on SOC measures ignores
  ! ag_data except for clean_zeros; soc_map_bld is a procedure pointer.
  !
+  if (do_timings) call psb_tic(idx_map_bld)
  clean_zeros = ag%do_clean_zeros
  call ag%soc_map_bld(parms%aggr_ord,parms%aggr_thresh,clean_zeros,a,desc_a,nlaggr,ilaggr,info)
+  if (do_timings) call psb_toc(idx_map_bld)
+  if (do_timings) call psb_tic(idx_map_tprol)

  if (info==psb_success_) call amg_map_to_tprol(desc_a,ilaggr,nlaggr,t_prol,info)
+  if (do_timings) call psb_toc(idx_map_tprol)
  if (info /= psb_success_) then
    info=psb_err_from_subroutine_
    call psb_errpush(info,name,a_err='soc_map_bld/map_to_tprol')
@@ -140,6 +140,9 @@ subroutine amg_saggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,&
  real(psb_spk_)     :: anorm, omega, tmp, dg, theta
  logical, parameter :: debug_new=.false.
  character(len=80) :: filename
+  logical, parameter :: do_timings=.false.
+  integer(psb_ipk_), save :: idx_spspmm=-1, idx_phase1=-1, idx_gtrans=-1, idx_phase2=-1, idx_refine=-1
+  integer(psb_ipk_), save :: idx_phase3=-1, idx_cdasb=-1, idx_ptap=-1

  name='amg_aggrmat_smth_bld'
  info=psb_success_
@@ -153,6 +156,23 @@ subroutine amg_saggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,&
  ctxt = desc_a%get_context()

  call psb_info(ctxt, me, np)
+  if ((do_timings).and.(idx_spspmm==-1)) &
+       & idx_spspmm = psb_get_timer_idx("DEC_SMTH_BLD: par_spspmm")
+  if ((do_timings).and.(idx_phase1==-1)) &
+       & idx_phase1 = psb_get_timer_idx("DEC_SMTH_BLD: phase1    ")
+  if ((do_timings).and.(idx_phase2==-1)) &
+       & idx_phase2 = psb_get_timer_idx("DEC_SMTH_BLD: phase2    ")
+  if ((do_timings).and.(idx_phase3==-1)) &
+       & idx_phase3 = psb_get_timer_idx("DEC_SMTH_BLD: phase3    ")
+  if ((do_timings).and.(idx_gtrans==-1)) &
+       & idx_gtrans = psb_get_timer_idx("DEC_SMTH_BLD: gtrans    ")
+  if ((do_timings).and.(idx_refine==-1)) &
+       & idx_refine = psb_get_timer_idx("DEC_SMTH_BLD: refine    ")
+  if ((do_timings).and.(idx_cdasb==-1)) &
+       & idx_cdasb = psb_get_timer_idx("DEC_SMTH_BLD: cdasb     ")
+  if ((do_timings).and.(idx_ptap==-1)) &
+       & idx_ptap = psb_get_timer_idx("DEC_SMTH_BLD: ptap_bld  ")
+

  nglob = desc_a%get_global_rows()
  nrow  = desc_a%get_local_rows()
@@ -171,6 +191,7 @@ subroutine amg_saggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,&
  ! naggr: number of local aggregates
  ! nrow: local rows. 
  ! 
+  if (do_timings) call psb_tic(idx_phase1)

  ! Get the diagonal D
  adiag = a%get_diag(info)
@@ -196,7 +217,7 @@ subroutine amg_saggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,&
    !
    ! Build the filtered matrix Af from A
    ! 
-
+    !$OMP parallel do private(i,j,tmp,jd) schedule(static)
    do i=1, nrow
      tmp = szero
      jd  = -1 
@@ -214,11 +235,13 @@ subroutine amg_saggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,&
        acsrf%val(jd)=acsrf%val(jd)-tmp
      end if
    enddo
+    !$OMP end parallel do 
    ! Take out zeroed terms 
    call acsrf%clean_zeros(info)
  end if


+  !$OMP parallel do private(i) schedule(static)
  do i=1,size(adiag)
    if (adiag(i) /= szero) then
      adiag(i) = sone / adiag(i)
@@ -226,7 +249,7 @@ subroutine amg_saggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,&
      adiag(i) = sone
    end if
  end do
-
+  !$OMP end parallel do 
  if (parms%aggr_omega_alg == amg_eig_est_) then 

    if (parms%aggr_eig == amg_max_norm_) then 
@@ -252,8 +275,9 @@ subroutine amg_saggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,&
    call psb_errpush(info,name,a_err='invalid amg_aggr_omega_alg_')
    goto 9999
  end if
+  if (do_timings) call psb_toc(idx_phase1)

-  
+  if (do_timings) call psb_tic(idx_phase2)
  call acsrf%scal(adiag,info)
  if (info /= psb_success_) goto 9999

@@ -267,6 +291,8 @@ subroutine amg_saggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,&

  call psb_cdasb(desc_ac,info)
  call psb_cd_reinit(desc_ac,info)
+  if (do_timings) call psb_toc(idx_phase2)
+  if (do_timings) call psb_tic(idx_phase3)
  !
  ! Build the smoothed prolongator using either A or Af
  !    acsr1 = (I-w*D*A) Prol      acsr1 = (I-w*D*Af) Prol 
@@ -279,8 +305,8 @@ subroutine amg_saggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,&
    call psb_errpush(psb_err_from_subroutine_,name,a_err='spspmm 1')
    goto 9999
  end if
-
-
+  if (do_timings) call psb_toc(idx_phase3)
+  if (do_timings) call psb_tic(idx_ptap)
  if (debug_level >= psb_debug_outer_) &
       & write(debug_unit,*) me,' ',trim(name),&
       & 'Done SPSPMM 1'
@@ -292,7 +318,7 @@ subroutine amg_saggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,&

  call op_prol%mv_from(coo_prol)
  call op_restr%mv_from(coo_restr)
-
+  if (do_timings) call psb_toc(idx_ptap)
  if (debug_level >= psb_debug_outer_) &
       & write(debug_unit,*) me,' ',trim(name),&
       & 'Done smooth_aggregate '
@@ -97,6 +97,8 @@ subroutine  amg_z_dec_aggregator_build_tprol(ag,parms,ag_data,&
  integer(psb_lpk_)   :: ntaggr
  integer(psb_ipk_)   :: debug_level, debug_unit
  logical             :: clean_zeros
+  integer(psb_ipk_), save :: idx_map_bld=-1, idx_map_tprol=-1
+  logical, parameter      :: do_timings=.false.

  name='amg_z_dec_aggregator_tprol'
  call psb_erractionsave(err_act)
@@ -108,6 +110,10 @@ subroutine  amg_z_dec_aggregator_build_tprol(ag,parms,ag_data,&
  info  = psb_success_
  ctxt = desc_a%get_context()
  call psb_info(ctxt,me,np)
+  if ((do_timings).and.(idx_map_bld==-1))       &
+       & idx_map_bld = psb_get_timer_idx("DEC_TPROL: map_bld")
+  if ((do_timings).and.(idx_map_tprol==-1))     &
+       & idx_map_tprol = psb_get_timer_idx("DEC_TPROL: map_tprol")

  call amg_check_def(parms%ml_cycle,'Multilevel cycle',&
       &   amg_mult_ml_,is_legal_ml_cycle)
@@ -121,10 +127,14 @@ subroutine  amg_z_dec_aggregator_build_tprol(ag,parms,ag_data,&
  ! The decoupled aggregator based on SOC measures ignores
  ! ag_data except for clean_zeros; soc_map_bld is a procedure pointer.
  !
+  if (do_timings) call psb_tic(idx_map_bld)
  clean_zeros = ag%do_clean_zeros
  call ag%soc_map_bld(parms%aggr_ord,parms%aggr_thresh,clean_zeros,a,desc_a,nlaggr,ilaggr,info)
+  if (do_timings) call psb_toc(idx_map_bld)
+  if (do_timings) call psb_tic(idx_map_tprol)

  if (info==psb_success_) call amg_map_to_tprol(desc_a,ilaggr,nlaggr,t_prol,info)
+  if (do_timings) call psb_toc(idx_map_tprol)
  if (info /= psb_success_) then
    info=psb_err_from_subroutine_
    call psb_errpush(info,name,a_err='soc_map_bld/map_to_tprol')
@@ -140,6 +140,9 @@ subroutine amg_zaggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,&
  real(psb_dpk_)     :: anorm, omega, tmp, dg, theta
  logical, parameter :: debug_new=.false.
  character(len=80) :: filename
+  logical, parameter :: do_timings=.false.
+  integer(psb_ipk_), save :: idx_spspmm=-1, idx_phase1=-1, idx_gtrans=-1, idx_phase2=-1, idx_refine=-1
+  integer(psb_ipk_), save :: idx_phase3=-1, idx_cdasb=-1, idx_ptap=-1

  name='amg_aggrmat_smth_bld'
  info=psb_success_
@@ -153,6 +156,23 @@ subroutine amg_zaggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,&
  ctxt = desc_a%get_context()

  call psb_info(ctxt, me, np)
+  if ((do_timings).and.(idx_spspmm==-1)) &
+       & idx_spspmm = psb_get_timer_idx("DEC_SMTH_BLD: par_spspmm")
+  if ((do_timings).and.(idx_phase1==-1)) &
+       & idx_phase1 = psb_get_timer_idx("DEC_SMTH_BLD: phase1    ")
+  if ((do_timings).and.(idx_phase2==-1)) &
+       & idx_phase2 = psb_get_timer_idx("DEC_SMTH_BLD: phase2    ")
+  if ((do_timings).and.(idx_phase3==-1)) &
+       & idx_phase3 = psb_get_timer_idx("DEC_SMTH_BLD: phase3    ")
+  if ((do_timings).and.(idx_gtrans==-1)) &
+       & idx_gtrans = psb_get_timer_idx("DEC_SMTH_BLD: gtrans    ")
+  if ((do_timings).and.(idx_refine==-1)) &
+       & idx_refine = psb_get_timer_idx("DEC_SMTH_BLD: refine    ")
+  if ((do_timings).and.(idx_cdasb==-1)) &
+       & idx_cdasb = psb_get_timer_idx("DEC_SMTH_BLD: cdasb     ")
+  if ((do_timings).and.(idx_ptap==-1)) &
+       & idx_ptap = psb_get_timer_idx("DEC_SMTH_BLD: ptap_bld  ")
+

  nglob = desc_a%get_global_rows()
  nrow  = desc_a%get_local_rows()
@@ -171,6 +191,7 @@ subroutine amg_zaggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,&
  ! naggr: number of local aggregates
  ! nrow: local rows. 
  ! 
+  if (do_timings) call psb_tic(idx_phase1)

  ! Get the diagonal D
  adiag = a%get_diag(info)
@@ -196,7 +217,7 @@ subroutine amg_zaggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,&
    !
    ! Build the filtered matrix Af from A
    ! 
-
+    !$OMP parallel do private(i,j,tmp,jd) schedule(static)
    do i=1, nrow
      tmp = zzero
      jd  = -1 
@@ -214,11 +235,13 @@ subroutine amg_zaggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,&
        acsrf%val(jd)=acsrf%val(jd)-tmp
      end if
    enddo
+    !$OMP end parallel do 
    ! Take out zeroed terms 
    call acsrf%clean_zeros(info)
  end if


+  !$OMP parallel do private(i) schedule(static)
  do i=1,size(adiag)
    if (adiag(i) /= zzero) then
      adiag(i) = zone / adiag(i)
@@ -226,7 +249,7 @@ subroutine amg_zaggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,&
      adiag(i) = zone
    end if
  end do
-
+  !$OMP end parallel do 
  if (parms%aggr_omega_alg == amg_eig_est_) then 

    if (parms%aggr_eig == amg_max_norm_) then 
@@ -252,8 +275,9 @@ subroutine amg_zaggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,&
    call psb_errpush(info,name,a_err='invalid amg_aggr_omega_alg_')
    goto 9999
  end if
+  if (do_timings) call psb_toc(idx_phase1)

-  
+  if (do_timings) call psb_tic(idx_phase2)
  call acsrf%scal(adiag,info)
  if (info /= psb_success_) goto 9999

@@ -267,6 +291,8 @@ subroutine amg_zaggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,&

  call psb_cdasb(desc_ac,info)
  call psb_cd_reinit(desc_ac,info)
+  if (do_timings) call psb_toc(idx_phase2)
+  if (do_timings) call psb_tic(idx_phase3)
  !
  ! Build the smoothed prolongator using either A or Af
  !    acsr1 = (I-w*D*A) Prol      acsr1 = (I-w*D*Af) Prol 
@@ -279,8 +305,8 @@ subroutine amg_zaggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,&
    call psb_errpush(psb_err_from_subroutine_,name,a_err='spspmm 1')
    goto 9999
  end if
-
-
+  if (do_timings) call psb_toc(idx_phase3)
+  if (do_timings) call psb_tic(idx_ptap)
  if (debug_level >= psb_debug_outer_) &
       & write(debug_unit,*) me,' ',trim(name),&
       & 'Done SPSPMM 1'
@@ -292,7 +318,7 @@ subroutine amg_zaggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,&

  call op_prol%mv_from(coo_prol)
  call op_restr%mv_from(coo_restr)
-
+  if (do_timings) call psb_toc(idx_ptap)
  if (debug_level >= psb_debug_outer_) &
       & write(debug_unit,*) me,' ',trim(name),&
       & 'Done smooth_aggregate '
@@ -0,0 +1,91 @@
+#include "MatchBoxPC.h"
+
+// TODO comment
+
+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)
+{
+    // Cleanup Phase
+
+#pragma omp parallel
+    {
+#pragma omp master
+        {
+#pragma omp task
+            {
+
+#ifdef PRINT_DEBUG_INFO_
+                cout << "\n(" << myRank << ") Waitall= " << endl;
+                fflush(stdout);
+#endif
+#ifdef DEBUG_HANG_
+                cout << "\n(" << myRank << ") Waitall " << endl;
+                fflush(stdout);
+#endif
+                //return;
+
+                MPI_Waitall(MessageIndex, &SRequest[0], &SStatus[0]);
+
+                // MPI_Buffer_attach(&Buffer, BufferSize); //Attach the Buffer
+                if (BufferSize > 0)
+                {
+                    MPI_Buffer_detach(&Buffer, &BufferSize); // Detach the Buffer
+                    free(Buffer);                            // Free the memory that was allocated
+                }
+            }
+
+#ifdef PRINT_DEBUG_INFO_
+            cout << "\n(" << myRank << ")End of function to compute matching: " << endl;
+            fflush(stdout);
+            cout << "\n(" << myRank << ")myCardinality: " << myCard << endl;
+            fflush(stdout);
+            cout << "\n(" << myRank << ")Matching took " << finishTime - startTime << "seconds" << endl;
+            fflush(stdout);
+            cout << "\n(" << myRank << ")** Getting out of the matching function **" << endl;
+            fflush(stdout);
+#endif
+#ifdef PRINT_DEBUG_INFO_
+            cout << "\n(" << myRank << ") Number of Ghost edges = " << numGhostEdges;
+            cout << "\n(" << myRank << ") Total number of potential message X 2 = " << numGhostEdges * 2;
+            cout << "\n(" << myRank << ") Number messages bundled = " << NumMessagesBundled;
+            cout << "\n(" << myRank << ") Total Individual Messages sent = " << msgInd;
+            if (msgInd > 0)
+            {
+                cout << "\n(" << myRank << ") Percentage of messages bundled = " << ((double)NumMessagesBundled / (double)(msgInd)) * 100.0 << "% \n";
+            }
+            fflush(stdout);
+#endif
+
+#pragma omp task
+            {
+                *msgActualSent = msgActual;
+                *msgIndSent = msgInd;
+                if (msgInd > 0)
+                {
+                    *msgPercent = ((double)NumMessagesBundled / (double)(msgInd)) * 100.0;
+                }
+                else
+                {
+                    *msgPercent = 0;
+                }
+            }
+
+#ifdef DEBUG_HANG_
+            if (myRank == 0)
+                cout << "\n(" << myRank << ") Done" << endl;
+            fflush(stdout);
+#endif
+        }
+    }
+}
@@ -0,0 +1,73 @@
+#include "MatchBoxPC.h"
+
+/**
+ * Execute the research fr the Candidate Mate without controlling if the vertices are already matched.
+ * Returns the vertices with the highest weight
+ * @param adj1
+ * @param adj2
+ * @param verLocInd
+ * @param edgeLocWeight
+ * @return
+ */
+MilanLongInt firstComputeCandidateMate(MilanLongInt adj1,
+                                       MilanLongInt adj2,
+                                       MilanLongInt *verLocInd,
+                                       MilanReal *edgeLocWeight)
+{
+    MilanInt w = -1;
+    MilanReal heaviestEdgeWt = MilanRealMin; // Assign the smallest Value possible first LDBL_MIN
+    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 computeCandidateMate(MilanLongInt adj1,
+                                  MilanLongInt adj2,
+                                  MilanReal *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;
+    MilanReal heaviestEdgeWt = MilanRealMin; // Assign the smallest Value possible first LDBL_MIN
+    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;
+}
@@ -80,9 +80,11 @@ class staticQueue
 		MilanLongInt squeueTail;
 		MilanLongInt NumNodes;

+		//FIXME I had to comment this piece of code in order to make everything work.
+		//		why?
 		//Prevent Assignment and Pass by Value:
-		staticQueue(const staticQueue& src);
-		staticQueue& operator=(const staticQueue& rhs);
+		//staticQueue(const staticQueue& src);
+		//staticQueue& operator=(const staticQueue& rhs);

 	public:
 		//Constructors and Destructors
@@ -0,0 +1,31 @@
+#include "MatchBoxPC.h"
+
+void extractUChunk(
+    vector<MilanLongInt> &UChunkBeingProcessed,
+    vector<MilanLongInt> &U,
+    vector<MilanLongInt> &privateU)
+{
+
+    UChunkBeingProcessed.clear();
+#pragma omp critical(U)
+    {
+
+        if (U.empty() && !privateU.empty()) // If U is empty but there are nodes in private U
+        {
+            while (!privateU.empty())
+                UChunkBeingProcessed.push_back(privateU.back());
+            privateU.pop_back();
+        }
+        else
+        {
+            for (int i = 0; i < UCHUNK; i++)
+            { // Pop the new nodes
+                if (U.empty())
+                    break;
+                UChunkBeingProcessed.push_back(U.back());
+                U.pop_back();
+            }
+        }
+
+    } // End of critical U // End of critical U
+}
@@ -0,0 +1,29 @@
+#include "MatchBoxPC.h"
+
+/// Find the owner of a ghost node:
+MilanInt findOwnerOfGhost(MilanLongInt vtxIndex, MilanLongInt *mVerDistance,
+                          MilanInt myRank, MilanInt numProcs)
+{
+
+  MilanLongInt mStartInd = mVerDistance[myRank];
+  MilanInt Start = 0;
+  MilanInt End = numProcs;
+  MilanInt Current = 0;
+
+  while (Start <= End)
+  {
+    Current = (End + Start) / 2;
+    // CASE-1:
+    if (mVerDistance[Current] == vtxIndex) return Current;
+    else // CASE 2:
+      if (mVerDistance[Current] > vtxIndex)
+        End = Current - 1;
+      else // CASE 3:
+        Start = Current + 1;
+  } // End of While()
+
+ if (mVerDistance[Current] > vtxIndex)
+      return (Current - 1);
+
+  return Current;
+} // End of findOwnerOfGhost()
@@ -0,0 +1,304 @@
+#include "MatchBoxPC.h"
+
+void initialize(MilanLongInt NLVer, MilanLongInt NLEdge,
+                MilanLongInt StartIndex, MilanLongInt EndIndex,
+                MilanLongInt *numGhostEdges,
+                MilanLongInt *numGhostVertices,
+                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)
+{
+
+    MilanLongInt insertMe = 0;
+    MilanLongInt adj1, adj2;
+    int i, v, k, w;
+    // index that starts with zero to |Vg|  - 1
+    map<MilanLongInt, MilanLongInt>::iterator storedAlready;
+
+#pragma omp parallel private(insertMe, k, w, v, adj1, adj2) firstprivate(StartIndex, EndIndex) default(shared) num_threads(NUM_THREAD)
+    {
+
+#pragma omp single
+        {
+
+#ifdef TIME_TRACKER
+            double Ghost2LocalInitialization = MPI_Wtime();
+#endif
+
+            /*
+             * OMP Ghost2LocalInitialization
+             * This loop analyzes all the edges and when finds a ghost edge
+             * puts it in the Ghost2LocalMap.
+             * A critical region is needed when inserting data in the map.
+             *
+             * Despite the critical region it is still productive to
+             * parallelize this cycle because the critical region is exeuted
+             * only when a ghost edge is found and ghost edges are a minority,
+             * circa 3.5% during the tests.
+             */
+#pragma omp task depend(out \
+                        : *numGhostEdges, Counter, Ghost2LocalMap, insertMe, storedAlready, *numGhostVertices)
+            {
+#pragma omp taskloop num_tasks(NUM_THREAD) reduction(+ \
+                                                     : numGhostEdges[:1])
+                for (i = 0; i < NLEdge; i++)
+                { // O(m) - Each edge stored twice
+                    insertMe = verLocInd[i];
+                    if ((insertMe < StartIndex) || (insertMe > EndIndex))
+                    { // Find a ghost
+                        (*numGhostEdges)++;
+#pragma omp critical
+                        {
+                            storedAlready = Ghost2LocalMap.find(insertMe);
+                            if (storedAlready != Ghost2LocalMap.end())
+                            {                                     // Has already been added
+                                Counter[storedAlready->second]++; // Increment the counter
+                            }
+                            else
+                            {                                                 // Insert an entry for the ghost:
+                                Ghost2LocalMap[insertMe] = *numGhostVertices; // Add a map entry
+                                Counter.push_back(1);                         // Initialize the counter
+                                (*numGhostVertices)++;                        // Increment the number of ghost vertices
+                            }                                                 // End of else()
+                        }
+                    } // End of if ( (insertMe < StartIndex) || (insertMe > EndIndex) )
+                }     // End of for(ghost vertices)
+            }         // end of task depend
+
+            // *numGhostEdges = atomicNumGhostEdges;
+#ifdef TIME_TRACKER
+            Ghost2LocalInitialization = MPI_Wtime() - Ghost2LocalInitialization;
+            fprintf(stderr, "Ghost2LocalInitialization time: %f\n", Ghost2LocalInitialization);
+#endif
+
+#ifdef PRINT_DEBUG_INFO_
+            cout << "\n(" << myRank << ")NGhosts:" << *numGhostVertices << " GhostEdges: " << *numGhostEdges;
+            if (!Ghost2LocalMap.empty())
+            {
+                cout << "\n(" << myRank << ")Final Map : on process ";
+                cout << "\n(" << myRank << ")Key \t Value \t Counter \n";
+                fflush(stdout);
+                storedAlready = Ghost2LocalMap.begin();
+                do
+                {
+                    cout << storedAlready->second << " - " << storedAlready->first << " : " << Counter[storedAlready->second] << endl;
+                    fflush(stdout);
+                    storedAlready++;
+                } while (storedAlready != Ghost2LocalMap.end());
+            }
+#endif
+
+#pragma omp task depend(out                                                       \
+                        : verGhostPtr, tempCounter, verGhostInd, GMate) depend(in \
+                                                                               : *numGhostVertices, *numGhostEdges)
+            {
+
+                // Initialize adjacency Lists for Ghost Vertices:
+                try
+                {
+                    verGhostPtr.reserve(*numGhostVertices + 1); // Pointer Vector
+                    tempCounter.reserve(*numGhostVertices);     // Pointer Vector
+                    verGhostInd.reserve(*numGhostEdges);        // Index Vector
+                    GMate.reserve(*numGhostVertices);           // Ghost Mate Vector
+                }
+                catch (length_error)
+                {
+                    cout << "Error in function algoDistEdgeApproxDominatingEdgesLinearSearch: \n";
+                    cout << "Not enough memory to allocate the internal variables \n";
+                    exit(1);
+                }
+                // Initialize the Vectors:
+                verGhostPtr.resize(*numGhostVertices + 1, 0); // Pointer Vector
+                tempCounter.resize(*numGhostVertices, 0);     // Temporary Counter
+                verGhostInd.resize(*numGhostEdges, -1);       // Index Vector
+                GMate.resize(*numGhostVertices, -1);          // Temporary Counter
+                verGhostPtr[0] = 0;                           // The first value
+#ifdef PRINT_DEBUG_INFO_
+                cout << "\n(" << myRank << ")Ghost Vertex Pointer: ";
+                fflush(stdout);
+#endif
+
+            } // End of task
+
+#pragma omp task depend(out                      \
+                        : verGhostPtr) depend(in \
+                                              : Counter, *numGhostVertices)
+            {
+
+#ifdef TIME_TRACKER
+                double verGhostPtrInitialization = MPI_Wtime();
+#endif
+                for (i = 0; i < *numGhostVertices; i++)
+                { // O(|Ghost Vertices|)
+                    verGhostPtr[i + 1] = verGhostPtr[i] + Counter[i];
+#ifdef PRINT_DEBUG_INFO_
+                    cout << verGhostPtr[i] << "\t";
+                    fflush(stdout);
+#endif
+                }
+
+#ifdef TIME_TRACKER
+                verGhostPtrInitialization = MPI_Wtime() - verGhostPtrInitialization;
+                fprintf(stderr, "verGhostPtrInitialization time: %f\n", verGhostPtrInitialization);
+#endif
+            } // End of task
+
+#ifdef PRINT_DEBUG_INFO_
+            if (*numGhostVertices > 0)
+                cout << verGhostPtr[*numGhostVertices] << "\n";
+            fflush(stdout);
+#endif
+
+#ifdef TIME_TRACKER
+            double verGhostIndInitialization = MPI_Wtime();
+#endif
+
+            /*
+             * OMP verGhostIndInitialization
+             *
+             * In this cycle the verGhostInd is initialized
+             * with the datas related to ghost edges.
+             * The check to see if a node is a ghost node is
+             * executed in paralle and when a ghost node
+             * is found a critical region is started.
+             *
+             * Despite the critical region it's still useful to
+             * parallelize the for cause the ghost nodes
+             * are a minority hence the critical region is executed
+             * few times, circa 3.5% of the times in the tests.
+             */
+#pragma omp task depend(in                                                               \
+                        : insertMe, Ghost2LocalMap, tempCounter, verGhostPtr) depend(out \
+                                                                                     : verGhostInd)
+            {
+#pragma omp taskloop num_tasks(NUM_THREAD)
+                for (v = 0; v < NLVer; v++)
+                {
+                    adj1 = verLocPtr[v]; // Vertex Pointer
+                    adj2 = verLocPtr[v + 1];
+                    for (k = adj1; k < adj2; k++)
+                    {
+                        w = verLocInd[k]; // Get the adjacent vertex
+                        if ((w < StartIndex) || (w > EndIndex))
+                        { // Find a ghost
+#pragma omp critical
+                            {
+                                insertMe = verGhostPtr[Ghost2LocalMap[w]] + tempCounter[Ghost2LocalMap[w]]; // Where to insert
+                                tempCounter[Ghost2LocalMap[w]]++;                                           // Increment the counter
+                            }
+                            verGhostInd[insertMe] = v + StartIndex; // Add the adjacency
+                        }                                           // End of if((w < StartIndex) || (w > EndIndex))
+                    }                                               // End of for(k)
+                }                                                   // End of for (v)
+            }                                                       // end of tasklopp
+
+#ifdef TIME_TRACKER
+            verGhostIndInitialization = MPI_Wtime() - verGhostIndInitialization;
+            fprintf(stderr, "verGhostIndInitialization time: %f\n", verGhostIndInitialization);
+#endif
+
+#ifdef PRINT_DEBUG_INFO_
+            cout << "\n(" << myRank << ")Ghost Vertex Index: ";
+            for (v = 0; v < *numGhostEdges; v++)
+                cout << verGhostInd[v] << "\t";
+            cout << endl;
+            fflush(stdout);
+#endif
+
+#pragma omp task depend(in                           \
+                        : *numGhostEdges) depend(out \
+                                                 : QLocalVtx, QGhostVtx, QMsgType, QOwner)
+            {
+                try
+                {
+                    QLocalVtx.reserve(*numGhostEdges); // Local Vertex
+                    QGhostVtx.reserve(*numGhostEdges); // Ghost Vertex
+                    QMsgType.reserve(*numGhostEdges);  // Message Type (Request/Failure)
+                    QOwner.reserve(*numGhostEdges);    // Owner of the ghost: COmpute once and use later
+                }
+                catch (length_error)
+                {
+                    cout << "Error in function algoDistEdgeApproxDominatingEdgesMessageBundling: \n";
+                    cout << "Not enough memory to allocate the internal variables \n";
+                    exit(1);
+                }
+            } // end of task
+
+#ifdef PRINT_DEBUG_INFO_
+            cout << "\n(" << myRank << ")Allocating CandidateMate.. ";
+            fflush(stdout);
+#endif
+
+#ifdef PRINT_DEBUG_INFO_
+            cout << "\n(" << myRank << "=========================************===============================" << endl;
+            fflush(stdout);
+            fflush(stdout);
+#endif
+
+#ifdef PRINT_DEBUG_INFO_
+            cout << "\n(" << myRank << ") Setup Time :" << *ph0_time << endl;
+            fflush(stdout);
+            fflush(stdout);
+#endif
+
+#ifdef DEBUG_HANG_
+            if (myRank == 0)
+                cout << "\n(" << myRank << ") Setup Time :" << *ph0_time << endl;
+            fflush(stdout);
+#endif
+
+#pragma omp task depend(in                              \
+                        : *numGhostVertices) depend(out \
+                                                    : candidateMate, S, U, privateU, privateQLocalVtx, privateQGhostVtx, privateQMsgType, privateQOwner)
+            {
+
+                // Allocate Data Structures:
+                /*
+                 * candidateMate was a vector and has been replaced with an array
+                 * there is no point in using the vector (or maybe there is (???))
+                 * so I replaced it with an array wich is slightly faster
+                 */
+                candidateMate = new MilanLongInt[NLVer + (*numGhostVertices)];
+
+                *S = (*numGhostVertices); // Initialize S with number of Ghost Vertices
+
+                /*
+                 * Create the Queue Data Structure for the Dominating Set
+                 *
+                 * I had to declare the staticuQueue U before the parallel region
+                 * to have it in the correct scope. Since we can't change the dimension
+                 * of a staticQueue I had to destroy the previous object and instantiate
+                 * a new one of the correct size.
+                 */
+                //new (&U) staticQueue(NLVer + (*numGhostVertices));
+                U.reserve(NLVer + (*numGhostVertices));
+
+                // Initialize the private vectors                
+                privateQLocalVtx.reserve(*numGhostVertices);
+                privateQGhostVtx.reserve(*numGhostVertices);
+                privateQMsgType.reserve(*numGhostVertices);
+                privateQOwner.reserve(*numGhostVertices);
+                privateU.reserve(*numGhostVertices);
+            } // end of task
+
+        } // End of single region
+    }     // End of parallel region
+}
@@ -0,0 +1,46 @@
+#include "MatchBoxPC.h"
+
+/**
+ * //TODO documentation
+ * @param k
+ * @param verLocInd
+ * @param StartIndex
+ * @param EndIndex
+ * @param GMate
+ * @param Mate
+ * @param Ghost2LocalMap
+ * @return
+ */
+bool isAlreadyMatched(MilanLongInt node,
+                      MilanLongInt StartIndex,
+                      MilanLongInt EndIndex,
+                      vector<MilanLongInt> &GMate,
+                      MilanLongInt *Mate,
+                      map<MilanLongInt, MilanLongInt> &Ghost2LocalMap)
+{
+
+    /*
+#pragma omp critical(Mate)
+    {
+        if ((node < StartIndex) || (node > EndIndex)) { //Is it a ghost vertex?
+            result = GMate[Ghost2LocalMap[node]] >= 0;// Already matched
+        } else { //A local vertex
+            result = (Mate[node - StartIndex] >= 0); // Already matched
+        }
+
+    }
+    */
+    MilanLongInt val;
+    if ((node < StartIndex) || (node > EndIndex)) // if ghost vertex
+    { 
+#pragma omp atomic read
+        val = GMate[Ghost2LocalMap[node]];
+        return val >= 0; // Already matched
+    }
+
+    // If not ghost vertex
+#pragma omp atomic read
+    val = Mate[node - StartIndex];
+
+    return val >= 0; // Already matched
+}
@@ -0,0 +1,27 @@
+#include "MatchBoxPC.h"
+
+void PARALLEL_COMPUTE_CANDIDATE_MATE_B(MilanLongInt NLVer,
+                                              MilanLongInt *verLocPtr,
+                                              MilanLongInt *verLocInd,
+                                              MilanInt myRank,
+                                              MilanReal *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] = firstComputeCandidateMate(verLocPtr[v], verLocPtr[v + 1], verLocInd, edgeLocWeight);
+            // End: PARALLEL_COMPUTE_CANDIDATE_MATE_B(v)
+        }
+    }
+}
@@ -0,0 +1,24 @@
+#include "MatchBoxPC.h"
+
+void PROCESS_CROSS_EDGE(MilanLongInt *edge,
+                        MilanLongInt *S)
+{
+    // Start: PARALLEL_PROCESS_CROSS_EDGE_B
+    MilanLongInt captureCounter;
+
+#pragma omp atomic capture
+    captureCounter = --(*edge); // Decrement
+
+    //assert(captureCounter >= 0);
+
+    if (captureCounter == 0)
+#pragma omp atomic
+        (*S)--; // Decrement S
+
+#ifdef PRINT_DEBUG_INFO_
+    cout << "\n(" << myRank << ")Decrementing S: Ghost vertex " << edge << " has received all its messages";
+    fflush(stdout);
+#endif
+
+    // End: PARALLEL_PROCESS_CROSS_EDGE_B
+}
@@ -0,0 +1,195 @@
+#include "MatchBoxPC.h"
+
+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 *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(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
+			  
+#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 critical processExposed
+		
+            } // 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
+
+#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);
+		  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);
+                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
+}
@@ -0,0 +1,294 @@
+#include "MatchBoxPC.h"
+
+void processMatchedVertices(
+    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,
+    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)
+{
+
+    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
+		    {
+		      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;
+			
+#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);
+		    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);
+			
+			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());
+	  }
+	  
+	  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());
+	  }
+	  
+	  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
+}
@@ -0,0 +1,308 @@
+#include "MatchBoxPC.h"
+//#define DEBUG_HANG_
+void processMatchedVerticesAndSendMessages(
+					   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,
+					   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)
+{
+
+    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
+		    {
+		      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;
+			
+#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,
+               omp_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
+}
@@ -0,0 +1,315 @@
+#include "MatchBoxPC.h"
+//#define DEBUG_HANG_
+
+void processMessages(
+    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,
+    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 *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 = computeCandidateMate(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;
+}
@@ -0,0 +1,36 @@
+#include "MatchBoxPC.h"
+
+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)
+{
+
+#pragma omp critical(U)
+    {
+        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());
+    }
+
+    privateQLocalVtx.clear();
+    privateQGhostVtx.clear();
+    privateQMsgType.clear();
+    privateQOwner.clear();
+
+}
@@ -0,0 +1,209 @@
+#include "MatchBoxPC.h"
+
+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)
+{
+
+    MilanLongInt myIndex = 0, numMessagesToSend;
+    MilanInt i = 0, OneMessageSize = 0;
+
+#ifdef DEBUG_HANG_
+    if (myRank == 0)
+        cout << "\n(" << myRank << ") Send Bundles" << endl;
+    fflush(stdout);
+#endif
+
+#pragma omp parallel private(i) default(shared) num_threads(NUM_THREAD)
+    {
+#pragma omp master
+        {
+// Data structures for Bundled Messages:
+#pragma omp task depend(inout                                                       \
+                        : PCumulative, PMessageBundle, PSizeInfoMessages) depend(in \
+                                                                                 : NumMessagesBundled, numProcs)
+            {
+	      try {
+                PMessageBundle.reserve(NumMessagesBundled * 3); // Three integers per message
+		PCumulative.reserve(numProcs + 1);                          // Similar to Row Pointer vector in CSR data structure
+		PSizeInfoMessages.reserve(numProcs * 3);                    // Buffer to hold the Size info message packets
+	      }
+	      catch (length_error)
+		{
+		  cout << "Error in function algoDistEdgeApproxDominatingEdgesMessageBundling: \n";
+		  cout << "Not enough memory to allocate the internal variables \n";
+		  exit(1);
+		}
+	      PMessageBundle.resize(NumMessagesBundled * 3, -1); // Initialize
+	      PCumulative.resize(numProcs + 1, 0);               // Only initialize the counter variable
+	      PSizeInfoMessages.resize(numProcs * 3, 0);
+	    }
+
+#pragma omp task depend(inout                    \
+                        : PCumulative) depend(in		\
+                                              : PCounter)
+	    {
+	      for (i = 0; i < numProcs; i++)
+		PCumulative[i + 1] = PCumulative[i] + PCounter[i];
+	    }
+	    
+#pragma omp task depend(inout \
+                        : PCounter)
+	    {
+	      // Reuse PCounter to keep track of how many messages were inserted:
+	      for (MilanInt i = 0; i < numProcs; i++) // Changed by Fabio to be an integer, addresses needs to be integers!
+		PCounter[i] = 0;
+	    }
+
+// Build the Message Bundle packet:
+#pragma omp task depend(in                                                                                          \
+                        : PCounter, QLocalVtx, QGhostVtx, QMsgType, QOwner, PMessageBundle, PCumulative) depend(out \
+                                                                                                                : myIndex, PMessageBundle, PCounter)
+{
+  for (i = 0; i < NumMessagesBundled; i++)   {
+    myIndex = (PCumulative[QOwner[i]] + PCounter[QOwner[i]]) * 3;
+    PMessageBundle[myIndex + 0] = QLocalVtx[i];
+    PMessageBundle[myIndex + 1] = QGhostVtx[i];
+    PMessageBundle[myIndex + 2] = QMsgType[i];
+    PCounter[QOwner[i]]++;
+  }
+ }
+ 
+// Send the Bundled Messages: Use ISend
+#pragma omp task depend(out \
+                        : SRequest, SStatus)
+ {
+   try
+     {
+       SRequest.reserve(numProcs * 2); // At most two messages per processor
+       SStatus.reserve(numProcs * 2);  // At most two messages per processor
+     }
+   catch (length_error)
+     {
+       cout << "Error in function algoDistEdgeApproxDominatingEdgesLinearSearchImmediateSend: \n";
+       cout << "Not enough memory to allocate the internal variables \n";
+       exit(1);
+     }
+ }
+ 
+// Send the Messages
+#pragma omp task depend(inout                                                  \
+                        : SRequest, PSizeInfoMessages, PCumulative) depend(out \
+                                                                           : *msgActual, *msgInd)
+{
+  for (i = 0; i < numProcs; i++)   {                    // Changed by Fabio to be an integer, addresses needs to be integers!
+    if (i == myRank) // Do not send anything to yourself
+      continue;
+    // Send the Message with information about the size of next message:
+    // Build the Message Packet:
+    PSizeInfoMessages[i * 3 + 0] = (PCumulative[i + 1] - PCumulative[i]) * 3; // # of integers in the next message
+    PSizeInfoMessages[i * 3 + 1] = -1;                                        // Dummy packet
+    PSizeInfoMessages[i * 3 + 2] = SIZEINFO;                                  // TYPE
+    // Send a Request (Asynchronous)
+#ifdef PRINT_DEBUG_INFO_
+    cout << "\n(" << myRank << ")Sending bundled message to process " << i << " size: " << PSizeInfoMessages[i * 3 + 0] << endl;
+    fflush(stdout);
+#endif
+    if (PSizeInfoMessages[i * 3 + 0] > 0)
+      { // Send only if it is a nonempty packet
+	MPI_Isend(&PSizeInfoMessages[i * 3 + 0], 3, TypeMap<MilanLongInt>(), i, ComputeTag, comm,
+		  &SRequest[(*msgInd)]);
+	(*msgActual)++;
+	(*msgInd)++;
+	// Now Send the message with the data packet:
+#ifdef PRINT_DEBUG_INFO_
+	cout << "\n(" << myRank << ")SendiFFng Bundle to : " << i << endl;
+	for (k = (PCumulative[i] * 3); k < (PCumulative[i] * 3 + PSizeInfoMessages[i * 3 + 0]); k++)
+	  cout << PMessageBundle[k] << ",";
+	cout << endl;
+	fflush(stdout);
+#endif
+	MPI_Isend(&PMessageBundle[PCumulative[i] * 3], PSizeInfoMessages[i * 3 + 0],
+		  TypeMap<MilanLongInt>(), i, BundleTag, comm, &SRequest[(*msgInd)]);
+	(*msgInd)++;
+      } // End of if size > 0
+  }
+}
+
+#pragma omp task depend(inout \
+                        : PCumulative, QLocalVtx, QGhostVtx, QMsgType, QOwner)
+{
+
+    // Free up temporary memory:
+  PCumulative.clear();
+  QLocalVtx.clear();
+  QGhostVtx.clear();
+  QMsgType.clear();
+  QOwner.clear();
+}
+
+#pragma omp task depend(inout : OneMessageSize, *BufferSize) depend(out : numMessagesToSend) depend(in : *numGhostEdges)
+{
+
+#ifdef PRINT_DEBUG_INFO_
+  cout << "\n(" << myRank << ")Number of Ghost edges = " << *numGhostEdges;
+  cout << "\n(" << myRank << ")Total number of potential message X 2 = " << *numGhostEdges * 2;
+  cout << "\n(" << myRank << ")Number messages already sent in bundles = " << NumMessagesBundled;
+  if (*numGhostEdges > 0)
+    {
+      cout << "\n(" << myRank << ")Percentage of total = " << ((double)NumMessagesBundled / (double)(*numGhostEdges * 2)) * 100.0 << "% \n";
+    }
+  fflush(stdout);
+#endif
+
+  // Allocate memory for MPI Send messages:
+  /* WILL COME BACK HERE - NO NEED TO STORE ALL THIS MEMORY !! */
+  OneMessageSize = 0;
+  MPI_Pack_size(3, TypeMap<MilanLongInt>(), comm, &OneMessageSize); // Size of one message packet
+  // How many messages to send?
+  // Potentially three kinds of messages will be sent/received:
+  // Request, Success, Failure.
+  // But only two will be sent from a given processor.
+  // Substract the number of messages that have already been sent as bundled messages:
+  numMessagesToSend = (*numGhostEdges) * 2 - NumMessagesBundled;
+  *BufferSize = (OneMessageSize + MPI_BSEND_OVERHEAD) * numMessagesToSend;
+}
+
+#pragma omp task depend(out : Buffer) depend(in : *BufferSize)
+ {
+   Buffer = 0;
+#ifdef PRINT_DEBUG_INFO_
+   cout << "\n(" << myRank << ")Size of One Message from PACK= " << OneMessageSize;
+   cout << "\n(" << myRank << ")Size of Message overhead = " << MPI_BSEND_OVERHEAD;
+   cout << "\n(" << myRank << ")Number of Ghost edges = " << *numGhostEdges;
+   cout << "\n(" << myRank << ")Number of remaining message = " << numMessagesToSend;
+   cout << "\n(" << myRank << ")BufferSize = " << (*BufferSize);
+   cout << "\n(" << myRank << ")Attaching Buffer on.. ";
+   fflush(stdout);
+#endif
+   if ((*BufferSize) > 0)
+     {
+       Buffer = (MilanLongInt *)malloc((*BufferSize)); // Allocate memory
+       if (Buffer == 0)
+	 {
+	   cout << "Error in function algoDistEdgeApproxDominatingEdgesLinearSearch: \n";
+	   cout << "Not enough memory to allocate for send buffer on process " << myRank << "\n";
+	   exit(1);
+	 }
+       MPI_Buffer_attach(Buffer, *BufferSize); // Attach the Buffer
+     }
+ }
+}
+}
+}
@@ -109,6 +109,8 @@ subroutine amg_c_base_onelev_mat_asb(lv,a,desc_a,ilaggr,nlaggr,t_prol,info)
  type(psb_cspmat_type)    :: ac, op_restr, op_prol
  integer(psb_ipk_)        :: nzl, inl
  integer(psb_ipk_)        :: debug_level, debug_unit
+  integer(psb_ipk_), save  :: idx_matbld=-1, idx_matasb=-1, idx_mapbld=-1
+  logical, parameter :: do_timings=.false.

  name='amg_c_onelev_mat_asb'
  call psb_erractionsave(err_act)
@@ -120,6 +122,12 @@ subroutine amg_c_base_onelev_mat_asb(lv,a,desc_a,ilaggr,nlaggr,t_prol,info)
  info  = psb_success_
  ctxt = desc_a%get_context()
  call psb_info(ctxt,me,np)
+  if ((do_timings).and.(idx_matbld==-1))       &
+       & idx_matbld = psb_get_timer_idx("LEV_MASB: mat_bld")
+  if ((do_timings).and.(idx_matasb==-1))     &
+       & idx_matasb = psb_get_timer_idx("LEV_MASB: mat_asb")
+  if ((do_timings).and.(idx_mapbld==-1))       &
+       & idx_mapbld = psb_get_timer_idx("LEV_MASB: map_bld")

  call amg_check_def(lv%parms%aggr_prol,'Smoother',&
       &   amg_smooth_prol_,is_legal_ml_aggr_prol)
@@ -139,9 +147,10 @@ subroutine amg_c_base_onelev_mat_asb(lv,a,desc_a,ilaggr,nlaggr,t_prol,info)
  ! the mapping defined by amg_aggrmap_bld and applying the aggregation
  ! algorithm specified by lv%iprcparm(amg_aggr_prol_)
  !
+  if (do_timings) call psb_tic(idx_matbld)
  call lv%aggr%mat_bld(lv%parms,a,desc_a,ilaggr,nlaggr,&
       & lv%ac,lv%desc_ac,op_prol,op_restr,t_prol,info)
-
+  if (do_timings) call psb_toc(idx_matbld)
  if(info /= psb_success_) then
    call psb_errpush(psb_err_from_subroutine_,name,a_err='amg_aggrmat_asb')
    goto 9999
@@ -151,14 +160,17 @@ subroutine amg_c_base_onelev_mat_asb(lv,a,desc_a,ilaggr,nlaggr,t_prol,info)
  ! Now build its descriptor and convert global indices for
  ! ac, op_restr and op_prol
  !
+  if (do_timings) call psb_tic(idx_matasb)
  if (info == psb_success_) &
       & call lv%aggr%mat_asb(lv%parms,a,desc_a,&
       & lv%ac,lv%desc_ac,op_prol,op_restr,info)
-  
+  if (do_timings) call psb_toc(idx_matasb)
+  if (do_timings) call psb_tic(idx_mapbld)  
  if (info == psb_success_) call lv%ac%cscnv(info,type='csr',dupl=psb_dupl_add_)
  
  if (info == psb_success_) call lv%aggr%bld_map(desc_a, lv%desc_ac,&
       & ilaggr,nlaggr,op_restr,op_prol,lv%linmap,info)
+  if (do_timings) call psb_toc(idx_mapbld)  
  if(info /= psb_success_) then
    call psb_errpush(psb_err_from_subroutine_,name,a_err='mat_asb/map_bld')
    goto 9999
@@ -109,6 +109,8 @@ subroutine amg_d_base_onelev_mat_asb(lv,a,desc_a,ilaggr,nlaggr,t_prol,info)
  type(psb_dspmat_type)    :: ac, op_restr, op_prol
  integer(psb_ipk_)        :: nzl, inl
  integer(psb_ipk_)        :: debug_level, debug_unit
+  integer(psb_ipk_), save  :: idx_matbld=-1, idx_matasb=-1, idx_mapbld=-1
+  logical, parameter :: do_timings=.false.

  name='amg_d_onelev_mat_asb'
  call psb_erractionsave(err_act)
@@ -120,6 +122,12 @@ subroutine amg_d_base_onelev_mat_asb(lv,a,desc_a,ilaggr,nlaggr,t_prol,info)
  info  = psb_success_
  ctxt = desc_a%get_context()
  call psb_info(ctxt,me,np)
+  if ((do_timings).and.(idx_matbld==-1))       &
+       & idx_matbld = psb_get_timer_idx("LEV_MASB: mat_bld")
+  if ((do_timings).and.(idx_matasb==-1))     &
+       & idx_matasb = psb_get_timer_idx("LEV_MASB: mat_asb")
+  if ((do_timings).and.(idx_mapbld==-1))       &
+       & idx_mapbld = psb_get_timer_idx("LEV_MASB: map_bld")

  call amg_check_def(lv%parms%aggr_prol,'Smoother',&
       &   amg_smooth_prol_,is_legal_ml_aggr_prol)
@@ -139,9 +147,10 @@ subroutine amg_d_base_onelev_mat_asb(lv,a,desc_a,ilaggr,nlaggr,t_prol,info)
  ! the mapping defined by amg_aggrmap_bld and applying the aggregation
  ! algorithm specified by lv%iprcparm(amg_aggr_prol_)
  !
+  if (do_timings) call psb_tic(idx_matbld)
  call lv%aggr%mat_bld(lv%parms,a,desc_a,ilaggr,nlaggr,&
       & lv%ac,lv%desc_ac,op_prol,op_restr,t_prol,info)
-
+  if (do_timings) call psb_toc(idx_matbld)
  if(info /= psb_success_) then
    call psb_errpush(psb_err_from_subroutine_,name,a_err='amg_aggrmat_asb')
    goto 9999
@@ -151,14 +160,17 @@ subroutine amg_d_base_onelev_mat_asb(lv,a,desc_a,ilaggr,nlaggr,t_prol,info)
  ! Now build its descriptor and convert global indices for
  ! ac, op_restr and op_prol
  !
+  if (do_timings) call psb_tic(idx_matasb)
  if (info == psb_success_) &
       & call lv%aggr%mat_asb(lv%parms,a,desc_a,&
       & lv%ac,lv%desc_ac,op_prol,op_restr,info)
-  
+  if (do_timings) call psb_toc(idx_matasb)
+  if (do_timings) call psb_tic(idx_mapbld)  
  if (info == psb_success_) call lv%ac%cscnv(info,type='csr',dupl=psb_dupl_add_)
  
  if (info == psb_success_) call lv%aggr%bld_map(desc_a, lv%desc_ac,&
       & ilaggr,nlaggr,op_restr,op_prol,lv%linmap,info)
+  if (do_timings) call psb_toc(idx_mapbld)  
  if(info /= psb_success_) then
    call psb_errpush(psb_err_from_subroutine_,name,a_err='mat_asb/map_bld')
    goto 9999
@@ -109,6 +109,8 @@ subroutine amg_s_base_onelev_mat_asb(lv,a,desc_a,ilaggr,nlaggr,t_prol,info)
  type(psb_sspmat_type)    :: ac, op_restr, op_prol
  integer(psb_ipk_)        :: nzl, inl
  integer(psb_ipk_)        :: debug_level, debug_unit
+  integer(psb_ipk_), save  :: idx_matbld=-1, idx_matasb=-1, idx_mapbld=-1
+  logical, parameter :: do_timings=.false.

  name='amg_s_onelev_mat_asb'
  call psb_erractionsave(err_act)
@@ -120,6 +122,12 @@ subroutine amg_s_base_onelev_mat_asb(lv,a,desc_a,ilaggr,nlaggr,t_prol,info)
  info  = psb_success_
  ctxt = desc_a%get_context()
  call psb_info(ctxt,me,np)
+  if ((do_timings).and.(idx_matbld==-1))       &
+       & idx_matbld = psb_get_timer_idx("LEV_MASB: mat_bld")
+  if ((do_timings).and.(idx_matasb==-1))     &
+       & idx_matasb = psb_get_timer_idx("LEV_MASB: mat_asb")
+  if ((do_timings).and.(idx_mapbld==-1))       &
+       & idx_mapbld = psb_get_timer_idx("LEV_MASB: map_bld")

  call amg_check_def(lv%parms%aggr_prol,'Smoother',&
       &   amg_smooth_prol_,is_legal_ml_aggr_prol)
@@ -139,9 +147,10 @@ subroutine amg_s_base_onelev_mat_asb(lv,a,desc_a,ilaggr,nlaggr,t_prol,info)
  ! the mapping defined by amg_aggrmap_bld and applying the aggregation
  ! algorithm specified by lv%iprcparm(amg_aggr_prol_)
  !
+  if (do_timings) call psb_tic(idx_matbld)
  call lv%aggr%mat_bld(lv%parms,a,desc_a,ilaggr,nlaggr,&
       & lv%ac,lv%desc_ac,op_prol,op_restr,t_prol,info)
-
+  if (do_timings) call psb_toc(idx_matbld)
  if(info /= psb_success_) then
    call psb_errpush(psb_err_from_subroutine_,name,a_err='amg_aggrmat_asb')
    goto 9999
@@ -151,14 +160,17 @@ subroutine amg_s_base_onelev_mat_asb(lv,a,desc_a,ilaggr,nlaggr,t_prol,info)
  ! Now build its descriptor and convert global indices for
  ! ac, op_restr and op_prol
  !
+  if (do_timings) call psb_tic(idx_matasb)
  if (info == psb_success_) &
       & call lv%aggr%mat_asb(lv%parms,a,desc_a,&
       & lv%ac,lv%desc_ac,op_prol,op_restr,info)
-  
+  if (do_timings) call psb_toc(idx_matasb)
+  if (do_timings) call psb_tic(idx_mapbld)  
  if (info == psb_success_) call lv%ac%cscnv(info,type='csr',dupl=psb_dupl_add_)
  
  if (info == psb_success_) call lv%aggr%bld_map(desc_a, lv%desc_ac,&
       & ilaggr,nlaggr,op_restr,op_prol,lv%linmap,info)
+  if (do_timings) call psb_toc(idx_mapbld)  
  if(info /= psb_success_) then
    call psb_errpush(psb_err_from_subroutine_,name,a_err='mat_asb/map_bld')
    goto 9999
@@ -109,6 +109,8 @@ subroutine amg_z_base_onelev_mat_asb(lv,a,desc_a,ilaggr,nlaggr,t_prol,info)
  type(psb_zspmat_type)    :: ac, op_restr, op_prol
  integer(psb_ipk_)        :: nzl, inl
  integer(psb_ipk_)        :: debug_level, debug_unit
+  integer(psb_ipk_), save  :: idx_matbld=-1, idx_matasb=-1, idx_mapbld=-1
+  logical, parameter :: do_timings=.false.

  name='amg_z_onelev_mat_asb'
  call psb_erractionsave(err_act)
@@ -120,6 +122,12 @@ subroutine amg_z_base_onelev_mat_asb(lv,a,desc_a,ilaggr,nlaggr,t_prol,info)
  info  = psb_success_
  ctxt = desc_a%get_context()
  call psb_info(ctxt,me,np)
+  if ((do_timings).and.(idx_matbld==-1))       &
+       & idx_matbld = psb_get_timer_idx("LEV_MASB: mat_bld")
+  if ((do_timings).and.(idx_matasb==-1))     &
+       & idx_matasb = psb_get_timer_idx("LEV_MASB: mat_asb")
+  if ((do_timings).and.(idx_mapbld==-1))       &
+       & idx_mapbld = psb_get_timer_idx("LEV_MASB: map_bld")

  call amg_check_def(lv%parms%aggr_prol,'Smoother',&
       &   amg_smooth_prol_,is_legal_ml_aggr_prol)
@@ -139,9 +147,10 @@ subroutine amg_z_base_onelev_mat_asb(lv,a,desc_a,ilaggr,nlaggr,t_prol,info)
  ! the mapping defined by amg_aggrmap_bld and applying the aggregation
  ! algorithm specified by lv%iprcparm(amg_aggr_prol_)
  !
+  if (do_timings) call psb_tic(idx_matbld)
  call lv%aggr%mat_bld(lv%parms,a,desc_a,ilaggr,nlaggr,&
       & lv%ac,lv%desc_ac,op_prol,op_restr,t_prol,info)
-
+  if (do_timings) call psb_toc(idx_matbld)
  if(info /= psb_success_) then
    call psb_errpush(psb_err_from_subroutine_,name,a_err='amg_aggrmat_asb')
    goto 9999
@@ -151,14 +160,17 @@ subroutine amg_z_base_onelev_mat_asb(lv,a,desc_a,ilaggr,nlaggr,t_prol,info)
  ! Now build its descriptor and convert global indices for
  ! ac, op_restr and op_prol
  !
+  if (do_timings) call psb_tic(idx_matasb)
  if (info == psb_success_) &
       & call lv%aggr%mat_asb(lv%parms,a,desc_a,&
       & lv%ac,lv%desc_ac,op_prol,op_restr,info)
-  
+  if (do_timings) call psb_toc(idx_matasb)
+  if (do_timings) call psb_tic(idx_mapbld)  
  if (info == psb_success_) call lv%ac%cscnv(info,type='csr',dupl=psb_dupl_add_)
  
  if (info == psb_success_) call lv%aggr%bld_map(desc_a, lv%desc_ac,&
       & ilaggr,nlaggr,op_restr,op_prol,lv%linmap,info)
+  if (do_timings) call psb_toc(idx_mapbld)  
  if(info /= psb_success_) then
    call psb_errpush(psb_err_from_subroutine_,name,a_err='mat_asb/map_bld')
    goto 9999
@@ -0,0 +1,25 @@
+cd amgprec/impl/aggregator/
+rm MatchBoxPC.o 
+rm sendBundledMessages.o 
+rm initialize.o 
+rm extractUChunk.o 
+rm isAlreadyMatched.o 
+rm findOwnerOfGhost.o 
+rm computeCandidateMate.o 
+rm parallelComputeCandidateMateB.o 
+rm processMatchedVertices.o 
+rm processCrossEdge.o 
+rm queueTransfer.o 
+rm processMessages.o 
+rm processExposedVertex.o 
+rm algoDistEdgeApproxDomEdgesLinearSearchMesgBndlSmallMateC.o 
+rm algoDistEdgeApproxDomEdgesLinearSearchMesgBndlSmallMateCMP.o
+cd ../../../
+make all
+cd samples/advanced/pdegen
+make amg_d_pde3d
+cd runs
+mpirun -np 4 amg_d_pde3d amg_pde3d.inp
+
+
+
@@ -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_krylov -lamg_prec -lpsb_prec 
 FINCLUDES=$(FMFLAG). $(FMFLAG)$(AMGMODDIR) $(FMFLAG)$(AMGINCDIR) $(PSBLAS_INCLUDES) $(FIFLAG).

 LINKOPT=
@@ -93,6 +93,9 @@ contains
       & a1,a2,a3,b1,b2,b3,c,g,info,f,amold,vmold,partition, nrl,iv)
    use psb_base_mod
    use psb_util_mod
+#if defined(OPENMP)
+    use omp_lib
+#endif
    !
    !   Discretizes the partial differential equation
    !
@@ -128,7 +131,6 @@ contains
    type(psb_d_csc_sparse_mat)  :: acsc
    type(psb_d_coo_sparse_mat)  :: acoo
    type(psb_d_csr_sparse_mat)  :: acsr
-    real(psb_dpk_)           :: zt(nb),x,y,z,xph,xmh,yph,ymh,zph,zmh
    integer(psb_ipk_) :: nnz,nr,nlr,i,j,ii,ib,k, partition_
    integer(psb_lpk_) :: m,n,glob_row,nt
    integer(psb_ipk_) :: ix,iy,iz,ia,indx_owner
@@ -141,8 +143,7 @@ contains
    ! Process grid
    integer(psb_ipk_) :: np, iam
    integer(psb_ipk_) :: icoeff
-    integer(psb_lpk_), allocatable     :: irow(:),icol(:),myidx(:)
-    real(psb_dpk_), allocatable :: val(:)
+    integer(psb_lpk_), allocatable     :: myidx(:)
    ! deltah dimension of each grid cell
    ! deltat discretization time
    real(psb_dpk_)            :: deltah, sqdeltah, deltah2
@@ -368,119 +369,128 @@ contains
    call psb_barrier(ctxt)
    talc = psb_wtime()-t0

-    if (info /= psb_success_) then
-      info=psb_err_from_subroutine_
-      ch_err='allocation rout.'
-      call psb_errpush(info,name,a_err=ch_err)
-      goto 9999
-    end if
-
-    ! we build an auxiliary matrix consisting of one row at a
-    ! time; just a small matrix. might be extended to generate
-    ! a bunch of rows per call.
-    !
-    allocate(val(20*nb),irow(20*nb),&
-         &icol(20*nb),stat=info)
-    if (info /= psb_success_ ) then
-      info=psb_err_alloc_dealloc_
-      call psb_errpush(info,name)
-      goto 9999
-    endif
-
-
-    ! loop over rows belonging to current process in a block
-    ! distribution.
-
    call psb_barrier(ctxt)
    t1 = psb_wtime()
-    do ii=1, nlr,nb
-      ib = min(nb,nlr-ii+1)
-      icoeff = 1
-      do k=1,ib
-        i=ii+k-1
-        ! local matrix pointer
-        glob_row=myidx(i)
-        ! compute gridpoint coordinates
-        call idx2ijk(ix,iy,iz,glob_row,idim,idim,idim)
-        ! x, y, z coordinates
-        x = (ix-1)*deltah
-        y = (iy-1)*deltah
-        z = (iz-1)*deltah
-        zt(k) = f_(x,y,z)
-        ! internal point: build discretization
-        !
-        !  term depending on   (x-1,y,z)
-        !
-        val(icoeff) = -a1(x,y,z)/sqdeltah-b1(x,y,z)/deltah2
-        if (ix == 1) then
-          zt(k) = g(dzero,y,z)*(-val(icoeff)) + zt(k)
-        else
-          call ijk2idx(icol(icoeff),ix-1,iy,iz,idim,idim,idim)
-          irow(icoeff) = glob_row
-          icoeff       = icoeff+1
-        endif
-        !  term depending on     (x,y-1,z)
-        val(icoeff)  = -a2(x,y,z)/sqdeltah-b2(x,y,z)/deltah2
-        if (iy == 1) then
-          zt(k) = g(x,dzero,z)*(-val(icoeff))   + zt(k)
-        else
-          call ijk2idx(icol(icoeff),ix,iy-1,iz,idim,idim,idim)
-          irow(icoeff) = glob_row
-          icoeff       = icoeff+1
-        endif
-        !  term depending on     (x,y,z-1)
-        val(icoeff)=-a3(x,y,z)/sqdeltah-b3(x,y,z)/deltah2
-        if (iz == 1) then
-          zt(k) = g(x,y,dzero)*(-val(icoeff))   + zt(k)
-        else
-          call ijk2idx(icol(icoeff),ix,iy,iz-1,idim,idim,idim)
-          irow(icoeff) = glob_row
-          icoeff       = icoeff+1
-        endif
+    !$omp parallel shared(deltah,myidx,a,desc_a)
+    !
+    block 
+      integer(psb_ipk_) :: i,j,k,ii,ib,icoeff, ix,iy,iz, ith,nth
+      integer(psb_lpk_) :: glob_row
+      integer(psb_lpk_), allocatable :: irow(:),icol(:)
+      real(psb_dpk_), allocatable :: val(:)
+      real(psb_dpk_)     :: x,y,z, zt(nb)
+#if defined(OPENMP)
+      nth = omp_get_num_threads()
+      ith = omp_get_thread_num()
+#else
+      nth = 1
+      ith = 0
+#endif
+      allocate(val(20*nb),irow(20*nb),&
+           &icol(20*nb),stat=info)
+      if (info /= psb_success_ ) then
+        info=psb_err_alloc_dealloc_
+        call psb_errpush(info,name)
+        !goto 9999
+      endif

-        !  term depending on     (x,y,z)
-        val(icoeff)=(2*done)*(a1(x,y,z)+a2(x,y,z)+a3(x,y,z))/sqdeltah &
-             & + c(x,y,z)
-        call ijk2idx(icol(icoeff),ix,iy,iz,idim,idim,idim)
-        irow(icoeff) = glob_row
-        icoeff       = icoeff+1
-        !  term depending on     (x,y,z+1)
-        val(icoeff)=-a3(x,y,z)/sqdeltah+b3(x,y,z)/deltah2
-        if (iz == idim) then
-          zt(k) = g(x,y,done)*(-val(icoeff))   + zt(k)
-        else
-          call ijk2idx(icol(icoeff),ix,iy,iz+1,idim,idim,idim)
-          irow(icoeff) = glob_row
-          icoeff       = icoeff+1
-        endif
-        !  term depending on     (x,y+1,z)
-        val(icoeff)=-a2(x,y,z)/sqdeltah+b2(x,y,z)/deltah2
-        if (iy == idim) then
-          zt(k) = g(x,done,z)*(-val(icoeff))   + zt(k)
-        else
-          call ijk2idx(icol(icoeff),ix,iy+1,iz,idim,idim,idim)
-          irow(icoeff) = glob_row
-          icoeff       = icoeff+1
-        endif
-        !  term depending on     (x+1,y,z)
-        val(icoeff)=-a1(x,y,z)/sqdeltah+b1(x,y,z)/deltah2
-        if (ix==idim) then
-          zt(k) = g(done,y,z)*(-val(icoeff))   + zt(k)
-        else
-          call ijk2idx(icol(icoeff),ix+1,iy,iz,idim,idim,idim)
-          irow(icoeff) = glob_row
-          icoeff       = icoeff+1
-        endif
+      !$omp  do schedule(dynamic)
+      !     
+      do ii=1, nlr, nb
+        if (info /= psb_success_) cycle
+        ib = min(nb,nlr-ii+1)
+        icoeff = 1
+        do k=1,ib
+          i=ii+k-1
+          ! local matrix pointer
+          glob_row=myidx(i)
+          ! compute gridpoint coordinates
+          call idx2ijk(ix,iy,iz,glob_row,idim,idim,idim)
+          ! x, y, z coordinates
+          x = (ix-1)*deltah
+          y = (iy-1)*deltah
+          z = (iz-1)*deltah
+          zt(k) = f_(x,y,z)
+          ! internal point: build discretization
+          !
+          !  term depending on   (x-1,y,z)
+          !
+          val(icoeff) = -a1(x,y,z)/sqdeltah-b1(x,y,z)/deltah2
+          if (ix == 1) then
+            zt(k) = g(dzero,y,z)*(-val(icoeff)) + zt(k)
+          else
+            call ijk2idx(icol(icoeff),ix-1,iy,iz,idim,idim,idim)
+            irow(icoeff) = glob_row
+            icoeff       = icoeff+1
+          endif
+          !  term depending on     (x,y-1,z)
+          val(icoeff)  = -a2(x,y,z)/sqdeltah-b2(x,y,z)/deltah2
+          if (iy == 1) then
+            zt(k) = g(x,dzero,z)*(-val(icoeff))   + zt(k)
+          else
+            call ijk2idx(icol(icoeff),ix,iy-1,iz,idim,idim,idim)
+            irow(icoeff) = glob_row
+            icoeff       = icoeff+1
+          endif
+          !  term depending on     (x,y,z-1)
+          val(icoeff)=-a3(x,y,z)/sqdeltah-b3(x,y,z)/deltah2
+          if (iz == 1) then
+            zt(k) = g(x,y,dzero)*(-val(icoeff))   + zt(k)
+          else
+            call ijk2idx(icol(icoeff),ix,iy,iz-1,idim,idim,idim)
+            irow(icoeff) = glob_row
+            icoeff       = icoeff+1
+          endif

+          !  term depending on     (x,y,z)
+          val(icoeff)=(2*done)*(a1(x,y,z)+a2(x,y,z)+a3(x,y,z))/sqdeltah &
+               & + c(x,y,z)
+          call ijk2idx(icol(icoeff),ix,iy,iz,idim,idim,idim)
+          irow(icoeff) = glob_row
+          icoeff       = icoeff+1
+          !  term depending on     (x,y,z+1)
+          val(icoeff)=-a3(x,y,z)/sqdeltah+b3(x,y,z)/deltah2
+          if (iz == idim) then
+            zt(k) = g(x,y,done)*(-val(icoeff))   + zt(k)
+          else
+            call ijk2idx(icol(icoeff),ix,iy,iz+1,idim,idim,idim)
+            irow(icoeff) = glob_row
+            icoeff       = icoeff+1
+          endif
+          !  term depending on     (x,y+1,z)
+          val(icoeff)=-a2(x,y,z)/sqdeltah+b2(x,y,z)/deltah2
+          if (iy == idim) then
+            zt(k) = g(x,done,z)*(-val(icoeff))   + zt(k)
+          else
+            call ijk2idx(icol(icoeff),ix,iy+1,iz,idim,idim,idim)
+            irow(icoeff) = glob_row
+            icoeff       = icoeff+1
+          endif
+          !  term depending on     (x+1,y,z)
+          val(icoeff)=-a1(x,y,z)/sqdeltah+b1(x,y,z)/deltah2
+          if (ix==idim) then
+            zt(k) = g(done,y,z)*(-val(icoeff))   + zt(k)
+          else
+            call ijk2idx(icol(icoeff),ix+1,iy,iz,idim,idim,idim)
+            irow(icoeff) = glob_row
+            icoeff       = icoeff+1
+          endif
+
+        end do
+        !write(0,*) ' Outer in_parallel ',omp_in_parallel()
+        call psb_spins(icoeff-1,irow,icol,val,a,desc_a,info)
+        if(info /= psb_success_) cycle
+        call psb_geins(ib,myidx(ii:ii+ib-1),zt(1:ib),bv,desc_a,info)
+        if(info /= psb_success_) cycle
+        zt(:)=dzero
+        call psb_geins(ib,myidx(ii:ii+ib-1),zt(1:ib),xv,desc_a,info)
+        if(info /= psb_success_) cycle
      end do
-      call psb_spins(icoeff-1,irow,icol,val,a,desc_a,info)
-      if(info /= psb_success_) exit
-      call psb_geins(ib,myidx(ii:ii+ib-1),zt(1:ib),bv,desc_a,info)
-      if(info /= psb_success_) exit
-      zt(:)=dzero
-      call psb_geins(ib,myidx(ii:ii+ib-1),zt(1:ib),xv,desc_a,info)
-      if(info /= psb_success_) exit
-    end do
+      !$omp end do
+
+      deallocate(val,irow,icol)
+    end block
+    !$omp end parallel

    tgen = psb_wtime()-t1
    if(info /= psb_success_) then
@@ -490,7 +500,6 @@ contains
      goto 9999
    end if

-    deallocate(val,irow,icol)

    call psb_barrier(ctxt)
    t1 = psb_wtime()
@@ -557,6 +566,9 @@ contains
       & a1,a2,b1,b2,c,g,info,f,amold,vmold,partition, nrl,iv)
    use psb_base_mod
    use psb_util_mod
+#if defined(OPENMP)
+    use omp_lib
+#endif
    !
    !   Discretizes the partial differential equation
    !
@@ -591,7 +603,6 @@ contains
    type(psb_d_csc_sparse_mat)  :: acsc
    type(psb_d_coo_sparse_mat)  :: acoo
    type(psb_d_csr_sparse_mat)  :: acsr
-    real(psb_dpk_)           :: zt(nb),x,y,z,xph,xmh,yph,ymh,zph,zmh
    integer(psb_ipk_) :: nnz,nr,nlr,i,j,ii,ib,k, partition_
    integer(psb_lpk_) :: m,n,glob_row,nt
    integer(psb_ipk_) :: ix,iy,iz,ia,indx_owner
@@ -604,8 +615,7 @@ contains
    ! Process grid
    integer(psb_ipk_) :: np, iam
    integer(psb_ipk_) :: icoeff
-    integer(psb_lpk_), allocatable     :: irow(:),icol(:),myidx(:)
-    real(psb_dpk_), allocatable :: val(:)
+    integer(psb_lpk_), allocatable     :: myidx(:)
    ! deltah dimension of each grid cell
    ! deltat discretization time
    real(psb_dpk_)            :: deltah, sqdeltah, deltah2, dd
@@ -791,7 +801,7 @@ contains
          !write(0,*) iam,' Check on neighbours: ',desc_a%get_p_adjcncy()
        end if
      end block
-      
+
    case default
      write(psb_err_unit,*) iam, 'Initialization error: should not get here'
      info = -1
@@ -816,93 +826,109 @@ contains
      goto 9999
    end if

-    ! we build an auxiliary matrix consisting of one row at a
-    ! time; just a small matrix. might be extended to generate
-    ! a bunch of rows per call.
-    !
-    allocate(val(20*nb),irow(20*nb),&
-         &icol(20*nb),stat=info)
-    if (info /= psb_success_ ) then
-      info=psb_err_alloc_dealloc_
-      call psb_errpush(info,name)
-      goto 9999
-    endif
-
-
-    ! loop over rows belonging to current process in a block
-    ! distribution.
-
    call psb_barrier(ctxt)
    t1 = psb_wtime()
-    do ii=1, nlr,nb
-      ib = min(nb,nlr-ii+1)
-      icoeff = 1
-      do k=1,ib
-        i=ii+k-1
-        ! local matrix pointer
-        glob_row=myidx(i)
-        ! compute gridpoint coordinates
-        call idx2ijk(ix,iy,glob_row,idim,idim)
-        ! x, y coordinates
-        x = (ix-1)*deltah
-        y = (iy-1)*deltah
+    !$omp parallel shared(deltah,myidx,a,desc_a)
+    !
+    block 
+      integer(psb_ipk_) :: i,j,k,ii,ib,icoeff, ix,iy,iz, ith,nth
+      integer(psb_lpk_) :: glob_row
+      integer(psb_lpk_), allocatable :: irow(:),icol(:)
+      real(psb_dpk_), allocatable :: val(:)
+      real(psb_dpk_)     :: x,y,z, zt(nb)
+#if defined(OPENMP)
+      nth = omp_get_num_threads()
+      ith = omp_get_thread_num()
+#else
+      nth = 1
+      ith = 0
+#endif
+      allocate(val(20*nb),irow(20*nb),&
+           &icol(20*nb),stat=info)
+      if (info /= psb_success_ ) then
+        info=psb_err_alloc_dealloc_
+        call psb_errpush(info,name)
+        !goto 9999
+      endif

-        zt(k) = f_(x,y)
-        ! internal point: build discretization
-        !
-        !  term depending on   (x-1,y)
-        !
-        val(icoeff) = -a1(x,y)/sqdeltah-b1(x,y)/deltah2
-        if (ix == 1) then
-          zt(k) = g(dzero,y)*(-val(icoeff)) + zt(k)
-        else
-          call ijk2idx(icol(icoeff),ix-1,iy,idim,idim)
-          irow(icoeff) = glob_row
-          icoeff       = icoeff+1
-        endif
-        !  term depending on     (x,y-1)
-        val(icoeff)  = -a2(x,y)/sqdeltah-b2(x,y)/deltah2
-        if (iy == 1) then
-          zt(k) = g(x,dzero)*(-val(icoeff))   + zt(k)
-        else
-          call ijk2idx(icol(icoeff),ix,iy-1,idim,idim)
-          irow(icoeff) = glob_row
-          icoeff       = icoeff+1
-        endif
+      ! loop over rows belonging to current process in a block
+      ! distribution.
+      !$omp  do schedule(dynamic)
+      !     
+      do ii=1, nlr,nb
+        ib = min(nb,nlr-ii+1)
+        icoeff = 1
+        do k=1,ib
+          i=ii+k-1
+          ! local matrix pointer
+          glob_row=myidx(i)
+          ! compute gridpoint coordinates
+          call idx2ijk(ix,iy,glob_row,idim,idim)
+          ! x, y coordinates
+          x = (ix-1)*deltah
+          y = (iy-1)*deltah

-        !  term depending on     (x,y)
-        val(icoeff)=(2*done)*(a1(x,y) + a2(x,y))/sqdeltah + c(x,y)
-        call ijk2idx(icol(icoeff),ix,iy,idim,idim)
-        irow(icoeff) = glob_row
-        icoeff       = icoeff+1
-        !  term depending on     (x,y+1)
-        val(icoeff)=-a2(x,y)/sqdeltah+b2(x,y)/deltah2
-        if (iy == idim) then
-          zt(k) = g(x,done)*(-val(icoeff))   + zt(k)
-        else
-          call ijk2idx(icol(icoeff),ix,iy+1,idim,idim)
-          irow(icoeff) = glob_row
-          icoeff       = icoeff+1
-        endif
-        !  term depending on     (x+1,y)
-        val(icoeff)=-a1(x,y)/sqdeltah+b1(x,y)/deltah2
-        if (ix==idim) then
-          zt(k) = g(done,y)*(-val(icoeff))   + zt(k)
-        else
-          call ijk2idx(icol(icoeff),ix+1,iy,idim,idim)
-          irow(icoeff) = glob_row
-          icoeff       = icoeff+1
-        endif
+          zt(k) = f_(x,y)
+          ! internal point: build discretization
+          !
+          !  term depending on   (x-1,y)
+          !
+          val(icoeff) = -a1(x,y)/sqdeltah-b1(x,y)/deltah2
+          if (ix == 1) then
+            zt(k) = g(dzero,y)*(-val(icoeff)) + zt(k)
+          else
+            call ijk2idx(icol(icoeff),ix-1,iy,idim,idim)
+            irow(icoeff) = glob_row
+            icoeff       = icoeff+1
+          endif
+          !  term depending on     (x,y-1)
+          val(icoeff)  = -a2(x,y)/sqdeltah-b2(x,y)/deltah2
+          if (iy == 1) then
+            zt(k) = g(x,dzero)*(-val(icoeff))   + zt(k)
+          else
+            call ijk2idx(icol(icoeff),ix,iy-1,idim,idim)
+            irow(icoeff) = glob_row
+            icoeff       = icoeff+1
+          endif

+          !  term depending on     (x,y)
+          val(icoeff)=(2*done)*(a1(x,y) + a2(x,y))/sqdeltah + c(x,y)
+          call ijk2idx(icol(icoeff),ix,iy,idim,idim)
+          irow(icoeff) = glob_row
+          icoeff       = icoeff+1
+          !  term depending on     (x,y+1)
+          val(icoeff)=-a2(x,y)/sqdeltah+b2(x,y)/deltah2
+          if (iy == idim) then
+            zt(k) = g(x,done)*(-val(icoeff))   + zt(k)
+          else
+            call ijk2idx(icol(icoeff),ix,iy+1,idim,idim)
+            irow(icoeff) = glob_row
+            icoeff       = icoeff+1
+          endif
+          !  term depending on     (x+1,y)
+          val(icoeff)=-a1(x,y)/sqdeltah+b1(x,y)/deltah2
+          if (ix==idim) then
+            zt(k) = g(done,y)*(-val(icoeff))   + zt(k)
+          else
+            call ijk2idx(icol(icoeff),ix+1,iy,idim,idim)
+            irow(icoeff) = glob_row
+            icoeff       = icoeff+1
+          endif
+
+        end do
+        call psb_spins(icoeff-1,irow,icol,val,a,desc_a,info)
+        if(info /= psb_success_) cycle
+        call psb_geins(ib,myidx(ii:ii+ib-1),zt(1:ib),bv,desc_a,info)
+        if(info /= psb_success_) cycle
+        zt(:)=dzero
+        call psb_geins(ib,myidx(ii:ii+ib-1),zt(1:ib),xv,desc_a,info)
+        if(info /= psb_success_) cycle
      end do
-      call psb_spins(icoeff-1,irow,icol,val,a,desc_a,info)
-      if(info /= psb_success_) exit
-      call psb_geins(ib,myidx(ii:ii+ib-1),zt(1:ib),bv,desc_a,info)
-      if(info /= psb_success_) exit
-      zt(:)=dzero
-      call psb_geins(ib,myidx(ii:ii+ib-1),zt(1:ib),xv,desc_a,info)
-      if(info /= psb_success_) exit
-    end do
+      !$omp end do
+
+      deallocate(val,irow,icol)
+    end block
+    !$omp end parallel

    tgen = psb_wtime()-t1
    if(info /= psb_success_) then
@@ -912,8 +938,6 @@ contains
      goto 9999
    end if

-    deallocate(val,irow,icol)
-
    call psb_barrier(ctxt)
    t1 = psb_wtime()
    call psb_cdasb(desc_a,info)
@@ -93,6 +93,9 @@ contains
       & a1,a2,a3,b1,b2,b3,c,g,info,f,amold,vmold,partition, nrl,iv)
    use psb_base_mod
    use psb_util_mod
+#if defined(OPENMP)
+    use omp_lib
+#endif
    !
    !   Discretizes the partial differential equation
    !
@@ -128,7 +131,6 @@ contains
    type(psb_s_csc_sparse_mat)  :: acsc
    type(psb_s_coo_sparse_mat)  :: acoo
    type(psb_s_csr_sparse_mat)  :: acsr
-    real(psb_spk_)           :: zt(nb),x,y,z,xph,xmh,yph,ymh,zph,zmh
    integer(psb_ipk_) :: nnz,nr,nlr,i,j,ii,ib,k, partition_
    integer(psb_lpk_) :: m,n,glob_row,nt
    integer(psb_ipk_) :: ix,iy,iz,ia,indx_owner
@@ -141,8 +143,7 @@ contains
    ! Process grid
    integer(psb_ipk_) :: np, iam
    integer(psb_ipk_) :: icoeff
-    integer(psb_lpk_), allocatable     :: irow(:),icol(:),myidx(:)
-    real(psb_spk_), allocatable :: val(:)
+    integer(psb_lpk_), allocatable     :: myidx(:)
    ! deltah dimension of each grid cell
    ! deltat discretization time
    real(psb_spk_)            :: deltah, sqdeltah, deltah2
@@ -368,119 +369,128 @@ contains
    call psb_barrier(ctxt)
    talc = psb_wtime()-t0

-    if (info /= psb_success_) then
-      info=psb_err_from_subroutine_
-      ch_err='allocation rout.'
-      call psb_errpush(info,name,a_err=ch_err)
-      goto 9999
-    end if
-
-    ! we build an auxiliary matrix consisting of one row at a
-    ! time; just a small matrix. might be extended to generate
-    ! a bunch of rows per call.
-    !
-    allocate(val(20*nb),irow(20*nb),&
-         &icol(20*nb),stat=info)
-    if (info /= psb_success_ ) then
-      info=psb_err_alloc_dealloc_
-      call psb_errpush(info,name)
-      goto 9999
-    endif
-
-
-    ! loop over rows belonging to current process in a block
-    ! distribution.
-
    call psb_barrier(ctxt)
    t1 = psb_wtime()
-    do ii=1, nlr,nb
-      ib = min(nb,nlr-ii+1)
-      icoeff = 1
-      do k=1,ib
-        i=ii+k-1
-        ! local matrix pointer
-        glob_row=myidx(i)
-        ! compute gridpoint coordinates
-        call idx2ijk(ix,iy,iz,glob_row,idim,idim,idim)
-        ! x, y, z coordinates
-        x = (ix-1)*deltah
-        y = (iy-1)*deltah
-        z = (iz-1)*deltah
-        zt(k) = f_(x,y,z)
-        ! internal point: build discretization
-        !
-        !  term depending on   (x-1,y,z)
-        !
-        val(icoeff) = -a1(x,y,z)/sqdeltah-b1(x,y,z)/deltah2
-        if (ix == 1) then
-          zt(k) = g(szero,y,z)*(-val(icoeff)) + zt(k)
-        else
-          call ijk2idx(icol(icoeff),ix-1,iy,iz,idim,idim,idim)
-          irow(icoeff) = glob_row
-          icoeff       = icoeff+1
-        endif
-        !  term depending on     (x,y-1,z)
-        val(icoeff)  = -a2(x,y,z)/sqdeltah-b2(x,y,z)/deltah2
-        if (iy == 1) then
-          zt(k) = g(x,szero,z)*(-val(icoeff))   + zt(k)
-        else
-          call ijk2idx(icol(icoeff),ix,iy-1,iz,idim,idim,idim)
-          irow(icoeff) = glob_row
-          icoeff       = icoeff+1
-        endif
-        !  term depending on     (x,y,z-1)
-        val(icoeff)=-a3(x,y,z)/sqdeltah-b3(x,y,z)/deltah2
-        if (iz == 1) then
-          zt(k) = g(x,y,szero)*(-val(icoeff))   + zt(k)
-        else
-          call ijk2idx(icol(icoeff),ix,iy,iz-1,idim,idim,idim)
-          irow(icoeff) = glob_row
-          icoeff       = icoeff+1
-        endif
+    !$omp parallel shared(deltah,myidx,a,desc_a)
+    !
+    block 
+      integer(psb_ipk_) :: i,j,k,ii,ib,icoeff, ix,iy,iz, ith,nth
+      integer(psb_lpk_) :: glob_row
+      integer(psb_lpk_), allocatable :: irow(:),icol(:)
+      real(psb_spk_), allocatable :: val(:)
+      real(psb_spk_)     :: x,y,z, zt(nb)
+#if defined(OPENMP)
+      nth = omp_get_num_threads()
+      ith = omp_get_thread_num()
+#else
+      nth = 1
+      ith = 0
+#endif
+      allocate(val(20*nb),irow(20*nb),&
+           &icol(20*nb),stat=info)
+      if (info /= psb_success_ ) then
+        info=psb_err_alloc_dealloc_
+        call psb_errpush(info,name)
+        !goto 9999
+      endif

-        !  term depending on     (x,y,z)
-        val(icoeff)=(2*sone)*(a1(x,y,z)+a2(x,y,z)+a3(x,y,z))/sqdeltah &
-             & + c(x,y,z)
-        call ijk2idx(icol(icoeff),ix,iy,iz,idim,idim,idim)
-        irow(icoeff) = glob_row
-        icoeff       = icoeff+1
-        !  term depending on     (x,y,z+1)
-        val(icoeff)=-a3(x,y,z)/sqdeltah+b3(x,y,z)/deltah2
-        if (iz == idim) then
-          zt(k) = g(x,y,sone)*(-val(icoeff))   + zt(k)
-        else
-          call ijk2idx(icol(icoeff),ix,iy,iz+1,idim,idim,idim)
-          irow(icoeff) = glob_row
-          icoeff       = icoeff+1
-        endif
-        !  term depending on     (x,y+1,z)
-        val(icoeff)=-a2(x,y,z)/sqdeltah+b2(x,y,z)/deltah2
-        if (iy == idim) then
-          zt(k) = g(x,sone,z)*(-val(icoeff))   + zt(k)
-        else
-          call ijk2idx(icol(icoeff),ix,iy+1,iz,idim,idim,idim)
-          irow(icoeff) = glob_row
-          icoeff       = icoeff+1
-        endif
-        !  term depending on     (x+1,y,z)
-        val(icoeff)=-a1(x,y,z)/sqdeltah+b1(x,y,z)/deltah2
-        if (ix==idim) then
-          zt(k) = g(sone,y,z)*(-val(icoeff))   + zt(k)
-        else
-          call ijk2idx(icol(icoeff),ix+1,iy,iz,idim,idim,idim)
-          irow(icoeff) = glob_row
-          icoeff       = icoeff+1
-        endif
+      !$omp  do schedule(dynamic)
+      !     
+      do ii=1, nlr, nb
+        if (info /= psb_success_) cycle
+        ib = min(nb,nlr-ii+1)
+        icoeff = 1
+        do k=1,ib
+          i=ii+k-1
+          ! local matrix pointer
+          glob_row=myidx(i)
+          ! compute gridpoint coordinates
+          call idx2ijk(ix,iy,iz,glob_row,idim,idim,idim)
+          ! x, y, z coordinates
+          x = (ix-1)*deltah
+          y = (iy-1)*deltah
+          z = (iz-1)*deltah
+          zt(k) = f_(x,y,z)
+          ! internal point: build discretization
+          !
+          !  term depending on   (x-1,y,z)
+          !
+          val(icoeff) = -a1(x,y,z)/sqdeltah-b1(x,y,z)/deltah2
+          if (ix == 1) then
+            zt(k) = g(szero,y,z)*(-val(icoeff)) + zt(k)
+          else
+            call ijk2idx(icol(icoeff),ix-1,iy,iz,idim,idim,idim)
+            irow(icoeff) = glob_row
+            icoeff       = icoeff+1
+          endif
+          !  term depending on     (x,y-1,z)
+          val(icoeff)  = -a2(x,y,z)/sqdeltah-b2(x,y,z)/deltah2
+          if (iy == 1) then
+            zt(k) = g(x,szero,z)*(-val(icoeff))   + zt(k)
+          else
+            call ijk2idx(icol(icoeff),ix,iy-1,iz,idim,idim,idim)
+            irow(icoeff) = glob_row
+            icoeff       = icoeff+1
+          endif
+          !  term depending on     (x,y,z-1)
+          val(icoeff)=-a3(x,y,z)/sqdeltah-b3(x,y,z)/deltah2
+          if (iz == 1) then
+            zt(k) = g(x,y,szero)*(-val(icoeff))   + zt(k)
+          else
+            call ijk2idx(icol(icoeff),ix,iy,iz-1,idim,idim,idim)
+            irow(icoeff) = glob_row
+            icoeff       = icoeff+1
+          endif

+          !  term depending on     (x,y,z)
+          val(icoeff)=(2*sone)*(a1(x,y,z)+a2(x,y,z)+a3(x,y,z))/sqdeltah &
+               & + c(x,y,z)
+          call ijk2idx(icol(icoeff),ix,iy,iz,idim,idim,idim)
+          irow(icoeff) = glob_row
+          icoeff       = icoeff+1
+          !  term depending on     (x,y,z+1)
+          val(icoeff)=-a3(x,y,z)/sqdeltah+b3(x,y,z)/deltah2
+          if (iz == idim) then
+            zt(k) = g(x,y,sone)*(-val(icoeff))   + zt(k)
+          else
+            call ijk2idx(icol(icoeff),ix,iy,iz+1,idim,idim,idim)
+            irow(icoeff) = glob_row
+            icoeff       = icoeff+1
+          endif
+          !  term depending on     (x,y+1,z)
+          val(icoeff)=-a2(x,y,z)/sqdeltah+b2(x,y,z)/deltah2
+          if (iy == idim) then
+            zt(k) = g(x,sone,z)*(-val(icoeff))   + zt(k)
+          else
+            call ijk2idx(icol(icoeff),ix,iy+1,iz,idim,idim,idim)
+            irow(icoeff) = glob_row
+            icoeff       = icoeff+1
+          endif
+          !  term depending on     (x+1,y,z)
+          val(icoeff)=-a1(x,y,z)/sqdeltah+b1(x,y,z)/deltah2
+          if (ix==idim) then
+            zt(k) = g(sone,y,z)*(-val(icoeff))   + zt(k)
+          else
+            call ijk2idx(icol(icoeff),ix+1,iy,iz,idim,idim,idim)
+            irow(icoeff) = glob_row
+            icoeff       = icoeff+1
+          endif
+
+        end do
+        !write(0,*) ' Outer in_parallel ',omp_in_parallel()
+        call psb_spins(icoeff-1,irow,icol,val,a,desc_a,info)
+        if(info /= psb_success_) cycle
+        call psb_geins(ib,myidx(ii:ii+ib-1),zt(1:ib),bv,desc_a,info)
+        if(info /= psb_success_) cycle
+        zt(:)=szero
+        call psb_geins(ib,myidx(ii:ii+ib-1),zt(1:ib),xv,desc_a,info)
+        if(info /= psb_success_) cycle
      end do
-      call psb_spins(icoeff-1,irow,icol,val,a,desc_a,info)
-      if(info /= psb_success_) exit
-      call psb_geins(ib,myidx(ii:ii+ib-1),zt(1:ib),bv,desc_a,info)
-      if(info /= psb_success_) exit
-      zt(:)=szero
-      call psb_geins(ib,myidx(ii:ii+ib-1),zt(1:ib),xv,desc_a,info)
-      if(info /= psb_success_) exit
-    end do
+      !$omp end do
+
+      deallocate(val,irow,icol)
+    end block
+    !$omp end parallel

    tgen = psb_wtime()-t1
    if(info /= psb_success_) then
@@ -490,7 +500,6 @@ contains
      goto 9999
    end if

-    deallocate(val,irow,icol)

    call psb_barrier(ctxt)
    t1 = psb_wtime()
@@ -557,6 +566,9 @@ contains
       & a1,a2,b1,b2,c,g,info,f,amold,vmold,partition, nrl,iv)
    use psb_base_mod
    use psb_util_mod
+#if defined(OPENMP)
+    use omp_lib
+#endif
    !
    !   Discretizes the partial differential equation
    !
@@ -591,7 +603,6 @@ contains
    type(psb_s_csc_sparse_mat)  :: acsc
    type(psb_s_coo_sparse_mat)  :: acoo
    type(psb_s_csr_sparse_mat)  :: acsr
-    real(psb_spk_)           :: zt(nb),x,y,z,xph,xmh,yph,ymh,zph,zmh
    integer(psb_ipk_) :: nnz,nr,nlr,i,j,ii,ib,k, partition_
    integer(psb_lpk_) :: m,n,glob_row,nt
    integer(psb_ipk_) :: ix,iy,iz,ia,indx_owner
@@ -604,8 +615,7 @@ contains
    ! Process grid
    integer(psb_ipk_) :: np, iam
    integer(psb_ipk_) :: icoeff
-    integer(psb_lpk_), allocatable     :: irow(:),icol(:),myidx(:)
-    real(psb_spk_), allocatable :: val(:)
+    integer(psb_lpk_), allocatable     :: myidx(:)
    ! deltah dimension of each grid cell
    ! deltat discretization time
    real(psb_spk_)            :: deltah, sqdeltah, deltah2, dd
@@ -791,7 +801,7 @@ contains
          !write(0,*) iam,' Check on neighbours: ',desc_a%get_p_adjcncy()
        end if
      end block
-      
+
    case default
      write(psb_err_unit,*) iam, 'Initialization error: should not get here'
      info = -1
@@ -816,93 +826,109 @@ contains
      goto 9999
    end if

-    ! we build an auxiliary matrix consisting of one row at a
-    ! time; just a small matrix. might be extended to generate
-    ! a bunch of rows per call.
-    !
-    allocate(val(20*nb),irow(20*nb),&
-         &icol(20*nb),stat=info)
-    if (info /= psb_success_ ) then
-      info=psb_err_alloc_dealloc_
-      call psb_errpush(info,name)
-      goto 9999
-    endif
-
-
-    ! loop over rows belonging to current process in a block
-    ! distribution.
-
    call psb_barrier(ctxt)
    t1 = psb_wtime()
-    do ii=1, nlr,nb
-      ib = min(nb,nlr-ii+1)
-      icoeff = 1
-      do k=1,ib
-        i=ii+k-1
-        ! local matrix pointer
-        glob_row=myidx(i)
-        ! compute gridpoint coordinates
-        call idx2ijk(ix,iy,glob_row,idim,idim)
-        ! x, y coordinates
-        x = (ix-1)*deltah
-        y = (iy-1)*deltah
+    !$omp parallel shared(deltah,myidx,a,desc_a)
+    !
+    block 
+      integer(psb_ipk_) :: i,j,k,ii,ib,icoeff, ix,iy,iz, ith,nth
+      integer(psb_lpk_) :: glob_row
+      integer(psb_lpk_), allocatable :: irow(:),icol(:)
+      real(psb_spk_), allocatable :: val(:)
+      real(psb_spk_)     :: x,y,z, zt(nb)
+#if defined(OPENMP)
+      nth = omp_get_num_threads()
+      ith = omp_get_thread_num()
+#else
+      nth = 1
+      ith = 0
+#endif
+      allocate(val(20*nb),irow(20*nb),&
+           &icol(20*nb),stat=info)
+      if (info /= psb_success_ ) then
+        info=psb_err_alloc_dealloc_
+        call psb_errpush(info,name)
+        !goto 9999
+      endif

-        zt(k) = f_(x,y)
-        ! internal point: build discretization
-        !
-        !  term depending on   (x-1,y)
-        !
-        val(icoeff) = -a1(x,y)/sqdeltah-b1(x,y)/deltah2
-        if (ix == 1) then
-          zt(k) = g(szero,y)*(-val(icoeff)) + zt(k)
-        else
-          call ijk2idx(icol(icoeff),ix-1,iy,idim,idim)
-          irow(icoeff) = glob_row
-          icoeff       = icoeff+1
-        endif
-        !  term depending on     (x,y-1)
-        val(icoeff)  = -a2(x,y)/sqdeltah-b2(x,y)/deltah2
-        if (iy == 1) then
-          zt(k) = g(x,szero)*(-val(icoeff))   + zt(k)
-        else
-          call ijk2idx(icol(icoeff),ix,iy-1,idim,idim)
-          irow(icoeff) = glob_row
-          icoeff       = icoeff+1
-        endif
+      ! loop over rows belonging to current process in a block
+      ! distribution.
+      !$omp  do schedule(dynamic)
+      !     
+      do ii=1, nlr,nb
+        ib = min(nb,nlr-ii+1)
+        icoeff = 1
+        do k=1,ib
+          i=ii+k-1
+          ! local matrix pointer
+          glob_row=myidx(i)
+          ! compute gridpoint coordinates
+          call idx2ijk(ix,iy,glob_row,idim,idim)
+          ! x, y coordinates
+          x = (ix-1)*deltah
+          y = (iy-1)*deltah

-        !  term depending on     (x,y)
-        val(icoeff)=(2*sone)*(a1(x,y) + a2(x,y))/sqdeltah + c(x,y)
-        call ijk2idx(icol(icoeff),ix,iy,idim,idim)
-        irow(icoeff) = glob_row
-        icoeff       = icoeff+1
-        !  term depending on     (x,y+1)
-        val(icoeff)=-a2(x,y)/sqdeltah+b2(x,y)/deltah2
-        if (iy == idim) then
-          zt(k) = g(x,sone)*(-val(icoeff))   + zt(k)
-        else
-          call ijk2idx(icol(icoeff),ix,iy+1,idim,idim)
-          irow(icoeff) = glob_row
-          icoeff       = icoeff+1
-        endif
-        !  term depending on     (x+1,y)
-        val(icoeff)=-a1(x,y)/sqdeltah+b1(x,y)/deltah2
-        if (ix==idim) then
-          zt(k) = g(sone,y)*(-val(icoeff))   + zt(k)
-        else
-          call ijk2idx(icol(icoeff),ix+1,iy,idim,idim)
-          irow(icoeff) = glob_row
-          icoeff       = icoeff+1
-        endif
+          zt(k) = f_(x,y)
+          ! internal point: build discretization
+          !
+          !  term depending on   (x-1,y)
+          !
+          val(icoeff) = -a1(x,y)/sqdeltah-b1(x,y)/deltah2
+          if (ix == 1) then
+            zt(k) = g(szero,y)*(-val(icoeff)) + zt(k)
+          else
+            call ijk2idx(icol(icoeff),ix-1,iy,idim,idim)
+            irow(icoeff) = glob_row
+            icoeff       = icoeff+1
+          endif
+          !  term depending on     (x,y-1)
+          val(icoeff)  = -a2(x,y)/sqdeltah-b2(x,y)/deltah2
+          if (iy == 1) then
+            zt(k) = g(x,szero)*(-val(icoeff))   + zt(k)
+          else
+            call ijk2idx(icol(icoeff),ix,iy-1,idim,idim)
+            irow(icoeff) = glob_row
+            icoeff       = icoeff+1
+          endif

+          !  term depending on     (x,y)
+          val(icoeff)=(2*sone)*(a1(x,y) + a2(x,y))/sqdeltah + c(x,y)
+          call ijk2idx(icol(icoeff),ix,iy,idim,idim)
+          irow(icoeff) = glob_row
+          icoeff       = icoeff+1
+          !  term depending on     (x,y+1)
+          val(icoeff)=-a2(x,y)/sqdeltah+b2(x,y)/deltah2
+          if (iy == idim) then
+            zt(k) = g(x,sone)*(-val(icoeff))   + zt(k)
+          else
+            call ijk2idx(icol(icoeff),ix,iy+1,idim,idim)
+            irow(icoeff) = glob_row
+            icoeff       = icoeff+1
+          endif
+          !  term depending on     (x+1,y)
+          val(icoeff)=-a1(x,y)/sqdeltah+b1(x,y)/deltah2
+          if (ix==idim) then
+            zt(k) = g(sone,y)*(-val(icoeff))   + zt(k)
+          else
+            call ijk2idx(icol(icoeff),ix+1,iy,idim,idim)
+            irow(icoeff) = glob_row
+            icoeff       = icoeff+1
+          endif
+
+        end do
+        call psb_spins(icoeff-1,irow,icol,val,a,desc_a,info)
+        if(info /= psb_success_) cycle
+        call psb_geins(ib,myidx(ii:ii+ib-1),zt(1:ib),bv,desc_a,info)
+        if(info /= psb_success_) cycle
+        zt(:)=szero
+        call psb_geins(ib,myidx(ii:ii+ib-1),zt(1:ib),xv,desc_a,info)
+        if(info /= psb_success_) cycle
      end do
-      call psb_spins(icoeff-1,irow,icol,val,a,desc_a,info)
-      if(info /= psb_success_) exit
-      call psb_geins(ib,myidx(ii:ii+ib-1),zt(1:ib),bv,desc_a,info)
-      if(info /= psb_success_) exit
-      zt(:)=szero
-      call psb_geins(ib,myidx(ii:ii+ib-1),zt(1:ib),xv,desc_a,info)
-      if(info /= psb_success_) exit
-    end do
+      !$omp end do
+
+      deallocate(val,irow,icol)
+    end block
+    !$omp end parallel

    tgen = psb_wtime()-t1
    if(info /= psb_success_) then
@@ -912,8 +938,6 @@ contains
      goto 9999
    end if

-    deallocate(val,irow,icol)
-
    call psb_barrier(ctxt)
    t1 = psb_wtime()
    call psb_cdasb(desc_a,info)
@@ -1,6 +1,6 @@
 %%%%%%%%%%%  General  arguments % Lines starting with % are ignored.
 CSR                         ! Storage format CSR COO JAD
-0080                        ! IDIM; domain size. Linear system size is IDIM**3
+0200                       ! IDIM; domain size. Linear system size is IDIM**3
 CONST                       ! PDECOEFF: CONST, EXP, GAUSS Coefficients of the PDE
 BICGSTAB                    ! Iterative method: BiCGSTAB BiCGSTABL BiCG CG CGS FCG GCR RGMRES
 2                           ! ISTOPC
Author	SHA1	Message	Date
sfilippone	494b8b925f	OpenMP loop in samples data generation	2023-06-05 11:46:39 +02:00
sfilippone	73e5d49913	Added timers to build phases	2023-06-02 11:37:58 +02:00
Salvatore Filippone	a612cea167	Debug for matchboxp	2023-02-10 07:53:04 -05:00
Salvatore Filippone	ebe9b45177	Modify MATCHBOXP to fix OpenMP. Performance to be reviewed	2023-02-10 07:50:58 -05:00
Salvatore Filippone	32994c7ce8	Better parameters in matchboxp_mod	2022-12-13 06:44:12 -05:00
Salvatore Filippone	d59c9e6c0a	Updates towards OpenMP version.	2022-11-22 03:02:51 -05:00
Salvatore Filippone	0d624df346	Merge branch 'petrilli-m' into openmp-match	2022-11-14 04:49:50 -05:00
StefanoPetrilli	6414d3aef3	U and privateU are now vectors	2022-07-23 12:47:43 -05:00
StefanoPetrilli	a259e8ab53	extractUChunch optimization	2022-07-23 11:34:43 -05:00
StefanoPetrilli	500403dbda	Replaced some staticQueues with vectors for performance reasons	2022-07-23 11:13:21 -05:00
StefanoPetrilli	066c1a5e62	optimization processMatchedVerticesAndSendMessages.cpp	2022-07-23 09:27:35 -05:00
StefanoPetrilli	1ab166b38b	Improved performance of processMatchedVerticesAndSendMessages.cpp	2022-07-23 08:24:50 -05:00
StefanoPetrilli	5efee20041	Optimization, replaced all useless atomic with reduction	2022-07-23 05:52:27 -05:00
StefanoPetrilli	aa45e2fe93	processMatchedVerticesAndSendMessages.cpp unoptimized	2022-07-23 05:14:26 -05:00
StefanoPetrilli	e328f3969c	queueTransfer optimization in processMatchedVertices	2022-07-22 07:25:09 -05:00
StefanoPetrilli	9d1a416f99	add rm to exec.sh	2022-07-21 15:45:31 -05:00
StefanoPetrilli	9b065602a8	Fixed race condition in processExposedVertices	2022-07-20 16:24:37 -05:00
StefanoPetrilli	abf258e2e8	isAlreadyMatched is now atomic	2022-07-20 15:45:29 -05:00
StefanoPetrilli	cdf92ea2b2	processMatchedVerticess add send messages with error	2022-07-20 15:37:29 -05:00
StefanoPetrilli	22d9baf296	isAlreadyMatched substituted with atomic read in one place	2022-07-18 14:11:53 -05:00
StefanoPetrilli	44f174a571	findOwnerOfGhost optimization and refactor	2022-07-17 13:44:58 -05:00
StefanoPetrilli	3e945c75b4	Refactoring, removed all useless Pointer passed in functions	2022-07-17 13:20:49 -05:00
StefanoPetrilli	a71fe82752	PROCESS_CROSS_EDGE refactoring	2022-07-17 12:03:48 -05:00
StefanoPetrilli	4f07a70ed1	initialize refactoring	2022-07-17 11:48:52 -05:00
StefanoPetrilli	cb660e044d	Remoe MateLock	2022-07-17 11:27:17 -05:00
StefanoPetrilli	d24c8c2d46	processCrossEdges is now atomic	2022-07-17 09:43:48 -05:00
StefanoPetrilli	9ab54adf3f	processMatchedVertices parallelized	2022-07-17 08:59:23 -05:00
StefanoPetrilli	71d4cdc319	processMatchedVertices rollback to critical regions	2022-07-17 06:11:11 -05:00
StefanoPetrilli	1374f21ba8	refactor increment on variables passed by reference in processMatchedVertices.cpp	2022-07-16 13:54:40 -05:00
StefanoPetrilli	a9bb6b26fa	processMatchedVertices partially working mixed critical and lock version	2022-07-16 11:20:39 -05:00
StefanoPetrilli	561cadee0f	parallelQueues working	2022-07-15 07:27:30 -05:00
StefanoPetrilli	5ca78fb871	Refactoring isAlreadyMatched and processCrossEdge	2022-07-14 17:10:36 -05:00
StefanoPetrilli	f17082b337	Refactoring: eliminatino of SPtr inside processMessages	2022-07-14 15:27:53 -05:00
StefanoPetrilli	1ea1be33ba	Refactoring, eliminated useless passed variables	2022-07-14 15:23:32 -05:00
StefanoPetrilli	47c6f4f2f8	comments	2022-07-13 16:19:52 -05:00
StefanoPetrilli	dc1675766f	processMessages.cpp further refactoring	2022-07-13 16:19:38 -05:00
StefanoPetrilli	ccac816f52	processCrossEdge small refactoring	2022-07-12 13:24:12 -05:00
StefanoPetrilli	c7e8193514	omp task in clean.cpp, lock destroy	2022-07-12 12:12:15 -05:00
StefanoPetrilli	36bd3a51a2	Makefile fix	2022-07-11 16:31:58 -05:00
StefanoPetrilli	32777cc15c	clean partial refactoring	2022-07-10 11:09:10 -05:00
StefanoPetrilli	64c23f93f8	processMessags partial refactoring, message const refactoring	2022-07-10 10:01:50 -05:00
StefanoPetrilli	d19443052d	Insert private queue error in processMatchedVertices.cpp	2022-07-10 05:24:31 -05:00
StefanoPetrilli	df1e4a4616	PROCESS_CROSS_EDGE refactoring	2022-07-10 04:31:51 -05:00
StefanoPetrilli	3de1e607eb	sendBundledMessages refactoring	2022-07-10 03:39:58 -05:00
StefanoPetrilli	9b13aef1ce	processMathedVertices refactoring	2022-07-08 13:32:24 -05:00
StefanoPetrilli	6dcae6d0c1	fix private queues in PARALLEL_PROCESS_EXPOSED_VERTEX_B	2022-07-06 15:33:29 -05:00
StefanoPetrilli	63b7602d3a	refactoring queueTransfer	2022-07-06 13:12:31 -05:00
StefanoPetrilli	b66de7f25c	Refactoring PARALLEL_PROCESS_EXPOSED_VERTEX_B	2022-07-06 12:58:00 -05:00
StefanoPetrilli	46047b2202	refactoring parallelComputeCandidateMateB	2022-06-30 16:48:18 -05:00
StefanoPetrilli	7cfe198d0f	Format	2022-06-26 10:45:06 -05:00
StefanoPetrilli	1aca17cd44	initialize fix	2022-06-26 10:02:11 -05:00
StefanoPetrilli	ea040ae5ee	Reformat initialize, refactoring of initialize completed	2022-06-26 04:48:49 -05:00
StefanoPetrilli	7741abd45d	Initialize parallelized with task	2022-06-26 04:40:13 -05:00
StefanoPetrilli	b5e52d31f5	Refactoring private queues, still not working	2022-06-25 15:25:13 -05:00
StefanoPetrilli	deab695294	Refactoring Initialization	2022-06-25 12:10:14 -05:00
StefanoPetrilli	a54f084ffb	refactoring, initialization	2022-06-25 10:16:30 -05:00
StefanoPetrilli	bf0532867d	Functions in different files	2022-06-25 08:48:49 -05:00
StefanoPetrilli	2044c5c8eb	Merge fix, lock error	2022-06-14 14:47:45 -05:00
StefanoPetrilli	f38f3cf09a	Merge branch 'tmp' into ompmpi_aggregator_stefano_petrilli	2022-06-14 14:44:24 -05:00
StefanoPetrilli	6fd571ecb2	Lock error	2022-06-14 14:33:31 -05:00
StefanoPetrilli	bf35c1659b	Further improved critical region U	2022-06-13 16:53:12 -05:00
StefanoPetrilli	b2230a6d6d	Improved critical region U	2022-06-13 16:09:00 -05:00
StefanoPetrilli	6c20cd7819	PROCESS MATCHED VERTICES draft of parallelization	2022-06-10 15:34:29 -05:00
StefanoPetrilli	f921aa47c4	Master region for tempCounter.clear() (Might have solved stucked runs)	2022-06-08 15:19:57 -05:00
StefanoPetrilli	532701031e	Extendend parallel region after SEND PACKET BUNDLE Nothing parallelizable founded	2022-06-02 09:15:31 -05:00
StefanoPetrilli	b079d71f30	Further optimizations PARALLEL_PROCESS_EXPOSED_VERTEX_B	2022-06-02 07:29:21 -05:00
StefanoPetrilli	e2ca97ca47	Removed one critical region from PARALLEL_PROCESS_EXPOSED_VERTEX_B	2022-05-31 16:04:56 -05:00
StefanoPetrilli	5bc4f2a080	PROCESS MATCHED VERTICES parallelization improvement	2022-05-30 14:27:26 -05:00
StefanoPetrilli	2c8dc2ffdd	PROCESS MATCHED VERTICES parallelization draft	2022-05-30 13:49:34 -05:00
StefanoPetrilli	f3d7b3ab5e	False sharing fix	2022-05-29 12:01:28 -05:00
StefanoPetrilli	766ef320c2	Refactoring + critical(Mate)	2022-05-29 12:01:24 -05:00
StefanoPetrilli	002239f5b6	False sharing fix	2022-05-22 17:35:08 -05:00
StefanoPetrilli	70b7c4db55	PARALLEL_PROCESS_EXPOSED_VERTEX_B named critical sections	2022-05-22 16:50:07 -05:00
StefanoPetrilli	2cac21b345	fix and reformatting	2022-05-21 11:46:40 -05:00
StefanoPetrilli	6180f29f39	PARALLEL_COMPUTE_CANDIDATE_MATE_B is now paralle and correct	2022-05-21 11:23:39 -05:00
StefanoPetrilli	b4bfdd83e5	computeCandidateMate and isAlreadyMatched	2022-05-21 10:22:58 -05:00
StefanoPetrilli	1140669ea7	firstComputeCandidateMate	2022-05-21 07:01:42 -05:00
StefanoPetrilli	919e2a2918	PARALLEL_PROCESS_EXPOSED_VERTEX_B is actually not parallelizable. Atleast not as I was doing.	2022-05-21 05:56:05 -05:00
StefanoPetrilli	baffff3d93	Instable PARALLEL_PROCESS_EXPOSED_VERTEX_B	2022-05-09 16:52:03 -05:00
StefanoPetrilli	25a603debe	PARALLEL_COMPUTE_CANDIDATE_MATE_B	2022-05-08 15:11:56 -05:00
StefanoPetrilli	a20f0d47e7	Solved the static queue out of scope problem	2022-05-08 12:15:46 -05:00
StefanoPetrilli	76e04ee997	The OMP and MPI version is now separated in two different files	2022-05-05 15:57:58 -05:00
StefanoPetrilli	0a8debe43a	Single parallel regions with multiple for cycles Added OMP for testing	2022-05-01 15:26:47 -05:00
StefanoPetrilli	8f6dc5fac2	verGhostPtrInitialization is now parallelized	2022-05-01 06:05:16 -05:00
StefanoPetrilli	7d40fde21d	verGhostIndInitialization and Ghost2LocalInitialization cycles parallelization	2022-05-01 05:42:42 -05:00
StefanoPetrilli	1760afbe97	Time tracking in algoDistEdge	2022-05-01 04:47:03 -05:00
StefanoPetrilli	60f90804d5	Time tracking in MatchBox	2022-05-01 04:42:33 -05:00