IFPACK  Development
 All Classes Namespaces Files Functions Variables Enumerations Friends Pages
Ifpack_Hypre.cpp
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, int (*pt2Func)(HYPRE_Solver, int**), int** parameter){
256  RCP<FunctionParameter> temp = rcp(new FunctionParameter(chooser, pt2Func, parameter));
257  IFPACK_CHK_ERR(AddFunToList(temp));
258  return 0;
259 } //SetParameter() - int** function pointer
260 
261 //==============================================================================
262 int Ifpack_Hypre::SetParameter(Hypre_Chooser chooser, Hypre_Solver solver){
263  if(chooser == Solver){
264  SolverType_ = solver;
265  } else {
266  PrecondType_ = solver;
267  }
268  return 0;
269 } //SetParameter() - set type of solver
270 
271 //==============================================================================
272 int Ifpack_Hypre::Compute(){
273  if(IsInitialized() == false){
274  IFPACK_CHK_ERR(Initialize());
275  }
276  Time_.ResetStartTime();
277  CopyEpetraToHypre();
278  // Hypre Setup must be called after matrix has values
279  if(SolveOrPrec_ == Solver){
280  IFPACK_CHK_ERR(SolverSetupPtr_(Solver_, ParMatrix_, ParX_, ParY_));
281  IsSolverSetup_[0] = true;
282  } else {
283  IFPACK_CHK_ERR(PrecondSetupPtr_(Preconditioner_, ParMatrix_, ParX_, ParY_));
284  IsPrecondSetup_[0] = true;
285  }
286  IsComputed_ = true;
287  NumCompute_ = NumCompute_ + 1;
288  ComputeTime_ = ComputeTime_ + Time_.ElapsedTime();
289  return 0;
290 } //Compute()
291 
292 //==============================================================================
293 int Ifpack_Hypre::CallFunctions() const{
294  for(int i = 0; i < NumFunsToCall_; i++){
295  IFPACK_CHK_ERR(FunsToCall_[i]->CallFunction(Solver_, Preconditioner_));
296  }
297  return 0;
298 } //CallFunctions()
299 
300 //==============================================================================
301 int Ifpack_Hypre::ApplyInverse(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const{
302  if(IsComputed() == false){
303  IFPACK_CHK_ERR(-1);
304  }
305  Time_.ResetStartTime();
306  bool SameVectors = false;
307  int NumVectors = X.NumVectors();
308  if (NumVectors != Y.NumVectors()) IFPACK_CHK_ERR(-1); // X and Y must have same number of vectors
309  if(X.Pointers() == Y.Pointers() || (NumVectors == 1 && X[0] == Y[0])){
310  SameVectors = true;
311  }
312  for(int VecNum = 0; VecNum < NumVectors; VecNum++) {
313  //Get values for current vector in multivector.
314  // FIXME amk Nov 23, 2015: This will not work for funky data layouts
315  double * XValues;
316  IFPACK_CHK_ERR((*X(VecNum)).ExtractView(&XValues));
317  double * YValues;
318  if(!SameVectors){
319  IFPACK_CHK_ERR((*Y(VecNum)).ExtractView(&YValues));
320  } else {
321  YValues = new double[X.MyLength()];
322  }
323  // Temporarily make a pointer to data in Hypre for end
324  double *XTemp = XLocal_->data;
325  // Replace data in Hypre vectors with Epetra data
326  XLocal_->data = XValues;
327  double *YTemp = YLocal_->data;
328  YLocal_->data = YValues;
329 
330  IFPACK_CHK_ERR(HYPRE_ParVectorSetConstantValues(ParY_, 0.0));
331  if(SolveOrPrec_ == Solver){
332  // Use the solver methods
333  IFPACK_CHK_ERR(SolverSolvePtr_(Solver_, ParMatrix_, ParX_, ParY_));
334  } else {
335  // Apply the preconditioner
336  IFPACK_CHK_ERR(PrecondSolvePtr_(Preconditioner_, ParMatrix_, ParX_, ParY_));
337  }
338  if(SameVectors){
339  int NumEntries = Y.MyLength();
340  std::vector<double> new_values; new_values.resize(NumEntries);
341  std::vector<int> new_indices; new_indices.resize(NumEntries);
342  for(int i = 0; i < NumEntries; i++){
343  new_values[i] = YValues[i];
344  new_indices[i] = i;
345  }
346  IFPACK_CHK_ERR((*Y(VecNum)).ReplaceMyValues(NumEntries, &new_values[0], &new_indices[0]));
347  delete[] YValues;
348  }
349  XLocal_->data = XTemp;
350  YLocal_->data = YTemp;
351  }
352  NumApplyInverse_ = NumApplyInverse_ + 1;
353  ApplyInverseTime_ = ApplyInverseTime_ + Time_.ElapsedTime();
354  return 0;
355 } //ApplyInverse()
356 
357 //==============================================================================
358 int Ifpack_Hypre::Multiply(bool TransA, const Epetra_MultiVector& X, Epetra_MultiVector& Y) const{
359  if(IsInitialized() == false){
360  IFPACK_CHK_ERR(-1);
361  }
362  bool SameVectors = false;
363  int NumVectors = X.NumVectors();
364  if (NumVectors != Y.NumVectors()) IFPACK_CHK_ERR(-1); // X and Y must have same number of vectors
365  if(X.Pointers() == Y.Pointers() || (NumVectors == 1 && X[0] == Y[0])){
366  SameVectors = true;
367  }
368  for(int VecNum = 0; VecNum < NumVectors; VecNum++) {
369  //Get values for current vector in multivector.
370  double * XValues;
371  double * YValues;
372  IFPACK_CHK_ERR((*X(VecNum)).ExtractView(&XValues));
373  double *XTemp = XLocal_->data;
374  double *YTemp = YLocal_->data;
375  if(!SameVectors){
376  IFPACK_CHK_ERR((*Y(VecNum)).ExtractView(&YValues));
377  } else {
378  YValues = new double[X.MyLength()];
379  }
380  YLocal_->data = YValues;
381  IFPACK_CHK_ERR(HYPRE_ParVectorSetConstantValues(ParY_,0.0));
382  // Temporarily make a pointer to data in Hypre for end
383  // Replace data in Hypre vectors with epetra values
384  XLocal_->data = XValues;
385  // Do actual computation.
386  if(TransA) {
387  // Use transpose of A in multiply
388  IFPACK_CHK_ERR(HYPRE_ParCSRMatrixMatvecT(1.0, ParMatrix_, ParX_, 1.0, ParY_));
389  } else {
390  IFPACK_CHK_ERR(HYPRE_ParCSRMatrixMatvec(1.0, ParMatrix_, ParX_, 1.0, ParY_));
391  }
392  if(SameVectors){
393  int NumEntries = Y.MyLength();
394  std::vector<double> new_values; new_values.resize(NumEntries);
395  std::vector<int> new_indices; new_indices.resize(NumEntries);
396  for(int i = 0; i < NumEntries; i++){
397  new_values[i] = YValues[i];
398  new_indices[i] = i;
399  }
400  IFPACK_CHK_ERR((*Y(VecNum)).ReplaceMyValues(NumEntries, new_values.data(), new_indices.data()));
401  delete[] YValues;
402  }
403  XLocal_->data = XTemp;
404  YLocal_->data = YTemp;
405  }
406  return 0;
407 } //Multiply()
408 
409 //==============================================================================
410 std::ostream& Ifpack_Hypre::Print(std::ostream& os) const{
411  using std::endl;
412  if (!Comm().MyPID()) {
413  os << endl;
414  os << "================================================================================" << endl;
415  os << "Ifpack_Hypre: " << Label() << endl << endl;
416  os << "Using " << Comm().NumProc() << " processors." << endl;
417  os << "Global number of rows = " << A_->NumGlobalRows() << endl;
418  os << "Global number of nonzeros = " << A_->NumGlobalNonzeros() << endl;
419  os << "Condition number estimate = " << Condest() << endl;
420  os << endl;
421  os << "Phase # calls Total Time (s) Total MFlops MFlops/s" << endl;
422  os << "----- ------- -------------- ------------ --------" << endl;
423  os << "Initialize() " << std::setw(5) << NumInitialize_
424  << " " << std::setw(15) << InitializeTime_
425  << " 0.0 0.0" << endl;
426  os << "Compute() " << std::setw(5) << NumCompute_
427  << " " << std::setw(15) << ComputeTime_
428  << " " << std::setw(15) << 1.0e-6 * ComputeFlops_;
429  if (ComputeTime_ != 0.0)
430  os << " " << std::setw(15) << 1.0e-6 * ComputeFlops_ / ComputeTime_ << endl;
431  else
432  os << " " << std::setw(15) << 0.0 << endl;
433  os << "ApplyInverse() " << std::setw(5) << NumApplyInverse_
434  << " " << std::setw(15) << ApplyInverseTime_
435  << " " << std::setw(15) << 1.0e-6 * ApplyInverseFlops_;
436  if (ApplyInverseTime_ != 0.0)
437  os << " " << std::setw(15) << 1.0e-6 * ApplyInverseFlops_ / ApplyInverseTime_ << endl;
438  else
439  os << " " << std::setw(15) << 0.0 << endl;
440  os << "================================================================================" << endl;
441  os << endl;
442  }
443  return os;
444 } //Print()
445 
446 //==============================================================================
447 double Ifpack_Hypre::Condest(const Ifpack_CondestType CT,
448  const int MaxIters,
449  const double Tol,
450  Epetra_RowMatrix* Matrix_in){
451  if (!IsComputed()) // cannot compute right now
452  return(-1.0);
453  Condest_ = Ifpack_Condest(*this, CT, MaxIters, Tol, Matrix_in);
454  return(Condest_);
455 } //Condest()
456 
457 //==============================================================================
458 int Ifpack_Hypre::SetSolverType(Hypre_Solver Solver){
459  switch(Solver) {
460  case BoomerAMG:
461  if(IsSolverSetup_[0]){
462  SolverDestroyPtr_(Solver_);
463  IsSolverSetup_[0] = false;
464  }
465  SolverCreatePtr_ = &Ifpack_Hypre::Hypre_BoomerAMGCreate;
466  SolverDestroyPtr_ = &HYPRE_BoomerAMGDestroy;
467  SolverSetupPtr_ = &HYPRE_BoomerAMGSetup;
468  SolverPrecondPtr_ = NULL;
469  SolverSolvePtr_ = &HYPRE_BoomerAMGSolve;
470  break;
471  case AMS:
472  if(IsSolverSetup_[0]){
473  SolverDestroyPtr_(Solver_);
474  IsSolverSetup_[0] = false;
475  }
476  SolverCreatePtr_ = &Ifpack_Hypre::Hypre_AMSCreate;
477  SolverDestroyPtr_ = &HYPRE_AMSDestroy;
478  SolverSetupPtr_ = &HYPRE_AMSSetup;
479  SolverSolvePtr_ = &HYPRE_AMSSolve;
480  SolverPrecondPtr_ = NULL;
481  break;
482  case Hybrid:
483  if(IsSolverSetup_[0]){
484  SolverDestroyPtr_(Solver_);
485  IsSolverSetup_[0] = false;
486  }
487  SolverCreatePtr_ = &Ifpack_Hypre::Hypre_ParCSRHybridCreate;
488  SolverDestroyPtr_ = &HYPRE_ParCSRHybridDestroy;
489  SolverSetupPtr_ = &HYPRE_ParCSRHybridSetup;
490  SolverSolvePtr_ = &HYPRE_ParCSRHybridSolve;
491  SolverPrecondPtr_ = &HYPRE_ParCSRHybridSetPrecond;
492  break;
493  case PCG:
494  if(IsSolverSetup_[0]){
495  SolverDestroyPtr_(Solver_);
496  IsSolverSetup_[0] = false;
497  }
498  SolverCreatePtr_ = &Ifpack_Hypre::Hypre_ParCSRPCGCreate;
499  SolverDestroyPtr_ = &HYPRE_ParCSRPCGDestroy;
500  SolverSetupPtr_ = &HYPRE_ParCSRPCGSetup;
501  SolverSolvePtr_ = &HYPRE_ParCSRPCGSolve;
502  SolverPrecondPtr_ = &HYPRE_ParCSRPCGSetPrecond;
503  break;
504  case GMRES:
505  if(IsSolverSetup_[0]){
506  SolverDestroyPtr_(Solver_);
507  IsSolverSetup_[0] = false;
508  }
509  SolverCreatePtr_ = &Ifpack_Hypre::Hypre_ParCSRGMRESCreate;
510  SolverDestroyPtr_ = &HYPRE_ParCSRGMRESDestroy;
511  SolverSetupPtr_ = &HYPRE_ParCSRGMRESSetup;
512  SolverPrecondPtr_ = &HYPRE_ParCSRGMRESSetPrecond;
513  break;
514  case FlexGMRES:
515  if(IsSolverSetup_[0]){
516  SolverDestroyPtr_(Solver_);
517  IsSolverSetup_[0] = false;
518  }
519  SolverCreatePtr_ = &Ifpack_Hypre::Hypre_ParCSRFlexGMRESCreate;
520  SolverDestroyPtr_ = &HYPRE_ParCSRFlexGMRESDestroy;
521  SolverSetupPtr_ = &HYPRE_ParCSRFlexGMRESSetup;
522  SolverSolvePtr_ = &HYPRE_ParCSRFlexGMRESSolve;
523  SolverPrecondPtr_ = &HYPRE_ParCSRFlexGMRESSetPrecond;
524  break;
525  case LGMRES:
526  if(IsSolverSetup_[0]){
527  SolverDestroyPtr_(Solver_);
528  IsSolverSetup_[0] = false;
529  }
530  SolverCreatePtr_ = &Ifpack_Hypre::Hypre_ParCSRLGMRESCreate;
531  SolverDestroyPtr_ = &HYPRE_ParCSRLGMRESDestroy;
532  SolverSetupPtr_ = &HYPRE_ParCSRLGMRESSetup;
533  SolverSolvePtr_ = &HYPRE_ParCSRLGMRESSolve;
534  SolverPrecondPtr_ = &HYPRE_ParCSRLGMRESSetPrecond;
535  break;
536  case BiCGSTAB:
537  if(IsSolverSetup_[0]){
538  SolverDestroyPtr_(Solver_);
539  IsSolverSetup_[0] = false;
540  }
541  SolverCreatePtr_ = &Ifpack_Hypre::Hypre_ParCSRBiCGSTABCreate;
542  SolverDestroyPtr_ = &HYPRE_ParCSRBiCGSTABDestroy;
543  SolverSetupPtr_ = &HYPRE_ParCSRBiCGSTABSetup;
544  SolverSolvePtr_ = &HYPRE_ParCSRBiCGSTABSolve;
545  SolverPrecondPtr_ = &HYPRE_ParCSRBiCGSTABSetPrecond;
546  break;
547  default:
548  return -1;
549  }
550  CreateSolver();
551  return 0;
552 } //SetSolverType()
553 
554 //==============================================================================
555 int Ifpack_Hypre::SetPrecondType(Hypre_Solver Precond){
556  switch(Precond) {
557  case BoomerAMG:
558  if(IsPrecondSetup_[0]){
559  PrecondDestroyPtr_(Preconditioner_);
560  IsPrecondSetup_[0] = false;
561  }
562  PrecondCreatePtr_ = &Ifpack_Hypre::Hypre_BoomerAMGCreate;
563  PrecondDestroyPtr_ = &HYPRE_BoomerAMGDestroy;
564  PrecondSetupPtr_ = &HYPRE_BoomerAMGSetup;
565  PrecondSolvePtr_ = &HYPRE_BoomerAMGSolve;
566  break;
567  case ParaSails:
568  if(IsPrecondSetup_[0]){
569  PrecondDestroyPtr_(Preconditioner_);
570  IsPrecondSetup_[0] = false;
571  }
572  PrecondCreatePtr_ = &Ifpack_Hypre::Hypre_ParaSailsCreate;
573  PrecondDestroyPtr_ = &HYPRE_ParaSailsDestroy;
574  PrecondSetupPtr_ = &HYPRE_ParaSailsSetup;
575  PrecondSolvePtr_ = &HYPRE_ParaSailsSolve;
576  break;
577  case Euclid:
578  if(IsPrecondSetup_[0]){
579  PrecondDestroyPtr_(Preconditioner_);
580  IsPrecondSetup_[0] = false;
581  }
582  PrecondCreatePtr_ = &Ifpack_Hypre::Hypre_EuclidCreate;
583  PrecondDestroyPtr_ = &HYPRE_EuclidDestroy;
584  PrecondSetupPtr_ = &HYPRE_EuclidSetup;
585  PrecondSolvePtr_ = &HYPRE_EuclidSolve;
586  break;
587  case AMS:
588  if(IsPrecondSetup_[0]){
589  PrecondDestroyPtr_(Preconditioner_);
590  IsPrecondSetup_[0] = false;
591  }
592  PrecondCreatePtr_ = &Ifpack_Hypre::Hypre_AMSCreate;
593  PrecondDestroyPtr_ = &HYPRE_AMSDestroy;
594  PrecondSetupPtr_ = &HYPRE_AMSSetup;
595  PrecondSolvePtr_ = &HYPRE_AMSSolve;
596  break;
597  default:
598  return -1;
599  }
600  CreatePrecond();
601  return 0;
602 
603 } //SetPrecondType()
604 
605 //==============================================================================
606 int Ifpack_Hypre::CreateSolver(){
607  MPI_Comm comm;
608  HYPRE_ParCSRMatrixGetComm(ParMatrix_, &comm);
609  return (this->*SolverCreatePtr_)(comm, &Solver_);
610 } //CreateSolver()
611 
612 //==============================================================================
613 int Ifpack_Hypre::CreatePrecond(){
614  MPI_Comm comm;
615  HYPRE_ParCSRMatrixGetComm(ParMatrix_, &comm);
616  return (this->*PrecondCreatePtr_)(comm, &Preconditioner_);
617 } //CreatePrecond()
618 
619 
620 //==============================================================================
621 int Ifpack_Hypre::CopyEpetraToHypre(){
622  for(int i = 0; i < A_->NumMyRows(); i++){
623  int numElements;
624  IFPACK_CHK_ERR(A_->NumMyRowEntries(i,numElements));
625  std::vector<int> indices; indices.resize(numElements);
626  std::vector<double> values; values.resize(numElements);
627  int numEntries;
628  IFPACK_CHK_ERR(A_->ExtractMyRowCopy(i, numElements, numEntries, values.data(), indices.data()));
629  for(int j = 0; j < numEntries; j++){
630  indices[j] = GloballyContiguousColMap_->GID(indices[j]);
631  }
632  int GlobalRow[1];
633  GlobalRow[0] = GloballyContiguousRowMap_->GID(i);
634  IFPACK_CHK_ERR(HYPRE_IJMatrixSetValues(HypreA_, 1, &numEntries, GlobalRow, indices.data(), values.data()));
635  }
636  IFPACK_CHK_ERR(HYPRE_IJMatrixAssemble(HypreA_));
637  IFPACK_CHK_ERR(HYPRE_IJMatrixGetObject(HypreA_, (void**)&ParMatrix_));
638  return 0;
639 } //CopyEpetraToHypre()
640 
641 #endif // HAVE_HYPRE && HAVE_MPI
Epetra_MultiVector
Epetra_Map
Epetra_IntVector
Epetra_Import
Ifpack_Utils.h
Epetra_RowMatrix
Generated by   doxygen 1.8.5