MueLu  Version of the Day
MueLu_RefMaxwell_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_REFMAXWELL_DEF_HPP
47 #define MUELU_REFMAXWELL_DEF_HPP
48 
49 #include "MueLu_ConfigDefs.hpp"
50 
51 #include "Xpetra_MatrixMatrix.hpp"
52 
53 #if defined(HAVE_MUELU_TPETRA) && defined(HAVE_MUELU_IFPACK2)
54 
56 #include "MueLu_Utilities.hpp"
57 #include "MueLu_Monitor.hpp"
58 #include "MueLu_MLParameterListInterpreter.hpp"
59 #include "MueLu_ParameterListInterpreter.hpp"
61 
62 namespace MueLu {
63 
64 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
65 Teuchos::RCP<const Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > RefMaxwell<Scalar,LocalOrdinal,GlobalOrdinal,Node>::getDomainMap() const {
66 
67  return Xpetra::toTpetraNonZero(SM_Matrix_->getDomainMap());
68 
69 }
70 
71 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
72 Teuchos::RCP<const Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > RefMaxwell<Scalar,LocalOrdinal,GlobalOrdinal,Node>::getRangeMap() const {
73 
74  return Xpetra::toTpetraNonZero(SM_Matrix_->getRangeMap());
75 
76 }
77 
78 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
80 
81  disable_addon_ = list.get("refmaxwell: disable add-on",true);
82  mode_ = list.get("refmaxwell: mode","additive");
83 
84  if(list.isSublist("refmaxwell: 11list"))
85  precList11_ = list.sublist("refmaxwell: 11list");
86 
87  if(list.isSublist("refmaxwell: 22list"))
88  precList22_ = list.sublist("refmaxwell: 22list");
89 
90  std::string ref("smoother:");
91  std::string replace("coarse:");
92  for(Teuchos::ParameterList::ConstIterator i=list.begin(); i !=list.end(); i++) {
93  const std::string & pname = list.name(i);
94  if(pname.find(ref)!=std::string::npos) {
95  smootherList_.setEntry(pname,list.entry(i));
96  std::string coarsename(pname);
97  coarsename.replace((size_t)0,(size_t)ref.length(),replace);
98  }
99  }
100  if(list.isSublist("smoother: params")) {
101  smootherList_.set("coarse: params",list.sublist("smoother: params"));
102  }
103  smootherList_.set("max levels",1);
104 
105 }
106 
107 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
109 
110  Teuchos::FancyOStream out(Teuchos::rcpFromRef(std::cout));
111  out.setOutputToRootOnly(0);
112  out.setShowProcRank(true);
113 
114  // clean rows associated with boundary conditions
115  findDirichletRows(SM_Matrix_,BCrows_);
116  findDirichletCols(D0_Matrix_,BCrows_,BCcols_);
117  D0_Matrix_->resumeFill();
118  Apply_BCsToMatrixRows(D0_Matrix_,BCrows_);
119  Apply_BCsToMatrixCols(D0_Matrix_,BCcols_);
120  D0_Matrix_->fillComplete(D0_Matrix_->getDomainMap(),D0_Matrix_->getRangeMap());
121  //D0_Matrix_->describe(out,Teuchos::VERB_EXTREME);
122 
123  // Form TMT_Matrix
124  Teuchos::RCP<XMat> C1 = MatrixFactory::Build(SM_Matrix_->getRowMap(),0);
125  TMT_Matrix_=MatrixFactory::Build(D0_Matrix_->getDomainMap(),0);
126  Xpetra::MatrixMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Multiply(*SM_Matrix_,false,*D0_Matrix_,false,*C1,true,true);
127  Xpetra::MatrixMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Multiply(*D0_Matrix_,true,*C1,false,*TMT_Matrix_,true,true);
128  TMT_Matrix_->resumeFill();
129  Remove_Zeroed_Rows(TMT_Matrix_,1.0e-16);
130  TMT_Matrix_->SetFixedBlockSize(1);
131  //TMT_Matrix_->describe(out,Teuchos::VERB_EXTREME);
132 
133  // build nullspace if necessary
134  if(Nullspace_ != Teuchos::null) {
135  // no need to do anything - nullspace is built
136  }
137  else if(Nullspace_ == Teuchos::null && Coords_ != Teuchos::null) {
138  Nullspace_ = MultiVectorFactory::Build(SM_Matrix_->getRowMap(),Coords_->getNumVectors());
139  D0_Matrix_->apply(*Coords_,*Nullspace_);
140  }
141  else {
142  std::cerr << "MueLu::RefMaxwell::compute(): either the nullspace or the nodal coordinates must be provided." << std::endl;
143  }
144 
145  // build special prolongator for (1,1)-block
146  if(P11_==Teuchos::null) {
147  buildProlongator();
148  }
149 
150  // build coarse grid operator for (1,1)-block
151  formCoarseMatrix();
152 
153  // build fine grid operator for (2,2)-block, D0* M1 D0
154  Teuchos::RCP<XMat> C = MatrixFactory::Build(M1_Matrix_->getRowMap(),0);
155  Xpetra::MatrixMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Multiply(*M1_Matrix_,false,*D0_Matrix_,false,*C,true,true);
156  A22_=MatrixFactory::Build(D0_Matrix_->getDomainMap(),0);
157  Xpetra::MatrixMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Multiply(*D0_Matrix_,true,*C,false,*A22_,true,true);
158  A22_->resumeFill();
159  Remove_Zeroed_Rows(A22_,1.0e-16);
160  A22_->SetFixedBlockSize(1);
161 
162  // Use HierarchyManagers to build 11 & 22 Hierarchies
164  RCP<HierarchyManager> Manager11, Manager22, ManagerSmoother;
165  std::string syntaxStr = "parameterlist: syntax";
166  if (parameterList_.isParameter(syntaxStr) && parameterList_.get<std::string>(syntaxStr) == "ml") {
167  parameterList_.remove(syntaxStr);
168  Manager11 = rcp(new MLParameterListInterpreter<SC,LO,GO,NO>(precList11_));
169  Manager22 = rcp(new MLParameterListInterpreter<SC,LO,GO,NO>(precList22_));
170  ManagerSmoother = rcp(new MLParameterListInterpreter<SC,LO,GO,NO>(smootherList_));
171  } else {
172  Manager11 = rcp(new ParameterListInterpreter <SC,LO,GO,NO>(precList11_,A11_->getDomainMap()->getComm()));
173  Manager22 = rcp(new ParameterListInterpreter <SC,LO,GO,NO>(precList22_,A22_->getDomainMap()->getComm()));
174  ManagerSmoother = rcp(new ParameterListInterpreter<SC,LO,GO,NO>(smootherList_,SM_Matrix_->getDomainMap()->getComm()));
175  }
176 
177  Hierarchy11_=Manager11->CreateHierarchy();
178  Hierarchy11_->setlib(Xpetra::UseTpetra);
179  Hierarchy11_->GetLevel(0)->Set("A", A11_);
180  Manager11->SetupHierarchy(*Hierarchy11_);
181 
182  Hierarchy22_=Manager22->CreateHierarchy();
183  Hierarchy22_->setlib(Xpetra::UseTpetra);
184  Hierarchy22_->GetLevel(0)->Set("A", A22_);
185  Manager22->SetupHierarchy(*Hierarchy22_);
186 
187  // build ifpack2 preconditioners for pre and post smoothing
188  HierarchySmoother_=ManagerSmoother->CreateHierarchy();
189  HierarchySmoother_->setlib(Xpetra::UseTpetra);
190  HierarchySmoother_->GetLevel(0)->Set("A", SM_Matrix_);
191  ManagerSmoother->SetupHierarchy(*HierarchySmoother_);
192 
193 }
194 
195 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
197 
198  Teuchos::FancyOStream out(Teuchos::rcpFromRef(std::cout));
199  out.setOutputToRootOnly(0);
200  out.setShowProcRank(true);
201 
202  // build prolongator: algorithm 1 in the reference paper
203  // First, aggregate nodal matrix by creating a 2-level hierarchy
204  Teuchos::RCP<Hierarchy> auxHierarchy
205  = Teuchos::rcp( new Hierarchy(TMT_Matrix_) );
206  Teuchos::RCP<FactoryManager> auxManager
207  = Teuchos::rcp( new FactoryManager );
208  Teuchos::RCP<TentativePFactory> TentPfact
209  = Teuchos::rcp( new TentativePFactory );
210  Teuchos::RCP<SaPFactory> Pfact
211  = Teuchos::rcp( new SaPFactory );
212  Teuchos::RCP<UncoupledAggregationFactory> Aggfact
213  = Teuchos::rcp( new UncoupledAggregationFactory() );
214  Teuchos::ParameterList params;
215  params.set("sa: damping factor",0.0);
216  Pfact -> SetParameterList(params);
217  auxManager -> SetFactory("P", Pfact);
218  auxManager -> SetFactory("Ptent", TentPfact);
219  auxManager -> SetFactory("Aggregates", Aggfact);
220  auxManager -> SetFactory("Smoother", Teuchos::null);
221  auxManager -> SetFactory("CoarseSolver", Teuchos::null);
222  auxHierarchy -> Keep("P", Pfact.get());
223  auxHierarchy -> SetMaxCoarseSize(1);
224  auxHierarchy -> Setup(*auxManager, 0, 2);
225 
226  // pull out tentative P
227  Teuchos::RCP<Level> Level1 = auxHierarchy -> GetLevel(1);
228  Teuchos::RCP<XMat> P = Level1 -> Get< Teuchos::RCP<XMat> >("P",Pfact.get());
229 
230  // make weighting matrix
231  Teuchos::RCP<XMat> D0_Matrix_Abs=MatrixFactory2::BuildCopy(D0_Matrix_);
232  D0_Matrix_Abs -> resumeFill();
233  D0_Matrix_Abs -> setAllToScalar((Scalar)0.5);
234  Apply_BCsToMatrixRows(D0_Matrix_Abs,BCrows_);
235  Apply_BCsToMatrixCols(D0_Matrix_Abs,BCcols_);
236  D0_Matrix_Abs -> fillComplete(D0_Matrix_->getDomainMap(),D0_Matrix_->getRangeMap());
237  Teuchos::RCP<XMat> Ptent = MatrixFactory::Build(D0_Matrix_Abs->getRowMap(),0);
238  Xpetra::MatrixMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Multiply(*D0_Matrix_Abs,false,*P,false,*Ptent,true,true);
239 
240  // put in entries to P11
241  size_t dim = Nullspace_->getNumVectors();
242  size_t numLocalRows = SM_Matrix_->getNodeNumRows();
243  Teuchos::RCP<XMap> BlockColMap
245  P11_ = Teuchos::rcp(new XCrsWrap(Ptent->getRowMap(),BlockColMap,0,Xpetra::StaticProfile));
246 
247  std::vector< Teuchos::ArrayRCP<const Scalar> > nullspace(dim);
248  for(size_t i=0; i<dim; i++) {
249  Teuchos::ArrayRCP<const Scalar> datavec = Nullspace_->getData(i);
250  nullspace[i]=datavec;
251  }
252 
253  size_t nnz=0;
254  std::vector<size_t> rowPtrs;
255  std::vector<LocalOrdinal> blockCols;
256  std::vector<Scalar> blockVals;
257  for(size_t i=0; i<numLocalRows; i++) {
258  rowPtrs.push_back(nnz);
259  Teuchos::ArrayView<const LocalOrdinal> localCols;
260  Teuchos::ArrayView<const Scalar> localVals;
261  Ptent->getLocalRowView(i,localCols,localVals);
262  size_t numCols = localCols.size();
263  for(size_t j=0; j<numCols; j++) {
264  for(size_t k=0; k<dim; k++) {
265  blockCols.push_back(localCols[j]*dim+k);
266  blockVals.push_back(localVals[j]*nullspace[k][i]);
267  nnz++;
268  }
269  }
270  }
271  rowPtrs.push_back(nnz);
272 
273  ArrayRCP<size_t> rcpRowPtr;
274  ArrayRCP<LocalOrdinal> rcpColumns;
275  ArrayRCP<Scalar> rcpValues;
276 
277  RCP<XCRS> TP11 = rcp_dynamic_cast<XCrsWrap>(P11_)->getCrsMatrix();
278  TP11->allocateAllValues(nnz, rcpRowPtr, rcpColumns, rcpValues);
279 
280  ArrayView<size_t> rows = rcpRowPtr();
281  ArrayView<LocalOrdinal> columns = rcpColumns();
282  ArrayView<Scalar> values = rcpValues();
283 
284  for (size_t ii = 0; ii < rowPtrs.size(); ii++) rows[ii] = rowPtrs[ii];
285  for (size_t ii = 0; ii < blockCols.size(); ii++) columns[ii] = blockCols[ii];
286  for (size_t ii = 0; ii < blockVals.size(); ii++) values[ii] = blockVals[ii];
287  TP11->setAllValues(rcpRowPtr, rcpColumns, rcpValues);
288  Teuchos::RCP<XMap> blockCoarseMap
290  TP11->expertStaticFillComplete(blockCoarseMap,SM_Matrix_->getDomainMap());
291 
292 }
293 
294 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
296 
297  // coarse matrix for P11* (M1 + D1* M2 D1) P11
298  Teuchos::RCP<XMat> C = MatrixFactory::Build(SM_Matrix_->getRowMap(),0);
299  Teuchos::RCP<XMat> Matrix1 = MatrixFactory::Build(P11_->getDomainMap(),0);
300 
301  // construct (M1 + D1* M2 D1) P11
302  Xpetra::MatrixMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Multiply(*SM_Matrix_,false,*P11_,false,*C,true,true);
303 
304  // construct P11* (M1 + D1* M2 D1) P11
305  Xpetra::MatrixMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Multiply(*P11_,true,*C,false,*Matrix1,true,true);
306 
307  if(disable_addon_==true) {
308  // if add-on is not chosen
309  A11_=Matrix1;
310  }
311  else {
312  // catch a failure
313  TEUCHOS_TEST_FOR_EXCEPTION(M0inv_Matrix_==Teuchos::null,std::invalid_argument,
314  "MueLu::RefMaxwell::formCoarseMatrix(): Inverse of "
315  "lumped mass matrix required for add-on (i.e. M0inv_Matrix is null)");
316 
317  // coarse matrix for add-on, i.e P11* (M1 D0 M0inv D0* M1) P11
318  Teuchos::RCP<XMat> Zaux = MatrixFactory::Build(M1_Matrix_->getRowMap(),0);
319  Teuchos::RCP<XMat> Z = MatrixFactory::Build(D0_Matrix_->getDomainMap(),0);
320  Teuchos::RCP<XMat> C2 = MatrixFactory::Build(M0inv_Matrix_->getRowMap(),0);
321  // construct M1 P11
322  Xpetra::MatrixMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Multiply(*M1_Matrix_,false,*P11_,false,*Zaux,true,true);
323  // construct Z = D0* M1 P11
324  Xpetra::MatrixMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Multiply(*D0_Matrix_,true,*Zaux,false,*Z,true,true);
325  // construct M0inv Z
326  Xpetra::MatrixMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Multiply(*M0inv_Matrix_,false,*Z,false,*C2,true,true);
327  // construct Z* M0inv Z
328  Teuchos::RCP<XMat> Matrix2 = MatrixFactory::Build(Z->getDomainMap(),0);
330  // add matrices together
331  RCP<Teuchos::FancyOStream> out = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
332  Xpetra::MatrixMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>::TwoMatrixAdd(*Matrix1,false,(Scalar)1.0,*Matrix2,false,(Scalar)1.0,A11_,*out);
333  A11_->fillComplete();
334  }
335 
336  // set fixed block size for vector nodal matrix
337  size_t dim = Nullspace_->getNumVectors();
338  A11_->SetFixedBlockSize(dim);
339 
340 }
341 
342 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
344 
345  // convert Tpetra matrices to Xpetra
346  Teuchos::RCP<XCRS> SM_tmp = Teuchos::rcp( new XTCRS(SM_Matrix_new) );
347  SM_Matrix_ = Teuchos::rcp( new XCrsWrap(SM_tmp) );
348 
349 }
350 
351 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
353 
354  // compute residuals
355  RCP<XMV> residual = Utilities::Residual(*SM_Matrix_, X, RHS);
356  RCP<XMV> P11res = MultiVectorFactory::Build(P11_->getDomainMap(),X.getNumVectors());
357  RCP<XMV> P11x = MultiVectorFactory::Build(P11_->getDomainMap(),X.getNumVectors());
358  RCP<XMV> D0res = MultiVectorFactory::Build(D0_Matrix_->getDomainMap(),X.getNumVectors());
359  RCP<XMV> D0x = MultiVectorFactory::Build(D0_Matrix_->getDomainMap(),X.getNumVectors());
360  P11_->apply(*residual,*P11res,Teuchos::TRANS);
361  D0_Matrix_->apply(*residual,*D0res,Teuchos::TRANS);
362 
363  // block diagonal preconditioner on 2x2 (V-cycle for diagonal blocks)
364  Hierarchy11_->Iterate(*P11res, *P11x, 1, true);
365  Hierarchy22_->Iterate(*D0res, *D0x, 1, true);
366 
367  // update current solution
368  P11_->apply(*P11x,*residual,Teuchos::NO_TRANS);
369  D0_Matrix_->apply(*D0x,*residual,Teuchos::NO_TRANS,(Scalar)1.0,(Scalar)1.0);
370  X.update((Scalar) 1.0, *residual, (Scalar) 1.0);
371 
372 }
373 
374 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
376 
377  RCP<XMV> P11res = MultiVectorFactory::Build(P11_->getDomainMap(),X.getNumVectors());
378  RCP<XMV> P11x = MultiVectorFactory::Build(P11_->getDomainMap(),X.getNumVectors());
379  RCP<XMV> D0res = MultiVectorFactory::Build(D0_Matrix_->getDomainMap(),X.getNumVectors());
380  RCP<XMV> D0x = MultiVectorFactory::Build(D0_Matrix_->getDomainMap(),X.getNumVectors());
381 
382  // precondition (1,1)-block
383  RCP<XMV> residual = Utilities::Residual(*SM_Matrix_, X, RHS);
384  P11_->apply(*residual,*P11res,Teuchos::TRANS);
385  Hierarchy11_->Iterate(*P11res, *P11x, 1, true);
386  P11_->apply(*P11x,*residual,Teuchos::NO_TRANS);
387  X.update((Scalar) 1.0, *residual, (Scalar) 1.0);
388 
389  // precondition (2,2)-block
390  residual = Utilities::Residual(*SM_Matrix_, X, RHS);
391  D0_Matrix_->apply(*residual,*D0res,Teuchos::TRANS);
392  Hierarchy22_->Iterate(*D0res, *D0x, 1, true);
393  D0_Matrix_->apply(*D0x,*residual,Teuchos::NO_TRANS);
394  X.update((Scalar) 1.0, *residual, (Scalar) 1.0);
395 
396  // precondition (1,1)-block
397  residual = Utilities::Residual(*SM_Matrix_, X, RHS);
398  P11_->apply(*residual,*P11res,Teuchos::TRANS);
399  Hierarchy11_->Iterate(*P11res, *P11x, 1, true);
400  P11_->apply(*P11x,*residual,Teuchos::NO_TRANS);
401  X.update((Scalar) 1.0, *residual, (Scalar) 1.0);
402 
403 }
404 
405 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
407 
408  RCP<XMV> P11res = MultiVectorFactory::Build(P11_->getDomainMap(),X.getNumVectors());
409  RCP<XMV> P11x = MultiVectorFactory::Build(P11_->getDomainMap(),X.getNumVectors());
410  RCP<XMV> D0res = MultiVectorFactory::Build(D0_Matrix_->getDomainMap(),X.getNumVectors());
411  RCP<XMV> D0x = MultiVectorFactory::Build(D0_Matrix_->getDomainMap(),X.getNumVectors());
412 
413  // precondition (2,2)-block
414  RCP<XMV> residual = Utilities::Residual(*SM_Matrix_, X, RHS);
415  D0_Matrix_->apply(*residual,*D0res,Teuchos::TRANS);
416  Hierarchy22_->Iterate(*D0res, *D0x, 1, true);
417  D0_Matrix_->apply(*D0x,*residual,Teuchos::NO_TRANS);
418  X.update((Scalar) 1.0, *residual, (Scalar) 1.0);
419 
420  // precondition (1,1)-block
421  residual = Utilities::Residual(*SM_Matrix_, X, RHS);
422  P11_->apply(*residual,*P11res,Teuchos::TRANS);
423  Hierarchy11_->Iterate(*P11res, *P11x, 1, true);
424  P11_->apply(*P11x,*residual,Teuchos::NO_TRANS);
425  X.update((Scalar) 1.0, *residual, (Scalar) 1.0);
426 
427  // precondition (2,2)-block
428  residual = Utilities::Residual(*SM_Matrix_, X, RHS);
429  D0_Matrix_->apply(*residual,*D0res,Teuchos::TRANS);
430  Hierarchy22_->Iterate(*D0res, *D0x, 1, true);
431  D0_Matrix_->apply(*D0x,*residual,Teuchos::NO_TRANS);
432  X.update((Scalar) 1.0, *residual, (Scalar) 1.0);
433 
434 }
435 
436 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
437 void RefMaxwell<Scalar,LocalOrdinal,GlobalOrdinal,Node>::apply(const Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>& X,
438  Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>& Y,
439  Teuchos::ETransp mode, Scalar alpha, Scalar beta) const {
440 
441  try {
442 
443  TMV& temp_x = const_cast<TMV &>(X);
444  const XTMV tX(rcpFromRef(temp_x));
445  XTMV tY(rcpFromRef(Y));
446  tY.putScalar(Teuchos::ScalarTraits<Scalar>::zero());
447 
448  // apply pre-smoothing
449  HierarchySmoother_->Iterate(tX,tY,1);
450 
451  // do solve for the 2x2 block system
452  if(mode_=="additive")
453  applyInverseAdditive(tX,tY);
454  else if(mode_=="121")
455  applyInverse121(tX,tY);
456  else if(mode_=="212")
457  applyInverse212(tX,tY);
458  else
459  applyInverseAdditive(tX,tY);
460 
461  // apply post-smoothing
462  HierarchySmoother_->Iterate(tX,tY,1);
463 
464  } catch (std::exception& e) {
465 
466  //FIXME add message and rethrow
467  std::cerr << "Caught an exception in MueLu::RefMaxwell::ApplyInverse():" << std::endl
468  << e.what() << std::endl;
469 
470  }
471 }
472 
473 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
475  return false;
476 }
477 
478 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
480 initialize(const Teuchos::RCP<TCRS> & D0_Matrix,
481  const Teuchos::RCP<TCRS> & M0inv_Matrix,
482  const Teuchos::RCP<TCRS> & M1_Matrix,
483  const Teuchos::RCP<TMV> & Nullspace,
484  const Teuchos::RCP<TMV> & Coords,
485  Teuchos::ParameterList& List)
486 {
487  // some pre-conditions
488  TEUCHOS_ASSERT(D0_Matrix!=Teuchos::null);
489  TEUCHOS_ASSERT(M1_Matrix!=Teuchos::null);
490 
491  Hierarchy11_ = Teuchos::null;
492  Hierarchy22_ = Teuchos::null;
493  HierarchySmoother_ = Teuchos::null;
494  parameterList_ = List;
495  disable_addon_ = false;
496  mode_ = "additive";
497 
498  // set parameters
499  setParameters(List);
500 
501  // convert Tpetra matrices to Xpetra
502  Teuchos::RCP<XCRS> D0_tmp = Teuchos::rcp( new XTCRS(D0_Matrix) );
503  D0_Matrix_ = Teuchos::rcp( new XCrsWrap(D0_tmp) );
504 
505  if(M0inv_Matrix != Teuchos::null) {
506  Teuchos::RCP<XCRS> M0inv_tmp = Teuchos::rcp( new XTCRS(M0inv_Matrix) );
507  M0inv_Matrix_ = Teuchos::rcp( new XCrsWrap(M0inv_tmp) );
508  }
509 
510  Teuchos::RCP<XCRS> M1_tmp = Teuchos::rcp( new XTCRS(M1_Matrix) );
511  M1_Matrix_ = Teuchos::rcp( new XCrsWrap(M1_tmp) );
512 
513  // convert Tpetra MultiVector to Xpetra
514  if(Coords != Teuchos::null)
515  Coords_ = Xpetra::toXpetra(Coords);
516  if(Nullspace != Teuchos::null)
517  Nullspace_ = Xpetra::toXpetra(Nullspace);
518 }
519 
520 } // namespace
521 #endif //ifdef HAVE_MUELU_TPETRA
522 
523 #define MUELU_REFMAXWELL_SHORT
524 #endif //ifdef MUELU_REFMAXWELL_DEF_HPP
size_t getNumVectors() const
This class specifies the default factory that should generate some data on a Level if the data does n...
static Teuchos::RCP< Map< LocalOrdinal, GlobalOrdinal, Node > > Build(UnderlyingLib lib, global_size_t numGlobalElements, GlobalOrdinal indexBase, const Teuchos::RCP< const Teuchos::Comm< int > > &comm, LocalGlobal lg=Xpetra::GloballyDistributed, const Teuchos::RCP< Node > &node=defaultArgNode())
const RCP< const Tpetra::Map< LocalOrdinal, GlobalOrdinal, Node > > toTpetraNonZero(const RCP< const Map< LocalOrdinal, GlobalOrdinal, Node > > &map)
void applyInverse121(const XTMV &RHS, XTMV &X) const
apply 1-2-1 algorithm for 2x2 solve
Namespace for MueLu classes and methods.
void setParameters(Teuchos::ParameterList &list)
Set parameters.
void buildProlongator()
Setup the prolongator for the (1,1)-block.
Factory for building tentative prolongator.
Teuchos::RCP< const Tpetra::Map< LocalOrdinal, GlobalOrdinal, Node > > getRangeMap() const
Returns the Tpetra::Map object associated with the range of this operator.
void applyInverse212(const XTMV &RHS, XTMV &X) const
apply 2-1-2 algorithm for 2x2 solve
void putScalar(const Scalar &value)
static RCP< Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Residual(const Xpetra::Operator< Scalar, LocalOrdinal, GlobalOrdinal, Node > &Op, const Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > &X, const Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > &RHS)
void compute()
Setup the preconditioner.
void update(const Scalar &alpha, const MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A, const Scalar &beta)
Always keep data, even accross run. This flag is set by Level::Keep(). This flag is propagated to coa...
void applyInverseAdditive(const XTMV &RHS, XTMV &X) const
apply additive algorithm for 2x2 solve
static void TwoMatrixAdd(const Matrix &A, bool transposeA, SC alpha, Matrix &B, SC beta)
Class that accepts ML-style parameters and builds a MueLu preconditioner. This interpreter uses the s...
bool hasTransposeApply() const
Indicates whether this operator supports applying the adjoint operator.
RCP< const CrsGraph< int, GlobalOrdinal, Node > > toXpetra(const Epetra_CrsGraph &g)
Tpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > TMV
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)
void initialize(const Teuchos::RCP< TCRS > &D0_Matrix, const Teuchos::RCP< TCRS > &M0inv_Matrix, const Teuchos::RCP< TCRS > &M1_Matrix, const Teuchos::RCP< TMV > &Nullspace, const Teuchos::RCP< TMV > &Coords, Teuchos::ParameterList &List)
void apply(const Tpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > &X, Tpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > &Y, Teuchos::ETransp mode=Teuchos::NO_TRANS, Scalar alpha=Teuchos::ScalarTraits< Scalar >::one(), Scalar beta=Teuchos::ScalarTraits< Scalar >::one()) const
void formCoarseMatrix()
Compute P11^{T}*A*P11 efficiently.
void resetMatrix(Teuchos::RCP< TCRS > SM_Matrix_new)
Reset system matrix.
Factory for building Smoothed Aggregation prolongators.Input/output of SaPFactory
Factory for building uncoupled aggregates.
Teuchos::RCP< const Tpetra::Map< LocalOrdinal, GlobalOrdinal, Node > > getDomainMap() const
Returns the Tpetra::Map object associated with the domain of this operator.
Provides methods to build a multigrid hierarchy and apply multigrid cycles.