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_Map.hpp"
52 #include "Xpetra_MatrixMatrix.hpp"
55 #include "Xpetra_MatrixUtils.hpp"
56 
58 
59 #include "MueLu_AmalgamationFactory.hpp"
60 #include "MueLu_RAPFactory.hpp"
61 #include "MueLu_ThresholdAFilterFactory.hpp"
62 #include "MueLu_TransPFactory.hpp"
63 #include "MueLu_SmootherFactory.hpp"
64 
65 #include "MueLu_CoalesceDropFactory.hpp"
66 #include "MueLu_CoarseMapFactory.hpp"
67 #include "MueLu_CoordinatesTransferFactory.hpp"
68 #include "MueLu_UncoupledAggregationFactory.hpp"
69 #include "MueLu_TentativePFactory.hpp"
70 #include "MueLu_Utilities.hpp"
71 
72 #ifdef HAVE_MUELU_KOKKOS_REFACTOR
73 #include "MueLu_CoalesceDropFactory_kokkos.hpp"
74 #include "MueLu_CoarseMapFactory_kokkos.hpp"
75 #include "MueLu_CoordinatesTransferFactory_kokkos.hpp"
76 #include "MueLu_UncoupledAggregationFactory_kokkos.hpp"
77 #include "MueLu_TentativePFactory_kokkos.hpp"
78 #include "MueLu_Utilities_kokkos.hpp"
79 #include <Kokkos_Core.hpp>
80 #include <KokkosSparse_CrsMatrix.hpp>
81 #endif
82 
83 
84 #include "MueLu_Monitor.hpp"
85 #include "MueLu_MLParameterListInterpreter.hpp"
86 #include "MueLu_ParameterListInterpreter.hpp"
87 #include "MueLu_VerbosityLevel.hpp"
88 
90 
91 #ifdef HAVE_MUELU_CUDA
92 #include "cuda_profiler_api.h"
93 #endif
94 
95 
96 namespace MueLu {
97 
98  template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
100  return SM_Matrix_->getDomainMap();
101  }
102 
103 
104  template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
106  return SM_Matrix_->getRangeMap();
107  }
108 
109 
110  template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
112 
113  parameterList_ = list;
114  disable_addon_ = list.get("refmaxwell: disable addon",true);
115  mode_ = list.get("refmaxwell: mode","additive");
116  use_as_preconditioner_ = list.get<bool>("refmaxwell: use as preconditioner");
117  dump_matrices_ = list.get("refmaxwell: dump matrices",false);
118 
119  if(list.isSublist("refmaxwell: 11list"))
120  precList11_ = list.sublist("refmaxwell: 11list");
121 
122  if(list.isSublist("refmaxwell: 22list"))
123  precList22_ = list.sublist("refmaxwell: 22list");
124 
125  if(list.isSublist("smoother: params")) {
126  smootherList_ = list.sublist("smoother: params");
127  }
128 
129 #if !defined(HAVE_MUELU_KOKKOS_REFACTOR)
130  useKokkos_ = false;
131 #elif defined(HAVE_MUELU_KOKKOS_REFACTOR_USE_BY_DEFAULT)
132  useKokkos_ = list.get("use kokkos refactor",true);
133 #else
134  useKokkos_ = list.get("use kokkos refactor",false);
135 #endif
136 
137  }
138 
139 
140  template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
142 
143 
144 #ifdef HAVE_MUELU_CUDA
145  if (parameterList_.get<bool>("refmaxwell: cuda profile setup", false)) cudaProfilerStart();
146 #endif
147 
148  Teuchos::TimeMonitor tmCompute(*Teuchos::TimeMonitor::getNewTimer("MueLu RefMaxwell: compute"));
149 
150  bool defaultFilter = false;
151 
152  // Remove zero entries from D0 if necessary.
153  // In the construction of the prolongator we use the graph of the
154  // matrix, so zero entries mess it up.
155  if (parameterList_.get<bool>("refmaxwell: filter D0", true) && D0_Matrix_->getNodeMaxNumRowEntries()>2) {
156  Level fineLevel;
157  fineLevel.SetFactoryManager(Teuchos::null);
158  fineLevel.SetLevelID(0);
159  fineLevel.Set("A",D0_Matrix_);
160  fineLevel.setlib(D0_Matrix_->getDomainMap()->lib());
161  // We expect D0 to have entries +-1, so any threshold value will do.
162  RCP<ThresholdAFilterFactory> ThreshFact = rcp(new ThresholdAFilterFactory("A",1.0e-8,/*keepDiagonal=*/false,/*expectedNNZperRow=*/2));
163  fineLevel.Request("A",ThreshFact.get());
164  ThreshFact->Build(fineLevel);
165  D0_Matrix_ = fineLevel.Get< RCP<Matrix> >("A",ThreshFact.get());
166 
167  // If D0 has too many zeros, maybe SM and M1 do as well.
168  defaultFilter = true;
169  }
170 
171  if (parameterList_.get<bool>("refmaxwell: filter SM", defaultFilter)) {
172  RCP<Vector> diag = VectorFactory::Build(SM_Matrix_->getRowMap());
173  // find a reasonable absolute value threshold
174  SM_Matrix_->getLocalDiagCopy(*diag);
175  magnitudeType threshold = 1.0e-8 * diag->normInf();
176 
177  Level fineLevel;
178  fineLevel.SetFactoryManager(Teuchos::null);
179  fineLevel.SetLevelID(0);
180  fineLevel.Set("A",SM_Matrix_);
181  fineLevel.setlib(SM_Matrix_->getDomainMap()->lib());
182  RCP<ThresholdAFilterFactory> ThreshFact = rcp(new ThresholdAFilterFactory("A",threshold,/*keepDiagonal=*/true));
183  fineLevel.Request("A",ThreshFact.get());
184  ThreshFact->Build(fineLevel);
185  SM_Matrix_ = fineLevel.Get< RCP<Matrix> >("A",ThreshFact.get());
186  }
187 
188  if (parameterList_.get<bool>("refmaxwell: filter M1", defaultFilter)) {
189  RCP<Vector> diag = VectorFactory::Build(M1_Matrix_->getRowMap());
190  // find a reasonable absolute value threshold
191  M1_Matrix_->getLocalDiagCopy(*diag);
192  magnitudeType threshold = 1.0e-8 * diag->normInf();
193 
194  Level fineLevel;
195  fineLevel.SetFactoryManager(Teuchos::null);
196  fineLevel.SetLevelID(0);
197  fineLevel.Set("A",M1_Matrix_);
198  fineLevel.setlib(M1_Matrix_->getDomainMap()->lib());
199  RCP<ThresholdAFilterFactory> ThreshFact = rcp(new ThresholdAFilterFactory("A",threshold,/*keepDiagonal=*/true));
200  fineLevel.Request("A",ThreshFact.get());
201  ThreshFact->Build(fineLevel);
202  M1_Matrix_ = fineLevel.Get< RCP<Matrix> >("A",ThreshFact.get());
203  }
204 
205  // clean rows associated with boundary conditions
206  // Find rows with only 1 or 2 nonzero entries, record them in BCrows_.
207  // BCrows_[i] is true, iff i is a boundary row
208  // BCcols_[i] is true, iff i is a boundary column
209 #ifdef HAVE_MUELU_KOKKOS_REFACTOR
210  if (useKokkos_) {
211  BCrowsKokkos_ = Utilities_kokkos::DetectDirichletRows(*SM_Matrix_,Teuchos::ScalarTraits<magnitudeType>::eps(),/*count_twos_as_dirichlet=*/true);
212  BCcolsKokkos_ = Utilities_kokkos::DetectDirichletCols(*D0_Matrix_,BCrowsKokkos_);
213  } else {
214  BCrows_ = Utilities::DetectDirichletRows(*SM_Matrix_,Teuchos::ScalarTraits<magnitudeType>::eps(),/*count_twos_as_dirichlet=*/true);
215  BCcols_ = Utilities::DetectDirichletCols(*D0_Matrix_,BCrows_);
216  }
217 #else
218  BCrows_ = Utilities::DetectDirichletRows(*SM_Matrix_,Teuchos::ScalarTraits<magnitudeType>::eps(),/*count_twos_as_dirichlet=*/true);
219  BCcols_ = Utilities::DetectDirichletCols(*D0_Matrix_,BCrows_);
220 #endif
221 
222  // build nullspace if necessary
223  if(Nullspace_ != Teuchos::null) {
224  // no need to do anything - nullspace is built
225  }
226  else if(Nullspace_ == Teuchos::null && Coords_ != Teuchos::null) {
227  // normalize coordinates
228  typedef typename RealValuedMultiVector::scalar_type realScalarType;
229  typedef typename Teuchos::ScalarTraits<realScalarType>::magnitudeType realMagnitudeType;
230  Teuchos::Array<realMagnitudeType> norms(Coords_->getNumVectors());
231  Coords_->norm2(norms);
232  for (size_t i=0;i<Coords_->getNumVectors();i++)
233  norms[i] = ((realMagnitudeType)1.0)/norms[i];
234  Coords_->scale(norms());
235  Nullspace_ = MultiVectorFactory::Build(SM_Matrix_->getRowMap(),Coords_->getNumVectors());
236 
237  // Cast coordinates to Scalar so they can be multiplied against D0
238 #ifdef HAVE_MUELU_KOKKOS_REFACTOR
239  RCP<MultiVector> CoordsSC;
240  if (useKokkos_)
241  CoordsSC = Utilities_kokkos::RealValuedToScalarMultiVector(Coords_);
242  else
243  CoordsSC = Utilities::RealValuedToScalarMultiVector(Coords_);
244 #else
246 #endif
247  D0_Matrix_->apply(*CoordsSC,*Nullspace_);
248  }
249  else {
250  GetOStream(Errors) << "MueLu::RefMaxwell::compute(): either the nullspace or the nodal coordinates must be provided." << std::endl;
251  }
252 
253  // Nuke the BC edges in nullspace
254 #ifdef HAVE_MUELU_KOKKOS_REFACTOR
255  if (useKokkos_)
256  Utilities_kokkos::ZeroDirichletRows(Nullspace_,BCrowsKokkos_);
257  else
258  Utilities::ZeroDirichletRows(Nullspace_,BCrows_);
259 #else
260  Utilities::ZeroDirichletRows(Nullspace_,BCrows_);
261 #endif
262 
263  if (dump_matrices_)
264  Xpetra::IO<SC, LO, GlobalOrdinal, Node>::Write(std::string("D0_clean.mat"), *D0_Matrix_);
265 
266  // build special prolongator for (1,1)-block
267  if(P11_==Teuchos::null) {
268  // Form A_nodal = D0* M1 D0 (aka TMT_agg)
269  Level fineLevel, coarseLevel;
270  fineLevel.SetFactoryManager(Teuchos::null);
271  coarseLevel.SetFactoryManager(Teuchos::null);
272  coarseLevel.SetPreviousLevel(rcpFromRef(fineLevel));
273  fineLevel.SetLevelID(0);
274  coarseLevel.SetLevelID(1);
275  fineLevel.Set("A",M1_Matrix_);
276  coarseLevel.Set("P",D0_Matrix_);
277  coarseLevel.setlib(M1_Matrix_->getDomainMap()->lib());
278  fineLevel.setlib(M1_Matrix_->getDomainMap()->lib());
279  coarseLevel.setObjectLabel("RefMaxwell (1,1) A_nodal");
280  fineLevel.setObjectLabel("RefMaxwell (1,1) A_nodal");
281 
282  RCP<RAPFactory> rapFact = rcp(new RAPFactory());
283  Teuchos::ParameterList rapList = *(rapFact->GetValidParameterList());
284  rapList.set("transpose: use implicit", parameterList_.get<bool>("transpose: use implicit", false));
285  rapList.set("rap: fix zero diagonals", parameterList_.get<bool>("rap: fix zero diagonals", true));
286  rapList.set("rap: triple product", parameterList_.get<bool>("rap: triple product", false));
287  rapFact->SetParameterList(rapList);
288 
289  RCP<TransPFactory> transPFactory;
290  if (!parameterList_.get<bool>("transpose: use implicit", false)) {
291  transPFactory = rcp(new TransPFactory());
292  rapFact->SetFactory("R", transPFactory);
293  }
294 
295  coarseLevel.Request("A", rapFact.get());
296 
297  A_nodal_Matrix_ = coarseLevel.Get< RCP<Matrix> >("A", rapFact.get());
298 
299  // build special prolongator
300  GetOStream(Runtime0) << "RefMaxwell::compute(): building special prolongator" << std::endl;
301  buildProlongator();
302 
303 #ifdef HAVE_MUELU_KOKKOS_REFACTOR
304  if (useKokkos_)
305  R11_ = Utilities_kokkos::Transpose(*P11_);
306  else
307  R11_ = Utilities::Transpose(*P11_);
308 #else
309  R11_ = Utilities::Transpose(*P11_);
310 #endif
311  }
312 
313  {
314  // build coarse grid operator for (1,1)-block
315  formCoarseMatrix();
316 
317 #ifdef HAVE_MUELU_KOKKOS_REFACTOR
318  // This should be taken out again as soon as
319  // CoalesceDropFactory_kokkos supports BlockSize > 1 and
320  // drop tol != 0.0
321  if (useKokkos_ && precList11_.isParameter("aggregation: drop tol") && precList11_.get<double>("aggregation: drop tol") != 0.0) {
322  GetOStream(Warnings0) << "RefMaxwell::compute(): Setting \"aggregation: drop tol\". to 0.0, since CoalesceDropFactory_kokkos does not "
323  << "support BlockSize > 1 and drop tol != 0.0" << std::endl;
324  precList11_.set("aggregation: drop tol", 0.0);
325  }
326 #endif
327  HierarchyH_ = MueLu::CreateXpetraPreconditioner(AH_, precList11_, CoordsH_);
328  }
329 
330  {
331  GetOStream(Runtime0) << "RefMaxwell::compute(): nuking BC edges of D0" << std::endl;
332 
333  D0_Matrix_->resumeFill();
334  // Scalar replaceWith = Teuchos::ScalarTraits<SC>::eps();
335  Scalar replaceWith = Teuchos::ScalarTraits<SC>::zero();
336 #ifdef HAVE_MUELU_KOKKOS_REFACTOR
337  if (useKokkos_) {
338  Utilities_kokkos::ZeroDirichletRows(D0_Matrix_,BCrowsKokkos_,replaceWith);
339  Utilities_kokkos::ZeroDirichletCols(D0_Matrix_,BCcolsKokkos_,replaceWith);
340  } else {
341  Utilities::ZeroDirichletRows(D0_Matrix_,BCrows_,replaceWith);
342  Utilities::ZeroDirichletCols(D0_Matrix_,BCcols_,replaceWith);
343  }
344 #else
345  Utilities::ZeroDirichletRows(D0_Matrix_,BCrows_,replaceWith);
346  Utilities::ZeroDirichletCols(D0_Matrix_,BCcols_,replaceWith);
347 #endif
348  D0_Matrix_->fillComplete(D0_Matrix_->getDomainMap(),D0_Matrix_->getRangeMap());
349  }
350 
351  {
352  GetOStream(Runtime0) << "RefMaxwell::compute(): building MG for (2,2)-block" << std::endl;
353 
354  { // build fine grid operator for (2,2)-block, D0* SM D0 (aka TMT)
355  Teuchos::TimeMonitor tm(*Teuchos::TimeMonitor::getNewTimer("MueLu RefMaxwell: Build A22"));
356 
357  Level fineLevel, coarseLevel;
358  fineLevel.SetFactoryManager(Teuchos::null);
359  coarseLevel.SetFactoryManager(Teuchos::null);
360  coarseLevel.SetPreviousLevel(rcpFromRef(fineLevel));
361  fineLevel.SetLevelID(0);
362  coarseLevel.SetLevelID(1);
363  fineLevel.Set("A",SM_Matrix_);
364  coarseLevel.Set("P",D0_Matrix_);
365  coarseLevel.setlib(SM_Matrix_->getDomainMap()->lib());
366  fineLevel.setlib(SM_Matrix_->getDomainMap()->lib());
367  coarseLevel.setObjectLabel("RefMaxwell (2,2)");
368  fineLevel.setObjectLabel("RefMaxwell (2,2)");
369 
370  RCP<RAPFactory> rapFact = rcp(new RAPFactory());
371  Teuchos::ParameterList rapList = *(rapFact->GetValidParameterList());
372  rapList.set("transpose: use implicit", false);
373  rapList.set("rap: fix zero diagonals", parameterList_.get<bool>("rap: fix zero diagonals", true));
374  rapList.set("rap: triple product", parameterList_.get<bool>("rap: triple product", false));
375  rapFact->SetParameterList(rapList);
376 
377  RCP<TransPFactory> transPFactory;
378  transPFactory = rcp(new TransPFactory());
379  rapFact->SetFactory("R", transPFactory);
380 
381  coarseLevel.Request("R", transPFactory.get());
382  coarseLevel.Request("A", rapFact.get());
383  A22_ = coarseLevel.Get< RCP<Matrix> >("A", rapFact.get());
384  A22_->setObjectLabel("RefMaxwell (2,2)");
385  D0_T_Matrix_ = coarseLevel.Get< RCP<Matrix> >("R", transPFactory.get());
386  }
387 
388  Hierarchy22_ = MueLu::CreateXpetraPreconditioner(A22_, precList22_, Coords_);
389  }
390 
391  {
392  if (parameterList_.isType<std::string>("smoother: type") &&
393  parameterList_.get<std::string>("smoother: type") == "hiptmair" &&
394  SM_Matrix_->getDomainMap()->lib() == Xpetra::UseTpetra &&
395  A22_->getDomainMap()->lib() == Xpetra::UseTpetra &&
396  D0_Matrix_->getDomainMap()->lib() == Xpetra::UseTpetra) {
397 #if defined(HAVE_MUELU_IFPACK2) && (!defined(HAVE_MUELU_EPETRA) || (defined(HAVE_MUELU_INST_DOUBLE_INT_INT)))
398  Teuchos::ParameterList hiptmairPreList, hiptmairPostList, smootherPreList, smootherPostList;
399 
400  if (smootherList_.isSublist("smoother: pre params"))
401  smootherPreList = smootherList_.sublist("smoother: pre params");
402  else if (smootherList_.isSublist("smoother: params"))
403  smootherPreList = smootherList_.sublist("smoother: params");
404  hiptmairPreList.set("hiptmair: smoother type 1",
405  smootherPreList.get<std::string>("hiptmair: smoother type 1", "CHEBYSHEV"));
406  hiptmairPreList.set("hiptmair: smoother type 2",
407  smootherPreList.get<std::string>("hiptmair: smoother type 2", "CHEBYSHEV"));
408  if(smootherPreList.isSublist("hiptmair: smoother list 1"))
409  hiptmairPreList.set("hiptmair: smoother list 1", smootherPreList.sublist("hiptmair: smoother list 1"));
410  if(smootherPreList.isSublist("hiptmair: smoother list 2"))
411  hiptmairPreList.set("hiptmair: smoother list 2", smootherPreList.sublist("hiptmair: smoother list 2"));
412  hiptmairPreList.set("hiptmair: pre or post",
413  smootherPreList.get<std::string>("hiptmair: pre or post", "pre"));
414  hiptmairPreList.set("hiptmair: zero starting solution",
415  smootherPreList.get<bool>("hiptmair: zero starting solution", true));
416 
417  if (smootherList_.isSublist("smoother: post params"))
418  smootherPostList = smootherList_.sublist("smoother: post params");
419  else if (smootherList_.isSublist("smoother: params"))
420  smootherPostList = smootherList_.sublist("smoother: params");
421  hiptmairPostList.set("hiptmair: smoother type 1",
422  smootherPostList.get<std::string>("hiptmair: smoother type 1", "CHEBYSHEV"));
423  hiptmairPostList.set("hiptmair: smoother type 2",
424  smootherPostList.get<std::string>("hiptmair: smoother type 2", "CHEBYSHEV"));
425  if(smootherPostList.isSublist("hiptmair: smoother list 1"))
426  hiptmairPostList.set("hiptmair: smoother list 1", smootherPostList.sublist("hiptmair: smoother list 1"));
427  if(smootherPostList.isSublist("hiptmair: smoother list 2"))
428  hiptmairPostList.set("hiptmair: smoother list 2", smootherPostList.sublist("hiptmair: smoother list 2"));
429  hiptmairPostList.set("hiptmair: pre or post",
430  smootherPostList.get<std::string>("hiptmair: pre or post", "post"));
431  hiptmairPostList.set("hiptmair: zero starting solution",
432  smootherPostList.get<bool>("hiptmair: zero starting solution", false));
433 
434  typedef Tpetra::RowMatrix<SC, LO, GO, NO> TROW;
438 
439  hiptmairPreSmoother_ = Teuchos::rcp( new Ifpack2::Hiptmair<TROW>(EdgeMatrix,NodeMatrix,PMatrix) );
440  hiptmairPreSmoother_ -> setParameters(hiptmairPreList);
441  hiptmairPreSmoother_ -> initialize();
442  hiptmairPreSmoother_ -> compute();
443  hiptmairPostSmoother_ = Teuchos::rcp( new Ifpack2::Hiptmair<TROW>(EdgeMatrix,NodeMatrix,PMatrix) );
444  hiptmairPostSmoother_ -> setParameters(hiptmairPostList);
445  hiptmairPostSmoother_ -> initialize();
446  hiptmairPostSmoother_ -> compute();
447  useHiptmairSmoothing_ = true;
448 #else
449  throw(Xpetra::Exceptions::RuntimeError("MueLu must be compiled with Ifpack2 for Hiptmair smoothing."));
450 #endif // defined(HAVE_MUELU_IFPACK2) && (!defined(HAVE_MUELU_EPETRA) || defined(HAVE_MUELU_INST_DOUBLE_INT_INT))
451  } else {
452  if (parameterList_.isType<std::string>("smoother: pre type") && parameterList_.isType<std::string>("smoother: post type")) {
453  std::string preSmootherType = parameterList_.get<std::string>("smoother: pre type");
454  std::string postSmootherType = parameterList_.get<std::string>("smoother: post type");
455 
456  Teuchos::ParameterList preSmootherList, postSmootherList;
457  if (parameterList_.isSublist("smoother: pre params"))
458  preSmootherList = parameterList_.sublist("smoother: pre params");
459  if (parameterList_.isSublist("smoother: post params"))
460  postSmootherList = parameterList_.sublist("smoother: post params");
461 
462  Level level;
463  RCP<MueLu::FactoryManagerBase> factoryHandler = rcp(new FactoryManager());
464  level.SetFactoryManager(factoryHandler);
465  level.SetLevelID(0);
466  level.setObjectLabel("RefMaxwell (1,1)");
467  level.Set("A",SM_Matrix_);
468  level.setlib(SM_Matrix_->getDomainMap()->lib());
469 
470  RCP<SmootherPrototype> preSmootherPrototype = rcp(new TrilinosSmoother(preSmootherType, preSmootherList));
471  RCP<SmootherFactory> preSmootherFact = rcp(new SmootherFactory(preSmootherPrototype));
472 
473  RCP<SmootherPrototype> postSmootherPrototype = rcp(new TrilinosSmoother(postSmootherType, postSmootherList));
474  RCP<SmootherFactory> postSmootherFact = rcp(new SmootherFactory(postSmootherPrototype));
475 
476  level.Request("PreSmoother",preSmootherFact.get());
477  preSmootherFact->Build(level);
478  PreSmoother_ = level.Get<RCP<SmootherBase> >("PreSmoother",preSmootherFact.get());
479 
480  level.Request("PostSmoother",postSmootherFact.get());
481  postSmootherFact->Build(level);
482  PostSmoother_ = level.Get<RCP<SmootherBase> >("PostSmoother",postSmootherFact.get());
483  } else {
484  std::string smootherType = parameterList_.get<std::string>("smoother: type", "CHEBYSHEV");
485  Level level;
486  RCP<MueLu::FactoryManagerBase> factoryHandler = rcp(new FactoryManager());
487  level.SetFactoryManager(factoryHandler);
488  level.SetLevelID(0);
489  level.setObjectLabel("RefMaxwell (1,1)");
490  level.Set("A",SM_Matrix_);
491  level.setlib(SM_Matrix_->getDomainMap()->lib());
492  RCP<SmootherPrototype> smootherPrototype = rcp(new TrilinosSmoother(smootherType, smootherList_));
493  RCP<SmootherFactory> SmootherFact = rcp(new SmootherFactory(smootherPrototype));
494  level.Request("PreSmoother",SmootherFact.get());
495  SmootherFact->Build(level);
496  PreSmoother_ = level.Get<RCP<SmootherBase> >("PreSmoother",SmootherFact.get());
497  PostSmoother_ = PreSmoother_;
498  }
499  useHiptmairSmoothing_ = false;
500  }
501  }
502 
503  // Allocate temporary MultiVectors for solve
504  P11res_ = MultiVectorFactory::Build(P11_->getDomainMap(),1);
505  P11x_ = MultiVectorFactory::Build(P11_->getDomainMap(),1);
506  D0res_ = MultiVectorFactory::Build(D0_Matrix_->getDomainMap(),1);
507  D0x_ = MultiVectorFactory::Build(D0_Matrix_->getDomainMap(),1);
508  residual_ = MultiVectorFactory::Build(SM_Matrix_->getDomainMap(),1);
509 
510 #ifdef HAVE_MUELU_CUDA
511  if (parameterList_.get<bool>("refmaxwell: cuda profile setup", false)) cudaProfilerStop();
512 #endif
513 
514  describe(GetOStream(Runtime0));
515 
516  if (dump_matrices_) {
517  GetOStream(Runtime0) << "RefMaxwell::compute(): dumping data" << std::endl;
518  Xpetra::IO<SC, LO, GlobalOrdinal, Node>::Write(std::string("SM.mat"), *SM_Matrix_);
519  Xpetra::IO<SC, LO, GlobalOrdinal, Node>::Write(std::string("M1.mat"), *M1_Matrix_);
520  Xpetra::IO<SC, LO, GlobalOrdinal, Node>::Write(std::string("M0inv.mat"), *M0inv_Matrix_);
521 #ifndef HAVE_MUELU_KOKKOS_REFACTOR
522  std::ofstream outBCrows("BCrows.mat");
523  std::copy(BCrows_.begin(), BCrows_.end(), std::ostream_iterator<LO>(outBCrows, "\n"));
524  std::ofstream outBCcols("BCcols.mat");
525  std::copy(BCcols_.begin(), BCcols_.end(), std::ostream_iterator<LO>(outBCcols, "\n"));
526 #else
527  if (useKokkos_) {
528  std::ofstream outBCrows("BCrows.mat");
529  auto BCrows = Kokkos::create_mirror_view (BCrowsKokkos_);
530  Kokkos::deep_copy(BCrows , BCrowsKokkos_);
531  for (size_t i = 0; i < BCrows.size(); i++)
532  outBCrows << BCrows[i] << "\n";
533 
534  std::ofstream outBCcols("BCcols.mat");
535  auto BCcols = Kokkos::create_mirror_view (BCcolsKokkos_);
536  Kokkos::deep_copy(BCcols , BCcolsKokkos_);
537  for (size_t i = 0; i < BCcols.size(); i++)
538  outBCcols << BCcols[i] << "\n";
539  } else {
540  std::ofstream outBCrows("BCrows.mat");
541  std::copy(BCrows_.begin(), BCrows_.end(), std::ostream_iterator<LO>(outBCrows, "\n"));
542  std::ofstream outBCcols("BCcols.mat");
543  std::copy(BCcols_.begin(), BCcols_.end(), std::ostream_iterator<LO>(outBCcols, "\n"));
544  }
545 #endif
546  Xpetra::IO<SC, LO, GlobalOrdinal, Node>::Write(std::string("nullspace.mat"), *Nullspace_);
547  if (Coords_ != Teuchos::null)
548  Xpetra::IO<double, LO, GlobalOrdinal, Node>::Write(std::string("coords.mat"), *Coords_);
549  Xpetra::IO<SC, LO, GlobalOrdinal, Node>::Write(std::string("D0_nuked.mat"), *D0_Matrix_);
550  Xpetra::IO<SC, LO, GlobalOrdinal, Node>::Write(std::string("A_nodal.mat"), *A_nodal_Matrix_);
551  Xpetra::IO<SC, LO, GO, NO>::Write(std::string("P11.mat"), *P11_);
552  Xpetra::IO<SC, LO, GlobalOrdinal, Node>::Write(std::string("AH.mat"), *AH_);
553  Xpetra::IO<SC, LO, GlobalOrdinal, Node>::Write(std::string("A22.mat"), *A22_);
554  }
555  }
556 
557 
558  template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
560  // The P11 matrix maps node based aggregrates { A_j } to edges { e_i }.
561  //
562  // The old implementation used
563  // P11(i, j*dim+k) = sum_{nodes n_l in e_i intersected with A_j} 0.5 * phi_k(e_i) * P(n_l, A_j)
564  // yet the paper gives
565  // P11(i, j*dim+k) = sum_{nodes n_l in e_i intersected with A_j} 0.5 * phi_k(e_i)
566  // where phi_k is the k-th nullspace vector.
567  //
568  // The graph of D0 contains the incidence from nodes to edges.
569  // The nodal prolongator P maps aggregates to nodes.
570 
571  const SC SC_ZERO = Teuchos::ScalarTraits<SC>::zero();
572  const SC SC_ONE = Teuchos::ScalarTraits<SC>::one();
573  size_t dim = Nullspace_->getNumVectors();
574  size_t numLocalRows = SM_Matrix_->getNodeNumRows();
575 
576  // build prolongator: algorithm 1 in the reference paper
577  // First, build nodal unsmoothed prolongator using the matrix A_nodal
578  RCP<Matrix> P_nodal;
579  bool read_P_from_file = parameterList_.get("refmaxwell: read_P_from_file",false);
580  if (read_P_from_file) {
581  // This permits to read in an ML prolongator, so that we get the same hierarchy.
582  // (ML and MueLu typically produce different aggregates.)
583  std::string P_filename = parameterList_.get("refmaxwell: P_filename",std::string("P.mat"));
584  P_nodal = Xpetra::IO<SC, LO, GO, NO>::Read(P_filename, A_nodal_Matrix_->getDomainMap());
585  } else {
586  Level fineLevel, coarseLevel;
587  fineLevel.SetFactoryManager(Teuchos::null);
588  coarseLevel.SetFactoryManager(Teuchos::null);
589  coarseLevel.SetPreviousLevel(rcpFromRef(fineLevel));
590  fineLevel.SetLevelID(0);
591  coarseLevel.SetLevelID(1);
592  fineLevel.Set("A",A_nodal_Matrix_);
593  fineLevel.Set("Coordinates",Coords_);
594  fineLevel.Set("DofsPerNode",1);
595  coarseLevel.setlib(A_nodal_Matrix_->getDomainMap()->lib());
596  fineLevel.setlib(A_nodal_Matrix_->getDomainMap()->lib());
597  coarseLevel.setObjectLabel("RefMaxwell (1,1)");
598  fineLevel.setObjectLabel("RefMaxwell (1,1)");
599 
600  LocalOrdinal NSdim = 1;
601  RCP<MultiVector> nullSpace = MultiVectorFactory::Build(A_nodal_Matrix_->getRowMap(),NSdim);
602  nullSpace->putScalar(SC_ONE);
603  fineLevel.Set("Nullspace",nullSpace);
604 
606 #ifdef HAVE_MUELU_KOKKOS_REFACTOR
607  RCP<Factory> dropFact, UncoupledAggFact, coarseMapFact, TentativePFact, Tfact;
608  if (useKokkos_) {
609  dropFact = rcp(new CoalesceDropFactory_kokkos());
610  UncoupledAggFact = rcp(new UncoupledAggregationFactory_kokkos());
611  coarseMapFact = rcp(new CoarseMapFactory_kokkos());
612  TentativePFact = rcp(new TentativePFactory_kokkos());
613  Tfact = rcp(new CoordinatesTransferFactory_kokkos());
614  } else {
615  dropFact = rcp(new CoalesceDropFactory());
616  UncoupledAggFact = rcp(new UncoupledAggregationFactory());
617  coarseMapFact = rcp(new CoarseMapFactory());
618  TentativePFact = rcp(new TentativePFactory());
619  Tfact = rcp(new CoordinatesTransferFactory());
620  }
621 #else
624  RCP<CoarseMapFactory> coarseMapFact = rcp(new CoarseMapFactory());
625  RCP<TentativePFactory> TentativePFact = rcp(new TentativePFactory());
627 #endif
628  dropFact->SetFactory("UnAmalgamationInfo", amalgFact);
629  double dropTol = parameterList_.get("aggregation: drop tol",0.0);
630  dropFact->SetParameter("aggregation: drop tol",Teuchos::ParameterEntry(dropTol));
631 
632  UncoupledAggFact->SetFactory("Graph", dropFact);
633 
634  coarseMapFact->SetFactory("Aggregates", UncoupledAggFact);
635 
636  TentativePFact->SetFactory("Aggregates", UncoupledAggFact);
637  TentativePFact->SetFactory("UnAmalgamationInfo", amalgFact);
638  TentativePFact->SetFactory("CoarseMap", coarseMapFact);
639 
640  Tfact->SetFactory("Aggregates", UncoupledAggFact);
641  Tfact->SetFactory("CoarseMap", coarseMapFact);
642 
643  coarseLevel.Request("P",TentativePFact.get());
644  coarseLevel.Request("Coordinates",Tfact.get());
645 
646  coarseLevel.Get("P",P_nodal,TentativePFact.get());
647  coarseLevel.Get("Coordinates",CoordsH_,Tfact.get());
648  }
649  if (dump_matrices_)
650  Xpetra::IO<SC, LO, GO, NO>::Write(std::string("P_nodal.mat"), *P_nodal);
651 
652  RCP<CrsMatrix> D0Crs = rcp_dynamic_cast<CrsMatrixWrap>(D0_Matrix_)->getCrsMatrix();
653 
654  // Import off-rank rows of P_nodal into P_nodal_imported
655  RCP<CrsMatrix> P_nodal_imported;
656  int numProcs = P_nodal->getDomainMap()->getComm()->getSize();
657  if (numProcs > 1) {
658  RCP<CrsMatrixWrap> P_nodal_temp;
659  RCP<const Map> targetMap = D0Crs->getColMap();
660  P_nodal_temp = Teuchos::rcp(new CrsMatrixWrap(targetMap, 0));
661  RCP<const Import> importer = D0Crs->getCrsGraph()->getImporter();
662  P_nodal_temp->doImport(*P_nodal, *importer, Xpetra::INSERT);
663  P_nodal_temp->fillComplete(rcp_dynamic_cast<CrsMatrixWrap>(P_nodal)->getCrsMatrix()->getDomainMap(),
664  rcp_dynamic_cast<CrsMatrixWrap>(P_nodal)->getCrsMatrix()->getRangeMap());
665  P_nodal_imported = P_nodal_temp->getCrsMatrix();
666  if (dump_matrices_)
667  Xpetra::IO<SC, LO, GO, NO>::Write(std::string("P_nodal_imported.mat"), *P_nodal_temp);
668  } else
669  P_nodal_imported = rcp_dynamic_cast<CrsMatrixWrap>(P_nodal)->getCrsMatrix();
670 
671 #ifdef HAVE_MUELU_KOKKOS_REFACTOR
672  if (useKokkos_) {
673 
674  using ATS = Kokkos::ArithTraits<SC>;
676 
677  typedef typename Matrix::local_matrix_type KCRS;
678  typedef typename KCRS::device_type device_t;
679  typedef typename KCRS::StaticCrsGraphType graph_t;
680  typedef typename graph_t::row_map_type::non_const_type lno_view_t;
681  typedef typename graph_t::entries_type::non_const_type lno_nnz_view_t;
682  typedef typename KCRS::values_type::non_const_type scalar_view_t;
683 
684  // Get data out of P_nodal_imported and D0.
685  auto localP = P_nodal_imported->getLocalMatrix();
686  auto localD0 = D0_Matrix_->getLocalMatrix();
687 
688  // Which algorithm should we use for the construction of the special prolongator?
689  // Option "mat-mat":
690  // Multiply D0 * P_nodal, take graph, blow up the domain space and compute the entries.
691  std::string defaultAlgo = "mat-mat";
692  std::string algo = parameterList_.get("refmaxwell: prolongator compute algorithm",defaultAlgo);
693 
694  if (algo == "mat-mat") {
695  Teuchos::RCP<Matrix> D0_P_nodal = MatrixFactory::Build(SM_Matrix_->getRowMap(),0);
696  Xpetra::MatrixMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Multiply(*D0_Matrix_,false,*P_nodal,false,*D0_P_nodal,true,true);
697 
698  // Get data out of D0*P.
699  auto localD0P = D0_P_nodal->getLocalMatrix();
700 
701  // Create the matrix object
702  RCP<Map> blockColMap = Xpetra::MapFactory<LO,GO,NO>::Build(P_nodal_imported->getColMap(), dim);
703  RCP<Map> blockDomainMap = Xpetra::MapFactory<LO,GO,NO>::Build(P_nodal->getDomainMap(), dim);
704 
705  lno_view_t P11rowptr("P11_rowptr", numLocalRows+1);
706  lno_nnz_view_t P11colind("P11_colind",dim*localD0P.graph.entries.size());
707  scalar_view_t P11vals("P11_vals",dim*localD0P.graph.entries.size());
708 
709  // adjust rowpointer
710  Kokkos::parallel_for("MueLu:RefMaxwell::buildProlongator_adjustRowptr", range_type(0,numLocalRows+1),
711  KOKKOS_LAMBDA(const size_t i) {
712  P11rowptr(i) = dim*localD0P.graph.row_map(i);
713  });
714 
715  // adjust column indices
716  Kokkos::parallel_for("MueLu:RefMaxwell::buildProlongator_adjustColind", range_type(0,localD0P.graph.entries.size()),
717  KOKKOS_LAMBDA(const size_t jj) {
718  for (size_t k = 0; k < dim; k++) {
719  P11colind(dim*jj+k) = dim*localD0P.graph.entries(jj)+k;
720  P11vals(dim*jj+k) = SC_ZERO;
721  }
722  });
723 
724  auto localNullspace = Nullspace_->template getLocalView<device_t>();
725 
726  // enter values
727  if (D0_Matrix_->getNodeMaxNumRowEntries()>2) {
728  // The matrix D0 has too many entries per row.
729  // Therefore we need to check whether its entries are actually non-zero.
730  // This is the case for the matrices built by MiniEM.
731  GetOStream(Warnings0) << "RefMaxwell::buildProlongator(): D0 matrix has more than 2 entries per row. Taking inefficient code path." << std::endl;
732 
734 
735  Kokkos::parallel_for("MueLu:RefMaxwell::buildProlongator_enterValues_D0wZeros", range_type(0,numLocalRows),
736  KOKKOS_LAMBDA(const size_t i) {
737  for (size_t ll = localD0.graph.row_map(i); ll < localD0.graph.row_map(i+1); ll++) {
738  LO l = localD0.graph.entries(ll);
739  SC p = localD0.values(ll);
740  if (ATS::magnitude(p) < tol)
741  continue;
742  for (size_t jj = localP.graph.row_map(l); jj < localP.graph.row_map(l+1); jj++) {
743  LO j = localP.graph.entries(jj);
744  SC v = localP.values(jj);
745  for (size_t k = 0; k < dim; k++) {
746  LO jNew = dim*j+k;
747  SC n = localNullspace(i,k);
748  size_t m;
749  for (m = P11rowptr(i); m < P11rowptr(i+1); m++)
750  if (P11colind(m) == jNew)
751  break;
752 #if defined(HAVE_MUELU_DEBUG) && !defined(HAVE_MUELU_CUDA)
753  TEUCHOS_ASSERT_EQUALITY(P11colind(m),jNew);
754 #endif
755  P11vals(m) += 0.5 * v * n;
756  }
757  }
758  }
759  });
760 
761  } else {
762  Kokkos::parallel_for("MueLu:RefMaxwell::buildProlongator_enterValues", range_type(0,numLocalRows),
763  KOKKOS_LAMBDA(const size_t i) {
764  for (size_t ll = localD0.graph.row_map(i); ll < localD0.graph.row_map(i+1); ll++) {
765  LO l = localD0.graph.entries(ll);
766  for (size_t jj = localP.graph.row_map(l); jj < localP.graph.row_map(l+1); jj++) {
767  LO j = localP.graph.entries(jj);
768  SC v = localP.values(jj);
769  for (size_t k = 0; k < dim; k++) {
770  LO jNew = dim*j+k;
771  SC n = localNullspace(i,k);
772  size_t m;
773  for (m = P11rowptr(i); m < P11rowptr(i+1); m++)
774  if (P11colind(m) == jNew)
775  break;
776 #if defined(HAVE_MUELU_DEBUG) && !defined(HAVE_MUELU_CUDA)
777  TEUCHOS_ASSERT_EQUALITY(P11colind(m),jNew);
778 #endif
779  P11vals(m) += 0.5 * v * n;
780  }
781  }
782  }
783  });
784  }
785 
786  P11_ = Teuchos::rcp(new CrsMatrixWrap(SM_Matrix_->getRowMap(), blockColMap, 0, Xpetra::StaticProfile));
787  RCP<CrsMatrix> P11Crs = rcp_dynamic_cast<CrsMatrixWrap>(P11_)->getCrsMatrix();
788  P11Crs->setAllValues(P11rowptr, P11colind, P11vals);
789  P11Crs->expertStaticFillComplete(blockDomainMap, SM_Matrix_->getRangeMap());
790 
791  } else
792  TEUCHOS_TEST_FOR_EXCEPTION(false,std::invalid_argument,algo << " is not a valid option for \"refmaxwell: prolongator compute algorithm\"");
793  } else {
794 #endif // ifdef(HAVE_MUELU_KOKKOS_REFACTOR)
795 
796  // get nullspace vectors
797  ArrayRCP<ArrayRCP<const SC> > nullspaceRCP(dim);
798  ArrayRCP<ArrayView<const SC> > nullspace(dim);
799  for(size_t i=0; i<dim; i++) {
800  nullspaceRCP[i] = Nullspace_->getData(i);
801  nullspace[i] = nullspaceRCP[i]();
802  }
803 
804  // Get data out of P_nodal_imported and D0.
805  ArrayRCP<const size_t> Prowptr_RCP, D0rowptr_RCP;
806  ArrayRCP<const LO> Pcolind_RCP, D0colind_RCP;
807  ArrayRCP<const SC> Pvals_RCP, D0vals_RCP;
808  ArrayRCP<size_t> P11rowptr_RCP;
809  ArrayRCP<LO> P11colind_RCP;
810  ArrayRCP<SC> P11vals_RCP;
811 
812  P_nodal_imported->getAllValues(Prowptr_RCP, Pcolind_RCP, Pvals_RCP);
813  rcp_dynamic_cast<CrsMatrixWrap>(D0_Matrix_)->getCrsMatrix()->getAllValues(D0rowptr_RCP, D0colind_RCP, D0vals_RCP);
814 
815  // For efficiency
816  // Refers to an issue where Teuchos::ArrayRCP::operator[] may be
817  // slower than Teuchos::ArrayView::operator[].
818  ArrayView<const size_t> Prowptr, D0rowptr;
819  ArrayView<const LO> Pcolind, D0colind;
820  ArrayView<const SC> Pvals, D0vals;
821  Prowptr = Prowptr_RCP(); Pcolind = Pcolind_RCP(); Pvals = Pvals_RCP();
822  D0rowptr = D0rowptr_RCP(); D0colind = D0colind_RCP(); D0vals = D0vals_RCP();
823 
824  // Which algorithm should we use for the construction of the special prolongator?
825  // Option "mat-mat":
826  // Multiply D0 * P_nodal, take graph, blow up the domain space and compute the entries.
827  // Option "gustavson":
828  // Loop over D0, P and nullspace and allocate directly. (Gustavson-like)
829  // More efficient, but only available for serial node.
830  std::string defaultAlgo = "gustavson";
831  std::string algo = parameterList_.get("refmaxwell: prolongator compute algorithm",defaultAlgo);
832 
833  if (algo == "mat-mat") {
834  Teuchos::RCP<Matrix> D0_P_nodal = MatrixFactory::Build(SM_Matrix_->getRowMap(),0);
835  Xpetra::MatrixMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Multiply(*D0_Matrix_,false,*P_nodal,false,*D0_P_nodal,true,true);
836 
837  // Get data out of D0*P.
838  ArrayRCP<const size_t> D0Prowptr_RCP;
839  ArrayRCP<const LO> D0Pcolind_RCP;
840  ArrayRCP<const SC> D0Pvals_RCP;
841  rcp_dynamic_cast<CrsMatrixWrap>(D0_P_nodal)->getCrsMatrix()->getAllValues(D0Prowptr_RCP, D0Pcolind_RCP, D0Pvals_RCP);
842 
843  // For efficiency
844  // Refers to an issue where Teuchos::ArrayRCP::operator[] may be
845  // slower than Teuchos::ArrayView::operator[].
846  ArrayView<const size_t> D0Prowptr;
847  ArrayView<const LO> D0Pcolind;
848  D0Prowptr = D0Prowptr_RCP(); D0Pcolind = D0Pcolind_RCP();
849 
850  // Create the matrix object
851  RCP<Map> blockColMap = Xpetra::MapFactory<LO,GO,NO>::Build(P_nodal_imported->getColMap(), dim);
852  RCP<Map> blockDomainMap = Xpetra::MapFactory<LO,GO,NO>::Build(P_nodal->getDomainMap(), dim);
853  P11_ = Teuchos::rcp(new CrsMatrixWrap(SM_Matrix_->getRowMap(), blockColMap, 0, Xpetra::StaticProfile));
854  RCP<CrsMatrix> P11Crs = rcp_dynamic_cast<CrsMatrixWrap>(P11_)->getCrsMatrix();
855  P11Crs->allocateAllValues(dim*D0Pcolind.size(), P11rowptr_RCP, P11colind_RCP, P11vals_RCP);
856 
857  ArrayView<size_t> P11rowptr = P11rowptr_RCP();
858  ArrayView<LO> P11colind = P11colind_RCP();
859  ArrayView<SC> P11vals = P11vals_RCP();
860 
861  // adjust rowpointer
862  for (size_t i = 0; i < numLocalRows+1; i++) {
863  P11rowptr[i] = dim*D0Prowptr[i];
864  }
865 
866  // adjust column indices
867  size_t nnz = 0;
868  for (size_t jj = 0; jj < (size_t) D0Pcolind.size(); jj++)
869  for (size_t k = 0; k < dim; k++) {
870  P11colind[nnz] = dim*D0Pcolind[jj]+k;
871  P11vals[nnz] = SC_ZERO;
872  nnz++;
873  }
874 
875  // enter values
876  if (D0_Matrix_->getNodeMaxNumRowEntries()>2) {
877  // The matrix D0 has too many entries per row.
878  // Therefore we need to check whether its entries are actually non-zero.
879  // This is the case for the matrices built by MiniEM.
880  GetOStream(Warnings0) << "RefMaxwell::buildProlongator(): D0 matrix has more than 2 entries per row. Taking inefficient code path." << std::endl;
881 
883  for (size_t i = 0; i < numLocalRows; i++) {
884  for (size_t ll = D0rowptr[i]; ll < D0rowptr[i+1]; ll++) {
885  LO l = D0colind[ll];
886  SC p = D0vals[ll];
888  continue;
889  for (size_t jj = Prowptr[l]; jj < Prowptr[l+1]; jj++) {
890  LO j = Pcolind[jj];
891  SC v = Pvals[jj];
892  for (size_t k = 0; k < dim; k++) {
893  LO jNew = dim*j+k;
894  SC n = nullspace[k][i];
895  size_t m;
896  for (m = P11rowptr[i]; m < P11rowptr[i+1]; m++)
897  if (P11colind[m] == jNew)
898  break;
899 #ifdef HAVE_MUELU_DEBUG
900  TEUCHOS_ASSERT_EQUALITY(P11colind[m],jNew);
901 #endif
902  P11vals[m] += 0.5 * v * n;
903  }
904  }
905  }
906  }
907  } else {
908  // enter values
909  for (size_t i = 0; i < numLocalRows; i++) {
910  for (size_t ll = D0rowptr[i]; ll < D0rowptr[i+1]; ll++) {
911  LO l = D0colind[ll];
912  for (size_t jj = Prowptr[l]; jj < Prowptr[l+1]; jj++) {
913  LO j = Pcolind[jj];
914  SC v = Pvals[jj];
915  for (size_t k = 0; k < dim; k++) {
916  LO jNew = dim*j+k;
917  SC n = nullspace[k][i];
918  size_t m;
919  for (m = P11rowptr[i]; m < P11rowptr[i+1]; m++)
920  if (P11colind[m] == jNew)
921  break;
922 #ifdef HAVE_MUELU_DEBUG
923  TEUCHOS_ASSERT_EQUALITY(P11colind[m],jNew);
924 #endif
925  P11vals[m] += 0.5 * v * n;
926  }
927  }
928  }
929  }
930  }
931 
932  P11Crs->setAllValues(P11rowptr_RCP, P11colind_RCP, P11vals_RCP);
933  P11Crs->expertStaticFillComplete(blockDomainMap, SM_Matrix_->getRangeMap());
934 
935  } else if (algo == "gustavson") {
936 
937  LO maxP11col = dim * P_nodal_imported->getColMap()->getMaxLocalIndex();
938  const size_t ST_INVALID = Teuchos::OrdinalTraits<LO>::invalid();
939  Array<size_t> P11_status(dim*maxP11col, ST_INVALID);
940  // This is ad-hoc and should maybe be replaced with some better heuristics.
941  size_t nnz_alloc = dim*D0vals_RCP.size();
942 
943  // Create the matrix object
944  RCP<Map> blockColMap = Xpetra::MapFactory<LO,GO,NO>::Build(P_nodal_imported->getColMap(), dim);
945  RCP<Map> blockDomainMap = Xpetra::MapFactory<LO,GO,NO>::Build(P_nodal->getDomainMap(), dim);
946  P11_ = Teuchos::rcp(new CrsMatrixWrap(SM_Matrix_->getRowMap(), blockColMap, 0, Xpetra::StaticProfile));
947  RCP<CrsMatrix> P11Crs = rcp_dynamic_cast<CrsMatrixWrap>(P11_)->getCrsMatrix();
948  P11Crs->allocateAllValues(nnz_alloc, P11rowptr_RCP, P11colind_RCP, P11vals_RCP);
949 
950  ArrayView<size_t> P11rowptr = P11rowptr_RCP();
951  ArrayView<LO> P11colind = P11colind_RCP();
952  ArrayView<SC> P11vals = P11vals_RCP();
953 
954  size_t nnz;
955  if (D0_Matrix_->getNodeMaxNumRowEntries()>2) {
956  // The matrix D0 has too many entries per row.
957  // Therefore we need to check whether its entries are actually non-zero.
958  // This is the case for the matrices built by MiniEM.
959  GetOStream(Warnings0) << "RefMaxwell::buildProlongator(): D0 matrix has more than 2 entries per row. Taking inefficient code path." << std::endl;
960 
962  nnz = 0;
963  size_t nnz_old = 0;
964  for (size_t i = 0; i < numLocalRows; i++) {
965  P11rowptr[i] = nnz;
966  for (size_t ll = D0rowptr[i]; ll < D0rowptr[i+1]; ll++) {
967  LO l = D0colind[ll];
968  SC p = D0vals[ll];
970  continue;
971  for (size_t jj = Prowptr[l]; jj < Prowptr[l+1]; jj++) {
972  LO j = Pcolind[jj];
973  SC v = Pvals[jj];
974  for (size_t k = 0; k < dim; k++) {
975  LO jNew = dim*j+k;
976  SC n = nullspace[k][i];
977  // do we already have an entry for (i, jNew)?
978  if (P11_status[jNew] == ST_INVALID || P11_status[jNew] < nnz_old) {
979  P11_status[jNew] = nnz;
980  P11colind[nnz] = jNew;
981  P11vals[nnz] = 0.5 * v * n;
982  // or should it be
983  // P11vals[nnz] = 0.5 * n;
984  nnz++;
985  } else {
986  P11vals[P11_status[jNew]] += 0.5 * v * n;
987  // or should it be
988  // P11vals[P11_status[jNew]] += 0.5 * n;
989  }
990  }
991  }
992  }
993  nnz_old = nnz;
994  }
995  P11rowptr[numLocalRows] = nnz;
996  } else {
997  nnz = 0;
998  size_t nnz_old = 0;
999  for (size_t i = 0; i < numLocalRows; i++) {
1000  P11rowptr[i] = nnz;
1001  for (size_t ll = D0rowptr[i]; ll < D0rowptr[i+1]; ll++) {
1002  LO l = D0colind[ll];
1003  for (size_t jj = Prowptr[l]; jj < Prowptr[l+1]; jj++) {
1004  LO j = Pcolind[jj];
1005  SC v = Pvals[jj];
1006  for (size_t k = 0; k < dim; k++) {
1007  LO jNew = dim*j+k;
1008  SC n = nullspace[k][i];
1009  // do we already have an entry for (i, jNew)?
1010  if (P11_status[jNew] == ST_INVALID || P11_status[jNew] < nnz_old) {
1011  P11_status[jNew] = nnz;
1012  P11colind[nnz] = jNew;
1013  P11vals[nnz] = 0.5 * v * n;
1014  // or should it be
1015  // P11vals[nnz] = 0.5 * n;
1016  nnz++;
1017  } else {
1018  P11vals[P11_status[jNew]] += 0.5 * v * n;
1019  // or should it be
1020  // P11vals[P11_status[jNew]] += 0.5 * n;
1021  }
1022  }
1023  }
1024  }
1025  nnz_old = nnz;
1026  }
1027  P11rowptr[numLocalRows] = nnz;
1028  }
1029 
1030  if (blockDomainMap->lib() == Xpetra::UseTpetra) {
1031  // Downward resize
1032  // - Cannot resize for Epetra, as it checks for same pointers
1033  // - Need to resize for Tpetra, as it checks ().size() == P11rowptr[numLocalRows]
1034  P11vals_RCP.resize(nnz);
1035  P11colind_RCP.resize(nnz);
1036  }
1037 
1038  P11Crs->setAllValues(P11rowptr_RCP, P11colind_RCP, P11vals_RCP);
1039  P11Crs->expertStaticFillComplete(blockDomainMap, SM_Matrix_->getRangeMap());
1040  } else
1041  TEUCHOS_TEST_FOR_EXCEPTION(false,std::invalid_argument,algo << " is not a valid option for \"refmaxwell: prolongator compute algorithm\"");
1042 #ifdef HAVE_MUELU_KOKKOS_REFACTOR
1043  }
1044 #endif
1045  }
1046 
1047  template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1049  Teuchos::TimeMonitor tm(*Teuchos::TimeMonitor::getNewTimer("MueLu RefMaxwell: Build coarse (1,1) matrix"));
1050 
1051  // coarse matrix for P11* (M1 + D1* M2 D1) P11
1052  Teuchos::RCP<Matrix> Matrix1 = MatrixFactory::Build(P11_->getDomainMap(),0);
1053  if (parameterList_.get<bool>("rap: triple product", false) == false) {
1054  Teuchos::RCP<Matrix> C = MatrixFactory::Build(SM_Matrix_->getRowMap(),0);
1055  // construct (M1 + D1* M2 D1) P11
1056  Xpetra::MatrixMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Multiply(*SM_Matrix_,false,*P11_,false,*C,true,true);
1057  // construct P11* (M1 + D1* M2 D1) P11
1058  Xpetra::MatrixMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Multiply(*R11_,false,*C,false,*Matrix1,true,true);
1059  } else {
1061  MultiplyRAP(*P11_, true, *SM_Matrix_, false, *P11_, false, *Matrix1, true, true);
1062  }
1063  if (parameterList_.get<bool>("rap: fix zero diagonals", true))
1065 
1066  if(disable_addon_==true) {
1067  // if add-on is not chosen
1068  AH_=Matrix1;
1069  }
1070  else {
1071  Teuchos::TimeMonitor tmAddon(*Teuchos::TimeMonitor::getNewTimer("MueLu RefMaxwell: Build coarse addon matrix"));
1072  // catch a failure
1073  TEUCHOS_TEST_FOR_EXCEPTION(M0inv_Matrix_==Teuchos::null,std::invalid_argument,
1074  "MueLu::RefMaxwell::formCoarseMatrix(): Inverse of "
1075  "lumped mass matrix required for add-on (i.e. M0inv_Matrix is null)");
1076 
1077  // coarse matrix for add-on, i.e P11* (M1 D0 M0inv D0* M1) P11
1078  Teuchos::RCP<Matrix> Zaux = MatrixFactory::Build(M1_Matrix_->getRowMap(),0);
1079  Teuchos::RCP<Matrix> Z = MatrixFactory::Build(D0_Matrix_->getDomainMap(),0);
1080 
1081  // construct Zaux = M1 P11
1082  Xpetra::MatrixMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Multiply(*M1_Matrix_,false,*P11_,false,*Zaux,true,true);
1083  // construct Z = D0* M1 P11 = D0* Zaux
1084  Xpetra::MatrixMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Multiply(*D0_Matrix_,true,*Zaux,false,*Z,true,true);
1085 
1086  // construct Z* M0inv Z
1087  Teuchos::RCP<Matrix> Matrix2 = MatrixFactory::Build(Z->getDomainMap(),0);
1088  if (M0inv_Matrix_->getGlobalMaxNumRowEntries()<=1) {
1089  // We assume that if M0inv has at most one entry per row then
1090  // these are all diagonal entries.
1091 #ifdef HAVE_MUELU_KOKKOS_REFACTOR
1092  RCP<Matrix> ZT;
1093  if (useKokkos_)
1094  ZT = Utilities_kokkos::Transpose(*Z);
1095  else
1096  ZT = Utilities::Transpose(*Z);
1097 #else
1099 #endif
1100  RCP<Vector> diag = VectorFactory::Build(M0inv_Matrix_->getRowMap());
1101  M0inv_Matrix_->getLocalDiagCopy(*diag);
1102  if (Z->getRowMap()->isSameAs(*(diag->getMap())))
1103  Z->leftScale(*diag);
1104  else {
1105  RCP<Import> importer = ImportFactory::Build(diag->getMap(),Z->getRowMap());
1106  RCP<Vector> diag2 = VectorFactory::Build(Z->getRowMap());
1107  diag2->doImport(*diag,*importer,Xpetra::INSERT);
1108  Z->leftScale(*diag2);
1109  }
1110  Xpetra::MatrixMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Multiply(*ZT,false,*Z,false,*Matrix2,true,true);
1111  } else if (parameterList_.get<bool>("rap: triple product", false) == false) {
1112  Teuchos::RCP<Matrix> C2 = MatrixFactory::Build(M0inv_Matrix_->getRowMap(),0);
1113  // construct C2 = M0inv Z
1114  Xpetra::MatrixMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Multiply(*M0inv_Matrix_,false,*Z,false,*C2,true,true);
1115  // construct Matrix2 = Z* M0inv Z = Z* C2
1116  Xpetra::MatrixMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Multiply(*Z,true,*C2,false,*Matrix2,true,true);
1117  } else {
1118  // construct Matrix2 = Z* M0inv Z
1120  MultiplyRAP(*Z, true, *M0inv_Matrix_, false, *Z, false, *Matrix2, true, true);
1121  }
1122  // add matrices together
1123  Xpetra::MatrixMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>::TwoMatrixAdd(*Matrix1,false,(Scalar)1.0,*Matrix2,false,(Scalar)1.0,AH_,GetOStream(Runtime0));
1124  AH_->fillComplete();
1125  }
1126 
1127  // set fixed block size for vector nodal matrix
1128  size_t dim = Nullspace_->getNumVectors();
1129  AH_->SetFixedBlockSize(dim);
1130  AH_->setObjectLabel("RefMaxwell (1,1)");
1131 
1132  }
1133 
1134 
1135  template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1137  SM_Matrix_ = SM_Matrix_new;
1138  }
1139 
1140 
1141  template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1143 
1144  // compute residuals
1145  Scalar one = Teuchos::ScalarTraits<Scalar>::one(), negone = -one, zero = Teuchos::ScalarTraits<Scalar>::zero();
1146  SM_Matrix_->apply(X, *residual_, Teuchos::NO_TRANS, one, zero);
1147  residual_->update(one, RHS, negone);
1148  R11_->apply(*residual_,*P11res_,Teuchos::NO_TRANS);
1149  D0_T_Matrix_->apply(*residual_,*D0res_,Teuchos::NO_TRANS);
1150 
1151  // block diagonal preconditioner on 2x2 (V-cycle for diagonal blocks)
1152  HierarchyH_->Iterate(*P11res_, *P11x_, 1, true);
1153  Hierarchy22_->Iterate(*D0res_, *D0x_, 1, true);
1154 
1155  // update current solution
1156  P11_->apply(*P11x_,*residual_,Teuchos::NO_TRANS);
1157  D0_Matrix_->apply(*D0x_,*residual_,Teuchos::NO_TRANS,one,one);
1158  X.update(one, *residual_, one);
1159 
1160  }
1161 
1162 
1163  template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1165 
1166  // precondition (1,1)-block
1167  Scalar one = Teuchos::ScalarTraits<Scalar>::one(), negone = -one, zero = Teuchos::ScalarTraits<Scalar>::zero();
1168  SM_Matrix_->apply(X, *residual_, Teuchos::NO_TRANS, one, zero);
1169  residual_->update(one, RHS, negone);
1170  R11_->apply(*residual_,*P11res_,Teuchos::NO_TRANS);
1171  HierarchyH_->Iterate(*P11res_, *P11x_, 1, true);
1172  P11_->apply(*P11x_,*residual_,Teuchos::NO_TRANS);
1173  X.update(one, *residual_, one);
1174 
1175  // precondition (2,2)-block
1176  SM_Matrix_->apply(X, *residual_, Teuchos::NO_TRANS, one, zero);
1177  residual_->update(one, RHS, negone);
1178  D0_T_Matrix_->apply(*residual_,*D0res_,Teuchos::NO_TRANS);
1179  Hierarchy22_->Iterate(*D0res_, *D0x_, 1, true);
1180  D0_Matrix_->apply(*D0x_,*residual_,Teuchos::NO_TRANS);
1181  X.update(one, *residual_, one);
1182 
1183  // precondition (1,1)-block
1184  SM_Matrix_->apply(X, *residual_, Teuchos::NO_TRANS, one, zero);
1185  residual_->update(one, RHS, negone);
1186  R11_->apply(*residual_,*P11res_,Teuchos::NO_TRANS);
1187  HierarchyH_->Iterate(*P11res_, *P11x_, 1, true);
1188  P11_->apply(*P11x_,*residual_,Teuchos::NO_TRANS);
1189  X.update(one, *residual_, one);
1190 
1191  }
1192 
1193 
1194  template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1196 
1197  // precondition (2,2)-block
1198  Scalar one = Teuchos::ScalarTraits<Scalar>::one(), negone = -one, zero = Teuchos::ScalarTraits<Scalar>::zero();
1199  SM_Matrix_->apply(X, *residual_, Teuchos::NO_TRANS, one, zero);
1200  residual_->update(one, RHS, negone);
1201  D0_T_Matrix_->apply(*residual_,*D0res_,Teuchos::NO_TRANS);
1202  Hierarchy22_->Iterate(*D0res_, *D0x_, 1, true);
1203  D0_Matrix_->apply(*D0x_,*residual_,Teuchos::NO_TRANS);
1204  X.update(one, *residual_, one);
1205 
1206  // precondition (1,1)-block
1207  SM_Matrix_->apply(X, *residual_, Teuchos::NO_TRANS, one, zero);
1208  residual_->update(one, RHS, negone);
1209  R11_->apply(*residual_,*P11res_,Teuchos::NO_TRANS);
1210  HierarchyH_->Iterate(*P11res_, *P11x_, 1, true);
1211  P11_->apply(*P11x_,*residual_,Teuchos::NO_TRANS);
1212  X.update(one, *residual_, one);
1213 
1214  // precondition (2,2)-block
1215  SM_Matrix_->apply(X, *residual_, Teuchos::NO_TRANS, one, zero);
1216  residual_->update(one, RHS, negone);
1217  D0_T_Matrix_->apply(*residual_,*D0res_,Teuchos::NO_TRANS);
1218  Hierarchy22_->Iterate(*D0res_, *D0x_, 1, true);
1219  D0_Matrix_->apply(*D0x_,*residual_,Teuchos::NO_TRANS);
1220  X.update(one, *residual_, one);
1221 
1222  }
1223 
1224  template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1226 
1227  // compute residuals
1228  Scalar one = Teuchos::ScalarTraits<Scalar>::one(), negone = -one, zero = Teuchos::ScalarTraits<Scalar>::zero();
1229  SM_Matrix_->apply(X, *residual_, Teuchos::NO_TRANS, one, zero);
1230  residual_->update(one, RHS, negone);
1231  R11_->apply(*residual_,*P11res_,Teuchos::NO_TRANS);
1232 
1233  // block diagonal preconditioner on 2x2 (V-cycle for diagonal blocks)
1234  HierarchyH_->Iterate(*P11res_, *P11x_, 1, true);
1235 
1236  // update current solution
1237  P11_->apply(*P11x_,*residual_,Teuchos::NO_TRANS);
1238  X.update(one, *residual_, one);
1239 
1240  }
1241 
1242 
1243  template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1245  Teuchos::ETransp mode,
1246  Scalar alpha,
1247  Scalar beta) const {
1248 
1249  Teuchos::TimeMonitor tm(*Teuchos::TimeMonitor::getNewTimer("MueLuRefmaxwell: Solve"));
1250 
1251  // make sure that we have enough temporary memory
1252  if (X.getNumVectors() != P11res_->getNumVectors()) {
1253  P11res_ = MultiVectorFactory::Build(P11_->getDomainMap(),X.getNumVectors());
1254  P11x_ = MultiVectorFactory::Build(P11_->getDomainMap(),X.getNumVectors());
1255  D0res_ = MultiVectorFactory::Build(D0_Matrix_->getDomainMap(),X.getNumVectors());
1256  D0x_ = MultiVectorFactory::Build(D0_Matrix_->getDomainMap(),X.getNumVectors());
1257  residual_ = MultiVectorFactory::Build(SM_Matrix_->getDomainMap(),1);
1258  }
1259 
1260  // apply pre-smoothing
1261 #if defined(HAVE_MUELU_IFPACK2) && (!defined(HAVE_MUELU_EPETRA) || defined(HAVE_MUELU_INST_DOUBLE_INT_INT))
1262  if (useHiptmairSmoothing_) {
1263  Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> tX = Utilities::MV2NonConstTpetraMV(X);
1264  Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> tRHS = Utilities::MV2TpetraMV(RHS);
1265  hiptmairPreSmoother_->apply(tRHS, tX);
1266  }
1267  else
1268 #endif
1269  PreSmoother_->Apply(X, RHS, use_as_preconditioner_);
1270 
1271  // do solve for the 2x2 block system
1272  if(mode_=="additive")
1273  applyInverseAdditive(RHS,X);
1274  else if(mode_=="121")
1275  applyInverse121(RHS,X);
1276  else if(mode_=="212")
1277  applyInverse212(RHS,X);
1278  else if(mode_=="11only")
1279  applyInverse11only(RHS,X);
1280  else if(mode_=="none") {
1281  // do nothing
1282  }
1283  else
1284  applyInverseAdditive(RHS,X);
1285 
1286  // apply post-smoothing
1287 #if defined(HAVE_MUELU_IFPACK2) && (!defined(HAVE_MUELU_EPETRA) || defined(HAVE_MUELU_INST_DOUBLE_INT_INT))
1288  if (useHiptmairSmoothing_)
1289  {
1290  Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> tX = Utilities::MV2NonConstTpetraMV(X);
1291  Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> tRHS = Utilities::MV2TpetraMV(RHS);
1292  hiptmairPostSmoother_->apply(tRHS, tX);
1293  }
1294  else
1295 #endif
1296  PostSmoother_->Apply(X, RHS, false);
1297 
1298  }
1299 
1300 
1301  template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1303  return false;
1304  }
1305 
1306 
1307  template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1310  const Teuchos::RCP<Matrix> & M0inv_Matrix,
1311  const Teuchos::RCP<Matrix> & M1_Matrix,
1312  const Teuchos::RCP<MultiVector> & Nullspace,
1313  const Teuchos::RCP<RealValuedMultiVector> & Coords,
1314  Teuchos::ParameterList& List)
1315  {
1316  // some pre-conditions
1317  TEUCHOS_ASSERT(D0_Matrix!=Teuchos::null);
1318  TEUCHOS_ASSERT(M1_Matrix!=Teuchos::null);
1319 
1320  HierarchyH_ = Teuchos::null;
1321  Hierarchy22_ = Teuchos::null;
1322  PreSmoother_ = Teuchos::null;
1323  PostSmoother_ = Teuchos::null;
1324  disable_addon_ = false;
1325  mode_ = "additive";
1326 
1327  // set parameters
1328  setParameters(List);
1329 
1330  D0_Matrix_ = D0_Matrix;
1331  M0inv_Matrix_ = M0inv_Matrix;
1332  M1_Matrix_ = M1_Matrix;
1333  Coords_ = Coords;
1334  Nullspace_ = Nullspace;
1335  }
1336 
1337  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1340  out << "\n--------------------------------------------------------------------------------\n" <<
1341  "--- RefMaxwell Summary ---\n"
1342  "--------------------------------------------------------------------------------" << std::endl;
1343  out << std::endl;
1344 
1345  GlobalOrdinal numRows;
1346  GlobalOrdinal nnz;
1347 
1348  numRows = SM_Matrix_->getGlobalNumRows();
1349  nnz = SM_Matrix_->getGlobalNumEntries();
1350 
1351  Xpetra::global_size_t tt = numRows;
1352  int rowspacer = 3; while (tt != 0) { tt /= 10; rowspacer++; }
1353  tt = nnz;
1354  int nnzspacer = 2; while (tt != 0) { tt /= 10; nnzspacer++; }
1355 
1356  out << "block " << std::setw(rowspacer) << " rows " << std::setw(nnzspacer) << " nnz " << std::setw(9) << " nnz/row" << std::endl;
1357  out << "(1, 1)" << std::setw(rowspacer) << numRows << std::setw(nnzspacer) << nnz << std::setw(9) << as<double>(nnz) / numRows << std::endl;
1358 
1359  numRows = A22_->getGlobalNumRows();
1360  nnz = A22_->getGlobalNumEntries();
1361 
1362  out << "(2, 2)" << std::setw(rowspacer) << numRows << std::setw(nnzspacer) << nnz << std::setw(9) << as<double>(nnz) / numRows << std::endl;
1363 
1364  out << std::endl;
1365 
1366  }
1367 
1368 } // namespace
1369 
1370 #define MUELU_REFMAXWELL_SHORT
1371 #endif //ifdef MUELU_REFMAXWELL_DEF_HPP
Important warning messages (one line)
Generic Smoother Factory for generating the smoothers of the MG hierarchy.
This class specifies the default factory that should generate some data on a Level if the data does n...
static void CheckRepairMainDiagonal(RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node >> &Ac, bool const &repairZeroDiagonals, Teuchos::FancyOStream &fos)
Various adapters that will create a MueLu preconditioner that is an Xpetra::Matrix.
void parallel_for(const ExecPolicy &policy, const FunctorType &functor, const std::string &str="", typename Impl::enable_if< Kokkos::Impl::is_execution_policy< ExecPolicy >::value >::type *=0)
T & Get(const std::string &ename, const FactoryBase *factory=NoFactory::get())
Get data without decrementing associated storage counter (i.e., read-only access). Usage: Level->Get< RCP<Matrix> >("A", factory) if factory == NULL => use default factory.
static Teuchos::ArrayRCP< const bool > DetectDirichletCols(const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A, const Teuchos::ArrayRCP< const bool > &dirichletRows)
Detect Dirichlet columns based on Dirichlet rows.
Factory for generating coarse level map. Used by TentativePFactory.
static RCP< const Tpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > MV2TpetraMV(RCP< Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > const vec)
Helper utility to pull out the underlying Tpetra objects from an Xpetra object.
static magnitudeType eps()
size_type size() const
T & get(const std::string &name, T def_value)
Class that encapsulates external library smoothers.
ParameterList & set(std::string const &name, T const &value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
static void MultiplyRAP(const Matrix &R, bool transposeR, const Matrix &A, bool transposeA, const Matrix &P, bool transposeP, Matrix &Ac, bool call_FillComplete_on_result=true, bool doOptimizeStorage=true, const std::string &label=std::string(), const RCP< ParameterList > &params=null)
static Teuchos::ArrayRCP< const bool > DetectDirichletRows(const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A, const Magnitude &tol=Teuchos::ScalarTraits< Scalar >::magnitude(0.), const bool count_twos_as_dirichlet=false)
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
static void Write(const std::string &fileName, const Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node > &M)
bool isSublist(const std::string &name) const
void ZeroDirichletCols(RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > &A, const Kokkos::View< const bool *, typename Node::device_type > &dirichletCols, Scalar replaceWith)
Teuchos::RCP< MueLu::Hierarchy< Scalar, LocalOrdinal, GlobalOrdinal, Node > > CreateXpetraPreconditioner(Teuchos::RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > op, const Teuchos::ParameterList &inParamList, Teuchos::RCP< Xpetra::MultiVector< double, LocalOrdinal, GlobalOrdinal, Node > > coords=Teuchos::null, Teuchos::RCP< Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > nullspace=Teuchos::null)
Helper function to create a MueLu preconditioner that can be used by Xpetra.Given an Xpetra::Matrix...
void ZeroDirichletRows(RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > &A, const Kokkos::View< const bool *, typename Node::device_type > &dirichletRows, Scalar replaceWith)
LocalOrdinal LO
One-liner description of what is happening.
size_type size() const
Namespace for MueLu classes and methods.
void SetPreviousLevel(const RCP< Level > &previousLevel)
Definition: MueLu_Level.cpp:85
void SetFactoryManager(const RCP< const FactoryManagerBase > &factoryManager)
Set default factories (used internally by Hierarchy::SetLevel()).
Definition: MueLu_Level.cpp:92
void setParameters(Teuchos::ParameterList &list)
Set parameters.
static RCP< Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > RealValuedToScalarMultiVector(RCP< Xpetra::MultiVector< double, LocalOrdinal, GlobalOrdinal, Node > > X)
void applyInverseAdditive(const MultiVector &RHS, MultiVector &X) const
apply additive algorithm for 2x2 solve
void buildProlongator()
Setup the prolongator for the (1,1)-block.
static RCP< Tpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > MV2NonConstTpetraMV(RCP< Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > vec)
Teuchos::ScalarTraits< Scalar >::magnitudeType magnitudeType
T * get() const
Factory for building tentative prolongator.
virtual void update(const Scalar &alpha, const MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A, const Scalar &beta)=0
void setlib(Xpetra::UnderlyingLib lib2)
static RCP< Time > getNewTimer(const std::string &name)
RCP< const ParameterList > GetValidParameterList() const
Return a const parameter list of valid parameters that setParameterList() will accept.
Additional warnings.
void resize(const size_type n, const T &val=T())
void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel=Teuchos::VERB_HIGH) const
void applyInverse121(const MultiVector &RHS, MultiVector &X) const
apply 1-2-1 algorithm for 2x2 solve
void Build(Level &currentLevel) const
Build an object with this factory.
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
virtual void SetParameterList(const ParameterList &paramList)
Set parameters from a parameter list and return with default values.
Class that holds all level-specific information.
Definition: MueLu_Level.hpp:99
void compute()
Setup the preconditioner.
Teuchos::RCP< const Map > getRangeMap() const
Returns the Xpetra::Map object associated with the range of this operator.
static RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Transpose(Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &Op, bool optimizeTranspose=false, const std::string &label=std::string(), const Teuchos::RCP< Teuchos::ParameterList > &params=Teuchos::null)
virtual void setObjectLabel(const std::string &objectLabel)
AmalgamationFactory for subblocks of strided map based amalgamation data.
void initialize(const Teuchos::RCP< Matrix > &D0_Matrix, const Teuchos::RCP< Matrix > &M0inv_Matrix, const Teuchos::RCP< Matrix > &M1_Matrix, const Teuchos::RCP< MultiVector > &Nullspace, const Teuchos::RCP< RealValuedMultiVector > &Coords, Teuchos::ParameterList &List)
static void ZeroDirichletRows(Teuchos::RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > &A, const std::vector< LocalOrdinal > &dirichletRows, Scalar replaceWith=Teuchos::ScalarTraits< Scalar >::zero())
size_t global_size_t
void Set(const std::string &ename, const T &entry, const FactoryBase *factory=NoFactory::get())
static void TwoMatrixAdd(const Matrix &A, bool transposeA, SC alpha, Matrix &B, SC beta)
void SetParameter(const std::string &name, const ParameterEntry &entry)
Set a parameter directly as a ParameterEntry.
bool hasTransposeApply() const
Indicates whether this operator supports applying the adjoint operator.
void Build(Level &currentLevel) const
Creates pre and post smoothers.
void applyInverse212(const MultiVector &RHS, MultiVector &X) const
apply 2-1-2 algorithm for 2x2 solve
Factory for creating a graph base on a given matrix.
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 SetLevelID(int levelID)
Set level number.
Definition: MueLu_Level.cpp:78
void deep_copy(const View< DT, DP... > &dst, typename ViewTraits< DT, DP... >::const_value_type &value, typename std::enable_if< std::is_same< typename ViewTraits< DT, DP... >::specialize, void >::value >::type *=0)
static Teuchos::RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Read(const std::string &fileName, Xpetra::UnderlyingLib lib, const RCP< const Teuchos::Comm< int > > &comm, bool binary=false)
Scalar SC
Class for transferring coordinates from a finer level to a coarser one.
Factory for building restriction operators.
Kokkos::View< const bool *, typename NO::device_type > DetectDirichletCols(const Xpetra::Matrix< SC, LO, GO, NO > &A, const Kokkos::View< const bool *, typename NO::device_type > &dirichletRows)
ParameterList & sublist(const std::string &name, bool mustAlreadyExist=false, const std::string &docString="")
virtual void SetFactory(const std::string &varName, const RCP< const FactoryBase > &factory)
Configuration.
void formCoarseMatrix()
Compute P11^{T}*A*P11 efficiently.
void applyInverse11only(const MultiVector &RHS, MultiVector &X) const
apply solve to 1-1 block only
void apply(const MultiVector &X, MultiVector &Y, Teuchos::ETransp mode=Teuchos::NO_TRANS, Scalar alpha=Teuchos::ScalarTraits< Scalar >::one(), Scalar beta=Teuchos::ScalarTraits< Scalar >::zero()) const
#define TEUCHOS_ASSERT(assertion_test)
static void ZeroDirichletCols(Teuchos::RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > &A, const Teuchos::ArrayRCP< const bool > &dirichletCols, Scalar replaceWith=Teuchos::ScalarTraits< Scalar >::zero())
#define TEUCHOS_ASSERT_EQUALITY(val1, val2)
Factory for building coarse matrices.
virtual size_t getNumVectors() const=0
static RCP< Tpetra::RowMatrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Op2NonConstTpetraRow(RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Op)
Factory for building uncoupled aggregates.
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 > &=Teuchos::null)
Factory for building a thresholded operator.
Kokkos::View< const bool *, typename NO::device_type > DetectDirichletRows(const Xpetra::Matrix< SC, LO, GO, NO > &A, const typename Teuchos::ScalarTraits< SC >::magnitudeType &tol, const bool count_twos_as_dirichlet)
void resetMatrix(Teuchos::RCP< Matrix > SM_Matrix_new)
Reset system matrix.
Teuchos::RCP< const Map > getDomainMap() const
Returns the Xpetra::Map object associated with the domain of this operator.
void Request(const FactoryBase &factory)
Increment the storage counter for all the inputs of a factory.