EpetraExt Development
EpetraExt_HypreIJMatrix.cpp
Go to the documentation of this file.
00001 //@HEADER
00002 // ***********************************************************************
00003 //
00004 //     EpetraExt: Epetra Extended - Linear Algebra Services Package
00005 //                 Copyright (2011) Sandia Corporation
00006 //
00007 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
00008 // the U.S. Government retains certain rights in this software.
00009 //
00010 // Redistribution and use in source and binary forms, with or without
00011 // modification, are permitted provided that the following conditions are
00012 // met:
00013 //
00014 // 1. Redistributions of source code must retain the above copyright
00015 // notice, this list of conditions and the following disclaimer.
00016 //
00017 // 2. Redistributions in binary form must reproduce the above copyright
00018 // notice, this list of conditions and the following disclaimer in the
00019 // documentation and/or other materials provided with the distribution.
00020 //
00021 // 3. Neither the name of the Corporation nor the names of the
00022 // contributors may be used to endorse or promote products derived from
00023 // this software without specific prior written permission.
00024 //
00025 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00026 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00027 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00028 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00029 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00030 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00031 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00032 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00033 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00034 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00035 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00036 //
00037 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
00038 //
00039 // ***********************************************************************
00040 //@HEADER
00041 
00042 #include "EpetraExt_ConfigDefs.h"
00043 #ifdef HAVE_HYPRE
00044 
00045 #include "Epetra_Map.h"
00046 #include "Epetra_Vector.h"
00047 #include "Epetra_MultiVector.h"
00048 #include "EpetraExt_HypreIJMatrix.h"
00049 #include "Teuchos_RCP.hpp"
00050 #include "Teuchos_Assert.hpp"
00051 #include "Teuchos_Array.hpp"
00052 #include <set>
00053 
00054 //=======================================================
00055 EpetraExt_HypreIJMatrix::EpetraExt_HypreIJMatrix(HYPRE_IJMatrix matrix)
00056   : Epetra_BasicRowMatrix(Epetra_MpiComm(hypre_IJMatrixComm(matrix))),
00057     Matrix_(matrix),
00058     ParMatrix_(0),
00059     NumMyRows_(-1),
00060     NumGlobalRows_(-1),
00061     NumGlobalCols_(-1),
00062     MyRowStart_(-1),
00063     MyRowEnd_(-1),
00064     MatType_(-1), 
00065     TransposeSolve_(false),
00066     SolveOrPrec_(Solver)
00067 {
00068   IsSolverSetup_ = new bool[1];
00069   IsPrecondSetup_ = new bool[1];
00070   IsSolverSetup_[0] = false;
00071   IsPrecondSetup_[0] = false;
00072   // Initialize default values for global variables
00073   int ierr = 0;
00074   ierr += InitializeDefaults();
00075   TEUCHOS_TEST_FOR_EXCEPTION(ierr != 0, std::logic_error, "Couldn't initialize default values.");
00076   
00077   // Create array of global row ids
00078   Teuchos::Array<int> GlobalRowIDs;  GlobalRowIDs.resize(NumMyRows_);
00079   
00080   for (int i = MyRowStart_; i <= MyRowEnd_; i++) {
00081     GlobalRowIDs[i-MyRowStart_] = i;
00082   }
00083   
00084   // Create array of global column ids
00085   int new_value = 0; int entries = 0;
00086   std::set<int> Columns;
00087   int num_entries;
00088   double *values;
00089   int *indices;
00090   for(int i = 0; i < NumMyRows_; i++){
00091     ierr += HYPRE_ParCSRMatrixGetRow(ParMatrix_, i+MyRowStart_, &num_entries, &indices, &values);
00092     ierr += HYPRE_ParCSRMatrixRestoreRow(ParMatrix_, i+MyRowStart_,&num_entries,&indices,&values);
00093     TEUCHOS_TEST_FOR_EXCEPTION(ierr != 0, std::logic_error, "Couldn't get row of matrix.");
00094     entries = num_entries;
00095     for(int j = 0; j < num_entries; j++){
00096       // Insert column ids from this row into set
00097       new_value = indices[j];
00098       Columns.insert(new_value);
00099     }
00100   }
00101   int NumMyCols = Columns.size(); 
00102   Teuchos::Array<int> GlobalColIDs; GlobalColIDs.resize(NumMyCols);
00103   
00104   std::set<int>::iterator it;
00105   int counter = 0;
00106   for (it = Columns.begin(); it != Columns.end(); it++) {
00107     // Get column ids in order
00108     GlobalColIDs[counter] = *it;
00109     counter = counter + 1;
00110   }
00111   //printf("Proc[%d] Rows from %d to %d, num = %d\n", Comm().MyPID(), MyRowStart_,MyRowEnd_, NumMyRows_);
00112   
00113   Epetra_Map RowMap(-1, NumMyRows_, &GlobalRowIDs[0], 0, Comm());
00114   Epetra_Map ColMap(-1, NumMyCols, &GlobalColIDs[0], 0, Comm());
00115   
00116   //Need to call SetMaps()
00117   SetMaps(RowMap, ColMap);
00118  
00119   // Get an MPI_Comm to create vectors.
00120   // The vectors will be reused in Multiply(), so that they aren't recreated every time.   
00121   MPI_Comm comm;
00122   ierr += HYPRE_ParCSRMatrixGetComm(ParMatrix_, &comm);
00123   TEUCHOS_TEST_FOR_EXCEPTION(ierr != 0, std::logic_error, "Couldn't get communicator from Hypre Matrix.");
00124   
00125   ierr += HYPRE_IJVectorCreate(comm, MyRowStart_, MyRowEnd_, &X_hypre);
00126   ierr += HYPRE_IJVectorSetObjectType(X_hypre, HYPRE_PARCSR);
00127   ierr += HYPRE_IJVectorInitialize(X_hypre);
00128   ierr += HYPRE_IJVectorAssemble(X_hypre);
00129   ierr += HYPRE_IJVectorGetObject(X_hypre, (void**) &par_x);
00130   TEUCHOS_TEST_FOR_EXCEPTION(ierr != 0, std::logic_error, "Couldn't create Hypre X vector.");
00131 
00132   ierr += HYPRE_IJVectorCreate(comm, MyRowStart_, MyRowEnd_, &Y_hypre);
00133   ierr += HYPRE_IJVectorSetObjectType(Y_hypre, HYPRE_PARCSR);
00134   ierr += HYPRE_IJVectorInitialize(Y_hypre);
00135   ierr += HYPRE_IJVectorAssemble(Y_hypre);
00136   ierr += HYPRE_IJVectorGetObject(Y_hypre, (void**) &par_y);
00137   TEUCHOS_TEST_FOR_EXCEPTION(ierr != 0, std::logic_error, "Couldn't create Hypre Y vector.");
00138 
00139   x_vec = (hypre_ParVector *) hypre_IJVectorObject(((hypre_IJVector *) X_hypre));
00140   x_local = hypre_ParVectorLocalVector(x_vec);
00141 
00142   y_vec = (hypre_ParVector *) hypre_IJVectorObject(((hypre_IJVector *) Y_hypre));
00143   y_local = hypre_ParVectorLocalVector(y_vec);
00144 
00145   SolverCreatePtr_ = &EpetraExt_HypreIJMatrix::Hypre_ParCSRPCGCreate;
00146   SolverDestroyPtr_ = &HYPRE_ParCSRPCGDestroy;
00147   SolverSetupPtr_ = &HYPRE_ParCSRPCGSetup;
00148   SolverSolvePtr_ = &HYPRE_ParCSRPCGSolve;
00149   SolverPrecondPtr_ = &HYPRE_ParCSRPCGSetPrecond;
00150   CreateSolver();
00151 
00152   PrecondCreatePtr_ = &EpetraExt_HypreIJMatrix::Hypre_EuclidCreate;
00153   PrecondDestroyPtr_ = &HYPRE_EuclidDestroy;
00154   PrecondSetupPtr_ = &HYPRE_EuclidSetup;
00155   PrecondSolvePtr_ = &HYPRE_EuclidSolve;
00156   CreatePrecond();
00157   ComputeNumericConstants();
00158   ComputeStructureConstants();
00159 } //EpetraExt_HYPREIJMatrix(Hypre_IJMatrix) Constructor
00160 
00161 //=======================================================
00162 EpetraExt_HypreIJMatrix::~EpetraExt_HypreIJMatrix(){
00163   int ierr = 0;
00164   ierr += HYPRE_IJVectorDestroy(X_hypre);
00165   TEUCHOS_TEST_FOR_EXCEPTION(ierr != 0, std::logic_error, "Couldn't destroy the X Vector.");
00166   ierr += HYPRE_IJVectorDestroy(Y_hypre);
00167   TEUCHOS_TEST_FOR_EXCEPTION(ierr != 0, std::logic_error, "Couldn't destroy the Y Vector.");
00168 
00169   /* Destroy solver and preconditioner */
00170   if(IsSolverSetup_[0] == true){
00171     ierr += SolverDestroyPtr_(Solver_);
00172     TEUCHOS_TEST_FOR_EXCEPTION(ierr != 0, std::logic_error, "Couldn't destroy the Solver.");
00173   }
00174   if(IsPrecondSetup_[0] == true){
00175     ierr += PrecondDestroyPtr_(Preconditioner_);
00176     TEUCHOS_TEST_FOR_EXCEPTION(ierr != 0, std::logic_error, "Couldn't destroy the Preconditioner.");
00177   }
00178   delete[] IsSolverSetup_;
00179   delete[] IsPrecondSetup_;
00180 } //EpetraExt_HypreIJMatrix destructor
00181 
00182 //=======================================================
00183 int EpetraExt_HypreIJMatrix::ExtractMyRowCopy(int Row, int Length, int & NumEntries, 
00184                      double * Values, int * Indices) const 
00185 {
00186   // Get values and indices of ith row of matrix
00187   int *indices;
00188   double *values;
00189   int num_entries;
00190   EPETRA_CHK_ERR(HYPRE_ParCSRMatrixGetRow(ParMatrix_, Row+RowMatrixRowMap().MinMyGID(), &num_entries, &indices, &values));
00191   EPETRA_CHK_ERR(HYPRE_ParCSRMatrixRestoreRow(ParMatrix_, Row+RowMatrixRowMap().MinMyGID(), &num_entries, &indices, &values));
00192   
00193   NumEntries = num_entries;
00194   
00195   if(Length < NumEntries){
00196     printf("The arrays passed in are not large enough. Allocate more space.\n");
00197     return -2;
00198   }
00199 
00200   for(int i = 0; i < NumEntries; i++){
00201     Values[i] = values[i];
00202     Indices[i] = RowMatrixColMap().LID(indices[i]);
00203   }
00204 
00205   return 0;
00206 } //ExtractMyRowCopy()
00207 
00208 //=======================================================
00209 int EpetraExt_HypreIJMatrix::NumMyRowEntries(int Row, int & NumEntries) const 
00210 {
00211   int my_row[1];
00212   my_row[0] = Row+RowMatrixRowMap().MinMyGID();
00213   int nentries[1];
00214   EPETRA_CHK_ERR(HYPRE_IJMatrixGetRowCounts(Matrix_, 1, my_row, nentries));
00215   NumEntries = nentries[0];
00216   return 0;
00217 } //NumMyRowEntries()
00218 
00219 //=======================================================
00220 int EpetraExt_HypreIJMatrix::ExtractMyEntryView(int CurEntry, double *&Value, int &RowIndex, int &ColIndex)
00221 {/* 
00222   This gives a copy, not a view of the values, so is not implemented.
00223   if(CurEntry >= NumMyNonzeros() || CurEntry < 0){
00224     return -1;
00225   }
00226   int counter = 0;
00227   for(int i = 0; i < NumMyRows(); i++){
00228     int entries;
00229     NumMyRowEntries(i, entries);
00230     counter += entries;
00231     if(counter > CurEntry) {
00232       counter -= entries;
00233       RowIndex = i;
00234       break;
00235     }
00236   }
00237   int entries;
00238   NumMyRowEntries(RowIndex, entries);
00239   int *indices;
00240   double *values;
00241   int num_entries;
00242   HYPRE_ParCSRMatrixGetRow(ParMatrix_, RowIndex+MyRowStart_, &num_entries, &indices, &values);
00243   HYPRE_ParCSRMatrixRestoreRow(ParMatrix_, RowIndex+MyRowStart_, &num_entries, &indices, &values);
00244   ColIndex = RowMatrixColMap().LID(indices[CurEntry-counter]);
00245   Value = &values[CurEntry-counter];
00246   return 0;*/return -1;
00247 } //ExtractMyEntryView() - not implemented
00248  
00249  //=======================================================
00250  int EpetraExt_HypreIJMatrix::ExtractMyEntryView(int CurEntry, double const * & Value, int & RowIndex, int & ColIndex) const 
00251 {/*
00252    if(CurEntry >= NumMyNonzeros() || CurEntry < 0){
00253      return -1;
00254    }
00255    int counter = 0;
00256    for(int i = 0; i < NumMyRows(); i++){
00257      int entries;
00258      NumMyRowEntries(i, entries);
00259      counter += entries;
00260      if(counter > CurEntry) {
00261        counter -= entries;
00262        RowIndex = i;
00263        break;
00264      }
00265    }
00266    int entries;
00267    NumMyRowEntries(RowIndex, entries);
00268   int *indices;
00269   double *values;
00270   int num_entries;
00271   HYPRE_ParCSRMatrixGetRow(ParMatrix_, RowIndex+MyRowStart_, &num_entries, &indices, &values);
00272   HYPRE_ParCSRMatrixRestoreRow(ParMatrix_, RowIndex+MyRowStart_, &num_entries, &indices, &values);
00273   ColIndex = RowMatrixColMap().LID(indices[CurEntry-counter]);
00274   Value = &values[CurEntry-counter];
00275   return 0;*/ return -1;
00276 } //ExtractMyEntryView() - const version, not implemented
00277 
00278 //=======================================================
00279 int EpetraExt_HypreIJMatrix::Multiply(bool TransA,
00280                                const Epetra_MultiVector& X,
00281                                Epetra_MultiVector& Y) const
00282 {
00283   
00284   //printf("Proc[%d], Row start: %d, Row End: %d\n", Comm().MyPID(), MyRowStart_, MyRowEnd_);
00285   bool SameVectors = false; 
00286   int NumVectors = X.NumVectors();
00287   if (NumVectors != Y.NumVectors()) return -1;  // X and Y must have same number of vectors
00288   if(X.Pointers() == Y.Pointers()){
00289     SameVectors = true;
00290   }
00291   for(int VecNum = 0; VecNum < NumVectors; VecNum++) {
00292     //Get values for current vector in multivector.
00293     double * x_values;
00294     double * y_values;
00295     EPETRA_CHK_ERR((*X(VecNum)).ExtractView(&x_values));
00296     double *x_temp = x_local->data; 
00297     double *y_temp = y_local->data;
00298     if(!SameVectors){
00299       EPETRA_CHK_ERR((*Y(VecNum)).ExtractView(&y_values));
00300     } else {
00301       y_values = new double[X.MyLength()];
00302     }
00303     y_local->data = y_values;
00304     EPETRA_CHK_ERR(HYPRE_ParVectorSetConstantValues(par_y,0.0));
00305     // Temporarily make a pointer to data in Hypre for end
00306     // Replace data in Hypre vectors with epetra values
00307     x_local->data = x_values;
00308     // Do actual computation.
00309     if(TransA) {
00310       // Use transpose of A in multiply
00311       EPETRA_CHK_ERR(HYPRE_ParCSRMatrixMatvecT(1.0, ParMatrix_, par_x, 1.0, par_y));
00312     } else {
00313       EPETRA_CHK_ERR(HYPRE_ParCSRMatrixMatvec(1.0, ParMatrix_, par_x, 1.0, par_y));
00314     }
00315     if(SameVectors){
00316       int NumEntries = Y.MyLength();
00317       std::vector<double> new_values; new_values.resize(NumEntries);
00318       std::vector<int> new_indices; new_indices.resize(NumEntries);
00319       for(int i = 0; i < NumEntries; i++){
00320         new_values[i] = y_values[i];
00321         new_indices[i] = i;
00322       }
00323       EPETRA_CHK_ERR((*Y(VecNum)).ReplaceMyValues(NumEntries, &new_values[0], &new_indices[0]));
00324       delete[] y_values;
00325     }
00326     x_local->data = x_temp;
00327     y_local->data = y_temp;
00328   }
00329   double flops = (double) NumVectors * (double) NumGlobalNonzeros();
00330   UpdateFlops(flops);
00331   return 0;
00332 } //Multiply() 
00333 
00334 //=======================================================
00335 int EpetraExt_HypreIJMatrix::SetParameter(Hypre_Chooser chooser, int (*pt2Func)(HYPRE_Solver, int), int parameter){
00336   if(chooser == Preconditioner){
00337     EPETRA_CHK_ERR(pt2Func(Preconditioner_, parameter));
00338   }  else {
00339     EPETRA_CHK_ERR(pt2Func(Solver_, parameter));
00340   }
00341   return 0;
00342 } //SetParameter() - int function pointer
00343 
00344 //=======================================================
00345 int EpetraExt_HypreIJMatrix::SetParameter(Hypre_Chooser chooser, int (*pt2Func)(HYPRE_Solver, double), double parameter){
00346   if(chooser == Preconditioner){
00347     EPETRA_CHK_ERR(pt2Func(Preconditioner_, parameter));
00348   } else {
00349     EPETRA_CHK_ERR(pt2Func(Solver_, parameter));
00350   }
00351   return 0;
00352 } //SetParamater() - double function pointer
00353 
00354 //=======================================================
00355 int EpetraExt_HypreIJMatrix::SetParameter(Hypre_Chooser chooser, int (*pt2Func)(HYPRE_Solver, double, int), double parameter1, int parameter2){
00356   if(chooser == Preconditioner){
00357     EPETRA_CHK_ERR(pt2Func(Preconditioner_, parameter1, parameter2));
00358   } else {
00359     EPETRA_CHK_ERR(pt2Func(Solver_, parameter1, parameter2));
00360   }
00361   return 0;
00362 } //SetParameter() - double,int function pointer
00363 
00364 //=======================================================
00365 int EpetraExt_HypreIJMatrix::SetParameter(Hypre_Chooser chooser, int (*pt2Func)(HYPRE_Solver, int, int), int parameter1, int parameter2){
00366   if(chooser == Preconditioner){
00367     EPETRA_CHK_ERR(pt2Func(Preconditioner_, parameter1, parameter2));
00368   } else {
00369     EPETRA_CHK_ERR(pt2Func(Solver_, parameter1, parameter2));
00370   }
00371   return 0;
00372 } //SetParameter() - int,int function pointer
00373 
00374 //=======================================================
00375 int EpetraExt_HypreIJMatrix::SetParameter(Hypre_Chooser chooser, int (*pt2Func)(HYPRE_Solver, double*), double* parameter){
00376   if(chooser == Preconditioner){
00377     EPETRA_CHK_ERR(pt2Func(Preconditioner_, parameter));
00378   } else {
00379     EPETRA_CHK_ERR(pt2Func(Solver_, parameter));
00380   }
00381   return 0;
00382 } //SetParameter() - double* function pointer
00383 
00384 //=======================================================
00385 int EpetraExt_HypreIJMatrix::SetParameter(Hypre_Chooser chooser, int (*pt2Func)(HYPRE_Solver, int*), int* parameter){
00386   if(chooser == Preconditioner){
00387     EPETRA_CHK_ERR(pt2Func(Preconditioner_, parameter));
00388   } else {
00389     EPETRA_CHK_ERR(pt2Func(Solver_, parameter));
00390   }
00391   return 0;
00392 } //SetParameter() - int* function pointer
00393 
00394 //=======================================================
00395 int EpetraExt_HypreIJMatrix::SetParameter(Hypre_Chooser chooser, Hypre_Solver solver, bool transpose){
00396   if(chooser == Solver){
00397   if(transpose && solver != BoomerAMG){
00398     EPETRA_CHK_ERR(-1);
00399   }
00400   switch(solver) {
00401     case BoomerAMG:
00402       if(IsSolverSetup_[0]){
00403         SolverDestroyPtr_(Solver_);
00404         IsSolverSetup_[0] = false;
00405       }
00406       SolverCreatePtr_ = &EpetraExt_HypreIJMatrix::Hypre_BoomerAMGCreate;
00407       SolverDestroyPtr_ = &HYPRE_BoomerAMGDestroy;
00408       SolverSetupPtr_ = &HYPRE_BoomerAMGSetup;
00409       SolverPrecondPtr_ = NULL;
00410       if(transpose){
00411         TransposeSolve_ = true;
00412         SolverSolvePtr_ = &HYPRE_BoomerAMGSolveT;
00413       } else {
00414         SolverSolvePtr_ = &HYPRE_BoomerAMGSolve;
00415       }
00416       break;
00417     case AMS:
00418       if(IsSolverSetup_[0]){
00419         SolverDestroyPtr_(Solver_);
00420         IsSolverSetup_[0] = false;
00421       }
00422       SolverCreatePtr_ = &EpetraExt_HypreIJMatrix::Hypre_AMSCreate;
00423       SolverDestroyPtr_ = &HYPRE_AMSDestroy;
00424       SolverSetupPtr_ = &HYPRE_AMSSetup;
00425       SolverSolvePtr_ = &HYPRE_AMSSolve;
00426       SolverPrecondPtr_ = NULL;
00427       break;
00428     case Hybrid:
00429       if(IsSolverSetup_[0]){
00430         SolverDestroyPtr_(Solver_);
00431         IsSolverSetup_[0] = false;
00432       }
00433       SolverCreatePtr_ = &EpetraExt_HypreIJMatrix::Hypre_ParCSRHybridCreate;
00434       SolverDestroyPtr_ = &HYPRE_ParCSRHybridDestroy;
00435       SolverSetupPtr_ = &HYPRE_ParCSRHybridSetup;
00436       SolverSolvePtr_ = &HYPRE_ParCSRHybridSolve;
00437       SolverPrecondPtr_ = &HYPRE_ParCSRHybridSetPrecond;
00438       break;
00439     case PCG:
00440       if(IsSolverSetup_[0]){
00441         SolverDestroyPtr_(Solver_);
00442         IsSolverSetup_[0] = false;
00443       }
00444       SolverCreatePtr_ = &EpetraExt_HypreIJMatrix::Hypre_ParCSRPCGCreate;
00445       SolverDestroyPtr_ = &HYPRE_ParCSRPCGDestroy;
00446       SolverSetupPtr_ = &HYPRE_ParCSRPCGSetup;
00447       SolverSolvePtr_ = &HYPRE_ParCSRPCGSolve;
00448       SolverPrecondPtr_ = &HYPRE_ParCSRPCGSetPrecond;
00449       break;
00450     case GMRES:
00451       if(IsSolverSetup_[0]){
00452         SolverDestroyPtr_(Solver_);
00453         IsSolverSetup_[0] = false;
00454       }
00455       SolverCreatePtr_ = &EpetraExt_HypreIJMatrix::Hypre_ParCSRGMRESCreate;
00456       SolverDestroyPtr_ = &HYPRE_ParCSRGMRESDestroy;
00457       SolverSetupPtr_ = &HYPRE_ParCSRGMRESSetup;
00458       SolverSolvePtr_ = &HYPRE_ParCSRGMRESSolve;
00459       SolverPrecondPtr_ = &HYPRE_ParCSRGMRESSetPrecond;
00460       break;
00461     case FlexGMRES:
00462       if(IsSolverSetup_[0]){
00463         SolverDestroyPtr_(Solver_);
00464         IsSolverSetup_[0] = false;
00465       }
00466       SolverCreatePtr_ = &EpetraExt_HypreIJMatrix::Hypre_ParCSRFlexGMRESCreate;
00467       SolverDestroyPtr_ = &HYPRE_ParCSRFlexGMRESDestroy;
00468       SolverSetupPtr_ = &HYPRE_ParCSRFlexGMRESSetup;
00469       SolverSolvePtr_ = &HYPRE_ParCSRFlexGMRESSolve;
00470       SolverPrecondPtr_ = &HYPRE_ParCSRFlexGMRESSetPrecond;
00471       break;
00472     case LGMRES:
00473       if(IsSolverSetup_[0]){
00474         SolverDestroyPtr_(Solver_);
00475         IsSolverSetup_[0] = false;
00476       }
00477       SolverCreatePtr_ = &EpetraExt_HypreIJMatrix::Hypre_ParCSRLGMRESCreate;
00478       SolverDestroyPtr_ = &HYPRE_ParCSRLGMRESDestroy;
00479       SolverSetupPtr_ = &HYPRE_ParCSRLGMRESSetup;
00480       SolverSolvePtr_ = &HYPRE_ParCSRLGMRESSolve;
00481       SolverPrecondPtr_ = &HYPRE_ParCSRLGMRESSetPrecond;
00482       break;
00483     case BiCGSTAB:
00484       if(IsSolverSetup_[0]){
00485         SolverDestroyPtr_(Solver_);
00486         IsSolverSetup_[0] = false;
00487       }
00488       SolverCreatePtr_ = &EpetraExt_HypreIJMatrix::Hypre_ParCSRBiCGSTABCreate;
00489       SolverDestroyPtr_ = &HYPRE_ParCSRBiCGSTABDestroy;
00490       SolverSetupPtr_ = &HYPRE_ParCSRBiCGSTABSetup;
00491       SolverSolvePtr_ = &HYPRE_ParCSRBiCGSTABSolve;
00492       SolverPrecondPtr_ = &HYPRE_ParCSRBiCGSTABSetPrecond;
00493       break;
00494     default:
00495       EPETRA_CHK_ERR(-1);
00496     }
00497   CreateSolver();
00498   } else {
00499   // Preconditioner
00500   switch(solver) {
00501     case BoomerAMG:
00502       if(IsPrecondSetup_[0]){
00503         PrecondDestroyPtr_(Preconditioner_);
00504         IsPrecondSetup_[0] = false;
00505       }
00506       PrecondCreatePtr_ = &EpetraExt_HypreIJMatrix::Hypre_BoomerAMGCreate;
00507       PrecondDestroyPtr_ = &HYPRE_BoomerAMGDestroy;
00508       PrecondSetupPtr_ = &HYPRE_BoomerAMGSetup;
00509       PrecondSolvePtr_ = &HYPRE_BoomerAMGSolve;
00510       break;
00511     case ParaSails:
00512       if(IsPrecondSetup_[0]){
00513         PrecondDestroyPtr_(Preconditioner_);
00514         IsPrecondSetup_[0] = false;
00515       }
00516       PrecondCreatePtr_ = &EpetraExt_HypreIJMatrix::Hypre_ParaSailsCreate;
00517       PrecondDestroyPtr_ = &HYPRE_ParaSailsDestroy;
00518       PrecondSetupPtr_ = &HYPRE_ParaSailsSetup;
00519       PrecondSolvePtr_ = &HYPRE_ParaSailsSolve;
00520       break;
00521     case Euclid:
00522       if(IsPrecondSetup_[0]){
00523         PrecondDestroyPtr_(Preconditioner_);
00524         IsPrecondSetup_[0] = false;
00525       }
00526       PrecondCreatePtr_ = &EpetraExt_HypreIJMatrix::Hypre_EuclidCreate;
00527       PrecondDestroyPtr_ = &HYPRE_EuclidDestroy;
00528       PrecondSetupPtr_ = &HYPRE_EuclidSetup;
00529       PrecondSolvePtr_ = &HYPRE_EuclidSolve;
00530       break;
00531     case AMS:
00532       if(IsPrecondSetup_[0]){
00533         PrecondDestroyPtr_(Preconditioner_);
00534         IsPrecondSetup_[0] = false;
00535       }
00536       PrecondCreatePtr_ = &EpetraExt_HypreIJMatrix::Hypre_AMSCreate;
00537       PrecondDestroyPtr_ = &HYPRE_AMSDestroy;
00538       PrecondSetupPtr_ = &HYPRE_AMSSetup;
00539       PrecondSolvePtr_ = &HYPRE_AMSSolve;
00540       break;
00541     default:
00542       EPETRA_CHK_ERR(-1);
00543     }
00544   CreatePrecond();
00545 
00546   }
00547   return 0;
00548 } //SetParameter() - Choose solver or preconditioner type
00549 
00550 //=======================================================
00551 int EpetraExt_HypreIJMatrix::CreateSolver(){
00552   MPI_Comm comm;
00553   HYPRE_ParCSRMatrixGetComm(ParMatrix_, &comm);
00554   return (this->*SolverCreatePtr_)(comm, &Solver_);
00555 } //CreateSolver()
00556 
00557 //=======================================================
00558 int EpetraExt_HypreIJMatrix::CreatePrecond(){
00559   MPI_Comm comm;
00560   HYPRE_ParCSRMatrixGetComm(ParMatrix_, &comm);
00561   return (this->*PrecondCreatePtr_)(comm, &Preconditioner_);
00562 } //CreatePrecond()
00563 
00564 //=======================================================
00565 int EpetraExt_HypreIJMatrix::SetParameter(bool UsePreconditioner){
00566   if(UsePreconditioner == false){
00567     return 0;
00568   } else {
00569     if(SolverPrecondPtr_ == NULL){
00570       return -1;
00571     } else {
00572       SolverPrecondPtr_(Solver_, PrecondSolvePtr_, PrecondSetupPtr_, Preconditioner_);
00573       return 0;
00574     }
00575   }
00576 } //SetParameter() - Set the precondioner pointer for solver
00577 
00578 //=======================================================
00579 int EpetraExt_HypreIJMatrix::SetParameter(Hypre_Chooser answer){
00580   SolveOrPrec_ = answer;
00581   return 0;
00582 } //SetParameter() - Choose to solve or precondition
00583 
00584 //=======================================================
00585 int EpetraExt_HypreIJMatrix::SetupSolver() const{
00586   SolverSetupPtr_(Solver_, ParMatrix_, par_x, par_y);
00587   IsSolverSetup_[0] = true;
00588   return 0;
00589 } //SetupSolver()
00590 
00591 //=======================================================
00592 int EpetraExt_HypreIJMatrix::SetupPrecond() const{
00593   PrecondSetupPtr_(Preconditioner_, ParMatrix_, par_x, par_y);
00594   IsPrecondSetup_[0] = true;
00595   return 0;
00596 } //SetupPrecond()
00597 
00598 //=======================================================
00599 int EpetraExt_HypreIJMatrix::Solve(bool Upper, bool transpose, bool UnitDiagonal, const Epetra_MultiVector & X, Epetra_MultiVector & Y) const {
00600   bool SameVectors = false;
00601   int NumVectors = X.NumVectors();
00602   if (NumVectors != Y.NumVectors()) return -1;  // X and Y must have same number of vectors
00603   if(X.Pointers() == Y.Pointers()){
00604     SameVectors = true;
00605   }
00606   if(SolveOrPrec_ == Solver){
00607     if(IsSolverSetup_[0] == false){
00608       SetupSolver();
00609     }
00610   } else {
00611     if(IsPrecondSetup_[0] == false){
00612       SetupPrecond();
00613     }
00614   }
00615   for(int VecNum = 0; VecNum < NumVectors; VecNum++) {
00616     //Get values for current vector in multivector.
00617     double * x_values;
00618     EPETRA_CHK_ERR((*X(VecNum)).ExtractView(&x_values));
00619     double * y_values;
00620     if(!SameVectors){
00621       EPETRA_CHK_ERR((*Y(VecNum)).ExtractView(&y_values));
00622     } else {
00623       y_values = new double[X.MyLength()];
00624     }
00625     // Temporarily make a pointer to data in Hypre for end
00626     double *x_temp = x_local->data; 
00627     // Replace data in Hypre vectors with epetra values
00628     x_local->data = x_values;
00629     double *y_temp = y_local->data;
00630     y_local->data = y_values;
00631     
00632     EPETRA_CHK_ERR(HYPRE_ParVectorSetConstantValues(par_y, 0.0));
00633     if(transpose && !TransposeSolve_){
00634       // User requested a transpose solve, but the solver selected doesn't provide one
00635       EPETRA_CHK_ERR(-1);
00636     }
00637     if(SolveOrPrec_ == Solver){
00638       // Use the solver methods
00639       EPETRA_CHK_ERR(SolverSolvePtr_(Solver_, ParMatrix_, par_x, par_y));
00640     } else {
00641       // Apply the preconditioner
00642       EPETRA_CHK_ERR(PrecondSolvePtr_(Preconditioner_, ParMatrix_, par_x, par_y));
00643     }
00644     if(SameVectors){
00645       int NumEntries = Y.MyLength();
00646       std::vector<double> new_values; new_values.resize(NumEntries);
00647       std::vector<int> new_indices; new_indices.resize(NumEntries);
00648       for(int i = 0; i < NumEntries; i++){
00649         new_values[i] = y_values[i];
00650         new_indices[i] = i;
00651       }
00652       EPETRA_CHK_ERR((*Y(VecNum)).ReplaceMyValues(NumEntries, &new_values[0], &new_indices[0]));
00653       delete[] y_values;
00654     }
00655     x_local->data = x_temp;
00656     y_local->data = y_temp;
00657   }
00658   double flops = (double) NumVectors * (double) NumGlobalNonzeros();
00659   UpdateFlops(flops);
00660   return 0;
00661 } //Solve()
00662 
00663 //=======================================================
00664 int EpetraExt_HypreIJMatrix::LeftScale(const Epetra_Vector& X) {
00665   for(int i = 0; i < NumMyRows_; i++){
00666     //Vector-scalar mult on ith row
00667     int num_entries;
00668     int *indices;
00669     double *values;
00670     EPETRA_CHK_ERR(HYPRE_ParCSRMatrixGetRow(ParMatrix_,i+MyRowStart_, &num_entries, &indices, &values));
00671     EPETRA_CHK_ERR(HYPRE_ParCSRMatrixRestoreRow(ParMatrix_,i+MyRowStart_, &num_entries, &indices, &values));
00672     Teuchos::Array<double> new_values; new_values.resize(num_entries);
00673     Teuchos::Array<int> new_indices; new_indices.resize(num_entries);
00674     for(int j = 0; j < num_entries; j++){
00675       // Scale this row with the appropriate values from the vector
00676       new_values[j] = X[i]*values[j];
00677       new_indices[j] = indices[j];
00678     }
00679     int rows[1];
00680     rows[0] = i+MyRowStart_;
00681     EPETRA_CHK_ERR(HYPRE_IJMatrixSetValues(Matrix_, 1, &num_entries, rows, &new_indices[0], &new_values[0]));
00682     // Finally set values of the Matrix for this row
00683   }
00684   HaveNumericConstants_ = false;
00685   UpdateFlops(NumGlobalNonzeros());
00686   return 0;
00687 } //LeftScale()
00688 
00689 //=======================================================
00690 int EpetraExt_HypreIJMatrix::RightScale(const Epetra_Vector& X) {
00691   // First we need to import off-processor values of the vector
00692   Epetra_Import Importer(RowMatrixColMap(), RowMatrixRowMap());
00693   Epetra_Vector Import_Vector(RowMatrixColMap(), true);
00694   EPETRA_CHK_ERR(Import_Vector.Import(X, Importer, Insert, 0));
00695   
00696   for(int i = 0; i < NumMyRows_; i++){
00697     //Vector-scalar mult on ith column
00698     int num_entries;
00699     double *values;
00700     int *indices;
00701     // Get values and indices of ith row of matrix
00702     EPETRA_CHK_ERR(HYPRE_ParCSRMatrixGetRow(ParMatrix_,i+MyRowStart_, &num_entries, &indices, &values));
00703     EPETRA_CHK_ERR(HYPRE_ParCSRMatrixRestoreRow(ParMatrix_,i+MyRowStart_,&num_entries,&indices,&values));
00704     Teuchos::Array<int> new_indices; new_indices.resize(num_entries);
00705     Teuchos::Array<double> new_values; new_values.resize(num_entries);
00706     for(int j = 0; j < num_entries; j++){
00707       // Multiply column j with jth element
00708       int index = RowMatrixColMap().LID(indices[j]);
00709       TEUCHOS_TEST_FOR_EXCEPTION(index < 0, std::logic_error, "Index is negtive.");
00710       new_values[j] = values[j] * Import_Vector[index];
00711       new_indices[j] = indices[j];
00712     }
00713     // Finally set values of the Matrix for this row
00714     int rows[1];
00715     rows[0] = i+MyRowStart_;
00716     EPETRA_CHK_ERR(HYPRE_IJMatrixSetValues(Matrix_, 1, &num_entries, rows, &new_indices[0], &new_values[0]));
00717   }
00718   
00719   HaveNumericConstants_ = false;
00720   UpdateFlops(NumGlobalNonzeros());
00721   return 0;
00722 } //RightScale()
00723 
00724 //=======================================================
00725 int EpetraExt_HypreIJMatrix::InitializeDefaults(){
00726   int my_type;
00727   // Get type of Hypre IJ Matrix
00728   EPETRA_CHK_ERR(HYPRE_IJMatrixGetObjectType(Matrix_, &my_type));
00729   MatType_ = my_type;
00730   // Currently only ParCSR is supported
00731   TEUCHOS_TEST_FOR_EXCEPTION(MatType_ != HYPRE_PARCSR, std::logic_error, "Object is not type PARCSR");
00732 
00733   // Get the actual ParCSR object from the IJ matrix
00734   EPETRA_CHK_ERR(HYPRE_IJMatrixGetObject(Matrix_, (void**) &ParMatrix_));
00735   
00736   int numRows, numCols;
00737   
00738   // Get dimensions of the matrix and store as global variables
00739   EPETRA_CHK_ERR(HYPRE_ParCSRMatrixGetDims(ParMatrix_, &numRows, &numCols));
00740   
00741   NumGlobalRows_ = numRows;
00742   NumGlobalCols_ = numCols;
00743   
00744   // Get the local dimensions of the matrix
00745   int ColStart, ColEnd;
00746   EPETRA_CHK_ERR(HYPRE_ParCSRMatrixGetLocalRange(ParMatrix_, &MyRowStart_, &MyRowEnd_, &ColStart, &ColEnd));
00747   
00748   // Determine number of local rows
00749   NumMyRows_ = MyRowEnd_ - MyRowStart_+1;
00750   
00751   return 0;
00752 }  //InitializeDefaults() 
00753 
00754 #endif /*HAVE_HYPRE*/
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines