MueLu  Version of the Day
MueLu_RAPShiftFactory_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_RAPSHIFTFACTORY_DEF_HPP
47 #define MUELU_RAPSHIFTFACTORY_DEF_HPP
48 
49 #include <sstream>
50 
51 #include <Xpetra_Matrix.hpp>
52 #include <Xpetra_MatrixMatrix.hpp>
53 #include <Xpetra_Vector.hpp>
54 #include <Xpetra_VectorFactory.hpp>
55 
57 #include "MueLu_MasterList.hpp"
58 #include "MueLu_Monitor.hpp"
59 #include "MueLu_PerfUtils.hpp"
60 #include "MueLu_Utilities.hpp"
61 
62 namespace MueLu {
63 
64  /*********************************************************************************************************/
65  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
67  : implicitTranspose_(false), checkAc_(false), repairZeroDiagonals_(false) { }
68 
69 
70  /*********************************************************************************************************/
71  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
73  RCP<ParameterList> validParamList = rcp(new ParameterList());
74 
75 #define SET_VALID_ENTRY(name) validParamList->setEntry(name, MasterList::getEntry(name))
76  SET_VALID_ENTRY("transpose: use implicit");
77  SET_VALID_ENTRY("rap: shift");
78 #undef SET_VALID_ENTRY
79 
80  validParamList->set< RCP<const FactoryBase> >("M", Teuchos::null, "Generating factory of the matrix M used during the non-Galerkin RAP");
81  validParamList->set< RCP<const FactoryBase> >("K", Teuchos::null, "Generating factory of the matrix K used during the non-Galerkin RAP");
82  validParamList->set< RCP<const FactoryBase> >("P", Teuchos::null, "Prolongator factory");
83  validParamList->set< RCP<const FactoryBase> >("R", Teuchos::null, "Restrictor factory");
84 
85  validParamList->set< bool > ("CheckMainDiagonal", false, "Check main diagonal for zeros");
86  validParamList->set< bool > ("RepairMainDiagonal", false, "Repair zeros on main diagonal");
87 
88  // Make sure we don't recursively validate options for the matrixmatrix kernels
89  ParameterList norecurse;
90  norecurse.disableRecursiveValidation();
91  validParamList->set<ParameterList> ("matrixmatrix: kernel params", norecurse, "MatrixMatrix kernel parameters");
92 
93  return validParamList;
94  }
95 
96  /*********************************************************************************************************/
97  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
99  if (implicitTranspose_ == false) {
100  Input(coarseLevel, "R");
101  }
102 
103  Input(fineLevel, "K");
104  Input(fineLevel, "M");
105  Input(coarseLevel, "P");
106 
107  // call DeclareInput of all user-given transfer factories
108  for(std::vector<RCP<const FactoryBase> >::const_iterator it = transferFacts_.begin(); it!=transferFacts_.end(); ++it) {
109  (*it)->CallDeclareInput(coarseLevel);
110  }
111  }
112 
113  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
114  void RAPShiftFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Build(Level &fineLevel, Level &coarseLevel) const { // FIXME make fineLevel const
115  {
116  FactoryMonitor m(*this, "Computing Ac", coarseLevel);
117  const Teuchos::ParameterList& pL = GetParameterList();
118 
119  // Inputs: K, M, P
120  RCP<Matrix> K = Get< RCP<Matrix> >(fineLevel, "K");
121  RCP<Matrix> M = Get< RCP<Matrix> >(fineLevel, "M");
122  RCP<Matrix> P = Get< RCP<Matrix> >(coarseLevel, "P");
123 
124  // Build Kc = RKP, Mc = RMP
125  RCP<Matrix> KP, MP;
126 
127  // Reuse pattern if available (multiple solve)
128  // FIXME: Old style reuse doesn't work any kore
129  // if (IsAvailable(coarseLevel, "AP Pattern")) {
130  // KP = Get< RCP<Matrix> >(coarseLevel, "AP Pattern");
131  // MP = Get< RCP<Matrix> >(coarseLevel, "AP Pattern");
132  // }
133 
134  {
135  SubFactoryMonitor subM(*this, "MxM: K x P", coarseLevel);
136  KP = Xpetra::MatrixMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Multiply(*K, false, *P, false, KP, GetOStream(Statistics2));
137  MP = Xpetra::MatrixMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Multiply(*M, false, *P, false, MP, GetOStream(Statistics2));
138  Set(coarseLevel, "AP Pattern", KP);
139  }
140 
141  // Optimization storage option. If not modifying matrix later (inserting local values),
142  // allow optimization of storage. This is necessary for new faster Epetra MM kernels.
143  bool doOptimizedStorage = !checkAc_;
144 
145  RCP<Matrix> Ac, Kc, Mc;
146 
147  // Reuse pattern if available (multiple solve)
148  // if (IsAvailable(coarseLevel, "RAP Pattern"))
149  // Ac = Get< RCP<Matrix> >(coarseLevel, "RAP Pattern");
150 
151  bool doFillComplete=true;
152  if (implicitTranspose_) {
153  SubFactoryMonitor m2(*this, "MxM: P' x (KP) (implicit)", coarseLevel);
154  Kc = Xpetra::MatrixMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Multiply(*P, true, *KP, false, Kc, GetOStream(Statistics2), doFillComplete, doOptimizedStorage);
155  Mc = Xpetra::MatrixMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Multiply(*P, true, *MP, false, Mc, GetOStream(Statistics2), doFillComplete, doOptimizedStorage);
156  }
157  else {
158  RCP<Matrix> R = Get< RCP<Matrix> >(coarseLevel, "R");
159  SubFactoryMonitor m2(*this, "MxM: R x (KP) (explicit)", coarseLevel);
160  Kc = Xpetra::MatrixMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Multiply(*R, false, *KP, false, Kc, GetOStream(Statistics2), doFillComplete, doOptimizedStorage);
161  Mc = Xpetra::MatrixMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Multiply(*R, false, *MP, false, Mc, GetOStream(Statistics2), doFillComplete, doOptimizedStorage);
162  }
163 
164  // Get the shift
165  // FIXME - We should really get rid of the shifts array and drive this the same way everything else works
166  int level = coarseLevel.GetLevelID();
167  Scalar shift = Teuchos::ScalarTraits<Scalar>::zero();
168  if(level < (int)shifts_.size())
169  shift = shifts_[level];
170  else
171  shift = Teuchos::as<Scalar>(pL.get<double>("rap: shift"));
172 
173  // recombine to get K+shift*M
174  {
175  SubFactoryMonitor m2(*this, "Add: RKP + s*RMP", coarseLevel);
176  Xpetra::MatrixMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::TwoMatrixAdd(*Kc, false, (Scalar) 1.0, *Mc, false, shift, Ac, GetOStream(Statistics2));
177  Ac->fillComplete();
178  }
179 
180  if (checkAc_)
181  CheckMainDiagonal(Ac);
182 
183  RCP<ParameterList> params = rcp(new ParameterList());;
184  params->set("printLoadBalancingInfo", true);
185  GetOStream(Statistics0) << PerfUtils::PrintMatrixInfo(*Ac, "Ac", params);
186 
187  Set(coarseLevel, "A", Ac);
188  Set(coarseLevel, "K", Kc);
189  Set(coarseLevel, "M", Mc);
190  // Set(coarseLevel, "RAP Pattern", Ac);
191  }
192 
193  if (transferFacts_.begin() != transferFacts_.end()) {
194  SubFactoryMonitor m(*this, "Projections", coarseLevel);
195 
196  // call Build of all user-given transfer factories
197  for (std::vector<RCP<const FactoryBase> >::const_iterator it = transferFacts_.begin(); it != transferFacts_.end(); ++it) {
198  RCP<const FactoryBase> fac = *it;
199  GetOStream(Runtime0) << "RAPShiftFactory: call transfer factory: " << fac->description() << std::endl;
200  fac->CallBuild(coarseLevel);
201  // AP (11/11/13): I am not sure exactly why we need to call Release, but we do need it to get rid
202  // of dangling data for CoordinatesTransferFactory
203  coarseLevel.Release(*fac);
204  }
205  }
206  }
207 
208  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
210  // plausibility check: no zeros on diagonal
211  LO lZeroDiags = 0;
212  RCP<Vector> diagVec = VectorFactory::Build(Ac->getRowMap());
213  Ac->getLocalDiagCopy(*diagVec);
214  Teuchos::ArrayRCP< Scalar > diagVal = diagVec->getDataNonConst(0);
215  for (size_t r=0; r<Ac->getRowMap()->getNodeNumElements(); r++) {
216  if(diagVal[r]==0.0) {
217  lZeroDiags++;
218  if(repairZeroDiagonals_) {
219  GlobalOrdinal grid = Ac->getRowMap()->getGlobalElement(r);
220  LocalOrdinal lcid = Ac->getColMap()->getLocalElement(grid);
221  Teuchos::ArrayRCP<LocalOrdinal> indout(1,lcid);
222  Teuchos::ArrayRCP<Scalar> valout(1,Teuchos::ScalarTraits<Scalar>::one());
223  Ac->insertLocalValues(r, indout.view(0,indout.size()), valout.view(0,valout.size()));
224  }
225  }
226  }
227 
228  if(IsPrint(Warnings0)) {
229  const RCP<const Teuchos::Comm<int> > & comm = Ac->getRowMap()->getComm();
230  GO lZeroDiagsGO = lZeroDiags; /* LO->GO conversion */
231  GO gZeroDiags = 0;
232  MueLu_sumAll(comm, lZeroDiagsGO, gZeroDiags);
233  if(repairZeroDiagonals_) GetOStream(Warnings0) << "RAPShiftFactory (WARNING): repaired " << gZeroDiags << " zeros on main diagonal of Ac." << std::endl;
234  else GetOStream(Warnings0) << "RAPShiftFactory (WARNING): found " << gZeroDiags << " zeros on main diagonal of Ac." << std::endl;
235  }
236  }
237 
238  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
240  // check if it's a TwoLevelFactoryBase based transfer factory
241  TEUCHOS_TEST_FOR_EXCEPTION(Teuchos::rcp_dynamic_cast<const TwoLevelFactoryBase>(factory) == Teuchos::null, Exceptions::BadCast, "MueLu::RAPShiftFactory::AddTransferFactory: Transfer factory is not derived from TwoLevelFactoryBase. This is very strange. (Note: you can remove this exception if there's a good reason for)");
242  transferFacts_.push_back(factory);
243  }
244 
245 } //namespace MueLu
246 
247 #define MUELU_RAPSHIFTFACTORY_SHORT
248 #endif // MUELU_RAPSHIFTFACTORY_DEF_HPP
Important warning messages (one line)
Exception indicating invalid cast attempted.
#define MueLu_sumAll(rcpComm, in, out)
void CheckMainDiagonal(RCP< Matrix > &Ac) const
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.
One-liner description of what is happening.
void Build(Level &fineLevel, Level &coarseLevel) const
Build an object with this factory.
Namespace for MueLu classes and methods.
Print even more statistics.
Print statistics that do not involve significant additional computation.
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.
RCP< const ParameterList > GetValidParameterList() const
Return a const parameter list of valid parameters that setParameterList() will accept.
static std::string PrintMatrixInfo(const Matrix &A, const std::string &msgTag, RCP< const Teuchos::ParameterList > params=Teuchos::null)
void AddTransferFactory(const RCP< const FactoryBase > &factory)
Add transfer factory in the end of list of transfer factories in RepartitionAcFactory.
void DeclareInput(Level &fineLevel, Level &coarseLevel) const
Input.
#define SET_VALID_ENTRY(name)