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