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_double_func)(HYPRE_Solver, int, double);
70typedef int (*int_int_func)(HYPRE_Solver, int, int);
71typedef int (*int_star_func)(HYPRE_Solver, int*);
72typedef int (*int_star_star_func)(HYPRE_Solver, int**);
73typedef int (*double_star_func)(HYPRE_Solver, double*);
74typedef int (*int_int_double_double_func)(HYPRE_Solver, int, int, double, double);
75typedef int (*int_int_int_double_int_int_func)(HYPRE_Solver, int, int, int, double, int, int);
76typedef int (*char_star_func)(HYPRE_Solver, char*);
77
79class FunctionParameter{
80 public:
82 FunctionParameter(Hypre_Chooser chooser, int_func funct, int param1) :
83 chooser_(chooser),
84 option_(0),
85 int_func_(funct),
86 int_param1_(param1) {}
87
88 FunctionParameter(Hypre_Chooser chooser, std::string funct_name, int param1) :
89 chooser_(chooser),
90 option_(0),
91 int_func_(hypreMapIntFunc_.at(funct_name)),
92 int_param1_(param1) {}
93
95 FunctionParameter(Hypre_Chooser chooser, double_func funct, double param1):
96 chooser_(chooser),
97 option_(1),
98 double_func_(funct),
99 double_param1_(param1) {}
100
101 FunctionParameter(Hypre_Chooser chooser, std::string funct_name, double param1):
102 chooser_(chooser),
103 option_(1),
104 double_func_(hypreMapDoubleFunc_.at(funct_name)),
105 double_param1_(param1) {}
106
108 FunctionParameter(Hypre_Chooser chooser, double_int_func funct, double param1, int param2):
109 chooser_(chooser),
110 option_(2),
111 double_int_func_(funct),
112 int_param1_(param2),
113 double_param1_(param1) {}
114
116 FunctionParameter(Hypre_Chooser chooser, int_double_func funct, int param1, double param2):
117 chooser_(chooser),
118 option_(10),
119 int_double_func_(funct),
120 int_param1_(param1),
121 double_param1_(param2) {}
122
123 FunctionParameter(Hypre_Chooser chooser, std::string funct_name, double param1, int param2):
124 chooser_(chooser),
125 option_(2),
126 double_int_func_(hypreMapDoubleIntFunc_.at(funct_name)),
127 int_param1_(param2),
128 double_param1_(param1) {}
129
131 FunctionParameter(Hypre_Chooser chooser, int_int_func funct, int param1, int param2):
132 chooser_(chooser),
133 option_(3),
134 int_int_func_(funct),
135 int_param1_(param1),
136 int_param2_(param2) {}
137
138 FunctionParameter(Hypre_Chooser chooser, std::string funct_name, int param1, int param2):
139 chooser_(chooser),
140 option_(3),
141 int_int_func_(hypreMapIntIntFunc_.at(funct_name)),
142 int_param1_(param1),
143 int_param2_(param2) {}
144
146 FunctionParameter(Hypre_Chooser chooser, int_star_func funct, int *param1):
147 chooser_(chooser),
148 option_(4),
149 int_star_func_(funct),
150 int_star_param_(param1) {}
151
152 FunctionParameter(Hypre_Chooser chooser, std::string funct_name, int *param1):
153 chooser_(chooser),
154 option_(4),
155 int_star_func_(hypreMapIntStarFunc_.at(funct_name)),
156 int_star_param_(param1) {}
157
159 FunctionParameter(Hypre_Chooser chooser, double_star_func funct, double* param1):
160 chooser_(chooser),
161 option_(5),
162 double_star_func_(funct),
163 double_star_param_(param1) {}
164
165 FunctionParameter(Hypre_Chooser chooser, std::string funct_name, double* param1):
166 chooser_(chooser),
167 option_(5),
168 double_star_func_(hypreMapDoubleStarFunc_.at(funct_name)),
169 double_star_param_(param1) {}
170
172 FunctionParameter(Hypre_Chooser chooser, int_int_double_double_func funct, int param1, int param2, double param3, double param4):
173 chooser_(chooser),
174 option_(6),
175 int_int_double_double_func_(funct),
176 int_param1_(param1),
177 int_param2_(param2),
178 double_param1_(param3),
179 double_param2_(param4) {}
180
181 FunctionParameter(Hypre_Chooser chooser, std::string funct_name, int param1, int param2, double param3, double param4):
182 chooser_(chooser),
183 option_(6),
184 int_int_double_double_func_(hypreMapIntIntDoubleDoubleFunc_.at(funct_name)),
185 int_param1_(param1),
186 int_param2_(param2),
187 double_param1_(param3),
188 double_param2_(param4) {}
189
191 FunctionParameter(Hypre_Chooser chooser, int_star_star_func funct, int ** param1):
192 chooser_(chooser),
193 option_(7),
194 int_star_star_func_(funct),
195 int_star_star_param_(param1) {}
196
197 FunctionParameter(Hypre_Chooser chooser, std::string funct_name, int** param1):
198 chooser_(chooser),
199 option_(7),
200 int_star_star_func_(hypreMapIntStarStarFunc_.at(funct_name)),
201 int_star_star_param_(param1) {}
202
204 FunctionParameter(Hypre_Chooser chooser, int_int_int_double_int_int_func funct, int param1, int param2, int param3, double param4, int param5, int param6):
205 chooser_(chooser),
206 option_(8),
207 int_int_int_double_int_int_func_(funct),
208 int_param1_(param1),
209 int_param2_(param2),
210 int_param3_(param3),
211 int_param4_(param5),
212 int_param5_(param6),
213 double_param1_(param4) {}
214
215 FunctionParameter(Hypre_Chooser chooser, std::string funct_name, int param1, int param2, int param3, double param4, int param5, int param6):
216 chooser_(chooser),
217 option_(8),
218 int_int_int_double_int_int_func_(hypreMapIntIntIntDoubleIntIntFunc_.at(funct_name)),
219 int_param1_(param1),
220 int_param2_(param2),
221 int_param3_(param3),
222 int_param4_(param5),
223 int_param5_(param6),
224 double_param1_(param4) {}
225
227 FunctionParameter(Hypre_Chooser chooser, char_star_func funct, char *param1):
228 chooser_(chooser),
229 option_(9),
230 char_star_func_(funct),
231 char_star_param_(param1) {}
232
233 FunctionParameter(Hypre_Chooser chooser, std::string funct_name, char *param1):
234 chooser_(chooser),
235 option_(9),
236 char_star_func_(hypreMapCharStarFunc_.at(funct_name)),
237 char_star_param_(param1) {}
238
240 int CallFunction(HYPRE_Solver solver, HYPRE_Solver precond) {
241 if(chooser_ == Solver){
242 if(option_ == 0){
243 return int_func_(solver, int_param1_);
244 } else if(option_ == 1){
245 return double_func_(solver, double_param1_);
246 } else if(option_ == 2){
247 return double_int_func_(solver, double_param1_, int_param1_);
248 } else if (option_ == 3){
249 return int_int_func_(solver, int_param1_, int_param2_);
250 } else if (option_ == 4){
251 return int_star_func_(solver, int_star_param_);
252 } else if (option_ == 5){
253 return double_star_func_(solver, double_star_param_);
254 } else if (option_ == 6) {
255 return int_int_double_double_func_(solver, int_param1_, int_param2_, double_param1_, double_param2_);
256 } else if (option_ == 7) {
257 return int_star_star_func_(solver, int_star_star_param_);
258 } else if (option_ == 8) {
259 return int_int_int_double_int_int_func_(solver, int_param1_, int_param2_, int_param3_, double_param1_, int_param4_, int_param5_);
260 } else if (option_ == 9) {
261 return char_star_func_(solver, char_star_param_);
262 } else if (option_ == 10) {
263 return int_double_func_(solver, int_param1_, double_param1_);
264 } else {
265 IFPACK_CHK_ERR(-2);
266 }
267 } else {
268 if(option_ == 0){
269 return int_func_(precond, int_param1_);
270 } else if(option_ == 1){
271 return double_func_(precond, double_param1_);
272 } else if(option_ == 2){
273 return double_int_func_(precond, double_param1_, int_param1_);
274 } else if(option_ == 3) {
275 return int_int_func_(precond, int_param1_, int_param2_);
276 } else if(option_ == 4) {
277 return int_star_func_(precond, int_star_param_);
278 } else if(option_ == 5) {
279 return double_star_func_(precond, double_star_param_);
280 } else if (option_ == 6) {
281 return int_int_double_double_func_(precond, int_param1_, int_param2_, double_param1_, double_param2_);
282 } else if (option_ == 7) {
283 return int_star_star_func_(precond, int_star_star_param_);
284 } else if (option_ == 8) {
285 return int_int_int_double_int_int_func_(precond, int_param1_, int_param2_, int_param3_, double_param1_, int_param4_, int_param5_);
286 } else if (option_ == 9) {
287 return char_star_func_(solver, char_star_param_);
288 } else if (option_ == 10) {
289 return int_double_func_(precond, int_param1_, double_param1_);
290 } else {
291 IFPACK_CHK_ERR(-2);
292 }
293 }
294 }
295
296 static bool isFuncIntInt(std::string funct_name) {
297 return (hypreMapIntIntFunc_.find(funct_name) != hypreMapIntIntFunc_.end());
298 }
299
300 static bool isFuncIntIntDoubleDouble(std::string funct_name) {
301 return (hypreMapIntIntDoubleDoubleFunc_.find(funct_name) != hypreMapIntIntDoubleDoubleFunc_.end());
302 }
303
304 static bool isFuncIntDouble(std::string funct_name) {
305 return (hypreMapIntDoubleFunc_.find(funct_name) != hypreMapIntDoubleFunc_.end());
306 }
307
308 static bool isFuncIntIntIntDoubleIntInt(std::string funct_name) {
309 return (hypreMapIntIntIntDoubleIntIntFunc_.find(funct_name) != hypreMapIntIntIntDoubleIntIntFunc_.end());
310 }
311
312 static bool isFuncIntStarStar(std::string funct_name) {
313 return (hypreMapIntStarStarFunc_.find(funct_name) != hypreMapIntStarStarFunc_.end());
314 }
315
316 private:
317 Hypre_Chooser chooser_;
318 int option_;
319 int_func int_func_;
320 double_func double_func_;
321 double_int_func double_int_func_;
322 int_double_func int_double_func_;
323 int_int_func int_int_func_;
324 int_star_func int_star_func_;
325 double_star_func double_star_func_;
326 int_int_double_double_func int_int_double_double_func_;
327 int_int_int_double_int_int_func int_int_int_double_int_int_func_;
328 int_star_star_func int_star_star_func_;
329 char_star_func char_star_func_;
330 int int_param1_;
331 int int_param2_;
332 int int_param3_;
333 int int_param4_;
334 int int_param5_;
335 double double_param1_;
336 double double_param2_;
337 int *int_star_param_;
338 int **int_star_star_param_;
339 double *double_star_param_;
340 char *char_star_param_;
341
342 static const std::map<std::string, int_func> hypreMapIntFunc_;
343 static const std::map<std::string, double_func> hypreMapDoubleFunc_;
344 static const std::map<std::string, double_int_func> hypreMapDoubleIntFunc_;
345 static const std::map<std::string, int_double_func> hypreMapIntDoubleFunc_;
346 static const std::map<std::string, int_int_func> hypreMapIntIntFunc_;
347 static const std::map<std::string, int_star_func> hypreMapIntStarFunc_;
348 static const std::map<std::string, double_star_func> hypreMapDoubleStarFunc_;
349 static const std::map<std::string, int_int_double_double_func> hypreMapIntIntDoubleDoubleFunc_;
350 static const std::map<std::string, int_int_int_double_int_int_func> hypreMapIntIntIntDoubleIntIntFunc_;
351 static const std::map<std::string, int_star_star_func> hypreMapIntStarStarFunc_;
352 static const std::map<std::string, char_star_func> hypreMapCharStarFunc_;
353
354};
355
356// NOTE: This really, really needs to be here and not up above, so please don't move it
357#include "Ifpack_HypreParameterMap.h"
358
359Ifpack_Hypre::Ifpack_Hypre(Epetra_RowMatrix* A):
360 A_(rcp(A,false)),
361 UseTranspose_(false),
362 IsInitialized_(false),
363 IsComputed_(false),
364 Label_(),
365 NumInitialize_(0),
366 NumCompute_(0),
367 NumApplyInverse_(0),
368 InitializeTime_(0.0),
369 ComputeTime_(0.0),
370 ApplyInverseTime_(0.0),
371 ComputeFlops_(0.0),
372 ApplyInverseFlops_(0.0),
373 Time_(A_->Comm()),
374 HypreA_(0),
375 HypreG_(0),
376 xHypre_(0),
377 yHypre_(0),
378 zHypre_(0),
379 IsSolverCreated_(false),
380 IsPrecondCreated_(false),
381 SolveOrPrec_(Solver),
382 NumFunsToCall_(0),
383 SolverType_(PCG),
384 PrecondType_(Euclid),
385 UsePreconditioner_(false),
386 Dump_(false)
387{
388 MPI_Comm comm = GetMpiComm();
389 // Check that RowMap and RangeMap are the same. While this could handle the
390 // case where RowMap and RangeMap are permutations, other Ifpack PCs don't
391 // handle this either.
392 if (!A_->RowMatrixRowMap().SameAs(A_->OperatorRangeMap())) {
393 IFPACK_CHK_ERRV(-1);
394 }
395 // Hypre expects the RowMap to be Linear.
396 if (A_->RowMatrixRowMap().LinearMap()) {
397 // note these are non-owning pointers, they are deleted by A_'s destructor
398 GloballyContiguousRowMap_ = rcpFromRef(A_->RowMatrixRowMap());
399 GloballyContiguousColMap_ = rcpFromRef(A_->RowMatrixColMap());
400 } else {
401 // Must create GloballyContiguous Maps for Hypre
402 if(A_->OperatorDomainMap().SameAs(A_->RowMatrixRowMap())) {
403 Teuchos::RCP<const Epetra_RowMatrix> Aconst = A_;
404 GloballyContiguousColMap_ = MakeContiguousColumnMap(Aconst);
405 GloballyContiguousRowMap_ = rcp(new Epetra_Map(A_->RowMatrixRowMap().NumGlobalElements(),
406 A_->RowMatrixRowMap().NumMyElements(), 0, Comm()));
407 }
408 else {
409 throw std::runtime_error("Ifpack_Hypre: Unsupported map configuration: Row/Domain maps do not match");
410 }
411 }
412 // Next create vectors that will be used when ApplyInverse() is called
413 int ilower = GloballyContiguousRowMap_->MinMyGID();
414 int iupper = GloballyContiguousRowMap_->MaxMyGID();
415 // X in AX = Y
416 IFPACK_CHK_ERRV(HYPRE_IJVectorCreate(comm, ilower, iupper, &XHypre_));
417 IFPACK_CHK_ERRV(HYPRE_IJVectorSetObjectType(XHypre_, HYPRE_PARCSR));
418 IFPACK_CHK_ERRV(HYPRE_IJVectorInitialize(XHypre_));
419 IFPACK_CHK_ERRV(HYPRE_IJVectorAssemble(XHypre_));
420 IFPACK_CHK_ERRV(HYPRE_IJVectorGetObject(XHypre_, (void**) &ParX_));
421 XVec_ = Teuchos::rcp((hypre_ParVector *) hypre_IJVectorObject(((hypre_IJVector *) XHypre_)),false);
422
423 // Y in AX = Y
424 IFPACK_CHK_ERRV(HYPRE_IJVectorCreate(comm, ilower, iupper, &YHypre_));
425 IFPACK_CHK_ERRV(HYPRE_IJVectorSetObjectType(YHypre_, HYPRE_PARCSR));
426 IFPACK_CHK_ERRV(HYPRE_IJVectorInitialize(YHypre_));
427 IFPACK_CHK_ERRV(HYPRE_IJVectorAssemble(YHypre_));
428 IFPACK_CHK_ERRV(HYPRE_IJVectorGetObject(YHypre_, (void**) &ParY_));
429 YVec_ = Teuchos::rcp((hypre_ParVector *) hypre_IJVectorObject(((hypre_IJVector *) YHypre_)),false);
430
431 // Cache
432 VectorCache_.resize(A->RowMatrixRowMap().NumMyElements());
433} //Constructor
434
435//==============================================================================
436void Ifpack_Hypre::Destroy(){
437 if(IsInitialized()){
438 IFPACK_CHK_ERRV(HYPRE_IJMatrixDestroy(HypreA_));
439 }
440 IFPACK_CHK_ERRV(HYPRE_IJVectorDestroy(XHypre_));
441 IFPACK_CHK_ERRV(HYPRE_IJVectorDestroy(YHypre_));
442 if(IsSolverCreated_){
443 IFPACK_CHK_ERRV(SolverDestroyPtr_(Solver_));
444 }
445 if(IsPrecondCreated_){
446 IFPACK_CHK_ERRV(PrecondDestroyPtr_(Preconditioner_));
447 }
448
449 // Maxwell
450 if(HypreG_) {
451 IFPACK_CHK_ERRV(HYPRE_IJMatrixDestroy(HypreG_));
452 }
453 if(xHypre_) {
454 IFPACK_CHK_ERRV(HYPRE_IJVectorDestroy(xHypre_));
455 }
456 if(yHypre_) {
457 IFPACK_CHK_ERRV(HYPRE_IJVectorDestroy(yHypre_));
458 }
459 if(zHypre_) {
460 IFPACK_CHK_ERRV(HYPRE_IJVectorDestroy(zHypre_));
461 }
462} //Destroy()
463
464//==============================================================================
465int Ifpack_Hypre::Initialize(){
466 Time_.ResetStartTime();
467 if(IsInitialized_) return 0;
468 // set flags
469 IsInitialized_=true;
470 NumInitialize_ = NumInitialize_ + 1;
471 InitializeTime_ = InitializeTime_ + Time_.ElapsedTime();
472 return 0;
473} //Initialize()
474
475//==============================================================================
476int Ifpack_Hypre::SetParameters(Teuchos::ParameterList& list){
477
478 std::map<std::string, Hypre_Solver> solverMap;
479 solverMap["BoomerAMG"] = BoomerAMG;
480 solverMap["ParaSails"] = ParaSails;
481 solverMap["Euclid"] = Euclid;
482 solverMap["AMS"] = AMS;
483 solverMap["Hybrid"] = Hybrid;
484 solverMap["PCG"] = PCG;
485 solverMap["GMRES"] = GMRES;
486 solverMap["FlexGMRES"] = FlexGMRES;
487 solverMap["LGMRES"] = LGMRES;
488 solverMap["BiCGSTAB"] = BiCGSTAB;
489
490 std::map<std::string, Hypre_Chooser> chooserMap;
491 chooserMap["Solver"] = Solver;
492 chooserMap["Preconditioner"] = Preconditioner;
493
494 List_ = list;
495 Hypre_Solver solType;
496 if (list.isType<std::string>("hypre: Solver"))
497 solType = solverMap[list.get<std::string>("hypre: Solver")];
498 else
499 solType = list.get("hypre: Solver", PCG);
500 SolverType_ = solType;
501 Hypre_Solver precType;
502 if (list.isType<std::string>("hypre: Preconditioner"))
503 precType = solverMap[list.get<std::string>("hypre: Preconditioner")];
504 else
505 precType = list.get("hypre: Preconditioner", Euclid);
506 PrecondType_ = precType;
507 Hypre_Chooser chooser;
508 if (list.isType<std::string>("hypre: SolveOrPrecondition"))
509 chooser = chooserMap[list.get<std::string>("hypre: SolveOrPrecondition")];
510 else
511 chooser = list.get("hypre: SolveOrPrecondition", Solver);
512 SolveOrPrec_ = chooser;
513 bool SetPrecond = list.get("hypre: SetPreconditioner", false);
514 IFPACK_CHK_ERR(SetParameter(SetPrecond));
515 int NumFunctions = list.get("hypre: NumFunctions", 0);
516 FunsToCall_.clear();
517 NumFunsToCall_ = 0;
518 if(NumFunctions > 0){
519 RCP<FunctionParameter>* params = list.get<RCP<FunctionParameter>*>("hypre: Functions");
520 for(int i = 0; i < NumFunctions; i++){
521 IFPACK_CHK_ERR(AddFunToList(params[i]));
522 }
523 }
524
525 if (list.isSublist("hypre: Solver functions")) {
526 Teuchos::ParameterList solverList = list.sublist("hypre: Solver functions");
527 for (auto it = solverList.begin(); it != solverList.end(); ++it) {
528 std::string funct_name = it->first;
529 if (it->second.isType<int>()) {
530 IFPACK_CHK_ERR(AddFunToList(rcp(new FunctionParameter(Solver, funct_name , Teuchos::getValue<int>(it->second)))));
531 } else if (it->second.isType<double>()) {
532 IFPACK_CHK_ERR(AddFunToList(rcp(new FunctionParameter(Solver, funct_name , Teuchos::getValue<double>(it->second)))));
533 } else {
534 IFPACK_CHK_ERR(-1);
535 }
536 }
537 }
538
539 if (list.isSublist("hypre: Preconditioner functions")) {
540 Teuchos::ParameterList precList = list.sublist("hypre: Preconditioner functions");
541 for (auto it = precList.begin(); it != precList.end(); ++it) {
542 std::string funct_name = it->first;
543 if (it->second.isType<int>()) {
544 IFPACK_CHK_ERR(AddFunToList(rcp(new FunctionParameter(Preconditioner, funct_name , Teuchos::getValue<int>(it->second)))));
545 } else if (it->second.isType<double>()) {
546 IFPACK_CHK_ERR(AddFunToList(rcp(new FunctionParameter(Preconditioner, funct_name , Teuchos::getValue<double>(it->second)))));
547 } else if (it->second.isList()) {
548 Teuchos::ParameterList pl = Teuchos::getValue<Teuchos::ParameterList>(it->second);
549 if (FunctionParameter::isFuncIntInt(funct_name)) {
550 int arg0 = pl.get<int>("arg 0");
551 int arg1 = pl.get<int>("arg 1");
552 IFPACK_CHK_ERR(AddFunToList(rcp(new FunctionParameter(Preconditioner, funct_name , arg0, arg1))));
553 } else if (FunctionParameter::isFuncIntDouble(funct_name)) {
554 int arg0 = pl.get<int>("arg 0");
555 double arg1 = pl.get<double>("arg 1");
556 IFPACK_CHK_ERR(AddFunToList(rcp(new FunctionParameter(Preconditioner, funct_name , arg0, arg1))));
557 } else if (FunctionParameter::isFuncIntIntDoubleDouble(funct_name)) {
558 int arg0 = pl.get<int>("arg 0");
559 int arg1 = pl.get<int>("arg 1");
560 double arg2 = pl.get<double>("arg 2");
561 double arg3 = pl.get<double>("arg 3");
562 IFPACK_CHK_ERR(AddFunToList(rcp(new FunctionParameter(Preconditioner, funct_name , arg0, arg1, arg2, arg3))));
563 } else if (FunctionParameter::isFuncIntIntIntDoubleIntInt(funct_name)) {
564 int arg0 = pl.get<int>("arg 0");
565 int arg1 = pl.get<int>("arg 1");
566 int arg2 = pl.get<int>("arg 2");
567 double arg3 = pl.get<double>("arg 3");
568 int arg4 = pl.get<int>("arg 4");
569 int arg5 = pl.get<int>("arg 5");
570 IFPACK_CHK_ERR(AddFunToList(rcp(new FunctionParameter(Preconditioner, funct_name , arg0, arg1, arg2, arg3, arg4, arg5))));
571 } else {
572 IFPACK_CHK_ERR(-1);
573 }
574 }
575 }
576 }
577
578 if (list.isSublist("Coordinates") && list.sublist("Coordinates").isType<Teuchos::RCP<Epetra_MultiVector> >("Coordinates"))
579 Coords_ = list.sublist("Coordinates").get<Teuchos::RCP<Epetra_MultiVector> >("Coordinates");
580 if (list.isSublist("Operators") && list.sublist("Operators").isType<Teuchos::RCP<const Epetra_CrsMatrix> >("G"))
581 G_ = list.sublist("Operators").get<Teuchos::RCP<const Epetra_CrsMatrix> >("G");
582
583 Dump_ = list.get("hypre: Dump", false);
584
585 return 0;
586} //SetParameters()
587
588//==============================================================================
589int Ifpack_Hypre::AddFunToList(RCP<FunctionParameter> NewFun){
590 NumFunsToCall_ = NumFunsToCall_+1;
591 FunsToCall_.resize(NumFunsToCall_);
592 FunsToCall_[NumFunsToCall_-1] = NewFun;
593 return 0;
594} //AddFunToList()
595
596//==============================================================================
597int Ifpack_Hypre::SetParameter(Hypre_Chooser chooser, int (*pt2Func)(HYPRE_Solver, int), int parameter){
598 RCP<FunctionParameter> temp = rcp(new FunctionParameter(chooser, pt2Func, parameter));
599 IFPACK_CHK_ERR(AddFunToList(temp));
600 return 0;
601} //SetParameter() - int function pointer
602
603//==============================================================================
604int Ifpack_Hypre::SetParameter(Hypre_Chooser chooser, int (*pt2Func)(HYPRE_Solver, double), double parameter){
605 RCP<FunctionParameter> temp = rcp(new FunctionParameter(chooser, pt2Func, parameter));
606 IFPACK_CHK_ERR(AddFunToList(temp));
607 return 0;
608} //SetParameter() - double function pointer
609
610//==============================================================================
611int Ifpack_Hypre::SetParameter(Hypre_Chooser chooser, int (*pt2Func)(HYPRE_Solver, double, int), double parameter1, int parameter2){
612 RCP<FunctionParameter> temp = rcp(new FunctionParameter(chooser, pt2Func, parameter1, parameter2));
613 IFPACK_CHK_ERR(AddFunToList(temp));
614 return 0;
615} //SetParameter() - double,int function pointer
616
617//==============================================================================
618int Ifpack_Hypre::SetParameter(Hypre_Chooser chooser, int (*pt2Func)(HYPRE_Solver, int, double), int parameter1, double parameter2){
619 RCP<FunctionParameter> temp = rcp(new FunctionParameter(chooser, pt2Func, parameter1, parameter2));
620 IFPACK_CHK_ERR(AddFunToList(temp));
621 return 0;
622} //SetParameter() - int,double function pointer
623
624//==============================================================================
625int Ifpack_Hypre::SetParameter(Hypre_Chooser chooser, int (*pt2Func)(HYPRE_Solver, int, int), int parameter1, int parameter2){
626 RCP<FunctionParameter> temp = rcp(new FunctionParameter(chooser, pt2Func, parameter1, parameter2));
627 IFPACK_CHK_ERR(AddFunToList(temp));
628 return 0;
629} //SetParameter() int,int function pointer
630
631//==============================================================================
632int Ifpack_Hypre::SetParameter(Hypre_Chooser chooser, int (*pt2Func)(HYPRE_Solver, double*), double* parameter){
633 RCP<FunctionParameter> temp = rcp(new FunctionParameter(chooser, pt2Func, parameter));
634 IFPACK_CHK_ERR(AddFunToList(temp));
635 return 0;
636} //SetParameter() - double* function pointer
637
638//==============================================================================
639int Ifpack_Hypre::SetParameter(Hypre_Chooser chooser, int (*pt2Func)(HYPRE_Solver, int*), int* parameter){
640 RCP<FunctionParameter> temp = rcp(new FunctionParameter(chooser, pt2Func, parameter));
641 IFPACK_CHK_ERR(AddFunToList(temp));
642 return 0;
643} //SetParameter() - int* function pointer
644
645//==============================================================================
646int Ifpack_Hypre::SetParameter(Hypre_Chooser chooser, int (*pt2Func)(HYPRE_Solver, int**), int** parameter){
647 RCP<FunctionParameter> temp = rcp(new FunctionParameter(chooser, pt2Func, parameter));
648 IFPACK_CHK_ERR(AddFunToList(temp));
649 return 0;
650} //SetParameter() - int** function pointer
651
652//==============================================================================
653int Ifpack_Hypre::SetParameter(Hypre_Chooser chooser, Hypre_Solver solver){
654 if(chooser == Solver){
655 SolverType_ = solver;
656 } else {
657 PrecondType_ = solver;
658 }
659 return 0;
660} //SetParameter() - set type of solver
661
662//==============================================================================
663int Ifpack_Hypre::SetDiscreteGradient(Teuchos::RCP<const Epetra_CrsMatrix> G){
664
665 // Sanity check
666 if(!A_->RowMatrixRowMap().SameAs(G->RowMap()))
667 throw std::runtime_error("Ifpack_Hypre: Edge map mismatch: A and discrete gradient");
668
669 // Get the maps for the nodes (assuming the edge map from A is OK);
670 GloballyContiguousNodeRowMap_ = rcp(new Epetra_Map(G->DomainMap().NumGlobalElements(),
671 G->DomainMap().NumMyElements(), 0, Comm()));
672 Teuchos::RCP<const Epetra_RowMatrix> Grow = Teuchos::rcp_dynamic_cast<const Epetra_RowMatrix>(G);
673 GloballyContiguousNodeColMap_ = MakeContiguousColumnMap(Grow);
674
675 // Start building G
676 MPI_Comm comm = GetMpiComm();
677 int ilower = GloballyContiguousRowMap_->MinMyGID();
678 int iupper = GloballyContiguousRowMap_->MaxMyGID();
679 int jlower = GloballyContiguousNodeRowMap_->MinMyGID();
680 int jupper = GloballyContiguousNodeRowMap_->MaxMyGID();
681 IFPACK_CHK_ERR(HYPRE_IJMatrixCreate(comm, ilower, iupper, jlower, jupper, &HypreG_));
682 IFPACK_CHK_ERR(HYPRE_IJMatrixSetObjectType(HypreG_, HYPRE_PARCSR));
683 IFPACK_CHK_ERR(HYPRE_IJMatrixInitialize(HypreG_));
684
685 std::vector<int> new_indices(G->MaxNumEntries());
686 for(int i = 0; i < G->NumMyRows(); i++){
687 int numEntries;
688 double * values;
689 int *indices;
690 IFPACK_CHK_ERR(G->ExtractMyRowView(i, numEntries, values, indices));
691 for(int j = 0; j < numEntries; j++){
692 new_indices[j] = GloballyContiguousNodeColMap_->GID(indices[j]);
693 }
694 int GlobalRow[1];
695 GlobalRow[0] = GloballyContiguousRowMap_->GID(i);
696 IFPACK_CHK_ERR(HYPRE_IJMatrixSetValues(HypreG_, 1, &numEntries, GlobalRow, new_indices.data(), values));
697 }
698 IFPACK_CHK_ERR(HYPRE_IJMatrixAssemble(HypreG_));
699 IFPACK_CHK_ERR(HYPRE_IJMatrixGetObject(HypreG_, (void**)&ParMatrixG_));
700
701 if (Dump_)
702 HYPRE_ParCSRMatrixPrint(ParMatrixG_,"G.mat");
703
704 if(SolverType_ == AMS)
705 HYPRE_AMSSetDiscreteGradient(Solver_, ParMatrixG_);
706 if(PrecondType_ == AMS)
707 HYPRE_AMSSetDiscreteGradient(Preconditioner_, ParMatrixG_);
708
709 return 0;
710} //SetDiscreteGradient()
711
712//==============================================================================
713int Ifpack_Hypre::SetCoordinates(Teuchos::RCP<Epetra_MultiVector> coords) {
714
715 if(!G_.is_null() && !G_->DomainMap().SameAs(coords->Map()))
716 throw std::runtime_error("Ifpack_Hypre: Node map mismatch: G->DomainMap() and coords");
717
718 if(SolverType_ != AMS && PrecondType_ != AMS)
719 return 0;
720
721 double *xPtr;
722 double *yPtr;
723 double *zPtr;
724
725 IFPACK_CHK_ERR(((*coords)(0))->ExtractView(&xPtr));
726 IFPACK_CHK_ERR(((*coords)(1))->ExtractView(&yPtr));
727 IFPACK_CHK_ERR(((*coords)(2))->ExtractView(&zPtr));
728
729 MPI_Comm comm = GetMpiComm();
730 int NumEntries = coords->MyLength();
731 int * indices = GloballyContiguousNodeRowMap_->MyGlobalElements();
732
733 int ilower = GloballyContiguousNodeRowMap_->MinMyGID();
734 int iupper = GloballyContiguousNodeRowMap_->MaxMyGID();
735
736 if( NumEntries != iupper-ilower+1) {
737 std::cout<<"Ifpack_Hypre::SetCoordinates(): Error on rank "<<Comm().MyPID()<<": MyLength = "<<coords->MyLength()<<" GID range = ["<<ilower<<","<<iupper<<"]"<<std::endl;
738 throw std::runtime_error("Ifpack_Hypre: SetCoordinates: Length mismatch");
739 }
740
741 IFPACK_CHK_ERR(HYPRE_IJVectorCreate(comm, ilower, iupper, &xHypre_));
742 IFPACK_CHK_ERR(HYPRE_IJVectorSetObjectType(xHypre_, HYPRE_PARCSR));
743 IFPACK_CHK_ERR(HYPRE_IJVectorInitialize(xHypre_));
744 IFPACK_CHK_ERR(HYPRE_IJVectorSetValues(xHypre_,NumEntries,indices,xPtr));
745 IFPACK_CHK_ERR(HYPRE_IJVectorAssemble(xHypre_));
746 IFPACK_CHK_ERR(HYPRE_IJVectorGetObject(xHypre_, (void**) &xPar_));
747
748 IFPACK_CHK_ERR(HYPRE_IJVectorCreate(comm, ilower, iupper, &yHypre_));
749 IFPACK_CHK_ERR(HYPRE_IJVectorSetObjectType(yHypre_, HYPRE_PARCSR));
750 IFPACK_CHK_ERR(HYPRE_IJVectorInitialize(yHypre_));
751 IFPACK_CHK_ERR(HYPRE_IJVectorSetValues(yHypre_,NumEntries,indices,yPtr));
752 IFPACK_CHK_ERR(HYPRE_IJVectorAssemble(yHypre_));
753 IFPACK_CHK_ERR(HYPRE_IJVectorGetObject(yHypre_, (void**) &yPar_));
754
755 IFPACK_CHK_ERR(HYPRE_IJVectorCreate(comm, ilower, iupper, &zHypre_));
756 IFPACK_CHK_ERR(HYPRE_IJVectorSetObjectType(zHypre_, HYPRE_PARCSR));
757 IFPACK_CHK_ERR(HYPRE_IJVectorInitialize(zHypre_));
758 IFPACK_CHK_ERR(HYPRE_IJVectorSetValues(zHypre_,NumEntries,indices,zPtr));
759 IFPACK_CHK_ERR(HYPRE_IJVectorAssemble(zHypre_));
760 IFPACK_CHK_ERR(HYPRE_IJVectorGetObject(zHypre_, (void**) &zPar_));
761
762 if (Dump_) {
763 HYPRE_ParVectorPrint(xPar_,"coordX.dat");
764 HYPRE_ParVectorPrint(yPar_,"coordY.dat");
765 HYPRE_ParVectorPrint(zPar_,"coordZ.dat");
766 }
767
768 if(SolverType_ == AMS)
769 HYPRE_AMSSetCoordinateVectors(Solver_, xPar_, yPar_, zPar_);
770 if(PrecondType_ == AMS)
771 HYPRE_AMSSetCoordinateVectors(Preconditioner_, xPar_, yPar_, zPar_);
772
773 return 0;
774
775} //SetCoordinates
776
777//==============================================================================
778int Ifpack_Hypre::Compute(){
779 if(IsInitialized() == false){
780 IFPACK_CHK_ERR(Initialize());
781 }
782 Time_.ResetStartTime();
783
784 // Create the Hypre matrix and copy values. Note this uses values (which
785 // Initialize() shouldn't do) but it doesn't care what they are (for
786 // instance they can be uninitialized data even). It should be possible to
787 // set the Hypre structure without copying values, but this is the easiest
788 // way to get the structure.
789 MPI_Comm comm = GetMpiComm();
790 int ilower = GloballyContiguousRowMap_->MinMyGID();
791 int iupper = GloballyContiguousRowMap_->MaxMyGID();
792 IFPACK_CHK_ERR(HYPRE_IJMatrixCreate(comm, ilower, iupper, ilower, iupper, &HypreA_));
793 IFPACK_CHK_ERR(HYPRE_IJMatrixSetObjectType(HypreA_, HYPRE_PARCSR));
794 IFPACK_CHK_ERR(HYPRE_IJMatrixInitialize(HypreA_));
795 CopyEpetraToHypre();
796 if(SolveOrPrec_ == Solver) {
797 IFPACK_CHK_ERR(SetSolverType(SolverType_));
798 if (SolverPrecondPtr_ != NULL && UsePreconditioner_) {
799 // both method allows a PC (first condition) and the user wants a PC (second)
800 IFPACK_CHK_ERR(SetPrecondType(PrecondType_));
801 CallFunctions();
802 IFPACK_CHK_ERR(SolverPrecondPtr_(Solver_, PrecondSolvePtr_, PrecondSetupPtr_, Preconditioner_));
803 } else {
804 CallFunctions();
805 }
806 } else {
807 IFPACK_CHK_ERR(SetPrecondType(PrecondType_));
808 CallFunctions();
809 }
810
811 if (!G_.is_null()) {
812 SetDiscreteGradient(G_);
813 }
814
815 if (!Coords_.is_null()) {
816 SetCoordinates(Coords_);
817 }
818
819 // Hypre Setup must be called after matrix has values
820 if(SolveOrPrec_ == Solver){
821 IFPACK_CHK_ERR(SolverSetupPtr_(Solver_, ParMatrix_, ParX_, ParY_));
822 } else {
823 IFPACK_CHK_ERR(PrecondSetupPtr_(Preconditioner_, ParMatrix_, ParX_, ParY_));
824 }
825
826 // Dump Hierarchy here for BoomerAMG Preconditioner
827 if(Dump_ && PrecondSolvePtr_ == &HYPRE_BoomerAMGSolve) {
828 hypre_ParAMGData *amg_data = (hypre_ParAMGData*) Preconditioner_;
829 hypre_ParCSRMatrix **A_array = hypre_ParAMGDataAArray(amg_data);
830 hypre_ParCSRMatrix **P_array = hypre_ParAMGDataPArray(amg_data);
831#if 0
832 HYPRE_Int **CF_marker_array = hypre_ParAMGDataCFMarkerArray(amg_data);
833#endif
834 HYPRE_Int num_levels = hypre_ParAMGDataNumLevels(amg_data);
835
836 char ofs[80];
837 for(int k=0; k<num_levels; k++) {
838 // A
839 sprintf(ofs,"A_matrix.bmg.%d.dat",k);
840 HYPRE_ParCSRMatrixPrint(A_array[k], ofs);
841 if(k!=num_levels-1) {
842 // P
843 sprintf(ofs,"P_matrix.bmg.%d.dat",k);
844 HYPRE_ParCSRMatrixPrint(P_array[k], ofs);
845
846#if 0
847 // CF
848 // Note: Hypre outputs "-1" for F Points and "1" for C Points
849 HYPRE_Int local_size = hypre_CSRMatrixNumRows(hypre_ParCSRMatrixDiag(A_array[k]));
850 sprintf(ofs,"cf_marker.bmg.%d.dat.%d",k,Comm().MyPID());
851 FILE * f = fopen(ofs,"w");
852 fprintf(f,"%%%%MatrixMarket matrix array real general\n");
853 fprintf(f,"%d 1\n",local_size);
854 for(int i=0; i<local_size; i++)
855 fprintf(f,"%d\n",(int)CF_marker_array[k][i]);
856 fclose(f);
857#endif
858 }
859
860 }
861 }//end dump for BoomerAMG
862
863
864 IsComputed_ = true;
865 NumCompute_ = NumCompute_ + 1;
866 ComputeTime_ = ComputeTime_ + Time_.ElapsedTime();
867 return 0;
868} //Compute()
869
870//==============================================================================
871int Ifpack_Hypre::CallFunctions() const{
872 for(int i = 0; i < NumFunsToCall_; i++){
873 IFPACK_CHK_ERR(FunsToCall_[i]->CallFunction(Solver_, Preconditioner_));
874 }
875 return 0;
876} //CallFunctions()
877
878//==============================================================================
879int Ifpack_Hypre::ApplyInverse(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const{
880 if(IsComputed() == false){
881 IFPACK_CHK_ERR(-1);
882 }
883 Time_.ResetStartTime();
884 hypre_Vector *XLocal_ = hypre_ParVectorLocalVector(XVec_);
885 hypre_Vector *YLocal_ = hypre_ParVectorLocalVector(YVec_);
886
887 bool SameVectors = false;
888 int NumVectors = X.NumVectors();
889 if (NumVectors != Y.NumVectors()) IFPACK_CHK_ERR(-1); // X and Y must have same number of vectors
890 if(X.Pointers() == Y.Pointers() || (NumVectors == 1 && X[0] == Y[0])){
891 SameVectors = true;
892 }
893
894 // NOTE: Here were assuming that the local ordering of Epetra's X/Y-vectors and
895 // Hypre's X/Y-vectors are the same. Seeing as as this is more or less how we constructed
896 // the Hypre matrices, this seems pretty reasoanble.
897
898 for(int VecNum = 0; VecNum < NumVectors; VecNum++) {
899 //Get values for current vector in multivector.
900 // FIXME amk Nov 23, 2015: This will not work for funky data layouts
901 double * XValues = const_cast<double*>(X[VecNum]);
902 double * YValues;
903 if(!SameVectors){
904 YValues = const_cast<double*>(Y[VecNum]);
905 } else {
906 YValues = VectorCache_.getRawPtr();
907 }
908 // Temporarily make a pointer to data in Hypre for end
909 double *XTemp = XLocal_->data;
910 // Replace data in Hypre vectors with Epetra data
911 XLocal_->data = XValues;
912 double *YTemp = YLocal_->data;
913 YLocal_->data = YValues;
914
915 IFPACK_CHK_ERR(HYPRE_ParVectorSetConstantValues(ParY_, 0.0));
916 if(SolveOrPrec_ == Solver){
917 // Use the solver methods
918 IFPACK_CHK_ERR(SolverSolvePtr_(Solver_, ParMatrix_, ParX_, ParY_));
919 } else {
920 // Apply the preconditioner
921 IFPACK_CHK_ERR(PrecondSolvePtr_(Preconditioner_, ParMatrix_, ParX_, ParY_));
922 }
923 if(SameVectors){
924 int NumEntries = Y.MyLength();
925 for(int i = 0; i < NumEntries; i++)
926 Y[VecNum][i] = YValues[i];
927 }
928 XLocal_->data = XTemp;
929 YLocal_->data = YTemp;
930 }
931 NumApplyInverse_ = NumApplyInverse_ + 1;
932 ApplyInverseTime_ = ApplyInverseTime_ + Time_.ElapsedTime();
933 return 0;
934} //ApplyInverse()
935
936//==============================================================================
937int Ifpack_Hypre::Multiply(bool TransA, const Epetra_MultiVector& X, Epetra_MultiVector& Y) const{
938 if(IsInitialized() == false){
939 IFPACK_CHK_ERR(-1);
940 }
941 hypre_Vector *XLocal_ = hypre_ParVectorLocalVector(XVec_);
942 hypre_Vector *YLocal_ = hypre_ParVectorLocalVector(YVec_);
943 bool SameVectors = false;
944 int NumVectors = X.NumVectors();
945 if (NumVectors != Y.NumVectors()) IFPACK_CHK_ERR(-1); // X and Y must have same number of vectors
946 if(X.Pointers() == Y.Pointers() || (NumVectors == 1 && X[0] == Y[0])){
947 SameVectors = true;
948 }
949
950 // NOTE: Here were assuming that the local ordering of Epetra's X/Y-vectors and
951 // Hypre's X/Y-vectors are the same. Seeing as as this is more or less how we constructed
952 // the Hypre matrices, this seems pretty reasoanble.
953
954 for(int VecNum = 0; VecNum < NumVectors; VecNum++) {
955 //Get values for current vector in multivector.
956 double * XValues=const_cast<double*>(X[VecNum]);
957 double * YValues;
958 double *XTemp = XLocal_->data;
959 double *YTemp = YLocal_->data;
960 if(!SameVectors){
961 YValues = const_cast<double*>(Y[VecNum]);
962 } else {
963 YValues = VectorCache_.getRawPtr();
964 }
965 YLocal_->data = YValues;
966 IFPACK_CHK_ERR(HYPRE_ParVectorSetConstantValues(ParY_,0.0));
967 // Temporarily make a pointer to data in Hypre for end
968 // Replace data in Hypre vectors with epetra values
969 XLocal_->data = XValues;
970 // Do actual computation.
971 if(TransA) {
972 // Use transpose of A in multiply
973 IFPACK_CHK_ERR(HYPRE_ParCSRMatrixMatvecT(1.0, ParMatrix_, ParX_, 1.0, ParY_));
974 } else {
975 IFPACK_CHK_ERR(HYPRE_ParCSRMatrixMatvec(1.0, ParMatrix_, ParX_, 1.0, ParY_));
976 }
977 if(SameVectors){
978 int NumEntries = Y.MyLength();
979 for(int i = 0; i < NumEntries; i++)
980 Y[VecNum][i] = YValues[i];
981 }
982 XLocal_->data = XTemp;
983 YLocal_->data = YTemp;
984 }
985 return 0;
986} //Multiply()
987
988//==============================================================================
989std::ostream& Ifpack_Hypre::Print(std::ostream& os) const{
990 using std::endl;
991 if (!Comm().MyPID()) {
992 os << endl;
993 os << "================================================================================" << endl;
994 os << "Ifpack_Hypre: " << Label() << endl << endl;
995 os << "Using " << Comm().NumProc() << " processors." << endl;
996 os << "Global number of rows = " << A_->NumGlobalRows() << endl;
997 os << "Global number of nonzeros = " << A_->NumGlobalNonzeros() << endl;
998 os << "Condition number estimate = " << Condest() << endl;
999 os << endl;
1000 os << "Phase # calls Total Time (s) Total MFlops MFlops/s" << endl;
1001 os << "----- ------- -------------- ------------ --------" << endl;
1002 os << "Initialize() " << std::setw(5) << NumInitialize_
1003 << " " << std::setw(15) << InitializeTime_
1004 << " 0.0 0.0" << endl;
1005 os << "Compute() " << std::setw(5) << NumCompute_
1006 << " " << std::setw(15) << ComputeTime_
1007 << " " << std::setw(15) << 1.0e-6 * ComputeFlops_;
1008 if (ComputeTime_ != 0.0)
1009 os << " " << std::setw(15) << 1.0e-6 * ComputeFlops_ / ComputeTime_ << endl;
1010 else
1011 os << " " << std::setw(15) << 0.0 << endl;
1012 os << "ApplyInverse() " << std::setw(5) << NumApplyInverse_
1013 << " " << std::setw(15) << ApplyInverseTime_
1014 << " " << std::setw(15) << 1.0e-6 * ApplyInverseFlops_;
1015 if (ApplyInverseTime_ != 0.0)
1016 os << " " << std::setw(15) << 1.0e-6 * ApplyInverseFlops_ / ApplyInverseTime_ << endl;
1017 else
1018 os << " " << std::setw(15) << 0.0 << endl;
1019 os << "================================================================================" << endl;
1020 os << endl;
1021 }
1022 return os;
1023} //Print()
1024
1025//==============================================================================
1026double Ifpack_Hypre::Condest(const Ifpack_CondestType CT,
1027 const int MaxIters,
1028 const double Tol,
1029 Epetra_RowMatrix* Matrix_in){
1030 if (!IsComputed()) // cannot compute right now
1031 return(-1.0);
1032 Condest_ = Ifpack_Condest(*this, CT, MaxIters, Tol, Matrix_in);
1033 return(Condest_);
1034} //Condest()
1035
1036//==============================================================================
1037int Ifpack_Hypre::SetSolverType(Hypre_Solver Solver){
1038 switch(Solver) {
1039 case BoomerAMG:
1040 if(IsSolverCreated_){
1041 SolverDestroyPtr_(Solver_);
1042 IsSolverCreated_ = false;
1043 }
1044 SolverCreatePtr_ = &Ifpack_Hypre::Hypre_BoomerAMGCreate;
1045 SolverDestroyPtr_ = &HYPRE_BoomerAMGDestroy;
1046 SolverSetupPtr_ = &HYPRE_BoomerAMGSetup;
1047 SolverPrecondPtr_ = NULL;
1048 SolverSolvePtr_ = &HYPRE_BoomerAMGSolve;
1049 break;
1050 case AMS:
1051 if(IsSolverCreated_){
1052 SolverDestroyPtr_(Solver_);
1053 IsSolverCreated_ = false;
1054 }
1055 SolverCreatePtr_ = &Ifpack_Hypre::Hypre_AMSCreate;
1056 SolverDestroyPtr_ = &HYPRE_AMSDestroy;
1057 SolverSetupPtr_ = &HYPRE_AMSSetup;
1058 SolverSolvePtr_ = &HYPRE_AMSSolve;
1059 SolverPrecondPtr_ = NULL;
1060 break;
1061 case Hybrid:
1062 if(IsSolverCreated_){
1063 SolverDestroyPtr_(Solver_);
1064 IsSolverCreated_ = false;
1065 }
1066 SolverCreatePtr_ = &Ifpack_Hypre::Hypre_ParCSRHybridCreate;
1067 SolverDestroyPtr_ = &HYPRE_ParCSRHybridDestroy;
1068 SolverSetupPtr_ = &HYPRE_ParCSRHybridSetup;
1069 SolverSolvePtr_ = &HYPRE_ParCSRHybridSolve;
1070 SolverPrecondPtr_ = &HYPRE_ParCSRHybridSetPrecond;
1071 break;
1072 case PCG:
1073 if(IsSolverCreated_){
1074 SolverDestroyPtr_(Solver_);
1075 IsSolverCreated_ = false;
1076 }
1077 SolverCreatePtr_ = &Ifpack_Hypre::Hypre_ParCSRPCGCreate;
1078 SolverDestroyPtr_ = &HYPRE_ParCSRPCGDestroy;
1079 SolverSetupPtr_ = &HYPRE_ParCSRPCGSetup;
1080 SolverSolvePtr_ = &HYPRE_ParCSRPCGSolve;
1081 SolverPrecondPtr_ = &HYPRE_ParCSRPCGSetPrecond;
1082 break;
1083 case GMRES:
1084 if(IsSolverCreated_){
1085 SolverDestroyPtr_(Solver_);
1086 IsSolverCreated_ = false;
1087 }
1088 SolverCreatePtr_ = &Ifpack_Hypre::Hypre_ParCSRGMRESCreate;
1089 SolverDestroyPtr_ = &HYPRE_ParCSRGMRESDestroy;
1090 SolverSetupPtr_ = &HYPRE_ParCSRGMRESSetup;
1091 SolverPrecondPtr_ = &HYPRE_ParCSRGMRESSetPrecond;
1092 break;
1093 case FlexGMRES:
1094 if(IsSolverCreated_){
1095 SolverDestroyPtr_(Solver_);
1096 IsSolverCreated_ = false;
1097 }
1098 SolverCreatePtr_ = &Ifpack_Hypre::Hypre_ParCSRFlexGMRESCreate;
1099 SolverDestroyPtr_ = &HYPRE_ParCSRFlexGMRESDestroy;
1100 SolverSetupPtr_ = &HYPRE_ParCSRFlexGMRESSetup;
1101 SolverSolvePtr_ = &HYPRE_ParCSRFlexGMRESSolve;
1102 SolverPrecondPtr_ = &HYPRE_ParCSRFlexGMRESSetPrecond;
1103 break;
1104 case LGMRES:
1105 if(IsSolverCreated_){
1106 SolverDestroyPtr_(Solver_);
1107 IsSolverCreated_ = false;
1108 }
1109 SolverCreatePtr_ = &Ifpack_Hypre::Hypre_ParCSRLGMRESCreate;
1110 SolverDestroyPtr_ = &HYPRE_ParCSRLGMRESDestroy;
1111 SolverSetupPtr_ = &HYPRE_ParCSRLGMRESSetup;
1112 SolverSolvePtr_ = &HYPRE_ParCSRLGMRESSolve;
1113 SolverPrecondPtr_ = &HYPRE_ParCSRLGMRESSetPrecond;
1114 break;
1115 case BiCGSTAB:
1116 if(IsSolverCreated_){
1117 SolverDestroyPtr_(Solver_);
1118 IsSolverCreated_ = false;
1119 }
1120 SolverCreatePtr_ = &Ifpack_Hypre::Hypre_ParCSRBiCGSTABCreate;
1121 SolverDestroyPtr_ = &HYPRE_ParCSRBiCGSTABDestroy;
1122 SolverSetupPtr_ = &HYPRE_ParCSRBiCGSTABSetup;
1123 SolverSolvePtr_ = &HYPRE_ParCSRBiCGSTABSolve;
1124 SolverPrecondPtr_ = &HYPRE_ParCSRBiCGSTABSetPrecond;
1125 break;
1126 default:
1127 return -1;
1128 }
1129 CreateSolver();
1130 return 0;
1131} //SetSolverType()
1132
1133//==============================================================================
1134int Ifpack_Hypre::SetPrecondType(Hypre_Solver Precond){
1135 switch(Precond) {
1136 case BoomerAMG:
1137 if(IsPrecondCreated_){
1138 PrecondDestroyPtr_(Preconditioner_);
1139 IsPrecondCreated_ = false;
1140 }
1141 PrecondCreatePtr_ = &Ifpack_Hypre::Hypre_BoomerAMGCreate;
1142 PrecondDestroyPtr_ = &HYPRE_BoomerAMGDestroy;
1143 PrecondSetupPtr_ = &HYPRE_BoomerAMGSetup;
1144 PrecondSolvePtr_ = &HYPRE_BoomerAMGSolve;
1145 break;
1146 case ParaSails:
1147 if(IsPrecondCreated_){
1148 PrecondDestroyPtr_(Preconditioner_);
1149 IsPrecondCreated_ = false;
1150 }
1151 PrecondCreatePtr_ = &Ifpack_Hypre::Hypre_ParaSailsCreate;
1152 PrecondDestroyPtr_ = &HYPRE_ParaSailsDestroy;
1153 PrecondSetupPtr_ = &HYPRE_ParaSailsSetup;
1154 PrecondSolvePtr_ = &HYPRE_ParaSailsSolve;
1155 break;
1156 case Euclid:
1157 if(IsPrecondCreated_){
1158 PrecondDestroyPtr_(Preconditioner_);
1159 IsPrecondCreated_ = false;
1160 }
1161 PrecondCreatePtr_ = &Ifpack_Hypre::Hypre_EuclidCreate;
1162 PrecondDestroyPtr_ = &HYPRE_EuclidDestroy;
1163 PrecondSetupPtr_ = &HYPRE_EuclidSetup;
1164 PrecondSolvePtr_ = &HYPRE_EuclidSolve;
1165 break;
1166 case AMS:
1167 if(IsPrecondCreated_){
1168 PrecondDestroyPtr_(Preconditioner_);
1169 IsPrecondCreated_ = false;
1170 }
1171 PrecondCreatePtr_ = &Ifpack_Hypre::Hypre_AMSCreate;
1172 PrecondDestroyPtr_ = &HYPRE_AMSDestroy;
1173 PrecondSetupPtr_ = &HYPRE_AMSSetup;
1174 PrecondSolvePtr_ = &HYPRE_AMSSolve;
1175 break;
1176 default:
1177 return -1;
1178 }
1179 CreatePrecond();
1180 return 0;
1181
1182} //SetPrecondType()
1183
1184//==============================================================================
1185int Ifpack_Hypre::CreateSolver(){
1186 MPI_Comm comm;
1187 HYPRE_ParCSRMatrixGetComm(ParMatrix_, &comm);
1188 int ierr = (this->*SolverCreatePtr_)(comm, &Solver_);
1189 IsSolverCreated_ = true;
1190 return ierr;
1191} //CreateSolver()
1192
1193//==============================================================================
1194int Ifpack_Hypre::CreatePrecond(){
1195 MPI_Comm comm;
1196 HYPRE_ParCSRMatrixGetComm(ParMatrix_, &comm);
1197 int ierr = (this->*PrecondCreatePtr_)(comm, &Preconditioner_);
1198 IsPrecondCreated_ = true;
1199 return ierr;
1200} //CreatePrecond()
1201
1202
1203//==============================================================================
1204int Ifpack_Hypre::CopyEpetraToHypre(){
1205 Teuchos::RCP<const Epetra_CrsMatrix> Matrix = Teuchos::rcp_dynamic_cast<const Epetra_CrsMatrix>(A_);
1206 if(Matrix.is_null())
1207 throw std::runtime_error("Ifpack_Hypre: Unsupported matrix configuration: Epetra_CrsMatrix required");
1208
1209 std::vector<int> new_indices(Matrix->MaxNumEntries());
1210 for(int i = 0; i < Matrix->NumMyRows(); i++){
1211 int numEntries;
1212 int *indices;
1213 double *values;
1214 IFPACK_CHK_ERR(Matrix->ExtractMyRowView(i, numEntries, values, indices));
1215 for(int j = 0; j < numEntries; j++){
1216 new_indices[j] = GloballyContiguousColMap_->GID(indices[j]);
1217 }
1218 int GlobalRow[1];
1219 GlobalRow[0] = GloballyContiguousRowMap_->GID(i);
1220 IFPACK_CHK_ERR(HYPRE_IJMatrixSetValues(HypreA_, 1, &numEntries, GlobalRow, new_indices.data(), values));
1221 }
1222 IFPACK_CHK_ERR(HYPRE_IJMatrixAssemble(HypreA_));
1223 IFPACK_CHK_ERR(HYPRE_IJMatrixGetObject(HypreA_, (void**)&ParMatrix_));
1224 if (Dump_)
1225 HYPRE_ParCSRMatrixPrint(ParMatrix_,"A.mat");
1226 return 0;
1227} //CopyEpetraToHypre()
1228
1229//==============================================================================
1230int Ifpack_Hypre::Hypre_BoomerAMGCreate(MPI_Comm /*comm*/, HYPRE_Solver *solver)
1231 { return HYPRE_BoomerAMGCreate(solver);}
1232
1233//==============================================================================
1234int Ifpack_Hypre::Hypre_ParaSailsCreate(MPI_Comm comm, HYPRE_Solver *solver)
1235 { return HYPRE_ParaSailsCreate(comm, solver);}
1236
1237//==============================================================================
1238int Ifpack_Hypre::Hypre_EuclidCreate(MPI_Comm comm, HYPRE_Solver *solver)
1239 { return HYPRE_EuclidCreate(comm, solver);}
1240
1241//==============================================================================
1242int Ifpack_Hypre::Hypre_AMSCreate(MPI_Comm /*comm*/, HYPRE_Solver *solver)
1243 { return HYPRE_AMSCreate(solver);}
1244
1245//==============================================================================
1246int Ifpack_Hypre::Hypre_ParCSRHybridCreate(MPI_Comm /*comm*/, HYPRE_Solver *solver)
1247 { return HYPRE_ParCSRHybridCreate(solver);}
1248
1249//==============================================================================
1250int Ifpack_Hypre::Hypre_ParCSRPCGCreate(MPI_Comm comm, HYPRE_Solver *solver)
1251 { return HYPRE_ParCSRPCGCreate(comm, solver);}
1252
1253//==============================================================================
1254int Ifpack_Hypre::Hypre_ParCSRGMRESCreate(MPI_Comm comm, HYPRE_Solver *solver)
1255 { return HYPRE_ParCSRGMRESCreate(comm, solver);}
1256
1257//==============================================================================
1258int Ifpack_Hypre::Hypre_ParCSRFlexGMRESCreate(MPI_Comm comm, HYPRE_Solver *solver)
1259 { return HYPRE_ParCSRFlexGMRESCreate(comm, solver);}
1260
1261//==============================================================================
1262int Ifpack_Hypre::Hypre_ParCSRLGMRESCreate(MPI_Comm comm, HYPRE_Solver *solver)
1263 { return HYPRE_ParCSRLGMRESCreate(comm, solver);}
1264
1265//==============================================================================
1266int Ifpack_Hypre::Hypre_ParCSRBiCGSTABCreate(MPI_Comm comm, HYPRE_Solver *solver)
1267 { return HYPRE_ParCSRBiCGSTABCreate(comm, solver);}
1268
1269//==============================================================================
1270Teuchos::RCP<const Epetra_Map> Ifpack_Hypre::MakeContiguousColumnMap(Teuchos::RCP<const Epetra_RowMatrix> &MatrixRow) const{
1271 // Must create GloballyContiguous DomainMap (which is a permutation of Matrix_'s
1272 // DomainMap) and the corresponding permuted ColumnMap.
1273 // Epetra_GID ---------> LID ----------> HYPRE_GID
1274 // via DomainMap.LID() via GloballyContiguousDomainMap.GID()
1275 Teuchos::RCP<const Epetra_CrsMatrix> Matrix = Teuchos::rcp_dynamic_cast<const Epetra_CrsMatrix>(MatrixRow);
1276 if(Matrix.is_null())
1277 throw std::runtime_error("Ifpack_Hypre: Unsupported matrix configuration: Epetra_CrsMatrix required");
1278 const Epetra_Map & DomainMap = Matrix->DomainMap();
1279 const Epetra_Map & ColumnMap = Matrix->ColMap();
1280 const Epetra_Import * importer = Matrix->Importer();
1281
1282 if(DomainMap.LinearMap() ) {
1283 // If the domain map is linear, then we can just use the column map as is.
1284 return rcpFromRef(ColumnMap);
1285 }
1286 else {
1287 // The domain map isn't linear, so we need a new domain map
1288 Teuchos::RCP<Epetra_Map> ContiguousDomainMap = rcp(new Epetra_Map(DomainMap.NumGlobalElements(),
1289 DomainMap.NumMyElements(), 0, Comm()));
1290 if(importer) {
1291 // If there's an importer then we can use it to get a new column map
1292 Epetra_IntVector MyGIDsHYPRE(View,DomainMap,ContiguousDomainMap->MyGlobalElements());
1293
1294 // import the HYPRE GIDs
1295 Epetra_IntVector ColGIDsHYPRE(ColumnMap);
1296 ColGIDsHYPRE.Import(MyGIDsHYPRE, *importer, Insert);
1297
1298 // Make a HYPRE numbering-based column map.
1299 return Teuchos::rcp(new Epetra_Map(ColumnMap.NumGlobalElements(),ColGIDsHYPRE.MyLength(), &ColGIDsHYPRE[0], 0, Comm()));
1300 }
1301 else {
1302 // The problem has matching domain/column maps, and somehow the domain map isn't linear, so just use the new domain map
1303 return Teuchos::rcp(new Epetra_Map(ColumnMap.NumGlobalElements(),ColumnMap.NumMyElements(), ContiguousDomainMap->MyGlobalElements(), 0, Comm()));
1304 }
1305 }
1306}
1307
1308
1309
1310#endif // HAVE_HYPRE && HAVE_MPI
Ifpack_Utils.h
Epetra_BlockMap::LinearMap
bool LinearMap() const
Epetra_BlockMap::NumGlobalElements
int NumGlobalElements() const
Epetra_BlockMap::NumMyElements
int NumMyElements() const
Epetra_Import
Epetra_IntVector
Epetra_Map
Epetra_MultiVector
Epetra_MultiVector::NumVectors
int NumVectors() const
Epetra_MultiVector::MyLength
int MyLength() const
Epetra_MultiVector::Pointers
double ** Pointers() const
Epetra_RowMatrix
Epetra_RowMatrix::RowMatrixRowMap
virtual const Epetra_Map & RowMatrixRowMap() const=0
Generated by doxygen 1.9.8