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);
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);
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)
virtual void CallBuild(Level &requestedLevel) const =0
ParameterList & disableRecursiveValidation()
void CheckMainDiagonal(RCP< Matrix > &Ac) const
size_type size() const
void Release(const FactoryBase &factory)
Decrement the storage counter for all the inputs of a factory.
GlobalOrdinal GO
T & get(const std::string &name, T def_value)
ParameterList & set(std::string const &name, T const &value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
Timer to be used in factories. Similar to Monitor but with additional timers.
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
LocalOrdinal LO
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.
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
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)
static void TwoMatrixAdd(const Matrix &A, bool transposeA, SC alpha, Matrix &B, SC beta)
void AddTransferFactory(const RCP< const FactoryBase > &factory)
Add transfer factory in the end of list of transfer factories in RepartitionAcFactory.
static void Multiply(const Matrix &A, bool transposeA, const Matrix &B, bool transposeB, Matrix &C, bool call_FillComplete_on_result=true, bool doOptimizeStorage=true, const std::string &label=std::string(), const RCP< ParameterList > &params=null)
ArrayView< T > view(size_type lowerOffset, size_type size) const
void DeclareInput(Level &fineLevel, Level &coarseLevel) const
Input.
virtual std::string description() const
Return a simple one-line description of this object.
#define SET_VALID_ENTRY(name)