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