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,LocalOrdinal,GlobalOrdinal,Node>(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,LocalOrdinal,GlobalOrdinal,Node>(Graph_->getL_Graph()->getRowMap()) );
00104   setAllocated(true);
00105 }
00106 
00107 //==========================================================================
00108 template<class MatrixType>
00109 void RILUK<MatrixType>::setParameters(const Teuchos::ParameterList& parameterlist) {
00110   Ifpack2::getParameter(parameterlist, "fact: iluk level-of-fill", LevelOfFill_);
00111   Ifpack2::getParameter(parameterlist, "fact: iluk level-of-overlap", LevelOfOverlap_);
00112   double tmp = -1;
00113   Ifpack2::getParameter(parameterlist, "fact: absolute threshold", tmp);
00114   if (tmp != -1) Athresh_ = tmp;
00115   tmp = -1;
00116   Ifpack2::getParameter(parameterlist, "fact: relative threshold", tmp);
00117   if (tmp != -1) Rthresh_ = tmp;
00118   tmp = -1;
00119   Ifpack2::getParameter(parameterlist, "fact: relax value", tmp);
00120   if (tmp != -1) RelaxValue_ = tmp;
00121 }
00122 
00123 //==========================================================================
00124 template<class MatrixType>
00125 void RILUK<MatrixType>::initialize() {
00126 
00127   if (Graph_ != Teuchos::null) return;
00128 
00129   Graph_ = Teuchos::rcp(new Ifpack2::IlukGraph<LocalOrdinal,GlobalOrdinal,Node>(A_->getCrsGraph(), LevelOfFill_, LevelOfOverlap_));
00130 
00131   Graph_->constructFilledGraph();
00132 
00133   setInitialized(true);
00134 
00135   if (!isAllocated()) allocate_L_and_U();
00136 
00137   if (isOverlapped_) {
00138     throw std::runtime_error("Error, Ifpack2::RILUK::initialize: overlapping not yet supported.");
00139 //    OverlapA = Teuchos::rcp( new MatrixType(Graph_->OverlapGraph()) );
00140 //    EPETRA_CHK_ERR(OverlapA->Import(A, *Graph_->OverlapImporter(), Insert));
00141 //    EPETRA_CHK_ERR(OverlapA->FillComplete());
00142   }
00143   else {
00144     initAllValues(*A_);
00145   }
00146 
00147   ++numInitialize_;
00148 }
00149 
00150 //==========================================================================
00151 template<class MatrixType>
00152 void RILUK<MatrixType>::initAllValues(const Tpetra::RowMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> & OverlapA) {
00153 
00154   size_t NumIn = 0, NumL = 0, NumU = 0;
00155   bool DiagFound = false;
00156   size_t NumNonzeroDiags = 0;
00157   size_t MaxNumEntries = OverlapA.getGlobalMaxNumRowEntries();
00158 
00159 
00160   Teuchos::Array<GlobalOrdinal> InI(MaxNumEntries); // Allocate temp space
00161   Teuchos::Array<GlobalOrdinal> LI(MaxNumEntries);
00162   Teuchos::Array<GlobalOrdinal> UI(MaxNumEntries);
00163   Teuchos::Array<Scalar> InV(MaxNumEntries);
00164   Teuchos::Array<Scalar> LV(MaxNumEntries);
00165   Teuchos::Array<Scalar> UV(MaxNumEntries);
00166 
00167   bool ReplaceValues = (L_->isStaticGraph() || L_->isLocallyIndexed()); // Check if values should be inserted or replaced
00168 
00169   L_->resumeFill();
00170   U_->resumeFill();
00171   if (ReplaceValues) {
00172     L_->setAllToScalar(0.0); // Zero out L and U matrices
00173     U_->setAllToScalar(0.0);
00174   }
00175 
00176   D_->putScalar(0.0); // Set diagonal values to zero
00177   Teuchos::ArrayRCP<Scalar> DV = D_->get1dViewNonConst(); // Get view of diagonal
00178 
00179   const Teuchos::RCP<const Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> >& rowMap =
00180     L_->getRowMap();
00181 
00182   // First we copy the user's matrix into L and U, regardless of fill level
00183 
00184   for (GlobalOrdinal i=rowMap->getMinGlobalIndex(); i<=rowMap->getMaxGlobalIndex(); ++i) {
00185     GlobalOrdinal global_row = i;
00186     LocalOrdinal local_row = rowMap->getLocalElement(global_row);
00187 
00188     OverlapA.getGlobalRowCopy(global_row, InI(), InV(), NumIn); // Get Values and Indices
00189     
00190     // Split into L and U (we don't assume that indices are ordered).
00191     
00192     NumL = 0; 
00193     NumU = 0; 
00194     DiagFound = false;
00195     
00196     for (size_t j=0; j< NumIn; j++) {
00197       GlobalOrdinal k = InI[j];
00198 
00199       if (k==i) {
00200         DiagFound = true;
00201         DV[local_row] += Rthresh_ * InV[j] + IFPACK2_SGN(InV[j]) * Athresh_; // Store perturbed diagonal in Tpetra::Vector D_
00202       }
00203 
00204       else if (k < 0) { // Out of range
00205         throw std::runtime_error("out of range (k<0) in Ifpack2::RILUK::initAllValues");
00206       }
00207 
00208       else if (k < i) {
00209         LI[NumL] = k;
00210         LV[NumL] = InV[j];
00211         NumL++;
00212       }
00213       else if (k <= rowMap->getMaxGlobalIndex()) {
00214         UI[NumU] = k;
00215         UV[NumU] = InV[j];
00216         NumU++;
00217       }
00218 //      else {
00219 //        throw std::runtime_error("out of range in Ifpack2::RILUK::initAllValues");
00220 //      }
00221     }
00222     
00223     // Check in things for this row of L and U
00224 
00225     if (DiagFound) ++NumNonzeroDiags;
00226     else DV[local_row] = Athresh_;
00227 
00228     if (NumL) {
00229       if (ReplaceValues) {
00230         L_->replaceGlobalValues(global_row, LI(0, NumL), LV(0,NumL));
00231       }
00232       else {
00233         L_->insertGlobalValues(global_row, LI(0,NumL), LV(0,NumL));
00234       }
00235     }
00236 
00237     if (NumU) {
00238       if (ReplaceValues) {
00239         U_->replaceGlobalValues(global_row, UI(0,NumU), UV(0,NumU));
00240       }
00241       else {
00242         U_->insertGlobalValues(global_row, UI(0,NumU), UV(0,NumU));
00243       }
00244     }
00245 
00246   }
00247 
00248   // The domain of L and the range of U are exactly their own row maps (there is no communication).
00249   // The domain of U and the range of L must be the same as those of the original matrix,
00250   // However if the original matrix is a VbrMatrix, these two latter maps are translation from
00251   // a block map to a point map.
00252   L_->fillComplete(L_->getColMap(), A_->getRangeMap());
00253   U_->fillComplete(A_->getDomainMap(), U_->getRowMap());
00254 
00255   // At this point L and U have the values of A in the structure of L and U, and diagonal vector D
00256 
00257   setInitialized(true);
00258   setFactored(false);
00259 
00260   size_t TotalNonzeroDiags = 0;
00261   Teuchos::reduceAll(*L_->getRowMap()->getComm(),Teuchos::REDUCE_SUM,
00262                      1,&NumNonzeroDiags,&TotalNonzeroDiags);
00263   NumMyDiagonals_ = NumNonzeroDiags;
00264   if (NumNonzeroDiags != U_->getNodeNumRows()) {
00265     throw std::runtime_error("Error in Ifpack2::RILUK::initAllValues, wrong number of diagonals.");
00266   }
00267 }
00268 
00269 //==========================================================================
00270 template<class MatrixType>
00271 void RILUK<MatrixType>::compute() {
00272 
00273   L_->resumeFill();
00274   U_->resumeFill();
00275 
00276   TEUCHOS_TEST_FOR_EXCEPTION(!isInitialized(), std::runtime_error,
00277       "Ifpack2::RILUK::compute() ERROR: isInitialized() must be true.");
00278   TEUCHOS_TEST_FOR_EXCEPTION(isComputed() == true, std::runtime_error,
00279       "Ifpack2::RILUK::compute() ERROR: Can't have already computed factors.");
00280 
00281   // MinMachNum should be officially defined, for now pick something a little 
00282   // bigger than IEEE underflow value
00283 
00284   Scalar MinDiagonalValue = Teuchos::ScalarTraits<Scalar>::rmin();
00285   Scalar MaxDiagonalValue = Teuchos::ScalarTraits<Scalar>::one()/MinDiagonalValue;
00286 
00287   size_t NumIn, NumL, NumU;
00288 
00289   // Get Maximun Row length
00290   size_t MaxNumEntries = L_->getNodeMaxNumRowEntries() + U_->getNodeMaxNumRowEntries() + 1;
00291 
00292   Teuchos::Array<LocalOrdinal> InI(MaxNumEntries); // Allocate temp space
00293   Teuchos::Array<Scalar> InV(MaxNumEntries);
00294   size_t num_cols = U_->getColMap()->getNodeNumElements();
00295   Teuchos::Array<int> colflag(num_cols);
00296 
00297   Teuchos::ArrayRCP<Scalar> DV = D_->get1dViewNonConst(); // Get view of diagonal
00298 
00299   size_t current_madds = 0; // We will count multiply-add as they happen
00300 
00301   // Now start the factorization.
00302 
00303   // Need some integer workspace and pointers
00304   size_t NumUU; 
00305   Teuchos::ArrayView<const LocalOrdinal> UUI;
00306   Teuchos::ArrayView<const Scalar> UUV;
00307   for (size_t j=0; j<num_cols; j++) colflag[j] = - 1;
00308 
00309   for(size_t i=0; i<L_->getNodeNumRows(); i++) {
00310     LocalOrdinal local_row = i;
00311 
00312  // Fill InV, InI with current row of L, D and U combined
00313 
00314     NumIn = MaxNumEntries;
00315     L_->getLocalRowCopy(local_row, InI(), InV(), NumL);
00316 
00317     InV[NumL] = DV[i]; // Put in diagonal
00318     InI[NumL] = local_row;
00319     
00320     U_->getLocalRowCopy(local_row, InI(NumL+1,MaxNumEntries-NumL-1), InV(NumL+1,MaxNumEntries-NumL-1), NumU);
00321     NumIn = NumL+NumU+1;
00322 
00323     // Set column flags
00324     for (size_t j=0; j<NumIn; j++) colflag[InI[j]] = j;
00325 
00326     Scalar diagmod = 0.0; // Off-diagonal accumulator
00327 
00328     for (size_t jj=0; jj<NumL; jj++) {
00329       LocalOrdinal j = InI[jj];
00330       Scalar multiplier = InV[jj]; // current_mults++;
00331 
00332       InV[jj] *= DV[j];
00333  
00334       U_->getLocalRowView(j, UUI, UUV); // View of row above
00335       NumUU = UUI.size();
00336 
00337       if (RelaxValue_==0.0) {
00338         for (size_t k=0; k<NumUU; k++) {
00339           int kk = colflag[UUI[k]];
00340           if (kk>-1) {
00341             InV[kk] -= multiplier*UUV[k];
00342             current_madds++;
00343           }
00344         }
00345       }
00346       else {
00347         for (size_t k=0; k<NumUU; k++) {
00348           int kk = colflag[UUI[k]];
00349           if (kk>-1) InV[kk] -= multiplier*UUV[k];
00350           else diagmod -= multiplier*UUV[k];
00351           current_madds++;
00352         }
00353       }
00354     }
00355     if (NumL) {
00356       L_->replaceLocalValues(local_row, InI(0,NumL), InV(0,NumL));  // Replace current row of L
00357     }
00358 
00359     DV[i] = InV[NumL]; // Extract Diagonal value
00360 
00361     if (RelaxValue_!=0.0) {
00362       DV[i] += RelaxValue_*diagmod; // Add off diagonal modifications
00363       // current_madds++;
00364     }
00365 
00366     if (Teuchos::ScalarTraits<Scalar>::magnitude(DV[i]) > Teuchos::ScalarTraits<Scalar>::magnitude(MaxDiagonalValue)) {
00367       if (Teuchos::ScalarTraits<Scalar>::real(DV[i]) < 0) DV[i] = - MinDiagonalValue;
00368       else DV[i] = MinDiagonalValue;
00369     }
00370     else
00371       DV[i] = Teuchos::ScalarTraits<Scalar>::one()/DV[i]; // Invert diagonal value
00372 
00373     for (size_t j=0; j<NumU; j++) InV[NumL+1+j] *= DV[i]; // Scale U by inverse of diagonal
00374 
00375     if (NumU) {
00376       U_->replaceLocalValues(local_row, InI(NumL+1,NumU), InV(NumL+1,NumU));  // Replace current row of L and U
00377     }
00378 
00379     // Reset column flags
00380     for (size_t j=0; j<NumIn; j++) colflag[InI[j]] = -1;
00381   }
00382 
00383   L_->fillComplete(L_->getColMap(), A_->getRangeMap());
00384   U_->fillComplete(A_->getDomainMap(), U_->getRowMap());
00385 
00386   // Validate that the L and U factors are actually lower and upper triangular
00387 
00388   if( !L_->isLowerTriangular() ) 
00389     throw std::runtime_error("Ifpack2::RILUK::compute() ERROR, L isn't lower triangular.");
00390   if( !U_->isUpperTriangular() ) 
00391     throw std::runtime_error("Ifpack2::RILUK::compute() ERROR, U isn't lower triangular.");
00392   
00393   // Add up flops
00394  
00395   double current_flops = 2 * current_madds;
00396   double total_flops = 0;
00397     
00398   // Get total madds across all PEs
00399   Teuchos::reduceAll(*L_->getRowMap()->getComm(),Teuchos::REDUCE_SUM,
00400                      1,&current_flops,&total_flops);
00401 
00402   // Now count the rest
00403   total_flops += (double) L_->getGlobalNumEntries(); // Accounts for multiplier above
00404   total_flops += (double) D_->getGlobalLength(); // Accounts for reciprocal of diagonal
00405   if (RelaxValue_!=0.0) total_flops += 2 * (double)D_->getGlobalLength(); // Accounts for relax update of diag
00406 
00407   //UpdateFlops(total_flops); // Update flop count
00408 
00409   setFactored(true);
00410   ++numCompute_;
00411 }
00412 
00413 //=============================================================================
00414 template<class MatrixType>
00415 void RILUK<MatrixType>::apply(
00416        const Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>& X, 
00417              Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>& Y,
00418              Teuchos::ETransp mode, Scalar alpha, Scalar beta) const
00419 {
00420   TEUCHOS_TEST_FOR_EXCEPTION(!isComputed(), std::runtime_error,
00421     "Ifpack2::RILUK::apply() ERROR, compute() hasn't been called yet.");
00422 
00423 //
00424 // This function finds Y such that
00425 // LDU Y = X, or
00426 // U(trans) D L(trans) Y = X
00427 // for multiple RHS
00428 //
00429 
00430   // First generate X and Y as needed for this function
00431   Teuchos::RCP<const Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > X1;
00432   Teuchos::RCP<Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > Y1;
00433   generateXY(mode, X, Y, X1, Y1);
00434 
00435   Scalar one = Teuchos::ScalarTraits<Scalar>::one();
00436   Scalar zero = Teuchos::ScalarTraits<Scalar>::zero();
00437 
00438   if (mode == Teuchos::NO_TRANS) {
00439 
00440     L_->localSolve(*X1, *Y1,mode);
00441     Y1->elementWiseMultiply(one, *D_, *Y1, zero); // y = D*y (D_ has inverse of diagonal)
00442     U_->localSolve(*Y1, *Y1,mode); // Solve Uy = y
00443     if (isOverlapped_) {
00444       // Export computed Y values if needed
00445       Y.doExport(*Y1,*L_->getGraph()->getExporter(), OverlapMode_);
00446     }
00447   }
00448   else {
00449     U_->localSolve(*X1, *Y1,mode); // Solve Uy = y
00450     Y1->elementWiseMultiply(one, *D_, *Y1, zero); // y = D*y (D_ has inverse of diagonal)
00451     L_->localSolve(*Y1, *Y1,mode);
00452     if (isOverlapped_) {Y.doExport(*Y1,*U_->getGraph()->getImporter(), OverlapMode_);} // Export computed Y values if needed
00453   } 
00454 
00455   ++numApply_;
00456 }
00457 
00458 //=============================================================================
00459 template<class MatrixType>
00460 int RILUK<MatrixType>::Multiply(const Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>& X, 
00461             Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>& Y,
00462             Teuchos::ETransp mode) const {
00463 //
00464 // This function finds X such that LDU Y = X or U(trans) D L(trans) Y = X for multiple RHS
00465 //
00466     
00467   // First generate X and Y as needed for this function
00468   Teuchos::RCP<const Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > X1;
00469   Teuchos::RCP<Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > Y1;
00470   generateXY(mode, X, Y, X1, Y1);
00471 
00472 //  Epetra_Flops * counter = this->GetFlopCounter();
00473 //  if (counter!=0) {
00474 //    L_->SetFlopCounter(*counter);
00475 //    Y1->SetFlopCounter(*counter);
00476 //    U_->SetFlopCounter(*counter);
00477 //  }
00478 
00479   if (!mode == Teuchos::NO_TRANS) {
00480     U_->apply(*X1, *Y1,mode); // 
00481     Y1->update(1.0, *X1, 1.0); // Y1 = Y1 + X1 (account for implicit unit diagonal)
00482     Y1->elementWiseMultiply(1.0, *D_, *Y1, 0.0); // y = D*y (D_ has inverse of diagonal)
00483     Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> Y1temp(*Y1); // Need a temp copy of Y1
00484     L_->apply(Y1temp, *Y1,mode);
00485     Y1->update(1.0, Y1temp, 1.0); // (account for implicit unit diagonal)
00486     if (isOverlapped_) {Y.doExport(*Y1,*L_->getGraph()->getExporter(), OverlapMode_);} // Export computed Y values if needed
00487   }
00488   else {
00489 
00490     L_->apply(*X1, *Y1,mode);
00491     Y1->update(1.0, *X1, 1.0); // Y1 = Y1 + X1 (account for implicit unit diagonal)
00492     Y1->elementWiseMultiply(1, *D_, *Y1, 0); // y = D*y (D_ has inverse of diagonal)
00493     Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> Y1temp(*Y1); // Need a temp copy of Y1
00494     U_->apply(Y1temp, *Y1,mode);
00495     Y1->update(1.0, Y1temp, 1.0); // (account for implicit unit diagonal)
00496     if (isOverlapped_) {Y.doExport(*Y1,*L_->getGraph()->getExporter(), OverlapMode_);}
00497   } 
00498   return(0);
00499 }
00500 
00501 //=============================================================================
00502 template<class MatrixType>
00503 typename Teuchos::ScalarTraits<typename MatrixType::scalar_type>::magnitudeType
00504 RILUK<MatrixType>::computeCondEst(Teuchos::ETransp mode) const {
00505 
00506   if (Condest_>=0.0) {
00507     return Condest_;
00508   }
00509   // Create a vector with all values equal to one
00510   Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> Ones(U_->getDomainMap());
00511   Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> OnesResult(L_->getRangeMap());
00512   Ones.putScalar(1.0);
00513 
00514   apply(Ones, OnesResult,mode); // Compute the effect of the solve on the vector of ones
00515   OnesResult.abs(OnesResult); // Make all values non-negative
00516   Teuchos::Array<magnitudeType> norms(1);
00517   OnesResult.normInf(norms());
00518   Condest_ = norms[0];
00519   return Condest_;
00520 }
00521 
00522 
00523 //=========================================================================
00524 template<class MatrixType>
00525 void RILUK<MatrixType>::generateXY(Teuchos::ETransp mode, 
00526     const Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>& Xin,
00527     const Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>& Yin,
00528     Teuchos::RCP<const Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> >& Xout, 
00529     Teuchos::RCP<Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> >& Yout) const {
00530 
00531   // Generate an X and Y suitable for performing Solve() and Multiply() methods
00532 
00533   TEUCHOS_TEST_FOR_EXCEPTION(Xin.getNumVectors()!=Yin.getNumVectors(), std::runtime_error,
00534        "Ifpack2::RILUK::GenerateXY ERROR: X and Y not the same size");
00535 
00536   //cout << "Xin = " << Xin << endl;
00537   Xout = Teuchos::rcp( (const Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> *) &Xin, false );
00538   Yout = Teuchos::rcp( (Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> *) &Yin, false );
00539   if (!isOverlapped_) return; // Nothing more to do
00540 
00541   if (isOverlapped_) {
00542     // Make sure the number of vectors in the multivector is the same as before.
00543     if (OverlapX_!=Teuchos::null) {
00544       if (OverlapX_->getNumVectors()!=Xin.getNumVectors()) {
00545         OverlapX_ = Teuchos::null;
00546         OverlapY_ = Teuchos::null;
00547       }
00548     }
00549     if (OverlapX_==Teuchos::null) { // Need to allocate space for overlap X and Y
00550       OverlapX_ = Teuchos::rcp( new Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>(U_->getColMap(), Xout->getNumVectors()) );
00551       OverlapY_ = Teuchos::rcp( new Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>(L_->getRowMap(), Yout->getNumVectors()) );
00552     }
00553     if (mode == Teuchos::NO_TRANS) {
00554       OverlapX_->doImport(*Xout,*U_->getGraph()->getImporter(), Tpetra::INSERT); // Import X values for solve
00555     }
00556     else {
00557       OverlapX_->doImport(*Xout,*L_->getGraph()->getExporter(), Tpetra::INSERT); // Import X values for solve
00558     }
00559     Xout = OverlapX_;
00560     Yout = OverlapY_; // Set pointers for Xout and Yout to point to overlap space
00561     //cout << "OverlapX_ = " << *OverlapX_ << endl;
00562   }
00563 }
00564 
00565 }//namespace Ifpack2
00566 
00567 #endif
00568 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends