Ifpack2 Templated Preconditioning Package Version 1.0
Ifpack2_ILUT_def.hpp
00001 /*@HEADER
00002 // ***********************************************************************
00003 // 
00004 //       Ifpack2: Tempated Object-Oriented Algebraic Preconditioner Package
00005 //                 Copyright (2009) 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 
00030 //-----------------------------------------------------
00031 // Ifpack2::ILUT is a translation of the Aztec ILUT
00032 // implementation. The Aztec ILUT implementation was
00033 // written by Ray Tuminaro.
00034 // See notes below, in the Ifpack2::ILUT::Compute method.
00035 // ABW.
00036 //------------------------------------------------------
00037 
00038 #ifndef IFPACK2_ILUT_DEF_HPP
00039 #define IFPACK2_ILUT_DEF_HPP
00040 
00041 namespace Ifpack2 {
00042 
00043 //Definitions for the ILUT methods:
00044 
00045 //==============================================================================
00046 template <class MatrixType>
00047 ILUT<MatrixType>::ILUT(const Teuchos::RCP<const Tpetra::RowMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >& A) :
00048   A_(A),
00049   Comm_(A->getRowMap()->getComm()),
00050   Athresh_(0.0),
00051   Rthresh_(1.0),
00052   RelaxValue_(0.0),
00053   LevelOfFill_(1.0),
00054   DropTolerance_(1e-12),
00055   Condest_(-1.0),
00056   IsInitialized_(false),
00057   IsComputed_(false),
00058   NumInitialize_(0),
00059   NumCompute_(0),
00060   NumApply_(0),
00061   InitializeTime_(0.0),
00062   ComputeTime_(0.0),
00063   ApplyTime_(0.0),
00064   Time_("Ifpack2::ILUT"),
00065   NumMyRows_(-1),
00066   NumGlobalNonzeros_(0)
00067 { 
00068   TEUCHOS_TEST_FOR_EXCEPTION(A_ == Teuchos::null, std::runtime_error, 
00069       Teuchos::typeName(*this) << "::ILUT(): input matrix reference was null.");
00070 }
00071 
00072 //==========================================================================
00073 template <class MatrixType>
00074 ILUT<MatrixType>::~ILUT() {
00075 }
00076 
00077 //==========================================================================
00078 template <class MatrixType>
00079 void ILUT<MatrixType>::setParameters(const Teuchos::ParameterList& params) {
00080   Ifpack2::getParameter(params, "fact: ilut level-of-fill", LevelOfFill_);
00081   TEUCHOS_TEST_FOR_EXCEPTION(LevelOfFill_ <= 0.0, std::runtime_error,
00082     "Ifpack2::ILUT::SetParameters ERROR, level-of-fill must be >= 0.");
00083 
00084   double tmp = -1;
00085   Ifpack2::getParameter(params, "fact: absolute threshold", tmp);
00086   if (tmp != -1) Athresh_ = tmp;
00087   tmp = -1;
00088   Ifpack2::getParameter(params, "fact: relative threshold", tmp);
00089   if (tmp != -1) Rthresh_ = tmp;
00090   tmp = -1;
00091   Ifpack2::getParameter(params, "fact: relax value", tmp);
00092   if (tmp != -1) RelaxValue_ = tmp;
00093   tmp = -1;
00094   Ifpack2::getParameter(params, "fact: drop tolerance", tmp);
00095   if (tmp != -1) DropTolerance_ = tmp;
00096 }
00097 
00098 //==========================================================================
00099 template <class MatrixType>
00100 const Teuchos::RCP<const Teuchos::Comm<int> > & 
00101 ILUT<MatrixType>::getComm() const{
00102   return(Comm_);
00103 }
00104 
00105 //==========================================================================
00106 template <class MatrixType>
00107 Teuchos::RCP<const Tpetra::RowMatrix<typename MatrixType::scalar_type,typename MatrixType::local_ordinal_type,typename MatrixType::global_ordinal_type,typename MatrixType::node_type> >
00108 ILUT<MatrixType>::getMatrix() const {
00109   return(A_);
00110 }
00111 
00112 //==========================================================================
00113 template <class MatrixType>
00114 const Teuchos::RCP<const Tpetra::Map<typename MatrixType::local_ordinal_type,typename MatrixType::global_ordinal_type,typename MatrixType::node_type> >&
00115 ILUT<MatrixType>::getDomainMap() const
00116 {
00117   return A_->getDomainMap();
00118 }
00119 
00120 //==========================================================================
00121 template <class MatrixType>
00122 const Teuchos::RCP<const Tpetra::Map<typename MatrixType::local_ordinal_type,typename MatrixType::global_ordinal_type,typename MatrixType::node_type> >&
00123 ILUT<MatrixType>::getRangeMap() const
00124 {
00125   return A_->getRangeMap();
00126 }
00127 
00128 //==============================================================================
00129 template <class MatrixType>
00130 bool ILUT<MatrixType>::hasTransposeApply() const {
00131   return true;
00132 }
00133 
00134 //==========================================================================
00135 template <class MatrixType>
00136 int ILUT<MatrixType>::getNumInitialize() const {
00137   return(NumInitialize_);
00138 }
00139 
00140 //==========================================================================
00141 template <class MatrixType>
00142 int ILUT<MatrixType>::getNumCompute() const {
00143   return(NumCompute_);
00144 }
00145 
00146 //==========================================================================
00147 template <class MatrixType>
00148 int ILUT<MatrixType>::getNumApply() const {
00149   return(NumApply_);
00150 }
00151 
00152 //==========================================================================
00153 template <class MatrixType>
00154 double ILUT<MatrixType>::getInitializeTime() const {
00155   return(InitializeTime_);
00156 }
00157 
00158 //==========================================================================
00159 template<class MatrixType>
00160 double ILUT<MatrixType>::getComputeTime() const {
00161   return(ComputeTime_);
00162 }
00163 
00164 //==========================================================================
00165 template<class MatrixType>
00166 double ILUT<MatrixType>::getApplyTime() const {
00167   return(ApplyTime_);
00168 }
00169 
00170 //==========================================================================
00171 template<class MatrixType>
00172 global_size_t ILUT<MatrixType>::getGlobalNumEntries() const { 
00173   return(L_->getGlobalNumEntries() + U_->getGlobalNumEntries());
00174 }
00175 
00176 //==========================================================================
00177 template<class MatrixType>
00178 size_t ILUT<MatrixType>::getNodeNumEntries() const {
00179   return(L_->getNodeNumEntries() + U_->getNodeNumEntries());
00180 }
00181 
00182 //=============================================================================
00183 template<class MatrixType>
00184 typename Teuchos::ScalarTraits<typename MatrixType::scalar_type>::magnitudeType
00185 ILUT<MatrixType>::computeCondEst(
00186                      CondestType CT,
00187                      LocalOrdinal MaxIters, 
00188                      magnitudeType Tol,
00189                      const Teuchos::Ptr<const Tpetra::RowMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > &matrix) {
00190   if (!isComputed()) { // cannot compute right now
00191     return(-1.0);
00192   }
00193   // NOTE: this is computing the *local* condest
00194   if (Condest_ == -1.0) {
00195     Condest_ = Ifpack2::Condest(*this, CT, MaxIters, Tol, matrix);
00196   }
00197   return(Condest_);
00198 }
00199 
00200 //==========================================================================
00201 template<class MatrixType>
00202 void ILUT<MatrixType>::initialize() {
00203   // clear any previous allocation
00204   IsInitialized_ = false;
00205   IsComputed_ = false;
00206   L_ = Teuchos::null;
00207   U_ = Teuchos::null;
00208 
00209   Time_.start(true);
00210 
00211   // check only in serial
00212   TEUCHOS_TEST_FOR_EXCEPTION(Comm_->getSize() == 1 && A_->getNodeNumRows() != A_->getNodeNumCols(), std::runtime_error, "Ifpack2::ILUT::Initialize ERROR, matrix must be square");
00213 
00214   NumMyRows_ = A_->getNodeNumRows();
00215 
00216   // nothing else to do here
00217   IsInitialized_ = true;
00218   ++NumInitialize_;
00219   Time_.stop();
00220   InitializeTime_ += Time_.totalElapsedTime();
00221 }
00222 
00223 template<typename Scalar>
00224 typename Teuchos::ScalarTraits<Scalar>::magnitudeType scalar_mag(const Scalar& s)
00225 {
00226   return Teuchos::ScalarTraits<Scalar>::magnitude(s);
00227 }
00228 
00229 //==========================================================================
00230 template<class MatrixType>
00231 void ILUT<MatrixType>::compute() {
00232   //--------------------------------------------------------------------------
00233   // Ifpack2::ILUT is a translation of the Aztec ILUT implementation. The Aztec
00234   // ILUT implementation was written by Ray Tuminaro.
00235   //
00236   // This isn't an exact translation of the Aztec ILUT algorithm, for the
00237   // following reasons:
00238   // 1. Minor differences result from the fact that Aztec factors a MSR format
00239   // matrix in place, while the code below factors an input CrsMatrix which
00240   // remains untouched and stores the resulting factors in separate L and U
00241   // CrsMatrix objects.
00242   // Also, the Aztec code begins by shifting the matrix pointers back
00243   // by one, and the pointer contents back by one, and then using 1-based
00244   // Fortran-style indexing in the algorithm. This Ifpack2 code uses C-style
00245   // 0-based indexing throughout.
00246   // 2. Aztec stores the inverse of the diagonal of U. This Ifpack2 code
00247   // stores the non-inverted diagonal in U.
00248   // The triangular solves (in Ifpack2::ILUT::apply()) are performed by
00249   // calling the Tpetra::CrsMatrix::solve method on the L and U objects, and
00250   // this requires U to contain the non-inverted diagonal.
00251   //
00252   // ABW.
00253   //--------------------------------------------------------------------------
00254 
00255   if (!isInitialized()) {
00256     initialize();
00257   }
00258 
00259   Time_.start(true);
00260 
00261   NumMyRows_ = A_->getNodeNumRows();
00262 
00263   L_ = Teuchos::rcp(new MatrixType(A_->getRowMap(), A_->getColMap(), 0));
00264   U_ = Teuchos::rcp(new MatrixType(A_->getRowMap(), A_->getColMap(), 0));
00265 
00266   TEUCHOS_TEST_FOR_EXCEPTION(L_ == Teuchos::null || U_ == Teuchos::null, std::runtime_error,
00267      "Ifpack2::ILUT::Compute ERROR, failed to allocate L_ or U_");
00268 
00269   const Scalar zero = Teuchos::ScalarTraits<Scalar>::zero();
00270   const Scalar one  = Teuchos::ScalarTraits<Scalar>::one();
00271 
00272   // CGB: note, this caching approach may not be necessary anymore
00273   // We will store ArrayView objects that are views of the rows of U, so that
00274   // we don't have to repeatedly retrieve the view for each row. These will
00275   // be populated row by row as the factorization proceeds.
00276   Teuchos::Array<Teuchos::ArrayView<const LocalOrdinal> > Uindices(NumMyRows_);
00277   Teuchos::Array<Teuchos::ArrayView<const Scalar> >       Ucoefs(NumMyRows_);
00278 
00279   // If this macro is defined, files containing the L and U factors
00280   // will be written. DON'T CHECK IN THE CODE WITH THIS MACRO ENABLED!!!
00281   // #define IFPACK2_WRITE_FACTORS
00282 #ifdef IFPACK2_WRITE_FACTORS
00283   std::ofstream ofsL("L.tif.mtx", std::ios::out);
00284   std::ofstream ofsU("U.tif.mtx", std::ios::out);
00285 #endif
00286 
00287   // Calculate how much fill will be allowed in addition to the space that
00288   // corresponds to the input matrix entries.
00289   double local_nnz = static_cast<double>(A_->getNodeNumEntries());
00290   double fill = ((getLevelOfFill()-1)*local_nnz)/(2*NumMyRows_);
00291 
00292   // std::ceil gives the smallest integer larger than the argument.
00293   // this may give a slightly different result than Aztec's fill value in
00294   // some cases.
00295   double fill_ceil=std::ceil(fill);
00296 
00297   typedef typename Teuchos::Array<LocalOrdinal>::size_type      Tsize_t;
00298 
00299   // Similarly to Aztec, we will allow the same amount of fill for each
00300   // row, half in L and half in U.
00301   Tsize_t fillL = static_cast<Tsize_t>(fill_ceil);
00302   Tsize_t fillU = static_cast<Tsize_t>(fill_ceil);
00303 
00304   Teuchos::Array<Scalar> InvDiagU(NumMyRows_,zero);
00305 
00306   Teuchos::Array<LocalOrdinal> tmp_idx;
00307   Teuchos::Array<Scalar> tmpv;
00308 
00309   enum { UNUSED, ORIG, FILL };
00310   LocalOrdinal max_col = NumMyRows_;
00311 
00312   Teuchos::Array<int> pattern(max_col, UNUSED);
00313   Teuchos::Array<Scalar> cur_row(max_col, zero);
00314   Teuchos::Array<magnitudeType> unorm(max_col);
00315   magnitudeType rownorm;
00316   Teuchos::Array<LocalOrdinal> L_cols_heap;
00317   Teuchos::Array<LocalOrdinal> U_cols;
00318   Teuchos::Array<LocalOrdinal> L_vals_heap;
00319   Teuchos::Array<LocalOrdinal> U_vals_heap;
00320 
00321   // A comparison object which will be used to create 'heaps' of indices
00322   // that are ordered according to the corresponding values in the
00323   // 'cur_row' array.
00324   greater_indirect<Scalar,LocalOrdinal> vals_comp(cur_row);
00325 
00326   // =================== //
00327   // start factorization //
00328   // =================== //
00329 
00330   for (LocalOrdinal row_i = 0 ; row_i < NumMyRows_ ; ++row_i) {
00331     Teuchos::ArrayView<const LocalOrdinal> ColIndicesA;
00332     Teuchos::ArrayView<const Scalar> ColValuesA;
00333 
00334     A_->getLocalRowView(row_i, ColIndicesA, ColValuesA);
00335     size_t RowNnz = ColIndicesA.size();
00336 
00337     // Always include the diagonal in the U factor. The value should get
00338     // set in the next loop below.
00339     U_cols.push_back(row_i);
00340     cur_row[row_i] = zero;
00341     pattern[row_i] = ORIG;
00342 
00343     Tsize_t L_cols_heaplen = 0;
00344     rownorm = (magnitudeType)0;
00345     for(size_t i=0; i<RowNnz; ++i) {
00346       if (ColIndicesA[i] < NumMyRows_) {
00347         if (ColIndicesA[i] < row_i) {
00348           add_to_heap(ColIndicesA[i], L_cols_heap, L_cols_heaplen);
00349         }
00350         else if (ColIndicesA[i] > row_i) {
00351           U_cols.push_back(ColIndicesA[i]);
00352         }
00353 
00354         cur_row[ColIndicesA[i]] = ColValuesA[i];
00355         pattern[ColIndicesA[i]] = ORIG;
00356         rownorm += scalar_mag(ColValuesA[i]);
00357       }
00358     }
00359 
00360     // Alter the diagonal according to the absolute-threshold and
00361     // relative-threshold values. If not set, those values default
00362     // to zero and one respectively.
00363     const typename Teuchos::ScalarTraits<Scalar>::magnitudeType rthresh = getRelativeThreshold();
00364     const Scalar& v = cur_row[row_i];
00365     cur_row[row_i] = (Scalar)(getAbsoluteThreshold()*IFPACK2_SGN(v)) + rthresh*v;
00366 
00367     Tsize_t orig_U_len = U_cols.size();
00368     RowNnz = L_cols_heap.size() + orig_U_len;
00369     rownorm = getDropTolerance() * rownorm/RowNnz;
00370 
00371     // The following while loop corresponds to the 'L30' goto's in Aztec.
00372     Tsize_t L_vals_heaplen = 0;
00373     while(L_cols_heaplen > 0) {
00374       LocalOrdinal row_k = L_cols_heap.front();
00375 
00376       Scalar multiplier = cur_row[row_k] * InvDiagU[row_k];
00377       cur_row[row_k] = multiplier;
00378       magnitudeType mag_mult = scalar_mag(multiplier);
00379       if (mag_mult*unorm[row_k] < rownorm) {
00380         pattern[row_k] = UNUSED;
00381         rm_heap_root(L_cols_heap, L_cols_heaplen);
00382         continue;
00383       }
00384       if (pattern[row_k] != ORIG) {
00385         if (L_vals_heaplen < fillL) {
00386           add_to_heap(row_k, L_vals_heap, L_vals_heaplen, vals_comp);
00387         }
00388         else if (L_vals_heaplen==0 ||
00389                  mag_mult < scalar_mag(cur_row[L_vals_heap.front()])) {
00390           pattern[row_k] = UNUSED;
00391           rm_heap_root(L_cols_heap, L_cols_heaplen);
00392           continue;
00393         }
00394         else {
00395           pattern[L_vals_heap.front()] = UNUSED;
00396           rm_heap_root(L_vals_heap, L_vals_heaplen, vals_comp);
00397           add_to_heap(row_k, L_vals_heap, L_vals_heaplen, vals_comp);
00398         }
00399       }
00400 
00401       /* Reduce current row */
00402 
00403       Teuchos::ArrayView<const LocalOrdinal>& ColIndicesU = Uindices[row_k];
00404       Teuchos::ArrayView<const Scalar>& ColValuesU = Ucoefs[row_k];
00405       Tsize_t ColNnzU = ColIndicesU.size();
00406 
00407       for(Tsize_t j=0; j<ColNnzU; ++j) {
00408         if (ColIndicesU[j] > row_k) {
00409           Scalar tmp = multiplier * ColValuesU[j];
00410           LocalOrdinal col_j = ColIndicesU[j];
00411           if (pattern[col_j] != UNUSED) {
00412             cur_row[col_j] -= tmp;
00413           }
00414           else if (scalar_mag(tmp) > rownorm) {
00415             cur_row[col_j] = -tmp;
00416             pattern[col_j] = FILL;
00417             if (col_j > row_i) {
00418               U_cols.push_back(col_j);
00419             }
00420             else {
00421               add_to_heap(col_j, L_cols_heap, L_cols_heaplen);
00422             }
00423           }
00424         }
00425       }
00426 
00427       rm_heap_root(L_cols_heap, L_cols_heaplen);
00428     }//end of while(L_cols_heaplen) loop
00429 
00430 
00431     // Put indices and values for L into arrays and then into the L_ matrix.
00432 
00433     //   first, the original entries from the L section of A:
00434     for(Tsize_t i=0; i<ColIndicesA.size(); ++i) {
00435       if (ColIndicesA[i] < row_i) {
00436         tmp_idx.push_back(ColIndicesA[i]);
00437         tmpv.push_back(cur_row[ColIndicesA[i]]);
00438         pattern[ColIndicesA[i]] = UNUSED;
00439       }
00440     }
00441 
00442     //   next, the L entries resulting from fill:
00443     for(Tsize_t j=0; j<L_vals_heaplen; ++j) {
00444       tmp_idx.push_back(L_vals_heap[j]);
00445       tmpv.push_back(cur_row[L_vals_heap[j]]);
00446       pattern[L_vals_heap[j]] = UNUSED;
00447     }
00448 
00449     // L has a one on the diagonal, but we don't explicitly store it.
00450     // If we don't store it, then the Tpetra/Kokkos kernel which performs
00451     // the triangular solve can assume a unit diagonal, take a short-cut
00452     // and perform faster.
00453 
00454     L_->insertLocalValues(row_i, tmp_idx(), tmpv());
00455 #ifdef IFPACK2_WRITE_FACTORS
00456     for(Tsize_t ii=0; ii<tmp_idx.size(); ++ii) {
00457       ofsL << row_i << " " << tmp_idx[ii] << " " << tmpv[ii] << std::endl;
00458     }
00459 #endif
00460 
00461     tmp_idx.clear();
00462     tmpv.clear();
00463 
00464     // Pick out the diagonal element, store its reciprocal.
00465     if (cur_row[row_i] == zero) {
00466       std::cerr << "Ifpack2::ILUT::Compute: zero pivot encountered! Replacing with rownorm and continuing...(You may need to set the parameter 'fact: absolute threshold'.)" << std::endl;
00467       cur_row[row_i] = rownorm;
00468     }
00469     InvDiagU[row_i] = one / cur_row[row_i];
00470 
00471     // Non-inverted diagonal is stored for U:
00472     tmp_idx.push_back(row_i);
00473     tmpv.push_back(cur_row[row_i]);
00474     unorm[row_i] = scalar_mag(cur_row[row_i]);
00475     pattern[row_i] = UNUSED;
00476 
00477     // Now put indices and values for U into arrays and then into the U_ matrix.
00478     // The first entry in U_cols is the diagonal, which we just handled, so we'll
00479     // start our loop at j=1.
00480 
00481     Tsize_t U_vals_heaplen = 0;
00482     for(Tsize_t j=1; j<U_cols.size(); ++j) {
00483       LocalOrdinal col = U_cols[j];
00484       if (pattern[col] != ORIG) {
00485         if (U_vals_heaplen < fillU) {
00486           add_to_heap(col, U_vals_heap, U_vals_heaplen, vals_comp);
00487         }
00488         else if (U_vals_heaplen!=0 && scalar_mag(cur_row[col]) >
00489                  scalar_mag(cur_row[U_vals_heap.front()])) {
00490           rm_heap_root(U_vals_heap, U_vals_heaplen, vals_comp);
00491           add_to_heap(col, U_vals_heap, U_vals_heaplen, vals_comp);
00492         }
00493       }
00494       else {
00495         tmp_idx.push_back(col);
00496         tmpv.push_back(cur_row[col]);
00497         unorm[row_i] += scalar_mag(cur_row[col]);
00498       }
00499       pattern[col] = UNUSED;
00500     }
00501 
00502     for(Tsize_t j=0; j<U_vals_heaplen; ++j) {
00503       tmp_idx.push_back(U_vals_heap[j]);
00504       tmpv.push_back(cur_row[U_vals_heap[j]]);
00505       unorm[row_i] += scalar_mag(cur_row[U_vals_heap[j]]);
00506     }
00507 
00508     unorm[row_i] /= (orig_U_len + U_vals_heaplen);
00509 
00510     U_->insertLocalValues(row_i, tmp_idx(), tmpv() );
00511 #ifdef IFPACK2_WRITE_FACTORS
00512     for(int ii=0; ii<tmp_idx.size(); ++ii) {
00513       ofsU <<row_i<< " " <<tmp_idx[ii]<< " " <<tmpv[ii]<< std::endl;
00514     }
00515 #endif
00516     tmp_idx.clear();
00517     tmpv.clear();
00518 
00519     U_->getLocalRowView(row_i, Uindices[row_i], Ucoefs[row_i] );
00520 
00521     L_cols_heap.clear();
00522     U_cols.clear();
00523     L_vals_heap.clear();
00524     U_vals_heap.clear();
00525   } // end of for(row_i) loop
00526 
00527   L_->fillComplete();
00528   U_->fillComplete();
00529 
00530   global_size_t MyNonzeros = L_->getGlobalNumEntries() + U_->getGlobalNumEntries();
00531   Teuchos::reduceAll(*Comm_,Teuchos::REDUCE_SUM,1,&MyNonzeros,&NumGlobalNonzeros_);
00532 
00533   IsComputed_ = true;
00534 
00535   ++NumCompute_;
00536   Time_.stop();
00537   ComputeTime_ += Time_.totalElapsedTime();
00538 }
00539 
00540 //==========================================================================
00541 template <class MatrixType>
00542 void ILUT<MatrixType>::apply(
00543            const Tpetra::MultiVector<typename MatrixType::scalar_type, typename MatrixType::local_ordinal_type, typename MatrixType::global_ordinal_type, typename MatrixType::node_type>& X,
00544                  Tpetra::MultiVector<typename MatrixType::scalar_type, typename MatrixType::local_ordinal_type, typename MatrixType::global_ordinal_type, typename MatrixType::node_type>& Y,
00545                  Teuchos::ETransp mode,
00546                typename MatrixType::scalar_type alpha,
00547                typename MatrixType::scalar_type beta) const
00548 {
00549   TEUCHOS_TEST_FOR_EXCEPTION(!isComputed(), std::runtime_error,
00550     "Ifpack2::ILUT::apply() ERROR, compute() hasn't been called yet.");
00551 
00552   TEUCHOS_TEST_FOR_EXCEPTION(X.getNumVectors() != Y.getNumVectors(), std::runtime_error,
00553     "Ifpack2::ILUT::apply() ERROR, X.getNumVectors() != Y.getNumVectors().");
00554 
00555   Time_.start(true);
00556 
00557   // If X and Y are pointing to the same memory location,
00558   // we need to create an auxiliary vector, Xcopy
00559   Teuchos::RCP< const Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > Xcopy;
00560   if (X.getLocalMV().getValues() == Y.getLocalMV().getValues())
00561     Xcopy = Teuchos::rcp( new Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>(X) );
00562   else
00563     Xcopy = Teuchos::rcp( &X, false );
00564 
00565   if (mode == Teuchos::NO_TRANS)
00566   {
00567     // solves LU Y = X
00568     L_->localSolve(*Xcopy,Y,Teuchos::NO_TRANS);
00569     U_->localSolve(Y,Y,Teuchos::NO_TRANS);
00570   }
00571   else
00572   {
00573     // solves U(trans) L(trans) Y = X
00574     U_->localSolve(*Xcopy,Y,mode);
00575     L_->localSolve(Y,Y,mode);
00576   }
00577 
00578   ++NumApply_;
00579   Time_.stop();
00580   ApplyTime_ += Time_.totalElapsedTime();
00581 }
00582 
00583 //=============================================================================
00584 template <class MatrixType>
00585 std::string ILUT<MatrixType>::description() const {
00586   std::ostringstream oss;
00587   oss << Teuchos::Describable::description();
00588   if (isInitialized()) {
00589     if (isComputed()) {
00590       oss << "{status = initialized, computed";
00591     }
00592     else {
00593       oss << "{status = initialized, not computed";
00594     }
00595   }
00596   else {
00597     oss << "{status = not initialized, not computed";
00598   }
00599   oss << ", global rows = " << A_->getGlobalNumRows()
00600       << ", global cols = " << A_->getGlobalNumCols()
00601       << "}";
00602   return oss.str();
00603 }
00604 
00605 //=============================================================================
00606 template <class MatrixType>
00607 void ILUT<MatrixType>::describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel) const {
00608   using std::endl;
00609   using std::setw;
00610   using Teuchos::VERB_DEFAULT;
00611   using Teuchos::VERB_NONE;
00612   using Teuchos::VERB_LOW;
00613   using Teuchos::VERB_MEDIUM;
00614   using Teuchos::VERB_HIGH;
00615   using Teuchos::VERB_EXTREME;
00616   Teuchos::EVerbosityLevel vl = verbLevel;
00617   if (vl == VERB_DEFAULT) vl = VERB_LOW;
00618   const int myImageID = Comm_->getRank();
00619   Teuchos::OSTab tab(out);
00620   //    none: print nothing
00621   //     low: print O(1) info from node 0
00622   //  medium: 
00623   //    high: 
00624   // extreme: 
00625   if (vl != VERB_NONE && myImageID == 0) {
00626     out << this->description() << endl;
00627     out << endl;
00628     out << "===============================================================================" << endl;
00629     out << "Level-of-fill      = " << getLevelOfFill()       << endl;
00630     out << "Absolute threshold = " << getAbsoluteThreshold() << endl;
00631     out << "Relative threshold = " << getRelativeThreshold() << endl;
00632     out << "Relax value        = " << getRelaxValue()        << endl;
00633     if   (Condest_ == -1.0) { out << "Condition number estimate       = N/A" << endl; }
00634     else                    { out << "Condition number estimate       = " << Condest_ << endl; }
00635     if (isComputed()) {
00636       out << "Number of nonzeros in A         = " << A_->getGlobalNumEntries() << endl;
00637       out << "Number of nonzeros in L + U     = " << getGlobalNumEntries()
00638           << " ( = " << 100.0 * (double)getGlobalNumEntries() / (double)A_->getGlobalNumEntries() << " % of A)" << endl;
00639       out << "nonzeros / rows                 = " << 1.0 * getGlobalNumEntries() / U_->getGlobalNumRows() << endl;
00640     }
00641     out << endl;
00642     out << "Phase           # calls    Total Time (s) " << endl;
00643     out << "------------    -------    ---------------" << endl;
00644     out << "initialize()    " << setw(7) << getNumInitialize() << "    " << setw(15) << getInitializeTime() << endl;
00645     out << "compute()       " << setw(7) << getNumCompute()    << "    " << setw(15) << getComputeTime()    << endl;
00646     out << "apply()         " << setw(7) << getNumApply()      << "    " << setw(15) << getApplyTime()      << endl;
00647     out << "==============================================================================="                << endl;
00648     out << endl;
00649   }
00650 }
00651 
00652 
00653 }//namespace Ifpack2
00654 
00655 #endif /* IFPACK2_ILUT_DEF_HPP */
00656 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends