42 #include "Ifpack_Hypre.h" 43 #if defined(HAVE_HYPRE) && defined(HAVE_MPI) 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" 55 Ifpack_Hypre::Ifpack_Hypre(Epetra_RowMatrix* A):
58 IsInitialized_(false),
66 ApplyInverseTime_(0.0),
68 ApplyInverseFlops_(0.0),
74 UsePreconditioner_(false),
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();
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]);
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;
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_));
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_));
115 XVec_ = (hypre_ParVector *) hypre_IJVectorObject(((hypre_IJVector *) XHypre_));
116 XLocal_ = hypre_ParVectorLocalVector(XVec_);
118 YVec_ = (hypre_ParVector *) hypre_IJVectorObject(((hypre_IJVector *) YHypre_));
119 YLocal_ = hypre_ParVectorLocalVector(YVec_);
120 std::vector<int> rows; rows.resize(iupper - ilower +1);
121 for(
int i = ilower; i <= iupper; i++){
124 MySimpleMap_ = rcp(
new Epetra_Map(-1, iupper-ilower+1, &rows[0], 0, Comm()));
128 void Ifpack_Hypre::Destroy(){
130 IFPACK_CHK_ERRV(HYPRE_IJMatrixDestroy(HypreA_));
132 IFPACK_CHK_ERRV(HYPRE_IJVectorDestroy(XHypre_));
133 IFPACK_CHK_ERRV(HYPRE_IJVectorDestroy(YHypre_));
134 if(IsSolverSetup_[0]){
135 IFPACK_CHK_ERRV(SolverDestroyPtr_(Solver_));
137 if(IsPrecondSetup_[0]){
138 IFPACK_CHK_ERRV(PrecondDestroyPtr_(Preconditioner_));
140 delete[] IsSolverSetup_;
141 delete[] IsPrecondSetup_;
145 int Ifpack_Hypre::Initialize(){
146 Time_.ResetStartTime();
147 MPI_Comm comm = GetMpiComm();
148 int ilower = MySimpleMap_->MinMyGID();
149 int iupper = MySimpleMap_->MaxMyGID();
150 IFPACK_CHK_ERR(HYPRE_IJMatrixCreate(comm, ilower, iupper, ilower, iupper, &HypreA_));
151 IFPACK_CHK_ERR(HYPRE_IJMatrixSetObjectType(HypreA_, HYPRE_PARCSR));
152 IFPACK_CHK_ERR(HYPRE_IJMatrixInitialize(HypreA_));
153 for(
int i = 0; i < A_->NumMyRows(); i++){
155 IFPACK_CHK_ERR(A_->NumMyRowEntries(i,numElements));
156 std::vector<int> indices; indices.resize(numElements);
157 std::vector<double> values; values.resize(numElements);
159 IFPACK_CHK_ERR(A_->ExtractMyRowCopy(i, numElements, numEntries, &values[0], &indices[0]));
160 for(
int j = 0; j < numEntries; j++){
161 indices[j] = A_->RowMatrixColMap().GID(indices[j]);
164 GlobalRow[0] = A_->RowMatrixRowMap().GID(i);
165 IFPACK_CHK_ERR(HYPRE_IJMatrixSetValues(HypreA_, 1, &numEntries, GlobalRow, &indices[0], &values[0]));
167 IFPACK_CHK_ERR(HYPRE_IJMatrixAssemble(HypreA_));
168 IFPACK_CHK_ERR(HYPRE_IJMatrixGetObject(HypreA_, (
void**)&ParMatrix_));
170 NumInitialize_ = NumInitialize_ + 1;
171 InitializeTime_ = InitializeTime_ + Time_.ElapsedTime();
176 int Ifpack_Hypre::SetParameters(Teuchos::ParameterList& list){
178 Hypre_Solver solType = list.get(
"Solver", PCG);
179 SolverType_ = solType;
180 Hypre_Solver precType = list.get(
"Preconditioner", Euclid);
181 PrecondType_ = precType;
182 Hypre_Chooser chooser = list.get(
"SolveOrPrecondition", Solver);
183 SolveOrPrec_ = chooser;
184 bool SetPrecond = list.get(
"SetPreconditioner",
false);
185 IFPACK_CHK_ERR(SetParameter(SetPrecond));
186 int NumFunctions = list.get(
"NumFunctions", 0);
189 if(NumFunctions > 0){
190 RCP<FunctionParameter>* params = list.get<RCP<FunctionParameter>*>(
"Functions");
191 for(
int i = 0; i < NumFunctions; i++){
192 IFPACK_CHK_ERR(AddFunToList(params[i]));
199 int Ifpack_Hypre::AddFunToList(RCP<FunctionParameter> NewFun){
200 NumFunsToCall_ = NumFunsToCall_+1;
201 FunsToCall_.resize(NumFunsToCall_);
202 FunsToCall_[NumFunsToCall_-1] = NewFun;
207 int Ifpack_Hypre::SetParameter(Hypre_Chooser chooser,
int (*pt2Func)(HYPRE_Solver,
int),
int parameter){
208 RCP<FunctionParameter> temp = rcp(
new FunctionParameter(chooser, pt2Func, parameter));
209 IFPACK_CHK_ERR(AddFunToList(temp));
214 int Ifpack_Hypre::SetParameter(Hypre_Chooser chooser,
int (*pt2Func)(HYPRE_Solver,
double),
double parameter){
215 RCP<FunctionParameter> temp = rcp(
new FunctionParameter(chooser, pt2Func, parameter));
216 IFPACK_CHK_ERR(AddFunToList(temp));
221 int Ifpack_Hypre::SetParameter(Hypre_Chooser chooser,
int (*pt2Func)(HYPRE_Solver,
double,
int),
double parameter1,
int parameter2){
222 RCP<FunctionParameter> temp = rcp(
new FunctionParameter(chooser, pt2Func, parameter1, parameter2));
223 IFPACK_CHK_ERR(AddFunToList(temp));
228 int Ifpack_Hypre::SetParameter(Hypre_Chooser chooser,
int (*pt2Func)(HYPRE_Solver,
int,
int),
int parameter1,
int parameter2){
229 RCP<FunctionParameter> temp = rcp(
new FunctionParameter(chooser, pt2Func, parameter1, parameter2));
230 IFPACK_CHK_ERR(AddFunToList(temp));
235 int Ifpack_Hypre::SetParameter(Hypre_Chooser chooser,
int (*pt2Func)(HYPRE_Solver,
double*),
double* parameter){
236 RCP<FunctionParameter> temp = rcp(
new FunctionParameter(chooser, pt2Func, parameter));
237 IFPACK_CHK_ERR(AddFunToList(temp));
242 int Ifpack_Hypre::SetParameter(Hypre_Chooser chooser,
int (*pt2Func)(HYPRE_Solver,
int*),
int* parameter){
243 RCP<FunctionParameter> temp = rcp(
new FunctionParameter(chooser, pt2Func, parameter));
244 IFPACK_CHK_ERR(AddFunToList(temp));
249 int Ifpack_Hypre::SetParameter(Hypre_Chooser chooser, Hypre_Solver solver){
250 if(chooser == Solver){
251 SolverType_ = solver;
253 PrecondType_ = solver;
259 int Ifpack_Hypre::Compute(){
260 if(IsInitialized() ==
false){
261 IFPACK_CHK_ERR(Initialize());
263 Time_.ResetStartTime();
264 IFPACK_CHK_ERR(SetSolverType(SolverType_));
265 IFPACK_CHK_ERR(SetPrecondType(PrecondType_));
267 if(UsePreconditioner_){
268 if(SolverPrecondPtr_ != NULL){
269 IFPACK_CHK_ERR(SolverPrecondPtr_(Solver_, PrecondSolvePtr_, PrecondSetupPtr_, Preconditioner_));
272 if(SolveOrPrec_ == Solver){
273 IFPACK_CHK_ERR(SolverSetupPtr_(Solver_, ParMatrix_, ParX_, ParY_));
274 IsSolverSetup_[0] =
true;
276 IFPACK_CHK_ERR(PrecondSetupPtr_(Preconditioner_, ParMatrix_, ParX_, ParY_));
277 IsPrecondSetup_[0] =
true;
280 NumCompute_ = NumCompute_ + 1;
281 ComputeTime_ = ComputeTime_ + Time_.ElapsedTime();
286 int Ifpack_Hypre::CallFunctions()
const{
287 for(
int i = 0; i < NumFunsToCall_; i++){
288 IFPACK_CHK_ERR(FunsToCall_[i]->CallFunction(Solver_, Preconditioner_));
294 int Ifpack_Hypre::ApplyInverse(
const Epetra_MultiVector& X, Epetra_MultiVector& Y)
const{
295 if(IsComputed() ==
false){
298 Time_.ResetStartTime();
299 bool SameVectors =
false;
300 int NumVectors = X.NumVectors();
301 if (NumVectors != Y.NumVectors()) IFPACK_CHK_ERR(-1);
302 if(X.Pointers() == Y.Pointers()){
305 for(
int VecNum = 0; VecNum < NumVectors; VecNum++) {
308 IFPACK_CHK_ERR((*X(VecNum)).ExtractView(&XValues));
311 IFPACK_CHK_ERR((*Y(VecNum)).ExtractView(&YValues));
313 YValues =
new double[X.MyLength()];
316 double *XTemp = XLocal_->data;
318 XLocal_->data = XValues;
319 double *YTemp = YLocal_->data;
320 YLocal_->data = YValues;
322 IFPACK_CHK_ERR(HYPRE_ParVectorSetConstantValues(ParY_, 0.0));
323 if(SolveOrPrec_ == Solver){
325 IFPACK_CHK_ERR(SolverSolvePtr_(Solver_, ParMatrix_, ParX_, ParY_));
328 IFPACK_CHK_ERR(PrecondSolvePtr_(Preconditioner_, ParMatrix_, ParX_, ParY_));
331 int NumEntries = Y.MyLength();
332 std::vector<double> new_values; new_values.resize(NumEntries);
333 std::vector<int> new_indices; new_indices.resize(NumEntries);
334 for(
int i = 0; i < NumEntries; i++){
335 new_values[i] = YValues[i];
338 IFPACK_CHK_ERR((*Y(VecNum)).ReplaceMyValues(NumEntries, &new_values[0], &new_indices[0]));
341 XLocal_->data = XTemp;
342 YLocal_->data = YTemp;
344 NumApplyInverse_ = NumApplyInverse_ + 1;
345 ApplyInverseTime_ = ApplyInverseTime_ + Time_.ElapsedTime();
350 int Ifpack_Hypre::Multiply(
bool TransA,
const Epetra_MultiVector& X, Epetra_MultiVector& Y)
const{
351 if(IsInitialized() ==
false){
354 bool SameVectors =
false;
355 int NumVectors = X.NumVectors();
356 if (NumVectors != Y.NumVectors()) IFPACK_CHK_ERR(-1);
357 if(X.Pointers() == Y.Pointers()){
360 for(
int VecNum = 0; VecNum < NumVectors; VecNum++) {
364 IFPACK_CHK_ERR((*X(VecNum)).ExtractView(&XValues));
365 double *XTemp = XLocal_->data;
366 double *YTemp = YLocal_->data;
368 IFPACK_CHK_ERR((*Y(VecNum)).ExtractView(&YValues));
370 YValues =
new double[X.MyLength()];
372 YLocal_->data = YValues;
373 IFPACK_CHK_ERR(HYPRE_ParVectorSetConstantValues(ParY_,0.0));
376 XLocal_->data = XValues;
380 IFPACK_CHK_ERR(HYPRE_ParCSRMatrixMatvecT(1.0, ParMatrix_, ParX_, 1.0, ParY_));
382 IFPACK_CHK_ERR(HYPRE_ParCSRMatrixMatvec(1.0, ParMatrix_, ParX_, 1.0, ParY_));
385 int NumEntries = Y.MyLength();
386 std::vector<double> new_values; new_values.resize(NumEntries);
387 std::vector<int> new_indices; new_indices.resize(NumEntries);
388 for(
int i = 0; i < NumEntries; i++){
389 new_values[i] = YValues[i];
392 IFPACK_CHK_ERR((*Y(VecNum)).ReplaceMyValues(NumEntries, &new_values[0], &new_indices[0]));
395 XLocal_->data = XTemp;
396 YLocal_->data = YTemp;
402 std::ostream& Ifpack_Hypre::Print(std::ostream& os)
const{
405 if (!Comm().MyPID()) {
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;
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;
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;
432 os <<
" " << std::setw(15) << 0.0 << endl;
433 os <<
"================================================================================" << endl;
440 double Ifpack_Hypre::Condest(
const Ifpack_CondestType CT,
443 Epetra_RowMatrix* Matrix_in){
446 Condest_ = Ifpack_Condest(*
this, CT, MaxIters, Tol, Matrix_in);
451 int Ifpack_Hypre::SetSolverType(Hypre_Solver Solver){
454 if(IsSolverSetup_[0]){
455 SolverDestroyPtr_(Solver_);
456 IsSolverSetup_[0] =
false;
458 SolverCreatePtr_ = &Ifpack_Hypre::Hypre_BoomerAMGCreate;
459 SolverDestroyPtr_ = &HYPRE_BoomerAMGDestroy;
460 SolverSetupPtr_ = &HYPRE_BoomerAMGSetup;
461 SolverPrecondPtr_ = NULL;
462 SolverSolvePtr_ = &HYPRE_BoomerAMGSolve;
465 if(IsSolverSetup_[0]){
466 SolverDestroyPtr_(Solver_);
467 IsSolverSetup_[0] =
false;
469 SolverCreatePtr_ = &Ifpack_Hypre::Hypre_AMSCreate;
470 SolverDestroyPtr_ = &HYPRE_AMSDestroy;
471 SolverSetupPtr_ = &HYPRE_AMSSetup;
472 SolverSolvePtr_ = &HYPRE_AMSSolve;
473 SolverPrecondPtr_ = NULL;
476 if(IsSolverSetup_[0]){
477 SolverDestroyPtr_(Solver_);
478 IsSolverSetup_[0] =
false;
480 SolverCreatePtr_ = &Ifpack_Hypre::Hypre_ParCSRHybridCreate;
481 SolverDestroyPtr_ = &HYPRE_ParCSRHybridDestroy;
482 SolverSetupPtr_ = &HYPRE_ParCSRHybridSetup;
483 SolverSolvePtr_ = &HYPRE_ParCSRHybridSolve;
484 SolverPrecondPtr_ = &HYPRE_ParCSRHybridSetPrecond;
487 if(IsSolverSetup_[0]){
488 SolverDestroyPtr_(Solver_);
489 IsSolverSetup_[0] =
false;
491 SolverCreatePtr_ = &Ifpack_Hypre::Hypre_ParCSRPCGCreate;
492 SolverDestroyPtr_ = &HYPRE_ParCSRPCGDestroy;
493 SolverSetupPtr_ = &HYPRE_ParCSRPCGSetup;
494 SolverSolvePtr_ = &HYPRE_ParCSRPCGSolve;
495 SolverPrecondPtr_ = &HYPRE_ParCSRPCGSetPrecond;
498 if(IsSolverSetup_[0]){
499 SolverDestroyPtr_(Solver_);
500 IsSolverSetup_[0] =
false;
502 SolverCreatePtr_ = &Ifpack_Hypre::Hypre_ParCSRGMRESCreate;
503 SolverDestroyPtr_ = &HYPRE_ParCSRGMRESDestroy;
504 SolverSetupPtr_ = &HYPRE_ParCSRGMRESSetup;
505 SolverPrecondPtr_ = &HYPRE_ParCSRGMRESSetPrecond;
508 if(IsSolverSetup_[0]){
509 SolverDestroyPtr_(Solver_);
510 IsSolverSetup_[0] =
false;
512 SolverCreatePtr_ = &Ifpack_Hypre::Hypre_ParCSRFlexGMRESCreate;
513 SolverDestroyPtr_ = &HYPRE_ParCSRFlexGMRESDestroy;
514 SolverSetupPtr_ = &HYPRE_ParCSRFlexGMRESSetup;
515 SolverSolvePtr_ = &HYPRE_ParCSRFlexGMRESSolve;
516 SolverPrecondPtr_ = &HYPRE_ParCSRFlexGMRESSetPrecond;
519 if(IsSolverSetup_[0]){
520 SolverDestroyPtr_(Solver_);
521 IsSolverSetup_[0] =
false;
523 SolverCreatePtr_ = &Ifpack_Hypre::Hypre_ParCSRLGMRESCreate;
524 SolverDestroyPtr_ = &HYPRE_ParCSRLGMRESDestroy;
525 SolverSetupPtr_ = &HYPRE_ParCSRLGMRESSetup;
526 SolverSolvePtr_ = &HYPRE_ParCSRLGMRESSolve;
527 SolverPrecondPtr_ = &HYPRE_ParCSRLGMRESSetPrecond;
530 if(IsSolverSetup_[0]){
531 SolverDestroyPtr_(Solver_);
532 IsSolverSetup_[0] =
false;
534 SolverCreatePtr_ = &Ifpack_Hypre::Hypre_ParCSRBiCGSTABCreate;
535 SolverDestroyPtr_ = &HYPRE_ParCSRBiCGSTABDestroy;
536 SolverSetupPtr_ = &HYPRE_ParCSRBiCGSTABSetup;
537 SolverSolvePtr_ = &HYPRE_ParCSRBiCGSTABSolve;
538 SolverPrecondPtr_ = &HYPRE_ParCSRBiCGSTABSetPrecond;
548 int Ifpack_Hypre::SetPrecondType(Hypre_Solver Precond){
551 if(IsPrecondSetup_[0]){
552 PrecondDestroyPtr_(Preconditioner_);
553 IsPrecondSetup_[0] =
false;
555 PrecondCreatePtr_ = &Ifpack_Hypre::Hypre_BoomerAMGCreate;
556 PrecondDestroyPtr_ = &HYPRE_BoomerAMGDestroy;
557 PrecondSetupPtr_ = &HYPRE_BoomerAMGSetup;
558 PrecondSolvePtr_ = &HYPRE_BoomerAMGSolve;
561 if(IsPrecondSetup_[0]){
562 PrecondDestroyPtr_(Preconditioner_);
563 IsPrecondSetup_[0] =
false;
565 PrecondCreatePtr_ = &Ifpack_Hypre::Hypre_ParaSailsCreate;
566 PrecondDestroyPtr_ = &HYPRE_ParaSailsDestroy;
567 PrecondSetupPtr_ = &HYPRE_ParaSailsSetup;
568 PrecondSolvePtr_ = &HYPRE_ParaSailsSolve;
571 if(IsPrecondSetup_[0]){
572 PrecondDestroyPtr_(Preconditioner_);
573 IsPrecondSetup_[0] =
false;
575 PrecondCreatePtr_ = &Ifpack_Hypre::Hypre_EuclidCreate;
576 PrecondDestroyPtr_ = &HYPRE_EuclidDestroy;
577 PrecondSetupPtr_ = &HYPRE_EuclidSetup;
578 PrecondSolvePtr_ = &HYPRE_EuclidSolve;
581 if(IsPrecondSetup_[0]){
582 PrecondDestroyPtr_(Preconditioner_);
583 IsPrecondSetup_[0] =
false;
585 PrecondCreatePtr_ = &Ifpack_Hypre::Hypre_AMSCreate;
586 PrecondDestroyPtr_ = &HYPRE_AMSDestroy;
587 PrecondSetupPtr_ = &HYPRE_AMSSetup;
588 PrecondSolvePtr_ = &HYPRE_AMSSolve;
599 int Ifpack_Hypre::CreateSolver(){
601 HYPRE_ParCSRMatrixGetComm(ParMatrix_, &comm);
602 return (this->*SolverCreatePtr_)(comm, &Solver_);
606 int Ifpack_Hypre::CreatePrecond(){
608 HYPRE_ParCSRMatrixGetComm(ParMatrix_, &comm);
609 return (this->*PrecondCreatePtr_)(comm, &Preconditioner_);
612 #endif // HAVE_HYPRE && HAVE_MPI