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#if defined(HAVE_IFPACK2_HYPRE_UNDERSCORE_KRYLOV_H)
24#include "_hypre_krylov.h"
25#else
26#include "krylov.h"
27#endif
28#include "_hypre_parcsr_mv.h"
29#include "_hypre_IJ_mv.h"
30#include "HYPRE_parcsr_mv.h"
31#include "HYPRE.h"
32
33using Teuchos::RCP;
34using Teuchos::rcp;
35using Teuchos::rcpFromRef;
36
37namespace Ifpack2 {
38
39template <class LocalOrdinal, class Node>
40Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_BigInt, Node> >::
41 Hypre(const Teuchos::RCP<const row_matrix_type> &A)
42 : A_(A)
43 , IsInitialized_(false)
44 , IsComputed_(false)
45 , NumInitialize_(0)
46 , NumCompute_(0)
47 , NumApply_(0)
48 , InitializeTime_(0.0)
49 , ComputeTime_(0.0)
50 , ApplyTime_(0.0)
51 , HypreA_(0)
52 , HypreG_(0)
53 , xHypre_(0)
54 , yHypre_(0)
55 , zHypre_(0)
56 , IsSolverCreated_(false)
57 , IsPrecondCreated_(false)
58 , SolveOrPrec_(Hypre_Is_Solver)
59 , NumFunsToCall_(0)
60 , SolverType_(PCG)
61 , PrecondType_(Euclid)
62 , UsePreconditioner_(false)
63 , Dump_(false) {}
64
65//==============================================================================
66template <class LocalOrdinal, class Node>
67Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_BigInt, Node> >::~Hypre() {
68 Destroy();
69}
70
71//==============================================================================
72template <class LocalOrdinal, class Node>
73void Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_BigInt, Node> >::Destroy() {
74 if (isInitialized()) {
75 IFPACK2_CHK_ERRV(HYPRE_IJMatrixDestroy(HypreA_));
76 IFPACK2_CHK_ERRV(HYPRE_IJVectorDestroy(XHypre_));
77 IFPACK2_CHK_ERRV(HYPRE_IJVectorDestroy(YHypre_));
78 }
79 if (IsSolverCreated_) {
80 IFPACK2_CHK_ERRV(SolverDestroyPtr_(Solver_));
81 }
82 if (IsPrecondCreated_) {
83 IFPACK2_CHK_ERRV(PrecondDestroyPtr_(Preconditioner_));
84 }
85
86 // Maxwell
87 if (HypreG_) {
88 IFPACK2_CHK_ERRV(HYPRE_IJMatrixDestroy(HypreG_));
89 }
90 if (xHypre_) {
91 IFPACK2_CHK_ERRV(HYPRE_IJVectorDestroy(xHypre_));
92 }
93 if (yHypre_) {
94 IFPACK2_CHK_ERRV(HYPRE_IJVectorDestroy(yHypre_));
95 }
96 if (zHypre_) {
97 IFPACK2_CHK_ERRV(HYPRE_IJVectorDestroy(zHypre_));
98 }
99} // Destroy()
100
101//==============================================================================
102template <class LocalOrdinal, class Node>
103void Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_BigInt, Node> >::initialize() {
104 const std::string timerName("Ifpack2::Hypre::initialize");
105 Teuchos::RCP<Teuchos::Time> timer = Teuchos::TimeMonitor::lookupCounter(timerName);
106 if (timer.is_null()) timer = Teuchos::TimeMonitor::getNewCounter(timerName);
107
108 if (IsInitialized_) return;
109 double startTime = timer->wallTime();
110 {
111 Teuchos::TimeMonitor timeMon(*timer);
112
113 MPI_Comm comm = *(Teuchos::rcp_dynamic_cast<const Teuchos::MpiComm<int> >(A_->getRowMap()->getComm())->getRawMpiComm());
114
115 // Check that RowMap and RangeMap are the same. While this could handle the
116 // case where RowMap and RangeMap are permutations, other Ifpack PCs don't
117 // handle this either.
118 if (!A_->getRowMap()->isSameAs(*A_->getRangeMap())) {
119 IFPACK2_CHK_ERRV(-1);
120 }
121 // Hypre expects the RowMap to be Linear.
122 if (A_->getRowMap()->isContiguous()) {
123 GloballyContiguousRowMap_ = A_->getRowMap();
124 GloballyContiguousColMap_ = A_->getColMap();
125 } else {
126 // Must create GloballyContiguous Maps for Hypre
127 if (A_->getDomainMap()->isSameAs(*A_->getRowMap())) {
128 Teuchos::RCP<const crs_matrix_type> Aconst = Teuchos::rcp_dynamic_cast<const crs_matrix_type>(A_);
129 GloballyContiguousColMap_ = MakeContiguousColumnMap(Aconst);
130 GloballyContiguousRowMap_ = rcp(new map_type(A_->getRowMap()->getGlobalNumElements(),
131 A_->getRowMap()->getLocalNumElements(), 0, A_->getRowMap()->getComm()));
132 } else {
133 throw std::runtime_error("Ifpack_Hypre: Unsupported map configuration: Row/Domain maps do not match");
134 }
135 }
136 // Next create vectors that will be used when ApplyInverse() is called
137 HYPRE_BigInt ilower = GloballyContiguousRowMap_->getMinGlobalIndex();
138 HYPRE_BigInt iupper = GloballyContiguousRowMap_->getMaxGlobalIndex();
139 // X in AX = Y
140 IFPACK2_CHK_ERRV(HYPRE_IJVectorCreate(comm, ilower, iupper, &XHypre_));
141 IFPACK2_CHK_ERRV(HYPRE_IJVectorSetObjectType(XHypre_, HYPRE_PARCSR));
142 IFPACK2_CHK_ERRV(HYPRE_IJVectorInitialize(XHypre_));
143 IFPACK2_CHK_ERRV(HYPRE_IJVectorAssemble(XHypre_));
144 IFPACK2_CHK_ERRV(HYPRE_IJVectorGetObject(XHypre_, (void **)&ParX_));
145 XVec_ = Teuchos::rcp((hypre_ParVector *)hypre_IJVectorObject(((hypre_IJVector *)XHypre_)), false);
146
147 // Y in AX = Y
148 IFPACK2_CHK_ERRV(HYPRE_IJVectorCreate(comm, ilower, iupper, &YHypre_));
149 IFPACK2_CHK_ERRV(HYPRE_IJVectorSetObjectType(YHypre_, HYPRE_PARCSR));
150 IFPACK2_CHK_ERRV(HYPRE_IJVectorInitialize(YHypre_));
151 IFPACK2_CHK_ERRV(HYPRE_IJVectorAssemble(YHypre_));
152 IFPACK2_CHK_ERRV(HYPRE_IJVectorGetObject(YHypre_, (void **)&ParY_));
153 YVec_ = Teuchos::rcp((hypre_ParVector *)hypre_IJVectorObject(((hypre_IJVector *)YHypre_)), false);
154
155 // Cache
156 VectorCache_.resize(A_->getRowMap()->getLocalNumElements());
157
158 // set flags
159 IsInitialized_ = true;
160 NumInitialize_++;
161 }
162 InitializeTime_ += (timer->wallTime() - startTime);
163} // Initialize()
164
165//==============================================================================
166template <class LocalOrdinal, class Node>
167void Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_BigInt, Node> >::setParameters(const Teuchos::ParameterList &list) {
168 std::map<std::string, Hypre_Solver> solverMap;
169 solverMap["BoomerAMG"] = BoomerAMG;
170 solverMap["ParaSails"] = ParaSails;
171 solverMap["Euclid"] = Euclid;
172 solverMap["AMS"] = AMS;
173 solverMap["Hybrid"] = Hybrid;
174 solverMap["PCG"] = PCG;
175 solverMap["GMRES"] = GMRES;
176 solverMap["FlexGMRES"] = FlexGMRES;
177 solverMap["LGMRES"] = LGMRES;
178 solverMap["BiCGSTAB"] = BiCGSTAB;
179
180 std::map<std::string, Hypre_Chooser> chooserMap;
181 chooserMap["Solver"] = Hypre_Is_Solver;
182 chooserMap["Preconditioner"] = Hypre_Is_Preconditioner;
183
184 List_ = list;
185 Hypre_Solver solType;
186 if (list.isType<std::string>("hypre: Solver"))
187 solType = solverMap[list.get<std::string>("hypre: Solver")];
188 else if (list.isParameter("hypre: Solver"))
189 solType = (Hypre_Solver)list.get<int>("hypre: Solver");
190 else
191 solType = PCG;
192 SolverType_ = solType;
193 Hypre_Solver precType;
194 if (list.isType<std::string>("hypre: Preconditioner"))
195 precType = solverMap[list.get<std::string>("hypre: Preconditioner")];
196 else if (list.isParameter("hypre: Preconditioner"))
197 precType = (Hypre_Solver)list.get<int>("hypre: Preconditioner");
198 else
199 precType = Euclid;
200 PrecondType_ = precType;
201 Hypre_Chooser chooser;
202 if (list.isType<std::string>("hypre: SolveOrPrecondition"))
203 chooser = chooserMap[list.get<std::string>("hypre: SolveOrPrecondition")];
204 else if (list.isParameter("hypre: SolveOrPrecondition"))
205 chooser = (Hypre_Chooser)list.get<int>("hypre: SolveOrPrecondition");
206 else
207 chooser = Hypre_Is_Solver;
208 SolveOrPrec_ = chooser;
209 bool SetPrecond = list.isParameter("hypre: SetPreconditioner") ? list.get<bool>("hypre: SetPreconditioner") : false;
210 IFPACK2_CHK_ERR(SetParameter(SetPrecond));
211 int NumFunctions = list.isParameter("hypre: NumFunctions") ? list.get<int>("hypre: NumFunctions") : 0;
212 FunsToCall_.clear();
213 NumFunsToCall_ = 0;
214 if (NumFunctions > 0) {
215 RCP<FunctionParameter> *params = list.get<RCP<FunctionParameter> *>("hypre: Functions");
216 for (int i = 0; i < NumFunctions; i++) {
217 IFPACK2_CHK_ERR(AddFunToList(params[i]));
218 }
219 }
220
221 if (list.isSublist("hypre: Solver functions")) {
222 Teuchos::ParameterList solverList = list.sublist("hypre: Solver functions");
223 for (auto it = solverList.begin(); it != solverList.end(); ++it) {
224 std::string funct_name = it->first;
225 if (it->second.isType<HYPRE_Int>()) {
226 IFPACK2_CHK_ERR(AddFunToList(rcp(new FunctionParameter(Hypre_Is_Solver, funct_name, Teuchos::getValue<HYPRE_Int>(it->second)))));
227 } else if (!std::is_same<HYPRE_Int, int>::value && it->second.isType<int>()) {
228 IFPACK2_CHK_ERR(AddFunToList(rcp(new FunctionParameter(Hypre_Is_Solver, funct_name, Teuchos::as<HYPRE_Int>(Teuchos::getValue<int>(it->second))))));
229 } else if (it->second.isType<HYPRE_Real>()) {
230 IFPACK2_CHK_ERR(AddFunToList(rcp(new FunctionParameter(Hypre_Is_Solver, funct_name, Teuchos::getValue<HYPRE_Real>(it->second)))));
231 } else {
232 IFPACK2_CHK_ERR(-1);
233 }
234 }
235 }
236
237 if (list.isSublist("hypre: Preconditioner functions")) {
238 Teuchos::ParameterList precList = list.sublist("hypre: Preconditioner functions");
239 for (auto it = precList.begin(); it != precList.end(); ++it) {
240 std::string funct_name = it->first;
241 if (it->second.isType<HYPRE_Int>()) {
242 IFPACK2_CHK_ERR(AddFunToList(rcp(new FunctionParameter(Hypre_Is_Preconditioner, funct_name, Teuchos::getValue<HYPRE_Int>(it->second)))));
243 } else if (!std::is_same<HYPRE_Int, int>::value && it->second.isType<int>()) {
244 IFPACK2_CHK_ERR(AddFunToList(rcp(new FunctionParameter(Hypre_Is_Preconditioner, funct_name, Teuchos::as<HYPRE_Int>(Teuchos::getValue<int>(it->second))))));
245 } else if (it->second.isType<HYPRE_Real>()) {
246 IFPACK2_CHK_ERR(AddFunToList(rcp(new FunctionParameter(Hypre_Is_Preconditioner, funct_name, Teuchos::getValue<HYPRE_Real>(it->second)))));
247 } else if (it->second.isList()) {
248 Teuchos::ParameterList pl = Teuchos::getValue<Teuchos::ParameterList>(it->second);
249 if (FunctionParameter::isFuncIntInt(funct_name)) {
250 HYPRE_Int arg0 = pl.get<HYPRE_Int>("arg 0");
251 HYPRE_Int arg1 = pl.get<HYPRE_Int>("arg 1");
252 IFPACK2_CHK_ERR(AddFunToList(rcp(new FunctionParameter(Hypre_Is_Preconditioner, funct_name, arg0, arg1))));
253 } else if (FunctionParameter::isFuncIntIntDoubleDouble(funct_name)) {
254 HYPRE_Int arg0 = pl.get<HYPRE_Int>("arg 0");
255 HYPRE_Int arg1 = pl.get<HYPRE_Int>("arg 1");
256 HYPRE_Real arg2 = pl.get<HYPRE_Real>("arg 2");
257 HYPRE_Real arg3 = pl.get<HYPRE_Real>("arg 3");
258 IFPACK2_CHK_ERR(AddFunToList(rcp(new FunctionParameter(Hypre_Is_Preconditioner, funct_name, arg0, arg1, arg2, arg3))));
259 } else if (FunctionParameter::isFuncIntIntIntDoubleIntInt(funct_name)) {
260 HYPRE_Int arg0 = pl.get<HYPRE_Int>("arg 0");
261 HYPRE_Int arg1 = pl.get<HYPRE_Int>("arg 1");
262 HYPRE_Int arg2 = pl.get<HYPRE_Int>("arg 2");
263 HYPRE_Real arg3 = pl.get<HYPRE_Real>("arg 3");
264 HYPRE_Int arg4 = pl.get<HYPRE_Int>("arg 4");
265 HYPRE_Int arg5 = pl.get<HYPRE_Int>("arg 5");
266 IFPACK2_CHK_ERR(AddFunToList(rcp(new FunctionParameter(Hypre_Is_Preconditioner, funct_name, arg0, arg1, arg2, arg3, arg4, arg5))));
267 } else {
268 IFPACK2_CHK_ERR(-1);
269 }
270 }
271 }
272 }
273
274 if (list.isSublist("Coordinates") && list.sublist("Coordinates").isType<Teuchos::RCP<multivector_type> >("Coordinates"))
275 Coords_ = list.sublist("Coordinates").get<Teuchos::RCP<multivector_type> >("Coordinates");
276 if (list.isSublist("Operators") && list.sublist("Operators").isType<Teuchos::RCP<const crs_matrix_type> >("G"))
277 G_ = list.sublist("Operators").get<Teuchos::RCP<const crs_matrix_type> >("G");
278
279 Dump_ = list.isParameter("hypre: Dump") ? list.get<bool>("hypre: Dump") : false;
280} // setParameters()
281
282//==============================================================================
283template <class LocalOrdinal, class Node>
284int Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_BigInt, Node> >::AddFunToList(RCP<FunctionParameter> NewFun) {
285 NumFunsToCall_ = NumFunsToCall_ + 1;
286 FunsToCall_.resize(NumFunsToCall_);
287 FunsToCall_[NumFunsToCall_ - 1] = NewFun;
288 return 0;
289} // AddFunToList()
290
291//==============================================================================
292template <class LocalOrdinal, class Node>
293int Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_BigInt, Node> >::SetParameter(Hypre_Chooser chooser, HYPRE_Int (*pt2Func)(HYPRE_Solver, HYPRE_Int), HYPRE_Int parameter) {
294 RCP<FunctionParameter> temp = rcp(new FunctionParameter(chooser, pt2Func, parameter));
295 IFPACK2_CHK_ERR(AddFunToList(temp));
296 return 0;
297} // SetParameter() - int function pointer
298
299//=============================================================================
300template <class LocalOrdinal, class Node>
301int Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_BigInt, Node> >::SetParameter(Hypre_Chooser chooser, HYPRE_Int (*pt2Func)(HYPRE_Solver, HYPRE_Real), HYPRE_Real parameter) {
302 RCP<FunctionParameter> temp = rcp(new FunctionParameter(chooser, pt2Func, parameter));
303 IFPACK2_CHK_ERR(AddFunToList(temp));
304 return 0;
305} // SetParameter() - double function pointer
306
307//==============================================================================
308template <class LocalOrdinal, class Node>
309int Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_BigInt, Node> >::SetParameter(Hypre_Chooser chooser, HYPRE_Int (*pt2Func)(HYPRE_Solver, HYPRE_Real, HYPRE_Int), HYPRE_Real parameter1, HYPRE_Int parameter2) {
310 RCP<FunctionParameter> temp = rcp(new FunctionParameter(chooser, pt2Func, parameter1, parameter2));
311 IFPACK2_CHK_ERR(AddFunToList(temp));
312 return 0;
313} // SetParameter() - double,int function pointer
314
315//==============================================================================
316template <class LocalOrdinal, class Node>
317int Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_BigInt, Node> >::SetParameter(Hypre_Chooser chooser, HYPRE_Int (*pt2Func)(HYPRE_Solver, HYPRE_Int, HYPRE_Real), HYPRE_Int parameter1, HYPRE_Real parameter2) {
318 RCP<FunctionParameter> temp = rcp(new FunctionParameter(chooser, pt2Func, parameter1, parameter2));
319 IFPACK2_CHK_ERR(AddFunToList(temp));
320 return 0;
321} // SetParameter() - int,double function pointer
322
323//==============================================================================
324template <class LocalOrdinal, class Node>
325int Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_BigInt, Node> >::SetParameter(Hypre_Chooser chooser, HYPRE_Int (*pt2Func)(HYPRE_Solver, HYPRE_Int, HYPRE_Int), HYPRE_Int parameter1, HYPRE_Int parameter2) {
326 RCP<FunctionParameter> temp = rcp(new FunctionParameter(chooser, pt2Func, parameter1, parameter2));
327 IFPACK2_CHK_ERR(AddFunToList(temp));
328 return 0;
329} // SetParameter() int,int function pointer
330
331//==============================================================================
332template <class LocalOrdinal, class Node>
333int Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_BigInt, Node> >::SetParameter(Hypre_Chooser chooser, HYPRE_Int (*pt2Func)(HYPRE_Solver, HYPRE_Real *), HYPRE_Real *parameter) {
334 RCP<FunctionParameter> temp = rcp(new FunctionParameter(chooser, pt2Func, parameter));
335 IFPACK2_CHK_ERR(AddFunToList(temp));
336 return 0;
337} // SetParameter() - double* function pointer
338
339//==============================================================================
340template <class LocalOrdinal, class Node>
341int Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_BigInt, Node> >::SetParameter(Hypre_Chooser chooser, HYPRE_Int (*pt2Func)(HYPRE_Solver, HYPRE_Int *), HYPRE_Int *parameter) {
342 RCP<FunctionParameter> temp = rcp(new FunctionParameter(chooser, pt2Func, parameter));
343 IFPACK2_CHK_ERR(AddFunToList(temp));
344 return 0;
345} // SetParameter() - int* function pointer
346
347//==============================================================================
348template <class LocalOrdinal, class Node>
349int Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_BigInt, Node> >::SetParameter(Hypre_Chooser chooser, HYPRE_Int (*pt2Func)(HYPRE_Solver, HYPRE_Int **), HYPRE_Int **parameter) {
350 RCP<FunctionParameter> temp = rcp(new FunctionParameter(chooser, pt2Func, parameter));
351 IFPACK2_CHK_ERR(AddFunToList(temp));
352 return 0;
353} // SetParameter() - int** function pointer
354
355//==============================================================================
356template <class LocalOrdinal, class Node>
357int Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_BigInt, Node> >::SetParameter(Hypre_Chooser chooser, Hypre_Solver solver) {
358 if (chooser == Hypre_Is_Solver) {
359 SolverType_ = solver;
360 } else {
361 PrecondType_ = solver;
362 }
363 return 0;
364} // SetParameter() - set type of solver
365
366//==============================================================================
367template <class LocalOrdinal, class Node>
368int Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_BigInt, Node> >::SetDiscreteGradient(Teuchos::RCP<const crs_matrix_type> G) {
369 using LO = local_ordinal_type;
370 using GO = global_ordinal_type;
371 // using SC = scalar_type;
372
373 // Sanity check
374 if (!A_->getRowMap()->isSameAs(*G->getRowMap()))
375 throw std::runtime_error("Hypre<Tpetra::RowMatrix<double, HYPRE_Int, long long, Node>: Edge map mismatch: A and discrete gradient");
376
377 // Get the maps for the nodes (assuming the edge map from A is OK);
378 GloballyContiguousNodeRowMap_ = rcp(new map_type(G->getDomainMap()->getGlobalNumElements(),
379 G->getDomainMap()->getLocalNumElements(), 0, A_->getRowMap()->getComm()));
380 GloballyContiguousNodeColMap_ = MakeContiguousColumnMap(G);
381
382 // Start building G
383 MPI_Comm comm = *(Teuchos::rcp_dynamic_cast<const Teuchos::MpiComm<int> >(A_->getRowMap()->getComm())->getRawMpiComm());
384 GO ilower = GloballyContiguousRowMap_->getMinGlobalIndex();
385 GO iupper = GloballyContiguousRowMap_->getMaxGlobalIndex();
386 GO jlower = GloballyContiguousNodeRowMap_->getMinGlobalIndex();
387 GO jupper = GloballyContiguousNodeRowMap_->getMaxGlobalIndex();
388 IFPACK2_CHK_ERR(HYPRE_IJMatrixCreate(comm, ilower, iupper, jlower, jupper, &HypreG_));
389 IFPACK2_CHK_ERR(HYPRE_IJMatrixSetObjectType(HypreG_, HYPRE_PARCSR));
390 IFPACK2_CHK_ERR(HYPRE_IJMatrixInitialize(HypreG_));
391
392 std::vector<GO> new_indices(G->getLocalMaxNumRowEntries());
393 for (LO i = 0; i < (LO)G->getLocalNumRows(); i++) {
394 typename crs_matrix_type::values_host_view_type values;
395 typename crs_matrix_type::local_inds_host_view_type indices;
396 G->getLocalRowView(i, indices, values);
397 for (LO j = 0; j < (LO)indices.extent(0); j++) {
398 new_indices[j] = GloballyContiguousNodeColMap_->getGlobalElement(indices(j));
399 }
400 HYPRE_BigInt GlobalRow[1];
401 HYPRE_Int numEntries = (HYPRE_Int)indices.extent(0);
402 GlobalRow[0] = GloballyContiguousRowMap_->getGlobalElement(i);
403 IFPACK2_CHK_ERR(HYPRE_IJMatrixSetValues(HypreG_, 1, &numEntries, GlobalRow, new_indices.data(), values.data()));
404 }
405 IFPACK2_CHK_ERR(HYPRE_IJMatrixAssemble(HypreG_));
406 IFPACK2_CHK_ERR(HYPRE_IJMatrixGetObject(HypreG_, (void **)&ParMatrixG_));
407
408 if (Dump_)
409 HYPRE_ParCSRMatrixPrint(ParMatrixG_, "G.mat");
410
411 if (SolverType_ == AMS)
412 HYPRE_AMSSetDiscreteGradient(Solver_, ParMatrixG_);
413 if (PrecondType_ == AMS)
414 HYPRE_AMSSetDiscreteGradient(Preconditioner_, ParMatrixG_);
415 return 0;
416} // SetDiscreteGradient()
417
418//==============================================================================
419template <class LocalOrdinal, class Node>
420int Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_BigInt, Node> >::SetCoordinates(Teuchos::RCP<multivector_type> coords) {
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_BigInt, 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_BigInt, 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_BigInt, 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 Tpetra'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 Tpetra 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//==============================================================================
626template <class LocalOrdinal, class Node>
627void Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_BigInt, Node> >::applyMat(const Tpetra::MultiVector<scalar_type, local_ordinal_type, global_ordinal_type, node_type> &X,
628 Tpetra::MultiVector<scalar_type, local_ordinal_type, global_ordinal_type, node_type> &Y,
629 Teuchos::ETransp mode) const {
630 A_->apply(X, Y, mode);
631} // applyMat()
632
633//==============================================================================
634template <class LocalOrdinal, class Node>
635std::string Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_BigInt, Node> >::description() const {
636 std::ostringstream out;
637
638 // Output is a valid YAML dictionary in flow style. If you don't
639 // like everything on a single line, you should call describe()
640 // instead.
641 out << "\"Ifpack2::Hypre\": {";
642 out << "Initialized: " << (isInitialized() ? "true" : "false") << ", "
643 << "Computed: " << (isComputed() ? "true" : "false") << ", ";
644
645 if (A_.is_null()) {
646 out << "Matrix: null";
647 } else {
648 out << "Global matrix dimensions: ["
649 << A_->getGlobalNumRows() << ", "
650 << A_->getGlobalNumCols() << "]"
651 << ", Global nnz: " << A_->getGlobalNumEntries();
652 }
653
654 out << "}";
655 return out.str();
656} // description()
657
658//==============================================================================
659template <class LocalOrdinal, class Node>
660void Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_BigInt, Node> >::describe(Teuchos::FancyOStream &os, const Teuchos::EVerbosityLevel verbLevel) const {
661 using std::endl;
662 os << endl;
663 os << "================================================================================" << endl;
664 os << "Ifpack2::Hypre: " << endl
665 << endl;
666 os << "Using " << A_->getComm()->getSize() << " processors." << endl;
667 os << "Global number of rows = " << A_->getGlobalNumRows() << endl;
668 os << "Global number of nonzeros = " << A_->getGlobalNumEntries() << endl;
669 // os << "Condition number estimate = " << Condest() << endl;
670 os << endl;
671 os << "Phase # calls Total Time (s)" << endl;
672 os << "----- ------- --------------" << endl;
673 os << "Initialize() " << std::setw(5) << NumInitialize_
674 << " " << std::setw(15) << InitializeTime_ << endl;
675 os << "Compute() " << std::setw(5) << NumCompute_
676 << " " << std::setw(15) << ComputeTime_ << endl;
677 os << "ApplyInverse() " << std::setw(5) << NumApply_
678 << " " << std::setw(15) << ApplyTime_ << endl;
679 os << "================================================================================" << endl;
680 os << endl;
681} // description
682
683//==============================================================================
684template <class LocalOrdinal, class Node>
685int Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_BigInt, Node> >::SetSolverType(Hypre_Solver Solver) {
686 switch (Solver) {
687 case BoomerAMG:
688 if (IsSolverCreated_) {
689 SolverDestroyPtr_(Solver_);
690 IsSolverCreated_ = false;
691 }
692 SolverCreatePtr_ = &Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_BigInt, Node> >::Hypre_BoomerAMGCreate;
693 SolverDestroyPtr_ = &HYPRE_BoomerAMGDestroy;
694 SolverSetupPtr_ = &HYPRE_BoomerAMGSetup;
695 SolverPrecondPtr_ = NULL;
696 SolverSolvePtr_ = &HYPRE_BoomerAMGSolve;
697 break;
698 case AMS:
699 if (IsSolverCreated_) {
700 SolverDestroyPtr_(Solver_);
701 IsSolverCreated_ = false;
702 }
703 SolverCreatePtr_ = &Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_BigInt, Node> >::Hypre_AMSCreate;
704 SolverDestroyPtr_ = &HYPRE_AMSDestroy;
705 SolverSetupPtr_ = &HYPRE_AMSSetup;
706 SolverSolvePtr_ = &HYPRE_AMSSolve;
707 SolverPrecondPtr_ = NULL;
708 break;
709 case Hybrid:
710 if (IsSolverCreated_) {
711 SolverDestroyPtr_(Solver_);
712 IsSolverCreated_ = false;
713 }
714 SolverCreatePtr_ = &Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_BigInt, Node> >::Hypre_ParCSRHybridCreate;
715 SolverDestroyPtr_ = &HYPRE_ParCSRHybridDestroy;
716 SolverSetupPtr_ = &HYPRE_ParCSRHybridSetup;
717 SolverSolvePtr_ = &HYPRE_ParCSRHybridSolve;
718 SolverPrecondPtr_ = &HYPRE_ParCSRHybridSetPrecond;
719 break;
720 case PCG:
721 if (IsSolverCreated_) {
722 SolverDestroyPtr_(Solver_);
723 IsSolverCreated_ = false;
724 }
725 SolverCreatePtr_ = &Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_BigInt, Node> >::Hypre_ParCSRPCGCreate;
726 SolverDestroyPtr_ = &HYPRE_ParCSRPCGDestroy;
727 SolverSetupPtr_ = &HYPRE_ParCSRPCGSetup;
728 SolverSolvePtr_ = &HYPRE_ParCSRPCGSolve;
729 SolverPrecondPtr_ = &HYPRE_ParCSRPCGSetPrecond;
730 break;
731 case GMRES:
732 if (IsSolverCreated_) {
733 SolverDestroyPtr_(Solver_);
734 IsSolverCreated_ = false;
735 }
736 SolverCreatePtr_ = &Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_BigInt, Node> >::Hypre_ParCSRGMRESCreate;
737 SolverDestroyPtr_ = &HYPRE_ParCSRGMRESDestroy;
738 SolverSetupPtr_ = &HYPRE_ParCSRGMRESSetup;
739 SolverPrecondPtr_ = &HYPRE_ParCSRGMRESSetPrecond;
740 break;
741 case FlexGMRES:
742 if (IsSolverCreated_) {
743 SolverDestroyPtr_(Solver_);
744 IsSolverCreated_ = false;
745 }
746 SolverCreatePtr_ = &Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_BigInt, Node> >::Hypre_ParCSRFlexGMRESCreate;
747 SolverDestroyPtr_ = &HYPRE_ParCSRFlexGMRESDestroy;
748 SolverSetupPtr_ = &HYPRE_ParCSRFlexGMRESSetup;
749 SolverSolvePtr_ = &HYPRE_ParCSRFlexGMRESSolve;
750 SolverPrecondPtr_ = &HYPRE_ParCSRFlexGMRESSetPrecond;
751 break;
752 case LGMRES:
753 if (IsSolverCreated_) {
754 SolverDestroyPtr_(Solver_);
755 IsSolverCreated_ = false;
756 }
757 SolverCreatePtr_ = &Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_BigInt, Node> >::Hypre_ParCSRLGMRESCreate;
758 SolverDestroyPtr_ = &HYPRE_ParCSRLGMRESDestroy;
759 SolverSetupPtr_ = &HYPRE_ParCSRLGMRESSetup;
760 SolverSolvePtr_ = &HYPRE_ParCSRLGMRESSolve;
761 SolverPrecondPtr_ = &HYPRE_ParCSRLGMRESSetPrecond;
762 break;
763 case BiCGSTAB:
764 if (IsSolverCreated_) {
765 SolverDestroyPtr_(Solver_);
766 IsSolverCreated_ = false;
767 }
768 SolverCreatePtr_ = &Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_BigInt, Node> >::Hypre_ParCSRBiCGSTABCreate;
769 SolverDestroyPtr_ = &HYPRE_ParCSRBiCGSTABDestroy;
770 SolverSetupPtr_ = &HYPRE_ParCSRBiCGSTABSetup;
771 SolverSolvePtr_ = &HYPRE_ParCSRBiCGSTABSolve;
772 SolverPrecondPtr_ = &HYPRE_ParCSRBiCGSTABSetPrecond;
773 break;
774 default:
775 return -1;
776 }
777 CreateSolver();
778 return 0;
779} // SetSolverType()
780
781//==============================================================================
782template <class LocalOrdinal, class Node>
783int Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_BigInt, Node> >::SetPrecondType(Hypre_Solver Precond) {
784 switch (Precond) {
785 case BoomerAMG:
786 if (IsPrecondCreated_) {
787 PrecondDestroyPtr_(Preconditioner_);
788 IsPrecondCreated_ = false;
789 }
790 PrecondCreatePtr_ = &Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_BigInt, Node> >::Hypre_BoomerAMGCreate;
791 PrecondDestroyPtr_ = &HYPRE_BoomerAMGDestroy;
792 PrecondSetupPtr_ = &HYPRE_BoomerAMGSetup;
793 PrecondSolvePtr_ = &HYPRE_BoomerAMGSolve;
794 break;
795 case ParaSails:
796 if (IsPrecondCreated_) {
797 PrecondDestroyPtr_(Preconditioner_);
798 IsPrecondCreated_ = false;
799 }
800 PrecondCreatePtr_ = &Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_BigInt, Node> >::Hypre_ParaSailsCreate;
801 PrecondDestroyPtr_ = &HYPRE_ParaSailsDestroy;
802 PrecondSetupPtr_ = &HYPRE_ParaSailsSetup;
803 PrecondSolvePtr_ = &HYPRE_ParaSailsSolve;
804 break;
805 case Euclid:
806 if (IsPrecondCreated_) {
807 PrecondDestroyPtr_(Preconditioner_);
808 IsPrecondCreated_ = false;
809 }
810 PrecondCreatePtr_ = &Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_BigInt, Node> >::Hypre_EuclidCreate;
811 PrecondDestroyPtr_ = &HYPRE_EuclidDestroy;
812 PrecondSetupPtr_ = &HYPRE_EuclidSetup;
813 PrecondSolvePtr_ = &HYPRE_EuclidSolve;
814 break;
815 case AMS:
816 if (IsPrecondCreated_) {
817 PrecondDestroyPtr_(Preconditioner_);
818 IsPrecondCreated_ = false;
819 }
820 PrecondCreatePtr_ = &Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_BigInt, Node> >::Hypre_AMSCreate;
821 PrecondDestroyPtr_ = &HYPRE_AMSDestroy;
822 PrecondSetupPtr_ = &HYPRE_AMSSetup;
823 PrecondSolvePtr_ = &HYPRE_AMSSolve;
824 break;
825 default:
826 return -1;
827 }
828 CreatePrecond();
829 return 0;
830
831} // SetPrecondType()
832
833//==============================================================================
834template <class LocalOrdinal, class Node>
835int Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_BigInt, Node> >::CreateSolver() {
836 MPI_Comm comm;
837 HYPRE_ParCSRMatrixGetComm(ParMatrix_, &comm);
838 int ierr = (this->*SolverCreatePtr_)(comm, &Solver_);
839 IsSolverCreated_ = true;
840 return ierr;
841} // CreateSolver()
842
843//==============================================================================
844template <class LocalOrdinal, class Node>
845int Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_BigInt, Node> >::CreatePrecond() {
846 MPI_Comm comm;
847 HYPRE_ParCSRMatrixGetComm(ParMatrix_, &comm);
848 int ierr = (this->*PrecondCreatePtr_)(comm, &Preconditioner_);
849 IsPrecondCreated_ = true;
850 return ierr;
851} // CreatePrecond()
852
853//==============================================================================
854template <class LocalOrdinal, class Node>
855int Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_BigInt, Node> >::CopyTpetraToHypre() {
856 using LO = local_ordinal_type;
857 using GO = global_ordinal_type;
858 // using SC = scalar_type;
859
860 Teuchos::RCP<const crs_matrix_type> Matrix = Teuchos::rcp_dynamic_cast<const crs_matrix_type>(A_);
861 if (Matrix.is_null())
862 throw std::runtime_error("Hypre<Tpetra::RowMatrix<double, LocalOrdinal, HYPRE_BigInt, Node>: Unsupported matrix configuration: Tpetra::CrsMatrix required");
863
864 std::vector<HYPRE_BigInt> new_indices(Matrix->getLocalMaxNumRowEntries());
865 for (LO i = 0; i < (LO)Matrix->getLocalNumRows(); i++) {
866 typename crs_matrix_type::values_host_view_type values;
867 typename crs_matrix_type::local_inds_host_view_type indices;
868 Matrix->getLocalRowView(i, indices, values);
869 for (LO j = 0; j < (LO)indices.extent(0); j++) {
870 new_indices[j] = GloballyContiguousColMap_->getGlobalElement(indices(j));
871 }
872 HYPRE_BigInt GlobalRow[1];
873 HYPRE_Int numEntries = (HYPRE_Int)indices.extent(0);
874 GlobalRow[0] = GloballyContiguousRowMap_->getGlobalElement(i);
875 IFPACK2_CHK_ERR(HYPRE_IJMatrixSetValues(HypreA_, 1, &numEntries, GlobalRow, new_indices.data(), values.data()));
876 }
877 IFPACK2_CHK_ERR(HYPRE_IJMatrixAssemble(HypreA_));
878 IFPACK2_CHK_ERR(HYPRE_IJMatrixGetObject(HypreA_, (void **)&ParMatrix_));
879 if (Dump_)
880 HYPRE_ParCSRMatrixPrint(ParMatrix_, "A.mat");
881 return 0;
882} // CopyTpetraToHypre()
883
884//==============================================================================
885template <class LocalOrdinal, class Node>
886HYPRE_Int Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_BigInt, Node> >::Hypre_BoomerAMGCreate(MPI_Comm /*comm*/, HYPRE_Solver *solver) { return HYPRE_BoomerAMGCreate(solver); }
887
888//==============================================================================
889template <class LocalOrdinal, class Node>
890HYPRE_Int Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_BigInt, Node> >::Hypre_ParaSailsCreate(MPI_Comm comm, HYPRE_Solver *solver) { return HYPRE_ParaSailsCreate(comm, solver); }
891
892//==============================================================================
893template <class LocalOrdinal, class Node>
894HYPRE_Int Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_BigInt, Node> >::Hypre_EuclidCreate(MPI_Comm comm, HYPRE_Solver *solver) { return HYPRE_EuclidCreate(comm, solver); }
895
896//==============================================================================
897template <class LocalOrdinal, class Node>
898HYPRE_Int Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_BigInt, Node> >::Hypre_AMSCreate(MPI_Comm /*comm*/, HYPRE_Solver *solver) { return HYPRE_AMSCreate(solver); }
899
900//==============================================================================
901template <class LocalOrdinal, class Node>
902HYPRE_Int Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_BigInt, Node> >::Hypre_ParCSRHybridCreate(MPI_Comm /*comm*/, HYPRE_Solver *solver) { return HYPRE_ParCSRHybridCreate(solver); }
903
904//==============================================================================
905template <class LocalOrdinal, class Node>
906HYPRE_Int Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_BigInt, Node> >::Hypre_ParCSRPCGCreate(MPI_Comm comm, HYPRE_Solver *solver) { return HYPRE_ParCSRPCGCreate(comm, solver); }
907
908//==============================================================================
909template <class LocalOrdinal, class Node>
910HYPRE_Int Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_BigInt, Node> >::Hypre_ParCSRGMRESCreate(MPI_Comm comm, HYPRE_Solver *solver) { return HYPRE_ParCSRGMRESCreate(comm, solver); }
911
912//==============================================================================
913template <class LocalOrdinal, class Node>
914HYPRE_Int Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_BigInt, Node> >::Hypre_ParCSRFlexGMRESCreate(MPI_Comm comm, HYPRE_Solver *solver) { return HYPRE_ParCSRFlexGMRESCreate(comm, solver); }
915
916//==============================================================================
917template <class LocalOrdinal, class Node>
918HYPRE_Int Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_BigInt, Node> >::Hypre_ParCSRLGMRESCreate(MPI_Comm comm, HYPRE_Solver *solver) { return HYPRE_ParCSRLGMRESCreate(comm, solver); }
919
920//==============================================================================
921template <class LocalOrdinal, class Node>
922HYPRE_Int Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_BigInt, Node> >::Hypre_ParCSRBiCGSTABCreate(MPI_Comm comm, HYPRE_Solver *solver) { return HYPRE_ParCSRBiCGSTABCreate(comm, solver); }
923
924//==============================================================================
925template <class LocalOrdinal, class Node>
926Teuchos::RCP<const Tpetra::Map<LocalOrdinal, HYPRE_BigInt, Node> >
927Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_BigInt, Node> >::MakeContiguousColumnMap(Teuchos::RCP<const crs_matrix_type> &Matrix) const {
928 using import_type = Tpetra::Import<local_ordinal_type, global_ordinal_type, node_type>;
929 using go_vector_type = Tpetra::Vector<global_ordinal_type, local_ordinal_type, global_ordinal_type, node_type>;
930
931 // Must create GloballyContiguous DomainMap (which is a permutation of Matrix_'s
932 // DomainMap) and the corresponding permuted ColumnMap.
933 // Tpetra_GID ---------> LID ----------> HYPRE_GID
934 // via DomainMap.LID() via GloballyContiguousDomainMap.GID()
935 if (Matrix.is_null())
936 throw std::runtime_error("Hypre<Tpetra::RowMatrix<HYPRE_Real, HYPRE_Int, long long, Node>: Unsupported matrix configuration: Tpetra::CrsMatrix required");
937 RCP<const map_type> DomainMap = Matrix->getDomainMap();
938 RCP<const map_type> ColumnMap = Matrix->getColMap();
939 RCP<const import_type> importer = Matrix->getGraph()->getImporter();
940
941 if (DomainMap->isContiguous()) {
942 // If the domain map is linear, then we can just use the column map as is.
943 return ColumnMap;
944 } else {
945 // The domain map isn't linear, so we need a new domain map
946 Teuchos::RCP<map_type> ContiguousDomainMap = rcp(new map_type(DomainMap->getGlobalNumElements(),
947 DomainMap->getLocalNumElements(), 0, DomainMap->getComm()));
948 if (importer) {
949 // If there's an importer then we can use it to get a new column map
950 go_vector_type MyGIDsHYPRE(DomainMap, ContiguousDomainMap->getLocalElementList());
951
952 // import the HYPRE GIDs
953 go_vector_type ColGIDsHYPRE(ColumnMap);
954 ColGIDsHYPRE.doImport(MyGIDsHYPRE, *importer, Tpetra::INSERT);
955
956 // Make a HYPRE numbering-based column map.
957 return Teuchos::rcp(new map_type(ColumnMap->getGlobalNumElements(), ColGIDsHYPRE.getDataNonConst()(), 0, ColumnMap->getComm()));
958 } else {
959 // The problem has matching domain/column maps, and somehow the domain map isn't linear, so just use the new domain map
960 return Teuchos::rcp(new map_type(ColumnMap->getGlobalNumElements(), ContiguousDomainMap->getLocalElementList(), 0, ColumnMap->getComm()));
961 }
962 }
963}
964
965template <class LocalOrdinal, class Node>
966int Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_BigInt, Node> >::getNumInitialize() const {
967 return NumInitialize_;
968}
969
970template <class LocalOrdinal, class Node>
971int Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_BigInt, Node> >::getNumCompute() const {
972 return NumCompute_;
973}
974
975template <class LocalOrdinal, class Node>
976int Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_BigInt, Node> >::getNumApply() const {
977 return NumApply_;
978}
979
980template <class LocalOrdinal, class Node>
981double Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_BigInt, Node> >::getInitializeTime() const {
982 return InitializeTime_;
983}
984
985template <class LocalOrdinal, class Node>
986double Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_BigInt, Node> >::getComputeTime() const {
987 return ComputeTime_;
988}
989
990template <class LocalOrdinal, class Node>
991double Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_BigInt, Node> >::getApplyTime() const {
992 return ApplyTime_;
993}
994
995template <class LocalOrdinal, class Node>
996Teuchos::RCP<const typename Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_BigInt, Node> >::map_type>
997Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_BigInt, Node> >::
998 getDomainMap() const {
999 Teuchos::RCP<const row_matrix_type> A = getMatrix();
1000 TEUCHOS_TEST_FOR_EXCEPTION(
1001 A.is_null(), std::runtime_error,
1002 "Ifpack2::Hypre::getDomainMap: The "
1003 "input matrix A is null. Please call setMatrix() with a nonnull input "
1004 "matrix before calling this method.");
1005 return A->getDomainMap();
1006}
1007
1008template <class LocalOrdinal, class Node>
1009Teuchos::RCP<const typename Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_BigInt, Node> >::map_type>
1010Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_BigInt, Node> >::
1011 getRangeMap() const {
1012 Teuchos::RCP<const row_matrix_type> A = getMatrix();
1013 TEUCHOS_TEST_FOR_EXCEPTION(
1014 A.is_null(), std::runtime_error,
1015 "Ifpack2::Hypre::getRangeMap: The "
1016 "input matrix A is null. Please call setMatrix() with a nonnull input "
1017 "matrix before calling this method.");
1018 return A->getRangeMap();
1019}
1020
1021template <class LocalOrdinal, class Node>
1022void Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_BigInt, Node> >::setMatrix(const Teuchos::RCP<const row_matrix_type> &A) {
1023 if (A.getRawPtr() != getMatrix().getRawPtr()) {
1024 IsInitialized_ = false;
1025 IsComputed_ = false;
1026 A_ = A;
1027 }
1028}
1029
1030template <class LocalOrdinal, class Node>
1031Teuchos::RCP<const typename Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_BigInt, Node> >::row_matrix_type>
1032Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_BigInt, Node> >::
1033 getMatrix() const {
1034 return A_;
1035}
1036
1037template <class LocalOrdinal, class Node>
1038bool Hypre<Tpetra::RowMatrix<HYPRE_Real, LocalOrdinal, HYPRE_BigInt, Node> >::hasTransposeApply() const {
1039 return false;
1040}
1041
1042} // namespace Ifpack2
1043
1044#define IFPACK2_HYPRE_INSTANT(S, LO, GO, N) \
1045 template class Ifpack2::Hypre<Tpetra::RowMatrix<S, LO, GO, N> >;
1046
1047#endif // HAVE_HYPRE && HAVE_MPI
1048#endif // IFPACK2_HYPRE_DEF_HPP
Ifpack2
Preconditioners and smoothers for Tpetra sparse matrices.
Definition Ifpack2_AdditiveSchwarz_decl.hpp:40
Ifpack2::GMRES
@ GMRES
Uses AztecOO's GMRES.
Definition Ifpack2_CondestType.hpp:20
Generated by doxygen 1.9.8