EpetraExt  Development
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 
42 #include "EpetraExt_ConfigDefs.h"
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 
152  PrecondCreatePtr_ = &EpetraExt_HypreIJMatrix::Hypre_EuclidCreate;
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){
175  ierr += PrecondDestroyPtr_(Preconditioner_);
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 //=======================================================
183 int 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 //=======================================================
209 int 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 //=======================================================
220 int 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 //=======================================================
279 int EpetraExt_HypreIJMatrix::Multiply(bool TransA,
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 //=======================================================
335 int 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 //=======================================================
345 int 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 //=======================================================
355 int 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 //=======================================================
365 int 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 //=======================================================
375 int 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 //=======================================================
385 int 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 //=======================================================
395 int 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]){
403  SolverDestroyPtr_(Solver_);
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]){
419  SolverDestroyPtr_(Solver_);
420  IsSolverSetup_[0] = false;
421  }
422  SolverCreatePtr_ = &EpetraExt_HypreIJMatrix::Hypre_AMSCreate;
423  SolverDestroyPtr_ = &HYPRE_AMSDestroy;
424  SolverSetupPtr_ = &HYPRE_AMSSetup;
425  SolverSolvePtr_ = &HYPRE_AMSSolve;
426  SolverPrecondPtr_ = NULL;
427  break;
428  case Hybrid:
429  if(IsSolverSetup_[0]){
430  SolverDestroyPtr_(Solver_);
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]){
441  SolverDestroyPtr_(Solver_);
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]){
452  SolverDestroyPtr_(Solver_);
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]){
463  SolverDestroyPtr_(Solver_);
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]){
474  SolverDestroyPtr_(Solver_);
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]){
485  SolverDestroyPtr_(Solver_);
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]){
503  PrecondDestroyPtr_(Preconditioner_);
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]){
513  PrecondDestroyPtr_(Preconditioner_);
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]){
523  PrecondDestroyPtr_(Preconditioner_);
524  IsPrecondSetup_[0] = false;
525  }
526  PrecondCreatePtr_ = &EpetraExt_HypreIJMatrix::Hypre_EuclidCreate;
527  PrecondDestroyPtr_ = &HYPRE_EuclidDestroy;
528  PrecondSetupPtr_ = &HYPRE_EuclidSetup;
529  PrecondSolvePtr_ = &HYPRE_EuclidSolve;
530  break;
531  case AMS:
532  if(IsPrecondSetup_[0]){
533  PrecondDestroyPtr_(Preconditioner_);
534  IsPrecondSetup_[0] = false;
535  }
536  PrecondCreatePtr_ = &EpetraExt_HypreIJMatrix::Hypre_AMSCreate;
537  PrecondDestroyPtr_ = &HYPRE_AMSDestroy;
538  PrecondSetupPtr_ = &HYPRE_AMSSetup;
539  PrecondSolvePtr_ = &HYPRE_AMSSolve;
540  break;
541  default:
542  EPETRA_CHK_ERR(-1);
543  }
544  CreatePrecond();
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 //=======================================================
565 int EpetraExt_HypreIJMatrix::SetParameter(bool UsePreconditioner){
566  if(UsePreconditioner == false){
567  return 0;
568  } else {
569  if(SolverPrecondPtr_ == NULL){
570  return -1;
571  } else {
572  SolverPrecondPtr_(Solver_, PrecondSolvePtr_, PrecondSetupPtr_, Preconditioner_);
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 //=======================================================
586  SolverSetupPtr_(Solver_, ParMatrix_, par_x, par_y);
587  IsSolverSetup_[0] = true;
588  return 0;
589 } //SetupSolver()
590 
591 //=======================================================
593  PrecondSetupPtr_(Preconditioner_, ParMatrix_, par_x, par_y);
594  IsPrecondSetup_[0] = true;
595  return 0;
596 } //SetupPrecond()
597 
598 //=======================================================
599 int 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
642  EPETRA_CHK_ERR(PrecondSolvePtr_(Preconditioner_, ParMatrix_, par_x, par_y));
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;
685  UpdateFlops(NumGlobalNonzeros());
686  return 0;
687 } //LeftScale()
688 
689 //=======================================================
691  // First we need to import off-processor values of the vector
692  Epetra_Import Importer(RowMatrixColMap(), RowMatrixRowMap());
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;
720  UpdateFlops(NumGlobalNonzeros());
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
749  NumMyRows_ = MyRowEnd_ - MyRowStart_+1;
750 
751  return 0;
752 } //InitializeDefaults()
753 
754 #endif /*HAVE_HYPRE*/
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 Hypre_ParCSRLGMRESCreate(MPI_Comm comm, HYPRE_Solver *solver)
LGMRES Create passing function.
virtual ~EpetraExt_HypreIJMatrix()
EpetraExt_HypreIJMatrix Destructor.
int LeftScale(const Epetra_Vector &X)
Scales the EpetraExt_HypreIJMatrix on the left with a Epetra_Vector x.
int InitializeDefaults()
Set global variables to default values.
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...
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 Hypre_ParCSRPCGCreate(MPI_Comm comm, HYPRE_Solver *solver)
PCG Create passing function.
EpetraExt_HypreIJMatrix(HYPRE_IJMatrix matrix)
Epetra_HypreIJMatrix constructor.
int NumMyRowEntries(int MyRow, int &NumEntries) const
Return the current number of values stored for the specified local row.
Hypre_Chooser
Enumerated type to choose to solve or precondition.
int Hypre_ParCSRGMRESCreate(MPI_Comm comm, HYPRE_Solver *solver)
GMRES Create passing function.
int RightScale(const Epetra_Vector &X)
Scales the EpetraExt_HypreIJMatrix on the right with a Epetra_Vector x.
int Hypre_ParCSRFlexGMRESCreate(MPI_Comm comm, HYPRE_Solver *solver)
FlexGMRES Create passing function.
int SetupPrecond() const
int CreatePrecond()
Create the preconditioner with selected type.
int Hypre_ParaSailsCreate(MPI_Comm comm, HYPRE_Solver *solver)
ParaSails Create passing function.
int Hypre_AMSCreate(MPI_Comm comm, HYPRE_Solver *solver)
AMS Create passing function.
int Hypre_BoomerAMGCreate(MPI_Comm comm, HYPRE_Solver *solver)
AMG Create passing function.
Hypre_Solver
Enumerated type for Hypre solvers.
int Hypre_EuclidCreate(MPI_Comm comm, HYPRE_Solver *solver)
Euclid Create passing function.
int Hypre_ParCSRHybridCreate(MPI_Comm comm, HYPRE_Solver *solver)
Hybrid Create passing function.
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 SetParameter(Hypre_Chooser chooser, int(*pt2Func)(HYPRE_Solver, int), int parameter)
Set a parameter that takes a single int.
int SetupSolver() const
int Hypre_ParCSRBiCGSTABCreate(MPI_Comm comm, HYPRE_Solver *solver)
BiCGSTAB Create passing function.
int CreateSolver()
Create the solver with selected type.