MueLu  Version of the Day
MueLu_AggregationPhase1Algorithm_kokkos_def.hpp
Go to the documentation of this file.
1 // @HEADER
2 //
3 // ***********************************************************************
4 //
5 // MueLu: A package for multigrid based preconditioning
6 // Copyright 2012 Sandia Corporation
7 //
8 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
9 // the U.S. Government retains certain rights in this software.
10 //
11 // Redistribution and use in source and binary forms, with or without
12 // modification, are permitted provided that the following conditions are
13 // met:
14 //
15 // 1. Redistributions of source code must retain the above copyright
16 // notice, this list of conditions and the following disclaimer.
17 //
18 // 2. Redistributions in binary form must reproduce the above copyright
19 // notice, this list of conditions and the following disclaimer in the
20 // documentation and/or other materials provided with the distribution.
21 //
22 // 3. Neither the name of the Corporation nor the names of the
23 // contributors may be used to endorse or promote products derived from
24 // this software without specific prior written permission.
25 //
26 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
27 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
28 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
29 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
30 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
31 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
32 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
33 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
34 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
35 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
36 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
37 //
38 // Questions? Contact
39 // Jonathan Hu (jhu@sandia.gov)
40 // Andrey Prokopenko (aprokop@sandia.gov)
41 // Ray Tuminaro (rstumin@sandia.gov)
42 //
43 // ***********************************************************************
44 //
45 // @HEADER
46 #ifndef MUELU_AGGREGATIONPHASE1ALGORITHM_KOKKOS_DEF_HPP
47 #define MUELU_AGGREGATIONPHASE1ALGORITHM_KOKKOS_DEF_HPP
48 
49 #ifdef HAVE_MUELU_KOKKOS_REFACTOR
50 
51 #include <queue>
52 
53 #include <Teuchos_Comm.hpp>
54 #include <Teuchos_CommHelpers.hpp>
55 
56 #include <Xpetra_Vector.hpp>
57 
59 
60 #include "MueLu_Aggregates_kokkos.hpp"
61 #include "MueLu_Exceptions.hpp"
62 #include "MueLu_LWGraph_kokkos.hpp"
63 #include "MueLu_Monitor.hpp"
64 
65 namespace MueLu {
66 
67  template <class LocalOrdinal, class GlobalOrdinal, class Node>
68  void AggregationPhase1Algorithm_kokkos<LocalOrdinal, GlobalOrdinal, Node>::
69  BuildAggregates(const ParameterList& params, const LWGraph_kokkos& graph, Aggregates_kokkos& aggregates, std::vector<unsigned>& aggStat,
70  LO& numNonAggregatedNodes) const {
71  Monitor m(*this, "BuildAggregates");
72 
73  std::string ordering = params.get<std::string>("aggregation: ordering");
74  int maxNeighAlreadySelected = params.get<int> ("aggregation: max selected neighbors");
75  int minNodesPerAggregate = params.get<int> ("aggregation: min agg size");
76  int maxNodesPerAggregate = params.get<int> ("aggregation: max agg size");
77 
78  TEUCHOS_TEST_FOR_EXCEPTION(maxNodesPerAggregate < minNodesPerAggregate, Exceptions::RuntimeError,
79  "MueLu::UncoupledAggregationAlgorithm::BuildAggregates: minNodesPerAggregate must be smaller or equal to MaxNodePerAggregate!");
80 
81  const LO numRows = graph.GetNodeNumVertices();
82  const int myRank = graph.GetComm()->getRank();
83 
84  ArrayRCP<LO> vertex2AggId = aggregates.GetVertex2AggId()->getDataNonConst(0);
85  ArrayRCP<LO> procWinner = aggregates.GetProcWinner() ->getDataNonConst(0);
86 
87  LO numLocalAggregates = aggregates.GetNumAggregates();
88 
89  ArrayRCP<LO> randomVector;
90  if (ordering == "random") {
91  randomVector = arcp<LO>(numRows);
92  for (LO i = 0; i < numRows; i++)
93  randomVector[i] = i;
94  RandomReorder(randomVector);
95  }
96 
97  int aggIndex = -1;
98  size_t aggSize = 0;
99  std::vector<int> aggList(graph.getNodeMaxNumRowEntries());
100 
101  std::queue<LO> graphOrderQueue;
102 
103  // Main loop over all local rows of graph(A)
104  for (LO i = 0; i < numRows; i++) {
105  // Step 1: pick the next node to aggregate
106  LO rootCandidate = 0;
107  if (ordering == "natural") rootCandidate = i;
108  else if (ordering == "random") rootCandidate = randomVector[i];
109  else if (ordering == "graph") {
110 
111  if (graphOrderQueue.size() == 0) {
112  // Current queue is empty for "graph" ordering, populate with one READY node
113  for (LO jnode = 0; jnode < numRows; jnode++)
114  if (aggStat[jnode] == READY) {
115  graphOrderQueue.push(jnode);
116  break;
117  }
118  }
119  if (graphOrderQueue.size() == 0) {
120  // There are no more ready nodes, end the phase
121  break;
122  }
123  rootCandidate = graphOrderQueue.front(); // take next node from graph ordering queue
124  graphOrderQueue.pop(); // delete this node in list
125  }
126 
127  if (aggStat[rootCandidate] != READY)
128  continue;
129 
130  // Step 2: build tentative aggregate
131  aggSize = 0;
132  aggList[aggSize++] = rootCandidate;
133 
134  typename LWGraph_kokkos::row_type neighOfINode = graph.getNeighborVertices(rootCandidate);
135 
136  // If the number of neighbors is less than the minimum number of nodes
137  // per aggregate, we know this is not going to be a valid root, and we
138  // may skip it, but only for "natural" and "random" (for "graph" we still
139  // need to fetch the list of local neighbors to continue)
140  if ((ordering == "natural" || ordering == "random") &&
141  neighOfINode.size() < minNodesPerAggregate) {
142  continue;
143  }
144 
145  LO numAggregatedNeighbours = 0;
146 
147  for (int j = 0; j < neighOfINode.size(); j++) {
148  LO neigh = neighOfINode(j);
149 
150  if (neigh != rootCandidate && graph.isLocalNeighborVertex(neigh)) {
151 
152  if (aggStat[neigh] == READY || aggStat[neigh] == NOTSEL) {
153  // If aggregate size does not exceed max size, add node to the
154  // tentative aggregate
155  // NOTE: We do not exit the loop over all neighbours since we have
156  // still to count all aggregated neighbour nodes for the
157  // aggregation criteria
158  // NOTE: We check here for the maximum aggregation size. If we
159  // would do it below with all the other check too big aggregates
160  // would not be accepted at all.
161  if (aggSize < as<size_t>(maxNodesPerAggregate))
162  aggList[aggSize++] = neigh;
163 
164  } else {
165  numAggregatedNeighbours++;
166  }
167  }
168  }
169 
170  // Step 3: check if tentative aggregate is acceptable
171  if ((numAggregatedNeighbours <= maxNeighAlreadySelected) && // too many connections to other aggregates
172  (aggSize >= as<size_t>(minNodesPerAggregate))) { // too few nodes in the tentative aggregate
173  // Accept new aggregate
174  // rootCandidate becomes the root of the newly formed aggregate
175  aggregates.SetIsRoot(rootCandidate);
176  aggIndex = numLocalAggregates++;
177 
178  for (size_t k = 0; k < aggSize; k++) {
179  aggStat [aggList[k]] = AGGREGATED;
180  vertex2AggId[aggList[k]] = aggIndex;
181  procWinner [aggList[k]] = myRank;
182  }
183 
184  numNonAggregatedNodes -= aggSize;
185 
186  } else {
187  // Aggregate is not accepted
188  aggStat[rootCandidate] = NOTSEL;
189 
190  // Need this for the "graph" ordering below
191  // The original candidate is always aggList[0]
192  aggSize = 1;
193  }
194 
195  if (ordering == "graph") {
196  // Add candidates to the list of nodes
197  // NOTE: the code have slightly different meanings depending on context:
198  // - if aggregate was accepted, we add neighbors of neighbors of the original candidate
199  // - if aggregate was not accepted, we add neighbors of the original candidate
200  for (size_t k = 0; k < aggSize; k++) {
201  typename LWGraph_kokkos::row_type neighOfJNode = graph.getNeighborVertices(aggList[k]);
202 
203  for (int j = 0; j < neighOfJNode.size(); j++) {
204  LO neigh = neighOfJNode(j);
205 
206  if (graph.isLocalNeighborVertex(neigh) && aggStat[neigh] == READY)
207  graphOrderQueue.push(neigh);
208  }
209  }
210  }
211  }
212 
213  // Reset all NOTSEL vertices to READY
214  // This simplifies other algorithms
215  for (LO i = 0; i < numRows; i++)
216  if (aggStat[i] == NOTSEL)
217  aggStat[i] = READY;
218 
219  // update aggregate object
220  aggregates.SetNumAggregates(numLocalAggregates);
221  }
222 
223  template <class LocalOrdinal, class GlobalOrdinal, class Node>
224  void AggregationPhase1Algorithm_kokkos<LocalOrdinal, GlobalOrdinal, Node>::RandomReorder(ArrayRCP<LO> list) const {
225  //TODO: replace int
226  int n = list.size();
227  for(int i = 0; i < n-1; i++)
228  std::swap(list[i], list[RandomOrdinal(i,n-1)]);
229  }
230 
231  template <class LocalOrdinal, class GlobalOrdinal, class Node>
232  int AggregationPhase1Algorithm_kokkos<LocalOrdinal, GlobalOrdinal, Node>::RandomOrdinal(int min, int max) const {
233  return min + as<int>((max-min+1) * (static_cast<double>(std::rand()) / (RAND_MAX + 1.0)));
234  }
235 
236 } // end namespace
237 
238 #endif // HAVE_MUELU_KOKKOS_REFACTOR
239 #endif // MUELU_AGGREGATIONPHASE1ALGORITHM_KOKKOS_DEF_HPP
LocalOrdinal LO
Namespace for MueLu classes and methods.