MueLu Version of the Day
Loading...
Searching...
No Matches
MueLu_AggregationPhase3Algorithm_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_AGGREGATIONPHASE3ALGORITHM_KOKKOS_DEF_HPP
47#define MUELU_AGGREGATIONPHASE3ALGORITHM_KOKKOS_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_kokkos.hpp"
57#include "MueLu_Exceptions.hpp"
58#include "MueLu_LWGraph_kokkos.hpp"
59#include "MueLu_Monitor.hpp"
60
61// #include "Kokkos_Core.hpp"
62
63namespace MueLu {
64
65 // Try to stick unaggregated nodes into a neighboring aggregate if they are
66 // not already too big. Otherwise, make a new aggregate
67 template <class LocalOrdinal, class GlobalOrdinal, class Node>
69 BuildAggregates(const ParameterList& params,
70 const LWGraph_kokkos& graph,
71 Aggregates_kokkos& aggregates,
72 Kokkos::View<unsigned*, typename LWGraph_kokkos::device_type>& aggStat,
73 LO& numNonAggregatedNodes) const {
74
75 // So far we only have the non-deterministic version of the algorithm...
76 if(params.get<bool>("aggregation: deterministic")) {
77 Monitor m(*this, "BuildAggregatesDeterministic");
78 BuildAggregatesRandom(params, graph, aggregates, aggStat, numNonAggregatedNodes);
79 } else {
80 Monitor m(*this, "BuildAggregatesRandom");
81 BuildAggregatesRandom(params, graph, aggregates, aggStat, numNonAggregatedNodes);
82 }
83
84 }
85
86 // Try to stick unaggregated nodes into a neighboring aggregate if they are
87 // not already too big. Otherwise, make a new aggregate
88 template <class LocalOrdinal, class GlobalOrdinal, class Node>
90 BuildAggregatesRandom(const ParameterList& params,
91 const LWGraph_kokkos& graph,
92 Aggregates_kokkos& aggregates,
93 Kokkos::View<unsigned*, typename LWGraph_kokkos::device_type>& aggStat,
94 LO& numNonAggregatedNodes) const {
95
96 bool error_on_isolated = params.get<bool>("aggregation: error on nodes with no on-rank neighbors");
97 bool makeNonAdjAggs = params.get<bool>("aggregation: phase3 avoid singletons");
98
99 const LO numRows = graph.GetNodeNumVertices();
100 const int myRank = graph.GetComm()->getRank();
101
102 auto vertex2AggId = aggregates.GetVertex2AggId()->getDeviceLocalView(Xpetra::Access::ReadWrite);
103 auto procWinner = aggregates.GetProcWinner() ->getDeviceLocalView(Xpetra::Access::ReadWrite);
104 auto colors = aggregates.GetGraphColors();
105 const LO numColors = aggregates.GetGraphNumColors();
106
107 auto lclLWGraph = graph.getLocalLWGraph();
108
109 Kokkos::View<LO, device_type> numAggregates("numAggregates");
110 Kokkos::deep_copy(numAggregates, aggregates.GetNumAggregates());
111
112 Kokkos::View<unsigned*, device_type> aggStatOld("Initial aggregation status", aggStat.extent(0));
113 Kokkos::deep_copy(aggStatOld, aggStat);
114 Kokkos::View<LO, device_type> numNonAggregated("numNonAggregated");
115 Kokkos::deep_copy(numNonAggregated, numNonAggregatedNodes);
116 for(int color = 1; color < numColors + 1; ++color) {
117 Kokkos::parallel_for("Aggregation Phase 3: aggregates clean-up",
118 Kokkos::RangePolicy<execution_space>(0, numRows),
119 KOKKOS_LAMBDA(const LO nodeIdx) {
120 // Check if node has already been treated?
121 if( (colors(nodeIdx) != color) ||
122 (aggStatOld(nodeIdx) == AGGREGATED) ||
123 (aggStatOld(nodeIdx) == IGNORED) ){ return; }
124
125 // Grab node neighbors
126 auto neighbors = lclLWGraph.getNeighborVertices(nodeIdx);
127 LO neighIdx;
128
129 // We don't want a singleton.
130 // So lets see if any neighbors can be used to form a new aggregate?
131 bool isNewAggregate = false;
132 for(int neigh = 0; neigh < neighbors.length; ++neigh) {
133 neighIdx = neighbors(neigh);
134
135 if((neighIdx != nodeIdx) &&
136 lclLWGraph.isLocalNeighborVertex(neighIdx) &&
137 (aggStatOld(neighIdx) == READY)) {
138 isNewAggregate = true;
139 break;
140 }
141 }
142
143 // We can form a new non singleton aggregate!
144 if(isNewAggregate) {
145 // If this is the aggregate root
146 // we need to process the nodes in the aggregate
147 const LO aggId = Kokkos::atomic_fetch_add(&numAggregates(), 1);
148 aggStat(nodeIdx) = AGGREGATED;
149 procWinner(nodeIdx, 0) = myRank;
150 vertex2AggId(nodeIdx, 0) = aggId;
151 // aggregates.SetIsRoot(nodeIdx);
152 Kokkos::atomic_decrement(&numNonAggregated());
153 for(int neigh = 0; neigh < neighbors.length; ++neigh) {
154 neighIdx = neighbors(neigh);
155 if((neighIdx != nodeIdx) &&
156 lclLWGraph.isLocalNeighborVertex(neighIdx) &&
157 (aggStatOld(neighIdx) == READY)) {
158 aggStat(neighIdx) = AGGREGATED;
159 procWinner(neighIdx, 0) = myRank;
160 vertex2AggId(neighIdx, 0) = aggId;
161 Kokkos::atomic_decrement(&numNonAggregated());
162 }
163 }
164 return;
165 }
166
167 // Getting a little desperate!
168 // Let us try to aggregate into a neighboring aggregate
169 for(int neigh = 0; neigh < neighbors.length; ++neigh) {
170 neighIdx = neighbors(neigh);
171 if (lclLWGraph.isLocalNeighborVertex(neighIdx) &&
172 (aggStatOld(neighIdx) == AGGREGATED)) {
173 aggStat(nodeIdx) = AGGREGATED;
174 procWinner(nodeIdx, 0) = myRank;
175 vertex2AggId(nodeIdx, 0) = vertex2AggId(neighIdx, 0);
176 Kokkos::atomic_decrement(&numNonAggregated());
177 return;
178 }
179 }
180
181 // Getting quite desperate!
182 // Let us try to make a non contiguous aggregate
183 if(makeNonAdjAggs) {
184 for(LO otherNodeIdx = 0; otherNodeIdx < numRows; ++otherNodeIdx) {
185 if((otherNodeIdx != nodeIdx) &&
186 (aggStatOld(otherNodeIdx) == AGGREGATED)) {
187 aggStat(nodeIdx) = AGGREGATED;
188 procWinner(nodeIdx, 0) = myRank;
189 vertex2AggId(nodeIdx, 0) = vertex2AggId(otherNodeIdx, 0);
190 Kokkos::atomic_decrement(&numNonAggregated());
191 return;
192 }
193 }
194 }
195
196 // Total deperation!
197 // Let us make a singleton
198 if(!error_on_isolated) {
199 const LO aggId = Kokkos::atomic_fetch_add(&numAggregates(), 1);
200 aggStat(nodeIdx) = AGGREGATED;
201 procWinner(nodeIdx, 0) = myRank;
202 vertex2AggId(nodeIdx, 0) = aggId;
203 Kokkos::atomic_decrement(&numNonAggregated());
204 }
205 });
206 // LBV on 09/27/19: here we could copy numNonAggregated to host
207 // and check for it to be equal to 0 in which case we can stop
208 // looping over the different colors...
209 Kokkos::deep_copy(aggStatOld, aggStat);
210 } // loop over colors
211
212 auto numNonAggregated_h = Kokkos::create_mirror_view(numNonAggregated);
213 Kokkos::deep_copy(numNonAggregated_h, numNonAggregated);
214 numNonAggregatedNodes = numNonAggregated_h();
215 if( (error_on_isolated) && (numNonAggregatedNodes > 0) ) {
216 // Error on this isolated node, as the user has requested
217 std::ostringstream oss;
218 oss<<"MueLu::AggregationPhase3Algorithm::BuildAggregates: MueLu has detected a non-Dirichlet node that has no on-rank neighbors and is terminating (by user request). "<<std::endl;
219 oss<<"If this error is being generated at level 0, this is due to an initial partitioning problem in your matrix."<<std::endl;
220 oss<<"If this error is being generated at any other level, try turning on repartitioning, which may fix this problem."<<std::endl;
221 throw Exceptions::RuntimeError(oss.str());
222 }
223
224 // update aggregate object
225 auto numAggregates_h = Kokkos::create_mirror_view(numAggregates);
226 Kokkos::deep_copy(numAggregates_h, numAggregates);
227 aggregates.SetNumAggregates(numAggregates_h());
228 }
229
230} // end namespace
231
232#endif // MUELU_AGGREGATIONPHASE3ALGORITHM_KOKKOS_DEF_HPP
void BuildAggregatesRandom(const ParameterList &params, const LWGraph_kokkos &graph, Aggregates_kokkos &aggregates, Kokkos::View< unsigned *, device_type > &aggStat, LO &numNonAggregatedNodes) const
void BuildAggregates(const ParameterList &params, const LWGraph_kokkos &graph, Aggregates_kokkos &aggregates, Kokkos::View< unsigned *, device_type > &aggStat, LO &numNonAggregatedNodes) const
Local aggregation.
Exception throws to report errors in the internal logical of the program.
Lightweight MueLu representation of a compressed row storage graph.
Timer to be used in non-factories.
Namespace for MueLu classes and methods.