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