MueLu Version of the Day
Loading...
Searching...
No Matches
MueLu_LocalAggregationAlgorithm_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_LOCALAGGREGATIONALGORITHM_DEF_HPP
47#define MUELU_LOCALAGGREGATIONALGORITHM_DEF_HPP
48
49#include <Teuchos_Comm.hpp>
50#include <Teuchos_CommHelpers.hpp>
51
52#include <Xpetra_Vector.hpp>
53
55
56#include "MueLu_Aggregates.hpp"
57#include "MueLu_Exceptions.hpp"
58#include "MueLu_GraphBase.hpp"
59#include "MueLu_LinkedList.hpp"
60#include "MueLu_Monitor.hpp"
61#include "MueLu_Utilities.hpp"
62
63namespace MueLu {
64
65 template <class LocalOrdinal, class GlobalOrdinal, class Node>
67 : ordering_("natural"), minNodesPerAggregate_(1), maxNeighAlreadySelected_(0)
68 { }
69
70 template <class LocalOrdinal, class GlobalOrdinal, class Node>
72 Monitor m(*this, "Coarsen Uncoupled");
73
74 GetOStream(Runtime1) << "Ordering: " << ordering_ << std::endl;
75 GetOStream(Runtime1) << "Min nodes per aggregate: " << minNodesPerAggregate_ << std::endl;
76 GetOStream(Runtime1) << "Max nbrs already selected: " << maxNeighAlreadySelected_ << std::endl;
77
78 /* Create Aggregation object */
79 my_size_t nAggregates = 0;
80
81 /* ============================================================= */
82 /* aggStat indicates whether this node has been aggreated, and */
83 /* vertex2AggId stores the aggregate number where this node has */
84 /* been aggregated into. */
85 /* ============================================================= */
86
87 Teuchos::ArrayRCP<CANodeState> aggStat;
88 const my_size_t nRows = graph.GetNodeNumVertices();
89 if (nRows > 0) aggStat = Teuchos::arcp<CANodeState>(nRows);
90 for ( my_size_t i = 0; i < nRows; ++i ) aggStat[i] = CA_READY;
91
92 /* ============================================================= */
93 /* Phase 1 : */
94 /* for all nodes, form a new aggregate with its neighbors */
95 /* if the number of its neighbors having been aggregated does */
96 /* not exceed a given threshold */
97 /* (GetMaxNeighAlreadySelected() = 0 ===> Vanek's scheme) */
98 /* ============================================================= */
99
100 /* some general variable declarations */
101 Teuchos::ArrayRCP<LO> randomVector;
102 RCP<MueLu::LinkedList> nodeList; /* list storing the next node to pick as a root point for ordering_ == "graph" */
103 MueLu_SuperNode *aggHead=NULL, *aggCurrent=NULL, *supernode=NULL;
104
105
106 if ( ordering_ == "random" ) /* random ordering */
107 {
108 //TODO: could be stored in a class that respect interface of LinkedList
109
110 randomVector = Teuchos::arcp<LO>(nRows); //size_t or int ?-> to be propagated
111 for (my_size_t i = 0; i < nRows; ++i) randomVector[i] = i;
112 RandomReorder(randomVector);
113 }
114 else if ( ordering_ == "graph" ) /* graph ordering */
115 {
116 nodeList = rcp(new MueLu::LinkedList());
117 nodeList->Add(0);
118 }
119
120 /* main loop */
121 {
122 LO iNode = 0;
123 LO iNode2 = 0;
124
125 Teuchos::ArrayRCP<LO> vertex2AggId = aggregates.GetVertex2AggId()->getDataNonConst(0); // output only: contents ignored
126
127 while (iNode2 < nRows)
128 {
129
130 /*------------------------------------------------------ */
131 /* pick the next node to aggregate */
132 /*------------------------------------------------------ */
133
134 if ( ordering_ == "natural" ) iNode = iNode2++;
135 else if ( ordering_ == "random" ) iNode = randomVector[iNode2++];
136 else if ( ordering_ == "graph" )
137 {
138 if ( nodeList->IsEmpty() )
139 {
140 for ( int jNode = 0; jNode < nRows; ++jNode )
141 {
142 if ( aggStat[jNode] == CA_READY )
143 {
144 nodeList->Add(jNode); //TODO optim: not necessary to create a node. Can just set iNode value and skip the end
145 break;
146 }
147 }
148 }
149 if ( nodeList->IsEmpty() ) break; /* end of the while loop */ //TODO: coding style :(
150
151 iNode = nodeList->Pop();
152 }
153 else {
154 throw(Exceptions::RuntimeError("CoarsenUncoupled: bad aggregation ordering option"));
155 }
156
157 /*------------------------------------------------------ */
158 /* consider further only if the node is in CA_READY mode */
159 /*------------------------------------------------------ */
160
161 if ( aggStat[iNode] == CA_READY )
162 {
163 // neighOfINode is the neighbor node list of node 'iNode'.
164 Teuchos::ArrayView<const LO> neighOfINode = graph.getNeighborVertices(iNode);
165 typename Teuchos::ArrayView<const LO>::size_type length = neighOfINode.size();
166
167 supernode = new MueLu_SuperNode;
168 try {
169 supernode->list = Teuchos::arcp<int>(length+1);
170 } catch (std::bad_alloc&) {
171 TEUCHOS_TEST_FOR_EXCEPTION(true, Exceptions::RuntimeError, "MueLu::LocalAggregationAlgorithm::CoarsenUncoupled(): Error: couldn't allocate memory for supernode! length=" + Teuchos::toString(length));
172 }
173
174 supernode->maxLength = length;
175 supernode->length = 1;
176 supernode->list[0] = iNode;
177
178 int selectFlag = 1;
179 {
180 /*--------------------------------------------------- */
181 /* count the no. of neighbors having been aggregated */
182 /*--------------------------------------------------- */
183
184 int count = 0;
185 for (typename Teuchos::ArrayView<const LO>::const_iterator it = neighOfINode.begin(); it != neighOfINode.end(); ++it)
186 {
187 int index = *it;
188 if ( index < nRows )
189 {
190 if ( aggStat[index] == CA_READY ||
191 aggStat[index] == CA_NOTSEL )
192 supernode->list[supernode->length++] = index;
193 else count++;
194
195 }
196 }
197
198 /*--------------------------------------------------- */
199 /* if there are too many neighbors aggregated or the */
200 /* number of nodes in the new aggregate is too few, */
201 /* don't do this one */
202 /*--------------------------------------------------- */
203
204 if ( count > GetMaxNeighAlreadySelected() ) selectFlag = 0;
205 }
206
207 // Note: the supernode length is actually 1 more than the
208 // number of nodes in the candidate aggregate. The
209 // root is counted twice. I'm not sure if this is
210 // a bug or a feature ... so I'll leave it and change
211 // < to <= in the if just below.
212
213 if (selectFlag != 1 ||
214 supernode->length <= GetMinNodesPerAggregate())
215 {
216 aggStat[iNode] = CA_NOTSEL;
217 delete supernode;
218 if ( ordering_ == "graph" ) /* if graph ordering */
219 {
220 for (typename Teuchos::ArrayView<const LO>::const_iterator it = neighOfINode.begin(); it != neighOfINode.end(); ++it)
221 {
222 int index = *it;
223 if ( index < nRows && aggStat[index] == CA_READY )
224 {
225 nodeList->Add(index);
226 }
227 }
228 }
229 }
230 else
231 {
232 aggregates.SetIsRoot(iNode);
233 for ( int j = 0; j < supernode->length; ++j )
234 {
235 int jNode = supernode->list[j];
236 aggStat[jNode] = CA_SELECTED;
237 vertex2AggId[jNode] = nAggregates;
238 if ( ordering_ == "graph" ) /* if graph ordering */
239 {
240
241 Teuchos::ArrayView<const LO> neighOfJNode = graph.getNeighborVertices(jNode);
242
243 for (typename Teuchos::ArrayView<const LO>::const_iterator it = neighOfJNode.begin(); it != neighOfJNode.end(); ++it)
244 {
245 int index = *it;
246 if ( index < nRows && aggStat[index] == CA_READY )
247 {
248 nodeList->Add(index);
249 }
250 }
251 }
252 }
253 supernode->next = NULL;
254 supernode->index = nAggregates;
255 if ( nAggregates == 0 )
256 {
257 aggHead = supernode;
258 aggCurrent = supernode;
259 }
260 else
261 {
262 aggCurrent->next = supernode;
263 aggCurrent = supernode;
264 }
265 nAggregates++;
266 // unused aggCntArray[nAggregates] = supernode->length;
267 }
268 }
269 } // end of 'for'
270
271 // views on distributed vectors are freed here.
272
273 } // end of 'main loop'
274
275 nodeList = Teuchos::null;
276
277 /* Update aggregate object */
278 aggregates.SetNumAggregates(nAggregates);
279
280 /* Verbose */
281 {
282 const RCP<const Teuchos::Comm<int> > & comm = graph.GetComm();
283
284 if (IsPrint(Warnings0)) {
285 GO localReady=0, globalReady;
286
287 // Compute 'localReady'
288 for ( my_size_t i = 0; i < nRows; ++i )
289 if (aggStat[i] == CA_READY) localReady++;
290
291 // Compute 'globalReady'
292 MueLu_sumAll(comm, localReady, globalReady);
293
294 if(globalReady > 0)
295 GetOStream(Warnings0) << globalReady << " CA_READY nodes left" << std::endl;
296 }
297
298 if (IsPrint(Statistics1)) {
299 // Compute 'localSelected'
300 LO localSelected=0;
301 for ( my_size_t i = 0; i < nRows; ++i )
302 if ( aggStat[i] == CA_SELECTED ) localSelected++;
303
304 // Compute 'globalSelected'
305 GO globalSelected; MueLu_sumAll(comm, (GO)localSelected, globalSelected);
306
307 // Compute 'globalNRows'
308 GO globalNRows; MueLu_sumAll(comm, (GO)nRows, globalNRows);
309
310 GetOStream(Statistics1) << "Nodes aggregated = " << globalSelected << " (" << globalNRows << ")" << std::endl;
311 }
312
313 if (IsPrint(Statistics1)) {
314 GO nAggregatesGlobal; MueLu_sumAll(comm, (GO)nAggregates, nAggregatesGlobal);
315 GetOStream(Statistics1) << "Total aggregates = " << nAggregatesGlobal << std::endl;
316 }
317
318 } // verbose
319
320 /* ------------------------------------------------------------- */
321 /* clean up */
322 /* ------------------------------------------------------------- */
323
324 aggCurrent = aggHead;
325 while ( aggCurrent != NULL )
326 {
327 supernode = aggCurrent;
328 aggCurrent = aggCurrent->next;
329 delete supernode;
330 }
331
332 } // CoarsenUncoupled
333
334 template <class LocalOrdinal, class GlobalOrdinal, class Node>
336 //TODO: replace int
337 int n = list.size();
338 for(int i=0; i<n-1; i++) {
339 std::swap(list[i], list[RandomOrdinal(i,n-1)]);
340 }
341 }
342
343 template <class LocalOrdinal, class GlobalOrdinal, class Node>
345 return min + Teuchos::as<int>((max-min+1) * (static_cast<double>(std::rand()) / (RAND_MAX + 1.0)));
346 }
347
348} // namespace MueLu
349
350#endif // MUELU_LOCALAGGREGATIONALGORITHM_DEF_HPP
#define MueLu_sumAll(rcpComm, in, out)
Container class for aggregation information.
void SetIsRoot(LO i, bool value=true)
Set root node information.
const RCP< LOMultiVector > & GetVertex2AggId() const
Returns constant vector that maps local node IDs to local aggregates IDs.
void SetNumAggregates(LO nAggregates)
Set number of local aggregates on current processor.
Exception throws to report errors in the internal logical of the program.
MueLu representation of a graph.
virtual const RCP< const Teuchos::Comm< int > > GetComm() const =0
virtual Teuchos::ArrayView< const LocalOrdinal > getNeighborVertices(LocalOrdinal v) const =0
Return the list of vertices adjacent to the vertex 'v'.
virtual size_t GetNodeNumVertices() const =0
Return number of vertices owned by the calling node.
int RandomOrdinal(int min, int max) const
Generate a random number in the range [min, max].
void RandomReorder(Teuchos::ArrayRCP< LO > list) const
Utility to take a list of integers and reorder them randomly (by using a local permutation).
void CoarsenUncoupled(GraphBase const &graph, Aggregates &aggregates) const
Local aggregation.
Timer to be used in non-factories.
Namespace for MueLu classes and methods.
@ Warnings0
Important warning messages (one line)
@ Statistics1
Print more statistics.
@ Runtime1
Description of what is happening (more verbose)
struct MueLu::MueLu_SuperNode_Struct MueLu_SuperNode