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