MueLu  Version of the Day
MueLu_StratimikosSmoother_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_STRATIMIKOSSMOOTHER_DEF_HPP
47 #define MUELU_STRATIMIKOSSMOOTHER_DEF_HPP
48 
49 #include "MueLu_ConfigDefs.hpp"
50 
51 #if defined(HAVE_MUELU_STRATIMIKOS) && defined(HAVE_MUELU_TPETRA) && defined(HAVE_MUELU_THYRA)
52 
53 #include <Teuchos_ParameterList.hpp>
54 
55 #include <Xpetra_CrsMatrix.hpp>
56 #include <Xpetra_CrsMatrixWrap.hpp>
57 #include <Xpetra_Matrix.hpp>
58 
60 #include "MueLu_Level.hpp"
62 #include "MueLu_Utilities.hpp"
63 #include "MueLu_Monitor.hpp"
64 #ifdef MUELU_RECURMG
66 #endif
67 
68 #include <Stratimikos_DefaultLinearSolverBuilder.hpp>
69 #include "Teuchos_AbstractFactoryStd.hpp"
70 #if defined(HAVE_MUELU_THYRATPETRAADAPTERS) && defined(HAVE_MUELU_IFPACK2)
71 #include <Thyra_Ifpack2PreconditionerFactory.hpp>
72 #endif
73 #include <Teuchos_ParameterList.hpp>
74 #include <unordered_map>
75 
76 
77 namespace MueLu {
78 
79  template <class LocalOrdinal, class GlobalOrdinal, class Node>
80  StratimikosSmoother<double, LocalOrdinal, GlobalOrdinal, Node>::StratimikosSmoother(const std::string type, const Teuchos::ParameterList& paramList)
81  : type_(type)
82  {
83  std::transform(type_.begin(), type_.end(), type_.begin(), ::toupper);
84  ParameterList& pList = const_cast<ParameterList&>(paramList);
85 
86  if (pList.isParameter("smoother: recurMgOnFilteredA")) {
87  recurMgOnFilteredA_ = true;
88  pList.remove("smoother: recurMgOnFilteredA");
89  }
90  bool isSupported = type_ == "STRATIMIKOS";
91  this->declareConstructionOutcome(!isSupported, "Stratimikos does not provide the smoother '" + type_ + "'.");
92  if (isSupported)
93  SetParameterList(paramList);
94  }
95 
96  template <class LocalOrdinal, class GlobalOrdinal, class Node>
97  void StratimikosSmoother<double, LocalOrdinal, GlobalOrdinal, Node>::SetParameterList(const Teuchos::ParameterList& paramList) {
98  Factory::SetParameterList(paramList);
99  }
100 
101  template <class LocalOrdinal, class GlobalOrdinal, class Node>
102  void StratimikosSmoother<double, LocalOrdinal, GlobalOrdinal, Node>::DeclareInput(Level& currentLevel) const {
103  this->Input(currentLevel, "A");
104  }
105 
106  template <class LocalOrdinal, class GlobalOrdinal, class Node>
107  void StratimikosSmoother<double, LocalOrdinal, GlobalOrdinal, Node>::Setup(Level& currentLevel) {
108  FactoryMonitor m(*this, "Setup Smoother", currentLevel);
109 
110  A_ = Factory::Get< RCP<Matrix> >(currentLevel, "A");
111  SetupStratimikos(currentLevel);
112  SmootherPrototype::IsSetup(true);
113  this->GetOStream(Statistics1) << description() << std::endl;
114  }
115 
116  template <class LocalOrdinal, class GlobalOrdinal, class Node>
117  void StratimikosSmoother<double, LocalOrdinal, GlobalOrdinal, Node>::SetupStratimikos(Level& currentLevel) {
118 
119  RCP<const Thyra::LinearOpBase<Scalar> > thyraA;
120  if (recurMgOnFilteredA_) {
121  RCP<Matrix> filteredA;
122  ExperimentalDropVertConnections(filteredA, currentLevel);
123  thyraA = Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node>::toThyra(Teuchos::rcp_dynamic_cast<CrsMatrixWrap>(filteredA)->getCrsMatrix());
124  }
125  else thyraA = Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node>::toThyra(Teuchos::rcp_dynamic_cast<CrsMatrixWrap>(A_)->getCrsMatrix());
126 
127  // Build Stratimikos solver
128  Stratimikos::DefaultLinearSolverBuilder linearSolverBuilder;
129  if (recurMgOnFilteredA_) {
130 #ifdef MUELU_RECURMG
131  Stratimikos::enableMueLu<LocalOrdinal,GlobalOrdinal,Node>(linearSolverBuilder);
132 #else
133  TEUCHOS_TEST_FOR_EXCEPTION( true , Exceptions::RuntimeError, "MueLu::StratimikosSmoother:: must compile with MUELU_RECURMG defined. Unfortunately, cmake does not always produce a proper link.txt file (which sometimes requires libmuelu.a before and after libmuelu-interface.a). After configuring, run script muelu/utils/misc/patchLinkForRecurMG to change link.txt files manually. If you want to create test example, add -DMUELU_RECURMG=ON to cmake arguments.");
134 #endif
135  }
136 
137  typedef Thyra::PreconditionerFactoryBase<Scalar> Base;
138 #if defined(HAVE_MUELU_THYRATPETRAADAPTERS) && defined(HAVE_MUELU_IFPACK2)
139  typedef Thyra::Ifpack2PreconditionerFactory<Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > Impl;
140  linearSolverBuilder.setPreconditioningStrategyFactory(Teuchos::abstractFactoryStd<Base, Impl>(), "Ifpack2");
141 #endif
142  linearSolverBuilder.setParameterList(rcpFromRef(const_cast<ParameterList&>(this->GetParameterList())));
143 
144 
145  // Build a new "solver factory" according to the previously specified parameter list.
146  RCP<Thyra::LinearOpWithSolveFactoryBase<Scalar> > solverFactory = Thyra::createLinearSolveStrategy(linearSolverBuilder);
147  solver_ = Thyra::linearOpWithSolve(*solverFactory, thyraA);
148 #ifdef dumpOutRecurMGDebug
149 char mystring[120];
150 sprintf(mystring,"for i in A_[0123456789].m P_[0123456789].m; do T=Xecho $i | sed Xs/.m$/%d.m/XX; mv $i $T; done", (int) currentLevel.GetLevelID()); fflush(stdout);
151 mystring[50]='`'; mystring[65]='"'; mystring[76]='"'; mystring[77]='`';
152 system(mystring);
153 #endif
154  }
155 
156  template <class LocalOrdinal, class GlobalOrdinal, class Node>
157  void StratimikosSmoother<double, LocalOrdinal, GlobalOrdinal, Node>::ExperimentalDropVertConnections(RCP<Matrix> & filteredA, Level& currentLevel) {
158 
159  // strip out the veritcal connections.
160  //
161  // There is some code, which is currently turned off (via sumDropped). That
162  // makes things more complicated. I want to keep it for now, and so it is
163  // explained here. The basic idea is to try and maintain the character
164  // of the horizontal stencil by summing dropped entries appropriately.
165  // However, this does not correspond to plane relexation, and so I am
166  // not sure it is really justified. Anyway, the picture below
167  // gives a situation with a 15-pt stencil
168  //
169  // a
170  // /
171  // /
172  // b ----c----- d
173  // /
174  // /
175  // e
176  // f dropped a & l summed into f
177  // / dropped b & m summed into g
178  // / dropped c & n summed into i
179  // g ----i----- j dropped d & o summed into j
180  // / dropped e & p summed into k
181  // /
182  // k
183  // l
184  // /
185  // /
186  // m ----n----- o
187  // /
188  // /
189  // p
190  // To do this, we use umap to record locations within the middle layer associated
191  // with each line ID (e.g. g corresponds to the 13th line). Then, when dropping (in a 2nd pass) we
192  // use lineId and umap to find where dropped entries should be summed (e.g., b corresponds to the
193  // 13th line and umap[13] points at location for g).
194  using TST = typename Teuchos::ScalarTraits<SC>;
195 
196  bool sumDropped = false;
197 
198  LO dofsPerNode = A_->GetFixedBlockSize();
199 
200 
201  RCP<ParameterList> fillCompleteParams(new ParameterList); // This code taken from Build method
202  fillCompleteParams->set("No Nonlocal Changes", true); // within MueLu_FilteredAFactory_def
203  filteredA = MatrixFactory::Build(A_->getCrsGraph());
204  filteredA->resumeFill();
205 
206  ArrayView<const LocalOrdinal> inds;
207  ArrayView<const Scalar> valsA;
208  ArrayView<Scalar> vals;
209  Teuchos::ArrayRCP<LocalOrdinal> TVertLineId = Factory::Get< Teuchos::ArrayRCP<LocalOrdinal> > (currentLevel, "LineDetection_VertLineIds");
210  Teuchos::ArrayRCP<LocalOrdinal> TLayerId = Factory::Get< Teuchos::ArrayRCP<LocalOrdinal> > (currentLevel, "LineDetection_Layers");
211  LocalOrdinal* VertLineId = TVertLineId.getRawPtr();
212  LocalOrdinal* LayerId = TLayerId.getRawPtr();
213  TEUCHOS_TEST_FOR_EXCEPTION( (LayerId == NULL) || (VertLineId == NULL), Exceptions::RuntimeError, "MueLu::StratimikosSmoother:: no line information found on this level. Cannot use recurMgOnFilteredA on this level.");
214 
215  Scalar ZERO = Teuchos::ScalarTraits<Scalar>::zero();
216  for (size_t i = 0; i < A_->getRowMap()->getNodeNumElements(); i++) {
217  A_->getLocalRowView(i, inds, valsA);
218  size_t nnz = inds.size();
219  ArrayView<const Scalar> vals1;
220  filteredA->getLocalRowView(i, inds, vals1);
221  vals = ArrayView<Scalar>(const_cast<Scalar*>(vals1.getRawPtr()), nnz);
222  memcpy(vals.getRawPtr(), valsA.getRawPtr(), nnz*sizeof(Scalar));
223  size_t inode, jdof, jnode, jdof_offset;
224  inode = i/dofsPerNode;
225 
226  std::unordered_map<LocalOrdinal, LocalOrdinal> umap; // umap[k] indicates where the dropped entry
227  // corresponding to kth line should be added
228  // within the row. See comments above.
229  if (sumDropped) {
230  for (size_t j = 0; j < nnz; j++) {
231  jdof = inds[j];
232  jnode = jdof/dofsPerNode;
233  jdof_offset = jdof - jnode*dofsPerNode;
234  if ( LayerId[jnode] == LayerId[inode] ) umap[dofsPerNode*VertLineId[jnode]+jdof_offset] = j;
235  }
236  }
237 
238  // drop non-middle layer entries. When sumDropped=true, sum dropped entries to corresponding mid-layer entry
239  for (size_t j = 0; j < nnz; j++) {
240  jdof = inds[j];
241  jnode = jdof/dofsPerNode;
242  jdof_offset = jdof - jnode*dofsPerNode;
243  if ( LayerId[jnode] != LayerId[inode] ) {
244  if (sumDropped) {
245  if (umap.find(dofsPerNode*VertLineId[jnode + jdof_offset]) != umap.end())
246  vals[umap[dofsPerNode*VertLineId[jnode + jdof_offset]]] += vals[j];
247  }
248  vals[j] = ZERO;
249  }
250  }
251  }
252  filteredA->fillComplete(fillCompleteParams);
253 
254  }
255 
256  template <class LocalOrdinal, class GlobalOrdinal, class Node>
257  void StratimikosSmoother<double, LocalOrdinal, GlobalOrdinal, Node>::Apply(MultiVector& X, const MultiVector& B, bool InitialGuessIsZero) const {
258  TEUCHOS_TEST_FOR_EXCEPTION(SmootherPrototype::IsSetup() == false, Exceptions::RuntimeError, "MueLu::StratimikosSmoother::Apply(): Setup() has not been called");
259 
260  // Apply
261  if (InitialGuessIsZero) {
262  X.putScalar(0.0);
263  RCP< Thyra::MultiVectorBase<Scalar> > thyraX = Teuchos::rcp_const_cast<Thyra::MultiVectorBase<Scalar> >(Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node>::toThyraMultiVector(rcpFromRef(X)));
264  RCP<const Thyra::MultiVectorBase<Scalar> > thyraB = Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node>::toThyraMultiVector(rcpFromRef(B));
265  Thyra::SolveStatus<Scalar> status = Thyra::solve<Scalar>(*solver_, Thyra::NOTRANS, *thyraB, thyraX.ptr());
266  RCP<MultiVector> thyXpX = Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node>::toXpetra(thyraX, X.getMap()->getComm());
267  X = *thyXpX;
268 
269  } else {
270  typedef Teuchos::ScalarTraits<Scalar> TST;
271  RCP<MultiVector> Residual = Utilities::Residual(*A_, X, B);
272 
273  RCP<MultiVector> Cor = Xpetra::MultiVectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(X.getMap(),X.getNumVectors(), true);
274  RCP< Thyra::MultiVectorBase<Scalar> > thyraCor = Teuchos::rcp_const_cast<Thyra::MultiVectorBase<Scalar> >(Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node>::toThyraMultiVector(Cor));
275  RCP<const Thyra::MultiVectorBase<Scalar> > thyraRes = Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node>::toThyraMultiVector(Residual);
276  Thyra::SolveStatus<Scalar> status = Thyra::solve<Scalar>(*solver_, Thyra::NOTRANS, *thyraRes, thyraCor.ptr());
277  RCP<MultiVector> thyXpCor = Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node>::toXpetra(thyraCor, X.getMap()->getComm());
278  X.update(TST::one(), *thyXpCor, TST::one());
279  }
280 
281  }
282 
283 
284  template <class LocalOrdinal, class GlobalOrdinal, class Node>
285  RCP<MueLu::SmootherPrototype<double, LocalOrdinal, GlobalOrdinal, Node> > StratimikosSmoother<double, LocalOrdinal, GlobalOrdinal, Node>::Copy() const {
286  RCP<StratimikosSmoother> smoother = rcp(new StratimikosSmoother(*this) );
287  smoother->SetParameterList(this->GetParameterList());
288  return smoother;
289  }
290 
291  template <class LocalOrdinal, class GlobalOrdinal, class Node>
292  std::string StratimikosSmoother<double, LocalOrdinal, GlobalOrdinal, Node>::description() const {
293  std::ostringstream out;
294  if (SmootherPrototype::IsSetup()) {
295  out << solver_->description();
296  } else {
297  out << "STRATIMIKOS {type = " << type_ << "}";
298  }
299  return out.str();
300  }
301 
302  template <class LocalOrdinal, class GlobalOrdinal, class Node>
303  void StratimikosSmoother<double, LocalOrdinal, GlobalOrdinal, Node>::print(Teuchos::FancyOStream &out, const VerbLevel verbLevel) const {
305 
306  if (verbLevel & Parameters1) {
307  out0 << "Parameter list: " << std::endl;
308  Teuchos::OSTab tab2(out);
309  out << this->GetParameterList();
310  }
311 
312  if (verbLevel & External)
313  if (solver_ != Teuchos::null) {
314  Teuchos::OSTab tab2(out);
315  out << *solver_ << std::endl;
316  }
317 
318  if (verbLevel & Debug) {
319  out0 << "IsSetup: " << Teuchos::toString(SmootherPrototype::IsSetup()) << std::endl
320  << "-" << std::endl
321  << "RCP<solver_>: " << solver_ << std::endl;
322  }
323  }
324 
325  template <class LocalOrdinal, class GlobalOrdinal, class Node>
326  size_t StratimikosSmoother<double, LocalOrdinal, GlobalOrdinal, Node>::getNodeSmootherComplexity() const {
327  return Teuchos::OrdinalTraits<size_t>::invalid();
328  }
329 
330 
331 } // namespace MueLu
332 
333 #endif // HAVE_MUELU_STRATIMIKOS
334 #endif // MUELU_STRATIMIKOSSMOOTHER_DEF_HPP
MueLu::DefaultLocalOrdinal LocalOrdinal
std::string toString(const T &what)
Little helper function to convert non-string types to strings.
Print external lib objects.
Print more statistics.
Print additional debugging information.
Namespace for MueLu classes and methods.
MueLu::DefaultScalar Scalar
#define MUELU_DESCRIBE
Helper macro for implementing Describable::describe() for BaseClass objects.
Print class parameters (more parameters, more verbose)