Ifpack2 Templated Preconditioning Package Version 1.0
Ifpack2_RILUK_def.hpp
00001 //@HEADER
00002 // ***********************************************************************
00003 //
00004 //       Ifpack2: Templated Object-Oriented Algebraic Preconditioner Package
00005 //                 Copyright (2010) 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 #ifndef IFPACK2_CRSRILUK_DEF_HPP
00030 #define IFPACK2_CRSRILUK_DEF_HPP
00031 
00032 namespace Ifpack2 {
00033 
00034 //==============================================================================
00035 template<class MatrixType>
00036 RILUK<MatrixType>::RILUK(const Teuchos::RCP<const MatrixType>& Matrix_in)
00037   : isOverlapped_(false),
00038     Graph_(),
00039     A_(Matrix_in),
00040     UseTranspose_(false),
00041     LevelOfFill_(0),
00042     LevelOfOverlap_(0),
00043     NumMyDiagonals_(0),
00044     isAllocated_(false),
00045     isInitialized_(false),
00046     numInitialize_(0),
00047     numCompute_(0),
00048     numApply_(0),
00049     Factored_(false),
00050     RelaxValue_(0.0),
00051     Athresh_(0.0),
00052     Rthresh_(1.0),
00053     Condest_(-1.0),
00054     OverlapMode_(Tpetra::REPLACE)
00055 {
00056 }
00057 
00058 //==============================================================================
00059 //template<class MatrixType>
00060 //RILUK<MatrixType>::RILUK(const RILUK<MatrixType>& src)
00061 //  : isOverlapped_(src.isOverlapped_),
00062 //    Graph_(src.Graph_),
00063 //    UseTranspose_(src.UseTranspose_),
00064 //    LevelOfFill_(src.LevelOfFill_),
00065 //    LevelOfOverlap_(src.LevelOfOverlap_),
00066 //    NumMyDiagonals_(src.NumMyDiagonals_),
00067 //    isAllocated_(src.isAllocated_),
00068 //    isInitialized_(src.isInitialized_),
00069 //    Factored_(src.Factored_),
00070 //    RelaxValue_(src.RelaxValue_),
00071 //    Athresh_(src.Athresh_),
00072 //    Rthresh_(src.Rthresh_),
00073 //    Condest_(src.Condest_),
00074 //    OverlapMode_(src.OverlapMode_)
00075 //{
00076 //  L_ = Teuchos::rcp( new MatrixType(src.getL()) );
00077 //  U_ = Teuchos::rcp( new MatrixType(src.getU()) );
00078 //  D_ = Teuchos::rcp( new Tpetra::Vector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>(src.getD()) );
00079 //}
00080 
00081 //==============================================================================
00082 template<class MatrixType>
00083 RILUK<MatrixType>::~RILUK() {
00084 }
00085 
00086 //==============================================================================
00087 template<class MatrixType>
00088 void RILUK<MatrixType>::allocate_L_and_U() {
00089 
00090   // Allocate Matrix using ILUK graphs
00091   L_ = Teuchos::rcp( new MatrixType(Graph_->getL_Graph()) );
00092   U_ = Teuchos::rcp( new MatrixType(Graph_->getU_Graph()) );
00093   L_->setAllToScalar(0.0); // Zero out L and U matrices
00094   U_->setAllToScalar(0.0);
00095   L_->fillComplete();
00096   U_->fillComplete();
00097   bool isLower = L_->isLowerTriangular();
00098   bool isUpper = U_->isUpperTriangular();
00099   if (!isLower || !isUpper) {
00100     std::cout << "error in triangular detection" << std::endl;
00101   }
00102 
00103   D_ = Teuchos::rcp( new Tpetra::Vector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>(Graph_->getL_Graph()->getRowMap()) );
00104   setAllocated(true);
00105 }
00106 
00107 //==========================================================================
00108 template<class MatrixType>
00109 void
00110 RILUK<MatrixType>::
00111 setParameters (const Teuchos::ParameterList& params)
00112 {
00113   using Teuchos::as;
00114   using Teuchos::Exceptions::InvalidParameterName;
00115   using Teuchos::Exceptions::InvalidParameterType;
00116   typedef Teuchos::ScalarTraits<magnitude_type> STM;
00117 
00118   // Default values of the various parameters.
00119   int fillLevel = 0;
00120   int overlapLevel = 0;
00121   magnitude_type absThresh = STM::zero ();
00122   magnitude_type relThresh = STM::one ();
00123   magnitude_type relaxValue = STM::zero ();
00124 
00125   //
00126   // "fact: iluk level-of-fill" parsing is more complicated, because
00127   // we want to allow as many types as make sense.  int is the native
00128   // type, but we also want to accept magnitude_type (for
00129   // compatibility with ILUT) and double (for backwards compatibilty
00130   // with ILUT).
00131   //
00132 
00133   bool gotFillLevel = false;
00134   try {
00135     fillLevel = params.get<int> ("fact: iluk level-of-fill");
00136     gotFillLevel = true;
00137   }
00138   catch (InvalidParameterType&) {
00139     // Throwing again in the catch block would just unwind the stack.
00140     // Instead, we do nothing here, and check the Boolean outside to
00141     // see if we got the value.
00142   }
00143   catch (InvalidParameterName&) {
00144     gotFillLevel = true; // Accept the default value.
00145   }
00146 
00147   if (! gotFillLevel) {
00148     try {
00149       // Try magnitude_type, for compatibility with ILUT.
00150       // The cast from magnitude_type to int must succeed.
00151       fillLevel = as<int> (params.get<magnitude_type> ("fact: iluk level-of-fill"));
00152       gotFillLevel = true;
00153     }
00154     catch (InvalidParameterType&) {
00155       // Try double next.
00156     }
00157     // Don't catch InvalidParameterName here; we've already done that above.
00158   }
00159 
00160   if (! gotFillLevel) {
00161     try {
00162       // Try double, for compatibility with ILUT.
00163       // The cast from double to int must succeed.
00164       fillLevel = as<int> (params.get<double> ("fact: iluk level-of-fill"));
00165       gotFillLevel = true;
00166     }
00167     catch (InvalidParameterType& e) {
00168       // We're out of options.  The user gave us the parameter, but it
00169       // doesn't have the right type.  The best thing for us to do in
00170       // that case is to throw, telling the user to use the right
00171       // type.
00172       throw e;
00173     }
00174     // Don't catch InvalidParameterName here; we've already done that above.
00175   }
00176 
00177   TEUCHOS_TEST_FOR_EXCEPTION(
00178     ! gotFillLevel,
00179     std::logic_error,
00180     "Ifpack2::RILUK::setParameters: We should never get here!  "
00181     "The method should either have read the \"fact: iluk level-of-fill\"  "
00182     "parameter by this point, or have thrown an exception.  "
00183     "Please let the Ifpack2 developers know about this bug.");
00184 
00185   //
00186   // overlapLevel was always int.  ILUT doesn't have this parameter.
00187   // However, some tests (as of 28 Nov 2012:
00188   // Ifpack2_RILUK_small_belos_MPI_1) depend on being able to set this
00189   // as a double instead of an int.  Thus, we go through the same
00190   // procedure as above with fill level.
00191   //
00192 
00193   bool gotOverlapLevel = false;
00194   try {
00195     overlapLevel = params.get<int> ("fact: iluk level-of-overlap");
00196     gotOverlapLevel = true;
00197   }
00198   catch (InvalidParameterType&) {
00199     // Throwing again in the catch block would just unwind the stack.
00200     // Instead, we do nothing here, and check the Boolean outside to
00201     // see if we got the value.
00202   }
00203   catch (InvalidParameterName&) {
00204     gotOverlapLevel = true; // Accept the default value.
00205   }
00206 
00207   if (! gotOverlapLevel) {
00208     try {
00209       // Try magnitude_type, for compatibility with ILUT.
00210       // The cast from magnitude_type to int must succeed.
00211       overlapLevel = as<int> (params.get<magnitude_type> ("fact: iluk level-of-overlap"));
00212       gotOverlapLevel = true;
00213     }
00214     catch (InvalidParameterType&) {
00215       // Try double next.
00216     }
00217     // Don't catch InvalidParameterName here; we've already done that above.
00218   }
00219 
00220   if (! gotOverlapLevel) {
00221     try {
00222       // Try double, for compatibility with ILUT.
00223       // The cast from double to int must succeed.
00224       overlapLevel = as<int> (params.get<double> ("fact: iluk level-of-overlap"));
00225       gotOverlapLevel = true;
00226     }
00227     catch (InvalidParameterType& e) {
00228       // We're out of options.  The user gave us the parameter, but it
00229       // doesn't have the right type.  The best thing for us to do in
00230       // that case is to throw, telling the user to use the right
00231       // type.
00232       throw e;
00233     }
00234     // Don't catch InvalidParameterName here; we've already done that above.
00235   }
00236 
00237   TEUCHOS_TEST_FOR_EXCEPTION(
00238     ! gotOverlapLevel,
00239     std::logic_error,
00240     "Ifpack2::RILUK::setParameters: We should never get here!  "
00241     "The method should either have read the \"fact: iluk level-of-overlap\"  "
00242     "parameter by this point, or have thrown an exception.  "
00243     "Please let the Ifpack2 developers know about this bug.");
00244 
00245   //
00246   // For the other parameters, we prefer magnitude_type, but allow
00247   // double for backwards compatibility.
00248   //
00249 
00250   try {
00251     absThresh = params.get<magnitude_type> ("fact: absolute threshold");
00252   }
00253   catch (InvalidParameterType&) {
00254     // Try double, for backwards compatibility.
00255     // The cast from double to magnitude_type must succeed.
00256     absThresh = as<magnitude_type> (params.get<double> ("fact: absolute threshold"));
00257   }
00258   catch (InvalidParameterName&) {
00259     // Accept the default value.
00260   }
00261 
00262   try {
00263     relThresh = params.get<magnitude_type> ("fact: relative threshold");
00264   }
00265   catch (InvalidParameterType&) {
00266     // Try double, for backwards compatibility.
00267     // The cast from double to magnitude_type must succeed.
00268     relThresh = as<magnitude_type> (params.get<double> ("fact: relative threshold"));
00269   }
00270   catch (InvalidParameterName&) {
00271     // Accept the default value.
00272   }
00273 
00274   try {
00275     relaxValue = params.get<magnitude_type> ("fact: relax value");
00276   }
00277   catch (InvalidParameterType&) {
00278     // Try double, for backwards compatibility.
00279     // The cast from double to magnitude_type must succeed.
00280     relaxValue = as<magnitude_type> (params.get<double> ("fact: relax value"));
00281   }
00282   catch (InvalidParameterName&) {
00283     // Accept the default value.
00284   }
00285 
00286   // "Commit" the values only after validating all of them.  This
00287   // ensures that there are no side effects if this routine throws an
00288   // exception.
00289 
00290   LevelOfFill_ = fillLevel;
00291   LevelOfOverlap_ = overlapLevel;
00292 
00293   // mfh 28 Nov 2012: The previous code would not assign Athresh_,
00294   // Rthresh_, or RelaxValue_, if the read-in value was -1.  I don't
00295   // know if keeping this behavior is correct, but I'll keep it just
00296   // so as not to change previous behavior.
00297 
00298   if (absThresh != -STM::one ()) {
00299     Athresh_ = absThresh;
00300   }
00301   if (relThresh != -STM::one ()) {
00302     Rthresh_ = relThresh;
00303   }
00304   if (relaxValue != -STM::one ()) {
00305     RelaxValue_ = relaxValue;
00306   }
00307 }
00308 
00309 //==========================================================================
00310 template<class MatrixType>
00311 void RILUK<MatrixType>::initialize() {
00312 
00313   if (Graph_ != Teuchos::null) return;
00314 
00315   Graph_ = Teuchos::rcp(new Ifpack2::IlukGraph<local_ordinal_type,global_ordinal_type,node_type>(A_->getCrsGraph(), LevelOfFill_, LevelOfOverlap_));
00316 
00317   Graph_->constructFilledGraph();
00318 
00319   setInitialized(true);
00320 
00321   if (!isAllocated()) allocate_L_and_U();
00322 
00323   if (isOverlapped_) {
00324     throw std::runtime_error("Error, Ifpack2::RILUK::initialize: overlapping not yet supported.");
00325 //    OverlapA = Teuchos::rcp( new MatrixType(Graph_->OverlapGraph()) );
00326 //    EPETRA_CHK_ERR(OverlapA->Import(A, *Graph_->OverlapImporter(), Insert));
00327 //    EPETRA_CHK_ERR(OverlapA->FillComplete());
00328   }
00329   else {
00330     initAllValues(*A_);
00331   }
00332 
00333   ++numInitialize_;
00334 }
00335 
00336 //==========================================================================
00337 template<class MatrixType>
00338 void RILUK<MatrixType>::initAllValues(const Tpetra::RowMatrix<scalar_type,local_ordinal_type,global_ordinal_type,node_type> & OverlapA) {
00339 
00340   size_t NumIn = 0, NumL = 0, NumU = 0;
00341   bool DiagFound = false;
00342   size_t NumNonzeroDiags = 0;
00343   size_t MaxNumEntries = OverlapA.getGlobalMaxNumRowEntries();
00344 
00345 
00346   Teuchos::Array<global_ordinal_type> InI(MaxNumEntries); // Allocate temp space
00347   Teuchos::Array<global_ordinal_type> LI(MaxNumEntries);
00348   Teuchos::Array<global_ordinal_type> UI(MaxNumEntries);
00349   Teuchos::Array<scalar_type> InV(MaxNumEntries);
00350   Teuchos::Array<scalar_type> LV(MaxNumEntries);
00351   Teuchos::Array<scalar_type> UV(MaxNumEntries);
00352 
00353   bool ReplaceValues = (L_->isStaticGraph() || L_->isLocallyIndexed()); // Check if values should be inserted or replaced
00354 
00355   L_->resumeFill();
00356   U_->resumeFill();
00357   if (ReplaceValues) {
00358     L_->setAllToScalar(0.0); // Zero out L and U matrices
00359     U_->setAllToScalar(0.0);
00360   }
00361 
00362   D_->putScalar(0.0); // Set diagonal values to zero
00363   Teuchos::ArrayRCP<scalar_type> DV = D_->get1dViewNonConst(); // Get view of diagonal
00364 
00365   const Teuchos::RCP<const Tpetra::Map<local_ordinal_type,global_ordinal_type,node_type> >& rowMap =
00366     L_->getRowMap();
00367 
00368   // First we copy the user's matrix into L and U, regardless of fill level
00369 
00370   for (global_ordinal_type i=rowMap->getMinGlobalIndex(); i<=rowMap->getMaxGlobalIndex(); ++i) {
00371     global_ordinal_type global_row = i;
00372     local_ordinal_type local_row = rowMap->getLocalElement(global_row);
00373 
00374     OverlapA.getGlobalRowCopy(global_row, InI(), InV(), NumIn); // Get Values and Indices
00375 
00376     // Split into L and U (we don't assume that indices are ordered).
00377 
00378     NumL = 0;
00379     NumU = 0;
00380     DiagFound = false;
00381 
00382     for (size_t j=0; j< NumIn; j++) {
00383       global_ordinal_type k = InI[j];
00384 
00385       if (k==i) {
00386         DiagFound = true;
00387         DV[local_row] += Rthresh_ * InV[j] + IFPACK2_SGN(InV[j]) * Athresh_; // Store perturbed diagonal in Tpetra::Vector D_
00388       }
00389 
00390       else if (k < 0) { // Out of range
00391         throw std::runtime_error("out of range (k<0) in Ifpack2::RILUK::initAllValues");
00392       }
00393 
00394       else if (k < i) {
00395         LI[NumL] = k;
00396         LV[NumL] = InV[j];
00397         NumL++;
00398       }
00399       else if (k <= rowMap->getMaxGlobalIndex()) {
00400         UI[NumU] = k;
00401         UV[NumU] = InV[j];
00402         NumU++;
00403       }
00404 //      else {
00405 //        throw std::runtime_error("out of range in Ifpack2::RILUK::initAllValues");
00406 //      }
00407     }
00408 
00409     // Check in things for this row of L and U
00410 
00411     if (DiagFound) ++NumNonzeroDiags;
00412     else DV[local_row] = Athresh_;
00413 
00414     if (NumL) {
00415       if (ReplaceValues) {
00416         L_->replaceGlobalValues(global_row, LI(0, NumL), LV(0,NumL));
00417       }
00418       else {
00419         L_->insertGlobalValues(global_row, LI(0,NumL), LV(0,NumL));
00420       }
00421     }
00422 
00423     if (NumU) {
00424       if (ReplaceValues) {
00425         U_->replaceGlobalValues(global_row, UI(0,NumU), UV(0,NumU));
00426       }
00427       else {
00428         U_->insertGlobalValues(global_row, UI(0,NumU), UV(0,NumU));
00429       }
00430     }
00431 
00432   }
00433 
00434   // The domain of L and the range of U are exactly their own row maps (there is no communication).
00435   // The domain of U and the range of L must be the same as those of the original matrix,
00436   // However if the original matrix is a VbrMatrix, these two latter maps are translation from
00437   // a block map to a point map.
00438   L_->fillComplete(L_->getColMap(), A_->getRangeMap());
00439   U_->fillComplete(A_->getDomainMap(), U_->getRowMap());
00440 
00441   // At this point L and U have the values of A in the structure of L and U, and diagonal vector D
00442 
00443   setInitialized(true);
00444   setFactored(false);
00445 
00446   size_t TotalNonzeroDiags = 0;
00447   Teuchos::reduceAll(*L_->getRowMap()->getComm(),Teuchos::REDUCE_SUM,
00448                      1,&NumNonzeroDiags,&TotalNonzeroDiags);
00449   NumMyDiagonals_ = NumNonzeroDiags;
00450   if (NumNonzeroDiags != U_->getNodeNumRows()) {
00451     throw std::runtime_error("Error in Ifpack2::RILUK::initAllValues, wrong number of diagonals.");
00452   }
00453 }
00454 
00455 //==========================================================================
00456 template<class MatrixType>
00457 void RILUK<MatrixType>::compute() {
00458 
00459   L_->resumeFill();
00460   U_->resumeFill();
00461 
00462   TEUCHOS_TEST_FOR_EXCEPTION(!isInitialized(), std::runtime_error,
00463       "Ifpack2::RILUK::compute() ERROR: isInitialized() must be true.");
00464   TEUCHOS_TEST_FOR_EXCEPTION(isComputed() == true, std::runtime_error,
00465       "Ifpack2::RILUK::compute() ERROR: Can't have already computed factors.");
00466 
00467   // MinMachNum should be officially defined, for now pick something a little
00468   // bigger than IEEE underflow value
00469 
00470   scalar_type MinDiagonalValue = Teuchos::ScalarTraits<scalar_type>::rmin();
00471   scalar_type MaxDiagonalValue = Teuchos::ScalarTraits<scalar_type>::one()/MinDiagonalValue;
00472 
00473   size_t NumIn, NumL, NumU;
00474 
00475   // Get Maximun Row length
00476   size_t MaxNumEntries = L_->getNodeMaxNumRowEntries() + U_->getNodeMaxNumRowEntries() + 1;
00477 
00478   Teuchos::Array<local_ordinal_type> InI(MaxNumEntries); // Allocate temp space
00479   Teuchos::Array<scalar_type> InV(MaxNumEntries);
00480   size_t num_cols = U_->getColMap()->getNodeNumElements();
00481   Teuchos::Array<int> colflag(num_cols);
00482 
00483   Teuchos::ArrayRCP<scalar_type> DV = D_->get1dViewNonConst(); // Get view of diagonal
00484 
00485   size_t current_madds = 0; // We will count multiply-add as they happen
00486 
00487   // Now start the factorization.
00488 
00489   // Need some integer workspace and pointers
00490   size_t NumUU;
00491   Teuchos::ArrayView<const local_ordinal_type> UUI;
00492   Teuchos::ArrayView<const scalar_type> UUV;
00493   for (size_t j=0; j<num_cols; j++) colflag[j] = - 1;
00494 
00495   for(size_t i=0; i<L_->getNodeNumRows(); i++) {
00496     local_ordinal_type local_row = i;
00497 
00498  // Fill InV, InI with current row of L, D and U combined
00499 
00500     NumIn = MaxNumEntries;
00501     L_->getLocalRowCopy(local_row, InI(), InV(), NumL);
00502 
00503     InV[NumL] = DV[i]; // Put in diagonal
00504     InI[NumL] = local_row;
00505 
00506     U_->getLocalRowCopy(local_row, InI(NumL+1,MaxNumEntries-NumL-1), InV(NumL+1,MaxNumEntries-NumL-1), NumU);
00507     NumIn = NumL+NumU+1;
00508 
00509     // Set column flags
00510     for (size_t j=0; j<NumIn; j++) colflag[InI[j]] = j;
00511 
00512     scalar_type diagmod = 0.0; // Off-diagonal accumulator
00513 
00514     for (size_t jj=0; jj<NumL; jj++) {
00515       local_ordinal_type j = InI[jj];
00516       scalar_type multiplier = InV[jj]; // current_mults++;
00517 
00518       InV[jj] *= DV[j];
00519 
00520       U_->getLocalRowView(j, UUI, UUV); // View of row above
00521       NumUU = UUI.size();
00522 
00523       if (RelaxValue_==0.0) {
00524         for (size_t k=0; k<NumUU; k++) {
00525           int kk = colflag[UUI[k]];
00526           if (kk>-1) {
00527             InV[kk] -= multiplier*UUV[k];
00528             current_madds++;
00529           }
00530         }
00531       }
00532       else {
00533         for (size_t k=0; k<NumUU; k++) {
00534           int kk = colflag[UUI[k]];
00535           if (kk>-1) InV[kk] -= multiplier*UUV[k];
00536           else diagmod -= multiplier*UUV[k];
00537           current_madds++;
00538         }
00539       }
00540     }
00541     if (NumL) {
00542       L_->replaceLocalValues(local_row, InI(0,NumL), InV(0,NumL));  // Replace current row of L
00543     }
00544 
00545     DV[i] = InV[NumL]; // Extract Diagonal value
00546 
00547     if (RelaxValue_!=0.0) {
00548       DV[i] += RelaxValue_*diagmod; // Add off diagonal modifications
00549       // current_madds++;
00550     }
00551 
00552     if (Teuchos::ScalarTraits<scalar_type>::magnitude(DV[i]) > Teuchos::ScalarTraits<scalar_type>::magnitude(MaxDiagonalValue)) {
00553       if (Teuchos::ScalarTraits<scalar_type>::real(DV[i]) < 0) DV[i] = - MinDiagonalValue;
00554       else DV[i] = MinDiagonalValue;
00555     }
00556     else
00557       DV[i] = Teuchos::ScalarTraits<scalar_type>::one()/DV[i]; // Invert diagonal value
00558 
00559     for (size_t j=0; j<NumU; j++) InV[NumL+1+j] *= DV[i]; // Scale U by inverse of diagonal
00560 
00561     if (NumU) {
00562       U_->replaceLocalValues(local_row, InI(NumL+1,NumU), InV(NumL+1,NumU));  // Replace current row of L and U
00563     }
00564 
00565     // Reset column flags
00566     for (size_t j=0; j<NumIn; j++) colflag[InI[j]] = -1;
00567   }
00568 
00569   L_->fillComplete(L_->getColMap(), A_->getRangeMap());
00570   U_->fillComplete(A_->getDomainMap(), U_->getRowMap());
00571 
00572   // Validate that the L and U factors are actually lower and upper triangular
00573 
00574   if( !L_->isLowerTriangular() )
00575     throw std::runtime_error("Ifpack2::RILUK::compute() ERROR, L isn't lower triangular.");
00576   if( !U_->isUpperTriangular() )
00577     throw std::runtime_error("Ifpack2::RILUK::compute() ERROR, U isn't lower triangular.");
00578 
00579   // Add up flops
00580 
00581   double current_flops = 2 * current_madds;
00582   double total_flops = 0;
00583 
00584   // Get total madds across all PEs
00585   Teuchos::reduceAll(*L_->getRowMap()->getComm(),Teuchos::REDUCE_SUM,
00586                      1,&current_flops,&total_flops);
00587 
00588   // Now count the rest
00589   total_flops += (double) L_->getGlobalNumEntries(); // Accounts for multiplier above
00590   total_flops += (double) D_->getGlobalLength(); // Accounts for reciprocal of diagonal
00591   if (RelaxValue_!=0.0) total_flops += 2 * (double)D_->getGlobalLength(); // Accounts for relax update of diag
00592 
00593   //UpdateFlops(total_flops); // Update flop count
00594 
00595   setFactored(true);
00596   ++numCompute_;
00597 }
00598 
00599 //=============================================================================
00600 template<class MatrixType>
00601 void RILUK<MatrixType>::apply(
00602        const Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& X,
00603              Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& Y,
00604              Teuchos::ETransp mode, scalar_type alpha, scalar_type beta) const
00605 {
00606   typedef Teuchos::ScalarTraits<scalar_type> STS;
00607 
00608   TEUCHOS_TEST_FOR_EXCEPTION(!isComputed(), std::runtime_error,
00609     "Ifpack2::RILUK::apply() ERROR, compute() hasn't been called yet.");
00610 
00611   TEUCHOS_TEST_FOR_EXCEPTION(
00612     alpha != STS::one (), 
00613     std::logic_error,
00614     "Ifpack2::RILUK::apply() does not currently allow alpha != 1.");
00615   TEUCHOS_TEST_FOR_EXCEPTION(
00616     beta != STS::zero (), 
00617     std::logic_error,
00618     "Ifpack2::RILUK::apply() does not currently allow zero != 0.");
00619 
00620 //
00621 // This function finds Y such that
00622 // LDU Y = X, or
00623 // U(trans) D L(trans) Y = X
00624 // for multiple RHS
00625 //
00626 
00627   // First generate X and Y as needed for this function
00628   Teuchos::RCP<const Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type> > X1;
00629   Teuchos::RCP<Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type> > Y1;
00630   generateXY(mode, X, Y, X1, Y1);
00631 
00632   scalar_type one = Teuchos::ScalarTraits<scalar_type>::one();
00633   scalar_type zero = Teuchos::ScalarTraits<scalar_type>::zero();
00634 
00635   if (mode == Teuchos::NO_TRANS) {
00636 
00637     L_->localSolve(*X1, *Y1,mode);
00638     Y1->elementWiseMultiply(one, *D_, *Y1, zero); // y = D*y (D_ has inverse of diagonal)
00639     U_->localSolve(*Y1, *Y1,mode); // Solve Uy = y
00640     if (isOverlapped_) {
00641       // Export computed Y values if needed
00642       Y.doExport(*Y1,*L_->getGraph()->getExporter(), OverlapMode_);
00643     }
00644   }
00645   else {
00646     U_->localSolve(*X1, *Y1,mode); // Solve Uy = y
00647     Y1->elementWiseMultiply(one, *D_, *Y1, zero); // y = D*y (D_ has inverse of diagonal)
00648     L_->localSolve(*Y1, *Y1,mode);
00649     if (isOverlapped_) {Y.doExport(*Y1,*U_->getGraph()->getImporter(), OverlapMode_);} // Export computed Y values if needed
00650   }
00651 
00652   ++numApply_;
00653 }
00654 
00655 //=============================================================================
00656 template<class MatrixType>
00657 int RILUK<MatrixType>::Multiply(const Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& X,
00658                               Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& Y,
00659             Teuchos::ETransp mode) const {
00660 //
00661 // This function finds X such that LDU Y = X or U(trans) D L(trans) Y = X for multiple RHS
00662 //
00663 
00664   // First generate X and Y as needed for this function
00665   Teuchos::RCP<const Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type> > X1;
00666   Teuchos::RCP<Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type> > Y1;
00667   generateXY(mode, X, Y, X1, Y1);
00668 
00669 //  Epetra_Flops * counter = this->GetFlopCounter();
00670 //  if (counter!=0) {
00671 //    L_->SetFlopCounter(*counter);
00672 //    Y1->SetFlopCounter(*counter);
00673 //    U_->SetFlopCounter(*counter);
00674 //  }
00675 
00676   if (!mode == Teuchos::NO_TRANS) {
00677     U_->apply(*X1, *Y1,mode); //
00678     Y1->update(1.0, *X1, 1.0); // Y1 = Y1 + X1 (account for implicit unit diagonal)
00679     Y1->elementWiseMultiply(1.0, *D_, *Y1, 0.0); // y = D*y (D_ has inverse of diagonal)
00680     Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type> Y1temp(*Y1); // Need a temp copy of Y1
00681     L_->apply(Y1temp, *Y1,mode);
00682     Y1->update(1.0, Y1temp, 1.0); // (account for implicit unit diagonal)
00683     if (isOverlapped_) {Y.doExport(*Y1,*L_->getGraph()->getExporter(), OverlapMode_);} // Export computed Y values if needed
00684   }
00685   else {
00686 
00687     L_->apply(*X1, *Y1,mode);
00688     Y1->update(1.0, *X1, 1.0); // Y1 = Y1 + X1 (account for implicit unit diagonal)
00689     Y1->elementWiseMultiply(1, *D_, *Y1, 0); // y = D*y (D_ has inverse of diagonal)
00690     Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type> Y1temp(*Y1); // Need a temp copy of Y1
00691     U_->apply(Y1temp, *Y1,mode);
00692     Y1->update(1.0, Y1temp, 1.0); // (account for implicit unit diagonal)
00693     if (isOverlapped_) {Y.doExport(*Y1,*L_->getGraph()->getExporter(), OverlapMode_);}
00694   }
00695   return(0);
00696 }
00697 
00698 //=============================================================================
00699 template<class MatrixType>
00700 typename Teuchos::ScalarTraits<typename MatrixType::scalar_type>::magnitudeType
00701 RILUK<MatrixType>::computeCondEst(Teuchos::ETransp mode) const {
00702 
00703   if (Condest_>=0.0) {
00704     return Condest_;
00705   }
00706   // Create a vector with all values equal to one
00707   Tpetra::Vector<scalar_type,local_ordinal_type,global_ordinal_type,node_type> Ones(U_->getDomainMap());
00708   Tpetra::Vector<scalar_type,local_ordinal_type,global_ordinal_type,node_type> OnesResult(L_->getRangeMap());
00709   Ones.putScalar(1.0);
00710 
00711   apply(Ones, OnesResult,mode); // Compute the effect of the solve on the vector of ones
00712   OnesResult.abs(OnesResult); // Make all values non-negative
00713   Teuchos::Array<magnitude_type> norms(1);
00714   OnesResult.normInf(norms());
00715   Condest_ = norms[0];
00716   return Condest_;
00717 }
00718 
00719 
00720 //=========================================================================
00721 template<class MatrixType>
00722 void RILUK<MatrixType>::generateXY(Teuchos::ETransp mode,
00723     const Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& Xin,
00724     const Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& Yin,
00725     Teuchos::RCP<const Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type> >& Xout,
00726     Teuchos::RCP<Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type> >& Yout) const {
00727 
00728   // Generate an X and Y suitable for performing Solve() and Multiply() methods
00729 
00730   TEUCHOS_TEST_FOR_EXCEPTION(Xin.getNumVectors()!=Yin.getNumVectors(), std::runtime_error,
00731        "Ifpack2::RILUK::GenerateXY ERROR: X and Y not the same size");
00732 
00733   //cout << "Xin = " << Xin << endl;
00734   Xout = Teuchos::rcp( (const Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type> *) &Xin, false );
00735   Yout = Teuchos::rcp( (Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type> *) &Yin, false );
00736   if (!isOverlapped_) return; // Nothing more to do
00737 
00738   if (isOverlapped_) {
00739     // Make sure the number of vectors in the multivector is the same as before.
00740     if (OverlapX_!=Teuchos::null) {
00741       if (OverlapX_->getNumVectors()!=Xin.getNumVectors()) {
00742         OverlapX_ = Teuchos::null;
00743         OverlapY_ = Teuchos::null;
00744       }
00745     }
00746     if (OverlapX_==Teuchos::null) { // Need to allocate space for overlap X and Y
00747       OverlapX_ = Teuchos::rcp( new Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>(U_->getColMap(), Xout->getNumVectors()) );
00748       OverlapY_ = Teuchos::rcp( new Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>(L_->getRowMap(), Yout->getNumVectors()) );
00749     }
00750     if (mode == Teuchos::NO_TRANS) {
00751       OverlapX_->doImport(*Xout,*U_->getGraph()->getImporter(), Tpetra::INSERT); // Import X values for solve
00752     }
00753     else {
00754       OverlapX_->doImport(*Xout,*L_->getGraph()->getExporter(), Tpetra::INSERT); // Import X values for solve
00755     }
00756     Xout = OverlapX_;
00757     Yout = OverlapY_; // Set pointers for Xout and Yout to point to overlap space
00758     //cout << "OverlapX_ = " << *OverlapX_ << endl;
00759   }
00760 }
00761 
00762 }//namespace Ifpack2
00763 
00764 #endif
00765 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends