MueLu  Version of the Day
MueLu_UncoupledAggregationFactory_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_UNCOUPLEDAGGREGATIONFACTORY_KOKKOS_DEF_HPP_
47 #define MUELU_UNCOUPLEDAGGREGATIONFACTORY_KOKKOS_DEF_HPP_
48 
49 #ifdef HAVE_MUELU_KOKKOS_REFACTOR
50 
51 #include <climits>
52 
53 #include <Xpetra_Map.hpp>
54 #include <Xpetra_Vector.hpp>
55 #include <Xpetra_VectorFactory.hpp>
56 
58 
59 #include "MueLu_OnePtAggregationAlgorithm_kokkos.hpp"
60 #include "MueLu_PreserveDirichletAggregationAlgorithm_kokkos.hpp"
61 #include "MueLu_IsolatedNodeAggregationAlgorithm_kokkos.hpp"
62 
63 #include "MueLu_AggregationPhase1Algorithm_kokkos.hpp"
64 #include "MueLu_AggregationPhase2aAlgorithm_kokkos.hpp"
65 #include "MueLu_AggregationPhase2bAlgorithm_kokkos.hpp"
66 #include "MueLu_AggregationPhase3Algorithm_kokkos.hpp"
67 
68 #include "MueLu_Level.hpp"
69 #include "MueLu_LWGraph_kokkos.hpp"
70 #include "MueLu_Aggregates_kokkos.hpp"
71 #include "MueLu_MasterList.hpp"
72 #include "MueLu_Monitor.hpp"
73 #include "MueLu_AmalgamationInfo.hpp"
74 #include "MueLu_Utilities.hpp" // for sum_all and similar stuff...
75 
76 namespace MueLu {
77 
78  template <class LocalOrdinal, class GlobalOrdinal, class Node>
79  UncoupledAggregationFactory_kokkos<LocalOrdinal, GlobalOrdinal, Node>::UncoupledAggregationFactory_kokkos()
80  : bDefinitionPhase_(true)
81  { }
82 
83  template <class LocalOrdinal, class GlobalOrdinal, class Node>
84  RCP<const ParameterList> UncoupledAggregationFactory_kokkos<LocalOrdinal, GlobalOrdinal, Node>::GetValidParameterList() const {
85  RCP<ParameterList> validParamList = rcp(new ParameterList());
86 
87  // Aggregation parameters (used in aggregation algorithms)
88  // TODO introduce local member function for each aggregation algorithm such that each aggregation algorithm can define its own parameters
89 
90  typedef Teuchos::StringToIntegralParameterEntryValidator<int> validatorType;
91 #define SET_VALID_ENTRY(name) validParamList->setEntry(name, MasterList::getEntry(name))
92  SET_VALID_ENTRY("aggregation: max agg size");
93  SET_VALID_ENTRY("aggregation: min agg size");
94  SET_VALID_ENTRY("aggregation: max selected neighbors");
95  SET_VALID_ENTRY("aggregation: ordering");
96  validParamList->getEntry("aggregation: ordering").setValidator(
97  rcp(new validatorType(Teuchos::tuple<std::string>("natural", "graph", "random"), "aggregation: ordering")));
98  SET_VALID_ENTRY("aggregation: enable phase 1");
99  SET_VALID_ENTRY("aggregation: enable phase 2a");
100  SET_VALID_ENTRY("aggregation: enable phase 2b");
101  SET_VALID_ENTRY("aggregation: enable phase 3");
102  SET_VALID_ENTRY("aggregation: preserve Dirichlet points");
103  SET_VALID_ENTRY("aggregation: allow user-specified singletons");
104 #undef SET_VALID_ENTRY
105 
106  // general variables needed in AggregationFactory
107  validParamList->set< RCP<const FactoryBase> >("Graph", null, "Generating factory of the graph");
108  validParamList->set< RCP<const FactoryBase> >("DofsPerNode", null, "Generating factory for variable \'DofsPerNode\', usually the same as for \'Graph\'");
109 
110  // special variables necessary for OnePtAggregationAlgorithm
111  validParamList->set< std::string > ("OnePt aggregate map name", "", "Name of input map for single node aggregates. (default='')");
112  validParamList->set< std::string > ("OnePt aggregate map factory", "", "Generating factory of (DOF) map for single node aggregates.");
113  //validParamList->set< RCP<const FactoryBase> >("OnePt aggregate map factory", NoFactory::getRCP(), "Generating factory of (DOF) map for single node aggregates.");
114 
115  return validParamList;
116  }
117 
118  template <class LocalOrdinal, class GlobalOrdinal, class Node>
119  void UncoupledAggregationFactory_kokkos<LocalOrdinal, GlobalOrdinal, Node>::DeclareInput(Level& currentLevel) const {
120  Input(currentLevel, "Graph");
121  Input(currentLevel, "DofsPerNode");
122 
123  const ParameterList& pL = GetParameterList();
124 
125  // request special data necessary for OnePtAggregationAlgorithm
126  std::string mapOnePtName = pL.get<std::string>("OnePt aggregate map name");
127  if (mapOnePtName.length() > 0) {
128  std::string mapOnePtFactName = pL.get<std::string>("OnePt aggregate map factory");
129  if (mapOnePtFactName == "" || mapOnePtFactName == "NoFactory") {
130  currentLevel.DeclareInput(mapOnePtName, NoFactory::get());
131  } else {
132  RCP<const FactoryBase> mapOnePtFact = GetFactory(mapOnePtFactName);
133  currentLevel.DeclareInput(mapOnePtName, mapOnePtFact.get());
134  }
135  }
136  }
137 
138  template <class LocalOrdinal, class GlobalOrdinal, class Node>
139  void UncoupledAggregationFactory_kokkos<LocalOrdinal, GlobalOrdinal, Node>::Build(Level &currentLevel) const {
140  FactoryMonitor m(*this, "Build", currentLevel);
141 
142  ParameterList pL = GetParameterList();
143  bDefinitionPhase_ = false; // definition phase is finished, now all aggregation algorithm information is fixed
144 
145  if (pL.get<int>("aggregation: max agg size") == -1)
146  pL.set("aggregation: max agg size", INT_MAX);
147 
148  // define aggregation algorithms
149  RCP<const FactoryBase> graphFact = GetFactory("Graph");
150 
151  // TODO Can we keep different aggregation algorithms over more Build calls?
152  algos_.clear();
153  algos_.push_back(rcp(new PreserveDirichletAggregationAlgorithm_kokkos(graphFact)));
154  if (pL.get<bool>("aggregation: allow user-specified singletons") == true) algos_.push_back(rcp(new OnePtAggregationAlgorithm_kokkos (graphFact)));
155  if (pL.get<bool>("aggregation: enable phase 1" ) == true) algos_.push_back(rcp(new AggregationPhase1Algorithm_kokkos (graphFact)));
156  if (pL.get<bool>("aggregation: enable phase 2a") == true) algos_.push_back(rcp(new AggregationPhase2aAlgorithm_kokkos (graphFact)));
157  if (pL.get<bool>("aggregation: enable phase 2b") == true) algos_.push_back(rcp(new AggregationPhase2bAlgorithm_kokkos (graphFact)));
158  if (pL.get<bool>("aggregation: enable phase 3" ) == true) algos_.push_back(rcp(new AggregationPhase3Algorithm_kokkos (graphFact)));
159 
160  std::string mapOnePtName = pL.get<std::string>("OnePt aggregate map name");
161  RCP<Map> OnePtMap = Teuchos::null;
162  if (mapOnePtName.length()) {
163  std::string mapOnePtFactName = pL.get<std::string>("OnePt aggregate map factory");
164  if (mapOnePtFactName == "" || mapOnePtFactName == "NoFactory") {
165  OnePtMap = currentLevel.Get<RCP<Map> >(mapOnePtName, NoFactory::get());
166  } else {
167  RCP<const FactoryBase> mapOnePtFact = GetFactory(mapOnePtFactName);
168  OnePtMap = currentLevel.Get<RCP<Map> >(mapOnePtName, mapOnePtFact.get());
169  }
170  }
171 
172  RCP<const LWGraph_kokkos> graph = Get< RCP<LWGraph_kokkos> >(currentLevel, "Graph");
173 
174  // Build
175  RCP<Aggregates_kokkos> aggregates = rcp(new Aggregates_kokkos(*graph));
176  aggregates->setObjectLabel("UC");
177 
178  const LO numRows = graph->GetNodeNumVertices();
179 
180  // construct aggStat information
181  std::vector<unsigned> aggStat(numRows, READY);
182 
183  // TODO
184  //ArrayRCP<const bool> dirichletBoundaryMap = graph->GetBoundaryNodeMap();
185  ArrayRCP<const bool> dirichletBoundaryMap;
186 
187  if (dirichletBoundaryMap != Teuchos::null)
188  for (LO i = 0; i < numRows; i++)
189  if (dirichletBoundaryMap[i] == true)
190  aggStat[i] = BOUNDARY;
191 
192  LO nDofsPerNode = Get<LO>(currentLevel, "DofsPerNode");
193  GO indexBase = graph->GetDomainMap()->getIndexBase();
194  if (OnePtMap != Teuchos::null) {
195  for (LO i = 0; i < numRows; i++) {
196  // reconstruct global row id (FIXME only works for contiguous maps)
197  GO grid = (graph->GetDomainMap()->getGlobalElement(i)-indexBase) * nDofsPerNode + indexBase;
198 
199  for (LO kr = 0; kr < nDofsPerNode; kr++)
200  if (OnePtMap->isNodeGlobalElement(grid + kr))
201  aggStat[i] = ONEPT;
202  }
203  }
204 
205 
206  const RCP<const Teuchos::Comm<int> > comm = graph->GetComm();
207  GO numGlobalRows = 0;
208  if (IsPrint(Statistics1))
209  MueLu_sumAll(comm, as<GO>(numRows), numGlobalRows);
210 
211  LO numNonAggregatedNodes = numRows;
212  GO numGlobalAggregatedPrev = 0, numGlobalAggsPrev = 0;
213  for (size_t a = 0; a < algos_.size(); a++) {
214  std::string phase = algos_[a]->description();
215  SubFactoryMonitor sfm(*this, "Algo \"" + phase + "\"", currentLevel);
216 
217  int oldRank = algos_[a]->SetProcRankVerbose(this->GetProcRankVerbose());
218  algos_[a]->BuildAggregates(pL, *graph, *aggregates, aggStat, numNonAggregatedNodes);
219  algos_[a]->SetProcRankVerbose(oldRank);
220 
221  if (IsPrint(Statistics1)) {
222  GO numLocalAggregated = numRows - numNonAggregatedNodes, numGlobalAggregated = 0;
223  GO numLocalAggs = aggregates->GetNumAggregates(), numGlobalAggs = 0;
224  MueLu_sumAll(comm, numLocalAggregated, numGlobalAggregated);
225  MueLu_sumAll(comm, numLocalAggs, numGlobalAggs);
226 
227  double aggPercent = 100*as<double>(numGlobalAggregated)/as<double>(numGlobalRows);
228  if (aggPercent > 99.99 && aggPercent < 100.00) {
229  // Due to round off (for instance, for 140465733/140466897), we could
230  // get 100.00% display even if there are some remaining nodes. This
231  // is bad from the users point of view. It is much better to change
232  // it to display 99.99%.
233  aggPercent = 99.99;
234  }
235  GetOStream(Statistics1) << " aggregated : " << (numGlobalAggregated - numGlobalAggregatedPrev) << " (phase), " << std::fixed
236  << std::setprecision(2) << numGlobalAggregated << "/" << numGlobalRows << " [" << aggPercent << "%] (total)\n"
237  << " remaining : " << numGlobalRows - numGlobalAggregated << "\n"
238  << " aggregates : " << numGlobalAggs-numGlobalAggsPrev << " (phase), " << numGlobalAggs << " (total)" << std::endl;
239  numGlobalAggregatedPrev = numGlobalAggregated;
240  numGlobalAggsPrev = numGlobalAggs;
241  }
242  }
243 
244  TEUCHOS_TEST_FOR_EXCEPTION(numNonAggregatedNodes, Exceptions::RuntimeError, "MueLu::UncoupledAggregationFactory::Build: Leftover nodes found! Error!");
245 
246  aggregates->AggregatesCrossProcessors(false);
247 
248  Set(currentLevel, "Aggregates", aggregates);
249 
250  GetOStream(Statistics0) << aggregates->description() << std::endl;
251  }
252 
253 } //namespace MueLu
254 
255 #endif // HAVE_MUELU_KOKKOS_REFACTOR
256 #endif /* MUELU_UNCOUPLEDAGGREGATIONFACTORY_DEF_HPP_ */
#define MueLu_sumAll(rcpComm, in, out)
Print more statistics.
Namespace for MueLu classes and methods.
Print statistics that do not involve significant additional computation.
#define SET_VALID_ENTRY(name)