Ifpack_ICT.cpp

00001 /*@HEADER
00002 // ***********************************************************************
00003 // 
00004 //       Ifpack: Object-Oriented Algebraic Preconditioner Package
00005 //                 Copyright (2002) 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_ConfigDefs.h"
00031 #include "Ifpack_Preconditioner.h"
00032 #include "Ifpack_ICT.h"
00033 #include "Ifpack_Condest.h"
00034 #include "Ifpack_Utils.h"
00035 #include "Ifpack_HashTable.h"
00036 #include "Epetra_SerialComm.h"
00037 #include "Epetra_Comm.h"
00038 #include "Epetra_Map.h"
00039 #include "Epetra_RowMatrix.h"
00040 #include "Epetra_CrsMatrix.h"
00041 #include "Epetra_Vector.h"
00042 #include "Epetra_MultiVector.h"
00043 #include "Epetra_Util.h"
00044 #include "Teuchos_ParameterList.hpp"
00045 #include "Teuchos_RefCountPtr.hpp"
00046 
00047 //==============================================================================
00048 // FIXME: allocate Comm_ and Time_ the first Initialize() call
00049 Ifpack_ICT::Ifpack_ICT(const Epetra_RowMatrix* A) :
00050   A_(*A),
00051   Comm_(A_.Comm()),
00052   Condest_(-1.0),
00053   Athresh_(0.0),
00054   Rthresh_(1.0),
00055   LevelOfFill_(1.0),
00056   DropTolerance_(0.0),
00057   Relax_(0.0),
00058   IsInitialized_(false),
00059   IsComputed_(false),
00060   UseTranspose_(false),
00061   NumMyRows_(0),
00062   NumInitialize_(0),
00063   NumCompute_(0),
00064   NumApplyInverse_(0),
00065   InitializeTime_(0.0),
00066   ComputeTime_(0.0),
00067   ApplyInverseTime_(0.0),
00068   ComputeFlops_(0.0),
00069   ApplyInverseFlops_(0.0),
00070   Time_(Comm()),
00071   GlobalNonzeros_(0)
00072 {
00073   // do nothing here
00074 }
00075 
00076 //==============================================================================
00077 Ifpack_ICT::~Ifpack_ICT()
00078 {
00079   Destroy();
00080 }
00081 
00082 //==============================================================================
00083 void Ifpack_ICT::Destroy()
00084 {
00085   IsInitialized_ = false;
00086   IsComputed_ = false;
00087 }
00088 
00089 //==========================================================================
00090 int Ifpack_ICT::SetParameters(Teuchos::ParameterList& List)
00091 {
00092 
00093   try
00094   {
00095     LevelOfFill_ = List.get("fact: ict level-of-fill",LevelOfFill_);
00096     Athresh_ = List.get("fact: absolute threshold", Athresh_);
00097     Rthresh_ = List.get("fact: relative threshold", Rthresh_);
00098     Relax_ = List.get("fact: relax value", Relax_);
00099     DropTolerance_ = List.get("fact: drop tolerance", DropTolerance_);
00100 
00101     // set label
00102     Label_ = "ICT (fill=" + Ifpack_toString(LevelOfFill())
00103       + ", athr=" + Ifpack_toString(AbsoluteThreshold()) 
00104       + ", rthr=" + Ifpack_toString(RelativeThreshold())
00105       + ", relax=" + Ifpack_toString(RelaxValue())
00106       + ", droptol=" + Ifpack_toString(DropTolerance())
00107       + ")";
00108 
00109     return(0);
00110   }
00111   catch (...)
00112   {
00113     cerr << "Caught an exception while parsing the parameter list" << endl;
00114     cerr << "This typically means that a parameter was set with the" << endl;
00115     cerr << "wrong type (for example, int instead of double). " << endl;
00116     cerr << "please check the documentation for the type required by each parameer." << endl;
00117     IFPACK_CHK_ERR(-1);
00118   }
00119 }
00120 
00121 //==========================================================================
00122 int Ifpack_ICT::Initialize()
00123 {
00124   // clean data if present
00125   Destroy();
00126 
00127   Time_.ResetStartTime();
00128 
00129   // matrix must be square. Check only on one processor
00130   if (Comm().NumProc() == 1 && Matrix().NumMyRows() != Matrix().NumMyCols())
00131     IFPACK_CHK_ERR(-2);
00132     
00133   NumMyRows_ = Matrix().NumMyRows();
00134 
00135   // nothing else to do here
00136   IsInitialized_ = true;
00137   ++NumInitialize_;
00138   InitializeTime_ += Time_.ElapsedTime();
00139 
00140   return(0);
00141 }
00142 
00143 //==========================================================================
00144 int Ifpack_ICT::Compute() 
00145 {
00146   if (!IsInitialized()) 
00147     IFPACK_CHK_ERR(Initialize());
00148 
00149   Time_.ResetStartTime();
00150   IsComputed_ = false;
00151 
00152   NumMyRows_ = A_.NumMyRows();
00153   int Length = A_.MaxNumEntries();
00154   vector<int>    RowIndices(Length);
00155   vector<double> RowValues(Length);
00156 
00157   bool distributed = (Comm().NumProc() > 1)?true:false;
00158 
00159   if (distributed)
00160   {
00161     SerialComm_ = Teuchos::rcp(new Epetra_SerialComm);
00162     SerialMap_ = Teuchos::rcp(new Epetra_Map(NumMyRows_, 0, *SerialComm_));
00163     assert (SerialComm_.get() != 0);
00164     assert (SerialMap_.get() != 0);
00165   }
00166   else
00167     SerialMap_ = Teuchos::rcp(const_cast<Epetra_Map*>(&A_.RowMatrixRowMap()), false);
00168 
00169   int RowNnz;
00170   double flops = 0.0;
00171 
00172   H_ = Teuchos::rcp(new Epetra_CrsMatrix(Copy,*SerialMap_,0));
00173   if (H_.get() == 0)
00174     IFPACK_CHK_ERR(-5); // memory allocation error
00175 
00176   // get A(0,0) element and insert it (after sqrt)
00177   IFPACK_CHK_ERR(A_.ExtractMyRowCopy(0,Length,RowNnz,
00178                                      &RowValues[0],&RowIndices[0]));
00179 
00180   // skip off-processor elements
00181   if (distributed)
00182   {
00183     int count = 0;
00184     for (int i = 0 ;i < RowNnz ; ++i) 
00185     {
00186       if (RowIndices[i] < NumMyRows_){
00187         RowIndices[count] = RowIndices[i];
00188         RowValues[count] = RowValues[i];
00189         ++count;
00190       }
00191       else
00192         continue;
00193     }
00194     RowNnz = count;
00195   }
00196 
00197   // modify diagonal
00198   double diag_val = 0.0;
00199   for (int i = 0 ;i < RowNnz ; ++i) {
00200     if (RowIndices[i] == 0) {
00201       double& v = RowValues[i];
00202       diag_val = AbsoluteThreshold() * EPETRA_SGN(v) +
00203         RelativeThreshold() * v;
00204       break;
00205     }
00206   }
00207 
00208   diag_val = sqrt(diag_val);
00209   int diag_idx = 0;
00210   EPETRA_CHK_ERR(H_->InsertGlobalValues(0,1,&diag_val, &diag_idx));
00211 
00212   int oldSize = RowNnz;
00213 
00214   // The 10 is just a small constant to limit collisons as the actual keys
00215   // we store are the indices and not integers
00216   // [0..A_.MaxNumEntries()*LevelofFill()].
00217   Ifpack_HashTable Hash( 10 * A_.MaxNumEntries() * LevelOfFill(), 1);
00218 
00219   // start factorization for line 1
00220   for (int row_i = 1 ; row_i < NumMyRows_ ; ++row_i) {
00221 
00222     // get row `row_i' of the matrix
00223     IFPACK_CHK_ERR(A_.ExtractMyRowCopy(row_i,Length,RowNnz,
00224                                        &RowValues[0],&RowIndices[0]));
00225 
00226     // skip off-processor elements
00227     if (distributed)
00228     {
00229       int count = 0;
00230       for (int i = 0 ;i < RowNnz ; ++i) 
00231       {
00232         if (RowIndices[i] < NumMyRows_){
00233           RowIndices[count] = RowIndices[i];
00234           RowValues[count] = RowValues[i];
00235           ++count;
00236         }
00237         else
00238           continue;
00239       }
00240       RowNnz = count;
00241     }
00242 
00243     // number of nonzeros in this row are defined as the nonzeros
00244     // of the matrix, plus the level of fill 
00245     int LOF = (int)(LevelOfFill() * RowNnz);
00246     if (LOF == 0) LOF = 1;
00247 
00248     // convert line `row_i' into hash for fast access
00249     Hash.reset();
00250 
00251     double h_ii = 0.0;
00252     for (int i = 0 ; i < RowNnz ; ++i) {
00253       if (RowIndices[i] == row_i) {
00254         double& v = RowValues[i];
00255         h_ii = AbsoluteThreshold() * EPETRA_SGN(v) + RelativeThreshold() * v;
00256       }
00257       else if (RowIndices[i] < row_i)
00258       {
00259         Hash.set(RowIndices[i], RowValues[i], true);
00260       }
00261     }
00262       
00263     // form element (row_i, col_j)
00264     // I start from the first row that has a nonzero column
00265     // index in row_i.
00266     for (int col_j = RowIndices[0] ; col_j < row_i ; ++col_j) {
00267 
00268       double h_ij = 0.0, h_jj = 0.0;
00269       // note: get() returns 0.0 if col_j is not found
00270       h_ij = Hash.get(col_j);
00271 
00272       // get pointers to row `col_j'
00273       int* ColIndices;
00274       double* ColValues;
00275       int ColNnz;
00276       H_->ExtractGlobalRowView(col_j, ColNnz, ColValues, ColIndices);
00277 
00278       for (int k = 0 ; k < ColNnz ; ++k) {
00279         int col_k = ColIndices[k];
00280 
00281         if (col_k == col_j)
00282           h_jj = ColValues[k];
00283         else {
00284           double xxx = Hash.get(col_k);
00285           if (xxx != 0.0)
00286           {
00287             h_ij -= ColValues[k] * xxx;
00288             flops += 2.0;
00289           }
00290         }
00291       }
00292 
00293       h_ij /= h_jj;
00294 
00295       if (IFPACK_ABS(h_ij) > DropTolerance_)
00296       {
00297         Hash.set(col_j, h_ij);
00298       }
00299     
00300       // only approx
00301       ComputeFlops_ += 2.0 * flops + 1.0;
00302     }
00303 
00304     int size = Hash.getNumEntries();
00305 
00306     vector<double> AbsRow(size);
00307     int count = 0;
00308     
00309     // +1 because I use the extra position for diagonal in insert
00310     vector<int> keys(size + 1);
00311     vector<double> values(size + 1);
00312 
00313     Hash.arrayify(&keys[0], &values[0]);
00314 
00315     for (int i = 0 ; i < size ; ++i)
00316     {
00317       AbsRow[i] = IFPACK_ABS(values[i]);
00318     }
00319     count = size;
00320 
00321     double cutoff = 0.0;
00322     if (count > LOF) {
00323       nth_element(AbsRow.begin(), AbsRow.begin() + LOF, AbsRow.begin() + count, 
00324                   greater<double>());
00325       cutoff = AbsRow[LOF];
00326     }
00327 
00328     for (int i = 0 ; i < size ; ++i)
00329     {
00330       h_ii -= values[i] * values[i];
00331     }
00332 
00333     if (h_ii < 0.0) h_ii = 1e-12;;
00334 
00335     h_ii = sqrt(h_ii);
00336 
00337     // only approx, + 1 == sqrt
00338     ComputeFlops_ += 2 * size + 1;
00339 
00340     double DiscardedElements = 0.0;
00341 
00342     count = 0;
00343     for (int i = 0 ; i < size ; ++i)    
00344     { 
00345       if (IFPACK_ABS(values[i]) > cutoff)
00346       {
00347         values[count] = values[i];
00348         keys[count] = keys[i];
00349         ++count;
00350       }
00351       else  
00352         DiscardedElements += values[i];
00353     }
00354 
00355     if (RelaxValue() != 0.0) {
00356       DiscardedElements *= RelaxValue();
00357       h_ii += DiscardedElements;
00358     }
00359 
00360     values[count] = h_ii;
00361     keys[count] = row_i;
00362     ++count;
00363 
00364     H_->InsertGlobalValues(row_i, count, &(values[0]), (int*)&(keys[0]));
00365 
00366     oldSize = size;
00367   }
00368 
00369   IFPACK_CHK_ERR(H_->FillComplete());
00370 
00371 #if 0
00372   // to check the complete factorization
00373   Epetra_Vector LHS(Matrix().RowMatrixRowMap());
00374   Epetra_Vector RHS1(Matrix().RowMatrixRowMap());
00375   Epetra_Vector RHS2(Matrix().RowMatrixRowMap());
00376   Epetra_Vector RHS3(Matrix().RowMatrixRowMap());
00377   LHS.Random();
00378 
00379   Matrix().Multiply(false,LHS,RHS1);
00380   H_->Multiply(true,LHS,RHS2);
00381   H_->Multiply(false,RHS2,RHS3);
00382 
00383   RHS1.Update(-1.0, RHS3, 1.0);
00384   cout << endl;
00385   cout << RHS1;
00386 #endif
00387   int MyNonzeros = H_->NumGlobalNonzeros();
00388   Comm().SumAll(&MyNonzeros, &GlobalNonzeros_, 1);
00389 
00390   IsComputed_ = true;
00391   double TotalFlops; // sum across all the processors
00392   A_.Comm().SumAll(&flops, &TotalFlops, 1);
00393   ComputeFlops_ += TotalFlops;
00394   ++NumCompute_;
00395   ComputeTime_ += Time_.ElapsedTime();
00396 
00397   return(0);
00398 
00399 }
00400 
00401 //=============================================================================
00402 int Ifpack_ICT::ApplyInverse(const Epetra_MultiVector& X, 
00403                  Epetra_MultiVector& Y) const
00404 {
00405 
00406   if (!IsComputed())
00407     IFPACK_CHK_ERR(-3); // compute preconditioner first
00408 
00409   if (X.NumVectors() != Y.NumVectors()) 
00410     IFPACK_CHK_ERR(-2); // Return error: X and Y not the same size
00411 
00412   Time_.ResetStartTime();
00413 
00414   // AztecOO gives X and Y pointing to the same memory location,
00415   // need to create an auxiliary vector, Xcopy
00416   Teuchos::RefCountPtr<const Epetra_MultiVector> Xcopy;
00417   if (X.Pointers()[0] == Y.Pointers()[0])
00418     Xcopy = Teuchos::rcp( new Epetra_MultiVector(X) );
00419   else
00420     Xcopy = Teuchos::rcp( &X, false );
00421 
00422   // NOTE: H_ is based on SerialMap_, while Xcopy is based
00423   // on A.Map()... which are in general different. However, Solve()
00424   // does not seem to care... which is fine with me.
00425   //
00426   EPETRA_CHK_ERR(H_->Solve(false,false,false,*Xcopy,Y));
00427   EPETRA_CHK_ERR(H_->Solve(false,true,false,Y,Y));
00428 
00429   // these are global flop count
00430   ApplyInverseFlops_ += 4.0 * GlobalNonzeros_;
00431 
00432   ++NumApplyInverse_;
00433   ApplyInverseTime_ += Time_.ElapsedTime();
00434 
00435   return(0);
00436 }
00437 //=============================================================================
00438 // This function finds X such that LDU Y = X or U(trans) D L(trans) Y = X for multiple RHS
00439 int Ifpack_ICT::Apply(const Epetra_MultiVector& X, 
00440               Epetra_MultiVector& Y) const 
00441 {
00442 
00443   IFPACK_CHK_ERR(-98);
00444 }
00445 
00446 //=============================================================================
00447 double Ifpack_ICT::Condest(const Ifpack_CondestType CT, 
00448                             const int MaxIters, const double Tol,
00449                 Epetra_RowMatrix* Matrix_in)
00450 {
00451   if (!IsComputed()) // cannot compute right now
00452     return(-1.0);
00453 
00454   // NOTE: this is computing the *local* condest
00455   if (Condest_ == -1.0)
00456     Condest_ = Ifpack_Condest(*this, CT, MaxIters, Tol, Matrix_in);
00457 
00458   return(Condest_);
00459 }
00460 
00461 //=============================================================================
00462 std::ostream&
00463 Ifpack_ICT::Print(std::ostream& os) const
00464 {
00465   if (!Comm().MyPID()) {
00466     os << endl;
00467     os << "================================================================================" << endl;
00468     os << "Ifpack_ICT: " << Label() << endl << endl;
00469     os << "Level-of-fill      = " << LevelOfFill() << endl;
00470     os << "Absolute threshold = " << AbsoluteThreshold() << endl;
00471     os << "Relative threshold = " << RelativeThreshold() << endl;
00472     os << "Relax value        = " << RelaxValue() << endl;
00473     os << "Condition number estimate = " << Condest() << endl;
00474     os << "Global number of rows            = " << Matrix().NumGlobalRows() << endl;
00475     if (IsComputed_) {
00476       os << "Number of nonzeros of H         = " << H_->NumGlobalNonzeros() << endl;
00477       os << "nonzeros / rows                 = " 
00478          << 1.0 * H_->NumGlobalNonzeros() / H_->NumGlobalRows() << endl;
00479     }
00480     os << endl;
00481     os << "Phase           # calls   Total Time (s)       Total MFlops     MFlops/s" << endl;
00482     os << "-----           -------   --------------       ------------     --------" << endl;
00483     os << "Initialize()    "   << std::setw(5) << NumInitialize() 
00484        << "  " << std::setw(15) << InitializeTime() 
00485        << "               0.0            0.0" << endl;
00486     os << "Compute()       "   << std::setw(5) << NumCompute() 
00487        << "  " << std::setw(15) << ComputeTime()
00488        << "  " << std::setw(15) << 1.0e-6 * ComputeFlops();
00489     if (ComputeTime() != 0.0) 
00490       os << "  " << std::setw(15) << 1.0e-6 * ComputeFlops() / ComputeTime() << endl;
00491     else
00492       os << "  " << std::setw(15) << 0.0 << endl;
00493     os << "ApplyInverse()  "   << std::setw(5) << NumApplyInverse() 
00494        << "  " << std::setw(15) << ApplyInverseTime()
00495        << "  " << std::setw(15) << 1.0e-6 * ApplyInverseFlops();
00496     if (ApplyInverseTime() != 0.0)
00497       os << "  " << std::setw(15) << 1.0e-6 * ApplyInverseFlops() / ApplyInverseTime() << endl;
00498     else
00499       os << "  " << std::setw(15) << 0.0 << endl;
00500     os << "================================================================================" << endl;
00501     os << endl;
00502   }
00503 
00504   
00505   return(os);
00506 }
 All Classes Files Functions Variables Enumerations Friends
Generated on Wed Apr 13 10:05:23 2011 for IFPACK by  doxygen 1.6.3