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