Ifpack_Euclid.cpp

00001 /*@HEADER
00002 // ***********************************************************************
00003 // 
00004 //       Ifpack: Object-Oriented Algebraic Preconditioner 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 
00030 #include "Ifpack_Euclid.h"
00031 #if defined(HAVE_EUCLID) && defined(HAVE_MPI)
00032 
00033 #include "Ifpack_Utils.h"
00034 #include <algorithm>
00035 #include "Epetra_MpiComm.h"
00036 #include "Epetra_IntVector.h"
00037 #include "Epetra_Import.h"
00038 #include "Teuchos_ParameterList.hpp"
00039 #include "Teuchos_RCP.hpp"
00040 #include "getRow_dh.h"
00041 
00042 using Teuchos::RCP;
00043 using Teuchos::rcp;
00044 
00045 Ifpack_Euclid::Ifpack_Euclid(Epetra_CrsMatrix* A):
00046   A_(rcp(A,false)),
00047   UseTranspose_(false),
00048   Condest_(-1),
00049   IsInitialized_(false),
00050   IsComputed_(false),
00051   Label_(),
00052   NumInitialize_(0),
00053   NumCompute_(0),
00054   NumApplyInverse_(0),
00055   InitializeTime_(0.0),
00056   ComputeTime_(0.0),
00057   ApplyInverseTime_(0.0),
00058   ComputeFlops_(0.0),
00059   ApplyInverseFlops_(0.0),
00060   Time_(A_->Comm()),
00061   SetLevel_(1),
00062   SetBJ_(0),
00063   SetStats_(0),
00064   SetMem_(0),
00065   SetSparse_(0.0),
00066   SetRowScale_(0),
00067   SetILUT_(0.0)
00068 {
00069   // Here we need to change the view of each row to have global indices. This is
00070   // because Euclid directly extracts a row view and expects global indices.
00071   for(int i = 0; i < A_->NumMyRows(); i++){
00072     int *indices;
00073     int len;
00074     A_->Graph().ExtractMyRowView(i, len, indices);
00075     for(int j = 0; j < len; j++){
00076       indices[j] = A_->GCID(indices[j]);
00077     }
00078   }
00079 } //Constructor
00080 
00081 //==============================================================================
00082 void Ifpack_Euclid::Destroy(){
00083   // Destroy the euclid solver, only if it was setup
00084   if(IsComputed()){
00085     Euclid_dhDestroy(eu);
00086   }
00087   // Delete these euclid varaiables if they were created
00088   if(IsInitialized()){
00089     Parser_dhDestroy(parser_dh);
00090     parser_dh = NULL;
00091     TimeLog_dhDestroy(tlog_dh);
00092     tlog_dh = NULL;
00093     Mem_dhDestroy(mem_dh);
00094     mem_dh = NULL;
00095   }
00096   // Now that Euclid is done with the matrix, we change it back to having local indices.
00097   for(int i = 0; i < A_->NumMyRows(); i++){
00098     int *indices;
00099     int len;
00100     A_->Graph().ExtractMyRowView(i, len, indices);
00101     for(int j = 0; j < len; j++){
00102       indices[j] = A_->LCID(indices[j]);
00103     }
00104   }
00105 } //Destroy()
00106 
00107 //==============================================================================
00108 int Ifpack_Euclid::Initialize(){
00109   //These are global variables in Euclid
00110   comm_dh = GetMpiComm();
00111   MPI_Comm_size(comm_dh, &np_dh);
00112   MPI_Comm_rank(comm_dh, &myid_dh);
00113   Time_.ResetStartTime();
00114   if(mem_dh == NULL){
00115     Mem_dhCreate(&mem_dh);
00116   }
00117   if (tlog_dh == NULL) {
00118     TimeLog_dhCreate(&tlog_dh);
00119   }
00120 
00121   if (parser_dh == NULL) {
00122     Parser_dhCreate(&parser_dh);
00123   }
00124   Parser_dhInit(parser_dh, 0, NULL);
00125   // Create the solver, this doesn't malloc anything yet, so it's only destroyed if Compute() is called.
00126   Euclid_dhCreate(&eu);
00127   IsInitialized_=true;
00128   NumInitialize_ = NumInitialize_ + 1;
00129   InitializeTime_ = InitializeTime_ + Time_.ElapsedTime();
00130   return 0;
00131 } //Initialize()
00132 
00133 //==============================================================================
00134 int Ifpack_Euclid::SetParameters(Teuchos::ParameterList& list){
00135   List_ = list;
00136   SetLevel_ = list.get("SetLevel", (int)1);
00137   SetBJ_ = list.get("SetBJ", (int)0);
00138   SetStats_ = list.get("SetStats", (int)0);
00139   SetMem_ = list.get("SetMem", (int)0);
00140   SetSparse_ = list.get("SetSparse", (double)0.0);
00141   SetRowScale_ = list.get("SetRowScale", (int)0);
00142   SetILUT_ = list.get("SetILUT", (double)0.0);
00143   return 0;
00144 } //SetParamters()
00145 
00146 //==============================================================================
00147 int Ifpack_Euclid::SetParameter(string name, int value){
00148   //Convert to lowercase (so it's case insensitive)
00149   locale loc;
00150   for(size_t i = 0; i < name.length(); i++){
00151     name[i] = (char) tolower(name[i], loc);
00152   }
00153   if(name.compare("setlevel") == 0){
00154     SetLevel_ = value;
00155   } else if(name.compare("setbj") == 0){
00156     SetBJ_ = value;
00157   } else if(name.compare("setstats") == 0){
00158     SetStats_ = value;
00159   } else if(name.compare("setmem") == 0){
00160     SetMem_ = value;
00161   } else if(name.compare("setrowscale") == 0){
00162     SetRowScale_ = value;
00163   } else {
00164     cout << "\nThe string " << name << " is not an available option.\n";
00165     IFPACK_CHK_ERR(-1);
00166   }
00167   return 0;
00168 } //SetParameter() (int)
00169 
00170 //==============================================================================
00171 int Ifpack_Euclid::SetParameter(string name, double value){
00172   //Convert to lowercase (so it's case insensitive)
00173   locale loc;
00174   for(size_t i; i < name.length(); i++){
00175     name[i] = (char) tolower(name[i], loc);
00176   }
00177   if(name.compare("setsparse") == 0){
00178     SetSparse_ = value;
00179   } else if(name.compare("setilut") == 0){
00180     SetILUT_ = value;
00181   } else {
00182     cout << "\nThe string " << name << " is not an available option.\n";
00183     IFPACK_CHK_ERR(-1);
00184   }
00185   return 0;
00186 } //SetParameter() (double)
00187 
00188 //==============================================================================
00189 int Ifpack_Euclid::Compute(){
00190   if(IsInitialized() == false){
00191     IFPACK_CHK_ERR(Initialize());
00192   }
00193   Time_.ResetStartTime();
00194   sprintf(Label_, "IFPACK_Euclid (level=%d, bj=%d, stats=%d, mem=%d, sparse=%f, rowscale=%d, ilut=%f)",
00195       SetLevel_, SetBJ_, SetStats_, SetMem_, SetSparse_, SetRowScale_, SetILUT_);
00196   // Set the parameters
00197   eu->level = SetLevel_;
00198   if(SetBJ_ != 0){
00199     strcpy("bj", eu->algo_par);
00200   }
00201   if(SetSparse_ != 0.0){
00202     eu->sparseTolA = SetSparse_;
00203   }
00204   if(SetRowScale_ != 0){
00205     eu->isScaled = true;
00206   }
00207   if(SetILUT_ != 0.0){
00208     eu->droptol = SetILUT_;
00209   }
00210   if(SetStats_ != 0 || SetMem_ != 0){
00211     eu->logging = true;
00212     Parser_dhInsert(parser_dh, "-eu_stats", "1");
00213   }
00214   // eu->A is the matrix as a void pointer, eu->m is local rows, eu->n is global rows
00215   eu->A = (void*) A_.get();
00216   eu->m = A_->NumMyRows();
00217   eu->n = A_->NumGlobalRows();
00218   Euclid_dhSetup(eu);
00219   IsComputed_ = true;
00220   NumCompute_ = NumCompute_ + 1;
00221   ComputeTime_ = ComputeTime_ + Time_.ElapsedTime();
00222   return 0;
00223 } //Compute()
00224 
00225 //==============================================================================
00226 int Ifpack_Euclid::ApplyInverse(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const{
00227   if(IsComputed() == false){
00228     IFPACK_CHK_ERR(-1);
00229   }
00230   int NumVectors = X.NumVectors();
00231   if(NumVectors != Y.NumVectors()){
00232     IFPACK_CHK_ERR(-2);
00233   }
00234   Time_.ResetStartTime();
00235   // Loop through the vectors
00236   for(int vecNum = 0; vecNum < NumVectors; vecNum++){
00237     CallEuclid(X[vecNum], Y[vecNum]);
00238   }
00239   if(SetStats_ != 0){
00240     Euclid_dhPrintTestData(eu, stdout);
00241   }
00242   NumApplyInverse_ = NumApplyInverse_ + 1;
00243   ApplyInverseTime_ = ApplyInverseTime_ + Time_.ElapsedTime();
00244   return 0;
00245 } //ApplyInverse()
00246 
00247 //==============================================================================
00248 int Ifpack_Euclid::CallEuclid(double *x, double *y) const{
00249   Euclid_dhApply(eu, x, y);
00250   return 0;
00251 } //CallEuclid()
00252 
00253 //==============================================================================
00254 ostream& operator << (ostream& os, const Ifpack_Euclid& A){
00255   if (!A.Comm().MyPID()) {
00256     os << endl;
00257     os << "================================================================================" << endl;
00258     os << "Ifpack_Euclid: " << A.Label () << endl << endl;
00259     os << "Using " << A.Comm().NumProc() << " processors." << endl;
00260     os << "Global number of rows            = " << A.Matrix().NumGlobalRows() << endl;
00261     os << "Global number of nonzeros        = " << A.Matrix().NumGlobalNonzeros() << endl;
00262     os << "Condition number estimate = " << A.Condest() << endl;
00263     os << endl;
00264     os << "Phase           # calls   Total Time (s)       Total MFlops     MFlops/s" << endl;
00265     os << "-----           -------   --------------       ------------     --------" << endl;
00266     os << "Initialize()    "   << std::setw(5) << A.NumInitialize()
00267        << "  " << std::setw(15) << A.InitializeTime()
00268        << "              0.0              0.0" << endl;
00269     os << "Compute()       "   << std::setw(5) << A.NumCompute()
00270        << "  " << std::setw(15) << A.ComputeTime()
00271        << "  " << std::setw(15) << 1.0e-6 * A.ComputeFlops();
00272     if (A.ComputeTime() != 0.0)
00273       os << "  " << std::setw(15) << 1.0e-6 * A.ComputeFlops() / A.ComputeTime() << endl;
00274     else
00275       os << "  " << std::setw(15) << 0.0 << endl;
00276     os << "ApplyInverse()  "   << std::setw(5) << A.NumApplyInverse()
00277        << "  " << std::setw(15) << A.ApplyInverseTime()
00278        << "  " << std::setw(15) << 1.0e-6 * A.ApplyInverseFlops();
00279     if (A.ApplyInverseTime() != 0.0)
00280       os << "  " << std::setw(15) << 1.0e-6 * A.ApplyInverseFlops() / A.ApplyInverseTime() << endl;
00281     else
00282       os << "  " << std::setw(15) << 0.0 << endl;
00283     os << "================================================================================" << endl;
00284     os << endl;
00285   }
00286   return os;
00287 } // <<
00288 
00289 //==============================================================================
00290 double Ifpack_Euclid::Condest(const Ifpack_CondestType CT, 
00291                              const int MaxIters,
00292                              const double Tol,
00293                              Epetra_RowMatrix* Matrix_in){
00294   if (!IsComputed()) // cannot compute right now
00295     return(-1.0);
00296   return(Condest_);
00297 } //Condest() - not implemented
00298 
00299 #endif // HAVE_EUCLID && HAVE_MPI

Generated on Tue Jul 13 09:27:13 2010 for IFPACK by  doxygen 1.4.7