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" 54 using Teuchos::rcpFromRef;
56 Ifpack_Hypre::Ifpack_Hypre(Epetra_RowMatrix* A):
59 IsInitialized_(false),
67 ApplyInverseTime_(0.0),
69 ApplyInverseFlops_(0.0),
75 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();
85 if (!A_->RowMatrixRowMap().SameAs(A_->OperatorRangeMap())) {
89 if (A_->RowMatrixRowMap().LinearMap()) {
91 GloballyContiguousRowMap_ = rcpFromRef(A_->RowMatrixRowMap());
92 GloballyContiguousColMap_ = rcpFromRef(A_->RowMatrixColMap());
98 GloballyContiguousRowMap_ = rcp(
new Epetra_Map(A_->RowMatrixRowMap().NumGlobalElements(),
99 A_->RowMatrixRowMap().NumMyElements(), 0, Comm()));
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);
106 Epetra_IntVector ColGIDsHYPRE(A_->RowMatrixColMap());
107 IFPACK_CHK_ERRV(ColGIDsHYPRE.Import(MyGIDsHYPRE, importer, Insert, 0));
109 GloballyContiguousColMap_ = rcp(
new Epetra_Map(
110 A_->RowMatrixColMap().NumGlobalElements(),
111 ColGIDsHYPRE.MyLength(), &ColGIDsHYPRE[0], 0, Comm()));
114 int ilower = GloballyContiguousRowMap_->MinMyGID();
115 int iupper = GloballyContiguousRowMap_->MaxMyGID();
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_);
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_);
135 void Ifpack_Hypre::Destroy(){
137 IFPACK_CHK_ERRV(HYPRE_IJMatrixDestroy(HypreA_));
139 IFPACK_CHK_ERRV(HYPRE_IJVectorDestroy(XHypre_));
140 IFPACK_CHK_ERRV(HYPRE_IJVectorDestroy(YHypre_));
141 if(IsSolverSetup_[0]){
142 IFPACK_CHK_ERRV(SolverDestroyPtr_(Solver_));
144 if(IsPrecondSetup_[0]){
145 IFPACK_CHK_ERRV(PrecondDestroyPtr_(Preconditioner_));
147 delete[] IsSolverSetup_;
148 delete[] IsPrecondSetup_;
152 int Ifpack_Hypre::Initialize(){
153 Time_.ResetStartTime();
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_));
166 IFPACK_CHK_ERR(SetSolverType(SolverType_));
167 IFPACK_CHK_ERR(SetPrecondType(PrecondType_));
169 if(UsePreconditioner_){
170 if(SolverPrecondPtr_ != NULL){
171 IFPACK_CHK_ERR(SolverPrecondPtr_(Solver_, PrecondSolvePtr_, PrecondSetupPtr_, Preconditioner_));
176 NumInitialize_ = NumInitialize_ + 1;
177 InitializeTime_ = InitializeTime_ + Time_.ElapsedTime();
182 int Ifpack_Hypre::SetParameters(Teuchos::ParameterList& 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);
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]));
205 int Ifpack_Hypre::AddFunToList(RCP<FunctionParameter> NewFun){
206 NumFunsToCall_ = NumFunsToCall_+1;
207 FunsToCall_.resize(NumFunsToCall_);
208 FunsToCall_[NumFunsToCall_-1] = NewFun;
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));
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));
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));
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));
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));
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));
255 int Ifpack_Hypre::SetParameter(Hypre_Chooser chooser, Hypre_Solver solver){
256 if(chooser == Solver){
257 SolverType_ = solver;
259 PrecondType_ = solver;
265 int Ifpack_Hypre::Compute(){
266 if(IsInitialized() ==
false){
267 IFPACK_CHK_ERR(Initialize());
269 Time_.ResetStartTime();
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++) {
309 IFPACK_CHK_ERR((*X(VecNum)).ExtractView(&XValues));
312 IFPACK_CHK_ERR((*Y(VecNum)).ExtractView(&YValues));
314 YValues =
new double[X.MyLength()];
317 double *XTemp = XLocal_->data;
319 XLocal_->data = XValues;
320 double *YTemp = YLocal_->data;
321 YLocal_->data = YValues;
323 IFPACK_CHK_ERR(HYPRE_ParVectorSetConstantValues(ParY_, 0.0));
324 if(SolveOrPrec_ == Solver){
326 IFPACK_CHK_ERR(SolverSolvePtr_(Solver_, ParMatrix_, ParX_, ParY_));
329 IFPACK_CHK_ERR(PrecondSolvePtr_(Preconditioner_, ParMatrix_, ParX_, ParY_));
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];
339 IFPACK_CHK_ERR((*Y(VecNum)).ReplaceMyValues(NumEntries, &new_values[0], &new_indices[0]));
342 XLocal_->data = XTemp;
343 YLocal_->data = YTemp;
345 NumApplyInverse_ = NumApplyInverse_ + 1;
346 ApplyInverseTime_ = ApplyInverseTime_ + Time_.ElapsedTime();
351 int Ifpack_Hypre::Multiply(
bool TransA,
const Epetra_MultiVector& X, Epetra_MultiVector& Y)
const{
352 if(IsInitialized() ==
false){
355 bool SameVectors =
false;
356 int NumVectors = X.NumVectors();
357 if (NumVectors != Y.NumVectors()) IFPACK_CHK_ERR(-1);
358 if(X.Pointers() == Y.Pointers()){
361 for(
int VecNum = 0; VecNum < NumVectors; VecNum++) {
365 IFPACK_CHK_ERR((*X(VecNum)).ExtractView(&XValues));
366 double *XTemp = XLocal_->data;
367 double *YTemp = YLocal_->data;
369 IFPACK_CHK_ERR((*Y(VecNum)).ExtractView(&YValues));
371 YValues =
new double[X.MyLength()];
373 YLocal_->data = YValues;
374 IFPACK_CHK_ERR(HYPRE_ParVectorSetConstantValues(ParY_,0.0));
377 XLocal_->data = XValues;
381 IFPACK_CHK_ERR(HYPRE_ParCSRMatrixMatvecT(1.0, ParMatrix_, ParX_, 1.0, ParY_));
383 IFPACK_CHK_ERR(HYPRE_ParCSRMatrixMatvec(1.0, ParMatrix_, ParX_, 1.0, ParY_));
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];
393 IFPACK_CHK_ERR((*Y(VecNum)).ReplaceMyValues(NumEntries, new_values.data(), new_indices.data()));
396 XLocal_->data = XTemp;
397 YLocal_->data = YTemp;
403 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_);
614 int Ifpack_Hypre::CopyEpetraToHypre(){
615 for(
int i = 0; i < A_->NumMyRows(); i++){
617 IFPACK_CHK_ERR(A_->NumMyRowEntries(i,numElements));
618 std::vector<int> indices; indices.resize(numElements);
619 std::vector<double> values; values.resize(numElements);
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]);
626 GlobalRow[0] = GloballyContiguousRowMap_->GID(i);
627 IFPACK_CHK_ERR(HYPRE_IJMatrixSetValues(HypreA_, 1, &numEntries, GlobalRow, indices.data(), values.data()));
629 IFPACK_CHK_ERR(HYPRE_IJMatrixAssemble(HypreA_));
630 IFPACK_CHK_ERR(HYPRE_IJMatrixGetObject(HypreA_, (
void**)&ParMatrix_));
634 #endif // HAVE_HYPRE && HAVE_MPI