MueLu  Version of the Day
MueLu_LocalOrdinalTransferFactory_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_LOCALORDINALTRANSFER_FACTORY_DEF_HPP
47 #define MUELU_LOCALORDINALTRANSFER_FACTORY_DEF_HPP
48 
49 #include "Xpetra_ImportFactory.hpp"
50 #include "Xpetra_VectorFactory.hpp"
51 #include "Xpetra_MapFactory.hpp"
52 #include "Xpetra_CrsGraph.hpp"
53 
54 #include "Xpetra_IO.hpp"
55 
56 #include "MueLu_CoarseMapFactory.hpp"
57 #include "MueLu_Aggregates.hpp"
59 
60 #include "MueLu_Level.hpp"
61 #include "MueLu_Monitor.hpp"
62 
63 namespace MueLu {
64 
65  template <class LocalOrdinal, class GlobalOrdinal, class Node>
67  RCP<ParameterList> validParamList = rcp(new ParameterList());
68 
69  validParamList->set<RCP<const FactoryBase> >(TransferVecName_, Teuchos::null, "Factory for TransferVec generation");
70  validParamList->set<RCP<const FactoryBase> >("P Graph", Teuchos::null, "Factory for P generation");
71  validParamList->set<RCP<const FactoryBase> >("Aggregates", Teuchos::null, "Factory for aggregates generation");
72  validParamList->set<RCP<const FactoryBase> >("CoarseMap", Teuchos::null, "Generating factory of the coarse map");
73 
74  return validParamList;
75  }
76 
77  template <class LocalOrdinal, class GlobalOrdinal, class Node>
79  static bool isAvailableXfer = false;
80  if (coarseLevel.GetRequestMode() == Level::REQUEST) {
81  isAvailableXfer = coarseLevel.IsAvailable(TransferVecName_, this);
82  if (isAvailableXfer == false) {
83  Input(fineLevel, TransferVecName_);
84  Input(fineLevel, "CoarseMap");
85 
86  if(useAggregatesMode_)
87  Input(fineLevel, "Aggregates");
88  else {
89  Input(coarseLevel, "P Graph");
90  }
91  }
92  }
93 
94  }
95 
96  template <class LocalOrdinal, class GlobalOrdinal, class Node>
98  if(useAggregatesMode_) BuildAggregates(fineLevel,coarseLevel);
99  else BuildFC(fineLevel,coarseLevel);
100  }
101 
102  template <class LocalOrdinal, class GlobalOrdinal, class Node>
104  FactoryMonitor m(*this, "Build", coarseLevel);
105 
106  GetOStream(Runtime0) << "Transferring " <<TransferVecName_ << std::endl;
107  LO LO_INVALID = Teuchos::OrdinalTraits<LO>::invalid();
108 
109  if (coarseLevel.IsAvailable(TransferVecName_, this)) {
110  GetOStream(Runtime0) << "Reusing "<<TransferVecName_ << std::endl;
111  return;
112  }
113 
114  // Get everything we need
115  RCP<const CrsGraph> P = Get< RCP<const CrsGraph> >(coarseLevel,"P Graph");
116  RCP<LocalOrdinalVector> fineTV = Get< RCP<LocalOrdinalVector> >(fineLevel, TransferVecName_);
117  RCP<const Map> coarseMap = Get< RCP<const Map> > (fineLevel, "CoarseMap");
118  RCP<const Map> uniqueMap = fineTV->getMap();
119  ArrayRCP<const LO> fineData = fineTV->getData(0);
120 
121  // Allocate new LO Vector
122  RCP<LocalOrdinalVector> coarseTV = LocalOrdinalVectorFactory::Build(coarseMap,1);
123  ArrayRCP<LO> coarseData = coarseTV->getDataNonConst(0);
124 
125  // Invalidate everything first, to check for errors
126  for(LO i=0; i<coarseData.size(); i++)
127  coarseData[i] = LO_INVALID;
128 
129  // Fill in coarse TV
130  LO domMapNumElements = P->getDomainMap()->getNodeNumElements();
131  for (LO row=0; row<(LO)P->getNodeNumRows(); row++) {
132  LO fineNumber = fineData[row];
133  ArrayView<const LO> indices;
134  P->getLocalRowView(row,indices);
135 
136  for(LO j=0; j<(LO)indices.size(); j++) {
137  LO col = indices[j];
138  if (col >= domMapNumElements) {
139  // skip off rank entries of P
140  } else {
141  coarseData[col] = fineNumber;
142  }
143  }
144  }
145 
146 #ifdef HAVE_MUELU_DEBUG
147  size_t error_count = 0;
148  {
149  RCP<LocalOrdinalVector> coarseTVghosted;
150  RCP<const Import> importer = P->getImporter();
151  if (!importer.is_null()) {
152  coarseTVghosted = LocalOrdinalVectorFactory::Build(P->getColMap(),1);
153  coarseTVghosted->doImport(*coarseTV, *importer, Xpetra::INSERT);
154  } else {
155  coarseTVghosted = coarseTV;
156  }
157  ArrayRCP<LO> coarseDataGhosted = coarseTVghosted->getDataNonConst(0);
158  for (LO col=0; col<(LO)P->getColMap()->getNodeNumElements(); col++) {
159  if (coarseDataGhosted[col] == LO_INVALID)
160  error_count++;
161  }
162  for (LO row=0; row<(LO)P->getNodeNumRows(); row++) {
163  LO fineNumber = fineData[row];
164  ArrayView<const LO> indices;
165  P->getLocalRowView(row,indices);
166  for(LO j=0; j<(LO)indices.size(); j++) {
167  if (coarseDataGhosted[indices[j]] != fineNumber)
168  error_count++;
169  }
170  }
171  }
172 
173  // Error checking: All nodes in an aggregate must share a local ordinal
174  if(error_count > 0) {
175  std::ostringstream ofs;
176  ofs << "LocalOrdinalTransferFactory("<<TransferVecName_<<"): ERROR: Each coarse dof must have a unique LO value. We had "<<std::to_string(error_count)<<" unknowns that did not match.";
177  throw std::runtime_error(ofs.str());
178  }
179 #endif
180 
181  Set<RCP<LocalOrdinalVector> >(coarseLevel, TransferVecName_, coarseTV);
182 
183  }
184 
185 
186 
187  template <class LocalOrdinal, class GlobalOrdinal, class Node>
189  FactoryMonitor m(*this, "Build", coarseLevel);
190 
191  GetOStream(Runtime0) << "Transferring " <<TransferVecName_ << std::endl;
192  RCP<LocalOrdinalVector> coarseTV;
193  RCP<LocalOrdinalVector> fineTV;
194  LO LO_INVALID = Teuchos::OrdinalTraits<LO>::invalid();
195 
196  if (coarseLevel.IsAvailable(TransferVecName_, this)) {
197  GetOStream(Runtime0) << "Reusing "<<TransferVecName_ << std::endl;
198  return;
199  }
200 
201  RCP<Aggregates> aggregates = Get< RCP<Aggregates> > (fineLevel, "Aggregates");
202  fineTV = Get< RCP<LocalOrdinalVector> >(fineLevel, TransferVecName_);
203  RCP<const Map> coarseMap = Get< RCP<const Map> > (fineLevel, "CoarseMap");
204  RCP<const Map> uniqueMap = fineTV->getMap();
205 
206  ArrayView<const GO> elementAList = coarseMap->getNodeElementList();
207 
208  coarseTV = LocalOrdinalVectorFactory::Build(coarseMap,1);
209 
210  // Create overlapped fine TV to reduce global communication
211  RCP<LocalOrdinalVector> ghostedTV = fineTV;
212  if (aggregates->AggregatesCrossProcessors()) {
213 
214  RCP<const Map> nonUniqueMap = aggregates->GetMap();
215  RCP<const Import> importer = ImportFactory::Build(uniqueMap, nonUniqueMap);
216 
217  ghostedTV = LocalOrdinalVectorFactory::Build(nonUniqueMap, 1);
218  ghostedTV->doImport(*fineTV, *importer, Xpetra::INSERT);
219  }
220 
221  // Get some info about aggregates
222  int myPID = uniqueMap->getComm()->getRank();
223  ArrayRCP<LO> aggSizes = aggregates->ComputeAggregateSizes();
224  const ArrayRCP<const LO> vertex2AggID = aggregates->GetVertex2AggId()->getData(0);
225  const ArrayRCP<const LO> procWinner = aggregates->GetProcWinner()->getData(0);
226 
227 
228  ArrayRCP<const LO> fineData = ghostedTV->getData(0);
229  ArrayRCP<LO> coarseData = coarseTV->getDataNonConst(0);
230 
231  // Invalidate everything first, to check for errors
232  for(LO i=0; i<coarseData.size(); i++)
233  coarseData[i] = LO_INVALID;
234 
235  // Fill in coarse TV
236  size_t error_count = 0;
237  for (LO lnode = 0; lnode < vertex2AggID.size(); lnode++) {
238  if (procWinner[lnode] == myPID &&
239  //lnode < vertex2AggID.size() &&
240  lnode < fineData.size() && // TAW do not access off-processor data
241  vertex2AggID[lnode] < coarseData.size()) {
242  if(coarseData[vertex2AggID[lnode]] == LO_INVALID)
243  coarseData[vertex2AggID[lnode]] = fineData[lnode];
244  if(coarseData[vertex2AggID[lnode]] != fineData[lnode])
245  error_count++;
246  }
247  }
248 
249  // Error checking: All nodes in an aggregate must share a local ordinal
250  if(error_count > 0) {
251  std::ostringstream ofs;
252  ofs << "LocalOrdinalTransferFactory: ERROR: Each aggregate must have a unique LO value. We had "<<std::to_string(error_count)<<" unknowns that did not match.";
253  throw std::runtime_error(ofs.str());
254  }
255 
256  Set<RCP<LocalOrdinalVector> >(coarseLevel, TransferVecName_, coarseTV);
257 
258  }
259 
260 } // namespace MueLu
261 
262 #endif // MUELU_LOCALORDINALTRANSFER_FACTORY_DEF_HPP
void Build(Level &fineLevel, Level &coarseLevel) const
Build an object with this factory.
RCP< const ParameterList > GetValidParameterList() const
Return a const parameter list of valid parameters that setParameterList() will accept.
RequestMode GetRequestMode() const
Timer to be used in factories. Similar to Monitor but with additional timers.
One-liner description of what is happening.
Namespace for MueLu classes and methods.
void DeclareInput(Level &finelevel, Level &coarseLevel) const
Specifies the data that this class needs, and the factories that generate that data.
Class that holds all level-specific information.
Definition: MueLu_Level.hpp:99
void BuildFC(Level &fineLevel, Level &coarseLevel) const
bool IsAvailable(const std::string &ename, const FactoryBase *factory=NoFactory::get()) const
Test whether a need&#39;s value has been saved.
void BuildAggregates(Level &fineLevel, Level &coarseLevel) const