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 
58 #include "MueLu_Monitor.hpp"
59 #include "MueLu_PerfUtils.hpp"
60 #include "MueLu_Utilities.hpp"
61 
62 namespace MueLu {
63 
64  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
66  : implicitTranspose_(false), checkAc_(false), repairZeroDiagonals_(false) { }
67 
68  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
70  if (implicitTranspose_ == false) {
71  Input(coarseLevel, "R");
72  }
73 
74  Input(fineLevel, "K");
75  Input(fineLevel, "M");
76  Input(coarseLevel, "P");
77 
78  // call DeclareInput of all user-given transfer factories
79  for(std::vector<RCP<const FactoryBase> >::const_iterator it = transferFacts_.begin(); it!=transferFacts_.end(); ++it) {
80  (*it)->CallDeclareInput(coarseLevel);
81  }
82  }
83 
84  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
85  void RAPShiftFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Build(Level &fineLevel, Level &coarseLevel) const { // FIXME make fineLevel const
86  {
87  FactoryMonitor m(*this, "Computing Ac", coarseLevel);
88 
89  // Inputs: K, M, P
90  RCP<Matrix> K = Get< RCP<Matrix> >(fineLevel, "K");
91  RCP<Matrix> M = Get< RCP<Matrix> >(fineLevel, "M");
92  RCP<Matrix> P = Get< RCP<Matrix> >(coarseLevel, "P");
93 
94  // Build Kc = RKP, Mc = RMP
95  RCP<Matrix> KP, MP;
96 
97  // Reuse pattern if available (multiple solve)
98  if (IsAvailable(coarseLevel, "AP Pattern")) {
99  KP = Get< RCP<Matrix> >(coarseLevel, "AP Pattern");
100  MP = Get< RCP<Matrix> >(coarseLevel, "AP Pattern");
101  }
102 
103  {
104  SubFactoryMonitor subM(*this, "MxM: K x P", coarseLevel);
105  KP = Xpetra::MatrixMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Multiply(*K, false, *P, false, KP, GetOStream(Statistics2));
106  MP = Xpetra::MatrixMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Multiply(*M, false, *P, false, MP, GetOStream(Statistics2));
107  Set(coarseLevel, "AP Pattern", KP);
108  }
109 
110  // Optimization storage option. If not modifying matrix later (inserting local values),
111  // allow optimization of storage. This is necessary for new faster Epetra MM kernels.
112  bool doOptimizedStorage = !checkAc_;
113 
114  RCP<Matrix> Ac, Kc, Mc;
115 
116  // Reuse pattern if available (multiple solve)
117  // if (IsAvailable(coarseLevel, "RAP Pattern"))
118  // Ac = Get< RCP<Matrix> >(coarseLevel, "RAP Pattern");
119 
120  bool doFillComplete=true;
121  if (implicitTranspose_) {
122  SubFactoryMonitor m2(*this, "MxM: P' x (KP) (implicit)", coarseLevel);
123  Kc = Xpetra::MatrixMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Multiply(*P, true, *KP, false, Kc, GetOStream(Statistics2), doFillComplete, doOptimizedStorage);
124  Mc = Xpetra::MatrixMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Multiply(*P, true, *MP, false, Mc, GetOStream(Statistics2), doFillComplete, doOptimizedStorage);
125  }
126  else {
127  RCP<Matrix> R = Get< RCP<Matrix> >(coarseLevel, "R");
128  SubFactoryMonitor m2(*this, "MxM: R x (KP) (explicit)", coarseLevel);
129  Kc = Xpetra::MatrixMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Multiply(*R, false, *KP, false, Kc, GetOStream(Statistics2), doFillComplete, doOptimizedStorage);
130  Mc = Xpetra::MatrixMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Multiply(*R, false, *MP, false, Mc, GetOStream(Statistics2), doFillComplete, doOptimizedStorage);
131  }
132 
133  // recombine to get K+shift*M
134  int level = coarseLevel.GetLevelID();
135  Scalar shift = shifts_[level];
136  Xpetra::MatrixMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::TwoMatrixAdd(*Kc, false, (Scalar) 1.0, *Mc, false, shift, Ac, GetOStream(Statistics2));
137  Ac->fillComplete();
138 
139  if (checkAc_)
140  CheckMainDiagonal(Ac);
141 
142  RCP<ParameterList> params = rcp(new ParameterList());;
143  params->set("printLoadBalancingInfo", true);
144  GetOStream(Statistics0) << PerfUtils::PrintMatrixInfo(*Ac, "Ac", params);
145 
146  Set(coarseLevel, "A", Ac);
147  Set(coarseLevel, "K", Kc);
148  Set(coarseLevel, "M", Mc);
149  Set(coarseLevel, "RAP Pattern", Ac);
150  }
151 
152  if (transferFacts_.begin() != transferFacts_.end()) {
153  SubFactoryMonitor m(*this, "Projections", coarseLevel);
154 
155  // call Build of all user-given transfer factories
156  for (std::vector<RCP<const FactoryBase> >::const_iterator it = transferFacts_.begin(); it != transferFacts_.end(); ++it) {
157  RCP<const FactoryBase> fac = *it;
158  GetOStream(Runtime0) << "RAPShiftFactory: call transfer factory: " << fac->description() << std::endl;
159  fac->CallBuild(coarseLevel);
160  // AP (11/11/13): I am not sure exactly why we need to call Release, but we do need it to get rid
161  // of dangling data for CoordinatesTransferFactory
162  coarseLevel.Release(*fac);
163  }
164  }
165  }
166 
167  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
169  // plausibility check: no zeros on diagonal
170  LO lZeroDiags = 0;
171  RCP<Vector> diagVec = VectorFactory::Build(Ac->getRowMap());
172  Ac->getLocalDiagCopy(*diagVec);
173  Teuchos::ArrayRCP< Scalar > diagVal = diagVec->getDataNonConst(0);
174  for (size_t r=0; r<Ac->getRowMap()->getNodeNumElements(); r++) {
175  if(diagVal[r]==0.0) {
176  lZeroDiags++;
177  if(repairZeroDiagonals_) {
178  GlobalOrdinal grid = Ac->getRowMap()->getGlobalElement(r);
179  LocalOrdinal lcid = Ac->getColMap()->getLocalElement(grid);
180  Teuchos::ArrayRCP<LocalOrdinal> indout(1,lcid);
181  Teuchos::ArrayRCP<Scalar> valout(1,Teuchos::ScalarTraits<Scalar>::one());
182  Ac->insertLocalValues(r, indout.view(0,indout.size()), valout.view(0,valout.size()));
183  }
184  }
185  }
186 
187  if(IsPrint(Warnings0)) {
188  const RCP<const Teuchos::Comm<int> > & comm = Ac->getRowMap()->getComm();
189  GO lZeroDiagsGO = lZeroDiags; /* LO->GO conversion */
190  GO gZeroDiags = 0;
191  MueLu_sumAll(comm, lZeroDiagsGO, gZeroDiags);
192  if(repairZeroDiagonals_) GetOStream(Warnings0) << "RAPShiftFactory (WARNING): repaired " << gZeroDiags << " zeros on main diagonal of Ac." << std::endl;
193  else GetOStream(Warnings0) << "RAPShiftFactory (WARNING): found " << gZeroDiags << " zeros on main diagonal of Ac." << std::endl;
194  }
195  }
196 
197  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
199  // check if it's a TwoLevelFactoryBase based transfer factory
200  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)");
201  transferFacts_.push_back(factory);
202  }
203 
204 } //namespace MueLu
205 
206 #define MUELU_RAPSHIFTFACTORY_SHORT
207 #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.
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.