Ifpack_ILU.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 #include "Ifpack_ConfigDefs.h"
00030 #ifdef HAVE_IFPACK_TEUCHOS
00031 #include "Ifpack_CondestType.h"
00032 #include "Ifpack_ILU.h"
00033 #include "Epetra_ConfigDefs.h"
00034 #include "Epetra_Comm.h"
00035 #include "Epetra_Map.h"
00036 #include "Epetra_RowMatrix.h"
00037 #include "Epetra_Vector.h"
00038 #include "Epetra_MultiVector.h"
00039 #include "Epetra_CrsGraph.h"
00040 #include "Epetra_CrsMatrix.h"
00041 #include "Teuchos_ParameterList.hpp"
00042 using namespace Teuchos;
00043 
00044 //==============================================================================
00045 Ifpack_ILU::Ifpack_ILU(Epetra_RowMatrix* Matrix) :
00046   A_(Matrix),
00047   Comm_(Matrix->Comm()),
00048   UseTranspose_(false),
00049   NumMyDiagonals_(0),
00050   RelaxValue_(0.0),
00051   Athresh_(0.0),
00052   Rthresh_(1.0),
00053   Condest_(-1.0),
00054   LevelOfFill_(0),
00055   IsInitialized_(false),
00056   IsComputed_(false),
00057   NumInitialize_(0),
00058   NumCompute_(0),
00059   NumApplyInverse_(0),
00060   InitializeTime_(0.0),
00061   ComputeTime_(0.0),
00062   ApplyInverseTime_(0.0),
00063   ComputeFlops_(0.0),
00064   ApplyInverseFlops_(0.0),
00065   Time_(Comm())
00066 {
00067 #ifdef HAVE_IFPACK_TEUCHOS
00068   Teuchos::ParameterList List;
00069   SetParameters(List);
00070 #endif
00071 }
00072 
00073 //==============================================================================
00074 void Ifpack_ILU::Destroy()
00075 {
00076 /*
00077   if (L_)             delete L_;             L_ = 0;
00078   if (U_)             delete U_;             U_ = 0;
00079   if (D_)             delete D_;             D_ = 0;
00080   if (Graph_)         delete Graph_;         Graph_ = 0;
00081   if (CrsGraph_)      delete CrsGraph_;      CrsGraph_ = 0;
00082   if (L_Graph_)       delete L_Graph_;       L_Graph_ = 0;
00083   if (U_Graph_)       delete U_Graph_;       U_Graph_ = 0;
00084   if (IlukRowMap_)    delete IlukRowMap_;    IlukRowMap_ = 0;
00085   if (IlukDomainMap_) delete IlukDomainMap_; IlukDomainMap_ = 0;
00086   if (IlukRangeMap_)  delete IlukRangeMap_;  IlukRangeMap_ = 0;
00087 */
00088   // reset pointers to already allocated stuff
00089   U_DomainMap_ = 0;
00090   L_RangeMap_ = 0;
00091 }
00092 
00093 //==========================================================================
00094 int Ifpack_ILU::SetParameters(Teuchos::ParameterList& List)
00095 {
00096   RelaxValue_ = List.get("fact: relax value", RelaxValue_);
00097   Athresh_ = List.get("fact: absolute threshold", Athresh_);
00098   Rthresh_ = List.get("fact: relative threshold", Rthresh_);
00099   LevelOfFill_ = List.get("fact: level-of-fill", LevelOfFill_);
00100 
00101   // set label
00102   sprintf(Label_, "IFPACK ILU (fill=%d, relax=%f, athr=%f, rthr=%f)",
00103       LevelOfFill(), RelaxValue(), AbsoluteThreshold(), 
00104       RelativeThreshold());
00105   return(0);
00106 }
00107 
00108 //==========================================================================
00109 int Ifpack_ILU::ComputeSetup() 
00110 {
00111 
00112   L_ = rcp(new Epetra_CrsMatrix(Copy, Graph().L_Graph()));
00113   U_ = rcp(new Epetra_CrsMatrix(Copy, Graph().U_Graph()));
00114   D_ = rcp(new Epetra_Vector(Graph().L_Graph().RowMap()));
00115   if ((L_.get() == 0) || (U_.get() == 0) || (D_.get() == 0))
00116     IFPACK_CHK_ERR(-5);
00117 
00118   // Get Maximun Row length
00119   int MaxNumEntries = Matrix().MaxNumEntries();
00120 
00121   // Set L range map and U domain map
00122   U_DomainMap_ = &(Matrix().OperatorDomainMap());
00123   L_RangeMap_ = &(Matrix().OperatorRangeMap());
00124  
00125   // this is the old InitAllValues()
00126   int ierr = 0;
00127   int i, j;
00128   int * InI=0, * LI=0, * UI = 0;
00129   double * InV=0, * LV=0, * UV = 0;
00130   int NumIn, NumL, NumU;
00131   bool DiagFound;
00132   int NumNonzeroDiags = 0;
00133 
00134   InI = new int[MaxNumEntries]; // Allocate temp space
00135   LI = new int[MaxNumEntries];
00136   UI = new int[MaxNumEntries];
00137   InV = new double[MaxNumEntries];
00138   LV = new double[MaxNumEntries];
00139   UV = new double[MaxNumEntries];
00140 
00141   bool ReplaceValues = (L_->StaticGraph() || L_->IndicesAreLocal()); // Check if values should be inserted or replaced
00142 
00143   if (ReplaceValues) {
00144     L_->PutScalar(0.0); // Zero out L and U matrices
00145     U_->PutScalar(0.0);
00146   }
00147 
00148   D_->PutScalar(0.0); // Set diagonal values to zero
00149   double *DV;
00150   IFPACK_CHK_ERR(D_->ExtractView(&DV)); // Get view of diagonal
00151 
00152   // First we copy the user's matrix into L and U, regardless of fill level
00153 
00154   for (i=0; i< NumMyRows(); i++) {
00155 
00156     IFPACK_CHK_ERR(Matrix().ExtractMyRowCopy(i, MaxNumEntries, NumIn, InV, InI)); // Get Values and Indices
00157     
00158     // Split into L and U (we don't assume that indices are ordered).
00159     
00160     NumL = 0; 
00161     NumU = 0; 
00162     DiagFound = false;
00163     
00164     for (j=0; j< NumIn; j++) {
00165       int k = InI[j];
00166 
00167       if (k==i) {
00168     DiagFound = true;
00169     DV[i] += Rthresh_ * InV[j] + EPETRA_SGN(InV[j]) * Athresh_; // Store perturbed diagonal in Epetra_Vector D_
00170       }
00171 
00172       else if (k < 0) {IFPACK_CHK_ERR(-4);} // Out of range
00173 
00174       else if (k < i) {
00175     LI[NumL] = k;
00176     LV[NumL] = InV[j];
00177     NumL++;
00178       }
00179       else if (k<NumMyRows()) {
00180     UI[NumU] = k;
00181     UV[NumU] = InV[j];
00182     NumU++;
00183       }
00184     }
00185     
00186     // Check in things for this row of L and U
00187 
00188     if (DiagFound) NumNonzeroDiags++;
00189     else DV[i] = Athresh_;
00190 
00191     if (NumL) {
00192       if (ReplaceValues) {
00193     (L_->ReplaceMyValues(i, NumL, LV, LI));
00194       }
00195       else {
00196     IFPACK_CHK_ERR(L_->InsertMyValues(i, NumL, LV, LI));
00197       }
00198     }
00199 
00200     if (NumU) {
00201       if (ReplaceValues) {
00202     (U_->ReplaceMyValues(i, NumU, UV, UI));
00203       }
00204       else {
00205     IFPACK_CHK_ERR(U_->InsertMyValues(i, NumU, UV, UI));
00206       }
00207     }
00208     
00209   }
00210 
00211   delete [] LI;
00212   delete [] UI;
00213   delete [] LV;
00214   delete [] UV;
00215   delete [] InI;
00216   delete [] InV;
00217 
00218   if (!ReplaceValues) {
00219     // The domain of L and the range of U are exactly their own row maps (there is no communication).
00220     // The domain of U and the range of L must be the same as those of the original matrix,
00221     // However if the original matrix is a VbrMatrix, these two latter maps are translation from
00222     // a block map to a point map.
00223     IFPACK_CHK_ERR(L_->FillComplete((L_->RowMatrixColMap()), *L_RangeMap_));
00224     IFPACK_CHK_ERR(U_->FillComplete(*U_DomainMap_, U_->RowMatrixRowMap()));
00225   }
00226 
00227   // At this point L and U have the values of A in the structure of L and U, and diagonal vector D
00228 
00229 
00230   int TotalNonzeroDiags = 0;
00231   IFPACK_CHK_ERR(Graph().L_Graph().RowMap().Comm().SumAll(&NumNonzeroDiags, &TotalNonzeroDiags, 1));
00232   NumMyDiagonals_ = NumNonzeroDiags;
00233   if (NumNonzeroDiags != NumMyRows()) ierr = 1; // Diagonals are not right, warn user
00234 
00235   IFPACK_CHK_ERR(ierr);
00236   return(ierr);
00237 }
00238 
00239 //==========================================================================
00240 int Ifpack_ILU::Initialize() 
00241 {
00242   Time_.ResetStartTime();
00243   IsInitialized_ = false;
00244 
00245   // reset this object
00246   Destroy();
00247 
00248   Epetra_CrsMatrix* CrsMatrix;
00249   CrsMatrix = dynamic_cast<Epetra_CrsMatrix*>(A_);
00250   if (CrsMatrix == 0) {
00251     // this means that we have to create
00252     // the graph from a given Epetra_RowMatrix. Note
00253     // that at this point we are ignoring any possible
00254     // graph coming from VBR matrices.
00255     int size = A_->MaxNumEntries();
00256     CrsGraph_ = rcp(new Epetra_CrsGraph(Copy,A_->RowMatrixRowMap(), size));
00257     if (CrsGraph_.get() == 0)
00258       IFPACK_CHK_ERR(-5); // memory allocation error
00259 
00260     vector<int> Indices(size);
00261     vector<double> Values(size);
00262 
00263     // extract each row at-a-time, and insert it into
00264     // the graph, ignore all off-process entries
00265     for (int i = 0 ; i < A_->NumMyRows() ; ++i) {
00266       int NumEntries;
00267       int GlobalRow = A_->RowMatrixRowMap().GID(i);
00268       IFPACK_CHK_ERR(A_->ExtractMyRowCopy(i, size, NumEntries, 
00269                       &Values[0], &Indices[0]));
00270       // convert to global indices
00271       for (int j = 0 ; j < NumEntries ; ++j) {
00272     Indices[j] = A_->RowMatrixColMap().GID(Indices[j]); 
00273       }
00274       IFPACK_CHK_ERR(CrsGraph_->InsertGlobalIndices(GlobalRow,NumEntries,
00275                            &Indices[0]));
00276     }
00277     
00278     IFPACK_CHK_ERR(CrsGraph_->FillComplete(A_->RowMatrixRowMap(),
00279                       A_->RowMatrixRowMap()));
00280 
00281     // always overlap zero, wider overlap will be handled
00282     // by the AdditiveSchwarz preconditioner.
00283     Graph_ = rcp(new Ifpack_IlukGraph(*CrsGraph_, LevelOfFill_, 0));
00284 
00285   }
00286   else {
00287     // see comment above for the overlap.
00288     Graph_ = rcp(new Ifpack_IlukGraph(CrsMatrix->Graph(), LevelOfFill_, 0));
00289   }
00290 
00291   if (Graph_.get() == 0)
00292     IFPACK_CHK_ERR(-5); // memory allocation error
00293   IFPACK_CHK_ERR(Graph_->ConstructFilledGraph());
00294 
00295   IsInitialized_ = true;
00296   NumInitialize_++;
00297   InitializeTime_ += Time_.ElapsedTime();
00298 
00299   return(0);
00300 }
00301 
00302 //==========================================================================
00303 int Ifpack_ILU::Compute() 
00304 {
00305 
00306   if (!IsInitialized()) 
00307     IFPACK_CHK_ERR(Initialize());
00308 
00309   Time_.ResetStartTime();
00310   IsComputed_ = false;
00311 
00312   // convert Matrix() into L and U factors.
00313   IFPACK_CHK_ERR(ComputeSetup());
00314 
00315   // MinMachNum should be officially defined, for now pick something a little 
00316   // bigger than IEEE underflow value
00317 
00318   double MinDiagonalValue = Epetra_MinDouble;
00319   double MaxDiagonalValue = 1.0/MinDiagonalValue;
00320 
00321   int ierr = 0;
00322   int i, j, k;
00323   int * LI=0, * UI = 0;
00324   double * LV=0, * UV = 0;
00325   int NumIn, NumL, NumU;
00326 
00327   // Get Maximun Row length
00328   int MaxNumEntries = L_->MaxNumEntries() + U_->MaxNumEntries() + 1;
00329 
00330   int * InI = new int[MaxNumEntries]; // Allocate temp space
00331   double * InV = new double[MaxNumEntries];
00332   int * colflag = new int[NumMyCols()];
00333 
00334   double *DV;
00335   ierr = D_->ExtractView(&DV); // Get view of diagonal
00336 
00337   int current_madds = 0; // We will count multiply-add as they happen
00338 
00339   // =========================== //
00340   // Now start the factorization //
00341   // =========================== //
00342 
00343   // Need some integer workspace and pointers
00344   int NumUU; 
00345   int * UUI;
00346   double * UUV;
00347   for (j = 0; j < NumMyCols(); ++j) colflag[j] = - 1;
00348 
00349   for (i = 0; i < NumMyRows(); ++i) {
00350 
00351     // Fill InV, InI with current row of L, D and U combined
00352 
00353     NumIn = MaxNumEntries;
00354     IFPACK_CHK_ERR(L_->ExtractMyRowCopy(i, NumIn, NumL, InV, InI));
00355     LV = InV;
00356     LI = InI;
00357 
00358     InV[NumL] = DV[i]; // Put in diagonal
00359     InI[NumL] = i;
00360     
00361     IFPACK_CHK_ERR(U_->ExtractMyRowCopy(i, NumIn-NumL-1, NumU, InV+NumL+1, InI+NumL+1));
00362     NumIn = NumL+NumU+1;
00363     UV = InV+NumL+1;
00364     UI = InI+NumL+1;
00365 
00366     // Set column flags
00367     for (j=0; j<NumIn; j++) colflag[InI[j]] = j;
00368 
00369     double diagmod = 0.0; // Off-diagonal accumulator
00370 
00371     for (int jj=0; jj<NumL; jj++) {
00372       j = InI[jj];
00373       double multiplier = InV[jj]; // current_mults++;
00374 
00375       InV[jj] *= DV[j];
00376       
00377       IFPACK_CHK_ERR(U_->ExtractMyRowView(j, NumUU, UUV, UUI)); // View of row above
00378 
00379       if (RelaxValue_==0.0) {
00380     for (k=0; k<NumUU; k++) {
00381       int kk = colflag[UUI[k]];
00382       if (kk>-1) {
00383         InV[kk] -= multiplier*UUV[k];
00384         current_madds++;
00385       }
00386     }
00387       }
00388       else {
00389     for (k=0; k<NumUU; k++) {
00390       int kk = colflag[UUI[k]];
00391       if (kk>-1) InV[kk] -= multiplier*UUV[k];
00392       else diagmod -= multiplier*UUV[k];
00393       current_madds++;
00394     }
00395       }
00396      }
00397     if (NumL) {
00398       IFPACK_CHK_ERR(L_->ReplaceMyValues(i, NumL, LV, LI));  // Replace current row of L
00399     }
00400 
00401     DV[i] = InV[NumL]; // Extract Diagonal value
00402 
00403     if (RelaxValue_!=0.0) {
00404       DV[i] += RelaxValue_*diagmod; // Add off diagonal modifications
00405       // current_madds++;
00406     }
00407 
00408     if (fabs(DV[i]) > MaxDiagonalValue) {
00409       if (DV[i] < 0) DV[i] = - MinDiagonalValue;
00410       else DV[i] = MinDiagonalValue;
00411     }
00412     else
00413       DV[i] = 1.0/DV[i]; // Invert diagonal value
00414 
00415     for (j=0; j<NumU; j++) UV[j] *= DV[i]; // Scale U by inverse of diagonal
00416 
00417     if (NumU) {
00418       IFPACK_CHK_ERR(U_->ReplaceMyValues(i, NumU, UV, UI));  // Replace current row of L and U
00419     }
00420 
00421     // Reset column flags
00422     for (j=0; j<NumIn; j++) colflag[InI[j]] = -1;
00423   }
00424 
00425   // Validate that the L and U factors are actually lower and upper triangular
00426 
00427   if (!L_->LowerTriangular()) 
00428     IFPACK_CHK_ERR(-4);
00429   if (!U_->UpperTriangular()) 
00430     IFPACK_CHK_ERR(-4);
00431   
00432   // Add up flops
00433  
00434   double current_flops = 2 * current_madds;
00435   double total_flops = 0;
00436     
00437   IFPACK_CHK_ERR(Graph().L_Graph().RowMap().Comm().SumAll(&current_flops, &total_flops, 1)); // Get total madds across all PEs
00438 
00439   // Now count the rest
00440   total_flops += (double) L_->NumGlobalNonzeros(); // Accounts for multiplier above
00441   total_flops += (double) D_->GlobalLength(); // Accounts for reciprocal of diagonal
00442   if (RelaxValue_!=0.0) total_flops += 2 * (double)D_->GlobalLength(); // Accounts for relax update of diag
00443 
00444   // add to this object's counter
00445   ComputeFlops_ += total_flops;
00446   
00447   delete [] InI;
00448   delete [] InV;
00449   delete [] colflag;
00450   
00451   IsComputed_ = true;
00452   NumCompute_++;
00453   ComputeTime_ += Time_.ElapsedTime();
00454 
00455   return(ierr);
00456 
00457 }
00458 
00459 //=============================================================================
00460 // This function finds Y such that LDU Y = X or U(trans) D L(trans) Y = X for multiple RHS
00461 int Ifpack_ILU::Solve(bool Trans, const Epetra_MultiVector& X, 
00462                       Epetra_MultiVector& Y) const 
00463 {
00464 
00465   // in this function the overlap is always zero
00466   bool Upper = true;
00467   bool Lower = false;
00468   bool UnitDiagonal = true;
00469 
00470   if (!Trans) {
00471 
00472     IFPACK_CHK_ERR(L_->Solve(Lower, Trans, UnitDiagonal, X, Y));
00473     // y = D*y (D_ has inverse of diagonal)
00474     IFPACK_CHK_ERR(Y.Multiply(1.0, *D_, Y, 0.0)); 
00475     // Solve Uy = y
00476     IFPACK_CHK_ERR(U_->Solve(Upper, Trans, UnitDiagonal, Y, Y)); 
00477   }
00478   else {
00479     // Solve Uy = y
00480     IFPACK_CHK_ERR(U_->Solve(Upper, Trans, UnitDiagonal, X, Y)); 
00481     // y = D*y (D_ has inverse of diagonal)
00482     IFPACK_CHK_ERR(Y.Multiply(1.0, *D_, Y, 0.0)); 
00483     IFPACK_CHK_ERR(L_->Solve(Lower, Trans, UnitDiagonal, Y, Y));
00484   } 
00485 
00486 
00487   return(0);
00488 }
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_ILU::Multiply(bool Trans, const Epetra_MultiVector& X, 
00493                 Epetra_MultiVector& Y) const 
00494 {
00495   if (!IsComputed())
00496     IFPACK_CHK_ERR(-3);
00497 
00498   if (!Trans) {
00499     IFPACK_CHK_ERR(U_->Multiply(Trans, X, Y)); 
00500     // Y1 = Y1 + X1 (account for implicit unit diagonal)
00501     IFPACK_CHK_ERR(Y.Update(1.0, X, 1.0)); 
00502     // y = D*y (D_ has inverse of diagonal)
00503     IFPACK_CHK_ERR(Y.ReciprocalMultiply(1.0, *D_, Y, 0.0)); 
00504     Epetra_MultiVector Y1temp(Y); // Need a temp copy of Y1
00505     IFPACK_CHK_ERR(L_->Multiply(Trans, Y1temp, Y));
00506     // (account for implicit unit diagonal)
00507     IFPACK_CHK_ERR(Y.Update(1.0, Y1temp, 1.0)); 
00508   }
00509   else {
00510 
00511     IFPACK_CHK_ERR(L_->Multiply(Trans, X, Y));
00512     // Y1 = Y1 + X1 (account for implicit unit diagonal)
00513     IFPACK_CHK_ERR(Y.Update(1.0, X, 1.0)); 
00514     // y = D*y (D_ has inverse of diagonal)
00515     IFPACK_CHK_ERR(Y.ReciprocalMultiply(1.0, *D_, Y, 0.0)); 
00516     Epetra_MultiVector Y1temp(Y); // Need a temp copy of Y1
00517     IFPACK_CHK_ERR(U_->Multiply(Trans, Y1temp, Y));
00518     // (account for implicit unit diagonal)
00519     IFPACK_CHK_ERR(Y.Update(1.0, Y1temp, 1.0)); 
00520   } 
00521 
00522   return(0);
00523 }
00524 
00525 //=============================================================================
00526 // This function finds X such that LDU Y = X or U(trans) D L(trans) Y = X for multiple RHS
00527 int Ifpack_ILU::ApplyInverse(const Epetra_MultiVector& X, 
00528                              Epetra_MultiVector& Y) const
00529 {
00530 
00531   if (!IsComputed())
00532     IFPACK_CHK_ERR(-3);
00533 
00534   if (X.NumVectors() != Y.NumVectors())
00535     IFPACK_CHK_ERR(-2);
00536 
00537   Time_.ResetStartTime();
00538 
00539   // AztecOO gives X and Y pointing to the same memory location,
00540   // need to create an auxiliary vector, Xcopy
00541   const Epetra_MultiVector* Xcopy;
00542   if (X.Pointers()[0] == Y.Pointers()[0])
00543     Xcopy = new Epetra_MultiVector(X);
00544   else
00545     Xcopy = &X;
00546 
00547   IFPACK_CHK_ERR(Solve(Ifpack_ILU::UseTranspose(), *Xcopy, Y));
00548 
00549   if (Xcopy != &X)
00550     delete Xcopy;
00551 
00552   // approx is the number of nonzeros in L and U
00553   ApplyInverseFlops_ += X.NumVectors() * 4 * 
00554     (L_->NumGlobalNonzeros() + U_->NumGlobalNonzeros());
00555 
00556   ++NumApplyInverse_;
00557   ApplyInverseTime_ += Time_.ElapsedTime();
00558 
00559   return(0);
00560 
00561 }
00562 
00563 //=============================================================================
00564 double Ifpack_ILU::Condest(const Ifpack_CondestType CT, 
00565                            const int MaxIters, const double Tol,
00566                               Epetra_RowMatrix* Matrix)
00567 {
00568   if (!IsComputed()) // cannot compute right now
00569     return(-1.0);
00570 
00571   Condest_ = Ifpack_Condest(*this, CT, MaxIters, Tol, Matrix);
00572 
00573   return(Condest_);
00574 }
00575 
00576 //=============================================================================
00577 std::ostream&
00578 Ifpack_ILU::Print(std::ostream& os) const
00579 {
00580   if (!Comm().MyPID()) {
00581     os << endl;
00582     os << "================================================================================" << endl;
00583     os << "Ifpack_ILU: " << Label() << endl << endl;
00584     os << "Level-of-fill      = " << LevelOfFill() << endl;
00585     os << "Absolute threshold = " << AbsoluteThreshold() << endl;
00586     os << "Relative threshold = " << RelativeThreshold() << endl;
00587     os << "Relax value        = " << RelaxValue() << endl;
00588     os << "Condition number estimate = " << Condest() << endl;
00589     os << "Global number of rows            = " << A_->NumGlobalRows() << endl;
00590     if (IsComputed_) {
00591       os << "Number of rows of L, D, U       = " << L_->NumGlobalRows() << endl;
00592       os << "Number of nonzeros of L + U     = " << NumGlobalNonzeros() << endl;
00593       os << "nonzeros / rows                 = " 
00594         << 1.0 * NumGlobalNonzeros() / U_->NumGlobalRows() << endl;
00595     }
00596     os << endl;
00597     os << "Phase           # calls   Total Time (s)       Total MFlops     MFlops/s" << endl;
00598     os << "-----           -------   --------------       ------------     --------" << endl;
00599     os << "Initialize()    "   << std::setw(5) << NumInitialize() 
00600        << "  " << std::setw(15) << InitializeTime() 
00601        << "               0.0            0.0" << endl;
00602     os << "Compute()       "   << std::setw(5) << NumCompute() 
00603        << "  " << std::setw(15) << ComputeTime()
00604        << "  " << std::setw(15) << 1.0e-6 * ComputeFlops();
00605     if (ComputeTime() != 0.0) 
00606       os << "  " << std::setw(15) << 1.0e-6 * ComputeFlops() / ComputeTime() << endl;
00607     else
00608       os << "  " << std::setw(15) << 0.0 << endl;
00609     os << "ApplyInverse()  "   << std::setw(5) << NumApplyInverse() 
00610        << "  " << std::setw(15) << ApplyInverseTime()
00611        << "  " << std::setw(15) << 1.0e-6 * ApplyInverseFlops();
00612     if (ApplyInverseTime() != 0.0) 
00613       os << "  " << std::setw(15) << 1.0e-6 * ApplyInverseFlops() / ApplyInverseTime() << endl;
00614     else
00615       os << "  " << std::setw(15) << 0.0 << endl;
00616     os << "================================================================================" << endl;
00617     os << endl;
00618   }
00619 
00620   return(os);
00621 }
00622 #endif // IFPACK_TEUCHOS

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