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