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

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