Ifpack Package Browser (Single Doxygen Collection)  Development
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Ifpack_Hypre.cpp
Go to the documentation of this file.
1 /*@HEADER
2 // ***********************************************************************
3 //
4 // Ifpack: Object-Oriented Algebraic Preconditioner Package
5 // Copyright (2002) 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 #include "Ifpack_Hypre.h"
43 #if defined(HAVE_HYPRE) && defined(HAVE_MPI)
44 #include <stdexcept>
45 
46 #include "Ifpack_Utils.h"
47 #include "Epetra_MpiComm.h"
48 #include "Epetra_IntVector.h"
49 #include "Epetra_Import.h"
50 #include "Teuchos_ParameterList.hpp"
51 #include "Teuchos_RCP.hpp"
52 #include "HYPRE_IJ_mv.h"
53 #include "HYPRE_parcsr_ls.h"
54 #include "krylov.h"
55 #include "_hypre_parcsr_mv.h"
56 #include "_hypre_parcsr_ls.h"
57 #include "_hypre_IJ_mv.h"
58 #include "HYPRE_parcsr_mv.h"
59 #include "HYPRE.h"
60 
61 using Teuchos::RCP;
62 using Teuchos::rcp;
63 using Teuchos::rcpFromRef;
64 
65 // The Python script that generates the ParameterMap needs to be after these typedefs
66 typedef int (*int_func)(HYPRE_Solver, int);
67 typedef int (*double_func)(HYPRE_Solver, double);
68 typedef int (*double_int_func)(HYPRE_Solver, double, int);
69 typedef int (*int_int_func)(HYPRE_Solver, int, int);
70 typedef int (*int_star_func)(HYPRE_Solver, int*);
71 typedef int (*int_star_star_func)(HYPRE_Solver, int**);
72 typedef int (*double_star_func)(HYPRE_Solver, double*);
73 typedef int (*int_int_double_double_func)(HYPRE_Solver, int, int, double, double);
74 typedef int (*int_int_int_double_int_int_func)(HYPRE_Solver, int, int, int, double, int, int);
75 typedef int (*char_star_func)(HYPRE_Solver, char*);
76 
78 class FunctionParameter{
79  public:
81  FunctionParameter(Hypre_Chooser chooser, int_func funct, int param1) :
82  chooser_(chooser),
83  option_(0),
84  int_func_(funct),
85  int_param1_(param1) {}
86 
87  FunctionParameter(Hypre_Chooser chooser, std::string funct_name, int param1) :
88  chooser_(chooser),
89  option_(0),
90  int_func_(hypreMapIntFunc_.at(funct_name)),
91  int_param1_(param1) {}
92 
94  FunctionParameter(Hypre_Chooser chooser, double_func funct, double param1):
95  chooser_(chooser),
96  option_(1),
97  double_func_(funct),
98  double_param1_(param1) {}
99 
100  FunctionParameter(Hypre_Chooser chooser, std::string funct_name, double param1):
101  chooser_(chooser),
102  option_(1),
103  double_func_(hypreMapDoubleFunc_.at(funct_name)),
104  double_param1_(param1) {}
105 
107  FunctionParameter(Hypre_Chooser chooser, double_int_func funct, double param1, int param2):
108  chooser_(chooser),
109  option_(2),
110  double_int_func_(funct),
111  int_param1_(param2),
112  double_param1_(param1) {}
113 
114  FunctionParameter(Hypre_Chooser chooser, std::string funct_name, double param1, int param2):
115  chooser_(chooser),
116  option_(2),
117  double_int_func_(hypreMapDoubleIntFunc_.at(funct_name)),
118  int_param1_(param2),
119  double_param1_(param1) {}
120 
122  FunctionParameter(Hypre_Chooser chooser, int_int_func funct, int param1, int param2):
123  chooser_(chooser),
124  option_(3),
125  int_int_func_(funct),
126  int_param1_(param1),
127  int_param2_(param2) {}
128 
129  FunctionParameter(Hypre_Chooser chooser, std::string funct_name, int param1, int param2):
130  chooser_(chooser),
131  option_(3),
132  int_int_func_(hypreMapIntIntFunc_.at(funct_name)),
133  int_param1_(param1),
134  int_param2_(param2) {}
135 
137  FunctionParameter(Hypre_Chooser chooser, int_star_func funct, int *param1):
138  chooser_(chooser),
139  option_(4),
140  int_star_func_(funct),
141  int_star_param_(param1) {}
142 
143  FunctionParameter(Hypre_Chooser chooser, std::string funct_name, int *param1):
144  chooser_(chooser),
145  option_(4),
146  int_star_func_(hypreMapIntStarFunc_.at(funct_name)),
147  int_star_param_(param1) {}
148 
150  FunctionParameter(Hypre_Chooser chooser, double_star_func funct, double* param1):
151  chooser_(chooser),
152  option_(5),
153  double_star_func_(funct),
154  double_star_param_(param1) {}
155 
156  FunctionParameter(Hypre_Chooser chooser, std::string funct_name, double* param1):
157  chooser_(chooser),
158  option_(5),
159  double_star_func_(hypreMapDoubleStarFunc_.at(funct_name)),
160  double_star_param_(param1) {}
161 
163  FunctionParameter(Hypre_Chooser chooser, int_int_double_double_func funct, int param1, int param2, double param3, double param4):
164  chooser_(chooser),
165  option_(6),
166  int_int_double_double_func_(funct),
167  int_param1_(param1),
168  int_param2_(param2),
169  double_param1_(param3),
170  double_param2_(param4) {}
171 
172  FunctionParameter(Hypre_Chooser chooser, std::string funct_name, int param1, int param2, double param3, double param4):
173  chooser_(chooser),
174  option_(6),
175  int_int_double_double_func_(hypreMapIntIntDoubleDoubleFunc_.at(funct_name)),
176  int_param1_(param1),
177  int_param2_(param2),
178  double_param1_(param3),
179  double_param2_(param4) {}
180 
182  FunctionParameter(Hypre_Chooser chooser, int_star_star_func funct, int ** param1):
183  chooser_(chooser),
184  option_(7),
185  int_star_star_func_(funct),
186  int_star_star_param_(param1) {}
187 
188  FunctionParameter(Hypre_Chooser chooser, std::string funct_name, int** param1):
189  chooser_(chooser),
190  option_(7),
191  int_star_star_func_(hypreMapIntStarStarFunc_.at(funct_name)),
192  int_star_star_param_(param1) {}
193 
195  FunctionParameter(Hypre_Chooser chooser, int_int_int_double_int_int_func funct, int param1, int param2, int param3, double param4, int param5, int param6):
196  chooser_(chooser),
197  option_(8),
198  int_int_int_double_int_int_func_(funct),
199  int_param1_(param1),
200  int_param2_(param2),
201  int_param3_(param3),
202  int_param4_(param5),
203  int_param5_(param6),
204  double_param1_(param4) {}
205 
206  FunctionParameter(Hypre_Chooser chooser, std::string funct_name, int param1, int param2, int param3, double param4, int param5, int param6):
207  chooser_(chooser),
208  option_(8),
209  int_int_int_double_int_int_func_(hypreMapIntIntIntDoubleIntIntFunc_.at(funct_name)),
210  int_param1_(param1),
211  int_param2_(param2),
212  int_param3_(param3),
213  int_param4_(param5),
214  int_param5_(param6),
215  double_param1_(param4) {}
216 
218  FunctionParameter(Hypre_Chooser chooser, char_star_func funct, char *param1):
219  chooser_(chooser),
220  option_(9),
221  char_star_func_(funct),
222  char_star_param_(param1) {}
223 
224  FunctionParameter(Hypre_Chooser chooser, std::string funct_name, char *param1):
225  chooser_(chooser),
226  option_(9),
227  char_star_func_(hypreMapCharStarFunc_.at(funct_name)),
228  char_star_param_(param1) {}
229 
231  int CallFunction(HYPRE_Solver solver, HYPRE_Solver precond) {
232  if(chooser_ == Solver){
233  if(option_ == 0){
234  return int_func_(solver, int_param1_);
235  } else if(option_ == 1){
236  return double_func_(solver, double_param1_);
237  } else if(option_ == 2){
238  return double_int_func_(solver, double_param1_, int_param1_);
239  } else if (option_ == 3){
240  return int_int_func_(solver, int_param1_, int_param2_);
241  } else if (option_ == 4){
242  return int_star_func_(solver, int_star_param_);
243  } else if (option_ == 5){
244  return double_star_func_(solver, double_star_param_);
245  } else if (option_ == 6) {
246  return int_int_double_double_func_(solver, int_param1_, int_param2_, double_param1_, double_param2_);
247  } else if (option_ == 7) {
248  return int_star_star_func_(solver, int_star_star_param_);
249  } else if (option_ == 8) {
250  return int_int_int_double_int_int_func_(solver, int_param1_, int_param2_, int_param3_, double_param1_, int_param4_, int_param5_);
251  } else if (option_ == 9) {
252  return char_star_func_(solver, char_star_param_);
253  } else {
254  IFPACK_CHK_ERR(-2);
255  }
256  } else {
257  if(option_ == 0){
258  return int_func_(precond, int_param1_);
259  } else if(option_ == 1){
260  return double_func_(precond, double_param1_);
261  } else if(option_ == 2){
262  return double_int_func_(precond, double_param1_, int_param1_);
263  } else if(option_ == 3) {
264  return int_int_func_(precond, int_param1_, int_param2_);
265  } else if(option_ == 4) {
266  return int_star_func_(precond, int_star_param_);
267  } else if(option_ == 5) {
268  return double_star_func_(precond, double_star_param_);
269  } else if (option_ == 6) {
270  return int_int_double_double_func_(precond, int_param1_, int_param2_, double_param1_, double_param2_);
271  } else if (option_ == 7) {
272  return int_star_star_func_(precond, int_star_star_param_);
273  } else if (option_ == 8) {
274  return int_int_int_double_int_int_func_(precond, int_param1_, int_param2_, int_param3_, double_param1_, int_param4_, int_param5_);
275  } else if (option_ == 9) {
276  return char_star_func_(solver, char_star_param_);
277  } else {
278  IFPACK_CHK_ERR(-2);
279  }
280  }
281  }
282 
283  static bool isFuncIntInt(std::string funct_name) {
284  return (hypreMapIntIntFunc_.find(funct_name) != hypreMapIntIntFunc_.end());
285  }
286 
287  static bool isFuncIntIntDoubleDouble(std::string funct_name) {
288  return (hypreMapIntIntDoubleDoubleFunc_.find(funct_name) != hypreMapIntIntDoubleDoubleFunc_.end());
289  }
290 
291  static bool isFuncIntIntIntDoubleIntInt(std::string funct_name) {
292  return (hypreMapIntIntIntDoubleIntIntFunc_.find(funct_name) != hypreMapIntIntIntDoubleIntIntFunc_.end());
293  }
294 
295  static bool isFuncIntStarStar(std::string funct_name) {
296  return (hypreMapIntStarStarFunc_.find(funct_name) != hypreMapIntStarStarFunc_.end());
297  }
298 
299  private:
300  Hypre_Chooser chooser_;
301  int option_;
302  int_func int_func_;
303  double_func double_func_;
304  double_int_func double_int_func_;
305  int_int_func int_int_func_;
306  int_star_func int_star_func_;
307  double_star_func double_star_func_;
308  int_int_double_double_func int_int_double_double_func_;
309  int_int_int_double_int_int_func int_int_int_double_int_int_func_;
310  int_star_star_func int_star_star_func_;
311  char_star_func char_star_func_;
312  int int_param1_;
313  int int_param2_;
314  int int_param3_;
315  int int_param4_;
316  int int_param5_;
317  double double_param1_;
318  double double_param2_;
319  int *int_star_param_;
320  int **int_star_star_param_;
321  double *double_star_param_;
322  char *char_star_param_;
323 
324  static const std::map<std::string, int_func> hypreMapIntFunc_;
325  static const std::map<std::string, double_func> hypreMapDoubleFunc_;
326  static const std::map<std::string, double_int_func> hypreMapDoubleIntFunc_;
327  static const std::map<std::string, int_int_func> hypreMapIntIntFunc_;
328  static const std::map<std::string, int_star_func> hypreMapIntStarFunc_;
329  static const std::map<std::string, double_star_func> hypreMapDoubleStarFunc_;
330  static const std::map<std::string, int_int_double_double_func> hypreMapIntIntDoubleDoubleFunc_;
331  static const std::map<std::string, int_int_int_double_int_int_func> hypreMapIntIntIntDoubleIntIntFunc_;
332  static const std::map<std::string, int_star_star_func> hypreMapIntStarStarFunc_;
333  static const std::map<std::string, char_star_func> hypreMapCharStarFunc_;
334 
335 };
336 
337 // NOTE: This really, really needs to be here and not up above, so please don't move it
338 #include "Ifpack_HypreParameterMap.h"
339 
340 Ifpack_Hypre::Ifpack_Hypre(Epetra_RowMatrix* A):
341  A_(rcp(A,false)),
342  UseTranspose_(false),
343  IsInitialized_(false),
344  IsComputed_(false),
345  Label_(),
346  NumInitialize_(0),
347  NumCompute_(0),
348  NumApplyInverse_(0),
349  InitializeTime_(0.0),
350  ComputeTime_(0.0),
351  ApplyInverseTime_(0.0),
352  ComputeFlops_(0.0),
353  ApplyInverseFlops_(0.0),
354  Time_(A_->Comm()),
355  HypreA_(0),
356  HypreG_(0),
357  xHypre_(0),
358  yHypre_(0),
359  zHypre_(0),
360  IsSolverCreated_(false),
361  IsPrecondCreated_(false),
362  SolveOrPrec_(Solver),
363  NumFunsToCall_(0),
364  SolverType_(PCG),
365  PrecondType_(Euclid),
366  UsePreconditioner_(false),
367  Dump_(false)
368 {
369  MPI_Comm comm = GetMpiComm();
370  // Check that RowMap and RangeMap are the same. While this could handle the
371  // case where RowMap and RangeMap are permutations, other Ifpack PCs don't
372  // handle this either.
373  if (!A_->RowMatrixRowMap().SameAs(A_->OperatorRangeMap())) {
374  IFPACK_CHK_ERRV(-1);
375  }
376  // Hypre expects the RowMap to be Linear.
377  if (A_->RowMatrixRowMap().LinearMap()) {
378  // note these are non-owning pointers, they are deleted by A_'s destructor
379  GloballyContiguousRowMap_ = rcpFromRef(A_->RowMatrixRowMap());
380  GloballyContiguousColMap_ = rcpFromRef(A_->RowMatrixColMap());
381  } else {
382  // Must create GloballyContiguous Maps for Hypre
383  if(A_->OperatorDomainMap().SameAs(A_->RowMatrixRowMap())) {
385  GloballyContiguousColMap_ = MakeContiguousColumnMap(Aconst);
386  GloballyContiguousRowMap_ = rcp(new Epetra_Map(A_->RowMatrixRowMap().NumGlobalElements(),
387  A_->RowMatrixRowMap().NumMyElements(), 0, Comm()));
388  }
389  else {
390  throw std::runtime_error("Ifpack_Hypre: Unsupported map configuration: Row/Domain maps do not match");
391  }
392  }
393  // Next create vectors that will be used when ApplyInverse() is called
394  int ilower = GloballyContiguousRowMap_->MinMyGID();
395  int iupper = GloballyContiguousRowMap_->MaxMyGID();
396  // X in AX = Y
397  IFPACK_CHK_ERRV(HYPRE_IJVectorCreate(comm, ilower, iupper, &XHypre_));
398  IFPACK_CHK_ERRV(HYPRE_IJVectorSetObjectType(XHypre_, HYPRE_PARCSR));
399  IFPACK_CHK_ERRV(HYPRE_IJVectorInitialize(XHypre_));
400  IFPACK_CHK_ERRV(HYPRE_IJVectorAssemble(XHypre_));
401  IFPACK_CHK_ERRV(HYPRE_IJVectorGetObject(XHypre_, (void**) &ParX_));
402  XVec_ = Teuchos::rcp((hypre_ParVector *) hypre_IJVectorObject(((hypre_IJVector *) XHypre_)),false);
403 
404  // Y in AX = Y
405  IFPACK_CHK_ERRV(HYPRE_IJVectorCreate(comm, ilower, iupper, &YHypre_));
406  IFPACK_CHK_ERRV(HYPRE_IJVectorSetObjectType(YHypre_, HYPRE_PARCSR));
407  IFPACK_CHK_ERRV(HYPRE_IJVectorInitialize(YHypre_));
408  IFPACK_CHK_ERRV(HYPRE_IJVectorAssemble(YHypre_));
409  IFPACK_CHK_ERRV(HYPRE_IJVectorGetObject(YHypre_, (void**) &ParY_));
410  YVec_ = Teuchos::rcp((hypre_ParVector *) hypre_IJVectorObject(((hypre_IJVector *) YHypre_)),false);
411 
412  // Cache
413  VectorCache_.resize(A->RowMatrixRowMap().NumMyElements());
414 } //Constructor
415 
416 //==============================================================================
417 void Ifpack_Hypre::Destroy(){
418  if(IsInitialized()){
419  IFPACK_CHK_ERRV(HYPRE_IJMatrixDestroy(HypreA_));
420  }
421  IFPACK_CHK_ERRV(HYPRE_IJVectorDestroy(XHypre_));
422  IFPACK_CHK_ERRV(HYPRE_IJVectorDestroy(YHypre_));
423  if(IsSolverCreated_){
424  IFPACK_CHK_ERRV(SolverDestroyPtr_(Solver_));
425  }
426  if(IsPrecondCreated_){
427  IFPACK_CHK_ERRV(PrecondDestroyPtr_(Preconditioner_));
428  }
429 
430  // Maxwell
431  if(HypreG_) {
432  IFPACK_CHK_ERRV(HYPRE_IJMatrixDestroy(HypreG_));
433  }
434  if(xHypre_) {
435  IFPACK_CHK_ERRV(HYPRE_IJVectorDestroy(xHypre_));
436  }
437  if(yHypre_) {
438  IFPACK_CHK_ERRV(HYPRE_IJVectorDestroy(yHypre_));
439  }
440  if(zHypre_) {
441  IFPACK_CHK_ERRV(HYPRE_IJVectorDestroy(zHypre_));
442  }
443 } //Destroy()
444 
445 //==============================================================================
446 int Ifpack_Hypre::Initialize(){
447  Time_.ResetStartTime();
448  if(IsInitialized_) return 0;
449  // set flags
450  IsInitialized_=true;
451  NumInitialize_ = NumInitialize_ + 1;
452  InitializeTime_ = InitializeTime_ + Time_.ElapsedTime();
453  return 0;
454 } //Initialize()
455 
456 //==============================================================================
457 int Ifpack_Hypre::SetParameters(Teuchos::ParameterList& list){
458 
459  std::map<std::string, Hypre_Solver> solverMap;
460  solverMap["BoomerAMG"] = BoomerAMG;
461  solverMap["ParaSails"] = ParaSails;
462  solverMap["Euclid"] = Euclid;
463  solverMap["AMS"] = AMS;
464  solverMap["Hybrid"] = Hybrid;
465  solverMap["PCG"] = PCG;
466  solverMap["GMRES"] = GMRES;
467  solverMap["FlexGMRES"] = FlexGMRES;
468  solverMap["LGMRES"] = LGMRES;
469  solverMap["BiCGSTAB"] = BiCGSTAB;
470 
471  std::map<std::string, Hypre_Chooser> chooserMap;
472  chooserMap["Solver"] = Solver;
473  chooserMap["Preconditioner"] = Preconditioner;
474 
475  List_ = list;
476  Hypre_Solver solType;
477  if (list.isType<std::string>("hypre: Solver"))
478  solType = solverMap[list.get<std::string>("hypre: Solver")];
479  else
480  solType = list.get("hypre: Solver", PCG);
481  SolverType_ = solType;
482  Hypre_Solver precType;
483  if (list.isType<std::string>("hypre: Preconditioner"))
484  precType = solverMap[list.get<std::string>("hypre: Preconditioner")];
485  else
486  precType = list.get("hypre: Preconditioner", Euclid);
487  PrecondType_ = precType;
488  Hypre_Chooser chooser;
489  if (list.isType<std::string>("hypre: SolveOrPrecondition"))
490  chooser = chooserMap[list.get<std::string>("hypre: SolveOrPrecondition")];
491  else
492  chooser = list.get("hypre: SolveOrPrecondition", Solver);
493  SolveOrPrec_ = chooser;
494  bool SetPrecond = list.get("hypre: SetPreconditioner", false);
495  IFPACK_CHK_ERR(SetParameter(SetPrecond));
496  int NumFunctions = list.get("hypre: NumFunctions", 0);
497  FunsToCall_.clear();
498  NumFunsToCall_ = 0;
499  if(NumFunctions > 0){
500  RCP<FunctionParameter>* params = list.get<RCP<FunctionParameter>*>("hypre: Functions");
501  for(int i = 0; i < NumFunctions; i++){
502  IFPACK_CHK_ERR(AddFunToList(params[i]));
503  }
504  }
505 
506  if (list.isSublist("hypre: Solver functions")) {
507  Teuchos::ParameterList solverList = list.sublist("hypre: Solver functions");
508  for (auto it = solverList.begin(); it != solverList.end(); ++it) {
509  std::string funct_name = it->first;
510  if (it->second.isType<int>()) {
511  IFPACK_CHK_ERR(AddFunToList(rcp(new FunctionParameter(Solver, funct_name , Teuchos::getValue<int>(it->second)))));
512  } else if (it->second.isType<double>()) {
513  IFPACK_CHK_ERR(AddFunToList(rcp(new FunctionParameter(Solver, funct_name , Teuchos::getValue<double>(it->second)))));
514  } else {
515  IFPACK_CHK_ERR(-1);
516  }
517  }
518  }
519 
520  if (list.isSublist("hypre: Preconditioner functions")) {
521  Teuchos::ParameterList precList = list.sublist("hypre: Preconditioner functions");
522  for (auto it = precList.begin(); it != precList.end(); ++it) {
523  std::string funct_name = it->first;
524  if (it->second.isType<int>()) {
525  IFPACK_CHK_ERR(AddFunToList(rcp(new FunctionParameter(Preconditioner, funct_name , Teuchos::getValue<int>(it->second)))));
526  } else if (it->second.isType<double>()) {
527  IFPACK_CHK_ERR(AddFunToList(rcp(new FunctionParameter(Preconditioner, funct_name , Teuchos::getValue<double>(it->second)))));
528  } else if (it->second.isList()) {
529  Teuchos::ParameterList pl = Teuchos::getValue<Teuchos::ParameterList>(it->second);
530  if (FunctionParameter::isFuncIntInt(funct_name)) {
531  int arg0 = pl.get<int>("arg 0");
532  int arg1 = pl.get<int>("arg 1");
533  IFPACK_CHK_ERR(AddFunToList(rcp(new FunctionParameter(Preconditioner, funct_name , arg0, arg1))));
534  } else if (FunctionParameter::isFuncIntIntDoubleDouble(funct_name)) {
535  int arg0 = pl.get<int>("arg 0");
536  int arg1 = pl.get<int>("arg 1");
537  double arg2 = pl.get<double>("arg 2");
538  double arg3 = pl.get<double>("arg 3");
539  IFPACK_CHK_ERR(AddFunToList(rcp(new FunctionParameter(Preconditioner, funct_name , arg0, arg1, arg2, arg3))));
540  } else if (FunctionParameter::isFuncIntIntIntDoubleIntInt(funct_name)) {
541  int arg0 = pl.get<int>("arg 0");
542  int arg1 = pl.get<int>("arg 1");
543  int arg2 = pl.get<int>("arg 2");
544  double arg3 = pl.get<double>("arg 3");
545  int arg4 = pl.get<int>("arg 4");
546  int arg5 = pl.get<int>("arg 5");
547  IFPACK_CHK_ERR(AddFunToList(rcp(new FunctionParameter(Preconditioner, funct_name , arg0, arg1, arg2, arg3, arg4, arg5))));
548  } else {
549  IFPACK_CHK_ERR(-1);
550  }
551  }
552  }
553  }
554 
555  if (list.isSublist("Coordinates") && list.sublist("Coordinates").isType<Teuchos::RCP<Epetra_MultiVector> >("Coordinates"))
556  Coords_ = list.sublist("Coordinates").get<Teuchos::RCP<Epetra_MultiVector> >("Coordinates");
557  if (list.isSublist("Operators") && list.sublist("Operators").isType<Teuchos::RCP<const Epetra_CrsMatrix> >("G"))
558  G_ = list.sublist("Operators").get<Teuchos::RCP<const Epetra_CrsMatrix> >("G");
559 
560  Dump_ = list.get("hypre: Dump", false);
561 
562  return 0;
563 } //SetParameters()
564 
565 //==============================================================================
566 int Ifpack_Hypre::AddFunToList(RCP<FunctionParameter> NewFun){
567  NumFunsToCall_ = NumFunsToCall_+1;
568  FunsToCall_.resize(NumFunsToCall_);
569  FunsToCall_[NumFunsToCall_-1] = NewFun;
570  return 0;
571 } //AddFunToList()
572 
573 //==============================================================================
574 int Ifpack_Hypre::SetParameter(Hypre_Chooser chooser, int (*pt2Func)(HYPRE_Solver, int), int parameter){
575  RCP<FunctionParameter> temp = rcp(new FunctionParameter(chooser, pt2Func, parameter));
576  IFPACK_CHK_ERR(AddFunToList(temp));
577  return 0;
578 } //SetParameter() - int function pointer
579 
580 //==============================================================================
581 int Ifpack_Hypre::SetParameter(Hypre_Chooser chooser, int (*pt2Func)(HYPRE_Solver, double), double parameter){
582  RCP<FunctionParameter> temp = rcp(new FunctionParameter(chooser, pt2Func, parameter));
583  IFPACK_CHK_ERR(AddFunToList(temp));
584  return 0;
585 } //SetParameter() - double function pointer
586 
587 //==============================================================================
588 int Ifpack_Hypre::SetParameter(Hypre_Chooser chooser, int (*pt2Func)(HYPRE_Solver, double, int), double parameter1, int parameter2){
589  RCP<FunctionParameter> temp = rcp(new FunctionParameter(chooser, pt2Func, parameter1, parameter2));
590  IFPACK_CHK_ERR(AddFunToList(temp));
591  return 0;
592 } //SetParameter() - double,int function pointer
593 
594 //==============================================================================
595 int Ifpack_Hypre::SetParameter(Hypre_Chooser chooser, int (*pt2Func)(HYPRE_Solver, int, int), int parameter1, int parameter2){
596  RCP<FunctionParameter> temp = rcp(new FunctionParameter(chooser, pt2Func, parameter1, parameter2));
597  IFPACK_CHK_ERR(AddFunToList(temp));
598  return 0;
599 } //SetParameter() int,int function pointer
600 
601 //==============================================================================
602 int Ifpack_Hypre::SetParameter(Hypre_Chooser chooser, int (*pt2Func)(HYPRE_Solver, double*), double* parameter){
603  RCP<FunctionParameter> temp = rcp(new FunctionParameter(chooser, pt2Func, parameter));
604  IFPACK_CHK_ERR(AddFunToList(temp));
605  return 0;
606 } //SetParameter() - double* function pointer
607 
608 //==============================================================================
609 int Ifpack_Hypre::SetParameter(Hypre_Chooser chooser, int (*pt2Func)(HYPRE_Solver, int*), int* parameter){
610  RCP<FunctionParameter> temp = rcp(new FunctionParameter(chooser, pt2Func, parameter));
611  IFPACK_CHK_ERR(AddFunToList(temp));
612  return 0;
613 } //SetParameter() - int* function pointer
614 
615 //==============================================================================
616 int Ifpack_Hypre::SetParameter(Hypre_Chooser chooser, int (*pt2Func)(HYPRE_Solver, int**), int** parameter){
617  RCP<FunctionParameter> temp = rcp(new FunctionParameter(chooser, pt2Func, parameter));
618  IFPACK_CHK_ERR(AddFunToList(temp));
619  return 0;
620 } //SetParameter() - int** function pointer
621 
622 //==============================================================================
623 int Ifpack_Hypre::SetParameter(Hypre_Chooser chooser, Hypre_Solver solver){
624  if(chooser == Solver){
625  SolverType_ = solver;
626  } else {
627  PrecondType_ = solver;
628  }
629  return 0;
630 } //SetParameter() - set type of solver
631 
632 //==============================================================================
633 int Ifpack_Hypre::SetDiscreteGradient(Teuchos::RCP<const Epetra_CrsMatrix> G){
634 
635  // Sanity check
636  if(!A_->RowMatrixRowMap().SameAs(G->RowMap()))
637  throw std::runtime_error("Ifpack_Hypre: Edge map mismatch: A and discrete gradient");
638 
639  // Get the maps for the nodes (assuming the edge map from A is OK);
640  GloballyContiguousNodeRowMap_ = rcp(new Epetra_Map(G->DomainMap().NumGlobalElements(),
641  G->DomainMap().NumMyElements(), 0, Comm()));
642  Teuchos::RCP<const Epetra_RowMatrix> Grow = Teuchos::rcp_dynamic_cast<const Epetra_RowMatrix>(G);
643  GloballyContiguousNodeColMap_ = MakeContiguousColumnMap(Grow);
644 
645  // Start building G
646  MPI_Comm comm = GetMpiComm();
647  int ilower = GloballyContiguousRowMap_->MinMyGID();
648  int iupper = GloballyContiguousRowMap_->MaxMyGID();
649  int jlower = GloballyContiguousNodeRowMap_->MinMyGID();
650  int jupper = GloballyContiguousNodeRowMap_->MaxMyGID();
651  IFPACK_CHK_ERR(HYPRE_IJMatrixCreate(comm, ilower, iupper, jlower, jupper, &HypreG_));
652  IFPACK_CHK_ERR(HYPRE_IJMatrixSetObjectType(HypreG_, HYPRE_PARCSR));
653  IFPACK_CHK_ERR(HYPRE_IJMatrixInitialize(HypreG_));
654 
655  std::vector<int> new_indices(G->MaxNumEntries());
656  for(int i = 0; i < G->NumMyRows(); i++){
657  int numEntries;
658  double * values;
659  int *indices;
660  IFPACK_CHK_ERR(G->ExtractMyRowView(i, numEntries, values, indices));
661  for(int j = 0; j < numEntries; j++){
662  new_indices[j] = GloballyContiguousNodeColMap_->GID(indices[j]);
663  }
664  int GlobalRow[1];
665  GlobalRow[0] = GloballyContiguousRowMap_->GID(i);
666  IFPACK_CHK_ERR(HYPRE_IJMatrixSetValues(HypreG_, 1, &numEntries, GlobalRow, new_indices.data(), values));
667  }
668  IFPACK_CHK_ERR(HYPRE_IJMatrixAssemble(HypreG_));
669  IFPACK_CHK_ERR(HYPRE_IJMatrixGetObject(HypreG_, (void**)&ParMatrixG_));
670 
671  if (Dump_)
672  HYPRE_ParCSRMatrixPrint(ParMatrixG_,"G.mat");
673 
674  if(SolverType_ == AMS)
675  HYPRE_AMSSetDiscreteGradient(Solver_, ParMatrixG_);
676  if(PrecondType_ == AMS)
677  HYPRE_AMSSetDiscreteGradient(Preconditioner_, ParMatrixG_);
678 
679  return 0;
680 } //SetDiscreteGradient()
681 
682 //==============================================================================
683 int Ifpack_Hypre::SetCoordinates(Teuchos::RCP<Epetra_MultiVector> coords) {
684 
685  if(!G_.is_null() && !G_->DomainMap().SameAs(coords->Map()))
686  throw std::runtime_error("Ifpack_Hypre: Node map mismatch: G->DomainMap() and coords");
687 
688  if(SolverType_ != AMS && PrecondType_ != AMS)
689  return 0;
690 
691  double *xPtr;
692  double *yPtr;
693  double *zPtr;
694 
695  IFPACK_CHK_ERR(((*coords)(0))->ExtractView(&xPtr));
696  IFPACK_CHK_ERR(((*coords)(1))->ExtractView(&yPtr));
697  IFPACK_CHK_ERR(((*coords)(2))->ExtractView(&zPtr));
698 
699  MPI_Comm comm = GetMpiComm();
700  int NumEntries = coords->MyLength();
701  int * indices = GloballyContiguousNodeRowMap_->MyGlobalElements();
702 
703  int ilower = GloballyContiguousNodeRowMap_->MinMyGID();
704  int iupper = GloballyContiguousNodeRowMap_->MaxMyGID();
705 
706  if( NumEntries != iupper-ilower+1) {
707  std::cout<<"Ifpack_Hypre::SetCoordinates(): Error on rank "<<Comm().MyPID()<<": MyLength = "<<coords->MyLength()<<" GID range = ["<<ilower<<","<<iupper<<"]"<<std::endl;
708  throw std::runtime_error("Ifpack_Hypre: SetCoordinates: Length mismatch");
709  }
710 
711  IFPACK_CHK_ERR(HYPRE_IJVectorCreate(comm, ilower, iupper, &xHypre_));
712  IFPACK_CHK_ERR(HYPRE_IJVectorSetObjectType(xHypre_, HYPRE_PARCSR));
713  IFPACK_CHK_ERR(HYPRE_IJVectorInitialize(xHypre_));
714  IFPACK_CHK_ERR(HYPRE_IJVectorSetValues(xHypre_,NumEntries,indices,xPtr));
715  IFPACK_CHK_ERR(HYPRE_IJVectorAssemble(xHypre_));
716  IFPACK_CHK_ERR(HYPRE_IJVectorGetObject(xHypre_, (void**) &xPar_));
717 
718  IFPACK_CHK_ERR(HYPRE_IJVectorCreate(comm, ilower, iupper, &yHypre_));
719  IFPACK_CHK_ERR(HYPRE_IJVectorSetObjectType(yHypre_, HYPRE_PARCSR));
720  IFPACK_CHK_ERR(HYPRE_IJVectorInitialize(yHypre_));
721  IFPACK_CHK_ERR(HYPRE_IJVectorSetValues(yHypre_,NumEntries,indices,yPtr));
722  IFPACK_CHK_ERR(HYPRE_IJVectorAssemble(yHypre_));
723  IFPACK_CHK_ERR(HYPRE_IJVectorGetObject(yHypre_, (void**) &yPar_));
724 
725  IFPACK_CHK_ERR(HYPRE_IJVectorCreate(comm, ilower, iupper, &zHypre_));
726  IFPACK_CHK_ERR(HYPRE_IJVectorSetObjectType(zHypre_, HYPRE_PARCSR));
727  IFPACK_CHK_ERR(HYPRE_IJVectorInitialize(zHypre_));
728  IFPACK_CHK_ERR(HYPRE_IJVectorSetValues(zHypre_,NumEntries,indices,zPtr));
729  IFPACK_CHK_ERR(HYPRE_IJVectorAssemble(zHypre_));
730  IFPACK_CHK_ERR(HYPRE_IJVectorGetObject(zHypre_, (void**) &zPar_));
731 
732  if (Dump_) {
733  HYPRE_ParVectorPrint(xPar_,"coordX.dat");
734  HYPRE_ParVectorPrint(yPar_,"coordY.dat");
735  HYPRE_ParVectorPrint(zPar_,"coordZ.dat");
736  }
737 
738  if(SolverType_ == AMS)
739  HYPRE_AMSSetCoordinateVectors(Solver_, xPar_, yPar_, zPar_);
740  if(PrecondType_ == AMS)
741  HYPRE_AMSSetCoordinateVectors(Preconditioner_, xPar_, yPar_, zPar_);
742 
743  return 0;
744 
745 } //SetCoordinates
746 
747 //==============================================================================
748 int Ifpack_Hypre::Compute(){
749  if(IsInitialized() == false){
750  IFPACK_CHK_ERR(Initialize());
751  }
752  Time_.ResetStartTime();
753 
754  // Create the Hypre matrix and copy values. Note this uses values (which
755  // Initialize() shouldn't do) but it doesn't care what they are (for
756  // instance they can be uninitialized data even). It should be possible to
757  // set the Hypre structure without copying values, but this is the easiest
758  // way to get the structure.
759  MPI_Comm comm = GetMpiComm();
760  int ilower = GloballyContiguousRowMap_->MinMyGID();
761  int iupper = GloballyContiguousRowMap_->MaxMyGID();
762  IFPACK_CHK_ERR(HYPRE_IJMatrixCreate(comm, ilower, iupper, ilower, iupper, &HypreA_));
763  IFPACK_CHK_ERR(HYPRE_IJMatrixSetObjectType(HypreA_, HYPRE_PARCSR));
764  IFPACK_CHK_ERR(HYPRE_IJMatrixInitialize(HypreA_));
765  CopyEpetraToHypre();
766  if(SolveOrPrec_ == Solver) {
767  IFPACK_CHK_ERR(SetSolverType(SolverType_));
768  if (SolverPrecondPtr_ != NULL && UsePreconditioner_) {
769  // both method allows a PC (first condition) and the user wants a PC (second)
770  IFPACK_CHK_ERR(SetPrecondType(PrecondType_));
771  CallFunctions();
772  IFPACK_CHK_ERR(SolverPrecondPtr_(Solver_, PrecondSolvePtr_, PrecondSetupPtr_, Preconditioner_));
773  } else {
774  CallFunctions();
775  }
776  } else {
777  IFPACK_CHK_ERR(SetPrecondType(PrecondType_));
778  CallFunctions();
779  }
780 
781  if (!G_.is_null()) {
782  SetDiscreteGradient(G_);
783  }
784 
785  if (!Coords_.is_null()) {
786  SetCoordinates(Coords_);
787  }
788 
789  // Hypre Setup must be called after matrix has values
790  if(SolveOrPrec_ == Solver){
791  IFPACK_CHK_ERR(SolverSetupPtr_(Solver_, ParMatrix_, ParX_, ParY_));
792  } else {
793  IFPACK_CHK_ERR(PrecondSetupPtr_(Preconditioner_, ParMatrix_, ParX_, ParY_));
794  }
795 
796  // Dump Hierarchy here for BoomerAMG Preconditioner
797  if(Dump_ && PrecondSolvePtr_ == &HYPRE_BoomerAMGSolve) {
798  hypre_ParAMGData *amg_data = (hypre_ParAMGData*) Preconditioner_;
799  hypre_ParCSRMatrix **A_array = hypre_ParAMGDataAArray(amg_data);
800  hypre_ParCSRMatrix **P_array = hypre_ParAMGDataPArray(amg_data);
801  HYPRE_Int **CF_marker_array = hypre_ParAMGDataCFMarkerArray(amg_data);
802  HYPRE_Int num_levels = hypre_ParAMGDataNumLevels(amg_data);
803 
804  char ofs[80];
805  for(int k=0; k<num_levels; k++) {
806  // A
807  sprintf(ofs,"A_matrix.bmg.%d.dat",k);
808  HYPRE_ParCSRMatrixPrint(A_array[k], ofs);
809  if(k!=num_levels-1) {
810  // P
811  sprintf(ofs,"P_matrix.bmg.%d.dat",k);
812  HYPRE_ParCSRMatrixPrint(P_array[k], ofs);
813 
814  // CF
815  // Note: Hypre outputs "-1" for F Points and "1" for C Points
816  HYPRE_Int local_size = hypre_CSRMatrixNumRows(hypre_ParCSRMatrixDiag(A_array[k]));
817  sprintf(ofs,"cf_marker.bmg.%d.dat.%d",k,Comm().MyPID());
818  FILE * f = fopen(ofs,"w");
819  fprintf(f,"%%%%MatrixMarket matrix array real general\n");
820  fprintf(f,"%d 1\n",local_size);
821  for(int i=0; i<local_size; i++)
822  fprintf(f,"%d\n",(int)CF_marker_array[k][i]);
823  fclose(f);
824  }
825 
826  }
827  }//end dump for BoomerAMG
828 
829 
830  IsComputed_ = true;
831  NumCompute_ = NumCompute_ + 1;
832  ComputeTime_ = ComputeTime_ + Time_.ElapsedTime();
833  return 0;
834 } //Compute()
835 
836 //==============================================================================
837 int Ifpack_Hypre::CallFunctions() const{
838  for(int i = 0; i < NumFunsToCall_; i++){
839  IFPACK_CHK_ERR(FunsToCall_[i]->CallFunction(Solver_, Preconditioner_));
840  }
841  return 0;
842 } //CallFunctions()
843 
844 //==============================================================================
845 int Ifpack_Hypre::ApplyInverse(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const{
846  if(IsComputed() == false){
847  IFPACK_CHK_ERR(-1);
848  }
849  Time_.ResetStartTime();
850  hypre_Vector *XLocal_ = hypre_ParVectorLocalVector(XVec_);
851  hypre_Vector *YLocal_ = hypre_ParVectorLocalVector(YVec_);
852 
853  bool SameVectors = false;
854  int NumVectors = X.NumVectors();
855  if (NumVectors != Y.NumVectors()) IFPACK_CHK_ERR(-1); // X and Y must have same number of vectors
856  if(X.Pointers() == Y.Pointers() || (NumVectors == 1 && X[0] == Y[0])){
857  SameVectors = true;
858  }
859 
860  // NOTE: Here were assuming that the local ordering of Epetra's X/Y-vectors and
861  // Hypre's X/Y-vectors are the same. Seeing as as this is more or less how we constructed
862  // the Hypre matrices, this seems pretty reasoanble.
863 
864  for(int VecNum = 0; VecNum < NumVectors; VecNum++) {
865  //Get values for current vector in multivector.
866  // FIXME amk Nov 23, 2015: This will not work for funky data layouts
867  double * XValues = const_cast<double*>(X[VecNum]);
868  double * YValues;
869  if(!SameVectors){
870  YValues = const_cast<double*>(Y[VecNum]);
871  } else {
872  YValues = VectorCache_.getRawPtr();
873  }
874  // Temporarily make a pointer to data in Hypre for end
875  double *XTemp = XLocal_->data;
876  // Replace data in Hypre vectors with Epetra data
877  XLocal_->data = XValues;
878  double *YTemp = YLocal_->data;
879  YLocal_->data = YValues;
880 
881  IFPACK_CHK_ERR(HYPRE_ParVectorSetConstantValues(ParY_, 0.0));
882  if(SolveOrPrec_ == Solver){
883  // Use the solver methods
884  IFPACK_CHK_ERR(SolverSolvePtr_(Solver_, ParMatrix_, ParX_, ParY_));
885  } else {
886  // Apply the preconditioner
887  IFPACK_CHK_ERR(PrecondSolvePtr_(Preconditioner_, ParMatrix_, ParX_, ParY_));
888  }
889  if(SameVectors){
890  int NumEntries = Y.MyLength();
891  for(int i = 0; i < NumEntries; i++)
892  Y[VecNum][i] = YValues[i];
893  }
894  XLocal_->data = XTemp;
895  YLocal_->data = YTemp;
896  }
897  NumApplyInverse_ = NumApplyInverse_ + 1;
898  ApplyInverseTime_ = ApplyInverseTime_ + Time_.ElapsedTime();
899  return 0;
900 } //ApplyInverse()
901 
902 //==============================================================================
903 int Ifpack_Hypre::Multiply(bool TransA, const Epetra_MultiVector& X, Epetra_MultiVector& Y) const{
904  if(IsInitialized() == false){
905  IFPACK_CHK_ERR(-1);
906  }
907  hypre_Vector *XLocal_ = hypre_ParVectorLocalVector(XVec_);
908  hypre_Vector *YLocal_ = hypre_ParVectorLocalVector(YVec_);
909  bool SameVectors = false;
910  int NumVectors = X.NumVectors();
911  if (NumVectors != Y.NumVectors()) IFPACK_CHK_ERR(-1); // X and Y must have same number of vectors
912  if(X.Pointers() == Y.Pointers() || (NumVectors == 1 && X[0] == Y[0])){
913  SameVectors = true;
914  }
915 
916  // NOTE: Here were assuming that the local ordering of Epetra's X/Y-vectors and
917  // Hypre's X/Y-vectors are the same. Seeing as as this is more or less how we constructed
918  // the Hypre matrices, this seems pretty reasoanble.
919 
920  for(int VecNum = 0; VecNum < NumVectors; VecNum++) {
921  //Get values for current vector in multivector.
922  double * XValues=const_cast<double*>(X[VecNum]);
923  double * YValues;
924  double *XTemp = XLocal_->data;
925  double *YTemp = YLocal_->data;
926  if(!SameVectors){
927  YValues = const_cast<double*>(Y[VecNum]);
928  } else {
929  YValues = VectorCache_.getRawPtr();
930  }
931  YLocal_->data = YValues;
932  IFPACK_CHK_ERR(HYPRE_ParVectorSetConstantValues(ParY_,0.0));
933  // Temporarily make a pointer to data in Hypre for end
934  // Replace data in Hypre vectors with epetra values
935  XLocal_->data = XValues;
936  // Do actual computation.
937  if(TransA) {
938  // Use transpose of A in multiply
939  IFPACK_CHK_ERR(HYPRE_ParCSRMatrixMatvecT(1.0, ParMatrix_, ParX_, 1.0, ParY_));
940  } else {
941  IFPACK_CHK_ERR(HYPRE_ParCSRMatrixMatvec(1.0, ParMatrix_, ParX_, 1.0, ParY_));
942  }
943  if(SameVectors){
944  int NumEntries = Y.MyLength();
945  for(int i = 0; i < NumEntries; i++)
946  Y[VecNum][i] = YValues[i];
947  }
948  XLocal_->data = XTemp;
949  YLocal_->data = YTemp;
950  }
951  return 0;
952 } //Multiply()
953 
954 //==============================================================================
955 std::ostream& Ifpack_Hypre::Print(std::ostream& os) const{
956  using std::endl;
957  if (!Comm().MyPID()) {
958  os << endl;
959  os << "================================================================================" << endl;
960  os << "Ifpack_Hypre: " << Label() << endl << endl;
961  os << "Using " << Comm().NumProc() << " processors." << endl;
962  os << "Global number of rows = " << A_->NumGlobalRows() << endl;
963  os << "Global number of nonzeros = " << A_->NumGlobalNonzeros() << endl;
964  os << "Condition number estimate = " << Condest() << endl;
965  os << endl;
966  os << "Phase # calls Total Time (s) Total MFlops MFlops/s" << endl;
967  os << "----- ------- -------------- ------------ --------" << endl;
968  os << "Initialize() " << std::setw(5) << NumInitialize_
969  << " " << std::setw(15) << InitializeTime_
970  << " 0.0 0.0" << endl;
971  os << "Compute() " << std::setw(5) << NumCompute_
972  << " " << std::setw(15) << ComputeTime_
973  << " " << std::setw(15) << 1.0e-6 * ComputeFlops_;
974  if (ComputeTime_ != 0.0)
975  os << " " << std::setw(15) << 1.0e-6 * ComputeFlops_ / ComputeTime_ << endl;
976  else
977  os << " " << std::setw(15) << 0.0 << endl;
978  os << "ApplyInverse() " << std::setw(5) << NumApplyInverse_
979  << " " << std::setw(15) << ApplyInverseTime_
980  << " " << std::setw(15) << 1.0e-6 * ApplyInverseFlops_;
981  if (ApplyInverseTime_ != 0.0)
982  os << " " << std::setw(15) << 1.0e-6 * ApplyInverseFlops_ / ApplyInverseTime_ << endl;
983  else
984  os << " " << std::setw(15) << 0.0 << endl;
985  os << "================================================================================" << endl;
986  os << endl;
987  }
988  return os;
989 } //Print()
990 
991 //==============================================================================
992 double Ifpack_Hypre::Condest(const Ifpack_CondestType CT,
993  const int MaxIters,
994  const double Tol,
995  Epetra_RowMatrix* Matrix_in){
996  if (!IsComputed()) // cannot compute right now
997  return(-1.0);
998  Condest_ = Ifpack_Condest(*this, CT, MaxIters, Tol, Matrix_in);
999  return(Condest_);
1000 } //Condest()
1001 
1002 //==============================================================================
1003 int Ifpack_Hypre::SetSolverType(Hypre_Solver Solver){
1004  switch(Solver) {
1005  case BoomerAMG:
1006  if(IsSolverCreated_){
1007  SolverDestroyPtr_(Solver_);
1008  IsSolverCreated_ = false;
1009  }
1010  SolverCreatePtr_ = &Ifpack_Hypre::Hypre_BoomerAMGCreate;
1011  SolverDestroyPtr_ = &HYPRE_BoomerAMGDestroy;
1012  SolverSetupPtr_ = &HYPRE_BoomerAMGSetup;
1013  SolverPrecondPtr_ = NULL;
1014  SolverSolvePtr_ = &HYPRE_BoomerAMGSolve;
1015  break;
1016  case AMS:
1017  if(IsSolverCreated_){
1018  SolverDestroyPtr_(Solver_);
1019  IsSolverCreated_ = false;
1020  }
1021  SolverCreatePtr_ = &Ifpack_Hypre::Hypre_AMSCreate;
1022  SolverDestroyPtr_ = &HYPRE_AMSDestroy;
1023  SolverSetupPtr_ = &HYPRE_AMSSetup;
1024  SolverSolvePtr_ = &HYPRE_AMSSolve;
1025  SolverPrecondPtr_ = NULL;
1026  break;
1027  case Hybrid:
1028  if(IsSolverCreated_){
1029  SolverDestroyPtr_(Solver_);
1030  IsSolverCreated_ = false;
1031  }
1032  SolverCreatePtr_ = &Ifpack_Hypre::Hypre_ParCSRHybridCreate;
1033  SolverDestroyPtr_ = &HYPRE_ParCSRHybridDestroy;
1034  SolverSetupPtr_ = &HYPRE_ParCSRHybridSetup;
1035  SolverSolvePtr_ = &HYPRE_ParCSRHybridSolve;
1036  SolverPrecondPtr_ = &HYPRE_ParCSRHybridSetPrecond;
1037  break;
1038  case PCG:
1039  if(IsSolverCreated_){
1040  SolverDestroyPtr_(Solver_);
1041  IsSolverCreated_ = false;
1042  }
1043  SolverCreatePtr_ = &Ifpack_Hypre::Hypre_ParCSRPCGCreate;
1044  SolverDestroyPtr_ = &HYPRE_ParCSRPCGDestroy;
1045  SolverSetupPtr_ = &HYPRE_ParCSRPCGSetup;
1046  SolverSolvePtr_ = &HYPRE_ParCSRPCGSolve;
1047  SolverPrecondPtr_ = &HYPRE_ParCSRPCGSetPrecond;
1048  break;
1049  case GMRES:
1050  if(IsSolverCreated_){
1051  SolverDestroyPtr_(Solver_);
1052  IsSolverCreated_ = false;
1053  }
1054  SolverCreatePtr_ = &Ifpack_Hypre::Hypre_ParCSRGMRESCreate;
1055  SolverDestroyPtr_ = &HYPRE_ParCSRGMRESDestroy;
1056  SolverSetupPtr_ = &HYPRE_ParCSRGMRESSetup;
1057  SolverPrecondPtr_ = &HYPRE_ParCSRGMRESSetPrecond;
1058  break;
1059  case FlexGMRES:
1060  if(IsSolverCreated_){
1061  SolverDestroyPtr_(Solver_);
1062  IsSolverCreated_ = false;
1063  }
1064  SolverCreatePtr_ = &Ifpack_Hypre::Hypre_ParCSRFlexGMRESCreate;
1065  SolverDestroyPtr_ = &HYPRE_ParCSRFlexGMRESDestroy;
1066  SolverSetupPtr_ = &HYPRE_ParCSRFlexGMRESSetup;
1067  SolverSolvePtr_ = &HYPRE_ParCSRFlexGMRESSolve;
1068  SolverPrecondPtr_ = &HYPRE_ParCSRFlexGMRESSetPrecond;
1069  break;
1070  case LGMRES:
1071  if(IsSolverCreated_){
1072  SolverDestroyPtr_(Solver_);
1073  IsSolverCreated_ = false;
1074  }
1075  SolverCreatePtr_ = &Ifpack_Hypre::Hypre_ParCSRLGMRESCreate;
1076  SolverDestroyPtr_ = &HYPRE_ParCSRLGMRESDestroy;
1077  SolverSetupPtr_ = &HYPRE_ParCSRLGMRESSetup;
1078  SolverSolvePtr_ = &HYPRE_ParCSRLGMRESSolve;
1079  SolverPrecondPtr_ = &HYPRE_ParCSRLGMRESSetPrecond;
1080  break;
1081  case BiCGSTAB:
1082  if(IsSolverCreated_){
1083  SolverDestroyPtr_(Solver_);
1084  IsSolverCreated_ = false;
1085  }
1086  SolverCreatePtr_ = &Ifpack_Hypre::Hypre_ParCSRBiCGSTABCreate;
1087  SolverDestroyPtr_ = &HYPRE_ParCSRBiCGSTABDestroy;
1088  SolverSetupPtr_ = &HYPRE_ParCSRBiCGSTABSetup;
1089  SolverSolvePtr_ = &HYPRE_ParCSRBiCGSTABSolve;
1090  SolverPrecondPtr_ = &HYPRE_ParCSRBiCGSTABSetPrecond;
1091  break;
1092  default:
1093  return -1;
1094  }
1095  CreateSolver();
1096  return 0;
1097 } //SetSolverType()
1098 
1099 //==============================================================================
1100 int Ifpack_Hypre::SetPrecondType(Hypre_Solver Precond){
1101  switch(Precond) {
1102  case BoomerAMG:
1103  if(IsPrecondCreated_){
1104  PrecondDestroyPtr_(Preconditioner_);
1105  IsPrecondCreated_ = false;
1106  }
1107  PrecondCreatePtr_ = &Ifpack_Hypre::Hypre_BoomerAMGCreate;
1108  PrecondDestroyPtr_ = &HYPRE_BoomerAMGDestroy;
1109  PrecondSetupPtr_ = &HYPRE_BoomerAMGSetup;
1110  PrecondSolvePtr_ = &HYPRE_BoomerAMGSolve;
1111  break;
1112  case ParaSails:
1113  if(IsPrecondCreated_){
1114  PrecondDestroyPtr_(Preconditioner_);
1115  IsPrecondCreated_ = false;
1116  }
1117  PrecondCreatePtr_ = &Ifpack_Hypre::Hypre_ParaSailsCreate;
1118  PrecondDestroyPtr_ = &HYPRE_ParaSailsDestroy;
1119  PrecondSetupPtr_ = &HYPRE_ParaSailsSetup;
1120  PrecondSolvePtr_ = &HYPRE_ParaSailsSolve;
1121  break;
1122  case Euclid:
1123  if(IsPrecondCreated_){
1124  PrecondDestroyPtr_(Preconditioner_);
1125  IsPrecondCreated_ = false;
1126  }
1127  PrecondCreatePtr_ = &Ifpack_Hypre::Hypre_EuclidCreate;
1128  PrecondDestroyPtr_ = &HYPRE_EuclidDestroy;
1129  PrecondSetupPtr_ = &HYPRE_EuclidSetup;
1130  PrecondSolvePtr_ = &HYPRE_EuclidSolve;
1131  break;
1132  case AMS:
1133  if(IsPrecondCreated_){
1134  PrecondDestroyPtr_(Preconditioner_);
1135  IsPrecondCreated_ = false;
1136  }
1137  PrecondCreatePtr_ = &Ifpack_Hypre::Hypre_AMSCreate;
1138  PrecondDestroyPtr_ = &HYPRE_AMSDestroy;
1139  PrecondSetupPtr_ = &HYPRE_AMSSetup;
1140  PrecondSolvePtr_ = &HYPRE_AMSSolve;
1141  break;
1142  default:
1143  return -1;
1144  }
1145  CreatePrecond();
1146  return 0;
1147 
1148 } //SetPrecondType()
1149 
1150 //==============================================================================
1151 int Ifpack_Hypre::CreateSolver(){
1152  MPI_Comm comm;
1153  HYPRE_ParCSRMatrixGetComm(ParMatrix_, &comm);
1154  int ierr = (this->*SolverCreatePtr_)(comm, &Solver_);
1155  IsSolverCreated_ = true;
1156  return ierr;
1157 } //CreateSolver()
1158 
1159 //==============================================================================
1160 int Ifpack_Hypre::CreatePrecond(){
1161  MPI_Comm comm;
1162  HYPRE_ParCSRMatrixGetComm(ParMatrix_, &comm);
1163  int ierr = (this->*PrecondCreatePtr_)(comm, &Preconditioner_);
1164  IsPrecondCreated_ = true;
1165  return ierr;
1166 } //CreatePrecond()
1167 
1168 
1169 //==============================================================================
1170 int Ifpack_Hypre::CopyEpetraToHypre(){
1171  Teuchos::RCP<const Epetra_CrsMatrix> Matrix = Teuchos::rcp_dynamic_cast<const Epetra_CrsMatrix>(A_);
1172  if(Matrix.is_null())
1173  throw std::runtime_error("Ifpack_Hypre: Unsupported matrix configuration: Epetra_CrsMatrix required");
1174 
1175  std::vector<int> new_indices(Matrix->MaxNumEntries());
1176  for(int i = 0; i < Matrix->NumMyRows(); i++){
1177  int numEntries;
1178  int *indices;
1179  double *values;
1180  IFPACK_CHK_ERR(Matrix->ExtractMyRowView(i, numEntries, values, indices));
1181  for(int j = 0; j < numEntries; j++){
1182  new_indices[j] = GloballyContiguousColMap_->GID(indices[j]);
1183  }
1184  int GlobalRow[1];
1185  GlobalRow[0] = GloballyContiguousRowMap_->GID(i);
1186  IFPACK_CHK_ERR(HYPRE_IJMatrixSetValues(HypreA_, 1, &numEntries, GlobalRow, new_indices.data(), values));
1187  }
1188  IFPACK_CHK_ERR(HYPRE_IJMatrixAssemble(HypreA_));
1189  IFPACK_CHK_ERR(HYPRE_IJMatrixGetObject(HypreA_, (void**)&ParMatrix_));
1190  if (Dump_)
1191  HYPRE_ParCSRMatrixPrint(ParMatrix_,"A.mat");
1192  return 0;
1193 } //CopyEpetraToHypre()
1194 
1195 //==============================================================================
1196 int Ifpack_Hypre::Hypre_BoomerAMGCreate(MPI_Comm /*comm*/, HYPRE_Solver *solver)
1197  { return HYPRE_BoomerAMGCreate(solver);}
1198 
1199 //==============================================================================
1200 int Ifpack_Hypre::Hypre_ParaSailsCreate(MPI_Comm comm, HYPRE_Solver *solver)
1201  { return HYPRE_ParaSailsCreate(comm, solver);}
1202 
1203 //==============================================================================
1204 int Ifpack_Hypre::Hypre_EuclidCreate(MPI_Comm comm, HYPRE_Solver *solver)
1205  { return HYPRE_EuclidCreate(comm, solver);}
1206 
1207 //==============================================================================
1208 int Ifpack_Hypre::Hypre_AMSCreate(MPI_Comm /*comm*/, HYPRE_Solver *solver)
1209  { return HYPRE_AMSCreate(solver);}
1210 
1211 //==============================================================================
1212 int Ifpack_Hypre::Hypre_ParCSRHybridCreate(MPI_Comm /*comm*/, HYPRE_Solver *solver)
1213  { return HYPRE_ParCSRHybridCreate(solver);}
1214 
1215 //==============================================================================
1216 int Ifpack_Hypre::Hypre_ParCSRPCGCreate(MPI_Comm comm, HYPRE_Solver *solver)
1217  { return HYPRE_ParCSRPCGCreate(comm, solver);}
1218 
1219 //==============================================================================
1220 int Ifpack_Hypre::Hypre_ParCSRGMRESCreate(MPI_Comm comm, HYPRE_Solver *solver)
1221  { return HYPRE_ParCSRGMRESCreate(comm, solver);}
1222 
1223 //==============================================================================
1224 int Ifpack_Hypre::Hypre_ParCSRFlexGMRESCreate(MPI_Comm comm, HYPRE_Solver *solver)
1225  { return HYPRE_ParCSRFlexGMRESCreate(comm, solver);}
1226 
1227 //==============================================================================
1228 int Ifpack_Hypre::Hypre_ParCSRLGMRESCreate(MPI_Comm comm, HYPRE_Solver *solver)
1229  { return HYPRE_ParCSRLGMRESCreate(comm, solver);}
1230 
1231 //==============================================================================
1232 int Ifpack_Hypre::Hypre_ParCSRBiCGSTABCreate(MPI_Comm comm, HYPRE_Solver *solver)
1233  { return HYPRE_ParCSRBiCGSTABCreate(comm, solver);}
1234 
1235 //==============================================================================
1236 Teuchos::RCP<const Epetra_Map> Ifpack_Hypre::MakeContiguousColumnMap(Teuchos::RCP<const Epetra_RowMatrix> &MatrixRow) const{
1237  // Must create GloballyContiguous DomainMap (which is a permutation of Matrix_'s
1238  // DomainMap) and the corresponding permuted ColumnMap.
1239  // Epetra_GID ---------> LID ----------> HYPRE_GID
1240  // via DomainMap.LID() via GloballyContiguousDomainMap.GID()
1241  Teuchos::RCP<const Epetra_CrsMatrix> Matrix = Teuchos::rcp_dynamic_cast<const Epetra_CrsMatrix>(MatrixRow);
1242  if(Matrix.is_null())
1243  throw std::runtime_error("Ifpack_Hypre: Unsupported matrix configuration: Epetra_CrsMatrix required");
1244  const Epetra_Map & DomainMap = Matrix->DomainMap();
1245  const Epetra_Map & ColumnMap = Matrix->ColMap();
1246  const Epetra_Import * importer = Matrix->Importer();
1247 
1248  if(DomainMap.LinearMap() ) {
1249  // If the domain map is linear, then we can just use the column map as is.
1250  return rcpFromRef(ColumnMap);
1251  }
1252  else {
1253  // The domain map isn't linear, so we need a new domain map
1254  Teuchos::RCP<Epetra_Map> ContiguousDomainMap = rcp(new Epetra_Map(DomainMap.NumGlobalElements(),
1255  DomainMap.NumMyElements(), 0, Comm()));
1256  if(importer) {
1257  // If there's an importer then we can use it to get a new column map
1258  Epetra_IntVector MyGIDsHYPRE(View,DomainMap,ContiguousDomainMap->MyGlobalElements());
1259 
1260  // import the HYPRE GIDs
1261  Epetra_IntVector ColGIDsHYPRE(ColumnMap);
1262  ColGIDsHYPRE.Import(MyGIDsHYPRE, *importer, Insert);
1263 
1264  // Make a HYPRE numbering-based column map.
1265  return Teuchos::rcp(new Epetra_Map(ColumnMap.NumGlobalElements(),ColGIDsHYPRE.MyLength(), &ColGIDsHYPRE[0], 0, Comm()));
1266  }
1267  else {
1268  // The problem has matching domain/column maps, and somehow the domain map isn't linear, so just use the new domain map
1269  return Teuchos::rcp(new Epetra_Map(ColumnMap.NumGlobalElements(),ColumnMap.NumMyElements(), ContiguousDomainMap->MyGlobalElements(), 0, Comm()));
1270  }
1271  }
1272 }
1273 
1274 
1275 
1276 #endif // HAVE_HYPRE && HAVE_MPI
int NumGlobalElements() const
ConstIterator end() const
void f()
virtual const Epetra_Map & RowMatrixRowMap() const =0
int MyGlobalElements(int *MyGlobalElementList) const
T & get(ParameterList &l, const std::string &name)
const int NumVectors
Definition: performance.cpp:71
RCP< ParameterList > sublist(const RCP< ParameterList > &paramList, const std::string &name, bool mustAlreadyExist=false, const std::string &docString="")
const Epetra_Map & ColMap() const
static bool Solver
Definition: performance.cpp:70
int ExtractMyRowView(int MyRow, int &NumEntries, double *&Values, int *&Indices) const
Ifpack_CondestType
Ifpack_CondestType: enum to define the type of condition number estimate.
const Epetra_Map & RowMap() const
int NumMyElements() const
#define IFPACK_CHK_ERRV(ifpack_err)
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
const Epetra_Import * Importer() const
bool is_null(const RCP< T > &p)
bool isSublist(const std::string &name) const
int MaxNumEntries() const
Hypre_Chooser
virtual const Epetra_BlockMap & Map() const =0
int NumMyRows() const
ConstIterator begin() const
#define false
bool LinearMap() const
void BiCGSTAB(Epetra_CrsMatrix &A, Epetra_Vector &x, Epetra_Vector &b, Ifpack_CrsRiluk *M, int Maxiter, double Tolerance, double *residual, bool verbose)
double Ifpack_Condest(const Ifpack_Preconditioner &IFP, const Ifpack_CondestType CT, const int MaxIters, const double Tol, Epetra_RowMatrix *Matrix)
const Epetra_Map & DomainMap() const
bool isType(const std::string &name) const
Hypre_Solver
#define IFPACK_CHK_ERR(ifpack_err)