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;
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_);
123 MySimpleMap_ =
rcp(
new Epetra_Map(A_->NumGlobalRows(), iupper-ilower+1, 0, Comm()));
127 void Ifpack_Hypre::Destroy(){
133 if(IsSolverSetup_[0]){
136 if(IsPrecondSetup_[0]){
139 delete[] IsSolverSetup_;
140 delete[] IsPrecondSetup_;
144 int Ifpack_Hypre::Initialize(){
145 Time_.ResetStartTime();
146 MPI_Comm comm = GetMpiComm();
147 int ilower = MySimpleMap_->MinMyGID();
148 int iupper = MySimpleMap_->MaxMyGID();
149 IFPACK_CHK_ERR(HYPRE_IJMatrixCreate(comm, ilower, iupper, ilower, iupper, &HypreA_));
150 IFPACK_CHK_ERR(HYPRE_IJMatrixSetObjectType(HypreA_, HYPRE_PARCSR));
152 for(
int i = 0; i < A_->NumMyRows(); i++){
155 std::vector<int> indices; indices.resize(numElements);
156 std::vector<double> values; values.resize(numElements);
158 IFPACK_CHK_ERR(A_->ExtractMyRowCopy(i, numElements, numEntries, &values[0], &indices[0]));
159 for(
int j = 0; j < numEntries; j++){
160 indices[j] = A_->RowMatrixColMap().GID(indices[j]);
163 GlobalRow[0] = A_->RowMatrixRowMap().GID(i);
164 IFPACK_CHK_ERR(HYPRE_IJMatrixAddToValues(HypreA_, 1, &numEntries, GlobalRow, &indices[0], &values[0]));
167 IFPACK_CHK_ERR(HYPRE_IJMatrixGetObject(HypreA_, (
void**)&ParMatrix_));
169 NumInitialize_ = NumInitialize_ + 1;
170 InitializeTime_ = InitializeTime_ + Time_.ElapsedTime();
178 SolverType_ = solType;
180 PrecondType_ = precType;
182 SolveOrPrec_ = chooser;
183 bool SetPrecond = list.
get(
"SetPreconditioner",
false);
185 int NumFunctions = list.
get(
"NumFunctions", 0);
188 if(NumFunctions > 0){
189 RCP<FunctionParameter>* params = list.
get<RCP<FunctionParameter>*>(
"Functions");
190 for(
int i = 0; i < NumFunctions; i++){
198 int Ifpack_Hypre::AddFunToList(RCP<FunctionParameter> NewFun){
199 NumFunsToCall_ = NumFunsToCall_+1;
200 FunsToCall_.resize(NumFunsToCall_);
201 FunsToCall_[NumFunsToCall_-1] = NewFun;
206 int Ifpack_Hypre::SetParameter(
Hypre_Chooser chooser,
int (*pt2Func)(HYPRE_Solver,
int),
int parameter){
207 RCP<FunctionParameter> temp =
rcp(
new FunctionParameter(chooser, pt2Func, parameter));
213 int Ifpack_Hypre::SetParameter(
Hypre_Chooser chooser,
int (*pt2Func)(HYPRE_Solver,
double),
double parameter){
214 RCP<FunctionParameter> temp =
rcp(
new FunctionParameter(chooser, pt2Func, parameter));
220 int Ifpack_Hypre::SetParameter(
Hypre_Chooser chooser,
int (*pt2Func)(HYPRE_Solver,
double,
int),
double parameter1,
int parameter2){
221 RCP<FunctionParameter> temp =
rcp(
new FunctionParameter(chooser, pt2Func, parameter1, parameter2));
227 int Ifpack_Hypre::SetParameter(
Hypre_Chooser chooser,
int (*pt2Func)(HYPRE_Solver,
int,
int),
int 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,
double*),
double* parameter){
235 RCP<FunctionParameter> temp =
rcp(
new FunctionParameter(chooser, pt2Func, parameter));
241 int Ifpack_Hypre::SetParameter(
Hypre_Chooser chooser,
int (*pt2Func)(HYPRE_Solver,
int*),
int* parameter){
242 RCP<FunctionParameter> temp =
rcp(
new FunctionParameter(chooser, pt2Func, parameter));
250 SolverType_ = solver;
252 PrecondType_ = solver;
258 int Ifpack_Hypre::Compute(){
259 if(IsInitialized() ==
false){
262 Time_.ResetStartTime();
266 if(UsePreconditioner_){
267 if(SolverPrecondPtr_ != NULL){
268 IFPACK_CHK_ERR(SolverPrecondPtr_(Solver_, PrecondSolvePtr_, PrecondSetupPtr_, Preconditioner_));
271 if(SolveOrPrec_ ==
Solver){
272 IFPACK_CHK_ERR(SolverSetupPtr_(Solver_, ParMatrix_, ParX_, ParY_));
273 IsSolverSetup_[0] =
true;
275 IFPACK_CHK_ERR(PrecondSetupPtr_(Preconditioner_, ParMatrix_, ParX_, ParY_));
276 IsPrecondSetup_[0] =
true;
279 NumCompute_ = NumCompute_ + 1;
280 ComputeTime_ = ComputeTime_ + Time_.ElapsedTime();
285 int Ifpack_Hypre::CallFunctions()
const{
286 for(
int i = 0; i < NumFunsToCall_; i++){
287 IFPACK_CHK_ERR(FunsToCall_[i]->CallFunction(Solver_, Preconditioner_));
293 int Ifpack_Hypre::ApplyInverse(
const Epetra_MultiVector& X, Epetra_MultiVector& Y)
const{
294 if(IsComputed() ==
false){
302 if(!X.Map().LinearMap() || !Y.Map().LinearMap()) {
303 std::cerr <<
"ERROR: X and Y must have contiguous maps.\n";
306 if(!X.Map().PointSameAs(*MySimpleMap_) ||
307 !Y.Map().PointSameAs(*MySimpleMap_)) {
308 std::cerr <<
"ERROR: X, Y, and A must have the same distribution.\n";
312 Time_.ResetStartTime();
313 bool SameVectors =
false;
316 if(X.Pointers() == Y.Pointers()){
319 for(
int VecNum = 0; VecNum <
NumVectors; VecNum++) {
328 YValues =
new double[X.MyLength()];
331 double *XTemp = XLocal_->data;
333 XLocal_->data = XValues;
334 double *YTemp = YLocal_->data;
335 YLocal_->data = YValues;
338 if(SolveOrPrec_ ==
Solver){
340 IFPACK_CHK_ERR(SolverSolvePtr_(Solver_, ParMatrix_, ParX_, ParY_));
343 IFPACK_CHK_ERR(PrecondSolvePtr_(Preconditioner_, ParMatrix_, ParX_, ParY_));
346 int NumEntries = Y.MyLength();
347 std::vector<double> new_values; new_values.resize(NumEntries);
348 std::vector<int> new_indices; new_indices.resize(NumEntries);
349 for(
int i = 0; i < NumEntries; i++){
350 new_values[i] = YValues[i];
353 IFPACK_CHK_ERR((*Y(VecNum)).ReplaceMyValues(NumEntries, &new_values[0], &new_indices[0]));
356 XLocal_->data = XTemp;
357 YLocal_->data = YTemp;
359 NumApplyInverse_ = NumApplyInverse_ + 1;
360 ApplyInverseTime_ = ApplyInverseTime_ + Time_.ElapsedTime();
365 int Ifpack_Hypre::Multiply(
bool TransA,
const Epetra_MultiVector& X, Epetra_MultiVector& Y)
const{
366 if(IsInitialized() ==
false){
369 bool SameVectors =
false;
372 if(X.Pointers() == Y.Pointers()){
375 for(
int VecNum = 0; VecNum <
NumVectors; VecNum++) {
380 double *XTemp = XLocal_->data;
381 double *YTemp = YLocal_->data;
385 YValues =
new double[X.MyLength()];
387 YLocal_->data = YValues;
391 XLocal_->data = XValues;
395 IFPACK_CHK_ERR(HYPRE_ParCSRMatrixMatvecT(1.0, ParMatrix_, ParX_, 1.0, ParY_));
397 IFPACK_CHK_ERR(HYPRE_ParCSRMatrixMatvec(1.0, ParMatrix_, ParX_, 1.0, ParY_));
400 int NumEntries = Y.MyLength();
401 std::vector<double> new_values; new_values.resize(NumEntries);
402 std::vector<int> new_indices; new_indices.resize(NumEntries);
403 for(
int i = 0; i < NumEntries; i++){
404 new_values[i] = YValues[i];
407 IFPACK_CHK_ERR((*Y(VecNum)).ReplaceMyValues(NumEntries, &new_values[0], &new_indices[0]));
410 XLocal_->data = XTemp;
411 YLocal_->data = YTemp;
417 std::ostream& Ifpack_Hypre::Print(std::ostream& os)
const{
420 if (!Comm().MyPID()) {
422 os <<
"================================================================================" << endl;
423 os <<
"Ifpack_Hypre: " << Label () << endl << endl;
424 os <<
"Using " << Comm().NumProc() <<
" processors." << endl;
425 os <<
"Global number of rows = " << A_->NumGlobalRows() << endl;
426 os <<
"Global number of nonzeros = " << A_->NumGlobalNonzeros() << endl;
427 os <<
"Condition number estimate = " << Condest() << endl;
429 os <<
"Phase # calls Total Time (s) Total MFlops MFlops/s" << endl;
430 os <<
"----- ------- -------------- ------------ --------" << endl;
431 os <<
"Initialize() " << std::setw(5) << NumInitialize_
432 <<
" " << std::setw(15) << InitializeTime_
433 <<
" 0.0 0.0" << endl;
434 os <<
"Compute() " << std::setw(5) << NumCompute_
435 <<
" " << std::setw(15) << ComputeTime_
436 <<
" " << std::setw(15) << 1.0e-6 * ComputeFlops_;
437 if (ComputeTime_ != 0.0)
438 os <<
" " << std::setw(15) << 1.0e-6 * ComputeFlops_ / ComputeTime_ << endl;
440 os <<
" " << std::setw(15) << 0.0 << endl;
441 os <<
"ApplyInverse() " << std::setw(5) << NumApplyInverse_
442 <<
" " << std::setw(15) << ApplyInverseTime_
443 <<
" " << std::setw(15) << 1.0e-6 * ApplyInverseFlops_;
444 if (ApplyInverseTime_ != 0.0)
445 os <<
" " << std::setw(15) << 1.0e-6 * ApplyInverseFlops_ / ApplyInverseTime_ << endl;
447 os <<
" " << std::setw(15) << 0.0 << endl;
448 os <<
"================================================================================" << endl;
458 Epetra_RowMatrix* Matrix_in){
469 if(IsSolverSetup_[0]){
470 SolverDestroyPtr_(Solver_);
471 IsSolverSetup_[0] =
false;
473 SolverCreatePtr_ = &Ifpack_Hypre::Hypre_BoomerAMGCreate;
474 SolverDestroyPtr_ = &HYPRE_BoomerAMGDestroy;
475 SolverSetupPtr_ = &HYPRE_BoomerAMGSetup;
476 SolverPrecondPtr_ = NULL;
477 SolverSolvePtr_ = &HYPRE_BoomerAMGSolve;
480 if(IsSolverSetup_[0]){
481 SolverDestroyPtr_(Solver_);
482 IsSolverSetup_[0] =
false;
484 SolverCreatePtr_ = &Ifpack_Hypre::Hypre_AMSCreate;
485 SolverDestroyPtr_ = &HYPRE_AMSDestroy;
486 SolverSetupPtr_ = &HYPRE_AMSSetup;
487 SolverSolvePtr_ = &HYPRE_AMSSolve;
488 SolverPrecondPtr_ = NULL;
491 if(IsSolverSetup_[0]){
492 SolverDestroyPtr_(Solver_);
493 IsSolverSetup_[0] =
false;
495 SolverCreatePtr_ = &Ifpack_Hypre::Hypre_ParCSRHybridCreate;
496 SolverDestroyPtr_ = &HYPRE_ParCSRHybridDestroy;
497 SolverSetupPtr_ = &HYPRE_ParCSRHybridSetup;
498 SolverSolvePtr_ = &HYPRE_ParCSRHybridSolve;
499 SolverPrecondPtr_ = &HYPRE_ParCSRHybridSetPrecond;
502 if(IsSolverSetup_[0]){
503 SolverDestroyPtr_(Solver_);
504 IsSolverSetup_[0] =
false;
506 SolverCreatePtr_ = &Ifpack_Hypre::Hypre_ParCSRPCGCreate;
507 SolverDestroyPtr_ = &HYPRE_ParCSRPCGDestroy;
508 SolverSetupPtr_ = &HYPRE_ParCSRPCGSetup;
509 SolverSolvePtr_ = &HYPRE_ParCSRPCGSolve;
510 SolverPrecondPtr_ = &HYPRE_ParCSRPCGSetPrecond;
513 if(IsSolverSetup_[0]){
514 SolverDestroyPtr_(Solver_);
515 IsSolverSetup_[0] =
false;
517 SolverCreatePtr_ = &Ifpack_Hypre::Hypre_ParCSRGMRESCreate;
518 SolverDestroyPtr_ = &HYPRE_ParCSRGMRESDestroy;
519 SolverSetupPtr_ = &HYPRE_ParCSRGMRESSetup;
520 SolverPrecondPtr_ = &HYPRE_ParCSRGMRESSetPrecond;
523 if(IsSolverSetup_[0]){
524 SolverDestroyPtr_(Solver_);
525 IsSolverSetup_[0] =
false;
527 SolverCreatePtr_ = &Ifpack_Hypre::Hypre_ParCSRFlexGMRESCreate;
528 SolverDestroyPtr_ = &HYPRE_ParCSRFlexGMRESDestroy;
529 SolverSetupPtr_ = &HYPRE_ParCSRFlexGMRESSetup;
530 SolverSolvePtr_ = &HYPRE_ParCSRFlexGMRESSolve;
531 SolverPrecondPtr_ = &HYPRE_ParCSRFlexGMRESSetPrecond;
534 if(IsSolverSetup_[0]){
535 SolverDestroyPtr_(Solver_);
536 IsSolverSetup_[0] =
false;
538 SolverCreatePtr_ = &Ifpack_Hypre::Hypre_ParCSRLGMRESCreate;
539 SolverDestroyPtr_ = &HYPRE_ParCSRLGMRESDestroy;
540 SolverSetupPtr_ = &HYPRE_ParCSRLGMRESSetup;
541 SolverSolvePtr_ = &HYPRE_ParCSRLGMRESSolve;
542 SolverPrecondPtr_ = &HYPRE_ParCSRLGMRESSetPrecond;
545 if(IsSolverSetup_[0]){
546 SolverDestroyPtr_(Solver_);
547 IsSolverSetup_[0] =
false;
549 SolverCreatePtr_ = &Ifpack_Hypre::Hypre_ParCSRBiCGSTABCreate;
550 SolverDestroyPtr_ = &HYPRE_ParCSRBiCGSTABDestroy;
551 SolverSetupPtr_ = &HYPRE_ParCSRBiCGSTABSetup;
552 SolverSolvePtr_ = &HYPRE_ParCSRBiCGSTABSolve;
553 SolverPrecondPtr_ = &HYPRE_ParCSRBiCGSTABSetPrecond;
566 if(IsPrecondSetup_[0]){
567 PrecondDestroyPtr_(Preconditioner_);
568 IsPrecondSetup_[0] =
false;
570 PrecondCreatePtr_ = &Ifpack_Hypre::Hypre_BoomerAMGCreate;
571 PrecondDestroyPtr_ = &HYPRE_BoomerAMGDestroy;
572 PrecondSetupPtr_ = &HYPRE_BoomerAMGSetup;
573 PrecondSolvePtr_ = &HYPRE_BoomerAMGSolve;
576 if(IsPrecondSetup_[0]){
577 PrecondDestroyPtr_(Preconditioner_);
578 IsPrecondSetup_[0] =
false;
580 PrecondCreatePtr_ = &Ifpack_Hypre::Hypre_ParaSailsCreate;
581 PrecondDestroyPtr_ = &HYPRE_ParaSailsDestroy;
582 PrecondSetupPtr_ = &HYPRE_ParaSailsSetup;
583 PrecondSolvePtr_ = &HYPRE_ParaSailsSolve;
586 if(IsPrecondSetup_[0]){
587 PrecondDestroyPtr_(Preconditioner_);
588 IsPrecondSetup_[0] =
false;
590 PrecondCreatePtr_ = &Ifpack_Hypre::Hypre_EuclidCreate;
591 PrecondDestroyPtr_ = &HYPRE_EuclidDestroy;
592 PrecondSetupPtr_ = &HYPRE_EuclidSetup;
593 PrecondSolvePtr_ = &HYPRE_EuclidSolve;
596 if(IsPrecondSetup_[0]){
597 PrecondDestroyPtr_(Preconditioner_);
598 IsPrecondSetup_[0] =
false;
600 PrecondCreatePtr_ = &Ifpack_Hypre::Hypre_AMSCreate;
601 PrecondDestroyPtr_ = &HYPRE_AMSDestroy;
602 PrecondSetupPtr_ = &HYPRE_AMSSetup;
603 PrecondSolvePtr_ = &HYPRE_AMSSolve;
614 int Ifpack_Hypre::CreateSolver(){
616 HYPRE_ParCSRMatrixGetComm(ParMatrix_, &comm);
617 return (this->*SolverCreatePtr_)(comm, &Solver_);
621 int Ifpack_Hypre::CreatePrecond(){
623 HYPRE_ParCSRMatrixGetComm(ParMatrix_, &comm);
624 return (this->*PrecondCreatePtr_)(comm, &Preconditioner_);
627 #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)