43#if defined(HAVE_HYPRE) && defined(HAVE_MPI)
47#include "Epetra_MpiComm.h"
48#include "Epetra_IntVector.h"
49#include "Epetra_Import.h"
50#include "Teuchos_ParameterList.hpp"
51#include "Teuchos_RCP.hpp"
52#include "HYPRE_IJ_mv.h"
53#include "HYPRE_parcsr_ls.h"
55#include "_hypre_parcsr_mv.h"
56#include "_hypre_parcsr_ls.h"
57#include "_hypre_IJ_mv.h"
58#include "HYPRE_parcsr_mv.h"
63using Teuchos::rcpFromRef;
66typedef int (*int_func)(HYPRE_Solver, int);
67typedef int (*double_func)(HYPRE_Solver, double);
68typedef int (*double_int_func)(HYPRE_Solver, double, int);
69typedef int (*int_int_func)(HYPRE_Solver, int, int);
70typedef int (*int_star_func)(HYPRE_Solver,
int*);
71typedef int (*int_star_star_func)(HYPRE_Solver,
int**);
72typedef int (*double_star_func)(HYPRE_Solver,
double*);
73typedef int (*int_int_double_double_func)(HYPRE_Solver, int, int, double, double);
74typedef int (*int_int_int_double_int_int_func)(HYPRE_Solver, int, int, int, double, int, int);
75typedef int (*char_star_func)(HYPRE_Solver,
char*);
78class FunctionParameter{
81 FunctionParameter(
Hypre_Chooser chooser, int_func funct,
int param1) :
85 int_param1_(param1) {}
87 FunctionParameter(
Hypre_Chooser chooser, std::string funct_name,
int param1) :
90 int_func_(hypreMapIntFunc_.at(funct_name)),
91 int_param1_(param1) {}
94 FunctionParameter(
Hypre_Chooser chooser, double_func funct,
double param1):
98 double_param1_(param1) {}
100 FunctionParameter(
Hypre_Chooser chooser, std::string funct_name,
double param1):
103 double_func_(hypreMapDoubleFunc_.at(funct_name)),
104 double_param1_(param1) {}
107 FunctionParameter(
Hypre_Chooser chooser, double_int_func funct,
double param1,
int param2):
110 double_int_func_(funct),
112 double_param1_(param1) {}
114 FunctionParameter(
Hypre_Chooser chooser, std::string funct_name,
double param1,
int param2):
117 double_int_func_(hypreMapDoubleIntFunc_.at(funct_name)),
119 double_param1_(param1) {}
122 FunctionParameter(
Hypre_Chooser chooser, int_int_func funct,
int param1,
int param2):
125 int_int_func_(funct),
127 int_param2_(param2) {}
129 FunctionParameter(
Hypre_Chooser chooser, std::string funct_name,
int param1,
int param2):
132 int_int_func_(hypreMapIntIntFunc_.at(funct_name)),
134 int_param2_(param2) {}
137 FunctionParameter(
Hypre_Chooser chooser, int_star_func funct,
int *param1):
140 int_star_func_(funct),
141 int_star_param_(param1) {}
143 FunctionParameter(
Hypre_Chooser chooser, std::string funct_name,
int *param1):
146 int_star_func_(hypreMapIntStarFunc_.at(funct_name)),
147 int_star_param_(param1) {}
150 FunctionParameter(
Hypre_Chooser chooser, double_star_func funct,
double* param1):
153 double_star_func_(funct),
154 double_star_param_(param1) {}
156 FunctionParameter(
Hypre_Chooser chooser, std::string funct_name,
double* param1):
159 double_star_func_(hypreMapDoubleStarFunc_.at(funct_name)),
160 double_star_param_(param1) {}
163 FunctionParameter(
Hypre_Chooser chooser, int_int_double_double_func funct,
int param1,
int param2,
double param3,
double param4):
166 int_int_double_double_func_(funct),
169 double_param1_(param3),
170 double_param2_(param4) {}
172 FunctionParameter(
Hypre_Chooser chooser, std::string funct_name,
int param1,
int param2,
double param3,
double param4):
175 int_int_double_double_func_(hypreMapIntIntDoubleDoubleFunc_.at(funct_name)),
178 double_param1_(param3),
179 double_param2_(param4) {}
182 FunctionParameter(
Hypre_Chooser chooser, int_star_star_func funct,
int ** param1):
185 int_star_star_func_(funct),
186 int_star_star_param_(param1) {}
188 FunctionParameter(
Hypre_Chooser chooser, std::string funct_name,
int** param1):
191 int_star_star_func_(hypreMapIntStarStarFunc_.at(funct_name)),
192 int_star_star_param_(param1) {}
195 FunctionParameter(
Hypre_Chooser chooser, int_int_int_double_int_int_func funct,
int param1,
int param2,
int param3,
double param4,
int param5,
int param6):
198 int_int_int_double_int_int_func_(funct),
204 double_param1_(param4) {}
206 FunctionParameter(
Hypre_Chooser chooser, std::string funct_name,
int param1,
int param2,
int param3,
double param4,
int param5,
int param6):
209 int_int_int_double_int_int_func_(hypreMapIntIntIntDoubleIntIntFunc_.at(funct_name)),
215 double_param1_(param4) {}
218 FunctionParameter(
Hypre_Chooser chooser, char_star_func funct,
char *param1):
221 char_star_func_(funct),
222 char_star_param_(param1) {}
224 FunctionParameter(
Hypre_Chooser chooser, std::string funct_name,
char *param1):
227 char_star_func_(hypreMapCharStarFunc_.at(funct_name)),
228 char_star_param_(param1) {}
231 int CallFunction(HYPRE_Solver solver, HYPRE_Solver precond) {
232 if(chooser_ == Solver){
234 return int_func_(solver, int_param1_);
235 }
else if(option_ == 1){
236 return double_func_(solver, double_param1_);
237 }
else if(option_ == 2){
238 return double_int_func_(solver, double_param1_, int_param1_);
239 }
else if (option_ == 3){
240 return int_int_func_(solver, int_param1_, int_param2_);
241 }
else if (option_ == 4){
242 return int_star_func_(solver, int_star_param_);
243 }
else if (option_ == 5){
244 return double_star_func_(solver, double_star_param_);
245 }
else if (option_ == 6) {
246 return int_int_double_double_func_(solver, int_param1_, int_param2_, double_param1_, double_param2_);
247 }
else if (option_ == 7) {
248 return int_star_star_func_(solver, int_star_star_param_);
249 }
else if (option_ == 8) {
250 return int_int_int_double_int_int_func_(solver, int_param1_, int_param2_, int_param3_, double_param1_, int_param4_, int_param5_);
251 }
else if (option_ == 9) {
252 return char_star_func_(solver, char_star_param_);
258 return int_func_(precond, int_param1_);
259 }
else if(option_ == 1){
260 return double_func_(precond, double_param1_);
261 }
else if(option_ == 2){
262 return double_int_func_(precond, double_param1_, int_param1_);
263 }
else if(option_ == 3) {
264 return int_int_func_(precond, int_param1_, int_param2_);
265 }
else if(option_ == 4) {
266 return int_star_func_(precond, int_star_param_);
267 }
else if(option_ == 5) {
268 return double_star_func_(precond, double_star_param_);
269 }
else if (option_ == 6) {
270 return int_int_double_double_func_(precond, int_param1_, int_param2_, double_param1_, double_param2_);
271 }
else if (option_ == 7) {
272 return int_star_star_func_(precond, int_star_star_param_);
273 }
else if (option_ == 8) {
274 return int_int_int_double_int_int_func_(precond, int_param1_, int_param2_, int_param3_, double_param1_, int_param4_, int_param5_);
275 }
else if (option_ == 9) {
276 return char_star_func_(solver, char_star_param_);
283 static bool isFuncIntInt(std::string funct_name) {
284 return (hypreMapIntIntFunc_.find(funct_name) != hypreMapIntIntFunc_.end());
287 static bool isFuncIntIntDoubleDouble(std::string funct_name) {
288 return (hypreMapIntIntDoubleDoubleFunc_.find(funct_name) != hypreMapIntIntDoubleDoubleFunc_.end());
291 static bool isFuncIntIntIntDoubleIntInt(std::string funct_name) {
292 return (hypreMapIntIntIntDoubleIntIntFunc_.find(funct_name) != hypreMapIntIntIntDoubleIntIntFunc_.end());
295 static bool isFuncIntStarStar(std::string funct_name) {
296 return (hypreMapIntStarStarFunc_.find(funct_name) != hypreMapIntStarStarFunc_.end());
303 double_func double_func_;
304 double_int_func double_int_func_;
305 int_int_func int_int_func_;
306 int_star_func int_star_func_;
307 double_star_func double_star_func_;
308 int_int_double_double_func int_int_double_double_func_;
309 int_int_int_double_int_int_func int_int_int_double_int_int_func_;
310 int_star_star_func int_star_star_func_;
311 char_star_func char_star_func_;
317 double double_param1_;
318 double double_param2_;
319 int *int_star_param_;
320 int **int_star_star_param_;
321 double *double_star_param_;
322 char *char_star_param_;
324 static const std::map<std::string, int_func> hypreMapIntFunc_;
325 static const std::map<std::string, double_func> hypreMapDoubleFunc_;
326 static const std::map<std::string, double_int_func> hypreMapDoubleIntFunc_;
327 static const std::map<std::string, int_int_func> hypreMapIntIntFunc_;
328 static const std::map<std::string, int_star_func> hypreMapIntStarFunc_;
329 static const std::map<std::string, double_star_func> hypreMapDoubleStarFunc_;
330 static const std::map<std::string, int_int_double_double_func> hypreMapIntIntDoubleDoubleFunc_;
331 static const std::map<std::string, int_int_int_double_int_int_func> hypreMapIntIntIntDoubleIntIntFunc_;
332 static const std::map<std::string, int_star_star_func> hypreMapIntStarStarFunc_;
333 static const std::map<std::string, char_star_func> hypreMapCharStarFunc_;
338#include "Ifpack_HypreParameterMap.h"
342 UseTranspose_(
false),
343 IsInitialized_(
false),
349 InitializeTime_(0.0),
351 ApplyInverseTime_(0.0),
353 ApplyInverseFlops_(0.0),
360 IsSolverCreated_(
false),
361 IsPrecondCreated_(
false),
366 UsePreconditioner_(
false),
369 MPI_Comm comm = GetMpiComm();
373 if (!A_->RowMatrixRowMap().SameAs(A_->OperatorRangeMap())) {
377 if (A_->RowMatrixRowMap().LinearMap()) {
379 GloballyContiguousRowMap_ = rcpFromRef(A_->RowMatrixRowMap());
380 GloballyContiguousColMap_ = rcpFromRef(A_->RowMatrixColMap());
383 if(A_->OperatorDomainMap().SameAs(A_->RowMatrixRowMap())) {
384 Teuchos::RCP<const Epetra_RowMatrix> Aconst = A_;
385 GloballyContiguousColMap_ = MakeContiguousColumnMap(Aconst);
386 GloballyContiguousRowMap_ = rcp(
new Epetra_Map(A_->RowMatrixRowMap().NumGlobalElements(),
387 A_->RowMatrixRowMap().NumMyElements(), 0, Comm()));
390 throw std::runtime_error(
"Ifpack_Hypre: Unsupported map configuration: Row/Domain maps do not match");
394 int ilower = GloballyContiguousRowMap_->MinMyGID();
395 int iupper = GloballyContiguousRowMap_->MaxMyGID();
402 XVec_ = Teuchos::rcp((hypre_ParVector *) hypre_IJVectorObject(((hypre_IJVector *) XHypre_)),
false);
410 YVec_ = Teuchos::rcp((hypre_ParVector *) hypre_IJVectorObject(((hypre_IJVector *) YHypre_)),
false);
417void Ifpack_Hypre::Destroy(){
423 if(IsSolverCreated_){
426 if(IsPrecondCreated_){
446int Ifpack_Hypre::Initialize(){
447 Time_.ResetStartTime();
448 if(IsInitialized_)
return 0;
451 NumInitialize_ = NumInitialize_ + 1;
452 InitializeTime_ = InitializeTime_ + Time_.ElapsedTime();
457int Ifpack_Hypre::SetParameters(Teuchos::ParameterList& list){
459 std::map<std::string, Hypre_Solver> solverMap;
462 solverMap[
"Euclid"] =
Euclid;
463 solverMap[
"AMS"] =
AMS;
464 solverMap[
"Hybrid"] =
Hybrid;
465 solverMap[
"PCG"] =
PCG;
466 solverMap[
"GMRES"] =
GMRES;
468 solverMap[
"LGMRES"] =
LGMRES;
471 std::map<std::string, Hypre_Chooser> chooserMap;
472 chooserMap[
"Solver"] =
Solver;
477 if (list.isType<std::string>(
"hypre: Solver"))
478 solType = solverMap[list.get<std::string>(
"hypre: Solver")];
480 solType = list.get(
"hypre: Solver", PCG);
481 SolverType_ = solType;
483 if (list.isType<std::string>(
"hypre: Preconditioner"))
484 precType = solverMap[list.get<std::string>(
"hypre: Preconditioner")];
486 precType = list.get(
"hypre: Preconditioner", Euclid);
487 PrecondType_ = precType;
489 if (list.isType<std::string>(
"hypre: SolveOrPrecondition"))
490 chooser = chooserMap[list.get<std::string>(
"hypre: SolveOrPrecondition")];
492 chooser = list.get(
"hypre: SolveOrPrecondition", Solver);
493 SolveOrPrec_ = chooser;
494 bool SetPrecond = list.get(
"hypre: SetPreconditioner",
false);
496 int NumFunctions = list.get(
"hypre: NumFunctions", 0);
499 if(NumFunctions > 0){
500 RCP<FunctionParameter>* params = list.get<RCP<FunctionParameter>*>(
"hypre: Functions");
501 for(
int i = 0; i < NumFunctions; i++){
506 if (list.isSublist(
"hypre: Solver functions")) {
507 Teuchos::ParameterList solverList = list.sublist(
"hypre: Solver functions");
508 for (
auto it = solverList.begin(); it != solverList.end(); ++it) {
509 std::string funct_name = it->first;
510 if (it->second.isType<
int>()) {
511 IFPACK_CHK_ERR(AddFunToList(rcp(
new FunctionParameter(Solver, funct_name , Teuchos::getValue<int>(it->second)))));
512 }
else if (it->second.isType<
double>()) {
513 IFPACK_CHK_ERR(AddFunToList(rcp(
new FunctionParameter(Solver, funct_name , Teuchos::getValue<double>(it->second)))));
520 if (list.isSublist(
"hypre: Preconditioner functions")) {
521 Teuchos::ParameterList precList = list.sublist(
"hypre: Preconditioner functions");
522 for (
auto it = precList.begin(); it != precList.end(); ++it) {
523 std::string funct_name = it->first;
524 if (it->second.isType<
int>()) {
525 IFPACK_CHK_ERR(AddFunToList(rcp(
new FunctionParameter(Preconditioner, funct_name , Teuchos::getValue<int>(it->second)))));
526 }
else if (it->second.isType<
double>()) {
527 IFPACK_CHK_ERR(AddFunToList(rcp(
new FunctionParameter(Preconditioner, funct_name , Teuchos::getValue<double>(it->second)))));
528 }
else if (it->second.isList()) {
529 Teuchos::ParameterList pl = Teuchos::getValue<Teuchos::ParameterList>(it->second);
530 if (FunctionParameter::isFuncIntInt(funct_name)) {
531 int arg0 = pl.get<
int>(
"arg 0");
532 int arg1 = pl.get<
int>(
"arg 1");
533 IFPACK_CHK_ERR(AddFunToList(rcp(
new FunctionParameter(Preconditioner, funct_name , arg0, arg1))));
534 }
else if (FunctionParameter::isFuncIntIntDoubleDouble(funct_name)) {
535 int arg0 = pl.get<
int>(
"arg 0");
536 int arg1 = pl.get<
int>(
"arg 1");
537 double arg2 = pl.get<
double>(
"arg 2");
538 double arg3 = pl.get<
double>(
"arg 3");
539 IFPACK_CHK_ERR(AddFunToList(rcp(
new FunctionParameter(Preconditioner, funct_name , arg0, arg1, arg2, arg3))));
540 }
else if (FunctionParameter::isFuncIntIntIntDoubleIntInt(funct_name)) {
541 int arg0 = pl.get<
int>(
"arg 0");
542 int arg1 = pl.get<
int>(
"arg 1");
543 int arg2 = pl.get<
int>(
"arg 2");
544 double arg3 = pl.get<
double>(
"arg 3");
545 int arg4 = pl.get<
int>(
"arg 4");
546 int arg5 = pl.get<
int>(
"arg 5");
547 IFPACK_CHK_ERR(AddFunToList(rcp(
new FunctionParameter(Preconditioner, funct_name , arg0, arg1, arg2, arg3, arg4, arg5))));
555 if (list.isSublist(
"Coordinates") && list.sublist(
"Coordinates").isType<Teuchos::RCP<Epetra_MultiVector> >(
"Coordinates"))
556 Coords_ = list.sublist(
"Coordinates").get<Teuchos::RCP<Epetra_MultiVector> >(
"Coordinates");
557 if (list.isSublist(
"Operators") && list.sublist(
"Operators").isType<Teuchos::RCP<const Epetra_CrsMatrix> >(
"G"))
558 G_ = list.sublist(
"Operators").get<Teuchos::RCP<const Epetra_CrsMatrix> >(
"G");
560 Dump_ = list.get(
"hypre: Dump",
false);
566int Ifpack_Hypre::AddFunToList(RCP<FunctionParameter> NewFun){
567 NumFunsToCall_ = NumFunsToCall_+1;
568 FunsToCall_.resize(NumFunsToCall_);
569 FunsToCall_[NumFunsToCall_-1] = NewFun;
574int Ifpack_Hypre::SetParameter(
Hypre_Chooser chooser,
int (*pt2Func)(HYPRE_Solver,
int),
int parameter){
575 RCP<FunctionParameter> temp = rcp(
new FunctionParameter(chooser, pt2Func, parameter));
581int Ifpack_Hypre::SetParameter(
Hypre_Chooser chooser,
int (*pt2Func)(HYPRE_Solver,
double),
double parameter){
582 RCP<FunctionParameter> temp = rcp(
new FunctionParameter(chooser, pt2Func, parameter));
588int Ifpack_Hypre::SetParameter(
Hypre_Chooser chooser,
int (*pt2Func)(HYPRE_Solver,
double,
int),
double parameter1,
int parameter2){
589 RCP<FunctionParameter> temp = rcp(
new FunctionParameter(chooser, pt2Func, parameter1, parameter2));
595int Ifpack_Hypre::SetParameter(
Hypre_Chooser chooser,
int (*pt2Func)(HYPRE_Solver,
int,
int),
int parameter1,
int parameter2){
596 RCP<FunctionParameter> temp = rcp(
new FunctionParameter(chooser, pt2Func, parameter1, parameter2));
602int Ifpack_Hypre::SetParameter(
Hypre_Chooser chooser,
int (*pt2Func)(HYPRE_Solver,
double*),
double* parameter){
603 RCP<FunctionParameter> temp = rcp(
new FunctionParameter(chooser, pt2Func, parameter));
609int Ifpack_Hypre::SetParameter(
Hypre_Chooser chooser,
int (*pt2Func)(HYPRE_Solver,
int*),
int* parameter){
610 RCP<FunctionParameter> temp = rcp(
new FunctionParameter(chooser, pt2Func, parameter));
616int Ifpack_Hypre::SetParameter(
Hypre_Chooser chooser,
int (*pt2Func)(HYPRE_Solver,
int**),
int** parameter){
617 RCP<FunctionParameter> temp = rcp(
new FunctionParameter(chooser, pt2Func, parameter));
624 if(chooser == Solver){
625 SolverType_ = solver;
627 PrecondType_ = solver;
633int Ifpack_Hypre::SetDiscreteGradient(Teuchos::RCP<const Epetra_CrsMatrix> G){
636 if(!A_->RowMatrixRowMap().SameAs(G->RowMap()))
637 throw std::runtime_error(
"Ifpack_Hypre: Edge map mismatch: A and discrete gradient");
640 GloballyContiguousNodeRowMap_ = rcp(
new Epetra_Map(G->DomainMap().NumGlobalElements(),
641 G->DomainMap().NumMyElements(), 0, Comm()));
642 Teuchos::RCP<const Epetra_RowMatrix> Grow = Teuchos::rcp_dynamic_cast<const Epetra_RowMatrix>(G);
643 GloballyContiguousNodeColMap_ = MakeContiguousColumnMap(Grow);
646 MPI_Comm comm = GetMpiComm();
647 int ilower = GloballyContiguousRowMap_->MinMyGID();
648 int iupper = GloballyContiguousRowMap_->MaxMyGID();
649 int jlower = GloballyContiguousNodeRowMap_->MinMyGID();
650 int jupper = GloballyContiguousNodeRowMap_->MaxMyGID();
651 IFPACK_CHK_ERR(HYPRE_IJMatrixCreate(comm, ilower, iupper, jlower, jupper, &HypreG_));
652 IFPACK_CHK_ERR(HYPRE_IJMatrixSetObjectType(HypreG_, HYPRE_PARCSR));
655 std::vector<int> new_indices(G->MaxNumEntries());
656 for(
int i = 0; i < G->NumMyRows(); i++){
660 IFPACK_CHK_ERR(G->ExtractMyRowView(i, numEntries, values, indices));
661 for(
int j = 0; j < numEntries; j++){
662 new_indices[j] = GloballyContiguousNodeColMap_->GID(indices[j]);
665 GlobalRow[0] = GloballyContiguousRowMap_->GID(i);
666 IFPACK_CHK_ERR(HYPRE_IJMatrixSetValues(HypreG_, 1, &numEntries, GlobalRow, new_indices.data(), values));
669 IFPACK_CHK_ERR(HYPRE_IJMatrixGetObject(HypreG_, (
void**)&ParMatrixG_));
672 HYPRE_ParCSRMatrixPrint(ParMatrixG_,
"G.mat");
674 if(SolverType_ == AMS)
675 HYPRE_AMSSetDiscreteGradient(Solver_, ParMatrixG_);
676 if(PrecondType_ == AMS)
677 HYPRE_AMSSetDiscreteGradient(Preconditioner_, ParMatrixG_);
683int Ifpack_Hypre::SetCoordinates(Teuchos::RCP<Epetra_MultiVector> coords) {
685 if(!G_.is_null() && !G_->DomainMap().SameAs(coords->Map()))
686 throw std::runtime_error(
"Ifpack_Hypre: Node map mismatch: G->DomainMap() and coords");
688 if(SolverType_ != AMS && PrecondType_ != AMS)
699 MPI_Comm comm = GetMpiComm();
700 int NumEntries = coords->MyLength();
701 int * indices = GloballyContiguousNodeRowMap_->MyGlobalElements();
703 int ilower = GloballyContiguousNodeRowMap_->MinMyGID();
704 int iupper = GloballyContiguousNodeRowMap_->MaxMyGID();
706 if( NumEntries != iupper-ilower+1) {
707 std::cout<<
"Ifpack_Hypre::SetCoordinates(): Error on rank "<<Comm().MyPID()<<
": MyLength = "<<coords->MyLength()<<
" GID range = ["<<ilower<<
","<<iupper<<
"]"<<std::endl;
708 throw std::runtime_error(
"Ifpack_Hypre: SetCoordinates: Length mismatch");
711 IFPACK_CHK_ERR(HYPRE_IJVectorCreate(comm, ilower, iupper, &xHypre_));
712 IFPACK_CHK_ERR(HYPRE_IJVectorSetObjectType(xHypre_, HYPRE_PARCSR));
714 IFPACK_CHK_ERR(HYPRE_IJVectorSetValues(xHypre_,NumEntries,indices,xPtr));
716 IFPACK_CHK_ERR(HYPRE_IJVectorGetObject(xHypre_, (
void**) &xPar_));
718 IFPACK_CHK_ERR(HYPRE_IJVectorCreate(comm, ilower, iupper, &yHypre_));
719 IFPACK_CHK_ERR(HYPRE_IJVectorSetObjectType(yHypre_, HYPRE_PARCSR));
721 IFPACK_CHK_ERR(HYPRE_IJVectorSetValues(yHypre_,NumEntries,indices,yPtr));
723 IFPACK_CHK_ERR(HYPRE_IJVectorGetObject(yHypre_, (
void**) &yPar_));
725 IFPACK_CHK_ERR(HYPRE_IJVectorCreate(comm, ilower, iupper, &zHypre_));
726 IFPACK_CHK_ERR(HYPRE_IJVectorSetObjectType(zHypre_, HYPRE_PARCSR));
728 IFPACK_CHK_ERR(HYPRE_IJVectorSetValues(zHypre_,NumEntries,indices,zPtr));
730 IFPACK_CHK_ERR(HYPRE_IJVectorGetObject(zHypre_, (
void**) &zPar_));
733 HYPRE_ParVectorPrint(xPar_,
"coordX.dat");
734 HYPRE_ParVectorPrint(yPar_,
"coordY.dat");
735 HYPRE_ParVectorPrint(zPar_,
"coordZ.dat");
738 if(SolverType_ == AMS)
739 HYPRE_AMSSetCoordinateVectors(Solver_, xPar_, yPar_, zPar_);
740 if(PrecondType_ == AMS)
741 HYPRE_AMSSetCoordinateVectors(Preconditioner_, xPar_, yPar_, zPar_);
748int Ifpack_Hypre::Compute(){
749 if(IsInitialized() ==
false){
752 Time_.ResetStartTime();
759 MPI_Comm comm = GetMpiComm();
760 int ilower = GloballyContiguousRowMap_->MinMyGID();
761 int iupper = GloballyContiguousRowMap_->MaxMyGID();
762 IFPACK_CHK_ERR(HYPRE_IJMatrixCreate(comm, ilower, iupper, ilower, iupper, &HypreA_));
763 IFPACK_CHK_ERR(HYPRE_IJMatrixSetObjectType(HypreA_, HYPRE_PARCSR));
766 if(SolveOrPrec_ == Solver) {
768 if (SolverPrecondPtr_ != NULL && UsePreconditioner_) {
772 IFPACK_CHK_ERR(SolverPrecondPtr_(Solver_, PrecondSolvePtr_, PrecondSetupPtr_, Preconditioner_));
782 SetDiscreteGradient(G_);
785 if (!Coords_.is_null()) {
786 SetCoordinates(Coords_);
790 if(SolveOrPrec_ == Solver){
791 IFPACK_CHK_ERR(SolverSetupPtr_(Solver_, ParMatrix_, ParX_, ParY_));
793 IFPACK_CHK_ERR(PrecondSetupPtr_(Preconditioner_, ParMatrix_, ParX_, ParY_));
797 if(Dump_ && PrecondSolvePtr_ == &HYPRE_BoomerAMGSolve) {
798 hypre_ParAMGData *amg_data = (hypre_ParAMGData*) Preconditioner_;
799 hypre_ParCSRMatrix **A_array = hypre_ParAMGDataAArray(amg_data);
800 hypre_ParCSRMatrix **P_array = hypre_ParAMGDataPArray(amg_data);
802 HYPRE_Int **CF_marker_array = hypre_ParAMGDataCFMarkerArray(amg_data);
804 HYPRE_Int num_levels = hypre_ParAMGDataNumLevels(amg_data);
807 for(
int k=0; k<num_levels; k++) {
809 sprintf(ofs,
"A_matrix.bmg.%d.dat",k);
810 HYPRE_ParCSRMatrixPrint(A_array[k], ofs);
811 if(k!=num_levels-1) {
813 sprintf(ofs,
"P_matrix.bmg.%d.dat",k);
814 HYPRE_ParCSRMatrixPrint(P_array[k], ofs);
819 HYPRE_Int local_size = hypre_CSRMatrixNumRows(hypre_ParCSRMatrixDiag(A_array[k]));
820 sprintf(ofs,
"cf_marker.bmg.%d.dat.%d",k,Comm().MyPID());
821 FILE * f = fopen(ofs,
"w");
822 fprintf(f,
"%%%%MatrixMarket matrix array real general\n");
823 fprintf(f,
"%d 1\n",local_size);
824 for(
int i=0; i<local_size; i++)
825 fprintf(f,
"%d\n",(
int)CF_marker_array[k][i]);
835 NumCompute_ = NumCompute_ + 1;
836 ComputeTime_ = ComputeTime_ + Time_.ElapsedTime();
841int Ifpack_Hypre::CallFunctions()
const{
842 for(
int i = 0; i < NumFunsToCall_; i++){
843 IFPACK_CHK_ERR(FunsToCall_[i]->CallFunction(Solver_, Preconditioner_));
850 if(IsComputed() ==
false){
853 Time_.ResetStartTime();
854 hypre_Vector *XLocal_ = hypre_ParVectorLocalVector(XVec_);
855 hypre_Vector *YLocal_ = hypre_ParVectorLocalVector(YVec_);
857 bool SameVectors =
false;
868 for(
int VecNum = 0; VecNum <
NumVectors; VecNum++) {
871 double * XValues =
const_cast<double*
>(X[VecNum]);
874 YValues =
const_cast<double*
>(Y[VecNum]);
876 YValues = VectorCache_.getRawPtr();
879 double *XTemp = XLocal_->data;
881 XLocal_->data = XValues;
882 double *YTemp = YLocal_->data;
883 YLocal_->data = YValues;
886 if(SolveOrPrec_ == Solver){
888 IFPACK_CHK_ERR(SolverSolvePtr_(Solver_, ParMatrix_, ParX_, ParY_));
891 IFPACK_CHK_ERR(PrecondSolvePtr_(Preconditioner_, ParMatrix_, ParX_, ParY_));
895 for(
int i = 0; i < NumEntries; i++)
896 Y[VecNum][i] = YValues[i];
898 XLocal_->data = XTemp;
899 YLocal_->data = YTemp;
901 NumApplyInverse_ = NumApplyInverse_ + 1;
902 ApplyInverseTime_ = ApplyInverseTime_ + Time_.ElapsedTime();
908 if(IsInitialized() ==
false){
911 hypre_Vector *XLocal_ = hypre_ParVectorLocalVector(XVec_);
912 hypre_Vector *YLocal_ = hypre_ParVectorLocalVector(YVec_);
913 bool SameVectors =
false;
924 for(
int VecNum = 0; VecNum <
NumVectors; VecNum++) {
926 double * XValues=
const_cast<double*
>(X[VecNum]);
928 double *XTemp = XLocal_->data;
929 double *YTemp = YLocal_->data;
931 YValues =
const_cast<double*
>(Y[VecNum]);
933 YValues = VectorCache_.getRawPtr();
935 YLocal_->data = YValues;
939 XLocal_->data = XValues;
943 IFPACK_CHK_ERR(HYPRE_ParCSRMatrixMatvecT(1.0, ParMatrix_, ParX_, 1.0, ParY_));
945 IFPACK_CHK_ERR(HYPRE_ParCSRMatrixMatvec(1.0, ParMatrix_, ParX_, 1.0, ParY_));
949 for(
int i = 0; i < NumEntries; i++)
950 Y[VecNum][i] = YValues[i];
952 XLocal_->data = XTemp;
953 YLocal_->data = YTemp;
959std::ostream& Ifpack_Hypre::Print(std::ostream& os)
const{
961 if (!Comm().MyPID()) {
963 os <<
"================================================================================" << endl;
964 os <<
"Ifpack_Hypre: " << Label() << endl << endl;
965 os <<
"Using " << Comm().NumProc() <<
" processors." << endl;
966 os <<
"Global number of rows = " << A_->NumGlobalRows() << endl;
967 os <<
"Global number of nonzeros = " << A_->NumGlobalNonzeros() << endl;
968 os <<
"Condition number estimate = " << Condest() << endl;
970 os <<
"Phase # calls Total Time (s) Total MFlops MFlops/s" << endl;
971 os <<
"----- ------- -------------- ------------ --------" << endl;
972 os <<
"Initialize() " << std::setw(5) << NumInitialize_
973 <<
" " << std::setw(15) << InitializeTime_
974 <<
" 0.0 0.0" << endl;
975 os <<
"Compute() " << std::setw(5) << NumCompute_
976 <<
" " << std::setw(15) << ComputeTime_
977 <<
" " << std::setw(15) << 1.0e-6 * ComputeFlops_;
978 if (ComputeTime_ != 0.0)
979 os <<
" " << std::setw(15) << 1.0e-6 * ComputeFlops_ / ComputeTime_ << endl;
981 os <<
" " << std::setw(15) << 0.0 << endl;
982 os <<
"ApplyInverse() " << std::setw(5) << NumApplyInverse_
983 <<
" " << std::setw(15) << ApplyInverseTime_
984 <<
" " << std::setw(15) << 1.0e-6 * ApplyInverseFlops_;
985 if (ApplyInverseTime_ != 0.0)
986 os <<
" " << std::setw(15) << 1.0e-6 * ApplyInverseFlops_ / ApplyInverseTime_ << endl;
988 os <<
" " << std::setw(15) << 0.0 << endl;
989 os <<
"================================================================================" << endl;
1010 if(IsSolverCreated_){
1011 SolverDestroyPtr_(Solver_);
1012 IsSolverCreated_ =
false;
1014 SolverCreatePtr_ = &Ifpack_Hypre::Hypre_BoomerAMGCreate;
1015 SolverDestroyPtr_ = &HYPRE_BoomerAMGDestroy;
1016 SolverSetupPtr_ = &HYPRE_BoomerAMGSetup;
1017 SolverPrecondPtr_ = NULL;
1018 SolverSolvePtr_ = &HYPRE_BoomerAMGSolve;
1021 if(IsSolverCreated_){
1022 SolverDestroyPtr_(Solver_);
1023 IsSolverCreated_ =
false;
1025 SolverCreatePtr_ = &Ifpack_Hypre::Hypre_AMSCreate;
1026 SolverDestroyPtr_ = &HYPRE_AMSDestroy;
1027 SolverSetupPtr_ = &HYPRE_AMSSetup;
1028 SolverSolvePtr_ = &HYPRE_AMSSolve;
1029 SolverPrecondPtr_ = NULL;
1032 if(IsSolverCreated_){
1033 SolverDestroyPtr_(Solver_);
1034 IsSolverCreated_ =
false;
1036 SolverCreatePtr_ = &Ifpack_Hypre::Hypre_ParCSRHybridCreate;
1037 SolverDestroyPtr_ = &HYPRE_ParCSRHybridDestroy;
1038 SolverSetupPtr_ = &HYPRE_ParCSRHybridSetup;
1039 SolverSolvePtr_ = &HYPRE_ParCSRHybridSolve;
1040 SolverPrecondPtr_ = &HYPRE_ParCSRHybridSetPrecond;
1043 if(IsSolverCreated_){
1044 SolverDestroyPtr_(Solver_);
1045 IsSolverCreated_ =
false;
1047 SolverCreatePtr_ = &Ifpack_Hypre::Hypre_ParCSRPCGCreate;
1048 SolverDestroyPtr_ = &HYPRE_ParCSRPCGDestroy;
1049 SolverSetupPtr_ = &HYPRE_ParCSRPCGSetup;
1050 SolverSolvePtr_ = &HYPRE_ParCSRPCGSolve;
1051 SolverPrecondPtr_ = &HYPRE_ParCSRPCGSetPrecond;
1054 if(IsSolverCreated_){
1055 SolverDestroyPtr_(Solver_);
1056 IsSolverCreated_ =
false;
1058 SolverCreatePtr_ = &Ifpack_Hypre::Hypre_ParCSRGMRESCreate;
1059 SolverDestroyPtr_ = &HYPRE_ParCSRGMRESDestroy;
1060 SolverSetupPtr_ = &HYPRE_ParCSRGMRESSetup;
1061 SolverPrecondPtr_ = &HYPRE_ParCSRGMRESSetPrecond;
1064 if(IsSolverCreated_){
1065 SolverDestroyPtr_(Solver_);
1066 IsSolverCreated_ =
false;
1068 SolverCreatePtr_ = &Ifpack_Hypre::Hypre_ParCSRFlexGMRESCreate;
1069 SolverDestroyPtr_ = &HYPRE_ParCSRFlexGMRESDestroy;
1070 SolverSetupPtr_ = &HYPRE_ParCSRFlexGMRESSetup;
1071 SolverSolvePtr_ = &HYPRE_ParCSRFlexGMRESSolve;
1072 SolverPrecondPtr_ = &HYPRE_ParCSRFlexGMRESSetPrecond;
1075 if(IsSolverCreated_){
1076 SolverDestroyPtr_(Solver_);
1077 IsSolverCreated_ =
false;
1079 SolverCreatePtr_ = &Ifpack_Hypre::Hypre_ParCSRLGMRESCreate;
1080 SolverDestroyPtr_ = &HYPRE_ParCSRLGMRESDestroy;
1081 SolverSetupPtr_ = &HYPRE_ParCSRLGMRESSetup;
1082 SolverSolvePtr_ = &HYPRE_ParCSRLGMRESSolve;
1083 SolverPrecondPtr_ = &HYPRE_ParCSRLGMRESSetPrecond;
1086 if(IsSolverCreated_){
1087 SolverDestroyPtr_(Solver_);
1088 IsSolverCreated_ =
false;
1090 SolverCreatePtr_ = &Ifpack_Hypre::Hypre_ParCSRBiCGSTABCreate;
1091 SolverDestroyPtr_ = &HYPRE_ParCSRBiCGSTABDestroy;
1092 SolverSetupPtr_ = &HYPRE_ParCSRBiCGSTABSetup;
1093 SolverSolvePtr_ = &HYPRE_ParCSRBiCGSTABSolve;
1094 SolverPrecondPtr_ = &HYPRE_ParCSRBiCGSTABSetPrecond;
1107 if(IsPrecondCreated_){
1108 PrecondDestroyPtr_(Preconditioner_);
1109 IsPrecondCreated_ =
false;
1111 PrecondCreatePtr_ = &Ifpack_Hypre::Hypre_BoomerAMGCreate;
1112 PrecondDestroyPtr_ = &HYPRE_BoomerAMGDestroy;
1113 PrecondSetupPtr_ = &HYPRE_BoomerAMGSetup;
1114 PrecondSolvePtr_ = &HYPRE_BoomerAMGSolve;
1117 if(IsPrecondCreated_){
1118 PrecondDestroyPtr_(Preconditioner_);
1119 IsPrecondCreated_ =
false;
1121 PrecondCreatePtr_ = &Ifpack_Hypre::Hypre_ParaSailsCreate;
1122 PrecondDestroyPtr_ = &HYPRE_ParaSailsDestroy;
1123 PrecondSetupPtr_ = &HYPRE_ParaSailsSetup;
1124 PrecondSolvePtr_ = &HYPRE_ParaSailsSolve;
1127 if(IsPrecondCreated_){
1128 PrecondDestroyPtr_(Preconditioner_);
1129 IsPrecondCreated_ =
false;
1131 PrecondCreatePtr_ = &Ifpack_Hypre::Hypre_EuclidCreate;
1132 PrecondDestroyPtr_ = &HYPRE_EuclidDestroy;
1133 PrecondSetupPtr_ = &HYPRE_EuclidSetup;
1134 PrecondSolvePtr_ = &HYPRE_EuclidSolve;
1137 if(IsPrecondCreated_){
1138 PrecondDestroyPtr_(Preconditioner_);
1139 IsPrecondCreated_ =
false;
1141 PrecondCreatePtr_ = &Ifpack_Hypre::Hypre_AMSCreate;
1142 PrecondDestroyPtr_ = &HYPRE_AMSDestroy;
1143 PrecondSetupPtr_ = &HYPRE_AMSSetup;
1144 PrecondSolvePtr_ = &HYPRE_AMSSolve;
1155int Ifpack_Hypre::CreateSolver(){
1157 HYPRE_ParCSRMatrixGetComm(ParMatrix_, &comm);
1158 int ierr = (this->*SolverCreatePtr_)(comm, &Solver_);
1159 IsSolverCreated_ =
true;
1164int Ifpack_Hypre::CreatePrecond(){
1166 HYPRE_ParCSRMatrixGetComm(ParMatrix_, &comm);
1167 int ierr = (this->*PrecondCreatePtr_)(comm, &Preconditioner_);
1168 IsPrecondCreated_ =
true;
1174int Ifpack_Hypre::CopyEpetraToHypre(){
1175 Teuchos::RCP<const Epetra_CrsMatrix> Matrix = Teuchos::rcp_dynamic_cast<const Epetra_CrsMatrix>(A_);
1176 if(Matrix.is_null())
1177 throw std::runtime_error(
"Ifpack_Hypre: Unsupported matrix configuration: Epetra_CrsMatrix required");
1179 std::vector<int> new_indices(Matrix->MaxNumEntries());
1180 for(
int i = 0; i < Matrix->NumMyRows(); i++){
1184 IFPACK_CHK_ERR(Matrix->ExtractMyRowView(i, numEntries, values, indices));
1185 for(
int j = 0; j < numEntries; j++){
1186 new_indices[j] = GloballyContiguousColMap_->GID(indices[j]);
1189 GlobalRow[0] = GloballyContiguousRowMap_->GID(i);
1190 IFPACK_CHK_ERR(HYPRE_IJMatrixSetValues(HypreA_, 1, &numEntries, GlobalRow, new_indices.data(), values));
1193 IFPACK_CHK_ERR(HYPRE_IJMatrixGetObject(HypreA_, (
void**)&ParMatrix_));
1195 HYPRE_ParCSRMatrixPrint(ParMatrix_,
"A.mat");
1200int Ifpack_Hypre::Hypre_BoomerAMGCreate(MPI_Comm , HYPRE_Solver *solver)
1201 {
return HYPRE_BoomerAMGCreate(solver);}
1204int Ifpack_Hypre::Hypre_ParaSailsCreate(MPI_Comm comm, HYPRE_Solver *solver)
1205 {
return HYPRE_ParaSailsCreate(comm, solver);}
1208int Ifpack_Hypre::Hypre_EuclidCreate(MPI_Comm comm, HYPRE_Solver *solver)
1209 {
return HYPRE_EuclidCreate(comm, solver);}
1212int Ifpack_Hypre::Hypre_AMSCreate(MPI_Comm , HYPRE_Solver *solver)
1213 {
return HYPRE_AMSCreate(solver);}
1216int Ifpack_Hypre::Hypre_ParCSRHybridCreate(MPI_Comm , HYPRE_Solver *solver)
1217 {
return HYPRE_ParCSRHybridCreate(solver);}
1220int Ifpack_Hypre::Hypre_ParCSRPCGCreate(MPI_Comm comm, HYPRE_Solver *solver)
1221 {
return HYPRE_ParCSRPCGCreate(comm, solver);}
1224int Ifpack_Hypre::Hypre_ParCSRGMRESCreate(MPI_Comm comm, HYPRE_Solver *solver)
1225 {
return HYPRE_ParCSRGMRESCreate(comm, solver);}
1228int Ifpack_Hypre::Hypre_ParCSRFlexGMRESCreate(MPI_Comm comm, HYPRE_Solver *solver)
1229 {
return HYPRE_ParCSRFlexGMRESCreate(comm, solver);}
1232int Ifpack_Hypre::Hypre_ParCSRLGMRESCreate(MPI_Comm comm, HYPRE_Solver *solver)
1233 {
return HYPRE_ParCSRLGMRESCreate(comm, solver);}
1236int Ifpack_Hypre::Hypre_ParCSRBiCGSTABCreate(MPI_Comm comm, HYPRE_Solver *solver)
1237 {
return HYPRE_ParCSRBiCGSTABCreate(comm, solver);}
1240Teuchos::RCP<const Epetra_Map> Ifpack_Hypre::MakeContiguousColumnMap(Teuchos::RCP<const Epetra_RowMatrix> &MatrixRow)
const{
1245 Teuchos::RCP<const Epetra_CrsMatrix> Matrix = Teuchos::rcp_dynamic_cast<const Epetra_CrsMatrix>(MatrixRow);
1246 if(Matrix.is_null())
1247 throw std::runtime_error(
"Ifpack_Hypre: Unsupported matrix configuration: Epetra_CrsMatrix required");
1248 const Epetra_Map & DomainMap = Matrix->DomainMap();
1249 const Epetra_Map & ColumnMap = Matrix->ColMap();
1254 return rcpFromRef(ColumnMap);
1262 Epetra_IntVector MyGIDsHYPRE(View,DomainMap,ContiguousDomainMap->MyGlobalElements());
1266 ColGIDsHYPRE.Import(MyGIDsHYPRE, *importer, Insert);
Ifpack_CondestType
Ifpack_CondestType: enum to define the type of condition number estimate.
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)
#define IFPACK_CHK_ERRV(ifpack_err)
int NumGlobalElements() const
int NumMyElements() const
double ** Pointers() const
virtual const Epetra_Map & RowMatrixRowMap() const=0