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