Ifpack2 Templated Preconditioning Package Version 1.0
Loading...
Searching...
No Matches
Ifpack2_Hypre_def.hpp
1// @HEADER
2// *****************************************************************************
3// Ifpack2: Templated Object-Oriented Algebraic Preconditioner Package
4//
5// Copyright 2009 NTESS and the Ifpack2 contributors.
6// SPDX-License-Identifier: BSD-3-Clause
7// *****************************************************************************
8// @HEADER
9
10#ifndef IFPACK2_HYPRE_DEF_HPP
11#define IFPACK2_HYPRE_DEF_HPP
12
13#include "Ifpack2_Hypre_decl.hpp"
14#if defined(HAVE_IFPACK2_HYPRE) && defined(HAVE_IFPACK2_MPI)
15#include <stdexcept>
16
17#include "Tpetra_Import.hpp"
18#include "Teuchos_ParameterList.hpp"
19#include "Teuchos_RCP.hpp"
20#include "Teuchos_DefaultMpiComm.hpp"
21#include "HYPRE_IJ_mv.h"
22#include "HYPRE_parcsr_ls.h"
23#include "krylov.h"
24#include "_hypre_parcsr_mv.h"
25#include "_hypre_IJ_mv.h"
26#include "HYPRE_parcsr_mv.h"
27#include "HYPRE.h"
28
29
30using Teuchos::RCP;
31using Teuchos::rcp;
32using Teuchos::rcpFromRef;
33
34
35namespace Ifpack2 {
36
37template<class LocalOrdinal, class Node>
38Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_Int, Node> >::
39Hypre(const Teuchos::RCP<const row_matrix_type>& A):
40 A_(A),
41 IsInitialized_(false),
42 IsComputed_(false), NumInitialize_(0),
43 NumCompute_(0),
44 NumApply_(0),
45 InitializeTime_(0.0),
46 ComputeTime_(0.0),
47 ApplyTime_(0.0),
48 HypreA_(0),
49 HypreG_(0),
50 xHypre_(0),
51 yHypre_(0),
52 zHypre_(0),
53 IsSolverCreated_(false),
54 IsPrecondCreated_(false),
55 SolveOrPrec_(Hypre_Is_Solver),
56 NumFunsToCall_(0),
57 SolverType_(PCG),
58 PrecondType_(Euclid),
59 UsePreconditioner_(false),
60 Dump_(false) { }
61
62//==============================================================================
63template<class LocalOrdinal, class Node>
64Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_Int, Node> >::~Hypre() {
65 Destroy();
66}
67
68//==============================================================================
69template<class LocalOrdinal, class Node>
70void Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_Int, Node> >::Destroy(){
71 if(isInitialized()){
72 IFPACK2_CHK_ERRV(HYPRE_IJMatrixDestroy(HypreA_));
73 IFPACK2_CHK_ERRV(HYPRE_IJVectorDestroy(XHypre_));
74 IFPACK2_CHK_ERRV(HYPRE_IJVectorDestroy(YHypre_));
75 }
76 if(IsSolverCreated_){
77 IFPACK2_CHK_ERRV(SolverDestroyPtr_(Solver_));
78 }
79 if(IsPrecondCreated_){
80 IFPACK2_CHK_ERRV(PrecondDestroyPtr_(Preconditioner_));
81 }
82
83 // Maxwell
84 if(HypreG_) {
85 IFPACK2_CHK_ERRV(HYPRE_IJMatrixDestroy(HypreG_));
86 }
87 if(xHypre_) {
88 IFPACK2_CHK_ERRV(HYPRE_IJVectorDestroy(xHypre_));
89 }
90 if(yHypre_) {
91 IFPACK2_CHK_ERRV(HYPRE_IJVectorDestroy(yHypre_));
92 }
93 if(zHypre_) {
94 IFPACK2_CHK_ERRV(HYPRE_IJVectorDestroy(zHypre_));
95 }
96} //Destroy()
97
98//==============================================================================
99template<class LocalOrdinal, class Node>
100void Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_Int, Node> >::initialize(){
101 const std::string timerName ("Ifpack2::Hypre::initialize");
102 Teuchos::RCP<Teuchos::Time> timer = Teuchos::TimeMonitor::lookupCounter (timerName);
103 if (timer.is_null ()) timer = Teuchos::TimeMonitor::getNewCounter (timerName);
104
105 if(IsInitialized_) return;
106 double startTime = timer->wallTime();
107 {
108 Teuchos::TimeMonitor timeMon (*timer);
109
110 MPI_Comm comm = * (Teuchos::rcp_dynamic_cast<const Teuchos::MpiComm<int> >(A_->getRowMap()->getComm())->getRawMpiComm());
111
112 // Check that RowMap and RangeMap are the same. While this could handle the
113 // case where RowMap and RangeMap are permutations, other Ifpack PCs don't
114 // handle this either.
115 if (!A_->getRowMap()->isSameAs(*A_->getRangeMap())) {
116 IFPACK2_CHK_ERRV(-1);
117 }
118 // Hypre expects the RowMap to be Linear.
119 if (A_->getRowMap()->isContiguous()) {
120 GloballyContiguousRowMap_ = A_->getRowMap();
121 GloballyContiguousColMap_ = A_->getColMap();
122 } else {
123 // Must create GloballyContiguous Maps for Hypre
124 if(A_->getDomainMap()->isSameAs(*A_->getRowMap())) {
125 Teuchos::RCP<const crs_matrix_type> Aconst = Teuchos::rcp_dynamic_cast<const crs_matrix_type>(A_);
126 GloballyContiguousColMap_ = MakeContiguousColumnMap(Aconst);
127 GloballyContiguousRowMap_ = rcp(new map_type(A_->getRowMap()->getGlobalNumElements(),
128 A_->getRowMap()->getLocalNumElements(), 0, A_->getRowMap()->getComm()));
129 }
130 else {
131 throw std::runtime_error("Ifpack_Hypre: Unsupported map configuration: Row/Domain maps do not match");
132 }
133 }
134 // Next create vectors that will be used when ApplyInverse() is called
135 HYPRE_Int ilower = GloballyContiguousRowMap_->getMinGlobalIndex();
136 HYPRE_Int iupper = GloballyContiguousRowMap_->getMaxGlobalIndex();
137 // X in AX = Y
138 IFPACK2_CHK_ERRV(HYPRE_IJVectorCreate(comm, ilower, iupper, &XHypre_));
139 IFPACK2_CHK_ERRV(HYPRE_IJVectorSetObjectType(XHypre_, HYPRE_PARCSR));
140 IFPACK2_CHK_ERRV(HYPRE_IJVectorInitialize(XHypre_));
141 IFPACK2_CHK_ERRV(HYPRE_IJVectorAssemble(XHypre_));
142 IFPACK2_CHK_ERRV(HYPRE_IJVectorGetObject(XHypre_, (void**) &ParX_));
143 XVec_ = Teuchos::rcp((hypre_ParVector *) hypre_IJVectorObject(((hypre_IJVector *) XHypre_)),false);
144
145 // Y in AX = Y
146 IFPACK2_CHK_ERRV(HYPRE_IJVectorCreate(comm, ilower, iupper, &YHypre_));
147 IFPACK2_CHK_ERRV(HYPRE_IJVectorSetObjectType(YHypre_, HYPRE_PARCSR));
148 IFPACK2_CHK_ERRV(HYPRE_IJVectorInitialize(YHypre_));
149 IFPACK2_CHK_ERRV(HYPRE_IJVectorAssemble(YHypre_));
150 IFPACK2_CHK_ERRV(HYPRE_IJVectorGetObject(YHypre_, (void**) &ParY_));
151 YVec_ = Teuchos::rcp((hypre_ParVector *) hypre_IJVectorObject(((hypre_IJVector *) YHypre_)),false);
152
153 // Cache
154 VectorCache_.resize(A_->getRowMap()->getLocalNumElements());
155
156 // set flags
157 IsInitialized_=true;
158 NumInitialize_++;
159 }
160 InitializeTime_ += (timer->wallTime() - startTime);
161} //Initialize()
162
163//==============================================================================
164template<class LocalOrdinal, class Node>
165void Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_Int, Node> >::setParameters(const Teuchos::ParameterList& list){
166
167 std::map<std::string, Hypre_Solver> solverMap;
168 solverMap["BoomerAMG"] = BoomerAMG;
169 solverMap["ParaSails"] = ParaSails;
170 solverMap["Euclid"] = Euclid;
171 solverMap["AMS"] = AMS;
172 solverMap["Hybrid"] = Hybrid;
173 solverMap["PCG"] = PCG;
174 solverMap["GMRES"] = GMRES;
175 solverMap["FlexGMRES"] = FlexGMRES;
176 solverMap["LGMRES"] = LGMRES;
177 solverMap["BiCGSTAB"] = BiCGSTAB;
178
179 std::map<std::string, Hypre_Chooser> chooserMap;
180 chooserMap["Solver"] = Hypre_Is_Solver;
181 chooserMap["Preconditioner"] = Hypre_Is_Preconditioner;
182
183 List_ = list;
184 Hypre_Solver solType;
185 if (list.isType<std::string>("hypre: Solver"))
186 solType = solverMap[list.get<std::string>("hypre: Solver")];
187 else if(list.isParameter("hypre: Solver"))
188 solType = (Hypre_Solver) list.get<int>("hypre: Solver");
189 else
190 solType = PCG;
191 SolverType_ = solType;
192 Hypre_Solver precType;
193 if (list.isType<std::string>("hypre: Preconditioner"))
194 precType = solverMap[list.get<std::string>("hypre: Preconditioner")];
195 else if(list.isParameter("hypre: Preconditioner"))
196 precType = (Hypre_Solver) list.get<int>("hypre: Preconditioner");
197 else
198 precType = Euclid;
199 PrecondType_ = precType;
200 Hypre_Chooser chooser;
201 if (list.isType<std::string>("hypre: SolveOrPrecondition"))
202 chooser = chooserMap[list.get<std::string>("hypre: SolveOrPrecondition")];
203 else if(list.isParameter("hypre: SolveOrPrecondition"))
204 chooser = (Hypre_Chooser) list.get<int>("hypre: SolveOrPrecondition");
205 else
206 chooser = Hypre_Is_Solver;
207 SolveOrPrec_ = chooser;
208 bool SetPrecond = list.isParameter("hypre: SetPreconditioner") ? list.get<bool>("hypre: SetPreconditioner") : false;
209 IFPACK2_CHK_ERR(SetParameter(SetPrecond));
210 int NumFunctions = list.isParameter("hypre: NumFunctions") ? list.get<int>("hypre: NumFunctions") : 0;
211 FunsToCall_.clear();
212 NumFunsToCall_ = 0;
213 if(NumFunctions > 0){
214 RCP<FunctionParameter>* params = list.get<RCP<FunctionParameter>*>("hypre: Functions");
215 for(int i = 0; i < NumFunctions; i++){
216 IFPACK2_CHK_ERR(AddFunToList(params[i]));
217 }
218 }
219
220 if (list.isSublist("hypre: Solver functions")) {
221 Teuchos::ParameterList solverList = list.sublist("hypre: Solver functions");
222 for (auto it = solverList.begin(); it != solverList.end(); ++it) {
223 std::string funct_name = it->first;
224 if (it->second.isType<HYPRE_Int>()) {
225 IFPACK2_CHK_ERR(AddFunToList(rcp(new FunctionParameter(Hypre_Is_Solver, funct_name , Teuchos::getValue<HYPRE_Int>(it->second)))));
226 } else if (!std::is_same<HYPRE_Int,int>::value && it->second.isType<int>()) {
227 IFPACK2_CHK_ERR(AddFunToList(rcp(new FunctionParameter(Hypre_Is_Solver, funct_name , Teuchos::as<HYPRE_Int>(Teuchos::getValue<int>(it->second))))));
228 } else if (it->second.isType<HYPRE_Real>()) {
229 IFPACK2_CHK_ERR(AddFunToList(rcp(new FunctionParameter(Hypre_Is_Solver, funct_name , Teuchos::getValue<HYPRE_Real>(it->second)))));
230 } else {
231 IFPACK2_CHK_ERR(-1);
232 }
233 }
234 }
235
236 if (list.isSublist("hypre: Preconditioner functions")) {
237 Teuchos::ParameterList precList = list.sublist("hypre: Preconditioner functions");
238 for (auto it = precList.begin(); it != precList.end(); ++it) {
239 std::string funct_name = it->first;
240 if (it->second.isType<HYPRE_Int>()) {
241 IFPACK2_CHK_ERR(AddFunToList(rcp(new FunctionParameter(Hypre_Is_Preconditioner, funct_name , Teuchos::getValue<HYPRE_Int>(it->second)))));
242 } else if (!std::is_same<HYPRE_Int,int>::value && it->second.isType<int>()) {
243 IFPACK2_CHK_ERR(AddFunToList(rcp(new FunctionParameter(Hypre_Is_Preconditioner, funct_name , Teuchos::as<HYPRE_Int>(Teuchos::getValue<int>(it->second))))));
244 } else if (it->second.isType<HYPRE_Real>()) {
245 IFPACK2_CHK_ERR(AddFunToList(rcp(new FunctionParameter(Hypre_Is_Preconditioner, funct_name , Teuchos::getValue<HYPRE_Real>(it->second)))));
246 } else if (it->second.isList()) {
247 Teuchos::ParameterList pl = Teuchos::getValue<Teuchos::ParameterList>(it->second);
248 if (FunctionParameter::isFuncIntInt(funct_name)) {
249 HYPRE_Int arg0 = pl.get<HYPRE_Int>("arg 0");
250 HYPRE_Int arg1 = pl.get<HYPRE_Int>("arg 1");
251 IFPACK2_CHK_ERR(AddFunToList(rcp(new FunctionParameter(Hypre_Is_Preconditioner, funct_name , arg0, arg1))));
252 } else if (FunctionParameter::isFuncIntIntDoubleDouble(funct_name)) {
253 HYPRE_Int arg0 = pl.get<HYPRE_Int>("arg 0");
254 HYPRE_Int arg1 = pl.get<HYPRE_Int>("arg 1");
255 HYPRE_Real arg2 = pl.get<HYPRE_Real>("arg 2");
256 HYPRE_Real arg3 = pl.get<HYPRE_Real>("arg 3");
257 IFPACK2_CHK_ERR(AddFunToList(rcp(new FunctionParameter(Hypre_Is_Preconditioner, funct_name , arg0, arg1, arg2, arg3))));
258 } else if (FunctionParameter::isFuncIntIntIntDoubleIntInt(funct_name)) {
259 HYPRE_Int arg0 = pl.get<HYPRE_Int>("arg 0");
260 HYPRE_Int arg1 = pl.get<HYPRE_Int>("arg 1");
261 HYPRE_Int arg2 = pl.get<HYPRE_Int>("arg 2");
262 HYPRE_Real arg3 = pl.get<HYPRE_Real>("arg 3");
263 HYPRE_Int arg4 = pl.get<HYPRE_Int>("arg 4");
264 HYPRE_Int arg5 = pl.get<HYPRE_Int>("arg 5");
265 IFPACK2_CHK_ERR(AddFunToList(rcp(new FunctionParameter(Hypre_Is_Preconditioner, funct_name , arg0, arg1, arg2, arg3, arg4, arg5))));
266 } else {
267 IFPACK2_CHK_ERR(-1);
268 }
269 }
270 }
271 }
272
273 if (list.isSublist("Coordinates") && list.sublist("Coordinates").isType<Teuchos::RCP<multivector_type> >("Coordinates"))
274 Coords_ = list.sublist("Coordinates").get<Teuchos::RCP<multivector_type> >("Coordinates");
275 if (list.isSublist("Operators") && list.sublist("Operators").isType<Teuchos::RCP<const crs_matrix_type> >("G"))
276 G_ = list.sublist("Operators").get<Teuchos::RCP<const crs_matrix_type> >("G");
277
278 Dump_ = list.isParameter("hypre: Dump") ? list.get<bool>("hypre: Dump") : false;
279} //setParameters()
280
281//==============================================================================
282template<class LocalOrdinal, class Node>
283int Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_Int, Node> >::AddFunToList(RCP<FunctionParameter> NewFun){
284 NumFunsToCall_ = NumFunsToCall_+1;
285 FunsToCall_.resize(NumFunsToCall_);
286 FunsToCall_[NumFunsToCall_-1] = NewFun;
287 return 0;
288} //AddFunToList()
289
290//==============================================================================
291template<class LocalOrdinal, class Node>
292int Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_Int, Node> >::SetParameter(Hypre_Chooser chooser, HYPRE_Int (*pt2Func)(HYPRE_Solver, HYPRE_Int), HYPRE_Int parameter){
293 RCP<FunctionParameter> temp = rcp(new FunctionParameter(chooser, pt2Func, parameter));
294 IFPACK2_CHK_ERR(AddFunToList(temp));
295 return 0;
296} //SetParameter() - int function pointer
297
298//=============================================================================
299template<class LocalOrdinal, class Node>
300int Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_Int, Node> >::SetParameter(Hypre_Chooser chooser, HYPRE_Int (*pt2Func)(HYPRE_Solver, HYPRE_Real), HYPRE_Real parameter){
301 RCP<FunctionParameter> temp = rcp(new FunctionParameter(chooser, pt2Func, parameter));
302 IFPACK2_CHK_ERR(AddFunToList(temp));
303 return 0;
304} //SetParameter() - double function pointer
305
306//==============================================================================
307template<class LocalOrdinal, class Node>
308int Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_Int, Node> >::SetParameter(Hypre_Chooser chooser, HYPRE_Int (*pt2Func)(HYPRE_Solver, HYPRE_Real, HYPRE_Int), HYPRE_Real parameter1, HYPRE_Int parameter2){
309 RCP<FunctionParameter> temp = rcp(new FunctionParameter(chooser, pt2Func, parameter1, parameter2));
310 IFPACK2_CHK_ERR(AddFunToList(temp));
311 return 0;
312} //SetParameter() - double,int function pointer
313
314//==============================================================================
315template<class LocalOrdinal, class Node>
316int Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_Int, Node> >::SetParameter(Hypre_Chooser chooser, HYPRE_Int (*pt2Func)(HYPRE_Solver, HYPRE_Int, HYPRE_Real), HYPRE_Int parameter1, HYPRE_Real parameter2){
317 RCP<FunctionParameter> temp = rcp(new FunctionParameter(chooser, pt2Func, parameter1, parameter2));
318 IFPACK2_CHK_ERR(AddFunToList(temp));
319 return 0;
320} //SetParameter() - int,double function pointer
321
322//==============================================================================
323template<class LocalOrdinal, class Node>
324int Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_Int, Node> >::SetParameter(Hypre_Chooser chooser, HYPRE_Int (*pt2Func)(HYPRE_Solver, HYPRE_Int, HYPRE_Int), HYPRE_Int parameter1, HYPRE_Int parameter2){
325 RCP<FunctionParameter> temp = rcp(new FunctionParameter(chooser, pt2Func, parameter1, parameter2));
326 IFPACK2_CHK_ERR(AddFunToList(temp));
327 return 0;
328} //SetParameter() int,int function pointer
329
330//==============================================================================
331template<class LocalOrdinal, class Node>
332int Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_Int, Node> >::SetParameter(Hypre_Chooser chooser, HYPRE_Int (*pt2Func)(HYPRE_Solver, HYPRE_Real*), HYPRE_Real* parameter){
333 RCP<FunctionParameter> temp = rcp(new FunctionParameter(chooser, pt2Func, parameter));
334 IFPACK2_CHK_ERR(AddFunToList(temp));
335 return 0;
336} //SetParameter() - double* function pointer
337
338//==============================================================================
339template<class LocalOrdinal, class Node>
340int Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_Int, Node> >::SetParameter(Hypre_Chooser chooser, HYPRE_Int (*pt2Func)(HYPRE_Solver, HYPRE_Int*), HYPRE_Int* parameter){
341 RCP<FunctionParameter> temp = rcp(new FunctionParameter(chooser, pt2Func, parameter));
342 IFPACK2_CHK_ERR(AddFunToList(temp));
343 return 0;
344} //SetParameter() - int* function pointer
345
346//==============================================================================
347template<class LocalOrdinal, class Node>
348int Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_Int, Node> >::SetParameter(Hypre_Chooser chooser, HYPRE_Int (*pt2Func)(HYPRE_Solver, HYPRE_Int**), HYPRE_Int** parameter){
349 RCP<FunctionParameter> temp = rcp(new FunctionParameter(chooser, pt2Func, parameter));
350 IFPACK2_CHK_ERR(AddFunToList(temp));
351 return 0;
352} //SetParameter() - int** function pointer
353
354//==============================================================================
355template<class LocalOrdinal, class Node>
356int Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_Int, Node> >::SetParameter(Hypre_Chooser chooser, Hypre_Solver solver){
357 if(chooser == Hypre_Is_Solver){
358 SolverType_ = solver;
359 } else {
360 PrecondType_ = solver;
361 }
362 return 0;
363} //SetParameter() - set type of solver
364
365//==============================================================================
366template<class LocalOrdinal, class Node>
367int Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_Int, Node> >::SetDiscreteGradient(Teuchos::RCP<const crs_matrix_type> G){
368 using LO = local_ordinal_type;
369 using GO = global_ordinal_type;
370 // using SC = scalar_type;
371
372 // Sanity check
373 if(!A_->getRowMap()->isSameAs(*G->getRowMap()))
374 throw std::runtime_error("Hypre<Tpetra::RowMatrix<double, HYPRE_Int, long long, Node>: Edge map mismatch: A and discrete gradient");
375
376 // Get the maps for the nodes (assuming the edge map from A is OK);
377 GloballyContiguousNodeRowMap_ = rcp(new map_type(G->getDomainMap()->getGlobalNumElements(),
378 G->getDomainMap()->getLocalNumElements(), 0, A_->getRowMap()->getComm()));
379 GloballyContiguousNodeColMap_ = MakeContiguousColumnMap(G);
380
381 // Start building G
382 MPI_Comm comm = * (Teuchos::rcp_dynamic_cast<const Teuchos::MpiComm<int> >(A_->getRowMap()->getComm())->getRawMpiComm());
383 GO ilower = GloballyContiguousRowMap_->getMinGlobalIndex();
384 GO iupper = GloballyContiguousRowMap_->getMaxGlobalIndex();
385 GO jlower = GloballyContiguousNodeRowMap_->getMinGlobalIndex();
386 GO jupper = GloballyContiguousNodeRowMap_->getMaxGlobalIndex();
387 IFPACK2_CHK_ERR(HYPRE_IJMatrixCreate(comm, ilower, iupper, jlower, jupper, &HypreG_));
388 IFPACK2_CHK_ERR(HYPRE_IJMatrixSetObjectType(HypreG_, HYPRE_PARCSR));
389 IFPACK2_CHK_ERR(HYPRE_IJMatrixInitialize(HypreG_));
390
391 std::vector<GO> new_indices(G->getLocalMaxNumRowEntries());
392 for(LO i = 0; i < (LO)G->getLocalNumRows(); i++){
393 typename crs_matrix_type::values_host_view_type values;
394 typename crs_matrix_type::local_inds_host_view_type indices;
395 G->getLocalRowView(i, indices, values);
396 for(LO j = 0; j < (LO) indices.extent(0); j++){
397 new_indices[j] = GloballyContiguousNodeColMap_->getGlobalElement(indices(j));
398 }
399 GO GlobalRow[1];
400 GO numEntries = (GO) indices.extent(0);
401 GlobalRow[0] = GloballyContiguousRowMap_->getGlobalElement(i);
402 IFPACK2_CHK_ERR(HYPRE_IJMatrixSetValues(HypreG_, 1, &numEntries, GlobalRow, new_indices.data(), values.data()));
403 }
404 IFPACK2_CHK_ERR(HYPRE_IJMatrixAssemble(HypreG_));
405 IFPACK2_CHK_ERR(HYPRE_IJMatrixGetObject(HypreG_, (void**)&ParMatrixG_));
406
407 if (Dump_)
408 HYPRE_ParCSRMatrixPrint(ParMatrixG_,"G.mat");
409
410 if(SolverType_ == AMS)
411 HYPRE_AMSSetDiscreteGradient(Solver_, ParMatrixG_);
412 if(PrecondType_ == AMS)
413 HYPRE_AMSSetDiscreteGradient(Preconditioner_, ParMatrixG_);
414 return 0;
415} //SetDiscreteGradient()
416
417//==============================================================================
418template<class LocalOrdinal, class Node>
419int Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_Int, Node> >::SetCoordinates(Teuchos::RCP<multivector_type> coords) {
420
421 if(!G_.is_null() && !G_->getDomainMap()->isSameAs(*coords->getMap()))
422 throw std::runtime_error("Hypre<Tpetra::RowMatrix<double, HYPRE_Int, long long, Node>: Node map mismatch: G->DomainMap() and coords");
423
424 if(SolverType_ != AMS && PrecondType_ != AMS)
425 return 0;
426
427 scalar_type *xPtr = coords->getDataNonConst(0).getRawPtr();
428 scalar_type *yPtr = coords->getDataNonConst(1).getRawPtr();
429 scalar_type *zPtr = coords->getDataNonConst(2).getRawPtr();
430
431 MPI_Comm comm = * (Teuchos::rcp_dynamic_cast<const Teuchos::MpiComm<int> >(A_->getRowMap()->getComm())->getRawMpiComm());
432 local_ordinal_type NumEntries = coords->getLocalLength();
433 global_ordinal_type * indices = const_cast<global_ordinal_type*>(GloballyContiguousNodeRowMap_->getLocalElementList().getRawPtr());
434
435 global_ordinal_type ilower = GloballyContiguousNodeRowMap_->getMinGlobalIndex();
436 global_ordinal_type iupper = GloballyContiguousNodeRowMap_->getMaxGlobalIndex();
437
438 if( NumEntries != iupper-ilower+1) {
439 std::cout<<"Ifpack2::Hypre::SetCoordinates(): Error on rank "<<A_->getRowMap()->getComm()->getRank()<<": MyLength = "<<coords->getLocalLength()<<" GID range = ["<<ilower<<","<<iupper<<"]"<<std::endl;
440 throw std::runtime_error("Hypre<Tpetra::RowMatrix<double, HYPRE_Int, long long, Node>: SetCoordinates: Length mismatch");
441 }
442
443 IFPACK2_CHK_ERR(HYPRE_IJVectorCreate(comm, ilower, iupper, &xHypre_));
444 IFPACK2_CHK_ERR(HYPRE_IJVectorSetObjectType(xHypre_, HYPRE_PARCSR));
445 IFPACK2_CHK_ERR(HYPRE_IJVectorInitialize(xHypre_));
446 IFPACK2_CHK_ERR(HYPRE_IJVectorSetValues(xHypre_,NumEntries,indices,xPtr));
447 IFPACK2_CHK_ERR(HYPRE_IJVectorAssemble(xHypre_));
448 IFPACK2_CHK_ERR(HYPRE_IJVectorGetObject(xHypre_, (void**) &xPar_));
449
450 IFPACK2_CHK_ERR(HYPRE_IJVectorCreate(comm, ilower, iupper, &yHypre_));
451 IFPACK2_CHK_ERR(HYPRE_IJVectorSetObjectType(yHypre_, HYPRE_PARCSR));
452 IFPACK2_CHK_ERR(HYPRE_IJVectorInitialize(yHypre_));
453 IFPACK2_CHK_ERR(HYPRE_IJVectorSetValues(yHypre_,NumEntries,indices,yPtr));
454 IFPACK2_CHK_ERR(HYPRE_IJVectorAssemble(yHypre_));
455 IFPACK2_CHK_ERR(HYPRE_IJVectorGetObject(yHypre_, (void**) &yPar_));
456
457 IFPACK2_CHK_ERR(HYPRE_IJVectorCreate(comm, ilower, iupper, &zHypre_));
458 IFPACK2_CHK_ERR(HYPRE_IJVectorSetObjectType(zHypre_, HYPRE_PARCSR));
459 IFPACK2_CHK_ERR(HYPRE_IJVectorInitialize(zHypre_));
460 IFPACK2_CHK_ERR(HYPRE_IJVectorSetValues(zHypre_,NumEntries,indices,zPtr));
461 IFPACK2_CHK_ERR(HYPRE_IJVectorAssemble(zHypre_));
462 IFPACK2_CHK_ERR(HYPRE_IJVectorGetObject(zHypre_, (void**) &zPar_));
463
464 if (Dump_) {
465 HYPRE_ParVectorPrint(xPar_,"coordX.dat");
466 HYPRE_ParVectorPrint(yPar_,"coordY.dat");
467 HYPRE_ParVectorPrint(zPar_,"coordZ.dat");
468 }
469
470 if(SolverType_ == AMS)
471 HYPRE_AMSSetCoordinateVectors(Solver_, xPar_, yPar_, zPar_);
472 if(PrecondType_ == AMS)
473 HYPRE_AMSSetCoordinateVectors(Preconditioner_, xPar_, yPar_, zPar_);
474
475 return 0;
476
477} //SetCoordinates
478
479//==============================================================================
480template<class LocalOrdinal, class Node>
481void Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_Int, Node> >::compute(){
482 const std::string timerName ("Ifpack2::Hypre::compute");
483 Teuchos::RCP<Teuchos::Time> timer = Teuchos::TimeMonitor::lookupCounter (timerName);
484 if (timer.is_null ()) timer = Teuchos::TimeMonitor::getNewCounter (timerName);
485 double startTime = timer->wallTime();
486 // Start timing here.
487 {
488 Teuchos::TimeMonitor timeMon (*timer);
489
490 if(isInitialized() == false){
491 initialize();
492 }
493
494 // Create the Hypre matrix and copy values. Note this uses values (which
495 // Initialize() shouldn't do) but it doesn't care what they are (for
496 // instance they can be uninitialized data even). It should be possible to
497 // set the Hypre structure without copying values, but this is the easiest
498 // way to get the structure.
499 MPI_Comm comm = * (Teuchos::rcp_dynamic_cast<const Teuchos::MpiComm<int> >(A_->getRowMap()->getComm())->getRawMpiComm());
500 global_ordinal_type ilower = GloballyContiguousRowMap_->getMinGlobalIndex();
501 global_ordinal_type iupper = GloballyContiguousRowMap_->getMaxGlobalIndex();
502 IFPACK2_CHK_ERR(HYPRE_IJMatrixCreate(comm, ilower, iupper, ilower, iupper, &HypreA_));
503 IFPACK2_CHK_ERR(HYPRE_IJMatrixSetObjectType(HypreA_, HYPRE_PARCSR));
504 IFPACK2_CHK_ERR(HYPRE_IJMatrixInitialize(HypreA_));
505 CopyTpetraToHypre();
506 if(SolveOrPrec_ == Hypre_Is_Solver) {
507 IFPACK2_CHK_ERR(SetSolverType(SolverType_));
508 if (SolverPrecondPtr_ != NULL && UsePreconditioner_) {
509 // both method allows a PC (first condition) and the user wants a PC (second)
510 IFPACK2_CHK_ERR(SetPrecondType(PrecondType_));
511 CallFunctions();
512 IFPACK2_CHK_ERR(SolverPrecondPtr_(Solver_, PrecondSolvePtr_, PrecondSetupPtr_, Preconditioner_));
513 } else {
514 CallFunctions();
515 }
516 } else {
517 IFPACK2_CHK_ERR(SetPrecondType(PrecondType_));
518 CallFunctions();
519 }
520
521 if (!G_.is_null()) {
522 SetDiscreteGradient(G_);
523 }
524
525 if (!Coords_.is_null()) {
526 SetCoordinates(Coords_);
527 }
528
529 // Hypre Setup must be called after matrix has values
530 if(SolveOrPrec_ == Hypre_Is_Solver){
531 IFPACK2_CHK_ERR(SolverSetupPtr_(Solver_, ParMatrix_, ParX_, ParY_));
532 } else {
533 IFPACK2_CHK_ERR(PrecondSetupPtr_(Preconditioner_, ParMatrix_, ParX_, ParY_));
534 }
535
536 IsComputed_ = true;
537 NumCompute_++;
538 }
539
540 ComputeTime_ += (timer->wallTime() - startTime);
541} //Compute()
542
543//==============================================================================
544template<class LocalOrdinal, class Node>
545int Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_Int, Node> >::CallFunctions() const{
546 for(int i = 0; i < NumFunsToCall_; i++){
547 IFPACK2_CHK_ERR(FunsToCall_[i]->CallFunction(Solver_, Preconditioner_));
548 }
549 return 0;
550} //CallFunctions()
551
552//==============================================================================
553template<class LocalOrdinal, class Node>
554void Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_Int, Node> >::apply (const Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& X,
555 Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& Y,
556 Teuchos::ETransp mode,
557 scalar_type alpha,
558 scalar_type beta) const {
559 using LO = local_ordinal_type;
560 using SC = scalar_type;
561 const std::string timerName ("Ifpack2::Hypre::apply");
562 Teuchos::RCP<Teuchos::Time> timer = Teuchos::TimeMonitor::lookupCounter (timerName);
563 if (timer.is_null ()) timer = Teuchos::TimeMonitor::getNewCounter (timerName);
564 double startTime = timer->wallTime();
565 // Start timing here.
566 {
567 Teuchos::TimeMonitor timeMon (*timer);
568
569 if(isComputed() == false){
570 IFPACK2_CHK_ERR(-1);
571 }
572 hypre_Vector *XLocal_ = hypre_ParVectorLocalVector(XVec_);
573 hypre_Vector *YLocal_ = hypre_ParVectorLocalVector(YVec_);
574 bool SameVectors = false;
575 size_t NumVectors = X.getNumVectors();
576 if (NumVectors != Y.getNumVectors()) IFPACK2_CHK_ERR(-1); // X and Y must have same number of vectors
577 if(&X == &Y) { //FIXME: Maybe not the right way to check this
578 SameVectors = true;
579 }
580
581 // NOTE: Here were assuming that the local ordering of Epetra's X/Y-vectors and
582 // Hypre's X/Y-vectors are the same. Seeing as as this is more or less how we constructed
583 // the Hypre matrices, this seems pretty reasoanble.
584
585 for(int VecNum = 0; VecNum < (int) NumVectors; VecNum++) {
586 //Get values for current vector in multivector.
587 // FIXME amk Nov 23, 2015: This will not work for funky data layouts
588 SC * XValues = const_cast<SC*>(X.getData(VecNum).getRawPtr());
589 SC * YValues;
590 if(!SameVectors){
591 YValues = const_cast<SC*>(Y.getData(VecNum).getRawPtr());
592 } else {
593 YValues = VectorCache_.getRawPtr();
594 }
595 // Temporarily make a pointer to data in Hypre for end
596 SC *XTemp = XLocal_->data;
597 // Replace data in Hypre vectors with Epetra data
598 XLocal_->data = XValues;
599 SC *YTemp = YLocal_->data;
600 YLocal_->data = YValues;
601
602 IFPACK2_CHK_ERR(HYPRE_ParVectorSetConstantValues(ParY_, 0.0));
603 if(SolveOrPrec_ == Hypre_Is_Solver){
604 // Use the solver methods
605 IFPACK2_CHK_ERR(SolverSolvePtr_(Solver_, ParMatrix_, ParX_, ParY_));
606 } else {
607 // Apply the preconditioner
608 IFPACK2_CHK_ERR(PrecondSolvePtr_(Preconditioner_, ParMatrix_, ParX_, ParY_));
609 }
610
611 if(SameVectors){
612 Teuchos::ArrayView<SC> Yv = Y.getDataNonConst(VecNum)();
613 LO NumEntries = Y.getLocalLength();
614 for(LO i = 0; i < NumEntries; i++)
615 Yv[i] = YValues[i];
616 }
617 XLocal_->data = XTemp;
618 YLocal_->data = YTemp;
619 }
620 NumApply_++;
621 }
622 ApplyTime_ += (timer->wallTime() - startTime);
623} //apply()
624
625
626//==============================================================================
627template<class LocalOrdinal, class Node>
628void Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_Int, Node> >::applyMat (const Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& X,
629 Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& Y,
630 Teuchos::ETransp mode) const {
631 A_->apply(X,Y,mode);
632} //applyMat()
633
634//==============================================================================
635template<class LocalOrdinal, class Node>
636std::string Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_Int, Node> >::description() const {
637 std::ostringstream out;
638
639 // Output is a valid YAML dictionary in flow style. If you don't
640 // like everything on a single line, you should call describe()
641 // instead.
642 out << "\"Ifpack2::Hypre\": {";
643 out << "Initialized: " << (isInitialized () ? "true" : "false") << ", "
644 << "Computed: " << (isComputed () ? "true" : "false") << ", ";
645
646 if (A_.is_null ()) {
647 out << "Matrix: null";
648 }
649 else {
650 out << "Global matrix dimensions: ["
651 << A_->getGlobalNumRows () << ", "
652 << A_->getGlobalNumCols () << "]"
653 << ", Global nnz: " << A_->getGlobalNumEntries();
654 }
655
656 out << "}";
657 return out.str ();
658} //description()
659
660//==============================================================================
661template<class LocalOrdinal, class Node>
662void Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_Int, Node> >::describe(Teuchos::FancyOStream &os, const Teuchos::EVerbosityLevel verbLevel) const {
663 using std::endl;
664 os << endl;
665 os << "================================================================================" << endl;
666 os << "Ifpack2::Hypre: " << endl << endl;
667 os << "Using " << A_->getComm()->getSize() << " processors." << endl;
668 os << "Global number of rows = " << A_->getGlobalNumRows() << endl;
669 os << "Global number of nonzeros = " << A_->getGlobalNumEntries() << endl;
670 // os << "Condition number estimate = " << Condest() << endl;
671 os << endl;
672 os << "Phase # calls Total Time (s)"<<endl;
673 os << "----- ------- --------------"<<endl;
674 os << "Initialize() " << std::setw(5) << NumInitialize_
675 << " " << std::setw(15) << InitializeTime_<<endl;
676 os << "Compute() " << std::setw(5) << NumCompute_
677 << " " << std::setw(15) << ComputeTime_ << endl;
678 os << "ApplyInverse() " << std::setw(5) << NumApply_
679 << " " << std::setw(15) << ApplyTime_ <<endl;
680 os << "================================================================================" << endl;
681 os << endl;
682} //description
683
684//==============================================================================
685template<class LocalOrdinal, class Node>
686int Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_Int, Node> >::SetSolverType(Hypre_Solver Solver){
687 switch(Solver) {
688 case BoomerAMG:
689 if(IsSolverCreated_){
690 SolverDestroyPtr_(Solver_);
691 IsSolverCreated_ = false;
692 }
693 SolverCreatePtr_ = &Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_Int, Node> >::Hypre_BoomerAMGCreate;
694 SolverDestroyPtr_ = &HYPRE_BoomerAMGDestroy;
695 SolverSetupPtr_ = &HYPRE_BoomerAMGSetup;
696 SolverPrecondPtr_ = NULL;
697 SolverSolvePtr_ = &HYPRE_BoomerAMGSolve;
698 break;
699 case AMS:
700 if(IsSolverCreated_){
701 SolverDestroyPtr_(Solver_);
702 IsSolverCreated_ = false;
703 }
704 SolverCreatePtr_ = &Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_Int, Node> >::Hypre_AMSCreate;
705 SolverDestroyPtr_ = &HYPRE_AMSDestroy;
706 SolverSetupPtr_ = &HYPRE_AMSSetup;
707 SolverSolvePtr_ = &HYPRE_AMSSolve;
708 SolverPrecondPtr_ = NULL;
709 break;
710 case Hybrid:
711 if(IsSolverCreated_){
712 SolverDestroyPtr_(Solver_);
713 IsSolverCreated_ = false;
714 }
715 SolverCreatePtr_ = &Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_Int, Node> >::Hypre_ParCSRHybridCreate;
716 SolverDestroyPtr_ = &HYPRE_ParCSRHybridDestroy;
717 SolverSetupPtr_ = &HYPRE_ParCSRHybridSetup;
718 SolverSolvePtr_ = &HYPRE_ParCSRHybridSolve;
719 SolverPrecondPtr_ = &HYPRE_ParCSRHybridSetPrecond;
720 break;
721 case PCG:
722 if(IsSolverCreated_){
723 SolverDestroyPtr_(Solver_);
724 IsSolverCreated_ = false;
725 }
726 SolverCreatePtr_ = &Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_Int, Node> >::Hypre_ParCSRPCGCreate;
727 SolverDestroyPtr_ = &HYPRE_ParCSRPCGDestroy;
728 SolverSetupPtr_ = &HYPRE_ParCSRPCGSetup;
729 SolverSolvePtr_ = &HYPRE_ParCSRPCGSolve;
730 SolverPrecondPtr_ = &HYPRE_ParCSRPCGSetPrecond;
731 break;
732 case GMRES:
733 if(IsSolverCreated_){
734 SolverDestroyPtr_(Solver_);
735 IsSolverCreated_ = false;
736 }
737 SolverCreatePtr_ = &Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_Int, Node> >::Hypre_ParCSRGMRESCreate;
738 SolverDestroyPtr_ = &HYPRE_ParCSRGMRESDestroy;
739 SolverSetupPtr_ = &HYPRE_ParCSRGMRESSetup;
740 SolverPrecondPtr_ = &HYPRE_ParCSRGMRESSetPrecond;
741 break;
742 case FlexGMRES:
743 if(IsSolverCreated_){
744 SolverDestroyPtr_(Solver_);
745 IsSolverCreated_ = false;
746 }
747 SolverCreatePtr_ = &Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_Int, Node> >::Hypre_ParCSRFlexGMRESCreate;
748 SolverDestroyPtr_ = &HYPRE_ParCSRFlexGMRESDestroy;
749 SolverSetupPtr_ = &HYPRE_ParCSRFlexGMRESSetup;
750 SolverSolvePtr_ = &HYPRE_ParCSRFlexGMRESSolve;
751 SolverPrecondPtr_ = &HYPRE_ParCSRFlexGMRESSetPrecond;
752 break;
753 case LGMRES:
754 if(IsSolverCreated_){
755 SolverDestroyPtr_(Solver_);
756 IsSolverCreated_ = false;
757 }
758 SolverCreatePtr_ = &Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_Int, Node> >::Hypre_ParCSRLGMRESCreate;
759 SolverDestroyPtr_ = &HYPRE_ParCSRLGMRESDestroy;
760 SolverSetupPtr_ = &HYPRE_ParCSRLGMRESSetup;
761 SolverSolvePtr_ = &HYPRE_ParCSRLGMRESSolve;
762 SolverPrecondPtr_ = &HYPRE_ParCSRLGMRESSetPrecond;
763 break;
764 case BiCGSTAB:
765 if(IsSolverCreated_){
766 SolverDestroyPtr_(Solver_);
767 IsSolverCreated_ = false;
768 }
769 SolverCreatePtr_ = &Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_Int, Node> >::Hypre_ParCSRBiCGSTABCreate;
770 SolverDestroyPtr_ = &HYPRE_ParCSRBiCGSTABDestroy;
771 SolverSetupPtr_ = &HYPRE_ParCSRBiCGSTABSetup;
772 SolverSolvePtr_ = &HYPRE_ParCSRBiCGSTABSolve;
773 SolverPrecondPtr_ = &HYPRE_ParCSRBiCGSTABSetPrecond;
774 break;
775 default:
776 return -1;
777 }
778 CreateSolver();
779 return 0;
780} //SetSolverType()
781
782//==============================================================================
783template<class LocalOrdinal, class Node>
784int Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_Int, Node> >::SetPrecondType(Hypre_Solver Precond){
785 switch(Precond) {
786 case BoomerAMG:
787 if(IsPrecondCreated_){
788 PrecondDestroyPtr_(Preconditioner_);
789 IsPrecondCreated_ = false;
790 }
791 PrecondCreatePtr_ = &Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_Int, Node> >::Hypre_BoomerAMGCreate;
792 PrecondDestroyPtr_ = &HYPRE_BoomerAMGDestroy;
793 PrecondSetupPtr_ = &HYPRE_BoomerAMGSetup;
794 PrecondSolvePtr_ = &HYPRE_BoomerAMGSolve;
795 break;
796 case ParaSails:
797 if(IsPrecondCreated_){
798 PrecondDestroyPtr_(Preconditioner_);
799 IsPrecondCreated_ = false;
800 }
801 PrecondCreatePtr_ = &Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_Int, Node> >::Hypre_ParaSailsCreate;
802 PrecondDestroyPtr_ = &HYPRE_ParaSailsDestroy;
803 PrecondSetupPtr_ = &HYPRE_ParaSailsSetup;
804 PrecondSolvePtr_ = &HYPRE_ParaSailsSolve;
805 break;
806 case Euclid:
807 if(IsPrecondCreated_){
808 PrecondDestroyPtr_(Preconditioner_);
809 IsPrecondCreated_ = false;
810 }
811 PrecondCreatePtr_ = &Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_Int, Node> >::Hypre_EuclidCreate;
812 PrecondDestroyPtr_ = &HYPRE_EuclidDestroy;
813 PrecondSetupPtr_ = &HYPRE_EuclidSetup;
814 PrecondSolvePtr_ = &HYPRE_EuclidSolve;
815 break;
816 case AMS:
817 if(IsPrecondCreated_){
818 PrecondDestroyPtr_(Preconditioner_);
819 IsPrecondCreated_ = false;
820 }
821 PrecondCreatePtr_ = &Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_Int, Node> >::Hypre_AMSCreate;
822 PrecondDestroyPtr_ = &HYPRE_AMSDestroy;
823 PrecondSetupPtr_ = &HYPRE_AMSSetup;
824 PrecondSolvePtr_ = &HYPRE_AMSSolve;
825 break;
826 default:
827 return -1;
828 }
829 CreatePrecond();
830 return 0;
831
832} //SetPrecondType()
833
834//==============================================================================
835template<class LocalOrdinal, class Node>
836int Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_Int, Node> >::CreateSolver(){
837 MPI_Comm comm;
838 HYPRE_ParCSRMatrixGetComm(ParMatrix_, &comm);
839 int ierr = (this->*SolverCreatePtr_)(comm, &Solver_);
840 IsSolverCreated_ = true;
841 return ierr;
842} //CreateSolver()
843
844//==============================================================================
845template<class LocalOrdinal, class Node>
846int Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_Int, Node> >::CreatePrecond(){
847 MPI_Comm comm;
848 HYPRE_ParCSRMatrixGetComm(ParMatrix_, &comm);
849 int ierr = (this->*PrecondCreatePtr_)(comm, &Preconditioner_);
850 IsPrecondCreated_ = true;
851 return ierr;
852} //CreatePrecond()
853
854
855//==============================================================================
856template<class LocalOrdinal, class Node>
857int Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_Int, Node> >::CopyTpetraToHypre(){
858 using LO = local_ordinal_type;
859 using GO = global_ordinal_type;
860 // using SC = scalar_type;
861
862 Teuchos::RCP<const crs_matrix_type> Matrix = Teuchos::rcp_dynamic_cast<const crs_matrix_type>(A_);
863 if(Matrix.is_null())
864 throw std::runtime_error("Hypre<Tpetra::RowMatrix<double, LocalOrdinal, HYPRE_Int, Node>: Unsupported matrix configuration: Tpetra::CrsMatrix required");
865
866 std::vector<HYPRE_Int> new_indices(Matrix->getLocalMaxNumRowEntries());
867 for(LO i = 0; i < (LO) Matrix->getLocalNumRows(); i++){
868 typename crs_matrix_type::values_host_view_type values;
869 typename crs_matrix_type::local_inds_host_view_type indices;
870 Matrix->getLocalRowView(i, indices, values);
871 for(LO j = 0; j < (LO)indices.extent(0); j++){
872 new_indices[j] = GloballyContiguousColMap_->getGlobalElement(indices(j));
873 }
874 HYPRE_Int GlobalRow[1];
875 HYPRE_Int numEntries = (GO) indices.extent(0);
876 GlobalRow[0] = GloballyContiguousRowMap_->getGlobalElement(i);
877 IFPACK2_CHK_ERR(HYPRE_IJMatrixSetValues(HypreA_, 1, &numEntries, GlobalRow, new_indices.data(), values.data()));
878 }
879 IFPACK2_CHK_ERR(HYPRE_IJMatrixAssemble(HypreA_));
880 IFPACK2_CHK_ERR(HYPRE_IJMatrixGetObject(HypreA_, (void**)&ParMatrix_));
881 if (Dump_)
882 HYPRE_ParCSRMatrixPrint(ParMatrix_,"A.mat");
883 return 0;
884} //CopyTpetraToHypre()
885
886//==============================================================================
887template<class LocalOrdinal, class Node>
888HYPRE_Int Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_Int, Node> >::Hypre_BoomerAMGCreate(MPI_Comm /*comm*/, HYPRE_Solver *solver)
889 { return HYPRE_BoomerAMGCreate(solver);}
890
891//==============================================================================
892template<class LocalOrdinal, class Node>
893HYPRE_Int Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_Int, Node> >::Hypre_ParaSailsCreate(MPI_Comm comm, HYPRE_Solver *solver)
894 { return HYPRE_ParaSailsCreate(comm, solver);}
895
896//==============================================================================
897template<class LocalOrdinal, class Node>
898HYPRE_Int Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_Int, Node> >::Hypre_EuclidCreate(MPI_Comm comm, HYPRE_Solver *solver)
899 { return HYPRE_EuclidCreate(comm, solver);}
900
901//==============================================================================
902template<class LocalOrdinal, class Node>
903HYPRE_Int Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_Int, Node> >::Hypre_AMSCreate(MPI_Comm /*comm*/, HYPRE_Solver *solver)
904 { return HYPRE_AMSCreate(solver);}
905
906//==============================================================================
907template<class LocalOrdinal, class Node>
908HYPRE_Int Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_Int, Node> >::Hypre_ParCSRHybridCreate(MPI_Comm /*comm*/, HYPRE_Solver *solver)
909 { return HYPRE_ParCSRHybridCreate(solver);}
910
911//==============================================================================
912template<class LocalOrdinal, class Node>
913HYPRE_Int Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_Int, Node> >::Hypre_ParCSRPCGCreate(MPI_Comm comm, HYPRE_Solver *solver)
914 { return HYPRE_ParCSRPCGCreate(comm, solver);}
915
916//==============================================================================
917template<class LocalOrdinal, class Node>
918HYPRE_Int Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_Int, Node> >::Hypre_ParCSRGMRESCreate(MPI_Comm comm, HYPRE_Solver *solver)
919 { return HYPRE_ParCSRGMRESCreate(comm, solver);}
920
921//==============================================================================
922template<class LocalOrdinal, class Node>
923HYPRE_Int Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_Int, Node> >::Hypre_ParCSRFlexGMRESCreate(MPI_Comm comm, HYPRE_Solver *solver)
924 { return HYPRE_ParCSRFlexGMRESCreate(comm, solver);}
925
926//==============================================================================
927template<class LocalOrdinal, class Node>
928HYPRE_Int Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_Int, Node> >::Hypre_ParCSRLGMRESCreate(MPI_Comm comm, HYPRE_Solver *solver)
929 { return HYPRE_ParCSRLGMRESCreate(comm, solver);}
930
931//==============================================================================
932template<class LocalOrdinal, class Node>
933HYPRE_Int Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_Int, Node> >::Hypre_ParCSRBiCGSTABCreate(MPI_Comm comm, HYPRE_Solver *solver)
934 { return HYPRE_ParCSRBiCGSTABCreate(comm, solver);}
935
936//==============================================================================
937template<class LocalOrdinal, class Node>
938Teuchos::RCP<const Tpetra::Map<LocalOrdinal, HYPRE_Int, Node> >
939Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_Int, Node> >::MakeContiguousColumnMap(Teuchos::RCP<const crs_matrix_type> &Matrix) const{
940 using import_type = Tpetra::Import<local_ordinal_type,global_ordinal_type,node_type>;
941 using go_vector_type = Tpetra::Vector<global_ordinal_type,local_ordinal_type,global_ordinal_type,node_type>;
942
943 // Must create GloballyContiguous DomainMap (which is a permutation of Matrix_'s
944 // DomainMap) and the corresponding permuted ColumnMap.
945 // Epetra_GID ---------> LID ----------> HYPRE_GID
946 // via DomainMap.LID() via GloballyContiguousDomainMap.GID()
947 if(Matrix.is_null())
948 throw std::runtime_error("Hypre<Tpetra::RowMatrix<HYPRE_Real, HYPRE_Int, long long, Node>: Unsupported matrix configuration: Tpetra::CrsMatrix required");
949 RCP<const map_type> DomainMap = Matrix->getDomainMap();
950 RCP<const map_type> ColumnMap = Matrix->getColMap();
951 RCP<const import_type> importer = Matrix->getGraph()->getImporter();
952
953 if(DomainMap->isContiguous() ) {
954 // If the domain map is linear, then we can just use the column map as is.
955 return ColumnMap;
956 }
957 else {
958 // The domain map isn't linear, so we need a new domain map
959 Teuchos::RCP<map_type> ContiguousDomainMap = rcp(new map_type(DomainMap->getGlobalNumElements(),
960 DomainMap->getLocalNumElements(), 0, DomainMap->getComm()));
961 if(importer) {
962 // If there's an importer then we can use it to get a new column map
963 go_vector_type MyGIDsHYPRE(DomainMap,ContiguousDomainMap->getLocalElementList());
964
965 // import the HYPRE GIDs
966 go_vector_type ColGIDsHYPRE(ColumnMap);
967 ColGIDsHYPRE.doImport(MyGIDsHYPRE, *importer, Tpetra::INSERT);
968
969 // Make a HYPRE numbering-based column map.
970 return Teuchos::rcp(new map_type(ColumnMap->getGlobalNumElements(),ColGIDsHYPRE.getDataNonConst()(),0, ColumnMap->getComm()));
971 }
972 else {
973 // The problem has matching domain/column maps, and somehow the domain map isn't linear, so just use the new domain map
974 return Teuchos::rcp(new map_type(ColumnMap->getGlobalNumElements(),ContiguousDomainMap->getLocalElementList(), 0, ColumnMap->getComm()));
975 }
976 }
977}
978
979
980
981template<class LocalOrdinal, class Node>
982int Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_Int, Node> >::getNumInitialize() const {
983 return NumInitialize_;
984}
985
986
987template<class LocalOrdinal, class Node>
988int Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_Int, Node> >::getNumCompute() const {
989 return NumCompute_;
990}
991
992
993template<class LocalOrdinal, class Node>
994int Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_Int, Node> >::getNumApply() const {
995 return NumApply_;
996}
997
998
999template<class LocalOrdinal, class Node>
1000double Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_Int, Node> >::getInitializeTime() const {
1001 return InitializeTime_;
1002}
1003
1004
1005template<class LocalOrdinal, class Node>
1006double Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_Int, Node> >::getComputeTime() const {
1007 return ComputeTime_;
1008}
1009
1010
1011template<class LocalOrdinal, class Node>
1012double Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_Int, Node> >::getApplyTime() const {
1013 return ApplyTime_;
1014}
1015
1016template<class LocalOrdinal, class Node>
1017Teuchos::RCP<const typename Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_Int, Node> >::map_type>
1018Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_Int, Node> >::
1019getDomainMap () const
1020{
1021 Teuchos::RCP<const row_matrix_type> A = getMatrix();
1022 TEUCHOS_TEST_FOR_EXCEPTION(
1023 A.is_null (), std::runtime_error, "Ifpack2::Hypre::getDomainMap: The "
1024 "input matrix A is null. Please call setMatrix() with a nonnull input "
1025 "matrix before calling this method.");
1026 return A->getDomainMap ();
1027}
1028
1029
1030template<class LocalOrdinal, class Node>
1031Teuchos::RCP<const typename Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_Int, Node> >::map_type>
1032Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_Int, Node> >::
1033getRangeMap () const
1034{
1035 Teuchos::RCP<const row_matrix_type> A = getMatrix();
1036 TEUCHOS_TEST_FOR_EXCEPTION(
1037 A.is_null (), std::runtime_error, "Ifpack2::Hypre::getRangeMap: The "
1038 "input matrix A is null. Please call setMatrix() with a nonnull input "
1039 "matrix before calling this method.");
1040 return A->getRangeMap ();
1041}
1042
1043
1044template<class LocalOrdinal, class Node>
1045void Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_Int, Node> >::setMatrix (const Teuchos::RCP<const row_matrix_type>& A)
1046{
1047 if (A.getRawPtr () != getMatrix().getRawPtr ()) {
1048 IsInitialized_ = false;
1049 IsComputed_ = false;
1050 A_ = A;
1051 }
1052}
1053
1054
1055template<class LocalOrdinal, class Node>
1056Teuchos::RCP<const typename Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_Int, Node> >::row_matrix_type>
1057Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_Int, Node> >::
1058getMatrix() const {
1059 return A_;
1060}
1061
1062
1063template<class LocalOrdinal, class Node>
1064bool Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_Int, Node> >::hasTransposeApply() const {
1065 return false;
1066}
1067
1068}// Ifpack2 namespace
1069
1070
1071#define IFPACK2_HYPRE_INSTANT(S,LO,GO,N) \
1072 template class Ifpack2::Hypre< Tpetra::RowMatrix<S, LO, GO, N> >;
1073
1074
1075#endif // HAVE_HYPRE && HAVE_MPI
1076#endif // IFPACK2_HYPRE_DEF_HPP
Ifpack2
Preconditioners and smoothers for Tpetra sparse matrices.
Definition Ifpack2_AdditiveSchwarz_decl.hpp:41
Ifpack2::GMRES
@ GMRES
Uses AztecOO's GMRES.
Definition Ifpack2_CondestType.hpp:20
Generated by doxygen 1.15.0