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