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