Ifpack Package Browser (Single Doxygen Collection)  Development
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 
45 #include "Ifpack_Utils.h"
46 #include "Epetra_MpiComm.h"
47 #include "Epetra_IntVector.h"
48 #include "Epetra_Import.h"
49 #include "Teuchos_ParameterList.hpp"
50 #include "Teuchos_RCP.hpp"
51 
52 using Teuchos::RCP;
53 using Teuchos::rcp;
54 
55 Ifpack_Hypre::Ifpack_Hypre(Epetra_RowMatrix* A):
56  A_(rcp(A,false)),
57  UseTranspose_(false),
58  IsInitialized_(false),
59  IsComputed_(false),
60  Label_(),
61  NumInitialize_(0),
62  NumCompute_(0),
63  NumApplyInverse_(0),
64  InitializeTime_(0.0),
65  ComputeTime_(0.0),
66  ApplyInverseTime_(0.0),
67  ComputeFlops_(0.0),
68  ApplyInverseFlops_(0.0),
69  Time_(A_->Comm()),
70  SolveOrPrec_(Solver),
71  NumFunsToCall_(0),
72  SolverType_(PCG),
73  PrecondType_(Euclid),
74  UsePreconditioner_(false),
75  NiceRowMap_(true)
76 {
77  IsSolverSetup_ = new bool[1];
78  IsPrecondSetup_ = new bool[1];
79  IsSolverSetup_[0] = false;
80  IsPrecondSetup_[0] = false;
81  MPI_Comm comm = GetMpiComm();
82  int ilower = A_->RowMatrixRowMap().MinMyGID();
83  int iupper = A_->RowMatrixRowMap().MaxMyGID();
84  // Need to check if the RowMap is the way Hypre expects (if not more difficult)
85  std::vector<int> ilowers; ilowers.resize(Comm().NumProc());
86  std::vector<int> iuppers; iuppers.resize(Comm().NumProc());
87  int myLower[1]; myLower[0] = ilower;
88  int myUpper[1]; myUpper[0] = iupper;
89  Comm().GatherAll(myLower, &ilowers[0], 1);
90  Comm().GatherAll(myUpper, &iuppers[0], 1);
91  for(int i = 0; i < Comm().NumProc()-1; i++){
92  NiceRowMap_ = (NiceRowMap_ && iuppers[i]+1 == ilowers[i+1]);
93  }
94  if(!NiceRowMap_){
95  ilower = (A_->NumGlobalRows() / Comm().NumProc())*Comm().MyPID();
96  iupper = (A_->NumGlobalRows() / Comm().NumProc())*(Comm().MyPID()+1)-1;
97  if(Comm().MyPID() == Comm().NumProc()-1){
98  iupper = A_-> NumGlobalRows()-1;
99  }
100  }
101 
102  // Next create vectors that will be used when ApplyInverse() is called
103  IFPACK_CHK_ERRV(HYPRE_IJVectorCreate(comm, ilower, iupper, &XHypre_));
104  IFPACK_CHK_ERRV(HYPRE_IJVectorSetObjectType(XHypre_, HYPRE_PARCSR));
105  IFPACK_CHK_ERRV(HYPRE_IJVectorInitialize(XHypre_));
106  IFPACK_CHK_ERRV(HYPRE_IJVectorAssemble(XHypre_));
107  IFPACK_CHK_ERRV(HYPRE_IJVectorGetObject(XHypre_, (void**) &ParX_));
108 
109  IFPACK_CHK_ERRV(HYPRE_IJVectorCreate(comm, ilower, iupper, &YHypre_));
110  IFPACK_CHK_ERRV(HYPRE_IJVectorSetObjectType(YHypre_, HYPRE_PARCSR));
111  IFPACK_CHK_ERRV(HYPRE_IJVectorInitialize(YHypre_));
112  IFPACK_CHK_ERRV(HYPRE_IJVectorAssemble(YHypre_));
113  IFPACK_CHK_ERRV(HYPRE_IJVectorGetObject(YHypre_, (void**) &ParY_));
114 
115  XVec_ = (hypre_ParVector *) hypre_IJVectorObject(((hypre_IJVector *) XHypre_));
116  XLocal_ = hypre_ParVectorLocalVector(XVec_);
117 
118  YVec_ = (hypre_ParVector *) hypre_IJVectorObject(((hypre_IJVector *) YHypre_));
119  YLocal_ = hypre_ParVectorLocalVector(YVec_);
120 
121  // amk November 24, 2015: This previously created a map that Epetra does not consider
122  // to be contiguous. hypre doesn't like that, so I changed it.
123  MySimpleMap_ = rcp(new Epetra_Map(A_->NumGlobalRows(), iupper-ilower+1, 0, Comm()));
124 } //Constructor
125 
126 //==============================================================================
127 void Ifpack_Hypre::Destroy(){
128  if(IsInitialized()){
129  IFPACK_CHK_ERRV(HYPRE_IJMatrixDestroy(HypreA_));
130  }
131  IFPACK_CHK_ERRV(HYPRE_IJVectorDestroy(XHypre_));
132  IFPACK_CHK_ERRV(HYPRE_IJVectorDestroy(YHypre_));
133  if(IsSolverSetup_[0]){
134  IFPACK_CHK_ERRV(SolverDestroyPtr_(Solver_));
135  }
136  if(IsPrecondSetup_[0]){
137  IFPACK_CHK_ERRV(PrecondDestroyPtr_(Preconditioner_));
138  }
139  delete[] IsSolverSetup_;
140  delete[] IsPrecondSetup_;
141 } //Destroy()
142 
143 //==============================================================================
144 int Ifpack_Hypre::Initialize(){
145  Time_.ResetStartTime();
146  MPI_Comm comm = GetMpiComm();
147  int ilower = MySimpleMap_->MinMyGID();
148  int iupper = MySimpleMap_->MaxMyGID();
149  IFPACK_CHK_ERR(HYPRE_IJMatrixCreate(comm, ilower, iupper, ilower, iupper, &HypreA_));
150  IFPACK_CHK_ERR(HYPRE_IJMatrixSetObjectType(HypreA_, HYPRE_PARCSR));
151  IFPACK_CHK_ERR(HYPRE_IJMatrixInitialize(HypreA_));
152  for(int i = 0; i < A_->NumMyRows(); i++){
153  int numElements;
154  IFPACK_CHK_ERR(A_->NumMyRowEntries(i,numElements));
155  std::vector<int> indices; indices.resize(numElements);
156  std::vector<double> values; values.resize(numElements);
157  int numEntries;
158  IFPACK_CHK_ERR(A_->ExtractMyRowCopy(i, numElements, numEntries, &values[0], &indices[0]));
159  for(int j = 0; j < numEntries; j++){
160  indices[j] = A_->RowMatrixColMap().GID(indices[j]);
161  }
162  int GlobalRow[1];
163  GlobalRow[0] = A_->RowMatrixRowMap().GID(i);
164  IFPACK_CHK_ERR(HYPRE_IJMatrixAddToValues(HypreA_, 1, &numEntries, GlobalRow, &indices[0], &values[0]));
165  }
166  IFPACK_CHK_ERR(HYPRE_IJMatrixAssemble(HypreA_));
167  IFPACK_CHK_ERR(HYPRE_IJMatrixGetObject(HypreA_, (void**)&ParMatrix_));
168  IsInitialized_=true;
169  NumInitialize_ = NumInitialize_ + 1;
170  InitializeTime_ = InitializeTime_ + Time_.ElapsedTime();
171  return 0;
172 } //Initialize()
173 
174 //==============================================================================
175 int Ifpack_Hypre::SetParameters(Teuchos::ParameterList& list){
176  List_ = list;
177  Hypre_Solver solType = list.get("Solver", PCG);
178  SolverType_ = solType;
179  Hypre_Solver precType = list.get("Preconditioner", Euclid);
180  PrecondType_ = precType;
181  Hypre_Chooser chooser = list.get("SolveOrPrecondition", Solver);
182  SolveOrPrec_ = chooser;
183  bool SetPrecond = list.get("SetPreconditioner", false);
184  IFPACK_CHK_ERR(SetParameter(SetPrecond));
185  int NumFunctions = list.get("NumFunctions", 0);
186  FunsToCall_.clear();
187  NumFunsToCall_ = 0;
188  if(NumFunctions > 0){
189  RCP<FunctionParameter>* params = list.get<RCP<FunctionParameter>*>("Functions");
190  for(int i = 0; i < NumFunctions; i++){
191  IFPACK_CHK_ERR(AddFunToList(params[i]));
192  }
193  }
194  return 0;
195 } //SetParameters()
196 
197 //==============================================================================
198 int Ifpack_Hypre::AddFunToList(RCP<FunctionParameter> NewFun){
199  NumFunsToCall_ = NumFunsToCall_+1;
200  FunsToCall_.resize(NumFunsToCall_);
201  FunsToCall_[NumFunsToCall_-1] = NewFun;
202  return 0;
203 } //AddFunToList()
204 
205 //==============================================================================
206 int Ifpack_Hypre::SetParameter(Hypre_Chooser chooser, int (*pt2Func)(HYPRE_Solver, int), int parameter){
207  RCP<FunctionParameter> temp = rcp(new FunctionParameter(chooser, pt2Func, parameter));
208  IFPACK_CHK_ERR(AddFunToList(temp));
209  return 0;
210 } //SetParameter() - int function pointer
211 
212 //==============================================================================
213 int Ifpack_Hypre::SetParameter(Hypre_Chooser chooser, int (*pt2Func)(HYPRE_Solver, double), double parameter){
214  RCP<FunctionParameter> temp = rcp(new FunctionParameter(chooser, pt2Func, parameter));
215  IFPACK_CHK_ERR(AddFunToList(temp));
216  return 0;
217 } //SetParameter() - double function pointer
218 
219 //==============================================================================
220 int Ifpack_Hypre::SetParameter(Hypre_Chooser chooser, int (*pt2Func)(HYPRE_Solver, double, int), double parameter1, int parameter2){
221  RCP<FunctionParameter> temp = rcp(new FunctionParameter(chooser, pt2Func, parameter1, parameter2));
222  IFPACK_CHK_ERR(AddFunToList(temp));
223  return 0;
224 } //SetParameter() - double,int function pointer
225 
226 //==============================================================================
227 int Ifpack_Hypre::SetParameter(Hypre_Chooser chooser, int (*pt2Func)(HYPRE_Solver, int, int), int parameter1, int parameter2){
228  RCP<FunctionParameter> temp = rcp(new FunctionParameter(chooser, pt2Func, parameter1, parameter2));
229  IFPACK_CHK_ERR(AddFunToList(temp));
230  return 0;
231 } //SetParameter() int,int function pointer
232 
233 //==============================================================================
234 int Ifpack_Hypre::SetParameter(Hypre_Chooser chooser, int (*pt2Func)(HYPRE_Solver, double*), double* parameter){
235  RCP<FunctionParameter> temp = rcp(new FunctionParameter(chooser, pt2Func, parameter));
236  IFPACK_CHK_ERR(AddFunToList(temp));
237  return 0;
238 } //SetParameter() - double* function pointer
239 
240 //==============================================================================
241 int Ifpack_Hypre::SetParameter(Hypre_Chooser chooser, int (*pt2Func)(HYPRE_Solver, int*), int* parameter){
242  RCP<FunctionParameter> temp = rcp(new FunctionParameter(chooser, pt2Func, parameter));
243  IFPACK_CHK_ERR(AddFunToList(temp));
244  return 0;
245 } //SetParameter() - int* function pointer
246 
247 //==============================================================================
248 int Ifpack_Hypre::SetParameter(Hypre_Chooser chooser, Hypre_Solver solver){
249  if(chooser == Solver){
250  SolverType_ = solver;
251  } else {
252  PrecondType_ = solver;
253  }
254  return 0;
255 } //SetParameter() - set type of solver
256 
257 //==============================================================================
258 int Ifpack_Hypre::Compute(){
259  if(IsInitialized() == false){
260  IFPACK_CHK_ERR(Initialize());
261  }
262  Time_.ResetStartTime();
263  IFPACK_CHK_ERR(SetSolverType(SolverType_));
264  IFPACK_CHK_ERR(SetPrecondType(PrecondType_));
265  CallFunctions();
266  if(UsePreconditioner_){
267  if(SolverPrecondPtr_ != NULL){
268  IFPACK_CHK_ERR(SolverPrecondPtr_(Solver_, PrecondSolvePtr_, PrecondSetupPtr_, Preconditioner_));
269  }
270  }
271  if(SolveOrPrec_ == Solver){
272  IFPACK_CHK_ERR(SolverSetupPtr_(Solver_, ParMatrix_, ParX_, ParY_));
273  IsSolverSetup_[0] = true;
274  } else {
275  IFPACK_CHK_ERR(PrecondSetupPtr_(Preconditioner_, ParMatrix_, ParX_, ParY_));
276  IsPrecondSetup_[0] = true;
277  }
278  IsComputed_ = true;
279  NumCompute_ = NumCompute_ + 1;
280  ComputeTime_ = ComputeTime_ + Time_.ElapsedTime();
281  return 0;
282 } //Compute()
283 
284 //==============================================================================
285 int Ifpack_Hypre::CallFunctions() const{
286  for(int i = 0; i < NumFunsToCall_; i++){
287  IFPACK_CHK_ERR(FunsToCall_[i]->CallFunction(Solver_, Preconditioner_));
288  }
289  return 0;
290 } //CallFunctions()
291 
292 //==============================================================================
293 int Ifpack_Hypre::ApplyInverse(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const{
294  if(IsComputed() == false){
295  IFPACK_CHK_ERR(-1);
296  }
297  // These are hypre requirements
298  // hypre needs A, X, and Y to have the same contiguous distribution
299  // NOTE: Maps are only considered to be contiguous if they were generated using a
300  // particular constructor. Otherwise, LinearMap() will not detect whether they are
301  // actually contiguous.
302  if(!X.Map().LinearMap() || !Y.Map().LinearMap()) {
303  std::cerr << "ERROR: X and Y must have contiguous maps.\n";
304  IFPACK_CHK_ERR(-1);
305  }
306  if(!X.Map().PointSameAs(*MySimpleMap_) ||
307  !Y.Map().PointSameAs(*MySimpleMap_)) {
308  std::cerr << "ERROR: X, Y, and A must have the same distribution.\n";
309  IFPACK_CHK_ERR(-1);
310  }
311 
312  Time_.ResetStartTime();
313  bool SameVectors = false;
314  int NumVectors = X.NumVectors();
315  if (NumVectors != Y.NumVectors()) IFPACK_CHK_ERR(-1); // X and Y must have same number of vectors
316  if(X.Pointers() == Y.Pointers()){
317  SameVectors = true;
318  }
319  for(int VecNum = 0; VecNum < NumVectors; VecNum++) {
320  //Get values for current vector in multivector.
321  // FIXME amk Nov 23, 2015: This will not work for funky data layouts
322  double * XValues;
323  IFPACK_CHK_ERR((*X(VecNum)).ExtractView(&XValues));
324  double * YValues;
325  if(!SameVectors){
326  IFPACK_CHK_ERR((*Y(VecNum)).ExtractView(&YValues));
327  } else {
328  YValues = new double[X.MyLength()];
329  }
330  // Temporarily make a pointer to data in Hypre for end
331  double *XTemp = XLocal_->data;
332  // Replace data in Hypre vectors with epetra values
333  XLocal_->data = XValues;
334  double *YTemp = YLocal_->data;
335  YLocal_->data = YValues;
336 
337  IFPACK_CHK_ERR(HYPRE_ParVectorSetConstantValues(ParY_, 0.0));
338  if(SolveOrPrec_ == Solver){
339  // Use the solver methods
340  IFPACK_CHK_ERR(SolverSolvePtr_(Solver_, ParMatrix_, ParX_, ParY_));
341  } else {
342  // Apply the preconditioner
343  IFPACK_CHK_ERR(PrecondSolvePtr_(Preconditioner_, ParMatrix_, ParX_, ParY_));
344  }
345  if(SameVectors){
346  int NumEntries = Y.MyLength();
347  std::vector<double> new_values; new_values.resize(NumEntries);
348  std::vector<int> new_indices; new_indices.resize(NumEntries);
349  for(int i = 0; i < NumEntries; i++){
350  new_values[i] = YValues[i];
351  new_indices[i] = i;
352  }
353  IFPACK_CHK_ERR((*Y(VecNum)).ReplaceMyValues(NumEntries, &new_values[0], &new_indices[0]));
354  delete[] YValues;
355  }
356  XLocal_->data = XTemp;
357  YLocal_->data = YTemp;
358  }
359  NumApplyInverse_ = NumApplyInverse_ + 1;
360  ApplyInverseTime_ = ApplyInverseTime_ + Time_.ElapsedTime();
361  return 0;
362 } //ApplyInverse()
363 
364 //==============================================================================
365 int Ifpack_Hypre::Multiply(bool TransA, const Epetra_MultiVector& X, Epetra_MultiVector& Y) const{
366  if(IsInitialized() == false){
367  IFPACK_CHK_ERR(-1);
368  }
369  bool SameVectors = false;
370  int NumVectors = X.NumVectors();
371  if (NumVectors != Y.NumVectors()) IFPACK_CHK_ERR(-1); // X and Y must have same number of vectors
372  if(X.Pointers() == Y.Pointers()){
373  SameVectors = true;
374  }
375  for(int VecNum = 0; VecNum < NumVectors; VecNum++) {
376  //Get values for current vector in multivector.
377  double * XValues;
378  double * YValues;
379  IFPACK_CHK_ERR((*X(VecNum)).ExtractView(&XValues));
380  double *XTemp = XLocal_->data;
381  double *YTemp = YLocal_->data;
382  if(!SameVectors){
383  IFPACK_CHK_ERR((*Y(VecNum)).ExtractView(&YValues));
384  } else {
385  YValues = new double[X.MyLength()];
386  }
387  YLocal_->data = YValues;
388  IFPACK_CHK_ERR(HYPRE_ParVectorSetConstantValues(ParY_,0.0));
389  // Temporarily make a pointer to data in Hypre for end
390  // Replace data in Hypre vectors with epetra values
391  XLocal_->data = XValues;
392  // Do actual computation.
393  if(TransA) {
394  // Use transpose of A in multiply
395  IFPACK_CHK_ERR(HYPRE_ParCSRMatrixMatvecT(1.0, ParMatrix_, ParX_, 1.0, ParY_));
396  } else {
397  IFPACK_CHK_ERR(HYPRE_ParCSRMatrixMatvec(1.0, ParMatrix_, ParX_, 1.0, ParY_));
398  }
399  if(SameVectors){
400  int NumEntries = Y.MyLength();
401  std::vector<double> new_values; new_values.resize(NumEntries);
402  std::vector<int> new_indices; new_indices.resize(NumEntries);
403  for(int i = 0; i < NumEntries; i++){
404  new_values[i] = YValues[i];
405  new_indices[i] = i;
406  }
407  IFPACK_CHK_ERR((*Y(VecNum)).ReplaceMyValues(NumEntries, &new_values[0], &new_indices[0]));
408  delete[] YValues;
409  }
410  XLocal_->data = XTemp;
411  YLocal_->data = YTemp;
412  }
413  return 0;
414 } //Multiply()
415 
416 //==============================================================================
417 std::ostream& Ifpack_Hypre::Print(std::ostream& os) const{
418  using std::endl;
419 
420  if (!Comm().MyPID()) {
421  os << endl;
422  os << "================================================================================" << endl;
423  os << "Ifpack_Hypre: " << Label () << endl << endl;
424  os << "Using " << Comm().NumProc() << " processors." << endl;
425  os << "Global number of rows = " << A_->NumGlobalRows() << endl;
426  os << "Global number of nonzeros = " << A_->NumGlobalNonzeros() << endl;
427  os << "Condition number estimate = " << Condest() << endl;
428  os << endl;
429  os << "Phase # calls Total Time (s) Total MFlops MFlops/s" << endl;
430  os << "----- ------- -------------- ------------ --------" << endl;
431  os << "Initialize() " << std::setw(5) << NumInitialize_
432  << " " << std::setw(15) << InitializeTime_
433  << " 0.0 0.0" << endl;
434  os << "Compute() " << std::setw(5) << NumCompute_
435  << " " << std::setw(15) << ComputeTime_
436  << " " << std::setw(15) << 1.0e-6 * ComputeFlops_;
437  if (ComputeTime_ != 0.0)
438  os << " " << std::setw(15) << 1.0e-6 * ComputeFlops_ / ComputeTime_ << endl;
439  else
440  os << " " << std::setw(15) << 0.0 << endl;
441  os << "ApplyInverse() " << std::setw(5) << NumApplyInverse_
442  << " " << std::setw(15) << ApplyInverseTime_
443  << " " << std::setw(15) << 1.0e-6 * ApplyInverseFlops_;
444  if (ApplyInverseTime_ != 0.0)
445  os << " " << std::setw(15) << 1.0e-6 * ApplyInverseFlops_ / ApplyInverseTime_ << endl;
446  else
447  os << " " << std::setw(15) << 0.0 << endl;
448  os << "================================================================================" << endl;
449  os << endl;
450  }
451  return os;
452 } //Print()
453 
454 //==============================================================================
455 double Ifpack_Hypre::Condest(const Ifpack_CondestType CT,
456  const int MaxIters,
457  const double Tol,
458  Epetra_RowMatrix* Matrix_in){
459  if (!IsComputed()) // cannot compute right now
460  return(-1.0);
461  Condest_ = Ifpack_Condest(*this, CT, MaxIters, Tol, Matrix_in);
462  return(Condest_);
463 } //Condest()
464 
465 //==============================================================================
466 int Ifpack_Hypre::SetSolverType(Hypre_Solver Solver){
467  switch(Solver) {
468  case BoomerAMG:
469  if(IsSolverSetup_[0]){
470  SolverDestroyPtr_(Solver_);
471  IsSolverSetup_[0] = false;
472  }
473  SolverCreatePtr_ = &Ifpack_Hypre::Hypre_BoomerAMGCreate;
474  SolverDestroyPtr_ = &HYPRE_BoomerAMGDestroy;
475  SolverSetupPtr_ = &HYPRE_BoomerAMGSetup;
476  SolverPrecondPtr_ = NULL;
477  SolverSolvePtr_ = &HYPRE_BoomerAMGSolve;
478  break;
479  case AMS:
480  if(IsSolverSetup_[0]){
481  SolverDestroyPtr_(Solver_);
482  IsSolverSetup_[0] = false;
483  }
484  SolverCreatePtr_ = &Ifpack_Hypre::Hypre_AMSCreate;
485  SolverDestroyPtr_ = &HYPRE_AMSDestroy;
486  SolverSetupPtr_ = &HYPRE_AMSSetup;
487  SolverSolvePtr_ = &HYPRE_AMSSolve;
488  SolverPrecondPtr_ = NULL;
489  break;
490  case Hybrid:
491  if(IsSolverSetup_[0]){
492  SolverDestroyPtr_(Solver_);
493  IsSolverSetup_[0] = false;
494  }
495  SolverCreatePtr_ = &Ifpack_Hypre::Hypre_ParCSRHybridCreate;
496  SolverDestroyPtr_ = &HYPRE_ParCSRHybridDestroy;
497  SolverSetupPtr_ = &HYPRE_ParCSRHybridSetup;
498  SolverSolvePtr_ = &HYPRE_ParCSRHybridSolve;
499  SolverPrecondPtr_ = &HYPRE_ParCSRHybridSetPrecond;
500  break;
501  case PCG:
502  if(IsSolverSetup_[0]){
503  SolverDestroyPtr_(Solver_);
504  IsSolverSetup_[0] = false;
505  }
506  SolverCreatePtr_ = &Ifpack_Hypre::Hypre_ParCSRPCGCreate;
507  SolverDestroyPtr_ = &HYPRE_ParCSRPCGDestroy;
508  SolverSetupPtr_ = &HYPRE_ParCSRPCGSetup;
509  SolverSolvePtr_ = &HYPRE_ParCSRPCGSolve;
510  SolverPrecondPtr_ = &HYPRE_ParCSRPCGSetPrecond;
511  break;
512  case GMRES:
513  if(IsSolverSetup_[0]){
514  SolverDestroyPtr_(Solver_);
515  IsSolverSetup_[0] = false;
516  }
517  SolverCreatePtr_ = &Ifpack_Hypre::Hypre_ParCSRGMRESCreate;
518  SolverDestroyPtr_ = &HYPRE_ParCSRGMRESDestroy;
519  SolverSetupPtr_ = &HYPRE_ParCSRGMRESSetup;
520  SolverPrecondPtr_ = &HYPRE_ParCSRGMRESSetPrecond;
521  break;
522  case FlexGMRES:
523  if(IsSolverSetup_[0]){
524  SolverDestroyPtr_(Solver_);
525  IsSolverSetup_[0] = false;
526  }
527  SolverCreatePtr_ = &Ifpack_Hypre::Hypre_ParCSRFlexGMRESCreate;
528  SolverDestroyPtr_ = &HYPRE_ParCSRFlexGMRESDestroy;
529  SolverSetupPtr_ = &HYPRE_ParCSRFlexGMRESSetup;
530  SolverSolvePtr_ = &HYPRE_ParCSRFlexGMRESSolve;
531  SolverPrecondPtr_ = &HYPRE_ParCSRFlexGMRESSetPrecond;
532  break;
533  case LGMRES:
534  if(IsSolverSetup_[0]){
535  SolverDestroyPtr_(Solver_);
536  IsSolverSetup_[0] = false;
537  }
538  SolverCreatePtr_ = &Ifpack_Hypre::Hypre_ParCSRLGMRESCreate;
539  SolverDestroyPtr_ = &HYPRE_ParCSRLGMRESDestroy;
540  SolverSetupPtr_ = &HYPRE_ParCSRLGMRESSetup;
541  SolverSolvePtr_ = &HYPRE_ParCSRLGMRESSolve;
542  SolverPrecondPtr_ = &HYPRE_ParCSRLGMRESSetPrecond;
543  break;
544  case BiCGSTAB:
545  if(IsSolverSetup_[0]){
546  SolverDestroyPtr_(Solver_);
547  IsSolverSetup_[0] = false;
548  }
549  SolverCreatePtr_ = &Ifpack_Hypre::Hypre_ParCSRBiCGSTABCreate;
550  SolverDestroyPtr_ = &HYPRE_ParCSRBiCGSTABDestroy;
551  SolverSetupPtr_ = &HYPRE_ParCSRBiCGSTABSetup;
552  SolverSolvePtr_ = &HYPRE_ParCSRBiCGSTABSolve;
553  SolverPrecondPtr_ = &HYPRE_ParCSRBiCGSTABSetPrecond;
554  break;
555  default:
556  return -1;
557  }
558  CreateSolver();
559  return 0;
560 } //SetSolverType()
561 
562 //==============================================================================
563 int Ifpack_Hypre::SetPrecondType(Hypre_Solver Precond){
564  switch(Precond) {
565  case BoomerAMG:
566  if(IsPrecondSetup_[0]){
567  PrecondDestroyPtr_(Preconditioner_);
568  IsPrecondSetup_[0] = false;
569  }
570  PrecondCreatePtr_ = &Ifpack_Hypre::Hypre_BoomerAMGCreate;
571  PrecondDestroyPtr_ = &HYPRE_BoomerAMGDestroy;
572  PrecondSetupPtr_ = &HYPRE_BoomerAMGSetup;
573  PrecondSolvePtr_ = &HYPRE_BoomerAMGSolve;
574  break;
575  case ParaSails:
576  if(IsPrecondSetup_[0]){
577  PrecondDestroyPtr_(Preconditioner_);
578  IsPrecondSetup_[0] = false;
579  }
580  PrecondCreatePtr_ = &Ifpack_Hypre::Hypre_ParaSailsCreate;
581  PrecondDestroyPtr_ = &HYPRE_ParaSailsDestroy;
582  PrecondSetupPtr_ = &HYPRE_ParaSailsSetup;
583  PrecondSolvePtr_ = &HYPRE_ParaSailsSolve;
584  break;
585  case Euclid:
586  if(IsPrecondSetup_[0]){
587  PrecondDestroyPtr_(Preconditioner_);
588  IsPrecondSetup_[0] = false;
589  }
590  PrecondCreatePtr_ = &Ifpack_Hypre::Hypre_EuclidCreate;
591  PrecondDestroyPtr_ = &HYPRE_EuclidDestroy;
592  PrecondSetupPtr_ = &HYPRE_EuclidSetup;
593  PrecondSolvePtr_ = &HYPRE_EuclidSolve;
594  break;
595  case AMS:
596  if(IsPrecondSetup_[0]){
597  PrecondDestroyPtr_(Preconditioner_);
598  IsPrecondSetup_[0] = false;
599  }
600  PrecondCreatePtr_ = &Ifpack_Hypre::Hypre_AMSCreate;
601  PrecondDestroyPtr_ = &HYPRE_AMSDestroy;
602  PrecondSetupPtr_ = &HYPRE_AMSSetup;
603  PrecondSolvePtr_ = &HYPRE_AMSSolve;
604  break;
605  default:
606  return -1;
607  }
608  CreatePrecond();
609  return 0;
610 
611 } //SetPrecondType()
612 
613 //==============================================================================
614 int Ifpack_Hypre::CreateSolver(){
615  MPI_Comm comm;
616  HYPRE_ParCSRMatrixGetComm(ParMatrix_, &comm);
617  return (this->*SolverCreatePtr_)(comm, &Solver_);
618 } //CreateSolver()
619 
620 //==============================================================================
621 int Ifpack_Hypre::CreatePrecond(){
622  MPI_Comm comm;
623  HYPRE_ParCSRMatrixGetComm(ParMatrix_, &comm);
624  return (this->*PrecondCreatePtr_)(comm, &Preconditioner_);
625 } //CreatePrecond()
626 
627 #endif // HAVE_HYPRE && HAVE_MPI
Hybrid
T & get(ParameterList &l, const std::string &name)
ParaSails
const int NumVectors
Definition: performance.cpp:71
LGMRES
static bool Solver
Definition: performance.cpp:70
BoomerAMG
Ifpack_CondestType
Ifpack_CondestType: enum to define the type of condition number estimate.
Euclid
#define NULL
#define IFPACK_CHK_ERRV(ifpack_err)
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Hypre_Chooser
#define false
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)
AMS
Hypre_Solver
PCG
FlexGMRES
GMRES
#define IFPACK_CHK_ERR(ifpack_err)
#define true