MueLu  Version of the Day
MueLu_RAPFactory_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_RAPFACTORY_DEF_HPP
47 #define MUELU_RAPFACTORY_DEF_HPP
48 
49 
50 #include <sstream>
51 
52 #include <Xpetra_Matrix.hpp>
53 #include <Xpetra_MatrixFactory.hpp>
54 #include <Xpetra_MatrixMatrix.hpp>
55 #include <Xpetra_Vector.hpp>
56 #include <Xpetra_VectorFactory.hpp>
57 
59 
60 #include "MueLu_MasterList.hpp"
61 #include "MueLu_Monitor.hpp"
62 #include "MueLu_PerfUtils.hpp"
64 //#include "MueLu_Utilities.hpp"
65 
66 namespace MueLu {
67 
68  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
70  : hasDeclaredInput_(false) { }
71 
72  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
74  RCP<ParameterList> validParamList = rcp(new ParameterList());
75 
76 #define SET_VALID_ENTRY(name) validParamList->setEntry(name, MasterList::getEntry(name))
77  SET_VALID_ENTRY("transpose: use implicit");
78 #undef SET_VALID_ENTRY
79  validParamList->set< RCP<const FactoryBase> >("A", null, "Generating factory of the matrix A used during the prolongator smoothing process");
80  validParamList->set< RCP<const FactoryBase> >("P", null, "Prolongator factory");
81  validParamList->set< RCP<const FactoryBase> >("R", null, "Restrictor factory");
82 
83  validParamList->set< bool > ("CheckMainDiagonal", false, "Check main diagonal for zeros");
84  validParamList->set< bool > ("RepairMainDiagonal", false, "Repair zeros on main diagonal");
85 
86  return validParamList;
87  }
88 
89  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
91  const Teuchos::ParameterList& pL = GetParameterList();
92  if (pL.get<bool>("transpose: use implicit") == false)
93  Input(coarseLevel, "R");
94 
95  Input(fineLevel, "A");
96  Input(coarseLevel, "P");
97 
98  // call DeclareInput of all user-given transfer factories
99  for (std::vector<RCP<const FactoryBase> >::const_iterator it = transferFacts_.begin(); it != transferFacts_.end(); ++it)
100  (*it)->CallDeclareInput(coarseLevel);
101 
102  hasDeclaredInput_ = true;
103  }
104 
105  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
107  const bool doTranspose = true;
108  const bool doFillComplete = true;
109  const bool doOptimizeStorage = true;
110  {
111  FactoryMonitor m(*this, "Computing Ac", coarseLevel);
112  std::ostringstream levelstr;
113  levelstr << coarseLevel.GetLevelID();
114 
115  TEUCHOS_TEST_FOR_EXCEPTION(hasDeclaredInput_ == false, Exceptions::RuntimeError,
116  "MueLu::RAPFactory::Build(): CallDeclareInput has not been called before Build!");
117 
118  const Teuchos::ParameterList& pL = GetParameterList();
119  RCP<Matrix> A = Get< RCP<Matrix> >(fineLevel, "A");
120  RCP<Matrix> P = Get< RCP<Matrix> >(coarseLevel, "P"), AP, Ac;
121 
122  // Reuse pattern if available (multiple solve)
123  RCP<ParameterList> APparams = rcp(new ParameterList);
124  if (coarseLevel.IsAvailable("AP reuse data", this)) {
125  GetOStream(static_cast<MsgType>(Runtime0 | Test)) << "Reusing previous AP data" << std::endl;
126 
127  APparams = coarseLevel.Get< RCP<ParameterList> >("AP reuse data", this);
128 
129  if (APparams->isParameter("graph"))
130  AP = APparams->get< RCP<Matrix> >("graph");
131  }
132 
133  {
134  SubFactoryMonitor subM(*this, "MxM: A x P", coarseLevel);
135 
136  AP = MatrixMatrix::Multiply(*A, !doTranspose, *P, !doTranspose, AP, GetOStream(Statistics2),
137  doFillComplete, doOptimizeStorage, std::string("MueLu::A*P-")+levelstr.str(), APparams);
138  }
139 
140  // Reuse coarse matrix memory if available (multiple solve)
141  RCP<ParameterList> RAPparams = rcp(new ParameterList);
142  if (coarseLevel.IsAvailable("RAP reuse data", this)) {
143  GetOStream(static_cast<MsgType>(Runtime0 | Test)) << "Reusing previous RAP data" << std::endl;
144 
145  RAPparams = coarseLevel.Get< RCP<ParameterList> >("RAP reuse data", this);
146 
147  if (RAPparams->isParameter("graph"))
148  Ac = RAPparams->get< RCP<Matrix> >("graph");
149 
150  // Some eigenvalue may have been cached with the matrix in the previous run.
151  // As the matrix values will be updated, we need to reset the eigenvalue.
152  Ac->SetMaxEigenvalueEstimate(-Teuchos::ScalarTraits<SC>::one());
153  }
154 
155  // Allow optimization of storage.
156  // This is necessary for new faster Epetra MM kernels.
157  // Seems to work with matrix modifications to repair diagonal entries.
158 
159  if (pL.get<bool>("transpose: use implicit") == true) {
160  SubFactoryMonitor m2(*this, "MxM: P' x (AP) (implicit)", coarseLevel);
161 
162  Ac = MatrixMatrix::Multiply(*P, doTranspose, *AP, !doTranspose, Ac, GetOStream(Statistics2),
163  doFillComplete, doOptimizeStorage, std::string("MueLu::R*(AP)-implicit-")+levelstr.str(), RAPparams);
164 
165  } else {
166  RCP<Matrix> R = Get< RCP<Matrix> >(coarseLevel, "R");
167 
168  SubFactoryMonitor m2(*this, "MxM: R x (AP) (explicit)", coarseLevel);
169 
170  Ac = MatrixMatrix::Multiply(*R, !doTranspose, *AP, !doTranspose, Ac, GetOStream(Statistics2),
171  doFillComplete, doOptimizeStorage, std::string("MueLu::R*(AP)-explicit-")+levelstr.str(), RAPparams);
172  }
173 
174  CheckRepairMainDiagonal(Ac);
175 
176  if (IsPrint(Statistics1)) {
177  RCP<ParameterList> params = rcp(new ParameterList());;
178  params->set("printLoadBalancingInfo", true);
179  params->set("printCommInfo", true);
180  GetOStream(Statistics1) << PerfUtils::PrintMatrixInfo(*Ac, "Ac", params);
181  }
182 
183  Set(coarseLevel, "A", Ac);
184 
185  APparams->set("graph", AP);
186  Set(coarseLevel, "AP reuse data", APparams);
187  RAPparams->set("graph", Ac);
188  Set(coarseLevel, "RAP reuse data", RAPparams);
189  }
190 
191  if (transferFacts_.begin() != transferFacts_.end()) {
192  SubFactoryMonitor m(*this, "Projections", coarseLevel);
193 
194  // call Build of all user-given transfer factories
195  for (std::vector<RCP<const FactoryBase> >::const_iterator it = transferFacts_.begin(); it != transferFacts_.end(); ++it) {
196  RCP<const FactoryBase> fac = *it;
197  GetOStream(Runtime0) << "RAPFactory: call transfer factory: " << fac->description() << std::endl;
198  fac->CallBuild(coarseLevel);
199  // Coordinates transfer is marginally different from all other operations
200  // because it is *optional*, and not required. For instance, we may need
201  // coordinates only on level 4 if we start repartitioning from that level,
202  // but we don't need them on level 1,2,3. As our current Hierarchy setup
203  // assumes propagation of dependencies only through three levels, this
204  // means that we need to rely on other methods to propagate optional data.
205  //
206  // The method currently used is through RAP transfer factories, which are
207  // simply factories which are called at the end of RAP with a single goal:
208  // transfer some fine data to coarser level. Because these factories are
209  // kind of outside of the mainline factories, they behave different. In
210  // particular, we call their Build method explicitly, rather than through
211  // Get calls. This difference is significant, as the Get call is smart
212  // enough to know when to release all factory dependencies, and Build is
213  // dumb. This led to the following CoordinatesTransferFactory sequence:
214  // 1. Request level 0
215  // 2. Request level 1
216  // 3. Request level 0
217  // 4. Release level 0
218  // 5. Release level 1
219  //
220  // The problem is missing "6. Release level 0". Because it was missing,
221  // we had outstanding request on "Coordinates", "Aggregates" and
222  // "CoarseMap" on level 0.
223  //
224  // This was fixed by explicitly calling Release on transfer factories in
225  // RAPFactory. I am still unsure how exactly it works, but now we have
226  // clear data requests for all levels.
227  coarseLevel.Release(*fac);
228  }
229  }
230 
231  }
232 
233  // Plausibility check: no zeros on diagonal
234  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
236  const Teuchos::ParameterList& pL = GetParameterList();
237  bool repairZeroDiagonals = pL.get<bool>("RepairMainDiagonal");
238  bool checkAc = pL.get<bool>("CheckMainDiagonal");
239 
240  if (!checkAc && !repairZeroDiagonals)
241  return;
242 
243  SC zero = Teuchos::ScalarTraits<SC>::zero(), one = Teuchos::ScalarTraits<SC>::one();
244 
245  Teuchos::RCP<Teuchos::ParameterList> p = Teuchos::rcp(new Teuchos::ParameterList());
246  p->set("DoOptimizeStorage", true);
247 
248  RCP<const Map> rowMap = Ac->getRowMap();
249  RCP<Vector> diagVec = VectorFactory::Build(rowMap);
250  Ac->getLocalDiagCopy(*diagVec);
251 
252  LO lZeroDiags = 0;
253  Teuchos::ArrayRCP< Scalar > diagVal = diagVec->getDataNonConst(0);
254 
255  for (size_t i = 0; i < rowMap->getNodeNumElements(); i++) {
256  if (diagVal[i] == zero) {
257  lZeroDiags++;
258  }
259  }
260  GO gZeroDiags;
261  MueLu_sumAll(rowMap->getComm(), Teuchos::as<GO>(lZeroDiags), gZeroDiags);
262 
263  if (repairZeroDiagonals && gZeroDiags > 0) {
264  // TAW: If Ac has empty rows, put a 1 on the diagonal of Ac. Be aware that Ac might have empty rows AND columns.
265  // The columns might not exist in the column map at all.
266  //
267  // It would be nice to add the entries to the original matrix Ac. But then we would have to use
268  // insertLocalValues. However we cannot add new entries for local column indices that do not exist in the column map
269  // of Ac (at least Epetra is not able to do this). With Tpetra it is also not possible to add new entries after the
270  // FillComplete call with a static map, since the column map already exists and the diagonal entries are completely missing.
271  //
272  // Here we build a diagonal matrix with zeros on the diagonal and ones on the diagonal for the rows where Ac has empty rows
273  // We have to build a new matrix to be able to use insertGlobalValues. Then we add the original matrix Ac to our new block
274  // diagonal matrix and use the result as new (non-singular) matrix Ac.
275  // This is very inefficient.
276  //
277  // If you know something better, please let me know.
278  RCP<Matrix> fixDiagMatrix = MatrixFactory::Build(rowMap, 1);
279  for (size_t r = 0; r < rowMap->getNodeNumElements(); r++) {
280  if (diagVal[r] == zero) {
281  GO grid = rowMap->getGlobalElement(r);
282  Teuchos::ArrayRCP<GO> indout(1,grid);
283  Teuchos::ArrayRCP<SC> valout(1, one);
284  fixDiagMatrix->insertGlobalValues(grid,indout.view(0, 1), valout.view(0, 1));
285  }
286  }
287  {
288  Teuchos::TimeMonitor m1(*Teuchos::TimeMonitor::getNewTimer("CheckRepairMainDiagonal: fillComplete1"));
289  if(rowMap->lib() == Xpetra::UseTpetra) Ac->resumeFill(); // TODO needed for refactored Tpetra because of the isFillActive flag???
290  Ac->fillComplete(p);
291  }
292  MatrixMatrix::TwoMatrixAdd(*Ac, false, 1.0, *fixDiagMatrix, 1.0);
293  if (Ac->IsView("stridedMaps"))
294  fixDiagMatrix->CreateView("stridedMaps", Ac);
295 
296  Ac = Teuchos::null; // free singular coarse level matrix
297  Ac = fixDiagMatrix; // set fixed non-singular coarse level matrix
298 
299  // call fillComplete with optimized storage option set to true
300  // This is necessary for new faster Epetra MM kernels.
301  {
302  Teuchos::TimeMonitor m1(*Teuchos::TimeMonitor::getNewTimer("CheckRepairMainDiagonal: fillComplete2"));
303  if(rowMap->lib() == Xpetra::UseTpetra) Ac->resumeFill(); // TODO needed for refactored Tpetra because of the isFillActive flag???
304  Ac->fillComplete(p);
305  }
306  } // end repair
307 
308 
309 
310  // print some output
311  if (IsPrint(Warnings0))
312  GetOStream(Warnings0) << "RAPFactory (WARNING): " << (repairZeroDiagonals ? "repaired " : "found ")
313  << gZeroDiags << " zeros on main diagonal of Ac." << std::endl;
314 
315 #ifdef HAVE_MUELU_DEBUG // only for debugging
316  // check whether Ac has been repaired...
317  Ac->getLocalDiagCopy(*diagVec);
318  Teuchos::ArrayRCP< Scalar > diagVal2 = diagVec->getDataNonConst(0);
319  for (size_t r = 0; r < Ac->getRowMap()->getNodeNumElements(); r++) {
320  if (diagVal2[r] == zero) {
321  GetOStream(Errors,-1) << "Error: there are zeros left on diagonal after repair..." << std::endl;
322  break;
323  }
324  }
325 #endif
326  }
327 
328  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
330  // check if it's a TwoLevelFactoryBase based transfer factory
331  TEUCHOS_TEST_FOR_EXCEPTION(Teuchos::rcp_dynamic_cast<const TwoLevelFactoryBase>(factory) == Teuchos::null, Exceptions::BadCast,
332  "MueLu::RAPFactory::AddTransferFactory: Transfer factory is not derived from TwoLevelFactoryBase. "
333  "This is very strange. (Note: you can remove this exception if there's a good reason for)");
334  TEUCHOS_TEST_FOR_EXCEPTION(hasDeclaredInput_, Exceptions::RuntimeError, "MueLu::RAPFactory::AddTransferFactory: Factory is being added after we have already declared input");
335  transferFacts_.push_back(factory);
336  }
337 
338 } //namespace MueLu
339 
340 #define MUELU_RAPFACTORY_SHORT
341 #endif // MUELU_RAPFACTORY_DEF_HPP
Important warning messages (one line)
#define SET_VALID_ENTRY(name)
Exception indicating invalid cast attempted.
#define MueLu_sumAll(rcpComm, in, out)
void Release(const FactoryBase &factory)
Decrement the storage counter for all the inputs of a factory.
GlobalOrdinal GO
Timer to be used in factories. Similar to Monitor but with additional timers.
Print more statistics.
LocalOrdinal LO
One-liner description of what is happening.
void AddTransferFactory(const RCP< const FactoryBase > &factory)
Add transfer factory in the end of list of transfer factories in RepartitionAcFactory.
Namespace for MueLu classes and methods.
void CheckRepairMainDiagonal(RCP< Matrix > &Ac) const
Print skeleton for the run, i.e. factory calls and used parameters.
Print even more statistics.
int GetLevelID() const
Return level number.
Definition: MueLu_Level.cpp:76
RCP< const ParameterList > GetValidParameterList() const
Return a const parameter list of valid parameters that setParameterList() will accept.
Class that holds all level-specific information.
Definition: MueLu_Level.hpp:99
Timer to be used in factories. Similar to SubMonitor but adds a timer level by level.
static std::string PrintMatrixInfo(const Matrix &A, const std::string &msgTag, RCP< const Teuchos::ParameterList > params=Teuchos::null)
Scalar SC
void DeclareInput(Level &fineLevel, Level &coarseLevel) const
Input.
Exception throws to report errors in the internal logical of the program.
void Build(Level &fineLevel, Level &coarseLevel) const
Build an object with this factory.