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