Ifpack_CrsRick.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_CrsRick.h"
00031 #include "Epetra_Comm.h"
00032 #include "Epetra_Map.h"
00033 #include "Epetra_CrsGraph.h"
00034 #include "Epetra_CrsMatrix.h"
00035 #include "Epetra_Vector.h"
00036 #include "Epetra_MultiVector.h"
00037 
00038 #ifdef HAVE_IFPACK_TEUCHOS
00039 #include <Teuchos_ParameterList.hpp>
00040 #include <ifp_parameters.h>
00041 #endif
00042 
00043 //==============================================================================
00044 Ifpack_CrsRick::Ifpack_CrsRick(const Epetra_CrsMatrix &A, const Ifpack_IlukGraph & Graph) 
00045   : A_(A),
00046     Graph_(Graph),
00047     UseTranspose_(false),
00048     Allocated_(false),
00049     ValuesInitialized_(false),
00050     Factored_(false),
00051     RelaxValue_(0.0),
00052     Condest_(-1.0),
00053     Athresh_(0.0),
00054     Rthresh_(1.0),
00055     OverlapX_(0),
00056     OverlapY_(0),
00057     OverlapMode_(Zero)
00058 {
00059   int ierr = Allocate();
00060 }
00061 
00062 //==============================================================================
00063 Ifpack_CrsRick::Ifpack_CrsRick(const Ifpack_CrsRick & FactoredMatrix) 
00064   : A_(FactoredMatrix.A_),
00065     Graph_(FactoredMatrix.Graph_),
00066     UseTranspose_(FactoredMatrix.UseTranspose_),
00067     Allocated_(FactoredMatrix.Allocated_),
00068     ValuesInitialized_(FactoredMatrix.ValuesInitialized_),
00069     Factored_(FactoredMatrix.Factored_),
00070     RelaxValue_(FactoredMatrix.RelaxValue_),
00071     Condest_(FactoredMatrix.Condest_),
00072     Athresh_(FactoredMatrix.Athresh_),
00073     Rthresh_(FactoredMatrix.Rthresh_),
00074     OverlapX_(0),
00075     OverlapY_(0),
00076     OverlapMode_(FactoredMatrix.OverlapMode_)
00077 {
00078   U_ = new Epetra_CrsMatrix(FactoredMatrix.U());
00079   D_ = new Epetra_Vector(Graph_.L_Graph().RowMap());
00080   
00081 }
00082 
00083 //==============================================================================
00084 int Ifpack_CrsRick::Allocate() {
00085 
00086   // Allocate Epetra_CrsMatrix using ILUK graphs
00087   U_ = new Epetra_CrsMatrix(Copy, Graph_.U_Graph());
00088   D_ = new Epetra_Vector(Graph_.L_Graph().RowMap());
00089   
00090   
00091     SetAllocated(true);
00092     return(0);
00093 }
00094 //==============================================================================
00095 Ifpack_CrsRick::~Ifpack_CrsRick(){
00096 
00097 
00098   delete U_;
00099   delete D_; // Diagonal is stored separately.  We store the inverse.
00100 
00101   if (OverlapX_!=0) delete OverlapX_;
00102   if (OverlapY_!=0) delete OverlapY_;
00103 
00104   ValuesInitialized_ = false;
00105   Factored_ = false;
00106   Allocated_ = false;
00107 }
00108 
00109 #ifdef HAVE_IFPACK_TEUCHOS
00110 //==========================================================================
00111 int Ifpack_CrsRick::SetParameters(const Teuchos::ParameterList& parameterlist,
00112                   bool cerr_warning_if_unused)
00113 {
00114   Ifpack::param_struct params;
00115   params.double_params[Ifpack::relax_value] = RelaxValue_;
00116   params.double_params[Ifpack::absolute_threshold] = Athresh_;
00117   params.double_params[Ifpack::relative_threshold] = Rthresh_;
00118   params.overlap_mode = OverlapMode_;
00119 
00120   Ifpack::set_parameters(parameterlist, params, cerr_warning_if_unused);
00121 
00122   RelaxValue_ = params.double_params[Ifpack::relax_value];
00123   Athresh_ = params.double_params[Ifpack::absolute_threshold];
00124   Rthresh_ = params.double_params[Ifpack::relative_threshold];
00125   OverlapMode_ = params.overlap_mode;
00126 
00127   return(0);
00128 }
00129 #endif
00130 
00131 //==========================================================================
00132 int Ifpack_CrsRick::InitValues() {
00133 
00134   // if (!Allocated()) return(-1); // This test is not needed at this time.  All constructors allocate.
00135 
00136   int ierr = 0;
00137   int i, j;
00138   int * InI=0, * LI=0, * UI = 0;
00139   double * InV=0, * LV=0, * UV = 0;
00140   int NumIn, NumL, NumU;
00141   bool DiagFound;
00142   int NumNonzeroDiags = 0;
00143 
00144   Epetra_CrsMatrix * OverlapA = (Epetra_CrsMatrix *) &A_;
00145 
00146   if (Graph_.LevelOverlap()>0 && Graph_.L_Graph().DomainMap().DistributedGlobal()) {
00147   
00148   OverlapA = new Epetra_CrsMatrix(Copy, *Graph_.OverlapGraph());
00149   OverlapA->Import(A_, *Graph_.OverlapImporter(), Insert);
00150   }
00151 
00152   // Get Maximun Row length
00153   int MaxNumEntries = OverlapA->MaxNumEntries();
00154 
00155   InI = new int[MaxNumEntries]; // Allocate temp space
00156   UI = new int[MaxNumEntries];
00157   InV = new double[MaxNumEntries];
00158   UV = new double[MaxNumEntries];
00159 
00160   double *DV;
00161   ierr = D_->ExtractView(&DV); // Get view of diagonal
00162     
00163 
00164   // First we copy the user's matrix into diagonal vector and U, regardless of fill level
00165 
00166   for (i=0; i< NumMyRows(); i++) {
00167 
00168     OverlapA->ExtractMyRowCopy(i, MaxNumEntries, NumIn, InV, InI); // Get Values and Indices
00169     
00170     // Split into L and U (we don't assume that indices are ordered).
00171     
00172     NumL = 0; 
00173     NumU = 0; 
00174     DiagFound = false;
00175     
00176     for (j=0; j< NumIn; j++) {
00177       int k = InI[j];
00178 
00179       if (k==i) {
00180     DiagFound = true;
00181     DV[i] += Rthresh_ * InV[j] + EPETRA_SGN(InV[j]) * Athresh_; // Store perturbed diagonal in Epetra_Vector D_
00182       }
00183 
00184       else if (k < 0) return(-1); // Out of range
00185       else if (k<NumMyRows()) {
00186     UI[NumU] = k;
00187     UV[NumU] = InV[j];
00188     NumU++;
00189       }
00190     }
00191     
00192     // Check in things for this row of L and U
00193 
00194     if (DiagFound) NumNonzeroDiags++;
00195     if (NumU) U_->ReplaceMyValues(i, NumU, UV, UI);
00196     
00197   }
00198 
00199   delete [] UI;
00200   delete [] UV;
00201   delete [] InI;
00202   delete [] InV;
00203 
00204   if (Graph_.LevelOverlap()>0 && Graph_.L_Graph().DomainMap().DistributedGlobal()) delete OverlapA;
00205 
00206 
00207   U_->FillComplete();
00208 
00209   // At this point L and U have the values of A in the structure of L and U, and diagonal vector D
00210 
00211   SetValuesInitialized(true);
00212   SetFactored(false);
00213 
00214   int TotalNonzeroDiags = 0;
00215   Graph_.U_Graph().RowMap().Comm().SumAll(&NumNonzeroDiags, &TotalNonzeroDiags, 1);
00216   if (Graph_.LevelOverlap()==0 &&
00217       ((TotalNonzeroDiags!=NumGlobalRows()) || 
00218        (TotalNonzeroDiags!=NumGlobalDiagonals()))) ierr = 1;
00219   if (NumNonzeroDiags != NumMyDiagonals()) ierr = 1; // Diagonals are not right, warn user
00220 
00221   return(ierr);
00222 }
00223 
00224 //==========================================================================
00225 int Ifpack_CrsRick::Factor() {
00226 
00227   // if (!Allocated()) return(-1); // This test is not needed at this time.  All constructors allocate.
00228   if (!ValuesInitialized()) return(-2); // Must have values initialized.
00229   if (Factored()) return(-3); // Can't have already computed factors.
00230 
00231   SetValuesInitialized(false);
00232 
00233   // MinMachNum should be officially defined, for now pick something a little 
00234   // bigger than IEEE underflow value
00235 
00236   double MinDiagonalValue = Epetra_MinDouble;
00237   double MaxDiagonalValue = 1.0/MinDiagonalValue;
00238 
00239   int ierr = 0;
00240   int i, j, k;
00241   int * UI = 0;
00242   double * UV = 0;
00243   int NumIn, NumU;
00244 
00245   // Get Maximun Row length
00246   int MaxNumEntries = U_->MaxNumEntries() + 1;
00247 
00248   int * InI = new int[MaxNumEntries]; // Allocate temp space
00249   double * InV = new double[MaxNumEntries];
00250   int * colflag = new int[NumMyCols()];
00251 
00252   double *DV;
00253   ierr = D_->ExtractView(&DV); // Get view of diagonal
00254 
00255   int current_madds = 0; // We will count multiply-add as they happen
00256 
00257   // Now start the factorization.
00258 
00259   // Need some integer workspace and pointers
00260   int NumUU; 
00261   int * UUI;
00262   double * UUV;
00263   for (j=0; j<NumMyCols(); j++) colflag[j] = - 1;
00264 
00265   for(i=0; i<NumMyRows(); i++) {
00266 
00267  // Fill InV, InI with current row of L, D and U combined
00268 
00269     NumIn = MaxNumEntries;
00270     IFPACK_CHK_ERR(L_->ExtractMyRowCopy(i, NumIn, NumL, InV, InI)==0);
00271     LV = InV;
00272     LI = InI;
00273 
00274     InV[NumL] = DV[i]; // Put in diagonal
00275     InI[NumL] = i;
00276     
00277     IFPACK_CHK_ERR(U_->ExtractMyRowCopy(i, NumIn-NumL-1, NumU, InV+NumL+1, InI+NumL+1));
00278     NumIn = NumL+NumU+1;
00279     UV = InV+NumL+1;
00280     UI = InI+NumL+1;
00281 
00282     // Set column flags
00283     for (j=0; j<NumIn; j++) colflag[InI[j]] = j;
00284 
00285     double diagmod = 0.0; // Off-diagonal accumulator
00286 
00287     for (int jj=0; jj<NumL; jj++) {
00288       j = InI[jj];
00289       double multiplier = InV[jj]; // current_mults++;
00290 
00291       InV[jj] *= DV[j];
00292       
00293       IFPACK_CHK_ERR(U_->ExtractMyRowView(j, NumUU, UUV, UUI)); // View of row above
00294 
00295       if (RelaxValue_==0.0) {
00296     for (k=0; k<NumUU; k++) {
00297       int kk = colflag[UUI[k]];
00298       if (kk>-1) {
00299         InV[kk] -= multiplier*UUV[k];
00300         current_madds++;
00301       }
00302     }
00303       }
00304       else {
00305     for (k=0; k<NumUU; k++) {
00306       int kk = colflag[UUI[k]];
00307       if (kk>-1) InV[kk] -= multiplier*UUV[k];
00308       else diagmod -= multiplier*UUV[k];
00309       current_madds++;
00310     }
00311       }
00312      }
00313     if (NumL) 
00314       IFPACK_CHK_ERR(L_->ReplaceMyValues(i, NumL, LV, LI));  // Replace current row of L
00315 
00316     DV[i] = InV[NumL]; // Extract Diagonal value
00317 
00318     if (RelaxValue_!=0.0) {
00319       DV[i] += RelaxValue_*diagmod; // Add off diagonal modifications
00320       // current_madds++;
00321     }
00322 
00323     if (fabs(DV[i]) > MaxDiagonalValue) {
00324       if (DV[i] < 0) DV[i] = - MinDiagonalValue;
00325       else DV[i] = MinDiagonalValue;
00326     }
00327     else
00328       DV[i] = 1.0/DV[i]; // Invert diagonal value
00329 
00330     for (j=0; j<NumU; j++) UV[j] *= DV[i]; // Scale U by inverse of diagonal
00331 
00332     if (NumU) 
00333       IFPACK_CHK_ERR(U_->ReplaceMyValues(i, NumU, UV, UI));  // Replace current row of L and U
00334 
00335 
00336     // Reset column flags
00337     for (j=0; j<NumIn; j++) colflag[InI[j]] = -1;
00338   }
00339 
00340   
00341   // Add up flops
00342  
00343   double current_flops = 2 * current_madds;
00344   double total_flops = 0;
00345     
00346   Graph_.L_Graph().RowMap().Comm().SumAll(&current_flops, &total_flops, 1); // Get total madds across all PEs
00347 
00348   // Now count the rest
00349   total_flops += (double) L_->NumGlobalNonzeros(); // Accounts for multiplier above
00350   total_flops += (double) D_->GlobalLength(); // Accounts for reciprocal of diagonal
00351   if (RelaxValue_!=0.0) total_flops += 2 * (double)D_->GlobalLength(); // Accounts for relax update of diag
00352 
00353   UpdateFlops(total_flops); // Update flop count
00354 
00355   delete [] InI;
00356   delete [] InV;
00357   delete [] colflag;
00358   
00359   SetFactored(true);
00360 
00361   return(ierr);
00362 
00363 }
00364 
00365 //=============================================================================
00366 int Ifpack_CrsRick::Solve(bool Trans, const Epetra_Vector& X, 
00367                 Epetra_Vector& Y) const {
00368 //
00369 // This function finds Y such that LDU Y = X or U(trans) D L(trans) Y = X for a single RHS
00370 //
00371 
00372   bool Upper = true;
00373   bool Lower = false;
00374   bool UnitDiagonal = true;
00375 
00376   Epetra_Vector * X1 = (Epetra_Vector *) &X;
00377   Epetra_Vector * Y1 = (Epetra_Vector *) &Y;
00378 
00379   if (Graph_.LevelOverlap()>0 && Graph_.L_Graph().DomainMap().DistributedGlobal()) {
00380     if (OverlapX_==0) { // Need to allocate space for overlap X and Y
00381       OverlapX_ = new Epetra_Vector(Graph_.OverlapGraph()->RowMap());
00382       OverlapY_ = new Epetra_Vector(Graph_.OverlapGraph()->RowMap());
00383     }
00384     OverlapX_->Import(X,*Graph_.OverlapImporter(), Insert); // Import X values for solve
00385     X1 = (Epetra_Vector *) OverlapX_;
00386     Y1 = (Epetra_Vector *) OverlapY_; // Set pointers for X1 and Y1 to point to overlap space
00387   }
00388 
00389   Epetra_Flops * counter = this->GetFlopCounter();
00390   if (counter!=0) {
00391     L_->SetFlopCounter(*counter);
00392     Y1->SetFlopCounter(*counter);
00393     U_->SetFlopCounter(*counter);
00394   }
00395 
00396   if (!Trans) {
00397 
00398     L_->Solve(Lower, Trans, UnitDiagonal, *X1, *Y1);
00399     Y1->Multiply(1.0, *D_, *Y1, 0.0); // y = D*y (D_ has inverse of diagonal)
00400     U_->Solve(Upper, Trans, UnitDiagonal, *Y1, *Y1); // Solve Uy = y
00401   }
00402   else
00403     {
00404       U_->Solve(Upper, Trans, UnitDiagonal, *X1, *Y1); // Solve Uy = y
00405       Y1->Multiply(1.0, *D_, *Y1, 0.0); // y = D*y (D_ has inverse of diagonal)
00406       L_->Solve(Lower, Trans, UnitDiagonal, *Y1, *Y1);
00407       
00408     } 
00409 
00410   // Export computed Y values as directed
00411   if (Graph_.LevelOverlap()>0 && Graph_.L_Graph().DomainMap().DistributedGlobal())
00412     Y.Export(*OverlapY_,*Graph_.OverlapImporter(), OverlapMode_);
00413   return(0);
00414 }
00415 
00416 
00417 //=============================================================================
00418 int Ifpack_CrsRick::Solve(bool Trans, const Epetra_MultiVector& X, 
00419                 Epetra_MultiVector& Y) const {
00420 //
00421 // This function finds Y such that LDU Y = X or U(trans) D L(trans) Y = X for multiple RHS
00422 //
00423 
00424   if (X.NumVectors()!=Y.NumVectors()) EPETRA_CHK_ERR(-1); // Return error: X and Y not the same size
00425 
00426   bool Upper = true;
00427   bool Lower = false;
00428   bool UnitDiagonal = true;
00429 
00430   Epetra_MultiVector * X1 = (Epetra_MultiVector *) &X;
00431   Epetra_MultiVector * Y1 = (Epetra_MultiVector *) &Y;
00432 
00433   if (Graph_.LevelOverlap()>0 && Graph_.L_Graph().DomainMap().DistributedGlobal()) {
00434     // Make sure the number of vectors in the multivector is the same as before.
00435     if (OverlapX_!=0) {
00436       if (OverlapX_->NumVectors()!=X.NumVectors()) {
00437     delete OverlapX_; OverlapX_ = 0;
00438     delete OverlapY_; OverlapY_ = 0;
00439       }
00440     }
00441     if (OverlapX_==0) { // Need to allocate space for overlap X and Y
00442       OverlapX_ = new Epetra_MultiVector(Graph_.OverlapGraph()->RowMap(), X.NumVectors());
00443       OverlapY_ = new Epetra_MultiVector(Graph_.OverlapGraph()->RowMap(), Y.NumVectors());
00444     }
00445     OverlapX_->Import(X,*Graph_.OverlapImporter(), Insert); // Import X values for solve
00446     X1 = OverlapX_;
00447     Y1 = OverlapY_; // Set pointers for X1 and Y1 to point to overlap space
00448   }
00449 
00450   Epetra_Flops * counter = this->GetFlopCounter();
00451   if (counter!=0) {
00452     L_->SetFlopCounter(*counter);
00453     Y1->SetFlopCounter(*counter);
00454     U_->SetFlopCounter(*counter);
00455   }
00456 
00457   if (!Trans) {
00458 
00459     L_->Solve(Lower, Trans, UnitDiagonal, *X1, *Y1);
00460     Y1->Multiply(1.0, *D_, *Y1, 0.0); // y = D*y (D_ has inverse of diagonal)
00461     U_->Solve(Upper, Trans, UnitDiagonal, *Y1, *Y1); // Solve Uy = y
00462   }
00463   else
00464     {
00465       U_->Solve(Upper, Trans, UnitDiagonal, *X1, *Y1); // Solve Uy = y
00466       Y1->Multiply(1.0, *D_, *Y1, 0.0); // y = D*y (D_ has inverse of diagonal)
00467       L_->Solve(Lower, Trans, UnitDiagonal, *Y1, *Y1);
00468       
00469     } 
00470 
00471   // Export computed Y values as directed
00472   if (Graph_.LevelOverlap()>0 && Graph_.L_Graph().DomainMap().DistributedGlobal())
00473     Y.Export(*OverlapY_,*Graph_.OverlapImporter(), OverlapMode_);
00474   return(0);
00475 }
00476 //=============================================================================
00477 int Ifpack_CrsRick::Multiply(bool Trans, const Epetra_MultiVector& X, 
00478                 Epetra_MultiVector& Y) const {
00479 //
00480 // This function finds X such that LDU Y = X or U(trans) D L(trans) Y = X for multiple RHS
00481 //
00482 
00483   if (X.NumVectors()!=Y.NumVectors()) EPETRA_CHK_ERR(-1); // Return error: X and Y not the same size
00484 
00485   bool Upper = true;
00486   bool Lower = false;
00487   bool UnitDiagonal = true;
00488 
00489   Epetra_MultiVector * X1 = (Epetra_MultiVector *) &X;
00490   Epetra_MultiVector * Y1 = (Epetra_MultiVector *) &Y;
00491 
00492   if (Graph_.LevelOverlap()>0 && Graph_.L_Graph().DomainMap().DistributedGlobal()) {
00493     // Make sure the number of vectors in the multivector is the same as before.
00494     if (OverlapX_!=0) {
00495       if (OverlapX_->NumVectors()!=X.NumVectors()) {
00496     delete OverlapX_; OverlapX_ = 0;
00497     delete OverlapY_; OverlapY_ = 0;
00498       }
00499     }
00500     if (OverlapX_==0) { // Need to allocate space for overlap X and Y
00501       OverlapX_ = new Epetra_MultiVector(Graph_.OverlapGraph()->RowMap(), X.NumVectors());
00502       OverlapY_ = new Epetra_MultiVector(Graph_.OverlapGraph()->RowMap(), Y.NumVectors());
00503     }
00504     OverlapX_->Import(X,*Graph_.OverlapImporter(), Insert); // Import X values for solve
00505     X1 = OverlapX_;
00506     Y1 = OverlapY_; // Set pointers for X1 and Y1 to point to overlap space
00507   }
00508 
00509   Epetra_Flops * counter = this->GetFlopCounter();
00510   if (counter!=0) {
00511     L_->SetFlopCounter(*counter);
00512     Y1->SetFlopCounter(*counter);
00513     U_->SetFlopCounter(*counter);
00514   }
00515 
00516   if (Trans) {
00517 
00518     L_->Multiply(Trans, *X1, *Y1);
00519     Y1->Update(1.0, *X1, 1.0); // Y1 = Y1 + X1 (account for implicit unit diagonal)
00520     Y1->ReciprocalMultiply(1.0, *D_, *Y1, 0.0); // y = D*y (D_ has inverse of diagonal)
00521     Epetra_MultiVector Y1temp(*Y1); // Need a temp copy of Y1
00522     U_->Multiply(Trans, Y1temp, *Y1);
00523     Y1->Update(1.0, Y1temp, 1.0); // (account for implicit unit diagonal)
00524   }
00525   else {
00526     U_->Multiply(Trans, *X1, *Y1); // 
00527     Y1->Update(1.0, *X1, 1.0); // Y1 = Y1 + X1 (account for implicit unit diagonal)
00528     Y1->ReciprocalMultiply(1.0, *D_, *Y1, 0.0); // y = D*y (D_ has inverse of diagonal)
00529     Epetra_MultiVector Y1temp(*Y1); // Need a temp copy of Y1
00530     L_->Multiply(Trans, Y1temp, *Y1);
00531     Y1->Update(1.0, Y1temp, 1.0); // (account for implicit unit diagonal)
00532     } 
00533 
00534   // Export computed Y values as directed
00535   if (Graph_.LevelOverlap()>0 && Graph_.L_Graph().DomainMap().DistributedGlobal())
00536     Y.Export(*OverlapY_,*Graph_.OverlapImporter(), OverlapMode_);
00537   return(0);
00538 }
00539 //=============================================================================
00540 int Ifpack_CrsRick::Condest(bool Trans, double & ConditionNumberEstimate) const {
00541 
00542   if (Condest_>=0.0) {
00543     ConditionNumberEstimate = Condest_;
00544     return(0);
00545   }
00546   // Create a vector with all values equal to one
00547   Epetra_Vector Ones(A_.RowMap());
00548   Epetra_Vector OnesResult(Ones);
00549   Ones.PutScalar(1.0);
00550 
00551   EPETRA_CHK_ERR(Solve(Trans, Ones, OnesResult)); // Compute the effect of the solve on the vector of ones
00552   EPETRA_CHK_ERR(OnesResult.Abs(OnesResult)); // Make all values non-negative
00553   EPETRA_CHK_ERR(OnesResult.MaxValue(&ConditionNumberEstimate)); // Get the maximum value across all processors
00554   Condest_ = ConditionNumberEstimate; // Save value for possible later calls
00555   return(0);
00556 }
00557 //=============================================================================
00558 // Non-member functions
00559 
00560 ostream& operator << (ostream& os, const Ifpack_CrsRick& A)
00561 {
00562 /*  Epetra_fmtflags olda = os.setf(ios::right,ios::adjustfield);
00563   Epetra_fmtflags oldf = os.setf(ios::scientific,ios::floatfield);
00564   int oldp = os.precision(12); */
00565   int LevelFill = A.Graph().LevelFill();
00566   int LevelOverlap = A.Graph().LevelOverlap();
00567   Epetra_CrsMatrix & L = (Epetra_CrsMatrix &) A.L();
00568   Epetra_CrsMatrix & U = (Epetra_CrsMatrix &) A.U();
00569   Epetra_Vector & D = (Epetra_Vector &) A.D();
00570 
00571   os.width(14);
00572   os << endl;
00573   os <<  "     Level of Fill = "; os << LevelFill;
00574   os << endl;
00575   os.width(14);
00576   os <<  "     Level of Overlap = "; os << LevelOverlap;
00577   os << endl;
00578 
00579   os.width(14);
00580   os <<  "     Lower Triangle = ";
00581   os << endl;
00582   os << L; // Let Epetra_CrsMatrix handle the rest.
00583   os << endl;
00584 
00585   os.width(14);
00586   os <<  "     Inverse of Diagonal = ";
00587   os << endl;
00588   os << D; // Let Epetra_Vector handle the rest.
00589   os << endl;
00590 
00591   os.width(14);
00592   os <<  "     Upper Triangle = ";
00593   os << endl;
00594   os << U; // Let Epetra_CrsMatrix handle the rest.
00595   os << endl;
00596  
00597   // Reset os flags
00598 
00599 /*  os.setf(olda,ios::adjustfield);
00600   os.setf(oldf,ios::floatfield);
00601   os.precision(oldp); */
00602 
00603   return os;
00604 }

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