IFPACK Development
Ifpack_IC.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_IC.h"
00046 #include "Ifpack_IC_Utils.h"
00047 #include "Ifpack_Condest.h"
00048 #include "Epetra_Comm.h"
00049 #include "Epetra_Map.h"
00050 #include "Epetra_RowMatrix.h"
00051 #include "Epetra_CrsMatrix.h"
00052 #include "Epetra_Vector.h"
00053 #include "Epetra_MultiVector.h"
00054 #include "Epetra_Util.h"
00055 #include "Teuchos_ParameterList.hpp"
00056 #include "Teuchos_RefCountPtr.hpp"
00057 using Teuchos::RefCountPtr;
00058 using Teuchos::rcp;
00059 
00060 //==============================================================================
00061 Ifpack_IC::Ifpack_IC(Epetra_RowMatrix* A) :
00062   A_(rcp(A,false)),
00063   Comm_(A->Comm()),
00064   UseTranspose_(false),
00065   Condest_(-1.0),
00066   Athresh_(0.0),
00067   Rthresh_(1.0),
00068   Droptol_(0.0),
00069   Lfil_(1.0),
00070   Aict_(0),
00071   Lict_(0),
00072   Ldiag_(0),
00073   IsInitialized_(false),
00074   IsComputed_(false),
00075   NumInitialize_(0),
00076   NumCompute_(0),
00077   NumApplyInverse_(0),
00078   InitializeTime_(0.0),
00079   ComputeTime_(0),
00080   ApplyInverseTime_(0),
00081   Time_(Comm()),
00082   ComputeFlops_(0.0),
00083   ApplyInverseFlops_(0.0)
00084 {
00085   Teuchos::ParameterList List;
00086   SetParameters(List);
00087 
00088 }
00089 //==============================================================================
00090 Ifpack_IC::~Ifpack_IC()
00091 {
00092   if (Lict_ != 0) {
00093     Ifpack_AIJMatrix * Lict = (Ifpack_AIJMatrix *) Lict_;
00094     delete [] Lict->ptr;
00095     delete [] Lict->col;
00096     delete [] Lict->val;
00097     delete Lict;
00098   }
00099   if (Aict_ != 0) {
00100     Ifpack_AIJMatrix * Aict = (Ifpack_AIJMatrix *) Aict_;
00101     delete Aict;
00102   }
00103   if (Ldiag_ != 0) delete [] Ldiag_;
00104 
00105   IsInitialized_ = false;
00106   IsComputed_ = false;
00107 }
00108 
00109 //==========================================================================
00110 int Ifpack_IC::SetParameters(Teuchos::ParameterList& List)
00111 {
00112 
00113   // Lfil_ = List.get("fact: level-of-fill",Lfil_); // Confusing parameter since Ifpack_IC can only do ICT not IC(k)
00114   Lfil_ = List.get("fact: ict level-of-fill", Lfil_); // Lfil_ is now the fill ratio as in ICT (1.0 means same #nonzeros as A)
00115   Athresh_ = List.get("fact: absolute threshold", Athresh_);
00116   Rthresh_ = List.get("fact: relative threshold", Rthresh_);
00117   Droptol_ = List.get("fact: drop tolerance", Droptol_);
00118 
00119   // set label
00120   sprintf(Label_, "IFPACK IC (fill=%f, drop=%f)",
00121       Lfil_, Droptol_);
00122   return(0);
00123 }
00124 
00125 //==========================================================================
00126 int Ifpack_IC::Initialize()
00127 {
00128 
00129   IsInitialized_ = false;
00130 
00131   // FIXME: construction of ILUK graph must be here
00132   
00133   ++NumInitialize_;
00134   IsInitialized_ = true;
00135   return(0);
00136 }
00137 
00138 //==========================================================================
00139 int Ifpack_IC::ComputeSetup()
00140 {
00141   // (re)allocate memory for ICT factors
00142   U_ = rcp(new Epetra_CrsMatrix(Copy, Matrix().RowMatrixRowMap(), 
00143                                 Matrix().RowMatrixRowMap(), 0));
00144   D_ = rcp(new Epetra_Vector(Matrix().RowMatrixRowMap()));
00145   
00146   if (U_.get() == 0 || D_.get() == 0)
00147     IFPACK_CHK_ERR(-5); // memory allocation error
00148 
00149   int ierr = 0;
00150   int i, j;
00151   int NumIn, NumU;
00152   bool DiagFound;
00153   int NumNonzeroDiags = 0;
00154 
00155   // Get Maximun Row length
00156   int MaxNumEntries = Matrix().MaxNumEntries();
00157 
00158   std::vector<int> InI(MaxNumEntries); // Allocate temp space
00159   std::vector<int> UI(MaxNumEntries);
00160   std::vector<double> InV(MaxNumEntries);
00161   std::vector<double> UV(MaxNumEntries);
00162 
00163   double *DV;
00164   ierr = D_->ExtractView(&DV); // Get view of diagonal
00165 
00166   // First we copy the user's matrix into diagonal vector and U, regardless of fill level
00167 
00168   int NumRows = Matrix().NumMyRows();
00169 
00170   for (i=0; i< NumRows; i++) {
00171 
00172     Matrix().ExtractMyRowCopy(i, MaxNumEntries, NumIn, &InV[0], &InI[0]); // Get Values and Indices
00173     
00174     // Split into L and U (we don't assume that indices are ordered).
00175     NumU = 0; 
00176     DiagFound = false;
00177     
00178     for (j=0; j< NumIn; j++) {
00179       int k = InI[j];
00180 
00181       if (k==i) {
00182     DiagFound = true;
00183     // Store perturbed diagonal in Epetra_Vector D_
00184     DV[i] += Rthresh_ * InV[j] + EPETRA_SGN(InV[j]) * Athresh_; 
00185       }
00186 
00187       else if (k < 0) return(-1); // Out of range
00188       else if (i<k && k<NumRows) {
00189     UI[NumU] = k;
00190     UV[NumU] = InV[j];
00191     NumU++;
00192       }
00193     }
00194     
00195     // Check in things for this row of L and U
00196 
00197     if (DiagFound) NumNonzeroDiags++;
00198     if (NumU) U_->InsertMyValues(i, NumU, &UV[0], &UI[0]);
00199     
00200   }
00201 
00202   U_->FillComplete(Matrix().OperatorDomainMap(), 
00203            Matrix().OperatorRangeMap());
00204   
00205   int ierr1 = 0;
00206   if (NumNonzeroDiags<U_->NumMyRows()) ierr1 = 1;
00207   Matrix().Comm().MaxAll(&ierr1, &ierr, 1);
00208   IFPACK_CHK_ERR(ierr);
00209   
00210   return(0);
00211 }
00212 
00213 //==========================================================================
00214 int Ifpack_IC::Compute() {
00215 
00216   if (!IsInitialized()) 
00217     IFPACK_CHK_ERR(Initialize());
00218 
00219   Time_.ResetStartTime();
00220   IsComputed_ = false;
00221   
00222   // copy matrix into L and U.
00223   IFPACK_CHK_ERR(ComputeSetup());
00224   
00225   int i;
00226   
00227   int m, n, nz, Nrhs, ldrhs, ldlhs;
00228   int * ptr=0, * ind;
00229   double * val, * rhs, * lhs;
00230   
00231   int ierr = Epetra_Util_ExtractHbData(U_.get(), 0, 0, m, n, nz, ptr, ind,
00232                        val, Nrhs, rhs, ldrhs, lhs, ldlhs);
00233   if (ierr < 0) 
00234     IFPACK_CHK_ERR(ierr);
00235   
00236   Ifpack_AIJMatrix * Aict;
00237   if (Aict_==0) {
00238     Aict = new Ifpack_AIJMatrix;
00239     Aict_ = (void *) Aict;
00240   }
00241   else Aict = (Ifpack_AIJMatrix *) Aict_;
00242   Ifpack_AIJMatrix * Lict;
00243   if (Lict_==0) {
00244     Lict = new Ifpack_AIJMatrix;
00245     Lict_ = (void *) Lict;
00246   }
00247   else {
00248     Lict = (Ifpack_AIJMatrix *) Lict_;
00249     Ifpack_AIJMatrix_dealloc( Lict );  // Delete storage, crout_ict will allocate it again.
00250   }
00251   if (Ldiag_ != 0) delete [] Ldiag_; // Delete storage, crout_ict will allocate it again.
00252   Aict->val = val;
00253   Aict->col = ind;
00254   Aict->ptr = ptr;
00255   double *DV;
00256   EPETRA_CHK_ERR(D_->ExtractView(&DV)); // Get view of diagonal
00257     
00258   // lfil is average number of nonzeros per row times fill ratio.
00259   // Currently crout_ict keeps a constant number of nonzeros per row.
00260   // TODO: Pass Lfil_ and make crout_ict keep variable #nonzeros per row.
00261   int lfil = (Lfil_ * nz)/n + .5;
00262 
00263   crout_ict(m, Aict, DV, Droptol_, lfil, Lict, &Ldiag_);
00264 
00265   // Get rid of unnecessary data
00266   delete [] ptr;
00267   
00268   // Create Epetra View of L from crout_ict
00269   U_ = rcp(new Epetra_CrsMatrix(View, A_->RowMatrixRowMap(), A_->RowMatrixRowMap(),0));
00270   D_ = rcp(new Epetra_Vector(View, A_->RowMatrixRowMap(), Ldiag_));
00271   
00272   ptr = Lict->ptr;
00273   ind = Lict->col;
00274   val = Lict->val;
00275   
00276   for (i=0; i< m; i++) {
00277     int NumEntries = ptr[i+1]-ptr[i];
00278     int * Indices = ind+ptr[i];
00279     double * Values = val+ptr[i];
00280     U_->InsertMyValues(i, NumEntries, Values, Indices);
00281   }
00282   
00283   U_->FillComplete(A_->OperatorDomainMap(), A_->OperatorRangeMap());
00284   D_->Reciprocal(*D_); // Put reciprocal of diagonal in this vector
00285   
00286 #ifdef IFPACK_FLOPCOUNTERS
00287   double current_flops = 2 * nz; // Just an estimate
00288   double total_flops = 0;
00289   
00290   A_->Comm().SumAll(&current_flops, &total_flops, 1); // Get total madds across all PEs
00291   
00292   ComputeFlops_ += total_flops; 
00293   // Now count the rest. NOTE: those flops are *global*
00294   ComputeFlops_ += (double) U_->NumGlobalNonzeros(); // Accounts for multiplier above
00295   ComputeFlops_ += (double) D_->GlobalLength(); // Accounts for reciprocal of diagonal
00296 #endif
00297   ++NumCompute_;
00298   ComputeTime_ += Time_.ElapsedTime();
00299  
00300  
00301   IsComputed_ = true;
00302   
00303   return(0);
00304   
00305 }
00306 
00307 //=============================================================================
00308 // This function finds Y such that LDU Y = X or U(trans) D L(trans) Y = X for multiple RHS
00309 int Ifpack_IC::ApplyInverse(const Epetra_MultiVector& X, 
00310                 Epetra_MultiVector& Y) const
00311 {
00312   
00313   if (!IsComputed())
00314     IFPACK_CHK_ERR(-3); // compute preconditioner first
00315   
00316   if (X.NumVectors() != Y.NumVectors()) 
00317     IFPACK_CHK_ERR(-2); // Return error: X and Y not the same size
00318  
00319   Time_.ResetStartTime(); 
00320 
00321   bool Upper = true;
00322   bool UnitDiagonal = true;
00323   
00324   // AztecOO gives X and Y pointing to the same memory location,
00325   // need to create an auxiliary vector, Xcopy
00326   RefCountPtr< const Epetra_MultiVector > Xcopy;
00327   if (X.Pointers()[0] == Y.Pointers()[0])
00328     Xcopy = rcp( new Epetra_MultiVector(X) );
00329   else
00330     Xcopy = rcp( &X, false );
00331   
00332   U_->Solve(Upper, true, UnitDiagonal, *Xcopy, Y);
00333   Y.Multiply(1.0, *D_, Y, 0.0); // y = D*y (D_ has inverse of diagonal)
00334   U_->Solve(Upper, false, UnitDiagonal, Y, Y); // Solve Uy = y
00335   
00336 #ifdef IFPACK_FLOPCOUNTERS
00337   ApplyInverseFlops_ += 4.0 * U_->NumGlobalNonzeros();
00338   ApplyInverseFlops_ += D_->GlobalLength();
00339 #endif
00340 
00341   ++NumApplyInverse_;
00342   ApplyInverseTime_ += Time_.ElapsedTime();
00343 
00344   return(0);
00345 
00346 }
00347 
00348 //=============================================================================
00349 // This function finds X such that LDU Y = X or U(trans) D L(trans) Y = X for multiple RHS
00350 int Ifpack_IC::Apply(const Epetra_MultiVector& X, 
00351               Epetra_MultiVector& Y) const 
00352 {
00353 
00354   if (X.NumVectors() != Y.NumVectors()) 
00355     IFPACK_CHK_ERR(-2); // Return error: X and Y not the same size
00356 
00357   Epetra_MultiVector * X1 = (Epetra_MultiVector *) &X;
00358   Epetra_MultiVector * Y1 = (Epetra_MultiVector *) &Y;
00359   
00360   U_->Multiply(false, *X1, *Y1);
00361   Y1->Update(1.0, *X1, 1.0); // Y1 = Y1 + X1 (account for implicit unit diagonal)
00362   Y1->ReciprocalMultiply(1.0, *D_, *Y1, 0.0); // y = D*y (D_ has inverse of diagonal)
00363   Epetra_MultiVector Y1temp(*Y1); // Need a temp copy of Y1
00364   U_->Multiply(true, Y1temp, *Y1);
00365   Y1->Update(1.0, Y1temp, 1.0); // (account for implicit unit diagonal)
00366   return(0);
00367 }
00368 
00369 //=============================================================================
00370 double Ifpack_IC::Condest(const Ifpack_CondestType CT, 
00371               const int MaxIters, const double Tol,
00372               Epetra_RowMatrix* Matrix_in)
00373 {
00374   if (!IsComputed()) // cannot compute right now
00375     return(-1.0);
00376   
00377   if (Condest_ == -1.0)
00378     Condest_ = Ifpack_Condest(*this, CT, MaxIters, Tol, Matrix_in);
00379   
00380   return(Condest_);
00381 }
00382 
00383 //=============================================================================
00384 std::ostream&
00385 Ifpack_IC::Print(std::ostream& os) const
00386 {
00387   if (!Comm().MyPID()) {
00388     os << endl;
00389     os << "================================================================================" << endl;
00390     os << "Ifpack_IC: " << Label() << endl << endl;
00391     os << "Level-of-fill      = " << LevelOfFill() << endl;
00392     os << "Absolute threshold = " << AbsoluteThreshold() << endl;
00393     os << "Relative threshold = " << RelativeThreshold() << endl;
00394     os << "Drop tolerance     = " << DropTolerance() << endl;
00395     os << "Condition number estimate = " << Condest() << endl;
00396     os << "Global number of rows            = " << A_->NumGlobalRows64() << endl;
00397     if (IsComputed_) {
00398       os << "Number of nonzeros of H         = " << U_->NumGlobalNonzeros64() << endl;
00399       os << "nonzeros / rows                 = " 
00400          << 1.0 * U_->NumGlobalNonzeros64() / U_->NumGlobalRows64() << endl;
00401     }
00402     os << endl;
00403     os << "Phase           # calls   Total Time (s)       Total MFlops     MFlops/s" << endl;
00404     os << "-----           -------   --------------       ------------     --------" << endl;
00405     os << "Initialize()    "   << std::setw(5) << NumInitialize() 
00406        << "  " << std::setw(15) << InitializeTime() 
00407        << "               0.0            0.0" << endl;
00408     os << "Compute()       "   << std::setw(5) << NumCompute() 
00409        << "  " << std::setw(15) << ComputeTime()
00410        << "  " << std::setw(15) << 1.0e-6 * ComputeFlops();
00411     if (ComputeTime() != 0.0)
00412       os << "  " << std::setw(15) << 1.0e-6 * ComputeFlops() / ComputeTime() << endl;
00413     else
00414       os << "  " << std::setw(15) << 0.0 << endl;
00415     os << "ApplyInverse()  "   << std::setw(5) << NumApplyInverse() 
00416        << "  " << std::setw(15) << ApplyInverseTime()
00417        << "  " << std::setw(15) << 1.0e-6 * ApplyInverseFlops();
00418     if (ApplyInverseTime() != 0.0) 
00419       os << "  " << std::setw(15) << 1.0e-6 * ApplyInverseFlops() / ApplyInverseTime() << endl;
00420     else
00421       os << "  " << std::setw(15) << 0.0 << endl;
00422     os << "================================================================================" << endl;
00423     os << endl;
00424   }
00425 
00426   
00427   return(os);
00428 } 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Friends