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