MueLu  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
MueLu_StructuredAggregationFactory_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_STRUCTUREDAGGREGATIONFACTORY_KOKKOS_DEF_HPP
47 #define MUELU_STRUCTUREDAGGREGATIONFACTORY_KOKKOS_DEF_HPP
48 
49 #ifdef HAVE_MUELU_KOKKOS_REFACTOR
50 
51 // Xpetra includes
52 #include <Xpetra_Map.hpp>
53 #include <Xpetra_CrsGraph.hpp>
54 
55 // MueLu generic includes
56 #include "MueLu_Level.hpp"
57 #include "MueLu_MasterList.hpp"
58 #include "MueLu_Monitor.hpp"
59 #include "MueLu_Utilities.hpp"
60 
61 // MueLu specific includes (kokkos version)
62 #include "MueLu_LWGraph_kokkos.hpp"
63 #include "MueLu_Aggregates_kokkos.hpp"
64 #include "MueLu_IndexManager_kokkos.hpp"
65 #include "MueLu_AggregationStructuredAlgorithm_kokkos.hpp"
66 
68 
69 namespace MueLu {
70 
71  template <class LocalOrdinal, class GlobalOrdinal, class Node>
72  StructuredAggregationFactory_kokkos<LocalOrdinal, GlobalOrdinal, Node>::
73  StructuredAggregationFactory_kokkos() : bDefinitionPhase_(true) { }
74 
75  template <class LocalOrdinal, class GlobalOrdinal, class Node>
76  RCP<const ParameterList> StructuredAggregationFactory_kokkos<LocalOrdinal, GlobalOrdinal, Node>::
77  GetValidParameterList() const {
78  RCP<ParameterList> validParamList = rcp(new ParameterList());
79 
80 #define SET_VALID_ENTRY(name) validParamList->setEntry(name, MasterList::getEntry(name))
81  SET_VALID_ENTRY("aggregation: preserve Dirichlet points");
82  SET_VALID_ENTRY("aggregation: allow user-specified singletons");
83  SET_VALID_ENTRY("aggregation: error on nodes with no on-rank neighbors");
84  SET_VALID_ENTRY("aggregation: phase3 avoid singletons");
85 #undef SET_VALID_ENTRY
86 
87  // general variables needed in StructuredAggregationFactory
88  validParamList->set<std::string> ("aggregation: output type", "Aggregates",
89  "Type of object holding the aggregation data: Aggregtes or CrsGraph");
90  validParamList->set<std::string> ("aggregation: coarsening rate", "{3}",
91  "Coarsening rate per spatial dimensions");
92  validParamList->set<int> ("aggregation: coarsening order", 0,
93  "The interpolation order used to construct grid transfer operators based off these aggregates.");
94 
95  validParamList->set<RCP<const FactoryBase> >("Graph", Teuchos::null,
96  "Graph of the matrix after amalgamation but without dropping.");
97  validParamList->set<RCP<const FactoryBase> >("DofsPerNode", Teuchos::null,
98  "Number of degrees of freedom per mesh node, provided by the coalsce drop factory.");
99  validParamList->set<RCP<const FactoryBase> >("numDimensions", Teuchos::null,
100  "Number of spatial dimension provided by CoordinatesTransferFactory.");
101  validParamList->set<RCP<const FactoryBase> >("lNodesPerDim", Teuchos::null,
102  "Number of nodes per spatial dimmension provided by CoordinatesTransferFactory.");
103 
104  return validParamList;
105  } // GetValidParameterList()
106 
107  template <class LocalOrdinal, class GlobalOrdinal, class Node>
108  void StructuredAggregationFactory_kokkos<LocalOrdinal, GlobalOrdinal, Node>::
109  DeclareInput(Level& currentLevel) const {
110  Input(currentLevel, "Graph");
111  Input(currentLevel, "DofsPerNode");
112 
113  // Request the local number of nodes per dimensions
114  if(currentLevel.GetLevelID() == 0) {
115  if(currentLevel.IsAvailable("numDimensions", NoFactory::get())) {
116  currentLevel.DeclareInput("numDimensions", NoFactory::get(), this);
117  } else {
118  TEUCHOS_TEST_FOR_EXCEPTION(currentLevel.IsAvailable("numDimensions", NoFactory::get()),
119  Exceptions::RuntimeError,
120  "numDimensions was not provided by the user on level0!");
121  }
122  if(currentLevel.IsAvailable("lNodesPerDim", NoFactory::get())) {
123  currentLevel.DeclareInput("lNodesPerDim", NoFactory::get(), this);
124  } else {
125  TEUCHOS_TEST_FOR_EXCEPTION(currentLevel.IsAvailable("lNodesPerDim", NoFactory::get()),
126  Exceptions::RuntimeError,
127  "lNodesPerDim was not provided by the user on level0!");
128  }
129  } else {
130  Input(currentLevel, "lNodesPerDim");
131  Input(currentLevel, "numDimensions");
132  }
133  } // DeclareInput()
134 
135  template <class LocalOrdinal, class GlobalOrdinal, class Node>
136  void StructuredAggregationFactory_kokkos<LocalOrdinal, GlobalOrdinal, Node>::
137  Build(Level &currentLevel) const {
138  FactoryMonitor m(*this, "Build", currentLevel);
139 
140  RCP<Teuchos::FancyOStream> out;
141  if(const char* dbg = std::getenv("MUELU_STRUCTUREDAGGREGATION_DEBUG")) {
142  out = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
143  out->setShowAllFrontMatter(false).setShowProcRank(true);
144  } else {
145  out = Teuchos::getFancyOStream(rcp(new Teuchos::oblackholestream()));
146  }
147 
148  typedef typename LWGraph_kokkos::local_graph_type::device_type::execution_space execution_space;
149  typedef typename LWGraph_kokkos::local_graph_type::device_type::memory_space memory_space;
150 
151  *out << "Entering structured aggregation" << std::endl;
152 
153  ParameterList pL = GetParameterList();
154  bDefinitionPhase_ = false; // definition phase is finished, now all aggregation algorithm information is fixed
155 
156  // General problem informations are gathered from data stored in the problem matix.
157  RCP<const LWGraph_kokkos> graph = Get<RCP<LWGraph_kokkos> >(currentLevel, "Graph");
158  RCP<const Map> fineMap = graph->GetDomainMap();
159  const int myRank = fineMap->getComm()->getRank();
160  const LO dofsPerNode = Get<LO>(currentLevel, "DofsPerNode");
161 
162  // Since we want to operate on nodes and not dof, we need to modify the rowMap in order to
163  // obtain a nodeMap.
164  const int interpolationOrder = pL.get<int>("aggregation: coarsening order");
165  std::string outputType = pL.get<std::string>("aggregation: output type");
166  const bool outputAggregates = (outputType == "Aggregates" ? true : false);
167  Array<LO> lFineNodesPerDir(3);
168  int numDimensions;
169  if(currentLevel.GetLevelID() == 0) {
170  // On level 0, data is provided by applications and has no associated factory.
171  lFineNodesPerDir = currentLevel.Get<Array<LO> >("lNodesPerDim", NoFactory::get());
172  numDimensions = currentLevel.Get<int>("numDimensions", NoFactory::get());
173  } else {
174  // On level > 0, data is provided directly by generating factories.
175  lFineNodesPerDir = Get<Array<LO> >(currentLevel, "lNodesPerDim");
176  numDimensions = Get<int>(currentLevel, "numDimensions");
177  }
178 
179 
180  // First make sure that input parameters are set logically based on dimension
181  for(int dim = 0; dim < 3; ++dim) {
182  if(dim >= numDimensions) {
183  lFineNodesPerDir[dim] = 1;
184  }
185  }
186 
187  // Get the coarsening rate
188  std::string coarseningRate = pL.get<std::string>("aggregation: coarsening rate");
189  Teuchos::Array<LO> coarseRate;
190  try {
191  coarseRate = Teuchos::fromStringToArray<LO>(coarseningRate);
193  GetOStream(Errors,-1) << " *** \"aggregation: coarsening rate\" must be a string convertible into an array! *** "
194  << std::endl;
195  throw e;
196  }
197  TEUCHOS_TEST_FOR_EXCEPTION((coarseRate.size() > 1) && (coarseRate.size() < numDimensions),
198  Exceptions::RuntimeError,
199  "\"aggregation: coarsening rate\" must have at least as many"
200  " components as the number of spatial dimensions in the problem.");
201 
202  // Now that we have extracted info from the level, create the IndexManager
203  RCP<IndexManager_kokkos> geoData = rcp(new IndexManager_kokkos(numDimensions,
204  interpolationOrder, myRank,
205  lFineNodesPerDir,
206  coarseRate));
207 
208  *out << "The index manager has now been built" << std::endl;
209  TEUCHOS_TEST_FOR_EXCEPTION(fineMap->getNodeNumElements()
210  != static_cast<size_t>(geoData->getNumLocalFineNodes()),
211  Exceptions::RuntimeError,
212  "The local number of elements in the graph's map is not equal to "
213  "the number of nodes given by: lNodesPerDim!");
214 
215  // Now we are ready for the big loop over the fine node that will assign each
216  // node on the fine grid to an aggregate and a processor.
217  RCP<AggregationStructuredAlgorithm_kokkos> myStructuredAlgorithm
218  = rcp(new AggregationStructuredAlgorithm_kokkos());
219 
220  if(interpolationOrder == 0 && outputAggregates){
221  RCP<Aggregates_kokkos> aggregates = rcp(new Aggregates_kokkos(graph->GetDomainMap()));
222  aggregates->setObjectLabel("ST");
223  aggregates->SetIndexManager(geoData);
224  aggregates->AggregatesCrossProcessors(false);
225  aggregates->SetNumAggregates(geoData->getNumCoarseNodes());
226 
227  LO numNonAggregatedNodes = geoData->getNumLocalFineNodes();
228  Kokkos::View<unsigned*, memory_space> aggStat("aggStat", numNonAggregatedNodes);
229  Kokkos::parallel_for("StructuredAggregation: initialize aggStat",
230  Kokkos::RangePolicy<execution_space>(0, numNonAggregatedNodes),
231  KOKKOS_LAMBDA(const LO nodeIdx) {aggStat(nodeIdx) = READY;});
232 
233  myStructuredAlgorithm->BuildAggregates(pL, *graph, *aggregates, aggStat,
234  numNonAggregatedNodes);
235 
236  *out << "numNonAggregatedNodes: " << numNonAggregatedNodes << std::endl;
237 
238  TEUCHOS_TEST_FOR_EXCEPTION(numNonAggregatedNodes, Exceptions::RuntimeError,
239  "MueLu::StructuredAggregationFactory::Build: Leftover nodes found! Error!");
240  aggregates->ComputeAggregateSizes(true/*forceRecompute*/);
241  GetOStream(Statistics1) << aggregates->description() << std::endl;
242  Set(currentLevel, "Aggregates", aggregates);
243 
244  } else {
245  // Create Coarse Data
246  RCP<CrsGraph> myGraph;
247  myStructuredAlgorithm->BuildGraph(*graph, geoData, dofsPerNode, myGraph);
248  Set(currentLevel, "prolongatorGraph", myGraph);
249  }
250 
251  Set(currentLevel, "lCoarseNodesPerDim", geoData->getCoarseNodesPerDirArray());
252  Set(currentLevel, "indexManager", geoData);
253  Set(currentLevel, "interpolationOrder", interpolationOrder);
254  Set(currentLevel, "numDimensions", numDimensions);
255 
256  } // Build()
257 
258 } //namespace MueLu
259 
260 #endif // HAVE_MUELU_KOKKOS_REFACTOR
261 #endif /* MUELU_STRUCTUREDAGGREGATIONFACTORY_KOKKOS_DEF_HPP */
void parallel_for(const ExecPolicy &policy, const FunctorType &functor, const std::string &str="", typename Impl::enable_if< Kokkos::Impl::is_execution_policy< ExecPolicy >::value >::type *=0)
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
Print more statistics.
LocalOrdinal LO
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
#define SET_VALID_ENTRY(name)