43 #ifndef IFPACK2_HYPRE_DEF_HPP 44 #define IFPACK2_HYPRE_DEF_HPP 46 #include "Ifpack2_Hypre_decl.hpp" 47 #if defined(HAVE_IFPACK2_HYPRE) && defined(HAVE_IFPACK2_MPI) 50 #include "Tpetra_Import.hpp" 51 #include "Teuchos_ParameterList.hpp" 52 #include "Teuchos_RCP.hpp" 53 #include "Teuchos_DefaultMpiComm.hpp" 54 #include "HYPRE_IJ_mv.h" 55 #include "HYPRE_parcsr_ls.h" 57 #include "_hypre_parcsr_mv.h" 58 #include "_hypre_IJ_mv.h" 59 #include "HYPRE_parcsr_mv.h" 65 using Teuchos::rcpFromRef;
70 template<
class MatrixType>
72 Hypre(
const Teuchos::RCP<const row_matrix_type>& A):
74 IsInitialized_(false),
75 IsComputed_(false), NumInitialize_(0),
86 IsSolverCreated_(false),
87 IsPrecondCreated_(false),
88 SolveOrPrec_(Hypre_Is_Solver),
92 UsePreconditioner_(false),
95 MPI_Comm comm = * (Teuchos::rcp_dynamic_cast<
const Teuchos::MpiComm<int> >(A->getRowMap()->getComm())->getRawMpiComm());
100 if (!A_->getRowMap()->isSameAs(*A_->getRangeMap())) {
101 IFPACK2_CHK_ERRV(-1);
104 if (A_->getRowMap()->isContiguous()) {
105 GloballyContiguousRowMap_ = A_->getRowMap();
106 GloballyContiguousColMap_ = A_->getColMap();
109 if(A_->getDomainMap()->isSameAs(*A_->getRowMap())) {
110 Teuchos::RCP<const crs_matrix_type> Aconst = Teuchos::rcp_dynamic_cast<
const crs_matrix_type>(A_);
111 GloballyContiguousColMap_ = MakeContiguousColumnMap(Aconst);
112 GloballyContiguousRowMap_ = rcp(
new map_type(A_->getRowMap()->getGlobalNumElements(),
113 A_->getRowMap()->getNodeNumElements(), 0, A_->getRowMap()->getComm()));
116 throw std::runtime_error(
"Ifpack_Hypre: Unsupported map configuration: Row/Domain maps do not match");
123 IFPACK2_CHK_ERRV(HYPRE_IJVectorCreate(comm, ilower, iupper, &XHypre_));
124 IFPACK2_CHK_ERRV(HYPRE_IJVectorSetObjectType(XHypre_, HYPRE_PARCSR));
125 IFPACK2_CHK_ERRV(HYPRE_IJVectorInitialize(XHypre_));
126 IFPACK2_CHK_ERRV(HYPRE_IJVectorAssemble(XHypre_));
127 IFPACK2_CHK_ERRV(HYPRE_IJVectorGetObject(XHypre_, (
void**) &ParX_));
128 XVec_ = Teuchos::rcp((hypre_ParVector *) hypre_IJVectorObject(((hypre_IJVector *) XHypre_)),
false);
131 IFPACK2_CHK_ERRV(HYPRE_IJVectorCreate(comm, ilower, iupper, &YHypre_));
132 IFPACK2_CHK_ERRV(HYPRE_IJVectorSetObjectType(YHypre_, HYPRE_PARCSR));
133 IFPACK2_CHK_ERRV(HYPRE_IJVectorInitialize(YHypre_));
134 IFPACK2_CHK_ERRV(HYPRE_IJVectorAssemble(YHypre_));
135 IFPACK2_CHK_ERRV(HYPRE_IJVectorGetObject(YHypre_, (
void**) &ParY_));
136 YVec_ = Teuchos::rcp((hypre_ParVector *) hypre_IJVectorObject(((hypre_IJVector *) YHypre_)),
false);
139 VectorCache_.resize(A->getRowMap()->getNodeNumElements());
143 template<
class MatrixType>
144 Hypre<MatrixType>::~Hypre() {
149 template<
class MatrixType>
150 void Hypre<MatrixType>::Destroy(){
152 IFPACK2_CHK_ERRV(HYPRE_IJMatrixDestroy(HypreA_));
154 IFPACK2_CHK_ERRV(HYPRE_IJVectorDestroy(XHypre_));
155 IFPACK2_CHK_ERRV(HYPRE_IJVectorDestroy(YHypre_));
156 if(IsSolverCreated_){
157 IFPACK2_CHK_ERRV(SolverDestroyPtr_(Solver_));
159 if(IsPrecondCreated_){
160 IFPACK2_CHK_ERRV(PrecondDestroyPtr_(Preconditioner_));
165 IFPACK2_CHK_ERRV(HYPRE_IJMatrixDestroy(HypreG_));
168 IFPACK2_CHK_ERRV(HYPRE_IJVectorDestroy(xHypre_));
171 IFPACK2_CHK_ERRV(HYPRE_IJVectorDestroy(yHypre_));
174 IFPACK2_CHK_ERRV(HYPRE_IJVectorDestroy(zHypre_));
179 template<
class MatrixType>
180 void Hypre<MatrixType>::initialize(){
181 const std::string timerName (
"Ifpack2::Hypre::initialize");
182 Teuchos::RCP<Teuchos::Time> timer = Teuchos::TimeMonitor::lookupCounter (timerName);
183 if (timer.is_null ()) timer = Teuchos::TimeMonitor::getNewCounter (timerName);
185 if(IsInitialized_)
return;
186 double startTime = timer->wallTime();
188 Teuchos::TimeMonitor timeMon (*timer);
194 InitializeTime_ += (timer->wallTime() - startTime);
198 template<
class MatrixType>
199 void Hypre<MatrixType>::setParameters(
const Teuchos::ParameterList& list){
201 std::map<std::string, Hypre_Solver> solverMap;
202 solverMap[
"BoomerAMG"] = BoomerAMG;
203 solverMap[
"ParaSails"] = ParaSails;
204 solverMap[
"Euclid"] = Euclid;
205 solverMap[
"AMS"] = AMS;
206 solverMap[
"Hybrid"] = Hybrid;
207 solverMap[
"PCG"] = PCG;
208 solverMap[
"GMRES"] =
GMRES;
209 solverMap[
"FlexGMRES"] = FlexGMRES;
210 solverMap[
"LGMRES"] = LGMRES;
211 solverMap[
"BiCGSTAB"] = BiCGSTAB;
213 std::map<std::string, Hypre_Chooser> chooserMap;
214 chooserMap[
"Solver"] = Hypre_Is_Solver;
215 chooserMap[
"Preconditioner"] = Hypre_Is_Preconditioner;
218 Hypre_Solver solType;
219 if (list.isType<std::string>(
"hypre: Solver"))
220 solType = solverMap[list.get<std::string>(
"hypre: Solver")];
221 else if(list.isParameter(
"hypre: Solver"))
222 solType = (Hypre_Solver) list.get<
int>(
"hypre: Solver");
225 SolverType_ = solType;
226 Hypre_Solver precType;
227 if (list.isType<std::string>(
"hypre: Preconditioner"))
228 precType = solverMap[list.get<std::string>(
"hypre: Preconditioner")];
229 else if(list.isParameter(
"hypre: Preconditioner"))
230 precType = (Hypre_Solver) list.get<
int>(
"hypre: Preconditioner");
233 PrecondType_ = precType;
234 Hypre_Chooser chooser;
235 if (list.isType<std::string>(
"hypre: SolveOrPrecondition"))
236 chooser = chooserMap[list.get<std::string>(
"hypre: SolveOrPrecondition")];
237 else if(list.isParameter(
"hypre: SolveOrPrecondition"))
238 chooser = (Hypre_Chooser) list.get<
int>(
"hypre: SolveOrPrecondition");
240 chooser = Hypre_Is_Solver;
241 SolveOrPrec_ = chooser;
242 bool SetPrecond = list.isParameter(
"hypre: SetPreconditioner") ? list.get<
bool>(
"hypre: SetPreconditioner") :
false;
243 IFPACK2_CHK_ERR(SetParameter(SetPrecond));
244 int NumFunctions = list.isParameter(
"hypre: NumFunctions") ? list.get<
int>(
"hypre: NumFunctions") : 0;
247 if(NumFunctions > 0){
248 RCP<FunctionParameter>* params = list.get<RCP<FunctionParameter>*>(
"hypre: Functions");
249 for(
int i = 0; i < NumFunctions; i++){
250 IFPACK2_CHK_ERR(AddFunToList(params[i]));
254 if (list.isSublist(
"hypre: Solver functions")) {
255 Teuchos::ParameterList solverList = list.sublist(
"hypre: Solver functions");
256 for (
auto it = solverList.begin(); it != solverList.end(); ++it) {
257 std::string funct_name = it->first;
258 if (it->second.isType<HYPRE_Int>()) {
259 IFPACK2_CHK_ERR(AddFunToList(rcp(
new FunctionParameter(Hypre_Is_Solver, funct_name , Teuchos::getValue<HYPRE_Int>(it->second)))));
260 }
else if (!std::is_same<HYPRE_Int,int>::value && it->second.isType<
int>()) {
261 IFPACK2_CHK_ERR(AddFunToList(rcp(
new FunctionParameter(Hypre_Is_Solver, funct_name , Teuchos::as<global_ordinal_type>(Teuchos::getValue<int>(it->second))))));
262 }
else if (it->second.isType<
double>()) {
263 IFPACK2_CHK_ERR(AddFunToList(rcp(
new FunctionParameter(Hypre_Is_Solver, funct_name , Teuchos::getValue<double>(it->second)))));
270 if (list.isSublist(
"hypre: Preconditioner functions")) {
271 Teuchos::ParameterList precList = list.sublist(
"hypre: Preconditioner functions");
272 for (
auto it = precList.begin(); it != precList.end(); ++it) {
273 std::string funct_name = it->first;
274 if (it->second.isType<HYPRE_Int>()) {
275 IFPACK2_CHK_ERR(AddFunToList(rcp(
new FunctionParameter(Hypre_Is_Preconditioner, funct_name , Teuchos::getValue<HYPRE_Int>(it->second)))));
276 }
else if (!std::is_same<HYPRE_Int,int>::value && it->second.isType<
int>()) {
277 IFPACK2_CHK_ERR(AddFunToList(rcp(
new FunctionParameter(Hypre_Is_Preconditioner, funct_name , Teuchos::as<global_ordinal_type>(Teuchos::getValue<int>(it->second))))));
278 }
else if (it->second.isType<
double>()) {
279 IFPACK2_CHK_ERR(AddFunToList(rcp(
new FunctionParameter(Hypre_Is_Preconditioner, funct_name , Teuchos::getValue<double>(it->second)))));
280 }
else if (it->second.isList()) {
281 Teuchos::ParameterList pl = Teuchos::getValue<Teuchos::ParameterList>(it->second);
282 if (FunctionParameter::isFuncIntInt(funct_name)) {
283 HYPRE_Int arg0 = pl.get<
int>(
"arg 0");
284 HYPRE_Int arg1 = pl.get<
int>(
"arg 1");
285 IFPACK2_CHK_ERR(AddFunToList(rcp(
new FunctionParameter(Hypre_Is_Preconditioner, funct_name , arg0, arg1))));
286 }
else if (FunctionParameter::isFuncIntIntDoubleDouble(funct_name)) {
287 HYPRE_Int arg0 = pl.get<
int>(
"arg 0");
288 HYPRE_Int arg1 = pl.get<
int>(
"arg 1");
289 double arg2 = pl.get<
double>(
"arg 2");
290 double arg3 = pl.get<
double>(
"arg 3");
291 IFPACK2_CHK_ERR(AddFunToList(rcp(
new FunctionParameter(Hypre_Is_Preconditioner, funct_name , arg0, arg1, arg2, arg3))));
292 }
else if (FunctionParameter::isFuncIntIntIntDoubleIntInt(funct_name)) {
293 HYPRE_Int arg0 = pl.get<
int>(
"arg 0");
294 HYPRE_Int arg1 = pl.get<
int>(
"arg 1");
295 HYPRE_Int arg2 = pl.get<
int>(
"arg 2");
296 double arg3 = pl.get<
double>(
"arg 3");
297 HYPRE_Int arg4 = pl.get<
int>(
"arg 4");
298 HYPRE_Int arg5 = pl.get<
int>(
"arg 5");
299 IFPACK2_CHK_ERR(AddFunToList(rcp(
new FunctionParameter(Hypre_Is_Preconditioner, funct_name , arg0, arg1, arg2, arg3, arg4, arg5))));
307 if (list.isSublist(
"Coordinates") && list.sublist(
"Coordinates").isType<Teuchos::RCP<multivector_type> >(
"Coordinates"))
308 Coords_ = list.sublist(
"Coordinates").get<Teuchos::RCP<multivector_type> >(
"Coordinates");
309 if (list.isSublist(
"Operators") && list.sublist(
"Operators").isType<Teuchos::RCP<const crs_matrix_type> >(
"G"))
310 G_ = list.sublist(
"Operators").get<Teuchos::RCP<const crs_matrix_type> >(
"G");
312 Dump_ = list.isParameter(
"hypre: Dump") ? list.get<
bool>(
"hypre: Dump") :
false;
316 template<
class MatrixType>
317 int Hypre<MatrixType>::AddFunToList(RCP<FunctionParameter> NewFun){
318 NumFunsToCall_ = NumFunsToCall_+1;
319 FunsToCall_.resize(NumFunsToCall_);
320 FunsToCall_[NumFunsToCall_-1] = NewFun;
325 template<
class MatrixType>
326 int Hypre<MatrixType>::SetParameter(Hypre_Chooser chooser, global_ordinal_type (*pt2Func)(HYPRE_Solver, global_ordinal_type), global_ordinal_type parameter){
327 RCP<FunctionParameter> temp = rcp(
new FunctionParameter(chooser, pt2Func, parameter));
328 IFPACK2_CHK_ERR(AddFunToList(temp));
333 template<
class MatrixType>
334 int Hypre<MatrixType>::SetParameter(Hypre_Chooser chooser, global_ordinal_type (*pt2Func)(HYPRE_Solver,
double),
double parameter){
335 RCP<FunctionParameter> temp = rcp(
new FunctionParameter(chooser, pt2Func, parameter));
336 IFPACK2_CHK_ERR(AddFunToList(temp));
341 template<
class MatrixType>
342 int Hypre<MatrixType>::SetParameter(Hypre_Chooser chooser, global_ordinal_type (*pt2Func)(HYPRE_Solver,
double, global_ordinal_type),
double parameter1, global_ordinal_type parameter2){
343 RCP<FunctionParameter> temp = rcp(
new FunctionParameter(chooser, pt2Func, parameter1, parameter2));
344 IFPACK2_CHK_ERR(AddFunToList(temp));
349 template<
class MatrixType>
350 int Hypre<MatrixType>::SetParameter(Hypre_Chooser chooser, global_ordinal_type (*pt2Func)(HYPRE_Solver, global_ordinal_type, global_ordinal_type), global_ordinal_type parameter1, global_ordinal_type parameter2){
351 RCP<FunctionParameter> temp = rcp(
new FunctionParameter(chooser, pt2Func, parameter1, parameter2));
352 IFPACK2_CHK_ERR(AddFunToList(temp));
357 template<
class MatrixType>
358 int Hypre<MatrixType>::SetParameter(Hypre_Chooser chooser, global_ordinal_type (*pt2Func)(HYPRE_Solver,
double*),
double* parameter){
359 RCP<FunctionParameter> temp = rcp(
new FunctionParameter(chooser, pt2Func, parameter));
360 IFPACK2_CHK_ERR(AddFunToList(temp));
365 template<
class MatrixType>
366 int Hypre<MatrixType>::SetParameter(Hypre_Chooser chooser, global_ordinal_type (*pt2Func)(HYPRE_Solver, global_ordinal_type*), global_ordinal_type* parameter){
367 RCP<FunctionParameter> temp = rcp(
new FunctionParameter(chooser, pt2Func, parameter));
368 IFPACK2_CHK_ERR(AddFunToList(temp));
373 template<
class MatrixType>
374 int Hypre<MatrixType>::SetParameter(Hypre_Chooser chooser, global_ordinal_type (*pt2Func)(HYPRE_Solver, global_ordinal_type**), global_ordinal_type** parameter){
375 RCP<FunctionParameter> temp = rcp(
new FunctionParameter(chooser, pt2Func, parameter));
376 IFPACK2_CHK_ERR(AddFunToList(temp));
381 template<
class MatrixType>
382 int Hypre<MatrixType>::SetParameter(Hypre_Chooser chooser, Hypre_Solver solver){
383 if(chooser == Hypre_Is_Solver){
384 SolverType_ = solver;
386 PrecondType_ = solver;
392 template<
class MatrixType>
393 int Hypre<MatrixType>::SetDiscreteGradient(Teuchos::RCP<const crs_matrix_type> G){
394 using LO = local_ordinal_type;
395 using GO = global_ordinal_type;
396 using SC = scalar_type;
399 if(!A_->getRowMap()->isSameAs(*G->getRowMap()))
400 throw std::runtime_error(
"Hypre<MatrixType>: Edge map mismatch: A and discrete gradient");
403 GloballyContiguousNodeRowMap_ = rcp(
new map_type(G->getDomainMap()->getGlobalNumElements(),
404 G->getDomainMap()->getNodeNumElements(), 0, A_->getRowMap()->getComm()));
405 GloballyContiguousNodeColMap_ = MakeContiguousColumnMap(G);
408 MPI_Comm comm = * (Teuchos::rcp_dynamic_cast<
const Teuchos::MpiComm<int> >(A_->getRowMap()->getComm())->getRawMpiComm());
409 GO ilower = GloballyContiguousRowMap_->getMinGlobalIndex();
410 GO iupper = GloballyContiguousRowMap_->getMaxGlobalIndex();
411 GO jlower = GloballyContiguousNodeRowMap_->getMinGlobalIndex();
412 GO jupper = GloballyContiguousNodeRowMap_->getMaxGlobalIndex();
413 IFPACK2_CHK_ERR(HYPRE_IJMatrixCreate(comm, ilower, iupper, jlower, jupper, &HypreG_));
414 IFPACK2_CHK_ERR(HYPRE_IJMatrixSetObjectType(HypreG_, HYPRE_PARCSR));
415 IFPACK2_CHK_ERR(HYPRE_IJMatrixInitialize(HypreG_));
417 std::vector<GO> new_indices(G->getNodeMaxNumRowEntries());
418 for(LO i = 0; i < (LO)G->getNodeNumRows(); i++){
419 Teuchos::ArrayView<const SC> values;
420 Teuchos::ArrayView<const LO> indices;
421 G->getLocalRowView(i, indices, values);
422 for(LO j = 0; j < (LO) indices.size(); j++){
423 new_indices[j] = GloballyContiguousNodeColMap_->getGlobalElement(indices[j]);
426 GO numEntries = (GO) indices.size();
427 GlobalRow[0] = GloballyContiguousRowMap_->getGlobalElement(i);
428 IFPACK2_CHK_ERR(HYPRE_IJMatrixSetValues(HypreG_, 1, &numEntries, GlobalRow, new_indices.data(), values.getRawPtr()));
430 IFPACK2_CHK_ERR(HYPRE_IJMatrixAssemble(HypreG_));
431 IFPACK2_CHK_ERR(HYPRE_IJMatrixGetObject(HypreG_, (
void**)&ParMatrixG_));
434 HYPRE_ParCSRMatrixPrint(ParMatrixG_,
"G.mat");
436 if(SolverType_ == AMS)
437 HYPRE_AMSSetDiscreteGradient(Solver_, ParMatrixG_);
438 if(PrecondType_ == AMS)
439 HYPRE_AMSSetDiscreteGradient(Preconditioner_, ParMatrixG_);
444 template<
class MatrixType>
445 int Hypre<MatrixType>::SetCoordinates(Teuchos::RCP<multivector_type> coords) {
447 if(!G_.is_null() && !G_->getDomainMap()->isSameAs(*coords->getMap()))
448 throw std::runtime_error(
"Hypre<MatrixType>: Node map mismatch: G->DomainMap() and coords");
450 if(SolverType_ != AMS && PrecondType_ != AMS)
453 scalar_type *xPtr = coords->getDataNonConst(0).getRawPtr();
454 scalar_type *yPtr = coords->getDataNonConst(1).getRawPtr();
455 scalar_type *zPtr = coords->getDataNonConst(2).getRawPtr();
457 MPI_Comm comm = * (Teuchos::rcp_dynamic_cast<
const Teuchos::MpiComm<int> >(A_->getRowMap()->getComm())->getRawMpiComm());
458 local_ordinal_type NumEntries = coords->getLocalLength();
459 global_ordinal_type * indices =
const_cast<global_ordinal_type*
>(GloballyContiguousNodeRowMap_->getNodeElementList().getRawPtr());
461 global_ordinal_type ilower = GloballyContiguousNodeRowMap_->getMinGlobalIndex();
462 global_ordinal_type iupper = GloballyContiguousNodeRowMap_->getMaxGlobalIndex();
464 if( NumEntries != iupper-ilower+1) {
465 std::cout<<
"Ifpack2::Hypre::SetCoordinates(): Error on rank "<<A_->getRowMap()->getComm()->getRank()<<
": MyLength = "<<coords->getLocalLength()<<
" GID range = ["<<ilower<<
","<<iupper<<
"]"<<std::endl;
466 throw std::runtime_error(
"Hypre<MatrixType>: SetCoordinates: Length mismatch");
469 IFPACK2_CHK_ERR(HYPRE_IJVectorCreate(comm, ilower, iupper, &xHypre_));
470 IFPACK2_CHK_ERR(HYPRE_IJVectorSetObjectType(xHypre_, HYPRE_PARCSR));
471 IFPACK2_CHK_ERR(HYPRE_IJVectorInitialize(xHypre_));
472 IFPACK2_CHK_ERR(HYPRE_IJVectorSetValues(xHypre_,NumEntries,indices,xPtr));
473 IFPACK2_CHK_ERR(HYPRE_IJVectorAssemble(xHypre_));
474 IFPACK2_CHK_ERR(HYPRE_IJVectorGetObject(xHypre_, (
void**) &xPar_));
476 IFPACK2_CHK_ERR(HYPRE_IJVectorCreate(comm, ilower, iupper, &yHypre_));
477 IFPACK2_CHK_ERR(HYPRE_IJVectorSetObjectType(yHypre_, HYPRE_PARCSR));
478 IFPACK2_CHK_ERR(HYPRE_IJVectorInitialize(yHypre_));
479 IFPACK2_CHK_ERR(HYPRE_IJVectorSetValues(yHypre_,NumEntries,indices,yPtr));
480 IFPACK2_CHK_ERR(HYPRE_IJVectorAssemble(yHypre_));
481 IFPACK2_CHK_ERR(HYPRE_IJVectorGetObject(yHypre_, (
void**) &yPar_));
483 IFPACK2_CHK_ERR(HYPRE_IJVectorCreate(comm, ilower, iupper, &zHypre_));
484 IFPACK2_CHK_ERR(HYPRE_IJVectorSetObjectType(zHypre_, HYPRE_PARCSR));
485 IFPACK2_CHK_ERR(HYPRE_IJVectorInitialize(zHypre_));
486 IFPACK2_CHK_ERR(HYPRE_IJVectorSetValues(zHypre_,NumEntries,indices,zPtr));
487 IFPACK2_CHK_ERR(HYPRE_IJVectorAssemble(zHypre_));
488 IFPACK2_CHK_ERR(HYPRE_IJVectorGetObject(zHypre_, (
void**) &zPar_));
491 HYPRE_ParVectorPrint(xPar_,
"coordX.dat");
492 HYPRE_ParVectorPrint(yPar_,
"coordY.dat");
493 HYPRE_ParVectorPrint(zPar_,
"coordZ.dat");
496 if(SolverType_ == AMS)
497 HYPRE_AMSSetCoordinateVectors(Solver_, xPar_, yPar_, zPar_);
498 if(PrecondType_ == AMS)
499 HYPRE_AMSSetCoordinateVectors(Preconditioner_, xPar_, yPar_, zPar_);
506 template<
class MatrixType>
507 void Hypre<MatrixType>::compute(){
508 const std::string timerName (
"Ifpack2::Hypre::compute");
509 Teuchos::RCP<Teuchos::Time> timer = Teuchos::TimeMonitor::lookupCounter (timerName);
510 if (timer.is_null ()) timer = Teuchos::TimeMonitor::getNewCounter (timerName);
511 double startTime = timer->wallTime();
514 Teuchos::TimeMonitor timeMon (*timer);
516 if(isInitialized() ==
false){
525 MPI_Comm comm = * (Teuchos::rcp_dynamic_cast<
const Teuchos::MpiComm<int> >(A_->getRowMap()->getComm())->getRawMpiComm());
526 global_ordinal_type ilower = GloballyContiguousRowMap_->getMinGlobalIndex();
527 global_ordinal_type iupper = GloballyContiguousRowMap_->getMaxGlobalIndex();
528 IFPACK2_CHK_ERR(HYPRE_IJMatrixCreate(comm, ilower, iupper, ilower, iupper, &HypreA_));
529 IFPACK2_CHK_ERR(HYPRE_IJMatrixSetObjectType(HypreA_, HYPRE_PARCSR));
530 IFPACK2_CHK_ERR(HYPRE_IJMatrixInitialize(HypreA_));
532 if(SolveOrPrec_ == Hypre_Is_Solver) {
533 IFPACK2_CHK_ERR(SetSolverType(SolverType_));
534 if (SolverPrecondPtr_ != NULL && UsePreconditioner_) {
536 IFPACK2_CHK_ERR(SetPrecondType(PrecondType_));
538 IFPACK2_CHK_ERR(SolverPrecondPtr_(Solver_, PrecondSolvePtr_, PrecondSetupPtr_, Preconditioner_));
543 IFPACK2_CHK_ERR(SetPrecondType(PrecondType_));
548 SetDiscreteGradient(G_);
551 if (!Coords_.is_null()) {
552 SetCoordinates(Coords_);
556 if(SolveOrPrec_ == Hypre_Is_Solver){
557 IFPACK2_CHK_ERR(SolverSetupPtr_(Solver_, ParMatrix_, ParX_, ParY_));
559 IFPACK2_CHK_ERR(PrecondSetupPtr_(Preconditioner_, ParMatrix_, ParX_, ParY_));
566 ComputeTime_ += (timer->wallTime() - startTime);
570 template<
class MatrixType>
571 int Hypre<MatrixType>::CallFunctions()
const{
572 for(
int i = 0; i < NumFunsToCall_; i++){
573 IFPACK2_CHK_ERR(FunsToCall_[i]->CallFunction(Solver_, Preconditioner_));
579 template<
class MatrixType>
580 void Hypre<MatrixType>::apply (
const Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& X,
581 Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& Y,
582 Teuchos::ETransp mode,
584 scalar_type beta)
const {
585 using LO = local_ordinal_type;
586 using SC = scalar_type;
587 const std::string timerName (
"Ifpack2::Hypre::apply");
588 Teuchos::RCP<Teuchos::Time> timer = Teuchos::TimeMonitor::lookupCounter (timerName);
589 if (timer.is_null ()) timer = Teuchos::TimeMonitor::getNewCounter (timerName);
590 double startTime = timer->wallTime();
593 Teuchos::TimeMonitor timeMon (*timer);
595 if(isComputed() ==
false){
598 hypre_Vector *XLocal_ = hypre_ParVectorLocalVector(XVec_);
599 hypre_Vector *YLocal_ = hypre_ParVectorLocalVector(YVec_);
601 bool SameVectors =
false;
602 size_t NumVectors = X.getNumVectors();
603 if (NumVectors != Y.getNumVectors()) IFPACK2_CHK_ERR(-1);
612 for(
int VecNum = 0; VecNum < (int) NumVectors; VecNum++) {
615 SC * XValues =
const_cast<SC*
>(X.getData(VecNum).getRawPtr());
618 YValues =
const_cast<SC*
>(Y.getData(VecNum).getRawPtr());
620 YValues = VectorCache_.getRawPtr();
623 SC *XTemp = XLocal_->data;
625 XLocal_->data = XValues;
626 SC *YTemp = YLocal_->data;
627 YLocal_->data = YValues;
629 IFPACK2_CHK_ERR(HYPRE_ParVectorSetConstantValues(ParY_, 0.0));
630 if(SolveOrPrec_ == Hypre_Is_Solver){
632 IFPACK2_CHK_ERR(SolverSolvePtr_(Solver_, ParMatrix_, ParX_, ParY_));
635 IFPACK2_CHK_ERR(PrecondSolvePtr_(Preconditioner_, ParMatrix_, ParX_, ParY_));
639 Teuchos::ArrayView<SC> Yv = Y.getDataNonConst(VecNum)();
640 LO NumEntries = Y.getLocalLength();
641 for(LO i = 0; i < NumEntries; i++)
644 XLocal_->data = XTemp;
645 YLocal_->data = YTemp;
649 ApplyTime_ += (timer->wallTime() - startTime);
654 template<
class MatrixType>
655 void Hypre<MatrixType>::applyMat (
const Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& X,
656 Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& Y,
657 Teuchos::ETransp mode)
const {
662 template<
class MatrixType>
663 std::string Hypre<MatrixType>::description()
const {
664 std::ostringstream out;
669 out <<
"\"Ifpack2::Hypre\": {";
670 out <<
"Initialized: " << (isInitialized () ?
"true" :
"false") <<
", " 671 <<
"Computed: " << (isComputed () ?
"true" :
"false") <<
", ";
674 out <<
"Matrix: null";
677 out <<
"Global matrix dimensions: [" 678 << A_->getGlobalNumRows () <<
", " 679 << A_->getGlobalNumCols () <<
"]" 680 <<
", Global nnz: " << A_->getGlobalNumEntries();
688 template<
class MatrixType>
689 void Hypre<MatrixType>::describe(Teuchos::FancyOStream &os,
const Teuchos::EVerbosityLevel verbLevel)
const {
692 os <<
"================================================================================" << endl;
693 os <<
"Ifpack2::Hypre: " << endl << endl;
694 os <<
"Using " << A_->getComm()->getSize() <<
" processors." << endl;
695 os <<
"Global number of rows = " << A_->getGlobalNumRows() << endl;
696 os <<
"Global number of nonzeros = " << A_->getGlobalNumEntries() << endl;
699 os <<
"Phase # calls Total Time (s)"<<endl;
700 os <<
"----- ------- --------------"<<endl;
701 os <<
"Initialize() " << std::setw(5) << NumInitialize_
702 <<
" " << std::setw(15) << InitializeTime_<<endl;
703 os <<
"Compute() " << std::setw(5) << NumCompute_
704 <<
" " << std::setw(15) << ComputeTime_ << endl;
705 os <<
"ApplyInverse() " << std::setw(5) << NumApply_
706 <<
" " << std::setw(15) << ApplyTime_ <<endl;
707 os <<
"================================================================================" << endl;
712 template<
class MatrixType>
713 int Hypre<MatrixType>::SetSolverType(Hypre_Solver Solver){
716 if(IsSolverCreated_){
717 SolverDestroyPtr_(Solver_);
718 IsSolverCreated_ =
false;
720 SolverCreatePtr_ = &Hypre<MatrixType>::Hypre_BoomerAMGCreate;
721 SolverDestroyPtr_ = &HYPRE_BoomerAMGDestroy;
722 SolverSetupPtr_ = &HYPRE_BoomerAMGSetup;
723 SolverPrecondPtr_ = NULL;
724 SolverSolvePtr_ = &HYPRE_BoomerAMGSolve;
727 if(IsSolverCreated_){
728 SolverDestroyPtr_(Solver_);
729 IsSolverCreated_ =
false;
731 SolverCreatePtr_ = &Hypre<MatrixType>::Hypre_AMSCreate;
732 SolverDestroyPtr_ = &HYPRE_AMSDestroy;
733 SolverSetupPtr_ = &HYPRE_AMSSetup;
734 SolverSolvePtr_ = &HYPRE_AMSSolve;
735 SolverPrecondPtr_ = NULL;
738 if(IsSolverCreated_){
739 SolverDestroyPtr_(Solver_);
740 IsSolverCreated_ =
false;
742 SolverCreatePtr_ = &Hypre<MatrixType>::Hypre_ParCSRHybridCreate;
743 SolverDestroyPtr_ = &HYPRE_ParCSRHybridDestroy;
744 SolverSetupPtr_ = &HYPRE_ParCSRHybridSetup;
745 SolverSolvePtr_ = &HYPRE_ParCSRHybridSolve;
746 SolverPrecondPtr_ = &HYPRE_ParCSRHybridSetPrecond;
749 if(IsSolverCreated_){
750 SolverDestroyPtr_(Solver_);
751 IsSolverCreated_ =
false;
753 SolverCreatePtr_ = &Hypre<MatrixType>::Hypre_ParCSRPCGCreate;
754 SolverDestroyPtr_ = &HYPRE_ParCSRPCGDestroy;
755 SolverSetupPtr_ = &HYPRE_ParCSRPCGSetup;
756 SolverSolvePtr_ = &HYPRE_ParCSRPCGSolve;
757 SolverPrecondPtr_ = &HYPRE_ParCSRPCGSetPrecond;
760 if(IsSolverCreated_){
761 SolverDestroyPtr_(Solver_);
762 IsSolverCreated_ =
false;
764 SolverCreatePtr_ = &Hypre<MatrixType>::Hypre_ParCSRGMRESCreate;
765 SolverDestroyPtr_ = &HYPRE_ParCSRGMRESDestroy;
766 SolverSetupPtr_ = &HYPRE_ParCSRGMRESSetup;
767 SolverPrecondPtr_ = &HYPRE_ParCSRGMRESSetPrecond;
770 if(IsSolverCreated_){
771 SolverDestroyPtr_(Solver_);
772 IsSolverCreated_ =
false;
774 SolverCreatePtr_ = &Hypre<MatrixType>::Hypre_ParCSRFlexGMRESCreate;
775 SolverDestroyPtr_ = &HYPRE_ParCSRFlexGMRESDestroy;
776 SolverSetupPtr_ = &HYPRE_ParCSRFlexGMRESSetup;
777 SolverSolvePtr_ = &HYPRE_ParCSRFlexGMRESSolve;
778 SolverPrecondPtr_ = &HYPRE_ParCSRFlexGMRESSetPrecond;
781 if(IsSolverCreated_){
782 SolverDestroyPtr_(Solver_);
783 IsSolverCreated_ =
false;
785 SolverCreatePtr_ = &Hypre<MatrixType>::Hypre_ParCSRLGMRESCreate;
786 SolverDestroyPtr_ = &HYPRE_ParCSRLGMRESDestroy;
787 SolverSetupPtr_ = &HYPRE_ParCSRLGMRESSetup;
788 SolverSolvePtr_ = &HYPRE_ParCSRLGMRESSolve;
789 SolverPrecondPtr_ = &HYPRE_ParCSRLGMRESSetPrecond;
792 if(IsSolverCreated_){
793 SolverDestroyPtr_(Solver_);
794 IsSolverCreated_ =
false;
796 SolverCreatePtr_ = &Hypre<MatrixType>::Hypre_ParCSRBiCGSTABCreate;
797 SolverDestroyPtr_ = &HYPRE_ParCSRBiCGSTABDestroy;
798 SolverSetupPtr_ = &HYPRE_ParCSRBiCGSTABSetup;
799 SolverSolvePtr_ = &HYPRE_ParCSRBiCGSTABSolve;
800 SolverPrecondPtr_ = &HYPRE_ParCSRBiCGSTABSetPrecond;
810 template<
class MatrixType>
811 int Hypre<MatrixType>::SetPrecondType(Hypre_Solver Precond){
814 if(IsPrecondCreated_){
815 PrecondDestroyPtr_(Preconditioner_);
816 IsPrecondCreated_ =
false;
818 PrecondCreatePtr_ = &Hypre<MatrixType>::Hypre_BoomerAMGCreate;
819 PrecondDestroyPtr_ = &HYPRE_BoomerAMGDestroy;
820 PrecondSetupPtr_ = &HYPRE_BoomerAMGSetup;
821 PrecondSolvePtr_ = &HYPRE_BoomerAMGSolve;
824 if(IsPrecondCreated_){
825 PrecondDestroyPtr_(Preconditioner_);
826 IsPrecondCreated_ =
false;
828 PrecondCreatePtr_ = &Hypre<MatrixType>::Hypre_ParaSailsCreate;
829 PrecondDestroyPtr_ = &HYPRE_ParaSailsDestroy;
830 PrecondSetupPtr_ = &HYPRE_ParaSailsSetup;
831 PrecondSolvePtr_ = &HYPRE_ParaSailsSolve;
834 if(IsPrecondCreated_){
835 PrecondDestroyPtr_(Preconditioner_);
836 IsPrecondCreated_ =
false;
838 PrecondCreatePtr_ = &Hypre<MatrixType>::Hypre_EuclidCreate;
839 PrecondDestroyPtr_ = &HYPRE_EuclidDestroy;
840 PrecondSetupPtr_ = &HYPRE_EuclidSetup;
841 PrecondSolvePtr_ = &HYPRE_EuclidSolve;
844 if(IsPrecondCreated_){
845 PrecondDestroyPtr_(Preconditioner_);
846 IsPrecondCreated_ =
false;
848 PrecondCreatePtr_ = &Hypre<MatrixType>::Hypre_AMSCreate;
849 PrecondDestroyPtr_ = &HYPRE_AMSDestroy;
850 PrecondSetupPtr_ = &HYPRE_AMSSetup;
851 PrecondSolvePtr_ = &HYPRE_AMSSolve;
862 template<
class MatrixType>
863 int Hypre<MatrixType>::CreateSolver(){
865 HYPRE_ParCSRMatrixGetComm(ParMatrix_, &comm);
866 int ierr = (this->*SolverCreatePtr_)(comm, &Solver_);
867 IsSolverCreated_ =
true;
872 template<
class MatrixType>
873 int Hypre<MatrixType>::CreatePrecond(){
875 HYPRE_ParCSRMatrixGetComm(ParMatrix_, &comm);
876 int ierr = (this->*PrecondCreatePtr_)(comm, &Preconditioner_);
877 IsPrecondCreated_ =
true;
883 template<
class MatrixType>
884 int Hypre<MatrixType>::CopyTpetraToHypre(){
885 using LO = local_ordinal_type;
886 using GO = global_ordinal_type;
887 using SC = scalar_type;
889 Teuchos::RCP<const crs_matrix_type> Matrix = Teuchos::rcp_dynamic_cast<
const crs_matrix_type>(A_);
891 throw std::runtime_error(
"Hypre<MatrixType>: Unsupported matrix configuration: Tpetra::CrsMatrix required");
893 std::vector<GO> new_indices(Matrix->getNodeMaxNumRowEntries());
894 for(LO i = 0; i < (LO) Matrix->getNodeNumRows(); i++){
895 Teuchos::ArrayView<const SC> values;
896 Teuchos::ArrayView<const LO> indices;
897 Matrix->getLocalRowView(i, indices, values);
898 for(LO j = 0; j < (LO)indices.size(); j++){
899 new_indices[j] = GloballyContiguousColMap_->getGlobalElement(indices[j]);
902 GO numEntries = (GO) indices.size();
903 GlobalRow[0] = GloballyContiguousRowMap_->getGlobalElement(i);
904 IFPACK2_CHK_ERR(HYPRE_IJMatrixSetValues(HypreA_, 1, &numEntries, GlobalRow, new_indices.data(), values.getRawPtr()));
906 IFPACK2_CHK_ERR(HYPRE_IJMatrixAssemble(HypreA_));
907 IFPACK2_CHK_ERR(HYPRE_IJMatrixGetObject(HypreA_, (
void**)&ParMatrix_));
909 HYPRE_ParCSRMatrixPrint(ParMatrix_,
"A.mat");
914 template<
class MatrixType>
915 int Hypre<MatrixType>::Hypre_BoomerAMGCreate(MPI_Comm , HYPRE_Solver *solver)
916 {
return HYPRE_BoomerAMGCreate(solver);}
919 template<
class MatrixType>
920 int Hypre<MatrixType>::Hypre_ParaSailsCreate(MPI_Comm comm, HYPRE_Solver *solver)
921 {
return HYPRE_ParaSailsCreate(comm, solver);}
924 template<
class MatrixType>
925 int Hypre<MatrixType>::Hypre_EuclidCreate(MPI_Comm comm, HYPRE_Solver *solver)
926 {
return HYPRE_EuclidCreate(comm, solver);}
929 template<
class MatrixType>
930 int Hypre<MatrixType>::Hypre_AMSCreate(MPI_Comm , HYPRE_Solver *solver)
931 {
return HYPRE_AMSCreate(solver);}
934 template<
class MatrixType>
935 int Hypre<MatrixType>::Hypre_ParCSRHybridCreate(MPI_Comm , HYPRE_Solver *solver)
936 {
return HYPRE_ParCSRHybridCreate(solver);}
939 template<
class MatrixType>
940 int Hypre<MatrixType>::Hypre_ParCSRPCGCreate(MPI_Comm comm, HYPRE_Solver *solver)
941 {
return HYPRE_ParCSRPCGCreate(comm, solver);}
944 template<
class MatrixType>
945 int Hypre<MatrixType>::Hypre_ParCSRGMRESCreate(MPI_Comm comm, HYPRE_Solver *solver)
946 {
return HYPRE_ParCSRGMRESCreate(comm, solver);}
949 template<
class MatrixType>
950 int Hypre<MatrixType>::Hypre_ParCSRFlexGMRESCreate(MPI_Comm comm, HYPRE_Solver *solver)
951 {
return HYPRE_ParCSRFlexGMRESCreate(comm, solver);}
954 template<
class MatrixType>
955 int Hypre<MatrixType>::Hypre_ParCSRLGMRESCreate(MPI_Comm comm, HYPRE_Solver *solver)
956 {
return HYPRE_ParCSRLGMRESCreate(comm, solver);}
959 template<
class MatrixType>
960 int Hypre<MatrixType>::Hypre_ParCSRBiCGSTABCreate(MPI_Comm comm, HYPRE_Solver *solver)
961 {
return HYPRE_ParCSRBiCGSTABCreate(comm, solver);}
964 template<
class MatrixType>
965 Teuchos::RCP<const Tpetra::Map<typename MatrixType::local_ordinal_type, typename MatrixType::global_ordinal_type, typename MatrixType::node_type> >
966 Hypre<MatrixType>::MakeContiguousColumnMap(Teuchos::RCP<const crs_matrix_type> &Matrix)
const{
967 using import_type = Tpetra::Import<local_ordinal_type,global_ordinal_type,node_type>;
968 using go_vector_type = Tpetra::Vector<global_ordinal_type,local_ordinal_type,global_ordinal_type,node_type>;
975 throw std::runtime_error(
"Hypre<MatrixType>: Unsupported matrix configuration: Tpetra::CrsMatrix required");
976 RCP<const map_type> DomainMap = Matrix->getDomainMap();
977 RCP<const map_type> ColumnMap = Matrix->getColMap();
978 RCP<const import_type> importer = Matrix->getGraph()->getImporter();
980 if(DomainMap->isContiguous() ) {
986 Teuchos::RCP<map_type> ContiguousDomainMap = rcp(
new map_type(DomainMap->getGlobalNumElements(),
987 DomainMap->getNodeNumElements(), 0, DomainMap->getComm()));
990 go_vector_type MyGIDsHYPRE(DomainMap,ContiguousDomainMap->getNodeElementList());
993 go_vector_type ColGIDsHYPRE(ColumnMap);
994 ColGIDsHYPRE.doImport(MyGIDsHYPRE, *importer, Tpetra::INSERT);
997 return Teuchos::rcp(
new map_type(ColumnMap->getGlobalNumElements(),ColGIDsHYPRE.getDataNonConst()(),0, ColumnMap->getComm()));
1001 return Teuchos::rcp(
new map_type(ColumnMap->getGlobalNumElements(),ContiguousDomainMap->getNodeElementList(), 0, ColumnMap->getComm()));
1008 template<
class MatrixType>
1009 int Hypre<MatrixType>::getNumInitialize()
const {
1010 return NumInitialize_;
1014 template<
class MatrixType>
1015 int Hypre<MatrixType>::getNumCompute()
const {
1020 template<
class MatrixType>
1021 int Hypre<MatrixType>::getNumApply()
const {
1026 template<
class MatrixType>
1027 double Hypre<MatrixType>::getInitializeTime()
const {
1028 return InitializeTime_;
1032 template<
class MatrixType>
1033 double Hypre<MatrixType>::getComputeTime()
const {
1034 return ComputeTime_;
1038 template<
class MatrixType>
1039 double Hypre<MatrixType>::getApplyTime()
const {
1043 template<
class MatrixType>
1044 Teuchos::RCP<const typename Hypre<MatrixType>::map_type>
1046 getDomainMap ()
const 1048 Teuchos::RCP<const row_matrix_type> A = getMatrix();
1049 TEUCHOS_TEST_FOR_EXCEPTION(
1050 A.is_null (), std::runtime_error,
"Ifpack2::Hypre::getDomainMap: The " 1051 "input matrix A is null. Please call setMatrix() with a nonnull input " 1052 "matrix before calling this method.");
1053 return A->getDomainMap ();
1057 template<
class MatrixType>
1058 Teuchos::RCP<const typename Hypre<MatrixType>::map_type>
1060 getRangeMap ()
const 1062 Teuchos::RCP<const row_matrix_type> A = getMatrix();
1063 TEUCHOS_TEST_FOR_EXCEPTION(
1064 A.is_null (), std::runtime_error,
"Ifpack2::Hypre::getRangeMap: The " 1065 "input matrix A is null. Please call setMatrix() with a nonnull input " 1066 "matrix before calling this method.");
1067 return A->getRangeMap ();
1071 template<
class MatrixType>
1072 void Hypre<MatrixType>::setMatrix (
const Teuchos::RCP<const row_matrix_type>& A)
1074 if (A.getRawPtr () != getMatrix().getRawPtr ()) {
1075 IsInitialized_ =
false;
1076 IsComputed_ =
false;
1082 template<
class MatrixType>
1083 Teuchos::RCP<const typename Hypre<MatrixType>::row_matrix_type>
1090 template<
class MatrixType>
1091 bool Hypre<MatrixType>::hasTransposeApply()
const {
1098 #define IFPACK2_HYPRE_INSTANT(S,LO,GO,N) \ 1099 template class Ifpack2::Hypre< Tpetra::RowMatrix<S, LO, GO, N> >; 1102 #endif // HAVE_HYPRE && HAVE_MPI 1103 #endif // IFPACK2_HYPRE_DEF_HPP MatrixType::global_ordinal_type global_ordinal_type
The type of global indices in the input MatrixType.
Definition: Ifpack2_Hiptmair_decl.hpp:88
Uses AztecOO's GMRES.
Definition: Ifpack2_CondestType.hpp:53
Preconditioners and smoothers for Tpetra sparse matrices.
Definition: Ifpack2_AdditiveSchwarz_decl.hpp:73