Ifpack_CrsRiluk.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 #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) 
00042   : UserMatrixIsVbr_(false),
00043     UserMatrixIsCrs_(false),
00044     Graph_(Graph),
00045     Comm_(Graph.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.LevelOverlap()>0 && Graph.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   int current_madds = 0; // We will count multiply-add as they happen
00376 
00377   // Now start the factorization.
00378 
00379   // Need some integer workspace and pointers
00380   int NumUU; 
00381   int * UUI;
00382   double * UUV;
00383   for (j=0; j<NumMyCols(); j++) colflag[j] = - 1;
00384 
00385   for(i=0; i<NumMyRows(); i++) {
00386 
00387  // Fill InV, InI with current row of L, D and U combined
00388 
00389     NumIn = MaxNumEntries;
00390     EPETRA_CHK_ERR(L_->ExtractMyRowCopy(i, NumIn, NumL, &InV[0], &InI[0]));
00391     LV = &InV[0];
00392     LI = &InI[0];
00393 
00394     InV[NumL] = DV[i]; // Put in diagonal
00395     InI[NumL] = i;
00396     
00397     EPETRA_CHK_ERR(U_->ExtractMyRowCopy(i, NumIn-NumL-1, NumU, &InV[NumL+1], &InI[NumL+1]));
00398     NumIn = NumL+NumU+1;
00399     UV = &InV[NumL+1];
00400     UI = &InI[NumL+1];
00401 
00402     // Set column flags
00403     for (j=0; j<NumIn; j++) colflag[InI[j]] = j;
00404 
00405     double diagmod = 0.0; // Off-diagonal accumulator
00406 
00407     for (int jj=0; jj<NumL; jj++) {
00408       j = InI[jj];
00409       double multiplier = InV[jj]; // current_mults++;
00410 
00411       InV[jj] *= DV[j];
00412       
00413       EPETRA_CHK_ERR(U_->ExtractMyRowView(j, NumUU, UUV, UUI)); // View of row above
00414 
00415       if (RelaxValue_==0.0) {
00416     for (k=0; k<NumUU; k++) {
00417       int kk = colflag[UUI[k]];
00418       if (kk>-1) {
00419         InV[kk] -= multiplier*UUV[k];
00420         current_madds++;
00421       }
00422     }
00423       }
00424       else {
00425     for (k=0; k<NumUU; k++) {
00426       int kk = colflag[UUI[k]];
00427       if (kk>-1) InV[kk] -= multiplier*UUV[k];
00428       else diagmod -= multiplier*UUV[k];
00429       current_madds++;
00430     }
00431       }
00432      }
00433     if (NumL) {
00434       EPETRA_CHK_ERR(L_->ReplaceMyValues(i, NumL, LV, LI));  // Replace current row of L
00435     }
00436 
00437     DV[i] = InV[NumL]; // Extract Diagonal value
00438 
00439     if (RelaxValue_!=0.0) {
00440       DV[i] += RelaxValue_*diagmod; // Add off diagonal modifications
00441       // current_madds++;
00442     }
00443 
00444     if (fabs(DV[i]) > MaxDiagonalValue) {
00445       if (DV[i] < 0) DV[i] = - MinDiagonalValue;
00446       else DV[i] = MinDiagonalValue;
00447     }
00448     else
00449       DV[i] = 1.0/DV[i]; // Invert diagonal value
00450 
00451     for (j=0; j<NumU; j++) UV[j] *= DV[i]; // Scale U by inverse of diagonal
00452 
00453     if (NumU) {
00454       EPETRA_CHK_ERR(U_->ReplaceMyValues(i, NumU, UV, UI));  // Replace current row of L and U
00455     }
00456 
00457     // Reset column flags
00458     for (j=0; j<NumIn; j++) colflag[InI[j]] = -1;
00459   }
00460 
00461   // Validate that the L and U factors are actually lower and upper triangular
00462 
00463   if( !L_->LowerTriangular() ) 
00464     EPETRA_CHK_ERR(-2);
00465   if( !U_->UpperTriangular() ) 
00466     EPETRA_CHK_ERR(-3);
00467   
00468   // Add up flops
00469  
00470   double current_flops = 2 * current_madds;
00471   double total_flops = 0;
00472     
00473   EPETRA_CHK_ERR(Graph_.L_Graph().RowMap().Comm().SumAll(&current_flops, &total_flops, 1)); // Get total madds across all PEs
00474 
00475   // Now count the rest
00476   total_flops += (double) L_->NumGlobalNonzeros(); // Accounts for multiplier above
00477   total_flops += (double) D_->GlobalLength(); // Accounts for reciprocal of diagonal
00478   if (RelaxValue_!=0.0) total_flops += 2 * (double)D_->GlobalLength(); // Accounts for relax update of diag
00479 
00480   UpdateFlops(total_flops); // Update flop count
00481 
00482   SetFactored(true);
00483 
00484   return(ierr);
00485 
00486 }
00487 
00488 //=============================================================================
00489 int Ifpack_CrsRiluk::Solve(bool Trans, const Epetra_MultiVector& X, 
00490                 Epetra_MultiVector& Y) const {
00491 //
00492 // This function finds Y such that LDU Y = X or U(trans) D L(trans) Y = X for multiple RHS
00493 //
00494 
00495   // First generate X and Y as needed for this function
00496   Teuchos::RefCountPtr<Epetra_MultiVector> X1;
00497   Teuchos::RefCountPtr<Epetra_MultiVector> Y1;
00498   EPETRA_CHK_ERR(GenerateXY(Trans, X, Y, &X1, &Y1));
00499 
00500   bool Upper = true;
00501   bool Lower = false;
00502   bool UnitDiagonal = true;
00503 
00504   Epetra_Flops * counter = this->GetFlopCounter();
00505   if (counter!=0) {
00506     L_->SetFlopCounter(*counter);
00507     Y1->SetFlopCounter(*counter);
00508     U_->SetFlopCounter(*counter);
00509   }
00510 
00511   if (!Trans) {
00512 
00513     EPETRA_CHK_ERR(L_->Solve(Lower, Trans, UnitDiagonal, *X1, *Y1));
00514     EPETRA_CHK_ERR(Y1->Multiply(1.0, *D_, *Y1, 0.0)); // y = D*y (D_ has inverse of diagonal)
00515     EPETRA_CHK_ERR(U_->Solve(Upper, Trans, UnitDiagonal, *Y1, *Y1)); // Solve Uy = y
00516     if (IsOverlapped_) {EPETRA_CHK_ERR(Y.Export(*Y1,*L_->Exporter(), OverlapMode_));} // Export computed Y values if needed
00517   }
00518   else {
00519     EPETRA_CHK_ERR(U_->Solve(Upper, Trans, UnitDiagonal, *X1, *Y1)); // Solve Uy = y
00520     EPETRA_CHK_ERR(Y1->Multiply(1.0, *D_, *Y1, 0.0)); // y = D*y (D_ has inverse of diagonal)
00521     EPETRA_CHK_ERR(L_->Solve(Lower, Trans, UnitDiagonal, *Y1, *Y1));
00522     if (IsOverlapped_) {EPETRA_CHK_ERR(Y.Export(*Y1,*U_->Importer(), OverlapMode_));} // Export computed Y values if needed
00523   } 
00524 
00525   return(0);
00526 }
00527 //=============================================================================
00528 int Ifpack_CrsRiluk::Multiply(bool Trans, const Epetra_MultiVector& X, 
00529                   Epetra_MultiVector& Y) const {
00530 //
00531 // This function finds X such that LDU Y = X or U(trans) D L(trans) Y = X for multiple RHS
00532 //
00533     
00534   // First generate X and Y as needed for this function
00535   Teuchos::RefCountPtr<Epetra_MultiVector> X1;
00536   Teuchos::RefCountPtr<Epetra_MultiVector> Y1;
00537   EPETRA_CHK_ERR(GenerateXY(Trans, X, Y, &X1, &Y1));
00538 
00539   Epetra_Flops * counter = this->GetFlopCounter();
00540   if (counter!=0) {
00541     L_->SetFlopCounter(*counter);
00542     Y1->SetFlopCounter(*counter);
00543     U_->SetFlopCounter(*counter);
00544   }
00545 
00546   if (!Trans) {
00547     EPETRA_CHK_ERR(U_->Multiply(Trans, *X1, *Y1)); // 
00548     EPETRA_CHK_ERR(Y1->Update(1.0, *X1, 1.0)); // Y1 = Y1 + X1 (account for implicit unit diagonal)
00549     EPETRA_CHK_ERR(Y1->ReciprocalMultiply(1.0, *D_, *Y1, 0.0)); // y = D*y (D_ has inverse of diagonal)
00550     Epetra_MultiVector Y1temp(*Y1); // Need a temp copy of Y1
00551     EPETRA_CHK_ERR(L_->Multiply(Trans, Y1temp, *Y1));
00552     EPETRA_CHK_ERR(Y1->Update(1.0, Y1temp, 1.0)); // (account for implicit unit diagonal)
00553     if (IsOverlapped_) {EPETRA_CHK_ERR(Y.Export(*Y1,*L_->Exporter(), OverlapMode_));} // Export computed Y values if needed
00554   }
00555   else {
00556 
00557     EPETRA_CHK_ERR(L_->Multiply(Trans, *X1, *Y1));
00558     EPETRA_CHK_ERR(Y1->Update(1.0, *X1, 1.0)); // Y1 = Y1 + X1 (account for implicit unit diagonal)
00559     EPETRA_CHK_ERR(Y1->ReciprocalMultiply(1.0, *D_, *Y1, 0.0)); // y = D*y (D_ has inverse of diagonal)
00560     Epetra_MultiVector Y1temp(*Y1); // Need a temp copy of Y1
00561     EPETRA_CHK_ERR(U_->Multiply(Trans, Y1temp, *Y1));
00562     EPETRA_CHK_ERR(Y1->Update(1.0, Y1temp, 1.0)); // (account for implicit unit diagonal)
00563     if (IsOverlapped_) {EPETRA_CHK_ERR(Y.Export(*Y1,*L_->Exporter(), OverlapMode_));}
00564   } 
00565   return(0);
00566 }
00567 //=============================================================================
00568 int Ifpack_CrsRiluk::Condest(bool Trans, double & ConditionNumberEstimate) const {
00569 
00570   if (Condest_>=0.0) {
00571     ConditionNumberEstimate = Condest_;
00572     return(0);
00573   }
00574   // Create a vector with all values equal to one
00575   Epetra_Vector Ones(U_->DomainMap());
00576   Epetra_Vector OnesResult(L_->RangeMap());
00577   Ones.PutScalar(1.0);
00578 
00579   EPETRA_CHK_ERR(Solve(Trans, Ones, OnesResult)); // Compute the effect of the solve on the vector of ones
00580   EPETRA_CHK_ERR(OnesResult.Abs(OnesResult)); // Make all values non-negative
00581   EPETRA_CHK_ERR(OnesResult.MaxValue(&ConditionNumberEstimate)); // Get the maximum value across all processors
00582   Condest_ = ConditionNumberEstimate; // Save value for possible later calls
00583   return(0);
00584 }
00585 //==============================================================================
00586 int Ifpack_CrsRiluk::BlockGraph2PointGraph(const Epetra_CrsGraph & BG, Epetra_CrsGraph & PG, bool Upper) {
00587 
00588   if (!BG.IndicesAreLocal()) {EPETRA_CHK_ERR(-1);} // Must have done FillComplete on BG
00589 
00590   int * ColFirstPointInElementList = BG.RowMap().FirstPointInElementList();
00591   int * ColElementSizeList = BG.RowMap().ElementSizeList();
00592   if (BG.Importer()!=0) {
00593     ColFirstPointInElementList = BG.ImportMap().FirstPointInElementList();
00594     ColElementSizeList = BG.ImportMap().ElementSizeList();
00595   }
00596 
00597   int Length = (BG.MaxNumIndices()+1) * BG.ImportMap().MaxMyElementSize();
00598   vector<int> tmpIndices(Length);
00599 
00600   int BlockRow, BlockOffset, NumEntries;
00601   int NumBlockEntries;
00602   int * BlockIndices;
00603 
00604   int NumMyRows = PG.NumMyRows();
00605 
00606   for (int i=0; i<NumMyRows; i++) {
00607     EPETRA_CHK_ERR(BG.RowMap().FindLocalElementID(i, BlockRow, BlockOffset));
00608     EPETRA_CHK_ERR(BG.ExtractMyRowView(BlockRow, NumBlockEntries, BlockIndices));
00609 
00610     int * ptr = &tmpIndices[0]; // Set pointer to beginning of buffer
00611 
00612     int RowDim = BG.RowMap().ElementSize(BlockRow);
00613     NumEntries = 0;
00614 
00615     // This next line make sure that the off-diagonal entries in the block diagonal of the 
00616     // original block entry matrix are included in the nonzero pattern of the point graph
00617     if (Upper) {
00618       int jstart = i+1;
00619       int jstop = EPETRA_MIN(NumMyRows,i+RowDim-BlockOffset);
00620       for (int j= jstart; j< jstop; j++) {*ptr++ = j; NumEntries++;}
00621     }
00622 
00623     for (int j=0; j<NumBlockEntries; j++) {
00624       int ColDim = ColElementSizeList[BlockIndices[j]];
00625       NumEntries += ColDim;
00626       assert(NumEntries<=Length); // Sanity test
00627       int Index = ColFirstPointInElementList[BlockIndices[j]];
00628       for (int k=0; k < ColDim; k++) *ptr++ = Index++;
00629     }
00630 
00631     // This next line make sure that the off-diagonal entries in the block diagonal of the 
00632     // original block entry matrix are included in the nonzero pattern of the point graph
00633     if (!Upper) {
00634       int jstart = EPETRA_MAX(0,i-RowDim+1);
00635       int jstop = i;
00636       for (int j = jstart; j < jstop; j++) {*ptr++ = j; NumEntries++;}
00637     }
00638 
00639     EPETRA_CHK_ERR(PG.InsertMyIndices(i, NumEntries, &tmpIndices[0]));
00640   }
00641 
00642   SetAllocated(true);
00643 
00644   return(0);
00645 }
00646 //=========================================================================
00647 int Ifpack_CrsRiluk::BlockMap2PointMap(const Epetra_BlockMap & BlockMap, Teuchos::RefCountPtr<Epetra_Map>* PointMap) {
00648     // Generate an Epetra_Map that has the same number and distribution of points
00649     // as the input Epetra_BlockMap object.  The global IDs for the output PointMap
00650     // are computed by using the MaxElementSize of the BlockMap.  For variable block
00651     // sizes this will create gaps in the GID space, but that is OK for Epetra_Maps.
00652 
00653     int MaxElementSize = BlockMap.MaxElementSize();
00654     int PtNumMyElements = BlockMap.NumMyPoints();
00655     vector<int> PtMyGlobalElements;
00656     if (PtNumMyElements>0) PtMyGlobalElements.resize(PtNumMyElements);
00657 
00658     int NumMyElements = BlockMap.NumMyElements();
00659 
00660     int curID = 0;
00661     for (int i=0; i<NumMyElements; i++) {
00662         int StartID = BlockMap.GID(i)*MaxElementSize;
00663         int ElementSize = BlockMap.ElementSize(i);
00664         for (int j=0; j<ElementSize; j++) PtMyGlobalElements[curID++] = StartID+j;
00665     }
00666     assert(curID==PtNumMyElements); // Sanity test
00667 
00668     (*PointMap) = Teuchos::rcp( new Epetra_Map(-1, PtNumMyElements, &PtMyGlobalElements[0], BlockMap.IndexBase(), BlockMap.Comm()) );
00669 
00670     if (!BlockMap.PointSameAs(*(*PointMap))) {EPETRA_CHK_ERR(-1);} // Maps not compatible
00671   return(0);
00672 }
00673 //=========================================================================
00674 int Ifpack_CrsRiluk::GenerateXY(bool Trans, 
00675                 const Epetra_MultiVector& Xin, const Epetra_MultiVector& Yin,
00676                 Teuchos::RefCountPtr<Epetra_MultiVector>* Xout, 
00677                 Teuchos::RefCountPtr<Epetra_MultiVector>* Yout) const {
00678 
00679   // Generate an X and Y suitable for performing Solve() and Multiply() methods
00680 
00681   if (Xin.NumVectors()!=Yin.NumVectors()) EPETRA_CHK_ERR(-1); // Return error: X and Y not the same size
00682 
00683   //cout << "Xin = " << Xin << endl;
00684   (*Xout) = Teuchos::rcp( (Epetra_MultiVector *) &Xin, false );
00685   (*Yout) = Teuchos::rcp( (Epetra_MultiVector *) &Yin, false );
00686   if (!IsOverlapped_ && UserMatrixIsCrs_) return(0); // Nothing more to do
00687 
00688   if (UserMatrixIsVbr_) {
00689     if (VbrX_!=Teuchos::null) {
00690       if (VbrX_->NumVectors()!=Xin.NumVectors()) {
00691     VbrX_ = Teuchos::null;
00692     VbrY_ = Teuchos::null;
00693       }
00694     }
00695     if (VbrX_==Teuchos::null) { // Need to allocate space for overlap X and Y
00696       VbrX_ = Teuchos::rcp( new Epetra_MultiVector(View, *U_DomainMap_, (*Xout)->Pointers(), (*Xout)->NumVectors()) );
00697       VbrY_ = Teuchos::rcp( new Epetra_MultiVector(View, *L_RangeMap_, (*Yout)->Pointers(), (*Yout)->NumVectors()) );
00698     }
00699     else {
00700       EPETRA_CHK_ERR(VbrX_->ResetView((*Xout)->Pointers()));
00701       EPETRA_CHK_ERR(VbrY_->ResetView((*Yout)->Pointers()));
00702     }
00703     (*Xout) = VbrX_;
00704     (*Yout) = VbrY_;
00705   }
00706     
00707   if (IsOverlapped_) {
00708     // Make sure the number of vectors in the multivector is the same as before.
00709     if (OverlapX_!=Teuchos::null) {
00710       if (OverlapX_->NumVectors()!=Xin.NumVectors()) {
00711     OverlapX_ = Teuchos::null;
00712     OverlapY_ = Teuchos::null;
00713       }
00714     }
00715     if (OverlapX_==Teuchos::null) { // Need to allocate space for overlap X and Y
00716       OverlapX_ = Teuchos::rcp( new Epetra_MultiVector(U_->RowMatrixColMap(), (*Xout)->NumVectors()) );
00717       OverlapY_ = Teuchos::rcp( new Epetra_MultiVector(L_->RowMatrixRowMap(), (*Yout)->NumVectors()) );
00718     }
00719     if (!Trans) {
00720       EPETRA_CHK_ERR(OverlapX_->Import(*(*Xout),*U_->Importer(), Insert)); // Import X values for solve
00721     }
00722     else {
00723       EPETRA_CHK_ERR(OverlapX_->Import(*(*Xout),*L_->Exporter(), Insert)); // Import X values for solve
00724     }
00725     (*Xout) = OverlapX_;
00726     (*Yout) = OverlapY_; // Set pointers for Xout and Yout to point to overlap space
00727     //cout << "OverlapX_ = " << *OverlapX_ << endl;
00728   }
00729   
00730   return(0);
00731 }
00732 //=============================================================================
00733 // Non-member functions
00734 
00735 ostream& operator << (ostream& os, const Ifpack_CrsRiluk& A)
00736 {
00737 /*  Epetra_fmtflags olda = os.setf(ios::right,ios::adjustfield);
00738   Epetra_fmtflags oldf = os.setf(ios::scientific,ios::floatfield);
00739   int oldp = os.precision(12); */
00740   int LevelFill = A.Graph().LevelFill();
00741   int LevelOverlap = A.Graph().LevelOverlap();
00742   Epetra_CrsMatrix & L = (Epetra_CrsMatrix &) A.L();
00743   Epetra_CrsMatrix & U = (Epetra_CrsMatrix &) A.U();
00744   Epetra_Vector & D = (Epetra_Vector &) A.D();
00745 
00746   os.width(14);
00747   os << endl;
00748   os <<  "     Level of Fill = "; os << LevelFill;
00749   os << endl;
00750   os.width(14);
00751   os <<  "     Level of Overlap = "; os << LevelOverlap;
00752   os << endl;
00753 
00754   os.width(14);
00755   os <<  "     Lower Triangle = ";
00756   os << endl;
00757   os << L; // Let Epetra_CrsMatrix handle the rest.
00758   os << endl;
00759 
00760   os.width(14);
00761   os <<  "     Inverse of Diagonal = ";
00762   os << endl;
00763   os << D; // Let Epetra_Vector handle the rest.
00764   os << endl;
00765 
00766   os.width(14);
00767   os <<  "     Upper Triangle = ";
00768   os << endl;
00769   os << U; // Let Epetra_CrsMatrix handle the rest.
00770   os << endl;
00771  
00772   // Reset os flags
00773 
00774 /*  os.setf(olda,ios::adjustfield);
00775   os.setf(oldf,ios::floatfield);
00776   os.precision(oldp); */
00777 
00778   return os;
00779 }

Generated on Tue Oct 20 12:48:53 2009 for IFPACK by doxygen 1.4.7