Ifpack Package Browser (Single Doxygen Collection) Development
Ifpack_CrsRiluk.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_CrsRiluk.h"
00030 #include "Epetra_Comm.h"
00031 #include "Epetra_Map.h"
00032 #include "Epetra_CrsGraph.h"
00033 #include "Epetra_VbrMatrix.h"
00034 #include "Epetra_RowMatrix.h"
00035 #include "Epetra_MultiVector.h"
00036 
00037 #include <Teuchos_ParameterList.hpp>
00038 #include <ifp_parameters.h>
00039 
00040 //==============================================================================
00041 Ifpack_CrsRiluk::Ifpack_CrsRiluk(const Ifpack_IlukGraph & Graph_in) 
00042   : UserMatrixIsVbr_(false),
00043     UserMatrixIsCrs_(false),
00044     Graph_(Graph_in),
00045     Comm_(Graph_in.Comm()),
00046     UseTranspose_(false),
00047     NumMyDiagonals_(0),
00048     Allocated_(false),
00049     ValuesInitialized_(false),
00050     Factored_(false),
00051     RelaxValue_(0.0),
00052     Athresh_(0.0),
00053     Rthresh_(1.0),
00054     Condest_(-1.0),
00055     OverlapMode_(Zero)
00056 {
00057   // Test for non-trivial overlap here so we can use it later.
00058   IsOverlapped_ = (Graph_in.LevelOverlap()>0 && Graph_in.DomainMap().DistributedGlobal());
00059 }
00060 
00061 //==============================================================================
00062 Ifpack_CrsRiluk::Ifpack_CrsRiluk(const Ifpack_CrsRiluk & FactoredMatrix) 
00063   : UserMatrixIsVbr_(FactoredMatrix.UserMatrixIsVbr_),
00064     UserMatrixIsCrs_(FactoredMatrix.UserMatrixIsCrs_),
00065     IsOverlapped_(FactoredMatrix.IsOverlapped_),
00066     Graph_(FactoredMatrix.Graph_),
00067     IlukRowMap_(FactoredMatrix.IlukRowMap_),
00068     IlukDomainMap_(FactoredMatrix.IlukDomainMap_),
00069     IlukRangeMap_(FactoredMatrix.IlukRangeMap_),
00070     Comm_(FactoredMatrix.Comm_),
00071     UseTranspose_(FactoredMatrix.UseTranspose_),
00072     NumMyDiagonals_(FactoredMatrix.NumMyDiagonals_),
00073     Allocated_(FactoredMatrix.Allocated_),
00074     ValuesInitialized_(FactoredMatrix.ValuesInitialized_),
00075     Factored_(FactoredMatrix.Factored_),
00076     RelaxValue_(FactoredMatrix.RelaxValue_),
00077     Athresh_(FactoredMatrix.Athresh_),
00078     Rthresh_(FactoredMatrix.Rthresh_),
00079     Condest_(FactoredMatrix.Condest_),
00080     OverlapMode_(FactoredMatrix.OverlapMode_)
00081 {
00082   L_ = Teuchos::rcp( new Epetra_CrsMatrix(FactoredMatrix.L()) );
00083   U_ = Teuchos::rcp( new Epetra_CrsMatrix(FactoredMatrix.U()) );
00084   D_ = Teuchos::rcp( new Epetra_Vector(FactoredMatrix.D()) );
00085   if (IlukRowMap_!=Teuchos::null) IlukRowMap_ = Teuchos::rcp( new Epetra_Map(*IlukRowMap_) );
00086   if (IlukDomainMap_!=Teuchos::null) IlukDomainMap_ = Teuchos::rcp( new Epetra_Map(*IlukDomainMap_) );
00087   if (IlukRangeMap_!=Teuchos::null) IlukRangeMap_ = Teuchos::rcp( new Epetra_Map(*IlukRangeMap_) );
00088   
00089 }
00090 //==============================================================================
00091 Ifpack_CrsRiluk::~Ifpack_CrsRiluk(){
00092 
00093   ValuesInitialized_ = false;
00094   Factored_ = false;
00095   Allocated_ = false;
00096 }
00097 //==============================================================================
00098 int Ifpack_CrsRiluk::AllocateCrs() {
00099 
00100   // Allocate Epetra_CrsMatrix using ILUK graphs
00101   L_ = Teuchos::rcp( new Epetra_CrsMatrix(Copy, Graph_.L_Graph()) );
00102   U_ = Teuchos::rcp( new Epetra_CrsMatrix(Copy, Graph_.U_Graph()) );
00103   D_ = Teuchos::rcp( new Epetra_Vector(Graph_.L_Graph().RowMap()) );
00104   L_Graph_ = Teuchos::null;
00105   U_Graph_ = Teuchos::null;
00106   SetAllocated(true);
00107   return(0);
00108 }
00109 //==============================================================================
00110 int Ifpack_CrsRiluk::AllocateVbr() {
00111 
00112   // First we need to create a set of Epetra_Maps that has the same number of points  as the
00113   // Epetra_BlockMaps associated with the Overlap Graph.
00114   EPETRA_CHK_ERR(BlockMap2PointMap(Graph_.L_Graph().RowMap(), &IlukRowMap_));
00115   EPETRA_CHK_ERR(BlockMap2PointMap(Graph_.U_Graph().DomainMap(), &IlukDomainMap_));
00116   EPETRA_CHK_ERR(BlockMap2PointMap(Graph_.L_Graph().RangeMap(), &IlukRangeMap_));
00117 
00118   // Set L range map and U domain map
00119   U_DomainMap_ = IlukDomainMap_;
00120   L_RangeMap_ = IlukRangeMap_;
00121   // If there is fill, then pre-build the L and U structures from the Block version of L and U.
00122   if (Graph().LevelFill()) { 
00123     L_Graph_ = Teuchos::rcp( new Epetra_CrsGraph(Copy, *IlukRowMap_, *IlukRowMap_, 0) );
00124     U_Graph_ = Teuchos::rcp( new Epetra_CrsGraph(Copy, *IlukRowMap_, *IlukRowMap_, 0) );
00125     EPETRA_CHK_ERR(BlockGraph2PointGraph(Graph_.L_Graph(), *L_Graph_, false));
00126     EPETRA_CHK_ERR(BlockGraph2PointGraph(Graph_.U_Graph(), *U_Graph_, true));
00127     
00128     L_Graph_->FillComplete(*IlukRowMap_, *IlukRangeMap_);
00129     U_Graph_->FillComplete(*IlukDomainMap_, *IlukRowMap_);
00130 
00131     L_ = Teuchos::rcp( new Epetra_CrsMatrix(Copy, *L_Graph_) );
00132     U_ = Teuchos::rcp( new Epetra_CrsMatrix(Copy, *U_Graph_) );
00133     D_ = Teuchos::rcp( new Epetra_Vector(*IlukRowMap_) );
00134   }
00135   else {
00136     // Allocate Epetra_CrsMatrix using ILUK graphs
00137     L_ = Teuchos::rcp( new Epetra_CrsMatrix(Copy, *IlukRowMap_, *IlukRowMap_, 0) );
00138     U_ = Teuchos::rcp( new Epetra_CrsMatrix(Copy, *IlukRowMap_, *IlukRowMap_, 0) );
00139     D_ = Teuchos::rcp( new Epetra_Vector(*IlukRowMap_) );
00140     L_Graph_ = Teuchos::null;
00141     U_Graph_ = Teuchos::null;
00142   }
00143   SetAllocated(true);
00144   return(0);
00145 }
00146 
00147 //==========================================================================
00148 int Ifpack_CrsRiluk::SetParameters(const Teuchos::ParameterList& parameterlist,
00149                                    bool cerr_warning_if_unused)
00150 {
00151   Ifpack::param_struct params;
00152   params.double_params[Ifpack::relax_value] = RelaxValue_;
00153   params.double_params[Ifpack::absolute_threshold] = Athresh_;
00154   params.double_params[Ifpack::relative_threshold] = Rthresh_;
00155   params.overlap_mode = OverlapMode_;
00156 
00157   Ifpack::set_parameters(parameterlist, params, cerr_warning_if_unused);
00158 
00159   RelaxValue_ = params.double_params[Ifpack::relax_value];
00160   Athresh_ = params.double_params[Ifpack::absolute_threshold];
00161   Rthresh_ = params.double_params[Ifpack::relative_threshold];
00162   OverlapMode_ = params.overlap_mode;
00163 
00164   return(0);
00165 }
00166 
00167 //==========================================================================
00168 int Ifpack_CrsRiluk::InitValues(const Epetra_CrsMatrix & A) {
00169 
00170   UserMatrixIsCrs_ = true;
00171 
00172   if (!Allocated()) AllocateCrs();
00173 
00174   Teuchos::RefCountPtr<Epetra_CrsMatrix> OverlapA = Teuchos::rcp( (Epetra_CrsMatrix *) &A, false );
00175 
00176   if (IsOverlapped_) {
00177   
00178     OverlapA = Teuchos::rcp( new Epetra_CrsMatrix(Copy, *Graph_.OverlapGraph()) );
00179     EPETRA_CHK_ERR(OverlapA->Import(A, *Graph_.OverlapImporter(), Insert));
00180     EPETRA_CHK_ERR(OverlapA->FillComplete());
00181   }
00182   
00183   // Get Maximun Row length
00184   int MaxNumEntries = OverlapA->MaxNumEntries();
00185 
00186   // Set L range map and U domain map
00187   U_DomainMap_ = Teuchos::rcp( &(A.DomainMap()), false );
00188   L_RangeMap_ = Teuchos::rcp( &(A.RangeMap()), false );
00189   // Do the rest using generic Epetra_RowMatrix interface
00190 
00191   EPETRA_CHK_ERR(InitAllValues(*OverlapA, MaxNumEntries));
00192 
00193   return(0);
00194 }
00195 
00196 //==========================================================================
00197 int Ifpack_CrsRiluk::InitValues(const Epetra_VbrMatrix & A) {
00198 
00199   UserMatrixIsVbr_ = true;
00200 
00201   if (!Allocated()) AllocateVbr();
00202 
00203   //cout << "Original Graph " << endl <<  A.Graph() << endl << flush;
00204   //A.Comm().Barrier(); 
00205   //if (A.Comm().MyPID()==0) cout << "*****************************************************" <<endl;
00206   //cout << "Original Matrix " << endl << A << endl << flush;
00207   //A.Comm().Barrier(); 
00208   //if (A.Comm().MyPID()==0) cout << "*****************************************************" <<endl;
00209   //cout << "Overlap Graph " << endl << *Graph_.OverlapGraph() << endl << flush;
00210   //A.Comm().Barrier(); 
00211   //if (A.Comm().MyPID()==0) cout << "*****************************************************" <<endl;
00212 
00213   Teuchos::RefCountPtr<Epetra_VbrMatrix> OverlapA = Teuchos::rcp( (Epetra_VbrMatrix *) &A, false );
00214 
00215   if (IsOverlapped_) {
00216   
00217     OverlapA = Teuchos::rcp( new Epetra_VbrMatrix(Copy, *Graph_.OverlapGraph()) );
00218     EPETRA_CHK_ERR(OverlapA->Import(A, *Graph_.OverlapImporter(), Insert));
00219     EPETRA_CHK_ERR(OverlapA->FillComplete());
00220   }
00221   
00222   //cout << "Overlap Matrix " << endl << *OverlapA << endl << flush;
00223 
00224   // Get Maximun Row length
00225   int MaxNumEntries = OverlapA->MaxNumNonzeros();
00226 
00227   // Do the rest using generic Epetra_RowMatrix interface
00228 
00229   EPETRA_CHK_ERR(InitAllValues(*OverlapA, MaxNumEntries));
00230 
00231   return(0);
00232 }
00233 //==========================================================================
00234 
00235 int Ifpack_CrsRiluk::InitAllValues(const Epetra_RowMatrix & OverlapA, int MaxNumEntries) {
00236 
00237   int ierr = 0;
00238   int i, j;
00239   int NumIn, NumL, NumU;
00240   bool DiagFound;
00241   int NumNonzeroDiags = 0;
00242 
00243 
00244   vector<int> InI(MaxNumEntries); // Allocate temp space
00245   vector<int> LI(MaxNumEntries);
00246   vector<int> UI(MaxNumEntries);
00247   vector<double> InV(MaxNumEntries);
00248   vector<double> LV(MaxNumEntries);
00249   vector<double> UV(MaxNumEntries);
00250 
00251   bool ReplaceValues = (L_->StaticGraph() || L_->IndicesAreLocal()); // Check if values should be inserted or replaced
00252 
00253   if (ReplaceValues) {
00254     L_->PutScalar(0.0); // Zero out L and U matrices
00255     U_->PutScalar(0.0);
00256   }
00257 
00258   D_->PutScalar(0.0); // Set diagonal values to zero
00259   double *DV;
00260   EPETRA_CHK_ERR(D_->ExtractView(&DV)); // Get view of diagonal
00261     
00262 
00263   // First we copy the user's matrix into L and U, regardless of fill level
00264 
00265   for (i=0; i< NumMyRows(); i++) {
00266 
00267     EPETRA_CHK_ERR(OverlapA.ExtractMyRowCopy(i, MaxNumEntries, NumIn, &InV[0], &InI[0])); // Get Values and Indices
00268     
00269     // Split into L and U (we don't assume that indices are ordered).
00270     
00271     NumL = 0; 
00272     NumU = 0; 
00273     DiagFound = false;
00274     
00275     for (j=0; j< NumIn; j++) {
00276       int k = InI[j];
00277 
00278       if (k==i) {
00279   DiagFound = true;
00280   DV[i] += Rthresh_ * InV[j] + EPETRA_SGN(InV[j]) * Athresh_; // Store perturbed diagonal in Epetra_Vector D_
00281       }
00282 
00283       else if (k < 0) {EPETRA_CHK_ERR(-1);} // Out of range
00284 
00285       else if (k < i) {
00286   LI[NumL] = k;
00287   LV[NumL] = InV[j];
00288   NumL++;
00289       }
00290       else if (k<NumMyRows()) {
00291   UI[NumU] = k;
00292   UV[NumU] = InV[j];
00293   NumU++;
00294       }
00295     }
00296     
00297     // Check in things for this row of L and U
00298 
00299     if (DiagFound) NumNonzeroDiags++;
00300     else DV[i] = Athresh_;
00301 
00302     if (NumL) {
00303       if (ReplaceValues) {
00304   EPETRA_CHK_ERR(L_->ReplaceMyValues(i, NumL, &LV[0], &LI[0]));
00305       }
00306       else {
00307   EPETRA_CHK_ERR(L_->InsertMyValues(i, NumL, &LV[0], &LI[0]));
00308       }
00309     }
00310 
00311     if (NumU) {
00312       if (ReplaceValues) {
00313   EPETRA_CHK_ERR(U_->ReplaceMyValues(i, NumU, &UV[0], &UI[0]));
00314       }
00315       else {
00316   EPETRA_CHK_ERR(U_->InsertMyValues(i, NumU, &UV[0], &UI[0]));
00317       }
00318     }
00319     
00320   }
00321 
00322   if (!ReplaceValues) {
00323     // The domain of L and the range of U are exactly their own row maps (there is no communication).
00324     // The domain of U and the range of L must be the same as those of the original matrix,
00325     // However if the original matrix is a VbrMatrix, these two latter maps are translation from
00326     // a block map to a point map.
00327     EPETRA_CHK_ERR(L_->FillComplete(L_->RowMatrixColMap(), *L_RangeMap_));
00328     EPETRA_CHK_ERR(U_->FillComplete(*U_DomainMap_, U_->RowMatrixRowMap()));
00329   }
00330 
00331   // At this point L and U have the values of A in the structure of L and U, and diagonal vector D
00332 
00333   SetValuesInitialized(true);
00334   SetFactored(false);
00335 
00336   int TotalNonzeroDiags = 0;
00337   EPETRA_CHK_ERR(Graph_.L_Graph().RowMap().Comm().SumAll(&NumNonzeroDiags, &TotalNonzeroDiags, 1));
00338   NumMyDiagonals_ = NumNonzeroDiags;
00339   if (NumNonzeroDiags != NumMyRows()) ierr = 1; // Diagonals are not right, warn user
00340 
00341   return(ierr);
00342 }
00343 
00344 //==========================================================================
00345 int Ifpack_CrsRiluk::Factor() {
00346 
00347   // if (!Allocated()) return(-1); // This test is not needed at this time.  All constructors allocate.
00348   if (!ValuesInitialized()) return(-2); // Must have values initialized.
00349   if (Factored()) return(-3); // Can't have already computed factors.
00350 
00351   SetValuesInitialized(false);
00352 
00353   // MinMachNum should be officially defined, for now pick something a little 
00354   // bigger than IEEE underflow value
00355 
00356   double MinDiagonalValue = Epetra_MinDouble;
00357   double MaxDiagonalValue = 1.0/MinDiagonalValue;
00358 
00359   int ierr = 0;
00360   int i, j, k;
00361   int * LI=0, * UI = 0;
00362   double * LV=0, * UV = 0;
00363   int NumIn, NumL, NumU;
00364 
00365   // Get Maximun Row length
00366   int MaxNumEntries = L_->MaxNumEntries() + U_->MaxNumEntries() + 1;
00367 
00368   vector<int> InI(MaxNumEntries); // Allocate temp space
00369   vector<double> InV(MaxNumEntries);
00370   vector<int> colflag(NumMyCols());
00371 
00372   double *DV;
00373   ierr = D_->ExtractView(&DV); // Get view of diagonal
00374 
00375 #ifdef IFPACK_FLOPCOUNTERS
00376   int current_madds = 0; // We will count multiply-add as they happen
00377 #endif
00378 
00379   // Now start the factorization.
00380 
00381   // Need some integer workspace and pointers
00382   int NumUU; 
00383   int * UUI;
00384   double * UUV;
00385   for (j=0; j<NumMyCols(); j++) colflag[j] = - 1;
00386 
00387   for(i=0; i<NumMyRows(); i++) {
00388 
00389  // Fill InV, InI with current row of L, D and U combined
00390 
00391     NumIn = MaxNumEntries;
00392     EPETRA_CHK_ERR(L_->ExtractMyRowCopy(i, NumIn, NumL, &InV[0], &InI[0]));
00393     LV = &InV[0];
00394     LI = &InI[0];
00395 
00396     InV[NumL] = DV[i]; // Put in diagonal
00397     InI[NumL] = i;
00398     
00399     EPETRA_CHK_ERR(U_->ExtractMyRowCopy(i, NumIn-NumL-1, NumU, &InV[NumL+1], &InI[NumL+1]));
00400     NumIn = NumL+NumU+1;
00401     UV = &InV[NumL+1];
00402     UI = &InI[NumL+1];
00403 
00404     // Set column flags
00405     for (j=0; j<NumIn; j++) colflag[InI[j]] = j;
00406 
00407     double diagmod = 0.0; // Off-diagonal accumulator
00408 
00409     for (int jj=0; jj<NumL; jj++) {
00410       j = InI[jj];
00411       double multiplier = InV[jj]; // current_mults++;
00412 
00413       InV[jj] *= DV[j];
00414       
00415       EPETRA_CHK_ERR(U_->ExtractMyRowView(j, NumUU, UUV, UUI)); // View of row above
00416 
00417       if (RelaxValue_==0.0) {
00418   for (k=0; k<NumUU; k++) {
00419     int kk = colflag[UUI[k]];
00420     if (kk>-1) {
00421       InV[kk] -= multiplier*UUV[k];
00422 #ifdef IFPACK_FLOPCOUNTERS
00423       current_madds++;
00424 #endif
00425     }
00426   }
00427       }
00428       else {
00429   for (k=0; k<NumUU; k++) {
00430     int kk = colflag[UUI[k]];
00431     if (kk>-1) InV[kk] -= multiplier*UUV[k];
00432     else diagmod -= multiplier*UUV[k];
00433 #ifdef IFPACK_FLOPCOUNTERS
00434     current_madds++;
00435 #endif
00436   }
00437       }
00438      }
00439     if (NumL) {
00440       EPETRA_CHK_ERR(L_->ReplaceMyValues(i, NumL, LV, LI));  // Replace current row of L
00441     }
00442 
00443     DV[i] = InV[NumL]; // Extract Diagonal value
00444 
00445     if (RelaxValue_!=0.0) {
00446       DV[i] += RelaxValue_*diagmod; // Add off diagonal modifications
00447       // current_madds++;
00448     }
00449 
00450     if (fabs(DV[i]) > MaxDiagonalValue) {
00451       if (DV[i] < 0) DV[i] = - MinDiagonalValue;
00452       else DV[i] = MinDiagonalValue;
00453     }
00454     else
00455       DV[i] = 1.0/DV[i]; // Invert diagonal value
00456 
00457     for (j=0; j<NumU; j++) UV[j] *= DV[i]; // Scale U by inverse of diagonal
00458 
00459     if (NumU) {
00460       EPETRA_CHK_ERR(U_->ReplaceMyValues(i, NumU, UV, UI));  // Replace current row of L and U
00461     }
00462 
00463     // Reset column flags
00464     for (j=0; j<NumIn; j++) colflag[InI[j]] = -1;
00465   }
00466 
00467   // Validate that the L and U factors are actually lower and upper triangular
00468 
00469   if( !L_->LowerTriangular() ) 
00470     EPETRA_CHK_ERR(-2);
00471   if( !U_->UpperTriangular() ) 
00472     EPETRA_CHK_ERR(-3);
00473   
00474 #ifdef IFPACK_FLOPCOUNTERS
00475   // Add up flops
00476  
00477   double current_flops = 2 * current_madds;
00478   double total_flops = 0;
00479     
00480   EPETRA_CHK_ERR(Graph_.L_Graph().RowMap().Comm().SumAll(&current_flops, &total_flops, 1)); // Get total madds across all PEs
00481 
00482   // Now count the rest
00483   total_flops += (double) L_->NumGlobalNonzeros(); // Accounts for multiplier above
00484   total_flops += (double) D_->GlobalLength(); // Accounts for reciprocal of diagonal
00485   if (RelaxValue_!=0.0) total_flops += 2 * (double)D_->GlobalLength(); // Accounts for relax update of diag
00486 
00487   UpdateFlops(total_flops); // Update flop count
00488 #endif
00489 
00490   SetFactored(true);
00491 
00492   return(ierr);
00493 
00494 }
00495 
00496 //=============================================================================
00497 int Ifpack_CrsRiluk::Solve(bool Trans, const Epetra_MultiVector& X, 
00498         Epetra_MultiVector& Y) const {
00499 //
00500 // This function finds Y such that LDU Y = X or U(trans) D L(trans) Y = X for multiple RHS
00501 //
00502 
00503   // First generate X and Y as needed for this function
00504   Teuchos::RefCountPtr<Epetra_MultiVector> X1;
00505   Teuchos::RefCountPtr<Epetra_MultiVector> Y1;
00506   EPETRA_CHK_ERR(GenerateXY(Trans, X, Y, &X1, &Y1));
00507 
00508   bool Upper = true;
00509   bool Lower = false;
00510   bool UnitDiagonal = true;
00511 
00512 #ifdef IFPACK_FLOPCOUNTERS
00513   Epetra_Flops * counter = this->GetFlopCounter();
00514   if (counter!=0) {
00515     L_->SetFlopCounter(*counter);
00516     Y1->SetFlopCounter(*counter);
00517     U_->SetFlopCounter(*counter);
00518   }
00519 #endif
00520 
00521   if (!Trans) {
00522 
00523     EPETRA_CHK_ERR(L_->Solve(Lower, Trans, UnitDiagonal, *X1, *Y1));
00524     EPETRA_CHK_ERR(Y1->Multiply(1.0, *D_, *Y1, 0.0)); // y = D*y (D_ has inverse of diagonal)
00525     EPETRA_CHK_ERR(U_->Solve(Upper, Trans, UnitDiagonal, *Y1, *Y1)); // Solve Uy = y
00526     if (IsOverlapped_) {EPETRA_CHK_ERR(Y.Export(*Y1,*L_->Exporter(), OverlapMode_));} // Export computed Y values if needed
00527   }
00528   else {
00529     EPETRA_CHK_ERR(U_->Solve(Upper, Trans, UnitDiagonal, *X1, *Y1)); // Solve Uy = y
00530     EPETRA_CHK_ERR(Y1->Multiply(1.0, *D_, *Y1, 0.0)); // y = D*y (D_ has inverse of diagonal)
00531     EPETRA_CHK_ERR(L_->Solve(Lower, Trans, UnitDiagonal, *Y1, *Y1));
00532     if (IsOverlapped_) {EPETRA_CHK_ERR(Y.Export(*Y1,*U_->Importer(), OverlapMode_));} // Export computed Y values if needed
00533   } 
00534 
00535   return(0);
00536 }
00537 //=============================================================================
00538 int Ifpack_CrsRiluk::Multiply(bool Trans, const Epetra_MultiVector& X, 
00539             Epetra_MultiVector& Y) const {
00540 //
00541 // This function finds X such that LDU Y = X or U(trans) D L(trans) Y = X for multiple RHS
00542 //
00543     
00544   // First generate X and Y as needed for this function
00545   Teuchos::RefCountPtr<Epetra_MultiVector> X1;
00546   Teuchos::RefCountPtr<Epetra_MultiVector> Y1;
00547   EPETRA_CHK_ERR(GenerateXY(Trans, X, Y, &X1, &Y1));
00548 
00549 #ifdef IFPACK_FLOPCOUNTERS
00550   Epetra_Flops * counter = this->GetFlopCounter();
00551   if (counter!=0) {
00552     L_->SetFlopCounter(*counter);
00553     Y1->SetFlopCounter(*counter);
00554     U_->SetFlopCounter(*counter);
00555   }
00556 #endif
00557 
00558   if (!Trans) {
00559     EPETRA_CHK_ERR(U_->Multiply(Trans, *X1, *Y1)); // 
00560     EPETRA_CHK_ERR(Y1->Update(1.0, *X1, 1.0)); // Y1 = Y1 + X1 (account for implicit unit diagonal)
00561     EPETRA_CHK_ERR(Y1->ReciprocalMultiply(1.0, *D_, *Y1, 0.0)); // y = D*y (D_ has inverse of diagonal)
00562     Epetra_MultiVector Y1temp(*Y1); // Need a temp copy of Y1
00563     EPETRA_CHK_ERR(L_->Multiply(Trans, Y1temp, *Y1));
00564     EPETRA_CHK_ERR(Y1->Update(1.0, Y1temp, 1.0)); // (account for implicit unit diagonal)
00565     if (IsOverlapped_) {EPETRA_CHK_ERR(Y.Export(*Y1,*L_->Exporter(), OverlapMode_));} // Export computed Y values if needed
00566   }
00567   else {
00568 
00569     EPETRA_CHK_ERR(L_->Multiply(Trans, *X1, *Y1));
00570     EPETRA_CHK_ERR(Y1->Update(1.0, *X1, 1.0)); // Y1 = Y1 + X1 (account for implicit unit diagonal)
00571     EPETRA_CHK_ERR(Y1->ReciprocalMultiply(1.0, *D_, *Y1, 0.0)); // y = D*y (D_ has inverse of diagonal)
00572     Epetra_MultiVector Y1temp(*Y1); // Need a temp copy of Y1
00573     EPETRA_CHK_ERR(U_->Multiply(Trans, Y1temp, *Y1));
00574     EPETRA_CHK_ERR(Y1->Update(1.0, Y1temp, 1.0)); // (account for implicit unit diagonal)
00575     if (IsOverlapped_) {EPETRA_CHK_ERR(Y.Export(*Y1,*L_->Exporter(), OverlapMode_));}
00576   } 
00577   return(0);
00578 }
00579 //=============================================================================
00580 int Ifpack_CrsRiluk::Condest(bool Trans, double & ConditionNumberEstimate) const {
00581 
00582   if (Condest_>=0.0) {
00583     ConditionNumberEstimate = Condest_;
00584     return(0);
00585   }
00586   // Create a vector with all values equal to one
00587   Epetra_Vector Ones(U_->DomainMap());
00588   Epetra_Vector OnesResult(L_->RangeMap());
00589   Ones.PutScalar(1.0);
00590 
00591   EPETRA_CHK_ERR(Solve(Trans, Ones, OnesResult)); // Compute the effect of the solve on the vector of ones
00592   EPETRA_CHK_ERR(OnesResult.Abs(OnesResult)); // Make all values non-negative
00593   EPETRA_CHK_ERR(OnesResult.MaxValue(&ConditionNumberEstimate)); // Get the maximum value across all processors
00594   Condest_ = ConditionNumberEstimate; // Save value for possible later calls
00595   return(0);
00596 }
00597 //==============================================================================
00598 int Ifpack_CrsRiluk::BlockGraph2PointGraph(const Epetra_CrsGraph & BG, Epetra_CrsGraph & PG, bool Upper) {
00599 
00600   if (!BG.IndicesAreLocal()) {EPETRA_CHK_ERR(-1);} // Must have done FillComplete on BG
00601 
00602   int * ColFirstPointInElementList = BG.RowMap().FirstPointInElementList();
00603   int * ColElementSizeList = BG.RowMap().ElementSizeList();
00604   if (BG.Importer()!=0) {
00605     ColFirstPointInElementList = BG.ImportMap().FirstPointInElementList();
00606     ColElementSizeList = BG.ImportMap().ElementSizeList();
00607   }
00608 
00609   int Length = (BG.MaxNumIndices()+1) * BG.ImportMap().MaxMyElementSize();
00610   vector<int> tmpIndices(Length);
00611 
00612   int BlockRow, BlockOffset, NumEntries;
00613   int NumBlockEntries;
00614   int * BlockIndices;
00615 
00616   int NumMyRows_tmp = PG.NumMyRows();
00617 
00618   for (int i=0; i<NumMyRows_tmp; i++) {
00619     EPETRA_CHK_ERR(BG.RowMap().FindLocalElementID(i, BlockRow, BlockOffset));
00620     EPETRA_CHK_ERR(BG.ExtractMyRowView(BlockRow, NumBlockEntries, BlockIndices));
00621 
00622     int * ptr = &tmpIndices[0]; // Set pointer to beginning of buffer
00623 
00624     int RowDim = BG.RowMap().ElementSize(BlockRow);
00625     NumEntries = 0;
00626 
00627     // This next line make sure that the off-diagonal entries in the block diagonal of the 
00628     // original block entry matrix are included in the nonzero pattern of the point graph
00629     if (Upper) {
00630       int jstart = i+1;
00631       int jstop = EPETRA_MIN(NumMyRows_tmp,i+RowDim-BlockOffset);
00632       for (int j= jstart; j< jstop; j++) {*ptr++ = j; NumEntries++;}
00633     }
00634 
00635     for (int j=0; j<NumBlockEntries; j++) {
00636       int ColDim = ColElementSizeList[BlockIndices[j]];
00637       NumEntries += ColDim;
00638       assert(NumEntries<=Length); // Sanity test
00639       int Index = ColFirstPointInElementList[BlockIndices[j]];
00640       for (int k=0; k < ColDim; k++) *ptr++ = Index++;
00641     }
00642 
00643     // This next line make sure that the off-diagonal entries in the block diagonal of the 
00644     // original block entry matrix are included in the nonzero pattern of the point graph
00645     if (!Upper) {
00646       int jstart = EPETRA_MAX(0,i-RowDim+1);
00647       int jstop = i;
00648       for (int j = jstart; j < jstop; j++) {*ptr++ = j; NumEntries++;}
00649     }
00650 
00651     EPETRA_CHK_ERR(PG.InsertMyIndices(i, NumEntries, &tmpIndices[0]));
00652   }
00653 
00654   SetAllocated(true);
00655 
00656   return(0);
00657 }
00658 //=========================================================================
00659 int Ifpack_CrsRiluk::BlockMap2PointMap(const Epetra_BlockMap & BlockMap, Teuchos::RefCountPtr<Epetra_Map>* PointMap) {
00660   // Generate an Epetra_Map that has the same number and distribution of points
00661   // as the input Epetra_BlockMap object.  The global IDs for the output PointMap
00662   // are computed by using the MaxElementSize of the BlockMap.  For variable block
00663   // sizes this will create gaps in the GID space, but that is OK for Epetra_Maps.
00664 
00665   int MaxElementSize = BlockMap.MaxElementSize();
00666   int PtNumMyElements = BlockMap.NumMyPoints();
00667   vector<int> PtMyGlobalElements;
00668   if (PtNumMyElements>0) PtMyGlobalElements.resize(PtNumMyElements);
00669 
00670   int NumMyElements = BlockMap.NumMyElements();
00671 
00672   int curID = 0;
00673   for (int i=0; i<NumMyElements; i++) {
00674     int StartID = BlockMap.GID(i)*MaxElementSize;
00675     int ElementSize = BlockMap.ElementSize(i);
00676     for (int j=0; j<ElementSize; j++) PtMyGlobalElements[curID++] = StartID+j;
00677   }
00678   assert(curID==PtNumMyElements); // Sanity test
00679 
00680   (*PointMap) = Teuchos::rcp( new Epetra_Map(-1, PtNumMyElements, &PtMyGlobalElements[0], BlockMap.IndexBase(), BlockMap.Comm()) );
00681 
00682   if (!BlockMap.PointSameAs(*(*PointMap))) {EPETRA_CHK_ERR(-1);} // Maps not compatible
00683   return(0);
00684 }
00685 //=========================================================================
00686 int Ifpack_CrsRiluk::GenerateXY(bool Trans, 
00687         const Epetra_MultiVector& Xin, const Epetra_MultiVector& Yin,
00688         Teuchos::RefCountPtr<Epetra_MultiVector>* Xout, 
00689         Teuchos::RefCountPtr<Epetra_MultiVector>* Yout) const {
00690 
00691   // Generate an X and Y suitable for performing Solve() and Multiply() methods
00692 
00693   if (Xin.NumVectors()!=Yin.NumVectors()) EPETRA_CHK_ERR(-1); // Return error: X and Y not the same size
00694 
00695   //cout << "Xin = " << Xin << endl;
00696   (*Xout) = Teuchos::rcp( (Epetra_MultiVector *) &Xin, false );
00697   (*Yout) = Teuchos::rcp( (Epetra_MultiVector *) &Yin, false );
00698   if (!IsOverlapped_ && UserMatrixIsCrs_) return(0); // Nothing more to do
00699 
00700   if (UserMatrixIsVbr_) {
00701     if (VbrX_!=Teuchos::null) {
00702       if (VbrX_->NumVectors()!=Xin.NumVectors()) {
00703   VbrX_ = Teuchos::null;
00704   VbrY_ = Teuchos::null;
00705       }
00706     }
00707     if (VbrX_==Teuchos::null) { // Need to allocate space for overlap X and Y
00708       VbrX_ = Teuchos::rcp( new Epetra_MultiVector(View, *U_DomainMap_, (*Xout)->Pointers(), (*Xout)->NumVectors()) );
00709       VbrY_ = Teuchos::rcp( new Epetra_MultiVector(View, *L_RangeMap_, (*Yout)->Pointers(), (*Yout)->NumVectors()) );
00710     }
00711     else {
00712       EPETRA_CHK_ERR(VbrX_->ResetView((*Xout)->Pointers()));
00713       EPETRA_CHK_ERR(VbrY_->ResetView((*Yout)->Pointers()));
00714     }
00715     (*Xout) = VbrX_;
00716     (*Yout) = VbrY_;
00717   }
00718     
00719   if (IsOverlapped_) {
00720     // Make sure the number of vectors in the multivector is the same as before.
00721     if (OverlapX_!=Teuchos::null) {
00722       if (OverlapX_->NumVectors()!=Xin.NumVectors()) {
00723   OverlapX_ = Teuchos::null;
00724   OverlapY_ = Teuchos::null;
00725       }
00726     }
00727     if (OverlapX_==Teuchos::null) { // Need to allocate space for overlap X and Y
00728       OverlapX_ = Teuchos::rcp( new Epetra_MultiVector(U_->RowMatrixColMap(), (*Xout)->NumVectors()) );
00729       OverlapY_ = Teuchos::rcp( new Epetra_MultiVector(L_->RowMatrixRowMap(), (*Yout)->NumVectors()) );
00730     }
00731     if (!Trans) {
00732       EPETRA_CHK_ERR(OverlapX_->Import(*(*Xout),*U_->Importer(), Insert)); // Import X values for solve
00733     }
00734     else {
00735       EPETRA_CHK_ERR(OverlapX_->Import(*(*Xout),*L_->Exporter(), Insert)); // Import X values for solve
00736     }
00737     (*Xout) = OverlapX_;
00738     (*Yout) = OverlapY_; // Set pointers for Xout and Yout to point to overlap space
00739     //cout << "OverlapX_ = " << *OverlapX_ << endl;
00740   }
00741   
00742   return(0);
00743 }
00744 //=============================================================================
00745 // Non-member functions
00746 
00747 ostream& operator << (ostream& os, const Ifpack_CrsRiluk& A)
00748 {
00749 /*  Epetra_fmtflags olda = os.setf(ios::right,ios::adjustfield);
00750   Epetra_fmtflags oldf = os.setf(ios::scientific,ios::floatfield);
00751   int oldp = os.precision(12); */
00752   int LevelFill = A.Graph().LevelFill();
00753   int LevelOverlap = A.Graph().LevelOverlap();
00754   Epetra_CrsMatrix & L = (Epetra_CrsMatrix &) A.L();
00755   Epetra_CrsMatrix & U = (Epetra_CrsMatrix &) A.U();
00756   Epetra_Vector & D = (Epetra_Vector &) A.D();
00757 
00758   os.width(14);
00759   os << endl;
00760   os <<  "     Level of Fill = "; os << LevelFill;
00761   os << endl;
00762   os.width(14);
00763   os <<  "     Level of Overlap = "; os << LevelOverlap;
00764   os << endl;
00765 
00766   os.width(14);
00767   os <<  "     Lower Triangle = ";
00768   os << endl;
00769   os << L; // Let Epetra_CrsMatrix handle the rest.
00770   os << endl;
00771 
00772   os.width(14);
00773   os <<  "     Inverse of Diagonal = ";
00774   os << endl;
00775   os << D; // Let Epetra_Vector handle the rest.
00776   os << endl;
00777 
00778   os.width(14);
00779   os <<  "     Upper Triangle = ";
00780   os << endl;
00781   os << U; // Let Epetra_CrsMatrix handle the rest.
00782   os << endl;
00783  
00784   // Reset os flags
00785 
00786 /*  os.setf(olda,ios::adjustfield);
00787   os.setf(oldf,ios::floatfield);
00788   os.precision(oldp); */
00789 
00790   return os;
00791 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines