EpetraExt Development
Loading...
Searching...
No Matches
EpetraExt_HypreIJMatrix.cpp
Go to the documentation of this file.
1//@HEADER
2// ***********************************************************************
3//
4// EpetraExt: Epetra Extended - Linear Algebra Services Package
5// Copyright (2011) Sandia Corporation
6//
7// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
8// the U.S. Government retains certain rights in this software.
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
43#ifdef HAVE_HYPRE
44
45#include "Epetra_Map.h"
46#include "Epetra_Vector.h"
47#include "Epetra_MultiVector.h"
49#include "Teuchos_RCP.hpp"
50#include "Teuchos_Assert.hpp"
51#include "Teuchos_Array.hpp"
52#include <set>
53
54//=======================================================
56 : Epetra_BasicRowMatrix(Epetra_MpiComm(hypre_IJMatrixComm(matrix))),
57 Matrix_(matrix),
58 ParMatrix_(0),
59 NumMyRows_(-1),
60 NumGlobalRows_(-1),
61 NumGlobalCols_(-1),
62 MyRowStart_(-1),
63 MyRowEnd_(-1),
64 MatType_(-1),
65 TransposeSolve_(false),
66 SolveOrPrec_(Solver)
67{
68 IsSolverSetup_ = new bool[1];
69 IsPrecondSetup_ = new bool[1];
70 IsSolverSetup_[0] = false;
71 IsPrecondSetup_[0] = false;
72 // Initialize default values for global variables
73 int ierr = 0;
74 ierr += InitializeDefaults();
75 TEUCHOS_TEST_FOR_EXCEPTION(ierr != 0, std::logic_error, "Couldn't initialize default values.");
76
77 // Create array of global row ids
78 Teuchos::Array<int> GlobalRowIDs; GlobalRowIDs.resize(NumMyRows_);
79
80 for (int i = MyRowStart_; i <= MyRowEnd_; i++) {
81 GlobalRowIDs[i-MyRowStart_] = i;
82 }
83
84 // Create array of global column ids
85 int new_value = 0; int entries = 0;
86 std::set<int> Columns;
87 int num_entries;
88 double *values;
89 int *indices;
90 for(int i = 0; i < NumMyRows_; i++){
91 ierr += HYPRE_ParCSRMatrixGetRow(ParMatrix_, i+MyRowStart_, &num_entries, &indices, &values);
92 ierr += HYPRE_ParCSRMatrixRestoreRow(ParMatrix_, i+MyRowStart_,&num_entries,&indices,&values);
93 TEUCHOS_TEST_FOR_EXCEPTION(ierr != 0, std::logic_error, "Couldn't get row of matrix.");
94 entries = num_entries;
95 for(int j = 0; j < num_entries; j++){
96 // Insert column ids from this row into set
97 new_value = indices[j];
98 Columns.insert(new_value);
99 }
100 }
101 int NumMyCols = Columns.size();
102 Teuchos::Array<int> GlobalColIDs; GlobalColIDs.resize(NumMyCols);
103
104 std::set<int>::iterator it;
105 int counter = 0;
106 for (it = Columns.begin(); it != Columns.end(); it++) {
107 // Get column ids in order
108 GlobalColIDs[counter] = *it;
109 counter = counter + 1;
110 }
111 //printf("Proc[%d] Rows from %d to %d, num = %d\n", Comm().MyPID(), MyRowStart_,MyRowEnd_, NumMyRows_);
112
113 Epetra_Map RowMap(-1, NumMyRows_, &GlobalRowIDs[0], 0, Comm());
114 Epetra_Map ColMap(-1, NumMyCols, &GlobalColIDs[0], 0, Comm());
115
116 //Need to call SetMaps()
117 SetMaps(RowMap, ColMap);
118
119 // Get an MPI_Comm to create vectors.
120 // The vectors will be reused in Multiply(), so that they aren't recreated every time.
121 MPI_Comm comm;
122 ierr += HYPRE_ParCSRMatrixGetComm(ParMatrix_, &comm);
123 TEUCHOS_TEST_FOR_EXCEPTION(ierr != 0, std::logic_error, "Couldn't get communicator from Hypre Matrix.");
124
125 ierr += HYPRE_IJVectorCreate(comm, MyRowStart_, MyRowEnd_, &X_hypre);
126 ierr += HYPRE_IJVectorSetObjectType(X_hypre, HYPRE_PARCSR);
127 ierr += HYPRE_IJVectorInitialize(X_hypre);
128 ierr += HYPRE_IJVectorAssemble(X_hypre);
129 ierr += HYPRE_IJVectorGetObject(X_hypre, (void**) &par_x);
130 TEUCHOS_TEST_FOR_EXCEPTION(ierr != 0, std::logic_error, "Couldn't create Hypre X vector.");
131
132 ierr += HYPRE_IJVectorCreate(comm, MyRowStart_, MyRowEnd_, &Y_hypre);
133 ierr += HYPRE_IJVectorSetObjectType(Y_hypre, HYPRE_PARCSR);
134 ierr += HYPRE_IJVectorInitialize(Y_hypre);
135 ierr += HYPRE_IJVectorAssemble(Y_hypre);
136 ierr += HYPRE_IJVectorGetObject(Y_hypre, (void**) &par_y);
137 TEUCHOS_TEST_FOR_EXCEPTION(ierr != 0, std::logic_error, "Couldn't create Hypre Y vector.");
138
139 x_vec = (hypre_ParVector *) hypre_IJVectorObject(((hypre_IJVector *) X_hypre));
140 x_local = hypre_ParVectorLocalVector(x_vec);
141
142 y_vec = (hypre_ParVector *) hypre_IJVectorObject(((hypre_IJVector *) Y_hypre));
143 y_local = hypre_ParVectorLocalVector(y_vec);
144
146 SolverDestroyPtr_ = &HYPRE_ParCSRPCGDestroy;
147 SolverSetupPtr_ = &HYPRE_ParCSRPCGSetup;
148 SolverSolvePtr_ = &HYPRE_ParCSRPCGSolve;
149 SolverPrecondPtr_ = &HYPRE_ParCSRPCGSetPrecond;
150 CreateSolver();
151
153 PrecondDestroyPtr_ = &HYPRE_EuclidDestroy;
154 PrecondSetupPtr_ = &HYPRE_EuclidSetup;
155 PrecondSolvePtr_ = &HYPRE_EuclidSolve;
156 CreatePrecond();
157 ComputeNumericConstants();
158 ComputeStructureConstants();
159} //EpetraExt_HYPREIJMatrix(Hypre_IJMatrix) Constructor
160
161//=======================================================
163 int ierr = 0;
164 ierr += HYPRE_IJVectorDestroy(X_hypre);
165 TEUCHOS_TEST_FOR_EXCEPTION(ierr != 0, std::logic_error, "Couldn't destroy the X Vector.");
166 ierr += HYPRE_IJVectorDestroy(Y_hypre);
167 TEUCHOS_TEST_FOR_EXCEPTION(ierr != 0, std::logic_error, "Couldn't destroy the Y Vector.");
168
169 /* Destroy solver and preconditioner */
170 if(IsSolverSetup_[0] == true){
171 ierr += SolverDestroyPtr_(Solver_);
172 TEUCHOS_TEST_FOR_EXCEPTION(ierr != 0, std::logic_error, "Couldn't destroy the Solver.");
173 }
174 if(IsPrecondSetup_[0] == true){
176 TEUCHOS_TEST_FOR_EXCEPTION(ierr != 0, std::logic_error, "Couldn't destroy the Preconditioner.");
177 }
178 delete[] IsSolverSetup_;
179 delete[] IsPrecondSetup_;
180} //EpetraExt_HypreIJMatrix destructor
181
182//=======================================================
183int EpetraExt_HypreIJMatrix::ExtractMyRowCopy(int Row, int Length, int & NumEntries,
184 double * Values, int * Indices) const
185{
186 // Get values and indices of ith row of matrix
187 int *indices;
188 double *values;
189 int num_entries;
190 EPETRA_CHK_ERR(HYPRE_ParCSRMatrixGetRow(ParMatrix_, Row+RowMatrixRowMap().MinMyGID(), &num_entries, &indices, &values));
191 EPETRA_CHK_ERR(HYPRE_ParCSRMatrixRestoreRow(ParMatrix_, Row+RowMatrixRowMap().MinMyGID(), &num_entries, &indices, &values));
192
193 NumEntries = num_entries;
194
195 if(Length < NumEntries){
196 printf("The arrays passed in are not large enough. Allocate more space.\n");
197 return -2;
198 }
199
200 for(int i = 0; i < NumEntries; i++){
201 Values[i] = values[i];
202 Indices[i] = RowMatrixColMap().LID(indices[i]);
203 }
204
205 return 0;
206} //ExtractMyRowCopy()
207
208//=======================================================
209int EpetraExt_HypreIJMatrix::NumMyRowEntries(int Row, int & NumEntries) const
210{
211 int my_row[1];
212 my_row[0] = Row+RowMatrixRowMap().MinMyGID();
213 int nentries[1];
214 EPETRA_CHK_ERR(HYPRE_IJMatrixGetRowCounts(Matrix_, 1, my_row, nentries));
215 NumEntries = nentries[0];
216 return 0;
217} //NumMyRowEntries()
218
219//=======================================================
220int EpetraExt_HypreIJMatrix::ExtractMyEntryView(int CurEntry, double *&Value, int &RowIndex, int &ColIndex)
221{/*
222 This gives a copy, not a view of the values, so is not implemented.
223 if(CurEntry >= NumMyNonzeros() || CurEntry < 0){
224 return -1;
225 }
226 int counter = 0;
227 for(int i = 0; i < NumMyRows(); i++){
228 int entries;
229 NumMyRowEntries(i, entries);
230 counter += entries;
231 if(counter > CurEntry) {
232 counter -= entries;
233 RowIndex = i;
234 break;
235 }
236 }
237 int entries;
238 NumMyRowEntries(RowIndex, entries);
239 int *indices;
240 double *values;
241 int num_entries;
242 HYPRE_ParCSRMatrixGetRow(ParMatrix_, RowIndex+MyRowStart_, &num_entries, &indices, &values);
243 HYPRE_ParCSRMatrixRestoreRow(ParMatrix_, RowIndex+MyRowStart_, &num_entries, &indices, &values);
244 ColIndex = RowMatrixColMap().LID(indices[CurEntry-counter]);
245 Value = &values[CurEntry-counter];
246 return 0;*/return -1;
247} //ExtractMyEntryView() - not implemented
248
249 //=======================================================
250 int EpetraExt_HypreIJMatrix::ExtractMyEntryView(int CurEntry, double const * & Value, int & RowIndex, int & ColIndex) const
251{/*
252 if(CurEntry >= NumMyNonzeros() || CurEntry < 0){
253 return -1;
254 }
255 int counter = 0;
256 for(int i = 0; i < NumMyRows(); i++){
257 int entries;
258 NumMyRowEntries(i, entries);
259 counter += entries;
260 if(counter > CurEntry) {
261 counter -= entries;
262 RowIndex = i;
263 break;
264 }
265 }
266 int entries;
267 NumMyRowEntries(RowIndex, entries);
268 int *indices;
269 double *values;
270 int num_entries;
271 HYPRE_ParCSRMatrixGetRow(ParMatrix_, RowIndex+MyRowStart_, &num_entries, &indices, &values);
272 HYPRE_ParCSRMatrixRestoreRow(ParMatrix_, RowIndex+MyRowStart_, &num_entries, &indices, &values);
273 ColIndex = RowMatrixColMap().LID(indices[CurEntry-counter]);
274 Value = &values[CurEntry-counter];
275 return 0;*/ return -1;
276} //ExtractMyEntryView() - const version, not implemented
277
278//=======================================================
280 const Epetra_MultiVector& X,
281 Epetra_MultiVector& Y) const
282{
283
284 //printf("Proc[%d], Row start: %d, Row End: %d\n", Comm().MyPID(), MyRowStart_, MyRowEnd_);
285 bool SameVectors = false;
286 int NumVectors = X.NumVectors();
287 if (NumVectors != Y.NumVectors()) return -1; // X and Y must have same number of vectors
288 if(X.Pointers() == Y.Pointers()){
289 SameVectors = true;
290 }
291 for(int VecNum = 0; VecNum < NumVectors; VecNum++) {
292 //Get values for current vector in multivector.
293 double * x_values;
294 double * y_values;
295 EPETRA_CHK_ERR((*X(VecNum)).ExtractView(&x_values));
296 double *x_temp = x_local->data;
297 double *y_temp = y_local->data;
298 if(!SameVectors){
299 EPETRA_CHK_ERR((*Y(VecNum)).ExtractView(&y_values));
300 } else {
301 y_values = new double[X.MyLength()];
302 }
303 y_local->data = y_values;
304 EPETRA_CHK_ERR(HYPRE_ParVectorSetConstantValues(par_y,0.0));
305 // Temporarily make a pointer to data in Hypre for end
306 // Replace data in Hypre vectors with epetra values
307 x_local->data = x_values;
308 // Do actual computation.
309 if(TransA) {
310 // Use transpose of A in multiply
311 EPETRA_CHK_ERR(HYPRE_ParCSRMatrixMatvecT(1.0, ParMatrix_, par_x, 1.0, par_y));
312 } else {
313 EPETRA_CHK_ERR(HYPRE_ParCSRMatrixMatvec(1.0, ParMatrix_, par_x, 1.0, par_y));
314 }
315 if(SameVectors){
316 int NumEntries = Y.MyLength();
317 std::vector<double> new_values; new_values.resize(NumEntries);
318 std::vector<int> new_indices; new_indices.resize(NumEntries);
319 for(int i = 0; i < NumEntries; i++){
320 new_values[i] = y_values[i];
321 new_indices[i] = i;
322 }
323 EPETRA_CHK_ERR((*Y(VecNum)).ReplaceMyValues(NumEntries, &new_values[0], &new_indices[0]));
324 delete[] y_values;
325 }
326 x_local->data = x_temp;
327 y_local->data = y_temp;
328 }
329 double flops = (double) NumVectors * (double) NumGlobalNonzeros();
330 UpdateFlops(flops);
331 return 0;
332} //Multiply()
333
334//=======================================================
335int EpetraExt_HypreIJMatrix::SetParameter(Hypre_Chooser chooser, int (*pt2Func)(HYPRE_Solver, int), int parameter){
336 if(chooser == Preconditioner){
337 EPETRA_CHK_ERR(pt2Func(Preconditioner_, parameter));
338 } else {
339 EPETRA_CHK_ERR(pt2Func(Solver_, parameter));
340 }
341 return 0;
342} //SetParameter() - int function pointer
343
344//=======================================================
345int EpetraExt_HypreIJMatrix::SetParameter(Hypre_Chooser chooser, int (*pt2Func)(HYPRE_Solver, double), double parameter){
346 if(chooser == Preconditioner){
347 EPETRA_CHK_ERR(pt2Func(Preconditioner_, parameter));
348 } else {
349 EPETRA_CHK_ERR(pt2Func(Solver_, parameter));
350 }
351 return 0;
352} //SetParamater() - double function pointer
353
354//=======================================================
355int EpetraExt_HypreIJMatrix::SetParameter(Hypre_Chooser chooser, int (*pt2Func)(HYPRE_Solver, double, int), double parameter1, int parameter2){
356 if(chooser == Preconditioner){
357 EPETRA_CHK_ERR(pt2Func(Preconditioner_, parameter1, parameter2));
358 } else {
359 EPETRA_CHK_ERR(pt2Func(Solver_, parameter1, parameter2));
360 }
361 return 0;
362} //SetParameter() - double,int function pointer
363
364//=======================================================
365int EpetraExt_HypreIJMatrix::SetParameter(Hypre_Chooser chooser, int (*pt2Func)(HYPRE_Solver, int, int), int parameter1, int parameter2){
366 if(chooser == Preconditioner){
367 EPETRA_CHK_ERR(pt2Func(Preconditioner_, parameter1, parameter2));
368 } else {
369 EPETRA_CHK_ERR(pt2Func(Solver_, parameter1, parameter2));
370 }
371 return 0;
372} //SetParameter() - int,int function pointer
373
374//=======================================================
375int EpetraExt_HypreIJMatrix::SetParameter(Hypre_Chooser chooser, int (*pt2Func)(HYPRE_Solver, double*), double* parameter){
376 if(chooser == Preconditioner){
377 EPETRA_CHK_ERR(pt2Func(Preconditioner_, parameter));
378 } else {
379 EPETRA_CHK_ERR(pt2Func(Solver_, parameter));
380 }
381 return 0;
382} //SetParameter() - double* function pointer
383
384//=======================================================
385int EpetraExt_HypreIJMatrix::SetParameter(Hypre_Chooser chooser, int (*pt2Func)(HYPRE_Solver, int*), int* parameter){
386 if(chooser == Preconditioner){
387 EPETRA_CHK_ERR(pt2Func(Preconditioner_, parameter));
388 } else {
389 EPETRA_CHK_ERR(pt2Func(Solver_, parameter));
390 }
391 return 0;
392} //SetParameter() - int* function pointer
393
394//=======================================================
395int EpetraExt_HypreIJMatrix::SetParameter(Hypre_Chooser chooser, Hypre_Solver solver, bool transpose){
396 if(chooser == Solver){
397 if(transpose && solver != BoomerAMG){
398 EPETRA_CHK_ERR(-1);
399 }
400 switch(solver) {
401 case BoomerAMG:
402 if(IsSolverSetup_[0]){
404 IsSolverSetup_[0] = false;
405 }
407 SolverDestroyPtr_ = &HYPRE_BoomerAMGDestroy;
408 SolverSetupPtr_ = &HYPRE_BoomerAMGSetup;
409 SolverPrecondPtr_ = NULL;
410 if(transpose){
411 TransposeSolve_ = true;
412 SolverSolvePtr_ = &HYPRE_BoomerAMGSolveT;
413 } else {
414 SolverSolvePtr_ = &HYPRE_BoomerAMGSolve;
415 }
416 break;
417 case AMS:
418 if(IsSolverSetup_[0]){
420 IsSolverSetup_[0] = false;
421 }
423 SolverDestroyPtr_ = &HYPRE_AMSDestroy;
424 SolverSetupPtr_ = &HYPRE_AMSSetup;
425 SolverSolvePtr_ = &HYPRE_AMSSolve;
426 SolverPrecondPtr_ = NULL;
427 break;
428 case Hybrid:
429 if(IsSolverSetup_[0]){
431 IsSolverSetup_[0] = false;
432 }
434 SolverDestroyPtr_ = &HYPRE_ParCSRHybridDestroy;
435 SolverSetupPtr_ = &HYPRE_ParCSRHybridSetup;
436 SolverSolvePtr_ = &HYPRE_ParCSRHybridSolve;
437 SolverPrecondPtr_ = &HYPRE_ParCSRHybridSetPrecond;
438 break;
439 case PCG:
440 if(IsSolverSetup_[0]){
442 IsSolverSetup_[0] = false;
443 }
445 SolverDestroyPtr_ = &HYPRE_ParCSRPCGDestroy;
446 SolverSetupPtr_ = &HYPRE_ParCSRPCGSetup;
447 SolverSolvePtr_ = &HYPRE_ParCSRPCGSolve;
448 SolverPrecondPtr_ = &HYPRE_ParCSRPCGSetPrecond;
449 break;
450 case GMRES:
451 if(IsSolverSetup_[0]){
453 IsSolverSetup_[0] = false;
454 }
456 SolverDestroyPtr_ = &HYPRE_ParCSRGMRESDestroy;
457 SolverSetupPtr_ = &HYPRE_ParCSRGMRESSetup;
458 SolverSolvePtr_ = &HYPRE_ParCSRGMRESSolve;
459 SolverPrecondPtr_ = &HYPRE_ParCSRGMRESSetPrecond;
460 break;
461 case FlexGMRES:
462 if(IsSolverSetup_[0]){
464 IsSolverSetup_[0] = false;
465 }
467 SolverDestroyPtr_ = &HYPRE_ParCSRFlexGMRESDestroy;
468 SolverSetupPtr_ = &HYPRE_ParCSRFlexGMRESSetup;
469 SolverSolvePtr_ = &HYPRE_ParCSRFlexGMRESSolve;
470 SolverPrecondPtr_ = &HYPRE_ParCSRFlexGMRESSetPrecond;
471 break;
472 case LGMRES:
473 if(IsSolverSetup_[0]){
475 IsSolverSetup_[0] = false;
476 }
478 SolverDestroyPtr_ = &HYPRE_ParCSRLGMRESDestroy;
479 SolverSetupPtr_ = &HYPRE_ParCSRLGMRESSetup;
480 SolverSolvePtr_ = &HYPRE_ParCSRLGMRESSolve;
481 SolverPrecondPtr_ = &HYPRE_ParCSRLGMRESSetPrecond;
482 break;
483 case BiCGSTAB:
484 if(IsSolverSetup_[0]){
486 IsSolverSetup_[0] = false;
487 }
489 SolverDestroyPtr_ = &HYPRE_ParCSRBiCGSTABDestroy;
490 SolverSetupPtr_ = &HYPRE_ParCSRBiCGSTABSetup;
491 SolverSolvePtr_ = &HYPRE_ParCSRBiCGSTABSolve;
492 SolverPrecondPtr_ = &HYPRE_ParCSRBiCGSTABSetPrecond;
493 break;
494 default:
495 EPETRA_CHK_ERR(-1);
496 }
497 CreateSolver();
498 } else {
499 // Preconditioner
500 switch(solver) {
501 case BoomerAMG:
502 if(IsPrecondSetup_[0]){
504 IsPrecondSetup_[0] = false;
505 }
507 PrecondDestroyPtr_ = &HYPRE_BoomerAMGDestroy;
508 PrecondSetupPtr_ = &HYPRE_BoomerAMGSetup;
509 PrecondSolvePtr_ = &HYPRE_BoomerAMGSolve;
510 break;
511 case ParaSails:
512 if(IsPrecondSetup_[0]){
514 IsPrecondSetup_[0] = false;
515 }
517 PrecondDestroyPtr_ = &HYPRE_ParaSailsDestroy;
518 PrecondSetupPtr_ = &HYPRE_ParaSailsSetup;
519 PrecondSolvePtr_ = &HYPRE_ParaSailsSolve;
520 break;
521 case Euclid:
522 if(IsPrecondSetup_[0]){
524 IsPrecondSetup_[0] = false;
525 }
527 PrecondDestroyPtr_ = &HYPRE_EuclidDestroy;
528 PrecondSetupPtr_ = &HYPRE_EuclidSetup;
529 PrecondSolvePtr_ = &HYPRE_EuclidSolve;
530 break;
531 case AMS:
532 if(IsPrecondSetup_[0]){
534 IsPrecondSetup_[0] = false;
535 }
537 PrecondDestroyPtr_ = &HYPRE_AMSDestroy;
538 PrecondSetupPtr_ = &HYPRE_AMSSetup;
539 PrecondSolvePtr_ = &HYPRE_AMSSolve;
540 break;
541 default:
542 EPETRA_CHK_ERR(-1);
543 }
545
546 }
547 return 0;
548} //SetParameter() - Choose solver or preconditioner type
549
550//=======================================================
552 MPI_Comm comm;
553 HYPRE_ParCSRMatrixGetComm(ParMatrix_, &comm);
554 return (this->*SolverCreatePtr_)(comm, &Solver_);
555} //CreateSolver()
556
557//=======================================================
559 MPI_Comm comm;
560 HYPRE_ParCSRMatrixGetComm(ParMatrix_, &comm);
561 return (this->*PrecondCreatePtr_)(comm, &Preconditioner_);
562} //CreatePrecond()
563
564//=======================================================
565int EpetraExt_HypreIJMatrix::SetParameter(bool UsePreconditioner){
566 if(UsePreconditioner == false){
567 return 0;
568 } else {
569 if(SolverPrecondPtr_ == NULL){
570 return -1;
571 } else {
573 return 0;
574 }
575 }
576} //SetParameter() - Set the precondioner pointer for solver
577
578//=======================================================
580 SolveOrPrec_ = answer;
581 return 0;
582} //SetParameter() - Choose to solve or precondition
583
584//=======================================================
587 IsSolverSetup_[0] = true;
588 return 0;
589} //SetupSolver()
590
591//=======================================================
594 IsPrecondSetup_[0] = true;
595 return 0;
596} //SetupPrecond()
597
598//=======================================================
599int EpetraExt_HypreIJMatrix::Solve(bool Upper, bool transpose, bool UnitDiagonal, const Epetra_MultiVector & X, Epetra_MultiVector & Y) const {
600 bool SameVectors = false;
601 int NumVectors = X.NumVectors();
602 if (NumVectors != Y.NumVectors()) return -1; // X and Y must have same number of vectors
603 if(X.Pointers() == Y.Pointers()){
604 SameVectors = true;
605 }
606 if(SolveOrPrec_ == Solver){
607 if(IsSolverSetup_[0] == false){
608 SetupSolver();
609 }
610 } else {
611 if(IsPrecondSetup_[0] == false){
612 SetupPrecond();
613 }
614 }
615 for(int VecNum = 0; VecNum < NumVectors; VecNum++) {
616 //Get values for current vector in multivector.
617 double * x_values;
618 EPETRA_CHK_ERR((*X(VecNum)).ExtractView(&x_values));
619 double * y_values;
620 if(!SameVectors){
621 EPETRA_CHK_ERR((*Y(VecNum)).ExtractView(&y_values));
622 } else {
623 y_values = new double[X.MyLength()];
624 }
625 // Temporarily make a pointer to data in Hypre for end
626 double *x_temp = x_local->data;
627 // Replace data in Hypre vectors with epetra values
628 x_local->data = x_values;
629 double *y_temp = y_local->data;
630 y_local->data = y_values;
631
632 EPETRA_CHK_ERR(HYPRE_ParVectorSetConstantValues(par_y, 0.0));
633 if(transpose && !TransposeSolve_){
634 // User requested a transpose solve, but the solver selected doesn't provide one
635 EPETRA_CHK_ERR(-1);
636 }
637 if(SolveOrPrec_ == Solver){
638 // Use the solver methods
639 EPETRA_CHK_ERR(SolverSolvePtr_(Solver_, ParMatrix_, par_x, par_y));
640 } else {
641 // Apply the preconditioner
643 }
644 if(SameVectors){
645 int NumEntries = Y.MyLength();
646 std::vector<double> new_values; new_values.resize(NumEntries);
647 std::vector<int> new_indices; new_indices.resize(NumEntries);
648 for(int i = 0; i < NumEntries; i++){
649 new_values[i] = y_values[i];
650 new_indices[i] = i;
651 }
652 EPETRA_CHK_ERR((*Y(VecNum)).ReplaceMyValues(NumEntries, &new_values[0], &new_indices[0]));
653 delete[] y_values;
654 }
655 x_local->data = x_temp;
656 y_local->data = y_temp;
657 }
658 double flops = (double) NumVectors * (double) NumGlobalNonzeros();
659 UpdateFlops(flops);
660 return 0;
661} //Solve()
662
663//=======================================================
665 for(int i = 0; i < NumMyRows_; i++){
666 //Vector-scalar mult on ith row
667 int num_entries;
668 int *indices;
669 double *values;
670 EPETRA_CHK_ERR(HYPRE_ParCSRMatrixGetRow(ParMatrix_,i+MyRowStart_, &num_entries, &indices, &values));
671 EPETRA_CHK_ERR(HYPRE_ParCSRMatrixRestoreRow(ParMatrix_,i+MyRowStart_, &num_entries, &indices, &values));
672 Teuchos::Array<double> new_values; new_values.resize(num_entries);
673 Teuchos::Array<int> new_indices; new_indices.resize(num_entries);
674 for(int j = 0; j < num_entries; j++){
675 // Scale this row with the appropriate values from the vector
676 new_values[j] = X[i]*values[j];
677 new_indices[j] = indices[j];
678 }
679 int rows[1];
680 rows[0] = i+MyRowStart_;
681 EPETRA_CHK_ERR(HYPRE_IJMatrixSetValues(Matrix_, 1, &num_entries, rows, &new_indices[0], &new_values[0]));
682 // Finally set values of the Matrix for this row
683 }
684 HaveNumericConstants_ = false;
686 return 0;
687} //LeftScale()
688
689//=======================================================
691 // First we need to import off-processor values of the vector
693 Epetra_Vector Import_Vector(RowMatrixColMap(), true);
694 EPETRA_CHK_ERR(Import_Vector.Import(X, Importer, Insert, 0));
695
696 for(int i = 0; i < NumMyRows_; i++){
697 //Vector-scalar mult on ith column
698 int num_entries;
699 double *values;
700 int *indices;
701 // Get values and indices of ith row of matrix
702 EPETRA_CHK_ERR(HYPRE_ParCSRMatrixGetRow(ParMatrix_,i+MyRowStart_, &num_entries, &indices, &values));
703 EPETRA_CHK_ERR(HYPRE_ParCSRMatrixRestoreRow(ParMatrix_,i+MyRowStart_,&num_entries,&indices,&values));
704 Teuchos::Array<int> new_indices; new_indices.resize(num_entries);
705 Teuchos::Array<double> new_values; new_values.resize(num_entries);
706 for(int j = 0; j < num_entries; j++){
707 // Multiply column j with jth element
708 int index = RowMatrixColMap().LID(indices[j]);
709 TEUCHOS_TEST_FOR_EXCEPTION(index < 0, std::logic_error, "Index is negtive.");
710 new_values[j] = values[j] * Import_Vector[index];
711 new_indices[j] = indices[j];
712 }
713 // Finally set values of the Matrix for this row
714 int rows[1];
715 rows[0] = i+MyRowStart_;
716 EPETRA_CHK_ERR(HYPRE_IJMatrixSetValues(Matrix_, 1, &num_entries, rows, &new_indices[0], &new_values[0]));
717 }
718
719 HaveNumericConstants_ = false;
721 return 0;
722} //RightScale()
723
724//=======================================================
726 int my_type;
727 // Get type of Hypre IJ Matrix
728 EPETRA_CHK_ERR(HYPRE_IJMatrixGetObjectType(Matrix_, &my_type));
729 MatType_ = my_type;
730 // Currently only ParCSR is supported
731 TEUCHOS_TEST_FOR_EXCEPTION(MatType_ != HYPRE_PARCSR, std::logic_error, "Object is not type PARCSR");
732
733 // Get the actual ParCSR object from the IJ matrix
734 EPETRA_CHK_ERR(HYPRE_IJMatrixGetObject(Matrix_, (void**) &ParMatrix_));
735
736 int numRows, numCols;
737
738 // Get dimensions of the matrix and store as global variables
739 EPETRA_CHK_ERR(HYPRE_ParCSRMatrixGetDims(ParMatrix_, &numRows, &numCols));
740
741 NumGlobalRows_ = numRows;
742 NumGlobalCols_ = numCols;
743
744 // Get the local dimensions of the matrix
745 int ColStart, ColEnd;
746 EPETRA_CHK_ERR(HYPRE_ParCSRMatrixGetLocalRange(ParMatrix_, &MyRowStart_, &MyRowEnd_, &ColStart, &ColEnd));
747
748 // Determine number of local rows
750
751 return 0;
752} //InitializeDefaults()
753
754#endif /*HAVE_HYPRE*/
Hypre_Chooser
Enumerated type to choose to solve or precondition.
@ Preconditioner
Hypre_Solver
Enumerated type for Hypre solvers.
int Multiply(bool TransA, const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
Returns the result of a EpetraExt_HypreIJMatrix multiplied by a Epetra_MultiVector X in Y.
int(* SolverPrecondPtr_)(HYPRE_Solver, HYPRE_PtrToParSolverFcn, HYPRE_PtrToParSolverFcn, HYPRE_Solver)
int Hypre_ParCSRBiCGSTABCreate(MPI_Comm comm, HYPRE_Solver *solver)
BiCGSTAB Create passing function.
bool * IsSolverSetup_
Flag to know if solver needs to be destoyed.
int(* SolverSolvePtr_)(HYPRE_Solver, HYPRE_ParCSRMatrix, HYPRE_ParVector, HYPRE_ParVector)
int(* SolverSetupPtr_)(HYPRE_Solver, HYPRE_ParCSRMatrix, HYPRE_ParVector, HYPRE_ParVector)
int MyRowStart_
First row on local processor.
int InitializeDefaults()
Set global variables to default values.
int CreatePrecond()
Create the preconditioner with selected type.
int Hypre_EuclidCreate(MPI_Comm comm, HYPRE_Solver *solver)
Euclid Create passing function.
bool TransposeSolve_
Do a transpose solve, only BoomerAMG.
int NumGlobalCols_
Number of columns across all processors.
int(* PrecondSolvePtr_)(HYPRE_Solver, HYPRE_ParCSRMatrix, HYPRE_ParVector, HYPRE_ParVector)
int(EpetraExt_HypreIJMatrix::* PrecondCreatePtr_)(MPI_Comm, HYPRE_Solver *)
int ExtractMyEntryView(int CurEntry, double *&Value, int &RowIndex, int &ColIndex)
Returns a reference to the ith entry in the matrix, along with its row and column index.
int CreateSolver()
Create the solver with selected type.
int Hypre_BoomerAMGCreate(MPI_Comm comm, HYPRE_Solver *solver)
AMG Create passing function.
int Hypre_ParCSRHybridCreate(MPI_Comm comm, HYPRE_Solver *solver)
Hybrid Create passing function.
int LeftScale(const Epetra_Vector &X)
Scales the EpetraExt_HypreIJMatrix on the left with a Epetra_Vector x.
int SetParameter(Hypre_Chooser chooser, int(*pt2Func)(HYPRE_Solver, int), int parameter)
Set a parameter that takes a single int.
int NumGlobalRows_
Number of rows across all processors.
int Hypre_AMSCreate(MPI_Comm comm, HYPRE_Solver *solver)
AMS Create passing function.
int Hypre_ParCSRFlexGMRESCreate(MPI_Comm comm, HYPRE_Solver *solver)
FlexGMRES Create passing function.
int(* SolverDestroyPtr_)(HYPRE_Solver)
EpetraExt_HypreIJMatrix(HYPRE_IJMatrix matrix)
Epetra_HypreIJMatrix constructor.
int(EpetraExt_HypreIJMatrix::* SolverCreatePtr_)(MPI_Comm, HYPRE_Solver *)
int(* PrecondDestroyPtr_)(HYPRE_Solver)
Hypre_Chooser SolveOrPrec_
Choose to solve or apply preconditioner.
int Solve(bool Upper, bool Trans, bool UnitDiagonal, const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
Returns the result of a EpetraExt_HypreIJMatrix solving a Epetra_MultiVector X in Y.
virtual ~EpetraExt_HypreIJMatrix()
EpetraExt_HypreIJMatrix Destructor.
int NumMyRows_
Number of rows on local processor.
int Hypre_ParCSRPCGCreate(MPI_Comm comm, HYPRE_Solver *solver)
PCG Create passing function.
int Hypre_ParCSRLGMRESCreate(MPI_Comm comm, HYPRE_Solver *solver)
LGMRES Create passing function.
int NumMyRowEntries(int MyRow, int &NumEntries) const
Return the current number of values stored for the specified local row.
int MyRowEnd_
Last row on local processor.
int Hypre_ParaSailsCreate(MPI_Comm comm, HYPRE_Solver *solver)
ParaSails Create passing function.
int Hypre_ParCSRGMRESCreate(MPI_Comm comm, HYPRE_Solver *solver)
GMRES Create passing function.
bool * IsPrecondSetup_
Flag to know if preconditioner needs to be destroyed.
int ExtractMyRowCopy(int MyRow, int Length, int &NumEntries, double *Values, int *Indices) const
Returns a copy of the specified local row in user-provided arrays.
int(* PrecondSetupPtr_)(HYPRE_Solver, HYPRE_ParCSRMatrix, HYPRE_ParVector, HYPRE_ParVector)
int MatType_
Hypre matrix type (parCSR).
int RightScale(const Epetra_Vector &X)
Scales the EpetraExt_HypreIJMatrix on the right with a Epetra_Vector x.
virtual const Epetra_Map & RowMatrixColMap() const
virtual const Epetra_Import * Importer() const
virtual const Epetra_Map & RowMatrixRowMap() const
virtual int NumGlobalNonzeros() const
int MinMyGID() const
int LID(int GID) const
void UpdateFlops(int Flops_in) const
int NumVectors() const
int MyLength() const
double ** Pointers() const