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;
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());
103 for (
int i=0; i!=A_->RowMatrixRowMap().NumMyElements(); ++i)
104 MyGIDsHYPRE[i] = GloballyContiguousRowMap_->GID(i);
107 IFPACK_CHK_ERRV(ColGIDsHYPRE.Import(MyGIDsHYPRE, importer, Insert, 0));
110 A_->RowMatrixColMap().NumGlobalElements(),
111 ColGIDsHYPRE.MyLength(), &ColGIDsHYPRE[0], 0, Comm()));
114 int ilower = GloballyContiguousRowMap_->MinMyGID();
115 int iupper = GloballyContiguousRowMap_->MaxMyGID();
122 XVec_ = (hypre_ParVector *) hypre_IJVectorObject(((hypre_IJVector *) XHypre_));
123 XLocal_ = hypre_ParVectorLocalVector(XVec_);
130 YVec_ = (hypre_ParVector *) hypre_IJVectorObject(((hypre_IJVector *) YHypre_));
131 YLocal_ = hypre_ParVectorLocalVector(YVec_);
135 void Ifpack_Hypre::Destroy(){
141 if(IsSolverSetup_[0]){
144 if(IsPrecondSetup_[0]){
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));
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();
185 SolverType_ = solType;
187 PrecondType_ = precType;
189 SolveOrPrec_ = chooser;
190 bool SetPrecond = list.
get(
"SetPreconditioner",
false);
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++){
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));
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));
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));
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));
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));
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));
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));
264 SolverType_ = solver;
266 PrecondType_ = solver;
272 int Ifpack_Hypre::Compute(){
273 if(IsInitialized() ==
false){
276 Time_.ResetStartTime();
279 if(SolveOrPrec_ ==
Solver){
280 IFPACK_CHK_ERR(SolverSetupPtr_(Solver_, ParMatrix_, ParX_, ParY_));
281 IsSolverSetup_[0] =
true;
283 IFPACK_CHK_ERR(PrecondSetupPtr_(Preconditioner_, ParMatrix_, ParX_, ParY_));
284 IsPrecondSetup_[0] =
true;
287 NumCompute_ = NumCompute_ + 1;
288 ComputeTime_ = ComputeTime_ + Time_.ElapsedTime();
293 int Ifpack_Hypre::CallFunctions()
const{
294 for(
int i = 0; i < NumFunsToCall_; i++){
295 IFPACK_CHK_ERR(FunsToCall_[i]->CallFunction(Solver_, Preconditioner_));
302 if(IsComputed() ==
false){
305 Time_.ResetStartTime();
306 bool SameVectors =
false;
309 if(X.Pointers() == Y.Pointers() || (NumVectors == 1 && X[0] == Y[0])){
312 for(
int VecNum = 0; VecNum <
NumVectors; VecNum++) {
321 YValues =
new double[X.MyLength()];
324 double *XTemp = XLocal_->data;
326 XLocal_->data = XValues;
327 double *YTemp = YLocal_->data;
328 YLocal_->data = YValues;
331 if(SolveOrPrec_ ==
Solver){
333 IFPACK_CHK_ERR(SolverSolvePtr_(Solver_, ParMatrix_, ParX_, ParY_));
336 IFPACK_CHK_ERR(PrecondSolvePtr_(Preconditioner_, ParMatrix_, ParX_, ParY_));
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];
346 IFPACK_CHK_ERR((*Y(VecNum)).ReplaceMyValues(NumEntries, &new_values[0], &new_indices[0]));
349 XLocal_->data = XTemp;
350 YLocal_->data = YTemp;
352 NumApplyInverse_ = NumApplyInverse_ + 1;
353 ApplyInverseTime_ = ApplyInverseTime_ + Time_.ElapsedTime();
359 if(IsInitialized() ==
false){
362 bool SameVectors =
false;
363 int NumVectors = X.NumVectors();
365 if(X.Pointers() == Y.Pointers() || (NumVectors == 1 && X[0] == Y[0])){
368 for(
int VecNum = 0; VecNum <
NumVectors; VecNum++) {
373 double *XTemp = XLocal_->data;
374 double *YTemp = YLocal_->data;
378 YValues =
new double[X.MyLength()];
380 YLocal_->data = YValues;
384 XLocal_->data = XValues;
388 IFPACK_CHK_ERR(HYPRE_ParCSRMatrixMatvecT(1.0, ParMatrix_, ParX_, 1.0, ParY_));
390 IFPACK_CHK_ERR(HYPRE_ParCSRMatrixMatvec(1.0, ParMatrix_, ParX_, 1.0, ParY_));
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];
400 IFPACK_CHK_ERR((*Y(VecNum)).ReplaceMyValues(NumEntries, new_values.data(), new_indices.data()));
403 XLocal_->data = XTemp;
404 YLocal_->data = YTemp;
410 std::ostream& Ifpack_Hypre::Print(std::ostream& os)
const{
412 if (!Comm().MyPID()) {
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;
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;
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;
439 os <<
" " << std::setw(15) << 0.0 << endl;
440 os <<
"================================================================================" << endl;
461 if(IsSolverSetup_[0]){
462 SolverDestroyPtr_(Solver_);
463 IsSolverSetup_[0] =
false;
465 SolverCreatePtr_ = &Ifpack_Hypre::Hypre_BoomerAMGCreate;
466 SolverDestroyPtr_ = &HYPRE_BoomerAMGDestroy;
467 SolverSetupPtr_ = &HYPRE_BoomerAMGSetup;
468 SolverPrecondPtr_ = NULL;
469 SolverSolvePtr_ = &HYPRE_BoomerAMGSolve;
472 if(IsSolverSetup_[0]){
473 SolverDestroyPtr_(Solver_);
474 IsSolverSetup_[0] =
false;
476 SolverCreatePtr_ = &Ifpack_Hypre::Hypre_AMSCreate;
477 SolverDestroyPtr_ = &HYPRE_AMSDestroy;
478 SolverSetupPtr_ = &HYPRE_AMSSetup;
479 SolverSolvePtr_ = &HYPRE_AMSSolve;
480 SolverPrecondPtr_ = NULL;
483 if(IsSolverSetup_[0]){
484 SolverDestroyPtr_(Solver_);
485 IsSolverSetup_[0] =
false;
487 SolverCreatePtr_ = &Ifpack_Hypre::Hypre_ParCSRHybridCreate;
488 SolverDestroyPtr_ = &HYPRE_ParCSRHybridDestroy;
489 SolverSetupPtr_ = &HYPRE_ParCSRHybridSetup;
490 SolverSolvePtr_ = &HYPRE_ParCSRHybridSolve;
491 SolverPrecondPtr_ = &HYPRE_ParCSRHybridSetPrecond;
494 if(IsSolverSetup_[0]){
495 SolverDestroyPtr_(Solver_);
496 IsSolverSetup_[0] =
false;
498 SolverCreatePtr_ = &Ifpack_Hypre::Hypre_ParCSRPCGCreate;
499 SolverDestroyPtr_ = &HYPRE_ParCSRPCGDestroy;
500 SolverSetupPtr_ = &HYPRE_ParCSRPCGSetup;
501 SolverSolvePtr_ = &HYPRE_ParCSRPCGSolve;
502 SolverPrecondPtr_ = &HYPRE_ParCSRPCGSetPrecond;
505 if(IsSolverSetup_[0]){
506 SolverDestroyPtr_(Solver_);
507 IsSolverSetup_[0] =
false;
509 SolverCreatePtr_ = &Ifpack_Hypre::Hypre_ParCSRGMRESCreate;
510 SolverDestroyPtr_ = &HYPRE_ParCSRGMRESDestroy;
511 SolverSetupPtr_ = &HYPRE_ParCSRGMRESSetup;
512 SolverPrecondPtr_ = &HYPRE_ParCSRGMRESSetPrecond;
515 if(IsSolverSetup_[0]){
516 SolverDestroyPtr_(Solver_);
517 IsSolverSetup_[0] =
false;
519 SolverCreatePtr_ = &Ifpack_Hypre::Hypre_ParCSRFlexGMRESCreate;
520 SolverDestroyPtr_ = &HYPRE_ParCSRFlexGMRESDestroy;
521 SolverSetupPtr_ = &HYPRE_ParCSRFlexGMRESSetup;
522 SolverSolvePtr_ = &HYPRE_ParCSRFlexGMRESSolve;
523 SolverPrecondPtr_ = &HYPRE_ParCSRFlexGMRESSetPrecond;
526 if(IsSolverSetup_[0]){
527 SolverDestroyPtr_(Solver_);
528 IsSolverSetup_[0] =
false;
530 SolverCreatePtr_ = &Ifpack_Hypre::Hypre_ParCSRLGMRESCreate;
531 SolverDestroyPtr_ = &HYPRE_ParCSRLGMRESDestroy;
532 SolverSetupPtr_ = &HYPRE_ParCSRLGMRESSetup;
533 SolverSolvePtr_ = &HYPRE_ParCSRLGMRESSolve;
534 SolverPrecondPtr_ = &HYPRE_ParCSRLGMRESSetPrecond;
537 if(IsSolverSetup_[0]){
538 SolverDestroyPtr_(Solver_);
539 IsSolverSetup_[0] =
false;
541 SolverCreatePtr_ = &Ifpack_Hypre::Hypre_ParCSRBiCGSTABCreate;
542 SolverDestroyPtr_ = &HYPRE_ParCSRBiCGSTABDestroy;
543 SolverSetupPtr_ = &HYPRE_ParCSRBiCGSTABSetup;
544 SolverSolvePtr_ = &HYPRE_ParCSRBiCGSTABSolve;
545 SolverPrecondPtr_ = &HYPRE_ParCSRBiCGSTABSetPrecond;
558 if(IsPrecondSetup_[0]){
559 PrecondDestroyPtr_(Preconditioner_);
560 IsPrecondSetup_[0] =
false;
562 PrecondCreatePtr_ = &Ifpack_Hypre::Hypre_BoomerAMGCreate;
563 PrecondDestroyPtr_ = &HYPRE_BoomerAMGDestroy;
564 PrecondSetupPtr_ = &HYPRE_BoomerAMGSetup;
565 PrecondSolvePtr_ = &HYPRE_BoomerAMGSolve;
568 if(IsPrecondSetup_[0]){
569 PrecondDestroyPtr_(Preconditioner_);
570 IsPrecondSetup_[0] =
false;
572 PrecondCreatePtr_ = &Ifpack_Hypre::Hypre_ParaSailsCreate;
573 PrecondDestroyPtr_ = &HYPRE_ParaSailsDestroy;
574 PrecondSetupPtr_ = &HYPRE_ParaSailsSetup;
575 PrecondSolvePtr_ = &HYPRE_ParaSailsSolve;
578 if(IsPrecondSetup_[0]){
579 PrecondDestroyPtr_(Preconditioner_);
580 IsPrecondSetup_[0] =
false;
582 PrecondCreatePtr_ = &Ifpack_Hypre::Hypre_EuclidCreate;
583 PrecondDestroyPtr_ = &HYPRE_EuclidDestroy;
584 PrecondSetupPtr_ = &HYPRE_EuclidSetup;
585 PrecondSolvePtr_ = &HYPRE_EuclidSolve;
588 if(IsPrecondSetup_[0]){
589 PrecondDestroyPtr_(Preconditioner_);
590 IsPrecondSetup_[0] =
false;
592 PrecondCreatePtr_ = &Ifpack_Hypre::Hypre_AMSCreate;
593 PrecondDestroyPtr_ = &HYPRE_AMSDestroy;
594 PrecondSetupPtr_ = &HYPRE_AMSSetup;
595 PrecondSolvePtr_ = &HYPRE_AMSSolve;
606 int Ifpack_Hypre::CreateSolver(){
608 HYPRE_ParCSRMatrixGetComm(ParMatrix_, &comm);
609 return (this->*SolverCreatePtr_)(comm, &Solver_);
613 int Ifpack_Hypre::CreatePrecond(){
615 HYPRE_ParCSRMatrixGetComm(ParMatrix_, &comm);
616 return (this->*PrecondCreatePtr_)(comm, &Preconditioner_);
621 int Ifpack_Hypre::CopyEpetraToHypre(){
622 for(
int i = 0; i < A_->NumMyRows(); i++){
625 std::vector<int> indices; indices.resize(numElements);
626 std::vector<double> values; values.resize(numElements);
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]);
633 GlobalRow[0] = GloballyContiguousRowMap_->GID(i);
634 IFPACK_CHK_ERR(HYPRE_IJMatrixSetValues(HypreA_, 1, &numEntries, GlobalRow, indices.data(), values.data()));
637 IFPACK_CHK_ERR(HYPRE_IJMatrixGetObject(HypreA_, (
void**)&ParMatrix_));
641 #endif // HAVE_HYPRE && HAVE_MPI
T & get(ParameterList &l, const std::string &name)
Ifpack_CondestType
Ifpack_CondestType: enum to define the type of condition number estimate.
#define IFPACK_CHK_ERRV(ifpack_err)
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
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)