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