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

Generated on Thu Sep 18 12:37:07 2008 for IFPACK by doxygen 1.3.9.1