Ifpack2 Templated Preconditioning Package  Version 1.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
Ifpack2_Hypre_def.hpp
1 // @HEADER
2 // *****************************************************************************
3 // Ifpack2: Templated Object-Oriented Algebraic Preconditioner Package
4 //
5 // Copyright 2009 NTESS and the Ifpack2 contributors.
6 // SPDX-License-Identifier: BSD-3-Clause
7 // *****************************************************************************
8 // @HEADER
9 
10 #ifndef IFPACK2_HYPRE_DEF_HPP
11 #define IFPACK2_HYPRE_DEF_HPP
12 
13 #include "Ifpack2_Hypre_decl.hpp"
14 #if defined(HAVE_IFPACK2_HYPRE) && defined(HAVE_IFPACK2_MPI)
15 #include <stdexcept>
16 
17 #include "Tpetra_Import.hpp"
18 #include "Teuchos_ParameterList.hpp"
19 #include "Teuchos_RCP.hpp"
20 #include "Teuchos_DefaultMpiComm.hpp"
21 #include "HYPRE_IJ_mv.h"
22 #include "HYPRE_parcsr_ls.h"
23 #include "krylov.h"
24 #include "_hypre_parcsr_mv.h"
25 #include "_hypre_IJ_mv.h"
26 #include "HYPRE_parcsr_mv.h"
27 #include "HYPRE.h"
28 
29 using Teuchos::RCP;
30 using Teuchos::rcp;
31 using Teuchos::rcpFromRef;
32 
33 namespace Ifpack2 {
34 
35 template <class LocalOrdinal, class Node>
36 Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_Int, Node> >::
37  Hypre(const Teuchos::RCP<const row_matrix_type> &A)
38  : A_(A)
39  , IsInitialized_(false)
40  , IsComputed_(false)
41  , NumInitialize_(0)
42  , NumCompute_(0)
43  , NumApply_(0)
44  , InitializeTime_(0.0)
45  , ComputeTime_(0.0)
46  , ApplyTime_(0.0)
47  , HypreA_(0)
48  , HypreG_(0)
49  , xHypre_(0)
50  , yHypre_(0)
51  , zHypre_(0)
52  , IsSolverCreated_(false)
53  , IsPrecondCreated_(false)
54  , SolveOrPrec_(Hypre_Is_Solver)
55  , NumFunsToCall_(0)
56  , SolverType_(PCG)
57  , PrecondType_(Euclid)
58  , UsePreconditioner_(false)
59  , Dump_(false) {}
60 
61 //==============================================================================
62 template <class LocalOrdinal, class Node>
63 Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_Int, Node> >::~Hypre() {
64  Destroy();
65 }
66 
67 //==============================================================================
68 template <class LocalOrdinal, class Node>
69 void Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_Int, Node> >::Destroy() {
70  if (isInitialized()) {
71  IFPACK2_CHK_ERRV(HYPRE_IJMatrixDestroy(HypreA_));
72  IFPACK2_CHK_ERRV(HYPRE_IJVectorDestroy(XHypre_));
73  IFPACK2_CHK_ERRV(HYPRE_IJVectorDestroy(YHypre_));
74  }
75  if (IsSolverCreated_) {
76  IFPACK2_CHK_ERRV(SolverDestroyPtr_(Solver_));
77  }
78  if (IsPrecondCreated_) {
79  IFPACK2_CHK_ERRV(PrecondDestroyPtr_(Preconditioner_));
80  }
81 
82  // Maxwell
83  if (HypreG_) {
84  IFPACK2_CHK_ERRV(HYPRE_IJMatrixDestroy(HypreG_));
85  }
86  if (xHypre_) {
87  IFPACK2_CHK_ERRV(HYPRE_IJVectorDestroy(xHypre_));
88  }
89  if (yHypre_) {
90  IFPACK2_CHK_ERRV(HYPRE_IJVectorDestroy(yHypre_));
91  }
92  if (zHypre_) {
93  IFPACK2_CHK_ERRV(HYPRE_IJVectorDestroy(zHypre_));
94  }
95 } // Destroy()
96 
97 //==============================================================================
98 template <class LocalOrdinal, class Node>
99 void Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_Int, Node> >::initialize() {
100  const std::string timerName("Ifpack2::Hypre::initialize");
101  Teuchos::RCP<Teuchos::Time> timer = Teuchos::TimeMonitor::lookupCounter(timerName);
102  if (timer.is_null()) timer = Teuchos::TimeMonitor::getNewCounter(timerName);
103 
104  if (IsInitialized_) return;
105  double startTime = timer->wallTime();
106  {
107  Teuchos::TimeMonitor timeMon(*timer);
108 
109  MPI_Comm comm = *(Teuchos::rcp_dynamic_cast<const Teuchos::MpiComm<int> >(A_->getRowMap()->getComm())->getRawMpiComm());
110 
111  // Check that RowMap and RangeMap are the same. While this could handle the
112  // case where RowMap and RangeMap are permutations, other Ifpack PCs don't
113  // handle this either.
114  if (!A_->getRowMap()->isSameAs(*A_->getRangeMap())) {
115  IFPACK2_CHK_ERRV(-1);
116  }
117  // Hypre expects the RowMap to be Linear.
118  if (A_->getRowMap()->isContiguous()) {
119  GloballyContiguousRowMap_ = A_->getRowMap();
120  GloballyContiguousColMap_ = A_->getColMap();
121  } else {
122  // Must create GloballyContiguous Maps for Hypre
123  if (A_->getDomainMap()->isSameAs(*A_->getRowMap())) {
124  Teuchos::RCP<const crs_matrix_type> Aconst = Teuchos::rcp_dynamic_cast<const crs_matrix_type>(A_);
125  GloballyContiguousColMap_ = MakeContiguousColumnMap(Aconst);
126  GloballyContiguousRowMap_ = rcp(new map_type(A_->getRowMap()->getGlobalNumElements(),
127  A_->getRowMap()->getLocalNumElements(), 0, A_->getRowMap()->getComm()));
128  } else {
129  throw std::runtime_error("Ifpack_Hypre: Unsupported map configuration: Row/Domain maps do not match");
130  }
131  }
132  // Next create vectors that will be used when ApplyInverse() is called
133  HYPRE_Int ilower = GloballyContiguousRowMap_->getMinGlobalIndex();
134  HYPRE_Int iupper = GloballyContiguousRowMap_->getMaxGlobalIndex();
135  // X in AX = Y
136  IFPACK2_CHK_ERRV(HYPRE_IJVectorCreate(comm, ilower, iupper, &XHypre_));
137  IFPACK2_CHK_ERRV(HYPRE_IJVectorSetObjectType(XHypre_, HYPRE_PARCSR));
138  IFPACK2_CHK_ERRV(HYPRE_IJVectorInitialize(XHypre_));
139  IFPACK2_CHK_ERRV(HYPRE_IJVectorAssemble(XHypre_));
140  IFPACK2_CHK_ERRV(HYPRE_IJVectorGetObject(XHypre_, (void **)&ParX_));
141  XVec_ = Teuchos::rcp((hypre_ParVector *)hypre_IJVectorObject(((hypre_IJVector *)XHypre_)), false);
142 
143  // Y in AX = Y
144  IFPACK2_CHK_ERRV(HYPRE_IJVectorCreate(comm, ilower, iupper, &YHypre_));
145  IFPACK2_CHK_ERRV(HYPRE_IJVectorSetObjectType(YHypre_, HYPRE_PARCSR));
146  IFPACK2_CHK_ERRV(HYPRE_IJVectorInitialize(YHypre_));
147  IFPACK2_CHK_ERRV(HYPRE_IJVectorAssemble(YHypre_));
148  IFPACK2_CHK_ERRV(HYPRE_IJVectorGetObject(YHypre_, (void **)&ParY_));
149  YVec_ = Teuchos::rcp((hypre_ParVector *)hypre_IJVectorObject(((hypre_IJVector *)YHypre_)), false);
150 
151  // Cache
152  VectorCache_.resize(A_->getRowMap()->getLocalNumElements());
153 
154  // set flags
155  IsInitialized_ = true;
156  NumInitialize_++;
157  }
158  InitializeTime_ += (timer->wallTime() - startTime);
159 } // Initialize()
160 
161 //==============================================================================
162 template <class LocalOrdinal, class Node>
163 void Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_Int, Node> >::setParameters(const Teuchos::ParameterList &list) {
164  std::map<std::string, Hypre_Solver> solverMap;
165  solverMap["BoomerAMG"] = BoomerAMG;
166  solverMap["ParaSails"] = ParaSails;
167  solverMap["Euclid"] = Euclid;
168  solverMap["AMS"] = AMS;
169  solverMap["Hybrid"] = Hybrid;
170  solverMap["PCG"] = PCG;
171  solverMap["GMRES"] = GMRES;
172  solverMap["FlexGMRES"] = FlexGMRES;
173  solverMap["LGMRES"] = LGMRES;
174  solverMap["BiCGSTAB"] = BiCGSTAB;
175 
176  std::map<std::string, Hypre_Chooser> chooserMap;
177  chooserMap["Solver"] = Hypre_Is_Solver;
178  chooserMap["Preconditioner"] = Hypre_Is_Preconditioner;
179 
180  List_ = list;
181  Hypre_Solver solType;
182  if (list.isType<std::string>("hypre: Solver"))
183  solType = solverMap[list.get<std::string>("hypre: Solver")];
184  else if (list.isParameter("hypre: Solver"))
185  solType = (Hypre_Solver)list.get<int>("hypre: Solver");
186  else
187  solType = PCG;
188  SolverType_ = solType;
189  Hypre_Solver precType;
190  if (list.isType<std::string>("hypre: Preconditioner"))
191  precType = solverMap[list.get<std::string>("hypre: Preconditioner")];
192  else if (list.isParameter("hypre: Preconditioner"))
193  precType = (Hypre_Solver)list.get<int>("hypre: Preconditioner");
194  else
195  precType = Euclid;
196  PrecondType_ = precType;
197  Hypre_Chooser chooser;
198  if (list.isType<std::string>("hypre: SolveOrPrecondition"))
199  chooser = chooserMap[list.get<std::string>("hypre: SolveOrPrecondition")];
200  else if (list.isParameter("hypre: SolveOrPrecondition"))
201  chooser = (Hypre_Chooser)list.get<int>("hypre: SolveOrPrecondition");
202  else
203  chooser = Hypre_Is_Solver;
204  SolveOrPrec_ = chooser;
205  bool SetPrecond = list.isParameter("hypre: SetPreconditioner") ? list.get<bool>("hypre: SetPreconditioner") : false;
206  IFPACK2_CHK_ERR(SetParameter(SetPrecond));
207  int NumFunctions = list.isParameter("hypre: NumFunctions") ? list.get<int>("hypre: NumFunctions") : 0;
208  FunsToCall_.clear();
209  NumFunsToCall_ = 0;
210  if (NumFunctions > 0) {
211  RCP<FunctionParameter> *params = list.get<RCP<FunctionParameter> *>("hypre: Functions");
212  for (int i = 0; i < NumFunctions; i++) {
213  IFPACK2_CHK_ERR(AddFunToList(params[i]));
214  }
215  }
216 
217  if (list.isSublist("hypre: Solver functions")) {
218  Teuchos::ParameterList solverList = list.sublist("hypre: Solver functions");
219  for (auto it = solverList.begin(); it != solverList.end(); ++it) {
220  std::string funct_name = it->first;
221  if (it->second.isType<HYPRE_Int>()) {
222  IFPACK2_CHK_ERR(AddFunToList(rcp(new FunctionParameter(Hypre_Is_Solver, funct_name, Teuchos::getValue<HYPRE_Int>(it->second)))));
223  } else if (!std::is_same<HYPRE_Int, int>::value && it->second.isType<int>()) {
224  IFPACK2_CHK_ERR(AddFunToList(rcp(new FunctionParameter(Hypre_Is_Solver, funct_name, Teuchos::as<HYPRE_Int>(Teuchos::getValue<int>(it->second))))));
225  } else if (it->second.isType<HYPRE_Real>()) {
226  IFPACK2_CHK_ERR(AddFunToList(rcp(new FunctionParameter(Hypre_Is_Solver, funct_name, Teuchos::getValue<HYPRE_Real>(it->second)))));
227  } else {
228  IFPACK2_CHK_ERR(-1);
229  }
230  }
231  }
232 
233  if (list.isSublist("hypre: Preconditioner functions")) {
234  Teuchos::ParameterList precList = list.sublist("hypre: Preconditioner functions");
235  for (auto it = precList.begin(); it != precList.end(); ++it) {
236  std::string funct_name = it->first;
237  if (it->second.isType<HYPRE_Int>()) {
238  IFPACK2_CHK_ERR(AddFunToList(rcp(new FunctionParameter(Hypre_Is_Preconditioner, funct_name, Teuchos::getValue<HYPRE_Int>(it->second)))));
239  } else if (!std::is_same<HYPRE_Int, int>::value && it->second.isType<int>()) {
240  IFPACK2_CHK_ERR(AddFunToList(rcp(new FunctionParameter(Hypre_Is_Preconditioner, funct_name, Teuchos::as<HYPRE_Int>(Teuchos::getValue<int>(it->second))))));
241  } else if (it->second.isType<HYPRE_Real>()) {
242  IFPACK2_CHK_ERR(AddFunToList(rcp(new FunctionParameter(Hypre_Is_Preconditioner, funct_name, Teuchos::getValue<HYPRE_Real>(it->second)))));
243  } else if (it->second.isList()) {
244  Teuchos::ParameterList pl = Teuchos::getValue<Teuchos::ParameterList>(it->second);
245  if (FunctionParameter::isFuncIntInt(funct_name)) {
246  HYPRE_Int arg0 = pl.get<HYPRE_Int>("arg 0");
247  HYPRE_Int arg1 = pl.get<HYPRE_Int>("arg 1");
248  IFPACK2_CHK_ERR(AddFunToList(rcp(new FunctionParameter(Hypre_Is_Preconditioner, funct_name, arg0, arg1))));
249  } else if (FunctionParameter::isFuncIntIntDoubleDouble(funct_name)) {
250  HYPRE_Int arg0 = pl.get<HYPRE_Int>("arg 0");
251  HYPRE_Int arg1 = pl.get<HYPRE_Int>("arg 1");
252  HYPRE_Real arg2 = pl.get<HYPRE_Real>("arg 2");
253  HYPRE_Real arg3 = pl.get<HYPRE_Real>("arg 3");
254  IFPACK2_CHK_ERR(AddFunToList(rcp(new FunctionParameter(Hypre_Is_Preconditioner, funct_name, arg0, arg1, arg2, arg3))));
255  } else if (FunctionParameter::isFuncIntIntIntDoubleIntInt(funct_name)) {
256  HYPRE_Int arg0 = pl.get<HYPRE_Int>("arg 0");
257  HYPRE_Int arg1 = pl.get<HYPRE_Int>("arg 1");
258  HYPRE_Int arg2 = pl.get<HYPRE_Int>("arg 2");
259  HYPRE_Real arg3 = pl.get<HYPRE_Real>("arg 3");
260  HYPRE_Int arg4 = pl.get<HYPRE_Int>("arg 4");
261  HYPRE_Int arg5 = pl.get<HYPRE_Int>("arg 5");
262  IFPACK2_CHK_ERR(AddFunToList(rcp(new FunctionParameter(Hypre_Is_Preconditioner, funct_name, arg0, arg1, arg2, arg3, arg4, arg5))));
263  } else {
264  IFPACK2_CHK_ERR(-1);
265  }
266  }
267  }
268  }
269 
270  if (list.isSublist("Coordinates") && list.sublist("Coordinates").isType<Teuchos::RCP<multivector_type> >("Coordinates"))
271  Coords_ = list.sublist("Coordinates").get<Teuchos::RCP<multivector_type> >("Coordinates");
272  if (list.isSublist("Operators") && list.sublist("Operators").isType<Teuchos::RCP<const crs_matrix_type> >("G"))
273  G_ = list.sublist("Operators").get<Teuchos::RCP<const crs_matrix_type> >("G");
274 
275  Dump_ = list.isParameter("hypre: Dump") ? list.get<bool>("hypre: Dump") : false;
276 } // setParameters()
277 
278 //==============================================================================
279 template <class LocalOrdinal, class Node>
280 int Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_Int, Node> >::AddFunToList(RCP<FunctionParameter> NewFun) {
281  NumFunsToCall_ = NumFunsToCall_ + 1;
282  FunsToCall_.resize(NumFunsToCall_);
283  FunsToCall_[NumFunsToCall_ - 1] = NewFun;
284  return 0;
285 } // AddFunToList()
286 
287 //==============================================================================
288 template <class LocalOrdinal, class Node>
289 int Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_Int, Node> >::SetParameter(Hypre_Chooser chooser, HYPRE_Int (*pt2Func)(HYPRE_Solver, HYPRE_Int), HYPRE_Int parameter) {
290  RCP<FunctionParameter> temp = rcp(new FunctionParameter(chooser, pt2Func, parameter));
291  IFPACK2_CHK_ERR(AddFunToList(temp));
292  return 0;
293 } // SetParameter() - int function pointer
294 
295 //=============================================================================
296 template <class LocalOrdinal, class Node>
297 int Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_Int, Node> >::SetParameter(Hypre_Chooser chooser, HYPRE_Int (*pt2Func)(HYPRE_Solver, HYPRE_Real), HYPRE_Real parameter) {
298  RCP<FunctionParameter> temp = rcp(new FunctionParameter(chooser, pt2Func, parameter));
299  IFPACK2_CHK_ERR(AddFunToList(temp));
300  return 0;
301 } // SetParameter() - double function pointer
302 
303 //==============================================================================
304 template <class LocalOrdinal, class Node>
305 int Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_Int, Node> >::SetParameter(Hypre_Chooser chooser, HYPRE_Int (*pt2Func)(HYPRE_Solver, HYPRE_Real, HYPRE_Int), HYPRE_Real parameter1, HYPRE_Int parameter2) {
306  RCP<FunctionParameter> temp = rcp(new FunctionParameter(chooser, pt2Func, parameter1, parameter2));
307  IFPACK2_CHK_ERR(AddFunToList(temp));
308  return 0;
309 } // SetParameter() - double,int function pointer
310 
311 //==============================================================================
312 template <class LocalOrdinal, class Node>
313 int Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_Int, Node> >::SetParameter(Hypre_Chooser chooser, HYPRE_Int (*pt2Func)(HYPRE_Solver, HYPRE_Int, HYPRE_Real), HYPRE_Int parameter1, HYPRE_Real parameter2) {
314  RCP<FunctionParameter> temp = rcp(new FunctionParameter(chooser, pt2Func, parameter1, parameter2));
315  IFPACK2_CHK_ERR(AddFunToList(temp));
316  return 0;
317 } // SetParameter() - int,double function pointer
318 
319 //==============================================================================
320 template <class LocalOrdinal, class Node>
321 int Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_Int, Node> >::SetParameter(Hypre_Chooser chooser, HYPRE_Int (*pt2Func)(HYPRE_Solver, HYPRE_Int, HYPRE_Int), HYPRE_Int parameter1, HYPRE_Int parameter2) {
322  RCP<FunctionParameter> temp = rcp(new FunctionParameter(chooser, pt2Func, parameter1, parameter2));
323  IFPACK2_CHK_ERR(AddFunToList(temp));
324  return 0;
325 } // SetParameter() int,int function pointer
326 
327 //==============================================================================
328 template <class LocalOrdinal, class Node>
329 int Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_Int, Node> >::SetParameter(Hypre_Chooser chooser, HYPRE_Int (*pt2Func)(HYPRE_Solver, HYPRE_Real *), HYPRE_Real *parameter) {
330  RCP<FunctionParameter> temp = rcp(new FunctionParameter(chooser, pt2Func, parameter));
331  IFPACK2_CHK_ERR(AddFunToList(temp));
332  return 0;
333 } // SetParameter() - double* function pointer
334 
335 //==============================================================================
336 template <class LocalOrdinal, class Node>
337 int Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_Int, Node> >::SetParameter(Hypre_Chooser chooser, HYPRE_Int (*pt2Func)(HYPRE_Solver, HYPRE_Int *), HYPRE_Int *parameter) {
338  RCP<FunctionParameter> temp = rcp(new FunctionParameter(chooser, pt2Func, parameter));
339  IFPACK2_CHK_ERR(AddFunToList(temp));
340  return 0;
341 } // SetParameter() - int* function pointer
342 
343 //==============================================================================
344 template <class LocalOrdinal, class Node>
345 int Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_Int, Node> >::SetParameter(Hypre_Chooser chooser, HYPRE_Int (*pt2Func)(HYPRE_Solver, HYPRE_Int **), HYPRE_Int **parameter) {
346  RCP<FunctionParameter> temp = rcp(new FunctionParameter(chooser, pt2Func, parameter));
347  IFPACK2_CHK_ERR(AddFunToList(temp));
348  return 0;
349 } // SetParameter() - int** function pointer
350 
351 //==============================================================================
352 template <class LocalOrdinal, class Node>
353 int Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_Int, Node> >::SetParameter(Hypre_Chooser chooser, Hypre_Solver solver) {
354  if (chooser == Hypre_Is_Solver) {
355  SolverType_ = solver;
356  } else {
357  PrecondType_ = solver;
358  }
359  return 0;
360 } // SetParameter() - set type of solver
361 
362 //==============================================================================
363 template <class LocalOrdinal, class Node>
364 int Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_Int, Node> >::SetDiscreteGradient(Teuchos::RCP<const crs_matrix_type> G) {
365  using LO = local_ordinal_type;
366  using GO = global_ordinal_type;
367  // using SC = scalar_type;
368 
369  // Sanity check
370  if (!A_->getRowMap()->isSameAs(*G->getRowMap()))
371  throw std::runtime_error("Hypre<Tpetra::RowMatrix<double, HYPRE_Int, long long, Node>: Edge map mismatch: A and discrete gradient");
372 
373  // Get the maps for the nodes (assuming the edge map from A is OK);
374  GloballyContiguousNodeRowMap_ = rcp(new map_type(G->getDomainMap()->getGlobalNumElements(),
375  G->getDomainMap()->getLocalNumElements(), 0, A_->getRowMap()->getComm()));
376  GloballyContiguousNodeColMap_ = MakeContiguousColumnMap(G);
377 
378  // Start building G
379  MPI_Comm comm = *(Teuchos::rcp_dynamic_cast<const Teuchos::MpiComm<int> >(A_->getRowMap()->getComm())->getRawMpiComm());
380  GO ilower = GloballyContiguousRowMap_->getMinGlobalIndex();
381  GO iupper = GloballyContiguousRowMap_->getMaxGlobalIndex();
382  GO jlower = GloballyContiguousNodeRowMap_->getMinGlobalIndex();
383  GO jupper = GloballyContiguousNodeRowMap_->getMaxGlobalIndex();
384  IFPACK2_CHK_ERR(HYPRE_IJMatrixCreate(comm, ilower, iupper, jlower, jupper, &HypreG_));
385  IFPACK2_CHK_ERR(HYPRE_IJMatrixSetObjectType(HypreG_, HYPRE_PARCSR));
386  IFPACK2_CHK_ERR(HYPRE_IJMatrixInitialize(HypreG_));
387 
388  std::vector<GO> new_indices(G->getLocalMaxNumRowEntries());
389  for (LO i = 0; i < (LO)G->getLocalNumRows(); i++) {
390  typename crs_matrix_type::values_host_view_type values;
391  typename crs_matrix_type::local_inds_host_view_type indices;
392  G->getLocalRowView(i, indices, values);
393  for (LO j = 0; j < (LO)indices.extent(0); j++) {
394  new_indices[j] = GloballyContiguousNodeColMap_->getGlobalElement(indices(j));
395  }
396  GO GlobalRow[1];
397  GO numEntries = (GO)indices.extent(0);
398  GlobalRow[0] = GloballyContiguousRowMap_->getGlobalElement(i);
399  IFPACK2_CHK_ERR(HYPRE_IJMatrixSetValues(HypreG_, 1, &numEntries, GlobalRow, new_indices.data(), values.data()));
400  }
401  IFPACK2_CHK_ERR(HYPRE_IJMatrixAssemble(HypreG_));
402  IFPACK2_CHK_ERR(HYPRE_IJMatrixGetObject(HypreG_, (void **)&ParMatrixG_));
403 
404  if (Dump_)
405  HYPRE_ParCSRMatrixPrint(ParMatrixG_, "G.mat");
406 
407  if (SolverType_ == AMS)
408  HYPRE_AMSSetDiscreteGradient(Solver_, ParMatrixG_);
409  if (PrecondType_ == AMS)
410  HYPRE_AMSSetDiscreteGradient(Preconditioner_, ParMatrixG_);
411  return 0;
412 } // SetDiscreteGradient()
413 
414 //==============================================================================
415 template <class LocalOrdinal, class Node>
416 int Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_Int, Node> >::SetCoordinates(Teuchos::RCP<multivector_type> coords) {
417  if (!G_.is_null() && !G_->getDomainMap()->isSameAs(*coords->getMap()))
418  throw std::runtime_error("Hypre<Tpetra::RowMatrix<double, HYPRE_Int, long long, Node>: Node map mismatch: G->DomainMap() and coords");
419 
420  if (SolverType_ != AMS && PrecondType_ != AMS)
421  return 0;
422 
423  scalar_type *xPtr = coords->getDataNonConst(0).getRawPtr();
424  scalar_type *yPtr = coords->getDataNonConst(1).getRawPtr();
425  scalar_type *zPtr = coords->getDataNonConst(2).getRawPtr();
426 
427  MPI_Comm comm = *(Teuchos::rcp_dynamic_cast<const Teuchos::MpiComm<int> >(A_->getRowMap()->getComm())->getRawMpiComm());
428  local_ordinal_type NumEntries = coords->getLocalLength();
429  global_ordinal_type *indices = const_cast<global_ordinal_type *>(GloballyContiguousNodeRowMap_->getLocalElementList().getRawPtr());
430 
431  global_ordinal_type ilower = GloballyContiguousNodeRowMap_->getMinGlobalIndex();
432  global_ordinal_type iupper = GloballyContiguousNodeRowMap_->getMaxGlobalIndex();
433 
434  if (NumEntries != iupper - ilower + 1) {
435  std::cout << "Ifpack2::Hypre::SetCoordinates(): Error on rank " << A_->getRowMap()->getComm()->getRank() << ": MyLength = " << coords->getLocalLength() << " GID range = [" << ilower << "," << iupper << "]" << std::endl;
436  throw std::runtime_error("Hypre<Tpetra::RowMatrix<double, HYPRE_Int, long long, Node>: SetCoordinates: Length mismatch");
437  }
438 
439  IFPACK2_CHK_ERR(HYPRE_IJVectorCreate(comm, ilower, iupper, &xHypre_));
440  IFPACK2_CHK_ERR(HYPRE_IJVectorSetObjectType(xHypre_, HYPRE_PARCSR));
441  IFPACK2_CHK_ERR(HYPRE_IJVectorInitialize(xHypre_));
442  IFPACK2_CHK_ERR(HYPRE_IJVectorSetValues(xHypre_, NumEntries, indices, xPtr));
443  IFPACK2_CHK_ERR(HYPRE_IJVectorAssemble(xHypre_));
444  IFPACK2_CHK_ERR(HYPRE_IJVectorGetObject(xHypre_, (void **)&xPar_));
445 
446  IFPACK2_CHK_ERR(HYPRE_IJVectorCreate(comm, ilower, iupper, &yHypre_));
447  IFPACK2_CHK_ERR(HYPRE_IJVectorSetObjectType(yHypre_, HYPRE_PARCSR));
448  IFPACK2_CHK_ERR(HYPRE_IJVectorInitialize(yHypre_));
449  IFPACK2_CHK_ERR(HYPRE_IJVectorSetValues(yHypre_, NumEntries, indices, yPtr));
450  IFPACK2_CHK_ERR(HYPRE_IJVectorAssemble(yHypre_));
451  IFPACK2_CHK_ERR(HYPRE_IJVectorGetObject(yHypre_, (void **)&yPar_));
452 
453  IFPACK2_CHK_ERR(HYPRE_IJVectorCreate(comm, ilower, iupper, &zHypre_));
454  IFPACK2_CHK_ERR(HYPRE_IJVectorSetObjectType(zHypre_, HYPRE_PARCSR));
455  IFPACK2_CHK_ERR(HYPRE_IJVectorInitialize(zHypre_));
456  IFPACK2_CHK_ERR(HYPRE_IJVectorSetValues(zHypre_, NumEntries, indices, zPtr));
457  IFPACK2_CHK_ERR(HYPRE_IJVectorAssemble(zHypre_));
458  IFPACK2_CHK_ERR(HYPRE_IJVectorGetObject(zHypre_, (void **)&zPar_));
459 
460  if (Dump_) {
461  HYPRE_ParVectorPrint(xPar_, "coordX.dat");
462  HYPRE_ParVectorPrint(yPar_, "coordY.dat");
463  HYPRE_ParVectorPrint(zPar_, "coordZ.dat");
464  }
465 
466  if (SolverType_ == AMS)
467  HYPRE_AMSSetCoordinateVectors(Solver_, xPar_, yPar_, zPar_);
468  if (PrecondType_ == AMS)
469  HYPRE_AMSSetCoordinateVectors(Preconditioner_, xPar_, yPar_, zPar_);
470 
471  return 0;
472 
473 } // SetCoordinates
474 
475 //==============================================================================
476 template <class LocalOrdinal, class Node>
477 void Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_Int, Node> >::compute() {
478  const std::string timerName("Ifpack2::Hypre::compute");
479  Teuchos::RCP<Teuchos::Time> timer = Teuchos::TimeMonitor::lookupCounter(timerName);
480  if (timer.is_null()) timer = Teuchos::TimeMonitor::getNewCounter(timerName);
481  double startTime = timer->wallTime();
482  // Start timing here.
483  {
484  Teuchos::TimeMonitor timeMon(*timer);
485 
486  if (isInitialized() == false) {
487  initialize();
488  }
489 
490  // Create the Hypre matrix and copy values. Note this uses values (which
491  // Initialize() shouldn't do) but it doesn't care what they are (for
492  // instance they can be uninitialized data even). It should be possible to
493  // set the Hypre structure without copying values, but this is the easiest
494  // way to get the structure.
495  MPI_Comm comm = *(Teuchos::rcp_dynamic_cast<const Teuchos::MpiComm<int> >(A_->getRowMap()->getComm())->getRawMpiComm());
496  global_ordinal_type ilower = GloballyContiguousRowMap_->getMinGlobalIndex();
497  global_ordinal_type iupper = GloballyContiguousRowMap_->getMaxGlobalIndex();
498  IFPACK2_CHK_ERR(HYPRE_IJMatrixCreate(comm, ilower, iupper, ilower, iupper, &HypreA_));
499  IFPACK2_CHK_ERR(HYPRE_IJMatrixSetObjectType(HypreA_, HYPRE_PARCSR));
500  IFPACK2_CHK_ERR(HYPRE_IJMatrixInitialize(HypreA_));
501  CopyTpetraToHypre();
502  if (SolveOrPrec_ == Hypre_Is_Solver) {
503  IFPACK2_CHK_ERR(SetSolverType(SolverType_));
504  if (SolverPrecondPtr_ != NULL && UsePreconditioner_) {
505  // both method allows a PC (first condition) and the user wants a PC (second)
506  IFPACK2_CHK_ERR(SetPrecondType(PrecondType_));
507  CallFunctions();
508  IFPACK2_CHK_ERR(SolverPrecondPtr_(Solver_, PrecondSolvePtr_, PrecondSetupPtr_, Preconditioner_));
509  } else {
510  CallFunctions();
511  }
512  } else {
513  IFPACK2_CHK_ERR(SetPrecondType(PrecondType_));
514  CallFunctions();
515  }
516 
517  if (!G_.is_null()) {
518  SetDiscreteGradient(G_);
519  }
520 
521  if (!Coords_.is_null()) {
522  SetCoordinates(Coords_);
523  }
524 
525  // Hypre Setup must be called after matrix has values
526  if (SolveOrPrec_ == Hypre_Is_Solver) {
527  IFPACK2_CHK_ERR(SolverSetupPtr_(Solver_, ParMatrix_, ParX_, ParY_));
528  } else {
529  IFPACK2_CHK_ERR(PrecondSetupPtr_(Preconditioner_, ParMatrix_, ParX_, ParY_));
530  }
531 
532  IsComputed_ = true;
533  NumCompute_++;
534  }
535 
536  ComputeTime_ += (timer->wallTime() - startTime);
537 } // Compute()
538 
539 //==============================================================================
540 template <class LocalOrdinal, class Node>
541 int Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_Int, Node> >::CallFunctions() const {
542  for (int i = 0; i < NumFunsToCall_; i++) {
543  IFPACK2_CHK_ERR(FunsToCall_[i]->CallFunction(Solver_, Preconditioner_));
544  }
545  return 0;
546 } // CallFunctions()
547 
548 //==============================================================================
549 template <class LocalOrdinal, class Node>
550 void Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_Int, Node> >::apply(const Tpetra::MultiVector<scalar_type, local_ordinal_type, global_ordinal_type, node_type> &X,
551  Tpetra::MultiVector<scalar_type, local_ordinal_type, global_ordinal_type, node_type> &Y,
552  Teuchos::ETransp mode,
553  scalar_type alpha,
554  scalar_type beta) const {
555  using LO = local_ordinal_type;
556  using SC = scalar_type;
557  const std::string timerName("Ifpack2::Hypre::apply");
558  Teuchos::RCP<Teuchos::Time> timer = Teuchos::TimeMonitor::lookupCounter(timerName);
559  if (timer.is_null()) timer = Teuchos::TimeMonitor::getNewCounter(timerName);
560  double startTime = timer->wallTime();
561  // Start timing here.
562  {
563  Teuchos::TimeMonitor timeMon(*timer);
564 
565  if (isComputed() == false) {
566  IFPACK2_CHK_ERR(-1);
567  }
568  hypre_Vector *XLocal_ = hypre_ParVectorLocalVector(XVec_);
569  hypre_Vector *YLocal_ = hypre_ParVectorLocalVector(YVec_);
570  bool SameVectors = false;
571  size_t NumVectors = X.getNumVectors();
572  if (NumVectors != Y.getNumVectors()) IFPACK2_CHK_ERR(-1); // X and Y must have same number of vectors
573  if (&X == &Y) { // FIXME: Maybe not the right way to check this
574  SameVectors = true;
575  }
576 
577  // NOTE: Here were assuming that the local ordering of Epetra's X/Y-vectors and
578  // Hypre's X/Y-vectors are the same. Seeing as as this is more or less how we constructed
579  // the Hypre matrices, this seems pretty reasoanble.
580 
581  for (int VecNum = 0; VecNum < (int)NumVectors; VecNum++) {
582  // Get values for current vector in multivector.
583  // FIXME amk Nov 23, 2015: This will not work for funky data layouts
584  SC *XValues = const_cast<SC *>(X.getData(VecNum).getRawPtr());
585  SC *YValues;
586  if (!SameVectors) {
587  YValues = const_cast<SC *>(Y.getData(VecNum).getRawPtr());
588  } else {
589  YValues = VectorCache_.getRawPtr();
590  }
591  // Temporarily make a pointer to data in Hypre for end
592  SC *XTemp = XLocal_->data;
593  // Replace data in Hypre vectors with Epetra data
594  XLocal_->data = XValues;
595  SC *YTemp = YLocal_->data;
596  YLocal_->data = YValues;
597 
598  IFPACK2_CHK_ERR(HYPRE_ParVectorSetConstantValues(ParY_, 0.0));
599  if (SolveOrPrec_ == Hypre_Is_Solver) {
600  // Use the solver methods
601  IFPACK2_CHK_ERR(SolverSolvePtr_(Solver_, ParMatrix_, ParX_, ParY_));
602  } else {
603  // Apply the preconditioner
604  IFPACK2_CHK_ERR(PrecondSolvePtr_(Preconditioner_, ParMatrix_, ParX_, ParY_));
605  }
606 
607  if (SameVectors) {
608  Teuchos::ArrayView<SC> Yv = Y.getDataNonConst(VecNum)();
609  LO NumEntries = Y.getLocalLength();
610  for (LO i = 0; i < NumEntries; i++)
611  Yv[i] = YValues[i];
612  }
613  XLocal_->data = XTemp;
614  YLocal_->data = YTemp;
615  }
616  NumApply_++;
617  }
618  ApplyTime_ += (timer->wallTime() - startTime);
619 } // apply()
620 
621 //==============================================================================
622 template <class LocalOrdinal, class Node>
623 void Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_Int, Node> >::applyMat(const Tpetra::MultiVector<scalar_type, local_ordinal_type, global_ordinal_type, node_type> &X,
624  Tpetra::MultiVector<scalar_type, local_ordinal_type, global_ordinal_type, node_type> &Y,
625  Teuchos::ETransp mode) const {
626  A_->apply(X, Y, mode);
627 } // applyMat()
628 
629 //==============================================================================
630 template <class LocalOrdinal, class Node>
631 std::string Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_Int, Node> >::description() const {
632  std::ostringstream out;
633 
634  // Output is a valid YAML dictionary in flow style. If you don't
635  // like everything on a single line, you should call describe()
636  // instead.
637  out << "\"Ifpack2::Hypre\": {";
638  out << "Initialized: " << (isInitialized() ? "true" : "false") << ", "
639  << "Computed: " << (isComputed() ? "true" : "false") << ", ";
640 
641  if (A_.is_null()) {
642  out << "Matrix: null";
643  } else {
644  out << "Global matrix dimensions: ["
645  << A_->getGlobalNumRows() << ", "
646  << A_->getGlobalNumCols() << "]"
647  << ", Global nnz: " << A_->getGlobalNumEntries();
648  }
649 
650  out << "}";
651  return out.str();
652 } // description()
653 
654 //==============================================================================
655 template <class LocalOrdinal, class Node>
656 void Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_Int, Node> >::describe(Teuchos::FancyOStream &os, const Teuchos::EVerbosityLevel verbLevel) const {
657  using std::endl;
658  os << endl;
659  os << "================================================================================" << endl;
660  os << "Ifpack2::Hypre: " << endl
661  << endl;
662  os << "Using " << A_->getComm()->getSize() << " processors." << endl;
663  os << "Global number of rows = " << A_->getGlobalNumRows() << endl;
664  os << "Global number of nonzeros = " << A_->getGlobalNumEntries() << endl;
665  // os << "Condition number estimate = " << Condest() << endl;
666  os << endl;
667  os << "Phase # calls Total Time (s)" << endl;
668  os << "----- ------- --------------" << endl;
669  os << "Initialize() " << std::setw(5) << NumInitialize_
670  << " " << std::setw(15) << InitializeTime_ << endl;
671  os << "Compute() " << std::setw(5) << NumCompute_
672  << " " << std::setw(15) << ComputeTime_ << endl;
673  os << "ApplyInverse() " << std::setw(5) << NumApply_
674  << " " << std::setw(15) << ApplyTime_ << endl;
675  os << "================================================================================" << endl;
676  os << endl;
677 } // description
678 
679 //==============================================================================
680 template <class LocalOrdinal, class Node>
681 int Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_Int, Node> >::SetSolverType(Hypre_Solver Solver) {
682  switch (Solver) {
683  case BoomerAMG:
684  if (IsSolverCreated_) {
685  SolverDestroyPtr_(Solver_);
686  IsSolverCreated_ = false;
687  }
688  SolverCreatePtr_ = &Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_Int, Node> >::Hypre_BoomerAMGCreate;
689  SolverDestroyPtr_ = &HYPRE_BoomerAMGDestroy;
690  SolverSetupPtr_ = &HYPRE_BoomerAMGSetup;
691  SolverPrecondPtr_ = NULL;
692  SolverSolvePtr_ = &HYPRE_BoomerAMGSolve;
693  break;
694  case AMS:
695  if (IsSolverCreated_) {
696  SolverDestroyPtr_(Solver_);
697  IsSolverCreated_ = false;
698  }
699  SolverCreatePtr_ = &Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_Int, Node> >::Hypre_AMSCreate;
700  SolverDestroyPtr_ = &HYPRE_AMSDestroy;
701  SolverSetupPtr_ = &HYPRE_AMSSetup;
702  SolverSolvePtr_ = &HYPRE_AMSSolve;
703  SolverPrecondPtr_ = NULL;
704  break;
705  case Hybrid:
706  if (IsSolverCreated_) {
707  SolverDestroyPtr_(Solver_);
708  IsSolverCreated_ = false;
709  }
710  SolverCreatePtr_ = &Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_Int, Node> >::Hypre_ParCSRHybridCreate;
711  SolverDestroyPtr_ = &HYPRE_ParCSRHybridDestroy;
712  SolverSetupPtr_ = &HYPRE_ParCSRHybridSetup;
713  SolverSolvePtr_ = &HYPRE_ParCSRHybridSolve;
714  SolverPrecondPtr_ = &HYPRE_ParCSRHybridSetPrecond;
715  break;
716  case PCG:
717  if (IsSolverCreated_) {
718  SolverDestroyPtr_(Solver_);
719  IsSolverCreated_ = false;
720  }
721  SolverCreatePtr_ = &Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_Int, Node> >::Hypre_ParCSRPCGCreate;
722  SolverDestroyPtr_ = &HYPRE_ParCSRPCGDestroy;
723  SolverSetupPtr_ = &HYPRE_ParCSRPCGSetup;
724  SolverSolvePtr_ = &HYPRE_ParCSRPCGSolve;
725  SolverPrecondPtr_ = &HYPRE_ParCSRPCGSetPrecond;
726  break;
727  case GMRES:
728  if (IsSolverCreated_) {
729  SolverDestroyPtr_(Solver_);
730  IsSolverCreated_ = false;
731  }
732  SolverCreatePtr_ = &Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_Int, Node> >::Hypre_ParCSRGMRESCreate;
733  SolverDestroyPtr_ = &HYPRE_ParCSRGMRESDestroy;
734  SolverSetupPtr_ = &HYPRE_ParCSRGMRESSetup;
735  SolverPrecondPtr_ = &HYPRE_ParCSRGMRESSetPrecond;
736  break;
737  case FlexGMRES:
738  if (IsSolverCreated_) {
739  SolverDestroyPtr_(Solver_);
740  IsSolverCreated_ = false;
741  }
742  SolverCreatePtr_ = &Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_Int, Node> >::Hypre_ParCSRFlexGMRESCreate;
743  SolverDestroyPtr_ = &HYPRE_ParCSRFlexGMRESDestroy;
744  SolverSetupPtr_ = &HYPRE_ParCSRFlexGMRESSetup;
745  SolverSolvePtr_ = &HYPRE_ParCSRFlexGMRESSolve;
746  SolverPrecondPtr_ = &HYPRE_ParCSRFlexGMRESSetPrecond;
747  break;
748  case LGMRES:
749  if (IsSolverCreated_) {
750  SolverDestroyPtr_(Solver_);
751  IsSolverCreated_ = false;
752  }
753  SolverCreatePtr_ = &Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_Int, Node> >::Hypre_ParCSRLGMRESCreate;
754  SolverDestroyPtr_ = &HYPRE_ParCSRLGMRESDestroy;
755  SolverSetupPtr_ = &HYPRE_ParCSRLGMRESSetup;
756  SolverSolvePtr_ = &HYPRE_ParCSRLGMRESSolve;
757  SolverPrecondPtr_ = &HYPRE_ParCSRLGMRESSetPrecond;
758  break;
759  case BiCGSTAB:
760  if (IsSolverCreated_) {
761  SolverDestroyPtr_(Solver_);
762  IsSolverCreated_ = false;
763  }
764  SolverCreatePtr_ = &Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_Int, Node> >::Hypre_ParCSRBiCGSTABCreate;
765  SolverDestroyPtr_ = &HYPRE_ParCSRBiCGSTABDestroy;
766  SolverSetupPtr_ = &HYPRE_ParCSRBiCGSTABSetup;
767  SolverSolvePtr_ = &HYPRE_ParCSRBiCGSTABSolve;
768  SolverPrecondPtr_ = &HYPRE_ParCSRBiCGSTABSetPrecond;
769  break;
770  default:
771  return -1;
772  }
773  CreateSolver();
774  return 0;
775 } // SetSolverType()
776 
777 //==============================================================================
778 template <class LocalOrdinal, class Node>
779 int Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_Int, Node> >::SetPrecondType(Hypre_Solver Precond) {
780  switch (Precond) {
781  case BoomerAMG:
782  if (IsPrecondCreated_) {
783  PrecondDestroyPtr_(Preconditioner_);
784  IsPrecondCreated_ = false;
785  }
786  PrecondCreatePtr_ = &Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_Int, Node> >::Hypre_BoomerAMGCreate;
787  PrecondDestroyPtr_ = &HYPRE_BoomerAMGDestroy;
788  PrecondSetupPtr_ = &HYPRE_BoomerAMGSetup;
789  PrecondSolvePtr_ = &HYPRE_BoomerAMGSolve;
790  break;
791  case ParaSails:
792  if (IsPrecondCreated_) {
793  PrecondDestroyPtr_(Preconditioner_);
794  IsPrecondCreated_ = false;
795  }
796  PrecondCreatePtr_ = &Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_Int, Node> >::Hypre_ParaSailsCreate;
797  PrecondDestroyPtr_ = &HYPRE_ParaSailsDestroy;
798  PrecondSetupPtr_ = &HYPRE_ParaSailsSetup;
799  PrecondSolvePtr_ = &HYPRE_ParaSailsSolve;
800  break;
801  case Euclid:
802  if (IsPrecondCreated_) {
803  PrecondDestroyPtr_(Preconditioner_);
804  IsPrecondCreated_ = false;
805  }
806  PrecondCreatePtr_ = &Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_Int, Node> >::Hypre_EuclidCreate;
807  PrecondDestroyPtr_ = &HYPRE_EuclidDestroy;
808  PrecondSetupPtr_ = &HYPRE_EuclidSetup;
809  PrecondSolvePtr_ = &HYPRE_EuclidSolve;
810  break;
811  case AMS:
812  if (IsPrecondCreated_) {
813  PrecondDestroyPtr_(Preconditioner_);
814  IsPrecondCreated_ = false;
815  }
816  PrecondCreatePtr_ = &Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_Int, Node> >::Hypre_AMSCreate;
817  PrecondDestroyPtr_ = &HYPRE_AMSDestroy;
818  PrecondSetupPtr_ = &HYPRE_AMSSetup;
819  PrecondSolvePtr_ = &HYPRE_AMSSolve;
820  break;
821  default:
822  return -1;
823  }
824  CreatePrecond();
825  return 0;
826 
827 } // SetPrecondType()
828 
829 //==============================================================================
830 template <class LocalOrdinal, class Node>
831 int Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_Int, Node> >::CreateSolver() {
832  MPI_Comm comm;
833  HYPRE_ParCSRMatrixGetComm(ParMatrix_, &comm);
834  int ierr = (this->*SolverCreatePtr_)(comm, &Solver_);
835  IsSolverCreated_ = true;
836  return ierr;
837 } // CreateSolver()
838 
839 //==============================================================================
840 template <class LocalOrdinal, class Node>
841 int Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_Int, Node> >::CreatePrecond() {
842  MPI_Comm comm;
843  HYPRE_ParCSRMatrixGetComm(ParMatrix_, &comm);
844  int ierr = (this->*PrecondCreatePtr_)(comm, &Preconditioner_);
845  IsPrecondCreated_ = true;
846  return ierr;
847 } // CreatePrecond()
848 
849 //==============================================================================
850 template <class LocalOrdinal, class Node>
851 int Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_Int, Node> >::CopyTpetraToHypre() {
852  using LO = local_ordinal_type;
853  using GO = global_ordinal_type;
854  // using SC = scalar_type;
855 
856  Teuchos::RCP<const crs_matrix_type> Matrix = Teuchos::rcp_dynamic_cast<const crs_matrix_type>(A_);
857  if (Matrix.is_null())
858  throw std::runtime_error("Hypre<Tpetra::RowMatrix<double, LocalOrdinal, HYPRE_Int, Node>: Unsupported matrix configuration: Tpetra::CrsMatrix required");
859 
860  std::vector<HYPRE_Int> new_indices(Matrix->getLocalMaxNumRowEntries());
861  for (LO i = 0; i < (LO)Matrix->getLocalNumRows(); i++) {
862  typename crs_matrix_type::values_host_view_type values;
863  typename crs_matrix_type::local_inds_host_view_type indices;
864  Matrix->getLocalRowView(i, indices, values);
865  for (LO j = 0; j < (LO)indices.extent(0); j++) {
866  new_indices[j] = GloballyContiguousColMap_->getGlobalElement(indices(j));
867  }
868  HYPRE_Int GlobalRow[1];
869  HYPRE_Int numEntries = (GO)indices.extent(0);
870  GlobalRow[0] = GloballyContiguousRowMap_->getGlobalElement(i);
871  IFPACK2_CHK_ERR(HYPRE_IJMatrixSetValues(HypreA_, 1, &numEntries, GlobalRow, new_indices.data(), values.data()));
872  }
873  IFPACK2_CHK_ERR(HYPRE_IJMatrixAssemble(HypreA_));
874  IFPACK2_CHK_ERR(HYPRE_IJMatrixGetObject(HypreA_, (void **)&ParMatrix_));
875  if (Dump_)
876  HYPRE_ParCSRMatrixPrint(ParMatrix_, "A.mat");
877  return 0;
878 } // CopyTpetraToHypre()
879 
880 //==============================================================================
881 template <class LocalOrdinal, class Node>
882 HYPRE_Int Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_Int, Node> >::Hypre_BoomerAMGCreate(MPI_Comm /*comm*/, HYPRE_Solver *solver) { return HYPRE_BoomerAMGCreate(solver); }
883 
884 //==============================================================================
885 template <class LocalOrdinal, class Node>
886 HYPRE_Int Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_Int, Node> >::Hypre_ParaSailsCreate(MPI_Comm comm, HYPRE_Solver *solver) { return HYPRE_ParaSailsCreate(comm, solver); }
887 
888 //==============================================================================
889 template <class LocalOrdinal, class Node>
890 HYPRE_Int Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_Int, Node> >::Hypre_EuclidCreate(MPI_Comm comm, HYPRE_Solver *solver) { return HYPRE_EuclidCreate(comm, solver); }
891 
892 //==============================================================================
893 template <class LocalOrdinal, class Node>
894 HYPRE_Int Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_Int, Node> >::Hypre_AMSCreate(MPI_Comm /*comm*/, HYPRE_Solver *solver) { return HYPRE_AMSCreate(solver); }
895 
896 //==============================================================================
897 template <class LocalOrdinal, class Node>
898 HYPRE_Int Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_Int, Node> >::Hypre_ParCSRHybridCreate(MPI_Comm /*comm*/, HYPRE_Solver *solver) { return HYPRE_ParCSRHybridCreate(solver); }
899 
900 //==============================================================================
901 template <class LocalOrdinal, class Node>
902 HYPRE_Int Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_Int, Node> >::Hypre_ParCSRPCGCreate(MPI_Comm comm, HYPRE_Solver *solver) { return HYPRE_ParCSRPCGCreate(comm, solver); }
903 
904 //==============================================================================
905 template <class LocalOrdinal, class Node>
906 HYPRE_Int Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_Int, Node> >::Hypre_ParCSRGMRESCreate(MPI_Comm comm, HYPRE_Solver *solver) { return HYPRE_ParCSRGMRESCreate(comm, solver); }
907 
908 //==============================================================================
909 template <class LocalOrdinal, class Node>
910 HYPRE_Int Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_Int, Node> >::Hypre_ParCSRFlexGMRESCreate(MPI_Comm comm, HYPRE_Solver *solver) { return HYPRE_ParCSRFlexGMRESCreate(comm, solver); }
911 
912 //==============================================================================
913 template <class LocalOrdinal, class Node>
914 HYPRE_Int Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_Int, Node> >::Hypre_ParCSRLGMRESCreate(MPI_Comm comm, HYPRE_Solver *solver) { return HYPRE_ParCSRLGMRESCreate(comm, solver); }
915 
916 //==============================================================================
917 template <class LocalOrdinal, class Node>
918 HYPRE_Int Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_Int, Node> >::Hypre_ParCSRBiCGSTABCreate(MPI_Comm comm, HYPRE_Solver *solver) { return HYPRE_ParCSRBiCGSTABCreate(comm, solver); }
919 
920 //==============================================================================
921 template <class LocalOrdinal, class Node>
922 Teuchos::RCP<const Tpetra::Map<LocalOrdinal, HYPRE_Int, Node> >
923 Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_Int, Node> >::MakeContiguousColumnMap(Teuchos::RCP<const crs_matrix_type> &Matrix) const {
924  using import_type = Tpetra::Import<local_ordinal_type, global_ordinal_type, node_type>;
925  using go_vector_type = Tpetra::Vector<global_ordinal_type, local_ordinal_type, global_ordinal_type, node_type>;
926 
927  // Must create GloballyContiguous DomainMap (which is a permutation of Matrix_'s
928  // DomainMap) and the corresponding permuted ColumnMap.
929  // Epetra_GID ---------> LID ----------> HYPRE_GID
930  // via DomainMap.LID() via GloballyContiguousDomainMap.GID()
931  if (Matrix.is_null())
932  throw std::runtime_error("Hypre<Tpetra::RowMatrix<HYPRE_Real, HYPRE_Int, long long, Node>: Unsupported matrix configuration: Tpetra::CrsMatrix required");
933  RCP<const map_type> DomainMap = Matrix->getDomainMap();
934  RCP<const map_type> ColumnMap = Matrix->getColMap();
935  RCP<const import_type> importer = Matrix->getGraph()->getImporter();
936 
937  if (DomainMap->isContiguous()) {
938  // If the domain map is linear, then we can just use the column map as is.
939  return ColumnMap;
940  } else {
941  // The domain map isn't linear, so we need a new domain map
942  Teuchos::RCP<map_type> ContiguousDomainMap = rcp(new map_type(DomainMap->getGlobalNumElements(),
943  DomainMap->getLocalNumElements(), 0, DomainMap->getComm()));
944  if (importer) {
945  // If there's an importer then we can use it to get a new column map
946  go_vector_type MyGIDsHYPRE(DomainMap, ContiguousDomainMap->getLocalElementList());
947 
948  // import the HYPRE GIDs
949  go_vector_type ColGIDsHYPRE(ColumnMap);
950  ColGIDsHYPRE.doImport(MyGIDsHYPRE, *importer, Tpetra::INSERT);
951 
952  // Make a HYPRE numbering-based column map.
953  return Teuchos::rcp(new map_type(ColumnMap->getGlobalNumElements(), ColGIDsHYPRE.getDataNonConst()(), 0, ColumnMap->getComm()));
954  } else {
955  // The problem has matching domain/column maps, and somehow the domain map isn't linear, so just use the new domain map
956  return Teuchos::rcp(new map_type(ColumnMap->getGlobalNumElements(), ContiguousDomainMap->getLocalElementList(), 0, ColumnMap->getComm()));
957  }
958  }
959 }
960 
961 template <class LocalOrdinal, class Node>
962 int Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_Int, Node> >::getNumInitialize() const {
963  return NumInitialize_;
964 }
965 
966 template <class LocalOrdinal, class Node>
967 int Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_Int, Node> >::getNumCompute() const {
968  return NumCompute_;
969 }
970 
971 template <class LocalOrdinal, class Node>
972 int Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_Int, Node> >::getNumApply() const {
973  return NumApply_;
974 }
975 
976 template <class LocalOrdinal, class Node>
977 double Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_Int, Node> >::getInitializeTime() const {
978  return InitializeTime_;
979 }
980 
981 template <class LocalOrdinal, class Node>
982 double Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_Int, Node> >::getComputeTime() const {
983  return ComputeTime_;
984 }
985 
986 template <class LocalOrdinal, class Node>
987 double Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_Int, Node> >::getApplyTime() const {
988  return ApplyTime_;
989 }
990 
991 template <class LocalOrdinal, class Node>
992 Teuchos::RCP<const typename Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_Int, Node> >::map_type>
993 Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_Int, Node> >::
994  getDomainMap() const {
995  Teuchos::RCP<const row_matrix_type> A = getMatrix();
996  TEUCHOS_TEST_FOR_EXCEPTION(
997  A.is_null(), std::runtime_error,
998  "Ifpack2::Hypre::getDomainMap: The "
999  "input matrix A is null. Please call setMatrix() with a nonnull input "
1000  "matrix before calling this method.");
1001  return A->getDomainMap();
1002 }
1003 
1004 template <class LocalOrdinal, class Node>
1005 Teuchos::RCP<const typename Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_Int, Node> >::map_type>
1006 Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_Int, Node> >::
1007  getRangeMap() const {
1008  Teuchos::RCP<const row_matrix_type> A = getMatrix();
1009  TEUCHOS_TEST_FOR_EXCEPTION(
1010  A.is_null(), std::runtime_error,
1011  "Ifpack2::Hypre::getRangeMap: The "
1012  "input matrix A is null. Please call setMatrix() with a nonnull input "
1013  "matrix before calling this method.");
1014  return A->getRangeMap();
1015 }
1016 
1017 template <class LocalOrdinal, class Node>
1018 void Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_Int, Node> >::setMatrix(const Teuchos::RCP<const row_matrix_type> &A) {
1019  if (A.getRawPtr() != getMatrix().getRawPtr()) {
1020  IsInitialized_ = false;
1021  IsComputed_ = false;
1022  A_ = A;
1023  }
1024 }
1025 
1026 template <class LocalOrdinal, class Node>
1027 Teuchos::RCP<const typename Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_Int, Node> >::row_matrix_type>
1028 Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_Int, Node> >::
1029  getMatrix() const {
1030  return A_;
1031 }
1032 
1033 template <class LocalOrdinal, class Node>
1034 bool Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_Int, Node> >::hasTransposeApply() const {
1035  return false;
1036 }
1037 
1038 } // namespace Ifpack2
1039 
1040 #define IFPACK2_HYPRE_INSTANT(S, LO, GO, N) \
1041  template class Ifpack2::Hypre<Tpetra::RowMatrix<S, LO, GO, N> >;
1042 
1043 #endif // HAVE_HYPRE && HAVE_MPI
1044 #endif // IFPACK2_HYPRE_DEF_HPP
Teuchos::ParameterList::end
ConstIterator end() const
Teuchos::ParameterList::get
T & get(const std::string &name, T def_value)
Teuchos::TimeMonitor::getNewCounter
static RCP< Time > getNewCounter(const std::string &name)
Teuchos::TimeMonitor::lookupCounter
static RCP< Time > lookupCounter(const std::string &name)
TEUCHOS_TEST_FOR_EXCEPTION
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
Teuchos::EVerbosityLevel
EVerbosityLevel
Teuchos::ArrayView::data
T * data() const
Teuchos::ParameterList::isParameter
bool isParameter(const std::string &name) const
Teuchos::RCP::getRawPtr
T * getRawPtr() const
Teuchos::rcp
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Teuchos_ParameterList.hpp
Teuchos::ParameterList::isSublist
bool isSublist(const std::string &name) const
Teuchos::basic_FancyOStream
Teuchos_DefaultMpiComm.hpp
Teuchos::ParameterList::begin
ConstIterator begin() const
Teuchos::ParameterList::isType
bool isType(const std::string &name) const
Ifpack2::GMRES
Uses AztecOO&#39;s GMRES.
Definition: Ifpack2_CondestType.hpp:20
Teuchos::ParameterList::sublist
ParameterList & sublist(const std::string &name, bool mustAlreadyExist=false, const std::string &docString="")
Teuchos::Time::wallTime
static double wallTime()
Teuchos_RCP.hpp
Teuchos::RCP::is_null
bool is_null() const
Generated by   doxygen 1.8.5