MueLu Version of the Day
Loading...
Searching...
No Matches
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
63namespace 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()->getLocalNumElements();
131 for (LO row=0; row<(LO)P->getLocalNumRows(); 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()->getLocalNumElements(); col++) {
159 if (coarseDataGhosted[col] == LO_INVALID)
160 error_count++;
161 }
162 for (LO row=0; row<(LO)P->getLocalNumRows(); 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->getLocalElementList();
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
Timer to be used in factories. Similar to Monitor but with additional timers.
Class that holds all level-specific information.
Definition: MueLu_Level.hpp:99
bool IsAvailable(const std::string &ename, const FactoryBase *factory=NoFactory::get()) const
Test whether a need's value has been saved.
RequestMode GetRequestMode() const
void BuildFC(Level &fineLevel, Level &coarseLevel) const
RCP< const ParameterList > GetValidParameterList() const
Return a const parameter list of valid parameters that setParameterList() will accept.
void Build(Level &fineLevel, Level &coarseLevel) const
Build an object with this factory.
void BuildAggregates(Level &fineLevel, Level &coarseLevel) const
void DeclareInput(Level &finelevel, Level &coarseLevel) const
Specifies the data that this class needs, and the factories that generate that data.
Namespace for MueLu classes and methods.
@ Runtime0
One-liner description of what is happening.