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