Get rid of empty aggregates at the end of PMIS algorithm

amgcl/mpi/coarsening/pmis.hpp

@@ -434,7 +434,9 @@ struct pmis {
         // 2. Apply PMIS algorithm to the symbolic square.
         ptrdiff_t n_undone = 0;
         std::vector<ptrdiff_t> rem_state(Sp.recv.count(), undone);
+        std::vector<int>       rem_owner(Sp.recv.count(), -1);
         std::vector<ptrdiff_t> send_state(Sp.send.count());
+        std::vector<int>       send_owner(Sp.send.count());
 
         // Remove lonely nodes.
 #pragma omp parallel for reduction(+:n_undone)
@@ -513,6 +515,7 @@ struct pmis {
                             std::tie(d,k) = Sp.remote_info(c);
 
                             rem_state[k] = id;
+
                             send_pts[d].push_back(c);
                             send_pts[d].push_back(id);
                         }
@@ -606,14 +609,85 @@ struct pmis {
                 MPI_Wait(&send_pts_req[i], MPI_STATUS_IGNORE);
             }
 
-
             for(ptrdiff_t i = 0, m = Sp.send.count(); i < m; ++i)
                 send_state[i] = loc_state[Sp.send.col[i]];
             Sp.exchange(&send_state[0], &rem_state[0]);
 
             if (0 == comm.reduce(MPI_SUM, n_undone))
                 break;
         }
+
+        // Some of the aggregates could potentially vanish during expansion
+        // step (*) above. We need to exclude those and renumber the rest.
+        AMGCL_TIC("drop empty aggregates");
+        for(ptrdiff_t i = 0, m = Sp.send.count(); i < m; ++i)
+            send_owner[i] = loc_owner[Sp.send.col[i]];
+        Sp.exchange(&send_owner[0], &rem_owner[0]);
+
+        std::vector<ptrdiff_t> new_id(naggr + 1, 0);
+        for(ptrdiff_t i = 0; i < n; ++i) {
+            if (loc_owner[i] == comm.rank && loc_state[i] >= 0)
+                new_id[loc_state[i] + 1] = 1;
+        }
+
+        for(size_t i = 0; i < Sp.recv.count(); ++i) {
+            if (rem_owner[i] == comm.rank && rem_state[i] >= 0)
+                new_id[rem_state[i] + 1] = 1;
+        }
+
+        std::partial_sum(new_id.begin(), new_id.end(), new_id.begin());
+
+        if (comm.reduce(MPI_SUM, naggr - new_id.back()) > 0) {
+            naggr = new_id.back();
+
+            for(ptrdiff_t i = 0; i < n; ++i) {
+                if (loc_owner[i] == comm.rank && loc_state[i] >= 0) {
+                    loc_state[i] = new_id[loc_state[i]];
+                }
+            }
+
+            for(size_t i = 0; i < Sp.recv.nbr.size(); ++i) {
+                send_pts[i].clear();
+            }
+
+            for (auto p = Sp.remote_begin(); p!= Sp.remote_end(); ++p) {
+                ptrdiff_t c = p->first;
+
+                int d, k;
+                std::tie(d, k) = p->second;
+
+                if (rem_owner[k] == comm.rank && rem_state[k] >= 0) {
+                    send_pts[d].push_back(c);
+                    send_pts[d].push_back(new_id[rem_state[k]]);
+                }
+            }
+
+            for(size_t i = 0; i < Sp.recv.nbr.size(); ++i) {
+                int npts = send_pts[i].size();
+                MPI_Isend(&npts, 1, MPI_INT, Sp.recv.nbr[i], tag_exc_cnt, comm, &send_cnt_req[i]);
+
+                if (!npts) continue;
+                MPI_Isend(&send_pts[i][0], npts, datatype<ptrdiff_t>(), Sp.recv.nbr[i], tag_exc_pts, comm, &send_pts_req[i]);
+            }
+
+            for(size_t i = 0; i < Sp.send.nbr.size(); ++i) {
+                int npts;
+                MPI_Recv(&npts, 1, MPI_INT, Sp.send.nbr[i], tag_exc_cnt, comm, MPI_STATUS_IGNORE);
+
+                if (!npts) continue;
+                recv_pts.resize(npts);
+                MPI_Recv(&recv_pts[0], npts, datatype<ptrdiff_t>(), Sp.send.nbr[i], tag_exc_pts, comm, MPI_STATUS_IGNORE);
+
+                for(int k = 0; k < npts; k += 2) {
+                    ptrdiff_t c  = recv_pts[k] - Sp.loc_col_shift();
+                    ptrdiff_t id = recv_pts[k+1];
+
+                    loc_state[c] = id;
+                }
+            }
+        }
+
+        AMGCL_TOC("drop empty aggregates");
         AMGCL_TOC("PMIS");
 
         return naggr;

amgcl/mpi/distributed_matrix.hpp

@@ -201,6 +201,16 @@ class comm_pattern {
             return idx.at(col);
         }
 
+        std::unordered_map<ptrdiff_t, std::tuple<int, int> >::const_iterator
+        remote_begin() const {
+            return idx.cbegin();
+        }
+
+        std::unordered_map<ptrdiff_t, std::tuple<int, int> >::const_iterator
+        remote_end() const {
+            return idx.cend();
+        }
+
         size_t renumber(size_t n, ptrdiff_t *col) const {
             for(size_t i = 0; i < n; ++i)
                 col[i] = std::get<1>(idx.at(col[i]));