IFPACK Development
Loading...
Searching...
No Matches
Ifpack_Hypre.cpp
1/*@HEADER
2// ***********************************************************************
3//
4// Ifpack: Object-Oriented Algebraic Preconditioner Package
5// Copyright (2002) Sandia Corporation
6//
7// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8// license for use of this work by or on behalf of the U.S. Government.
9//
10// Redistribution and use in source and binary forms, with or without
11// modification, are permitted provided that the following conditions are
12// met:
13//
14// 1. Redistributions of source code must retain the above copyright
15// notice, this list of conditions and the following disclaimer.
16//
17// 2. Redistributions in binary form must reproduce the above copyright
18// notice, this list of conditions and the following disclaimer in the
19// documentation and/or other materials provided with the distribution.
20//
21// 3. Neither the name of the Corporation nor the names of the
22// contributors may be used to endorse or promote products derived from
23// this software without specific prior written permission.
24//
25// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36//
37// Questions? Contact Michael A. Heroux (maherou@sandia.gov)
38//
39// ***********************************************************************
40//@HEADER
41*/
42#include "Ifpack_Hypre.h"
43#if defined(HAVE_HYPRE) && defined(HAVE_MPI)
44#include <stdexcept>
45
46#include "Ifpack_Utils.h"
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"
54#include "krylov.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"
59#include "HYPRE.h"
60
61using Teuchos::RCP;
62using Teuchos::rcp;
63using Teuchos::rcpFromRef;
64
65// The Python script that generates the ParameterMap needs to be after these typedefs
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*);
76
78class FunctionParameter{
79 public:
81 FunctionParameter(Hypre_Chooser chooser, int_func funct, int param1) :
82 chooser_(chooser),
83 option_(0),
84 int_func_(funct),
85 int_param1_(param1) {}
86
87 FunctionParameter(Hypre_Chooser chooser, std::string funct_name, int param1) :
88 chooser_(chooser),
89 option_(0),
90 int_func_(hypreMapIntFunc_.at(funct_name)),
91 int_param1_(param1) {}
92
94 FunctionParameter(Hypre_Chooser chooser, double_func funct, double param1):
95 chooser_(chooser),
96 option_(1),
97 double_func_(funct),
98 double_param1_(param1) {}
99
100 FunctionParameter(Hypre_Chooser chooser, std::string funct_name, double param1):
101 chooser_(chooser),
102 option_(1),
103 double_func_(hypreMapDoubleFunc_.at(funct_name)),
104 double_param1_(param1) {}
105
107 FunctionParameter(Hypre_Chooser chooser, double_int_func funct, double param1, int param2):
108 chooser_(chooser),
109 option_(2),
110 double_int_func_(funct),
111 int_param1_(param2),
112 double_param1_(param1) {}
113
114 FunctionParameter(Hypre_Chooser chooser, std::string funct_name, double param1, int param2):
115 chooser_(chooser),
116 option_(2),
117 double_int_func_(hypreMapDoubleIntFunc_.at(funct_name)),
118 int_param1_(param2),
119 double_param1_(param1) {}
120
122 FunctionParameter(Hypre_Chooser chooser, int_int_func funct, int param1, int param2):
123 chooser_(chooser),
124 option_(3),
125 int_int_func_(funct),
126 int_param1_(param1),
127 int_param2_(param2) {}
128
129 FunctionParameter(Hypre_Chooser chooser, std::string funct_name, int param1, int param2):
130 chooser_(chooser),
131 option_(3),
132 int_int_func_(hypreMapIntIntFunc_.at(funct_name)),
133 int_param1_(param1),
134 int_param2_(param2) {}
135
137 FunctionParameter(Hypre_Chooser chooser, int_star_func funct, int *param1):
138 chooser_(chooser),
139 option_(4),
140 int_star_func_(funct),
141 int_star_param_(param1) {}
142
143 FunctionParameter(Hypre_Chooser chooser, std::string funct_name, int *param1):
144 chooser_(chooser),
145 option_(4),
146 int_star_func_(hypreMapIntStarFunc_.at(funct_name)),
147 int_star_param_(param1) {}
148
150 FunctionParameter(Hypre_Chooser chooser, double_star_func funct, double* param1):
151 chooser_(chooser),
152 option_(5),
153 double_star_func_(funct),
154 double_star_param_(param1) {}
155
156 FunctionParameter(Hypre_Chooser chooser, std::string funct_name, double* param1):
157 chooser_(chooser),
158 option_(5),
159 double_star_func_(hypreMapDoubleStarFunc_.at(funct_name)),
160 double_star_param_(param1) {}
161
163 FunctionParameter(Hypre_Chooser chooser, int_int_double_double_func funct, int param1, int param2, double param3, double param4):
164 chooser_(chooser),
165 option_(6),
166 int_int_double_double_func_(funct),
167 int_param1_(param1),
168 int_param2_(param2),
169 double_param1_(param3),
170 double_param2_(param4) {}
171
172 FunctionParameter(Hypre_Chooser chooser, std::string funct_name, int param1, int param2, double param3, double param4):
173 chooser_(chooser),
174 option_(6),
175 int_int_double_double_func_(hypreMapIntIntDoubleDoubleFunc_.at(funct_name)),
176 int_param1_(param1),
177 int_param2_(param2),
178 double_param1_(param3),
179 double_param2_(param4) {}
180
182 FunctionParameter(Hypre_Chooser chooser, int_star_star_func funct, int ** param1):
183 chooser_(chooser),
184 option_(7),
185 int_star_star_func_(funct),
186 int_star_star_param_(param1) {}
187
188 FunctionParameter(Hypre_Chooser chooser, std::string funct_name, int** param1):
189 chooser_(chooser),
190 option_(7),
191 int_star_star_func_(hypreMapIntStarStarFunc_.at(funct_name)),
192 int_star_star_param_(param1) {}
193
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):
196 chooser_(chooser),
197 option_(8),
198 int_int_int_double_int_int_func_(funct),
199 int_param1_(param1),
200 int_param2_(param2),
201 int_param3_(param3),
202 int_param4_(param5),
203 int_param5_(param6),
204 double_param1_(param4) {}
205
206 FunctionParameter(Hypre_Chooser chooser, std::string funct_name, int param1, int param2, int param3, double param4, int param5, int param6):
207 chooser_(chooser),
208 option_(8),
209 int_int_int_double_int_int_func_(hypreMapIntIntIntDoubleIntIntFunc_.at(funct_name)),
210 int_param1_(param1),
211 int_param2_(param2),
212 int_param3_(param3),
213 int_param4_(param5),
214 int_param5_(param6),
215 double_param1_(param4) {}
216
218 FunctionParameter(Hypre_Chooser chooser, char_star_func funct, char *param1):
219 chooser_(chooser),
220 option_(9),
221 char_star_func_(funct),
222 char_star_param_(param1) {}
223
224 FunctionParameter(Hypre_Chooser chooser, std::string funct_name, char *param1):
225 chooser_(chooser),
226 option_(9),
227 char_star_func_(hypreMapCharStarFunc_.at(funct_name)),
228 char_star_param_(param1) {}
229
231 int CallFunction(HYPRE_Solver solver, HYPRE_Solver precond) {
232 if(chooser_ == Solver){
233 if(option_ == 0){
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_);
253 } else {
254 IFPACK_CHK_ERR(-2);
255 }
256 } else {
257 if(option_ == 0){
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_);
277 } else {
278 IFPACK_CHK_ERR(-2);
279 }
280 }
281 }
282
283 static bool isFuncIntInt(std::string funct_name) {
284 return (hypreMapIntIntFunc_.find(funct_name) != hypreMapIntIntFunc_.end());
285 }
286
287 static bool isFuncIntIntDoubleDouble(std::string funct_name) {
288 return (hypreMapIntIntDoubleDoubleFunc_.find(funct_name) != hypreMapIntIntDoubleDoubleFunc_.end());
289 }
290
291 static bool isFuncIntIntIntDoubleIntInt(std::string funct_name) {
292 return (hypreMapIntIntIntDoubleIntIntFunc_.find(funct_name) != hypreMapIntIntIntDoubleIntIntFunc_.end());
293 }
294
295 static bool isFuncIntStarStar(std::string funct_name) {
296 return (hypreMapIntStarStarFunc_.find(funct_name) != hypreMapIntStarStarFunc_.end());
297 }
298
299 private:
300 Hypre_Chooser chooser_;
301 int option_;
302 int_func int_func_;
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_;
312 int int_param1_;
313 int int_param2_;
314 int int_param3_;
315 int int_param4_;
316 int int_param5_;
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_;
323
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_;
334
335};
336
337// NOTE: This really, really needs to be here and not up above, so please don't move it
338#include "Ifpack_HypreParameterMap.h"
339
340Ifpack_Hypre::Ifpack_Hypre(Epetra_RowMatrix* A):
341 A_(rcp(A,false)),
342 UseTranspose_(false),
343 IsInitialized_(false),
344 IsComputed_(false),
345 Label_(),
346 NumInitialize_(0),
347 NumCompute_(0),
348 NumApplyInverse_(0),
349 InitializeTime_(0.0),
350 ComputeTime_(0.0),
351 ApplyInverseTime_(0.0),
352 ComputeFlops_(0.0),
353 ApplyInverseFlops_(0.0),
354 Time_(A_->Comm()),
355 HypreA_(0),
356 HypreG_(0),
357 xHypre_(0),
358 yHypre_(0),
359 zHypre_(0),
360 IsSolverCreated_(false),
361 IsPrecondCreated_(false),
362 SolveOrPrec_(Solver),
363 NumFunsToCall_(0),
364 SolverType_(PCG),
365 PrecondType_(Euclid),
366 UsePreconditioner_(false),
367 Dump_(false)
368{
369 MPI_Comm comm = GetMpiComm();
370 // Check that RowMap and RangeMap are the same. While this could handle the
371 // case where RowMap and RangeMap are permutations, other Ifpack PCs don't
372 // handle this either.
373 if (!A_->RowMatrixRowMap().SameAs(A_->OperatorRangeMap())) {
374 IFPACK_CHK_ERRV(-1);
375 }
376 // Hypre expects the RowMap to be Linear.
377 if (A_->RowMatrixRowMap().LinearMap()) {
378 // note these are non-owning pointers, they are deleted by A_'s destructor
379 GloballyContiguousRowMap_ = rcpFromRef(A_->RowMatrixRowMap());
380 GloballyContiguousColMap_ = rcpFromRef(A_->RowMatrixColMap());
381 } else {
382 // Must create GloballyContiguous Maps for Hypre
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()));
388 }
389 else {
390 throw std::runtime_error("Ifpack_Hypre: Unsupported map configuration: Row/Domain maps do not match");
391 }
392 }
393 // Next create vectors that will be used when ApplyInverse() is called
394 int ilower = GloballyContiguousRowMap_->MinMyGID();
395 int iupper = GloballyContiguousRowMap_->MaxMyGID();
396 // X in AX = Y
397 IFPACK_CHK_ERRV(HYPRE_IJVectorCreate(comm, ilower, iupper, &XHypre_));
398 IFPACK_CHK_ERRV(HYPRE_IJVectorSetObjectType(XHypre_, HYPRE_PARCSR));
399 IFPACK_CHK_ERRV(HYPRE_IJVectorInitialize(XHypre_));
400 IFPACK_CHK_ERRV(HYPRE_IJVectorAssemble(XHypre_));
401 IFPACK_CHK_ERRV(HYPRE_IJVectorGetObject(XHypre_, (void**) &ParX_));
402 XVec_ = Teuchos::rcp((hypre_ParVector *) hypre_IJVectorObject(((hypre_IJVector *) XHypre_)),false);
403
404 // Y in AX = Y
405 IFPACK_CHK_ERRV(HYPRE_IJVectorCreate(comm, ilower, iupper, &YHypre_));
406 IFPACK_CHK_ERRV(HYPRE_IJVectorSetObjectType(YHypre_, HYPRE_PARCSR));
407 IFPACK_CHK_ERRV(HYPRE_IJVectorInitialize(YHypre_));
408 IFPACK_CHK_ERRV(HYPRE_IJVectorAssemble(YHypre_));
409 IFPACK_CHK_ERRV(HYPRE_IJVectorGetObject(YHypre_, (void**) &ParY_));
410 YVec_ = Teuchos::rcp((hypre_ParVector *) hypre_IJVectorObject(((hypre_IJVector *) YHypre_)),false);
411
412 // Cache
413 VectorCache_.resize(A->RowMatrixRowMap().NumMyElements());
414} //Constructor
415
416//==============================================================================
417void Ifpack_Hypre::Destroy(){
418 if(IsInitialized()){
419 IFPACK_CHK_ERRV(HYPRE_IJMatrixDestroy(HypreA_));
420 }
421 IFPACK_CHK_ERRV(HYPRE_IJVectorDestroy(XHypre_));
422 IFPACK_CHK_ERRV(HYPRE_IJVectorDestroy(YHypre_));
423 if(IsSolverCreated_){
424 IFPACK_CHK_ERRV(SolverDestroyPtr_(Solver_));
425 }
426 if(IsPrecondCreated_){
427 IFPACK_CHK_ERRV(PrecondDestroyPtr_(Preconditioner_));
428 }
429
430 // Maxwell
431 if(HypreG_) {
432 IFPACK_CHK_ERRV(HYPRE_IJMatrixDestroy(HypreG_));
433 }
434 if(xHypre_) {
435 IFPACK_CHK_ERRV(HYPRE_IJVectorDestroy(xHypre_));
436 }
437 if(yHypre_) {
438 IFPACK_CHK_ERRV(HYPRE_IJVectorDestroy(yHypre_));
439 }
440 if(zHypre_) {
441 IFPACK_CHK_ERRV(HYPRE_IJVectorDestroy(zHypre_));
442 }
443} //Destroy()
444
445//==============================================================================
446int Ifpack_Hypre::Initialize(){
447 Time_.ResetStartTime();
448 if(IsInitialized_) return 0;
449 // set flags
450 IsInitialized_=true;
451 NumInitialize_ = NumInitialize_ + 1;
452 InitializeTime_ = InitializeTime_ + Time_.ElapsedTime();
453 return 0;
454} //Initialize()
455
456//==============================================================================
457int Ifpack_Hypre::SetParameters(Teuchos::ParameterList& list){
458
459 std::map<std::string, Hypre_Solver> solverMap;
460 solverMap["BoomerAMG"] = BoomerAMG;
461 solverMap["ParaSails"] = ParaSails;
462 solverMap["Euclid"] = Euclid;
463 solverMap["AMS"] = AMS;
464 solverMap["Hybrid"] = Hybrid;
465 solverMap["PCG"] = PCG;
466 solverMap["GMRES"] = GMRES;
467 solverMap["FlexGMRES"] = FlexGMRES;
468 solverMap["LGMRES"] = LGMRES;
469 solverMap["BiCGSTAB"] = BiCGSTAB;
470
471 std::map<std::string, Hypre_Chooser> chooserMap;
472 chooserMap["Solver"] = Solver;
473 chooserMap["Preconditioner"] = Preconditioner;
474
475 List_ = list;
476 Hypre_Solver solType;
477 if (list.isType<std::string>("hypre: Solver"))
478 solType = solverMap[list.get<std::string>("hypre: Solver")];
479 else
480 solType = list.get("hypre: Solver", PCG);
481 SolverType_ = solType;
482 Hypre_Solver precType;
483 if (list.isType<std::string>("hypre: Preconditioner"))
484 precType = solverMap[list.get<std::string>("hypre: Preconditioner")];
485 else
486 precType = list.get("hypre: Preconditioner", Euclid);
487 PrecondType_ = precType;
488 Hypre_Chooser chooser;
489 if (list.isType<std::string>("hypre: SolveOrPrecondition"))
490 chooser = chooserMap[list.get<std::string>("hypre: SolveOrPrecondition")];
491 else
492 chooser = list.get("hypre: SolveOrPrecondition", Solver);
493 SolveOrPrec_ = chooser;
494 bool SetPrecond = list.get("hypre: SetPreconditioner", false);
495 IFPACK_CHK_ERR(SetParameter(SetPrecond));
496 int NumFunctions = list.get("hypre: NumFunctions", 0);
497 FunsToCall_.clear();
498 NumFunsToCall_ = 0;
499 if(NumFunctions > 0){
500 RCP<FunctionParameter>* params = list.get<RCP<FunctionParameter>*>("hypre: Functions");
501 for(int i = 0; i < NumFunctions; i++){
502 IFPACK_CHK_ERR(AddFunToList(params[i]));
503 }
504 }
505
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)))));
514 } else {
515 IFPACK_CHK_ERR(-1);
516 }
517 }
518 }
519
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))));
548 } else {
549 IFPACK_CHK_ERR(-1);
550 }
551 }
552 }
553 }
554
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");
559
560 Dump_ = list.get("hypre: Dump", false);
561
562 return 0;
563} //SetParameters()
564
565//==============================================================================
566int Ifpack_Hypre::AddFunToList(RCP<FunctionParameter> NewFun){
567 NumFunsToCall_ = NumFunsToCall_+1;
568 FunsToCall_.resize(NumFunsToCall_);
569 FunsToCall_[NumFunsToCall_-1] = NewFun;
570 return 0;
571} //AddFunToList()
572
573//==============================================================================
574int Ifpack_Hypre::SetParameter(Hypre_Chooser chooser, int (*pt2Func)(HYPRE_Solver, int), int parameter){
575 RCP<FunctionParameter> temp = rcp(new FunctionParameter(chooser, pt2Func, parameter));
576 IFPACK_CHK_ERR(AddFunToList(temp));
577 return 0;
578} //SetParameter() - int function pointer
579
580//==============================================================================
581int Ifpack_Hypre::SetParameter(Hypre_Chooser chooser, int (*pt2Func)(HYPRE_Solver, double), double parameter){
582 RCP<FunctionParameter> temp = rcp(new FunctionParameter(chooser, pt2Func, parameter));
583 IFPACK_CHK_ERR(AddFunToList(temp));
584 return 0;
585} //SetParameter() - double function pointer
586
587//==============================================================================
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));
590 IFPACK_CHK_ERR(AddFunToList(temp));
591 return 0;
592} //SetParameter() - double,int function pointer
593
594//==============================================================================
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));
597 IFPACK_CHK_ERR(AddFunToList(temp));
598 return 0;
599} //SetParameter() int,int function pointer
600
601//==============================================================================
602int Ifpack_Hypre::SetParameter(Hypre_Chooser chooser, int (*pt2Func)(HYPRE_Solver, double*), double* parameter){
603 RCP<FunctionParameter> temp = rcp(new FunctionParameter(chooser, pt2Func, parameter));
604 IFPACK_CHK_ERR(AddFunToList(temp));
605 return 0;
606} //SetParameter() - double* function pointer
607
608//==============================================================================
609int Ifpack_Hypre::SetParameter(Hypre_Chooser chooser, int (*pt2Func)(HYPRE_Solver, int*), int* parameter){
610 RCP<FunctionParameter> temp = rcp(new FunctionParameter(chooser, pt2Func, parameter));
611 IFPACK_CHK_ERR(AddFunToList(temp));
612 return 0;
613} //SetParameter() - int* function pointer
614
615//==============================================================================
616int Ifpack_Hypre::SetParameter(Hypre_Chooser chooser, int (*pt2Func)(HYPRE_Solver, int**), int** parameter){
617 RCP<FunctionParameter> temp = rcp(new FunctionParameter(chooser, pt2Func, parameter));
618 IFPACK_CHK_ERR(AddFunToList(temp));
619 return 0;
620} //SetParameter() - int** function pointer
621
622//==============================================================================
623int Ifpack_Hypre::SetParameter(Hypre_Chooser chooser, Hypre_Solver solver){
624 if(chooser == Solver){
625 SolverType_ = solver;
626 } else {
627 PrecondType_ = solver;
628 }
629 return 0;
630} //SetParameter() - set type of solver
631
632//==============================================================================
633int Ifpack_Hypre::SetDiscreteGradient(Teuchos::RCP<const Epetra_CrsMatrix> G){
634
635 // Sanity check
636 if(!A_->RowMatrixRowMap().SameAs(G->RowMap()))
637 throw std::runtime_error("Ifpack_Hypre: Edge map mismatch: A and discrete gradient");
638
639 // Get the maps for the nodes (assuming the edge map from A is OK);
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);
644
645 // Start building G
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));
653 IFPACK_CHK_ERR(HYPRE_IJMatrixInitialize(HypreG_));
654
655 std::vector<int> new_indices(G->MaxNumEntries());
656 for(int i = 0; i < G->NumMyRows(); i++){
657 int numEntries;
658 double * values;
659 int *indices;
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]);
663 }
664 int GlobalRow[1];
665 GlobalRow[0] = GloballyContiguousRowMap_->GID(i);
666 IFPACK_CHK_ERR(HYPRE_IJMatrixSetValues(HypreG_, 1, &numEntries, GlobalRow, new_indices.data(), values));
667 }
668 IFPACK_CHK_ERR(HYPRE_IJMatrixAssemble(HypreG_));
669 IFPACK_CHK_ERR(HYPRE_IJMatrixGetObject(HypreG_, (void**)&ParMatrixG_));
670
671 if (Dump_)
672 HYPRE_ParCSRMatrixPrint(ParMatrixG_,"G.mat");
673
674 if(SolverType_ == AMS)
675 HYPRE_AMSSetDiscreteGradient(Solver_, ParMatrixG_);
676 if(PrecondType_ == AMS)
677 HYPRE_AMSSetDiscreteGradient(Preconditioner_, ParMatrixG_);
678
679 return 0;
680} //SetDiscreteGradient()
681
682//==============================================================================
683int Ifpack_Hypre::SetCoordinates(Teuchos::RCP<Epetra_MultiVector> coords) {
684
685 if(!G_.is_null() && !G_->DomainMap().SameAs(coords->Map()))
686 throw std::runtime_error("Ifpack_Hypre: Node map mismatch: G->DomainMap() and coords");
687
688 if(SolverType_ != AMS && PrecondType_ != AMS)
689 return 0;
690
691 double *xPtr;
692 double *yPtr;
693 double *zPtr;
694
695 IFPACK_CHK_ERR(((*coords)(0))->ExtractView(&xPtr));
696 IFPACK_CHK_ERR(((*coords)(1))->ExtractView(&yPtr));
697 IFPACK_CHK_ERR(((*coords)(2))->ExtractView(&zPtr));
698
699 MPI_Comm comm = GetMpiComm();
700 int NumEntries = coords->MyLength();
701 int * indices = GloballyContiguousNodeRowMap_->MyGlobalElements();
702
703 int ilower = GloballyContiguousNodeRowMap_->MinMyGID();
704 int iupper = GloballyContiguousNodeRowMap_->MaxMyGID();
705
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");
709 }
710
711 IFPACK_CHK_ERR(HYPRE_IJVectorCreate(comm, ilower, iupper, &xHypre_));
712 IFPACK_CHK_ERR(HYPRE_IJVectorSetObjectType(xHypre_, HYPRE_PARCSR));
713 IFPACK_CHK_ERR(HYPRE_IJVectorInitialize(xHypre_));
714 IFPACK_CHK_ERR(HYPRE_IJVectorSetValues(xHypre_,NumEntries,indices,xPtr));
715 IFPACK_CHK_ERR(HYPRE_IJVectorAssemble(xHypre_));
716 IFPACK_CHK_ERR(HYPRE_IJVectorGetObject(xHypre_, (void**) &xPar_));
717
718 IFPACK_CHK_ERR(HYPRE_IJVectorCreate(comm, ilower, iupper, &yHypre_));
719 IFPACK_CHK_ERR(HYPRE_IJVectorSetObjectType(yHypre_, HYPRE_PARCSR));
720 IFPACK_CHK_ERR(HYPRE_IJVectorInitialize(yHypre_));
721 IFPACK_CHK_ERR(HYPRE_IJVectorSetValues(yHypre_,NumEntries,indices,yPtr));
722 IFPACK_CHK_ERR(HYPRE_IJVectorAssemble(yHypre_));
723 IFPACK_CHK_ERR(HYPRE_IJVectorGetObject(yHypre_, (void**) &yPar_));
724
725 IFPACK_CHK_ERR(HYPRE_IJVectorCreate(comm, ilower, iupper, &zHypre_));
726 IFPACK_CHK_ERR(HYPRE_IJVectorSetObjectType(zHypre_, HYPRE_PARCSR));
727 IFPACK_CHK_ERR(HYPRE_IJVectorInitialize(zHypre_));
728 IFPACK_CHK_ERR(HYPRE_IJVectorSetValues(zHypre_,NumEntries,indices,zPtr));
729 IFPACK_CHK_ERR(HYPRE_IJVectorAssemble(zHypre_));
730 IFPACK_CHK_ERR(HYPRE_IJVectorGetObject(zHypre_, (void**) &zPar_));
731
732 if (Dump_) {
733 HYPRE_ParVectorPrint(xPar_,"coordX.dat");
734 HYPRE_ParVectorPrint(yPar_,"coordY.dat");
735 HYPRE_ParVectorPrint(zPar_,"coordZ.dat");
736 }
737
738 if(SolverType_ == AMS)
739 HYPRE_AMSSetCoordinateVectors(Solver_, xPar_, yPar_, zPar_);
740 if(PrecondType_ == AMS)
741 HYPRE_AMSSetCoordinateVectors(Preconditioner_, xPar_, yPar_, zPar_);
742
743 return 0;
744
745} //SetCoordinates
746
747//==============================================================================
748int Ifpack_Hypre::Compute(){
749 if(IsInitialized() == false){
750 IFPACK_CHK_ERR(Initialize());
751 }
752 Time_.ResetStartTime();
753
754 // Create the Hypre matrix and copy values. Note this uses values (which
755 // Initialize() shouldn't do) but it doesn't care what they are (for
756 // instance they can be uninitialized data even). It should be possible to
757 // set the Hypre structure without copying values, but this is the easiest
758 // way to get the structure.
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));
764 IFPACK_CHK_ERR(HYPRE_IJMatrixInitialize(HypreA_));
765 CopyEpetraToHypre();
766 if(SolveOrPrec_ == Solver) {
767 IFPACK_CHK_ERR(SetSolverType(SolverType_));
768 if (SolverPrecondPtr_ != NULL && UsePreconditioner_) {
769 // both method allows a PC (first condition) and the user wants a PC (second)
770 IFPACK_CHK_ERR(SetPrecondType(PrecondType_));
771 CallFunctions();
772 IFPACK_CHK_ERR(SolverPrecondPtr_(Solver_, PrecondSolvePtr_, PrecondSetupPtr_, Preconditioner_));
773 } else {
774 CallFunctions();
775 }
776 } else {
777 IFPACK_CHK_ERR(SetPrecondType(PrecondType_));
778 CallFunctions();
779 }
780
781 if (!G_.is_null()) {
782 SetDiscreteGradient(G_);
783 }
784
785 if (!Coords_.is_null()) {
786 SetCoordinates(Coords_);
787 }
788
789 // Hypre Setup must be called after matrix has values
790 if(SolveOrPrec_ == Solver){
791 IFPACK_CHK_ERR(SolverSetupPtr_(Solver_, ParMatrix_, ParX_, ParY_));
792 } else {
793 IFPACK_CHK_ERR(PrecondSetupPtr_(Preconditioner_, ParMatrix_, ParX_, ParY_));
794 }
795
796 // Dump Hierarchy here for BoomerAMG Preconditioner
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);
801#if 0
802 HYPRE_Int **CF_marker_array = hypre_ParAMGDataCFMarkerArray(amg_data);
803#endif
804 HYPRE_Int num_levels = hypre_ParAMGDataNumLevels(amg_data);
805
806 char ofs[80];
807 for(int k=0; k<num_levels; k++) {
808 // A
809 sprintf(ofs,"A_matrix.bmg.%d.dat",k);
810 HYPRE_ParCSRMatrixPrint(A_array[k], ofs);
811 if(k!=num_levels-1) {
812 // P
813 sprintf(ofs,"P_matrix.bmg.%d.dat",k);
814 HYPRE_ParCSRMatrixPrint(P_array[k], ofs);
815
816#if 0
817 // CF
818 // Note: Hypre outputs "-1" for F Points and "1" for C Points
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]);
826 fclose(f);
827#endif
828 }
829
830 }
831 }//end dump for BoomerAMG
832
833
834 IsComputed_ = true;
835 NumCompute_ = NumCompute_ + 1;
836 ComputeTime_ = ComputeTime_ + Time_.ElapsedTime();
837 return 0;
838} //Compute()
839
840//==============================================================================
841int Ifpack_Hypre::CallFunctions() const{
842 for(int i = 0; i < NumFunsToCall_; i++){
843 IFPACK_CHK_ERR(FunsToCall_[i]->CallFunction(Solver_, Preconditioner_));
844 }
845 return 0;
846} //CallFunctions()
847
848//==============================================================================
849int Ifpack_Hypre::ApplyInverse(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const{
850 if(IsComputed() == false){
851 IFPACK_CHK_ERR(-1);
852 }
853 Time_.ResetStartTime();
854 hypre_Vector *XLocal_ = hypre_ParVectorLocalVector(XVec_);
855 hypre_Vector *YLocal_ = hypre_ParVectorLocalVector(YVec_);
856
857 bool SameVectors = false;
858 int NumVectors = X.NumVectors();
859 if (NumVectors != Y.NumVectors()) IFPACK_CHK_ERR(-1); // X and Y must have same number of vectors
860 if(X.Pointers() == Y.Pointers() || (NumVectors == 1 && X[0] == Y[0])){
861 SameVectors = true;
862 }
863
864 // NOTE: Here were assuming that the local ordering of Epetra's X/Y-vectors and
865 // Hypre's X/Y-vectors are the same. Seeing as as this is more or less how we constructed
866 // the Hypre matrices, this seems pretty reasoanble.
867
868 for(int VecNum = 0; VecNum < NumVectors; VecNum++) {
869 //Get values for current vector in multivector.
870 // FIXME amk Nov 23, 2015: This will not work for funky data layouts
871 double * XValues = const_cast<double*>(X[VecNum]);
872 double * YValues;
873 if(!SameVectors){
874 YValues = const_cast<double*>(Y[VecNum]);
875 } else {
876 YValues = VectorCache_.getRawPtr();
877 }
878 // Temporarily make a pointer to data in Hypre for end
879 double *XTemp = XLocal_->data;
880 // Replace data in Hypre vectors with Epetra data
881 XLocal_->data = XValues;
882 double *YTemp = YLocal_->data;
883 YLocal_->data = YValues;
884
885 IFPACK_CHK_ERR(HYPRE_ParVectorSetConstantValues(ParY_, 0.0));
886 if(SolveOrPrec_ == Solver){
887 // Use the solver methods
888 IFPACK_CHK_ERR(SolverSolvePtr_(Solver_, ParMatrix_, ParX_, ParY_));
889 } else {
890 // Apply the preconditioner
891 IFPACK_CHK_ERR(PrecondSolvePtr_(Preconditioner_, ParMatrix_, ParX_, ParY_));
892 }
893 if(SameVectors){
894 int NumEntries = Y.MyLength();
895 for(int i = 0; i < NumEntries; i++)
896 Y[VecNum][i] = YValues[i];
897 }
898 XLocal_->data = XTemp;
899 YLocal_->data = YTemp;
900 }
901 NumApplyInverse_ = NumApplyInverse_ + 1;
902 ApplyInverseTime_ = ApplyInverseTime_ + Time_.ElapsedTime();
903 return 0;
904} //ApplyInverse()
905
906//==============================================================================
907int Ifpack_Hypre::Multiply(bool TransA, const Epetra_MultiVector& X, Epetra_MultiVector& Y) const{
908 if(IsInitialized() == false){
909 IFPACK_CHK_ERR(-1);
910 }
911 hypre_Vector *XLocal_ = hypre_ParVectorLocalVector(XVec_);
912 hypre_Vector *YLocal_ = hypre_ParVectorLocalVector(YVec_);
913 bool SameVectors = false;
914 int NumVectors = X.NumVectors();
915 if (NumVectors != Y.NumVectors()) IFPACK_CHK_ERR(-1); // X and Y must have same number of vectors
916 if(X.Pointers() == Y.Pointers() || (NumVectors == 1 && X[0] == Y[0])){
917 SameVectors = true;
918 }
919
920 // NOTE: Here were assuming that the local ordering of Epetra's X/Y-vectors and
921 // Hypre's X/Y-vectors are the same. Seeing as as this is more or less how we constructed
922 // the Hypre matrices, this seems pretty reasoanble.
923
924 for(int VecNum = 0; VecNum < NumVectors; VecNum++) {
925 //Get values for current vector in multivector.
926 double * XValues=const_cast<double*>(X[VecNum]);
927 double * YValues;
928 double *XTemp = XLocal_->data;
929 double *YTemp = YLocal_->data;
930 if(!SameVectors){
931 YValues = const_cast<double*>(Y[VecNum]);
932 } else {
933 YValues = VectorCache_.getRawPtr();
934 }
935 YLocal_->data = YValues;
936 IFPACK_CHK_ERR(HYPRE_ParVectorSetConstantValues(ParY_,0.0));
937 // Temporarily make a pointer to data in Hypre for end
938 // Replace data in Hypre vectors with epetra values
939 XLocal_->data = XValues;
940 // Do actual computation.
941 if(TransA) {
942 // Use transpose of A in multiply
943 IFPACK_CHK_ERR(HYPRE_ParCSRMatrixMatvecT(1.0, ParMatrix_, ParX_, 1.0, ParY_));
944 } else {
945 IFPACK_CHK_ERR(HYPRE_ParCSRMatrixMatvec(1.0, ParMatrix_, ParX_, 1.0, ParY_));
946 }
947 if(SameVectors){
948 int NumEntries = Y.MyLength();
949 for(int i = 0; i < NumEntries; i++)
950 Y[VecNum][i] = YValues[i];
951 }
952 XLocal_->data = XTemp;
953 YLocal_->data = YTemp;
954 }
955 return 0;
956} //Multiply()
957
958//==============================================================================
959std::ostream& Ifpack_Hypre::Print(std::ostream& os) const{
960 using std::endl;
961 if (!Comm().MyPID()) {
962 os << endl;
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;
969 os << 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;
980 else
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;
987 else
988 os << " " << std::setw(15) << 0.0 << endl;
989 os << "================================================================================" << endl;
990 os << endl;
991 }
992 return os;
993} //Print()
994
995//==============================================================================
996double Ifpack_Hypre::Condest(const Ifpack_CondestType CT,
997 const int MaxIters,
998 const double Tol,
999 Epetra_RowMatrix* Matrix_in){
1000 if (!IsComputed()) // cannot compute right now
1001 return(-1.0);
1002 Condest_ = Ifpack_Condest(*this, CT, MaxIters, Tol, Matrix_in);
1003 return(Condest_);
1004} //Condest()
1005
1006//==============================================================================
1007int Ifpack_Hypre::SetSolverType(Hypre_Solver Solver){
1008 switch(Solver) {
1009 case BoomerAMG:
1010 if(IsSolverCreated_){
1011 SolverDestroyPtr_(Solver_);
1012 IsSolverCreated_ = false;
1013 }
1014 SolverCreatePtr_ = &Ifpack_Hypre::Hypre_BoomerAMGCreate;
1015 SolverDestroyPtr_ = &HYPRE_BoomerAMGDestroy;
1016 SolverSetupPtr_ = &HYPRE_BoomerAMGSetup;
1017 SolverPrecondPtr_ = NULL;
1018 SolverSolvePtr_ = &HYPRE_BoomerAMGSolve;
1019 break;
1020 case AMS:
1021 if(IsSolverCreated_){
1022 SolverDestroyPtr_(Solver_);
1023 IsSolverCreated_ = false;
1024 }
1025 SolverCreatePtr_ = &Ifpack_Hypre::Hypre_AMSCreate;
1026 SolverDestroyPtr_ = &HYPRE_AMSDestroy;
1027 SolverSetupPtr_ = &HYPRE_AMSSetup;
1028 SolverSolvePtr_ = &HYPRE_AMSSolve;
1029 SolverPrecondPtr_ = NULL;
1030 break;
1031 case Hybrid:
1032 if(IsSolverCreated_){
1033 SolverDestroyPtr_(Solver_);
1034 IsSolverCreated_ = false;
1035 }
1036 SolverCreatePtr_ = &Ifpack_Hypre::Hypre_ParCSRHybridCreate;
1037 SolverDestroyPtr_ = &HYPRE_ParCSRHybridDestroy;
1038 SolverSetupPtr_ = &HYPRE_ParCSRHybridSetup;
1039 SolverSolvePtr_ = &HYPRE_ParCSRHybridSolve;
1040 SolverPrecondPtr_ = &HYPRE_ParCSRHybridSetPrecond;
1041 break;
1042 case PCG:
1043 if(IsSolverCreated_){
1044 SolverDestroyPtr_(Solver_);
1045 IsSolverCreated_ = false;
1046 }
1047 SolverCreatePtr_ = &Ifpack_Hypre::Hypre_ParCSRPCGCreate;
1048 SolverDestroyPtr_ = &HYPRE_ParCSRPCGDestroy;
1049 SolverSetupPtr_ = &HYPRE_ParCSRPCGSetup;
1050 SolverSolvePtr_ = &HYPRE_ParCSRPCGSolve;
1051 SolverPrecondPtr_ = &HYPRE_ParCSRPCGSetPrecond;
1052 break;
1053 case GMRES:
1054 if(IsSolverCreated_){
1055 SolverDestroyPtr_(Solver_);
1056 IsSolverCreated_ = false;
1057 }
1058 SolverCreatePtr_ = &Ifpack_Hypre::Hypre_ParCSRGMRESCreate;
1059 SolverDestroyPtr_ = &HYPRE_ParCSRGMRESDestroy;
1060 SolverSetupPtr_ = &HYPRE_ParCSRGMRESSetup;
1061 SolverPrecondPtr_ = &HYPRE_ParCSRGMRESSetPrecond;
1062 break;
1063 case FlexGMRES:
1064 if(IsSolverCreated_){
1065 SolverDestroyPtr_(Solver_);
1066 IsSolverCreated_ = false;
1067 }
1068 SolverCreatePtr_ = &Ifpack_Hypre::Hypre_ParCSRFlexGMRESCreate;
1069 SolverDestroyPtr_ = &HYPRE_ParCSRFlexGMRESDestroy;
1070 SolverSetupPtr_ = &HYPRE_ParCSRFlexGMRESSetup;
1071 SolverSolvePtr_ = &HYPRE_ParCSRFlexGMRESSolve;
1072 SolverPrecondPtr_ = &HYPRE_ParCSRFlexGMRESSetPrecond;
1073 break;
1074 case LGMRES:
1075 if(IsSolverCreated_){
1076 SolverDestroyPtr_(Solver_);
1077 IsSolverCreated_ = false;
1078 }
1079 SolverCreatePtr_ = &Ifpack_Hypre::Hypre_ParCSRLGMRESCreate;
1080 SolverDestroyPtr_ = &HYPRE_ParCSRLGMRESDestroy;
1081 SolverSetupPtr_ = &HYPRE_ParCSRLGMRESSetup;
1082 SolverSolvePtr_ = &HYPRE_ParCSRLGMRESSolve;
1083 SolverPrecondPtr_ = &HYPRE_ParCSRLGMRESSetPrecond;
1084 break;
1085 case BiCGSTAB:
1086 if(IsSolverCreated_){
1087 SolverDestroyPtr_(Solver_);
1088 IsSolverCreated_ = false;
1089 }
1090 SolverCreatePtr_ = &Ifpack_Hypre::Hypre_ParCSRBiCGSTABCreate;
1091 SolverDestroyPtr_ = &HYPRE_ParCSRBiCGSTABDestroy;
1092 SolverSetupPtr_ = &HYPRE_ParCSRBiCGSTABSetup;
1093 SolverSolvePtr_ = &HYPRE_ParCSRBiCGSTABSolve;
1094 SolverPrecondPtr_ = &HYPRE_ParCSRBiCGSTABSetPrecond;
1095 break;
1096 default:
1097 return -1;
1098 }
1099 CreateSolver();
1100 return 0;
1101} //SetSolverType()
1102
1103//==============================================================================
1104int Ifpack_Hypre::SetPrecondType(Hypre_Solver Precond){
1105 switch(Precond) {
1106 case BoomerAMG:
1107 if(IsPrecondCreated_){
1108 PrecondDestroyPtr_(Preconditioner_);
1109 IsPrecondCreated_ = false;
1110 }
1111 PrecondCreatePtr_ = &Ifpack_Hypre::Hypre_BoomerAMGCreate;
1112 PrecondDestroyPtr_ = &HYPRE_BoomerAMGDestroy;
1113 PrecondSetupPtr_ = &HYPRE_BoomerAMGSetup;
1114 PrecondSolvePtr_ = &HYPRE_BoomerAMGSolve;
1115 break;
1116 case ParaSails:
1117 if(IsPrecondCreated_){
1118 PrecondDestroyPtr_(Preconditioner_);
1119 IsPrecondCreated_ = false;
1120 }
1121 PrecondCreatePtr_ = &Ifpack_Hypre::Hypre_ParaSailsCreate;
1122 PrecondDestroyPtr_ = &HYPRE_ParaSailsDestroy;
1123 PrecondSetupPtr_ = &HYPRE_ParaSailsSetup;
1124 PrecondSolvePtr_ = &HYPRE_ParaSailsSolve;
1125 break;
1126 case Euclid:
1127 if(IsPrecondCreated_){
1128 PrecondDestroyPtr_(Preconditioner_);
1129 IsPrecondCreated_ = false;
1130 }
1131 PrecondCreatePtr_ = &Ifpack_Hypre::Hypre_EuclidCreate;
1132 PrecondDestroyPtr_ = &HYPRE_EuclidDestroy;
1133 PrecondSetupPtr_ = &HYPRE_EuclidSetup;
1134 PrecondSolvePtr_ = &HYPRE_EuclidSolve;
1135 break;
1136 case AMS:
1137 if(IsPrecondCreated_){
1138 PrecondDestroyPtr_(Preconditioner_);
1139 IsPrecondCreated_ = false;
1140 }
1141 PrecondCreatePtr_ = &Ifpack_Hypre::Hypre_AMSCreate;
1142 PrecondDestroyPtr_ = &HYPRE_AMSDestroy;
1143 PrecondSetupPtr_ = &HYPRE_AMSSetup;
1144 PrecondSolvePtr_ = &HYPRE_AMSSolve;
1145 break;
1146 default:
1147 return -1;
1148 }
1149 CreatePrecond();
1150 return 0;
1151
1152} //SetPrecondType()
1153
1154//==============================================================================
1155int Ifpack_Hypre::CreateSolver(){
1156 MPI_Comm comm;
1157 HYPRE_ParCSRMatrixGetComm(ParMatrix_, &comm);
1158 int ierr = (this->*SolverCreatePtr_)(comm, &Solver_);
1159 IsSolverCreated_ = true;
1160 return ierr;
1161} //CreateSolver()
1162
1163//==============================================================================
1164int Ifpack_Hypre::CreatePrecond(){
1165 MPI_Comm comm;
1166 HYPRE_ParCSRMatrixGetComm(ParMatrix_, &comm);
1167 int ierr = (this->*PrecondCreatePtr_)(comm, &Preconditioner_);
1168 IsPrecondCreated_ = true;
1169 return ierr;
1170} //CreatePrecond()
1171
1172
1173//==============================================================================
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");
1178
1179 std::vector<int> new_indices(Matrix->MaxNumEntries());
1180 for(int i = 0; i < Matrix->NumMyRows(); i++){
1181 int numEntries;
1182 int *indices;
1183 double *values;
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]);
1187 }
1188 int GlobalRow[1];
1189 GlobalRow[0] = GloballyContiguousRowMap_->GID(i);
1190 IFPACK_CHK_ERR(HYPRE_IJMatrixSetValues(HypreA_, 1, &numEntries, GlobalRow, new_indices.data(), values));
1191 }
1192 IFPACK_CHK_ERR(HYPRE_IJMatrixAssemble(HypreA_));
1193 IFPACK_CHK_ERR(HYPRE_IJMatrixGetObject(HypreA_, (void**)&ParMatrix_));
1194 if (Dump_)
1195 HYPRE_ParCSRMatrixPrint(ParMatrix_,"A.mat");
1196 return 0;
1197} //CopyEpetraToHypre()
1198
1199//==============================================================================
1200int Ifpack_Hypre::Hypre_BoomerAMGCreate(MPI_Comm /*comm*/, HYPRE_Solver *solver)
1201 { return HYPRE_BoomerAMGCreate(solver);}
1202
1203//==============================================================================
1204int Ifpack_Hypre::Hypre_ParaSailsCreate(MPI_Comm comm, HYPRE_Solver *solver)
1205 { return HYPRE_ParaSailsCreate(comm, solver);}
1206
1207//==============================================================================
1208int Ifpack_Hypre::Hypre_EuclidCreate(MPI_Comm comm, HYPRE_Solver *solver)
1209 { return HYPRE_EuclidCreate(comm, solver);}
1210
1211//==============================================================================
1212int Ifpack_Hypre::Hypre_AMSCreate(MPI_Comm /*comm*/, HYPRE_Solver *solver)
1213 { return HYPRE_AMSCreate(solver);}
1214
1215//==============================================================================
1216int Ifpack_Hypre::Hypre_ParCSRHybridCreate(MPI_Comm /*comm*/, HYPRE_Solver *solver)
1217 { return HYPRE_ParCSRHybridCreate(solver);}
1218
1219//==============================================================================
1220int Ifpack_Hypre::Hypre_ParCSRPCGCreate(MPI_Comm comm, HYPRE_Solver *solver)
1221 { return HYPRE_ParCSRPCGCreate(comm, solver);}
1222
1223//==============================================================================
1224int Ifpack_Hypre::Hypre_ParCSRGMRESCreate(MPI_Comm comm, HYPRE_Solver *solver)
1225 { return HYPRE_ParCSRGMRESCreate(comm, solver);}
1226
1227//==============================================================================
1228int Ifpack_Hypre::Hypre_ParCSRFlexGMRESCreate(MPI_Comm comm, HYPRE_Solver *solver)
1229 { return HYPRE_ParCSRFlexGMRESCreate(comm, solver);}
1230
1231//==============================================================================
1232int Ifpack_Hypre::Hypre_ParCSRLGMRESCreate(MPI_Comm comm, HYPRE_Solver *solver)
1233 { return HYPRE_ParCSRLGMRESCreate(comm, solver);}
1234
1235//==============================================================================
1236int Ifpack_Hypre::Hypre_ParCSRBiCGSTABCreate(MPI_Comm comm, HYPRE_Solver *solver)
1237 { return HYPRE_ParCSRBiCGSTABCreate(comm, solver);}
1238
1239//==============================================================================
1240Teuchos::RCP<const Epetra_Map> Ifpack_Hypre::MakeContiguousColumnMap(Teuchos::RCP<const Epetra_RowMatrix> &MatrixRow) const{
1241 // Must create GloballyContiguous DomainMap (which is a permutation of Matrix_'s
1242 // DomainMap) and the corresponding permuted ColumnMap.
1243 // Epetra_GID ---------> LID ----------> HYPRE_GID
1244 // via DomainMap.LID() via GloballyContiguousDomainMap.GID()
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();
1250 const Epetra_Import * importer = Matrix->Importer();
1251
1252 if(DomainMap.LinearMap() ) {
1253 // If the domain map is linear, then we can just use the column map as is.
1254 return rcpFromRef(ColumnMap);
1255 }
1256 else {
1257 // The domain map isn't linear, so we need a new domain map
1258 Teuchos::RCP<Epetra_Map> ContiguousDomainMap = rcp(new Epetra_Map(DomainMap.NumGlobalElements(),
1259 DomainMap.NumMyElements(), 0, Comm()));
1260 if(importer) {
1261 // If there's an importer then we can use it to get a new column map
1262 Epetra_IntVector MyGIDsHYPRE(View,DomainMap,ContiguousDomainMap->MyGlobalElements());
1263
1264 // import the HYPRE GIDs
1265 Epetra_IntVector ColGIDsHYPRE(ColumnMap);
1266 ColGIDsHYPRE.Import(MyGIDsHYPRE, *importer, Insert);
1267
1268 // Make a HYPRE numbering-based column map.
1269 return Teuchos::rcp(new Epetra_Map(ColumnMap.NumGlobalElements(),ColGIDsHYPRE.MyLength(), &ColGIDsHYPRE[0], 0, Comm()));
1270 }
1271 else {
1272 // The problem has matching domain/column maps, and somehow the domain map isn't linear, so just use the new domain map
1273 return Teuchos::rcp(new Epetra_Map(ColumnMap.NumGlobalElements(),ColumnMap.NumMyElements(), ContiguousDomainMap->MyGlobalElements(), 0, Comm()));
1274 }
1275 }
1276}
1277
1278
1279
1280#endif // HAVE_HYPRE && HAVE_MPI
bool LinearMap() const
int NumGlobalElements() const
int NumMyElements() const
int NumVectors() const
int MyLength() const
double ** Pointers() const
virtual const Epetra_Map & RowMatrixRowMap() const=0