Ifpack2 Templated Preconditioning Package Version 1.0
Ifpack2_RILUK_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 // Redistribution and use in source and binary forms, with or without
00011 // modification, are permitted provided that the following conditions are
00012 // met:
00013 //
00014 // 1. Redistributions of source code must retain the above copyright
00015 // notice, this list of conditions and the following disclaimer.
00016 //
00017 // 2. Redistributions in binary form must reproduce the above copyright
00018 // notice, this list of conditions and the following disclaimer in the
00019 // documentation and/or other materials provided with the distribution.
00020 //
00021 // 3. Neither the name of the Corporation nor the names of the
00022 // contributors may be used to endorse or promote products derived from
00023 // this software without specific prior written permission.
00024 //
00025 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00026 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00027 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00028 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00029 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00030 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00031 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00032 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00033 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00034 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00035 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00036 //
00037 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
00038 //
00039 // ***********************************************************************
00040 //@HEADER
00041 */
00042 
00043 #ifndef IFPACK2_CRSRILUK_DEF_HPP
00044 #define IFPACK2_CRSRILUK_DEF_HPP
00045 
00046 #include <Ifpack2_LocalFilter.hpp>
00047 
00048 namespace Ifpack2 {
00049 
00050 template<class MatrixType>
00051 RILUK<MatrixType>::RILUK (const Teuchos::RCP<const row_matrix_type>& Matrix_in)
00052   : A_ (Matrix_in),
00053     LevelOfFill_ (0),
00054     isAllocated_ (false),
00055     isInitialized_ (false),
00056     isComputed_ (false),
00057     numInitialize_ (0),
00058     numCompute_ (0),
00059     numApply_ (0),
00060     initializeTime_ (0.0),
00061     computeTime_ (0.0),
00062     applyTime_ (0.0),
00063     RelaxValue_ (Teuchos::ScalarTraits<magnitude_type>::zero ()),
00064     Athresh_ (Teuchos::ScalarTraits<magnitude_type>::zero ()),
00065     Rthresh_ (Teuchos::ScalarTraits<magnitude_type>::one ()),
00066     Condest_ (-Teuchos::ScalarTraits<magnitude_type>::one ())
00067 {}
00068 
00069 
00070 template<class MatrixType>
00071 RILUK<MatrixType>::RILUK (const Teuchos::RCP<const crs_matrix_type>& Matrix_in)
00072   : A_ (Matrix_in),
00073     LevelOfFill_ (0),
00074     isAllocated_ (false),
00075     isInitialized_ (false),
00076     isComputed_ (false),
00077     numInitialize_ (0),
00078     numCompute_ (0),
00079     numApply_ (0),
00080     initializeTime_ (0.0),
00081     computeTime_ (0.0),
00082     applyTime_ (0.0),
00083     RelaxValue_ (Teuchos::ScalarTraits<magnitude_type>::zero ()),
00084     Athresh_ (Teuchos::ScalarTraits<magnitude_type>::zero ()),
00085     Rthresh_ (Teuchos::ScalarTraits<magnitude_type>::one ()),
00086     Condest_ (-Teuchos::ScalarTraits<magnitude_type>::one ())
00087 {}
00088 
00089 
00090 template<class MatrixType>
00091 RILUK<MatrixType>::~RILUK() {}
00092 
00093 
00094 template<class MatrixType>
00095 void
00096 RILUK<MatrixType>::setMatrix (const Teuchos::RCP<const row_matrix_type>& A)
00097 {
00098   // It's legal for A to be null; in that case, you may not call
00099   // initialize() until calling setMatrix() with a nonnull input.
00100   // Regardless, setting the matrix invalidates any previous
00101   // factorization.
00102   if (A.getRawPtr () != A_.getRawPtr ()) {
00103     isAllocated_ = false;
00104     isInitialized_ = false;
00105     isComputed_ = false;
00106     A_local_crs_ = Teuchos::null;
00107     Graph_ = Teuchos::null;
00108     L_ = Teuchos::null;
00109     U_ = Teuchos::null;
00110     D_ = Teuchos::null;
00111     Condest_ = -Teuchos::ScalarTraits<magnitude_type>::one();
00112     A_ = A;
00113   }
00114 }
00115 
00116 
00117 
00118 template<class MatrixType>
00119 const typename RILUK<MatrixType>::crs_matrix_type&
00120 RILUK<MatrixType>::getL () const
00121 {
00122   TEUCHOS_TEST_FOR_EXCEPTION(
00123     L_.is_null (), std::runtime_error, "Ifpack2::RILUK::getL: The L factor "
00124     "is null.  Please call initialize() (and preferably also compute()) "
00125     "before calling this method.  If the input matrix has not yet been set, "
00126     "you must first call setMatrix() with a nonnull input matrix before you "
00127     "may call initialize() or compute().");
00128   return *L_;
00129 }
00130 
00131 
00132 template<class MatrixType>
00133 const Tpetra::Vector<typename RILUK<MatrixType>::scalar_type,
00134                      typename RILUK<MatrixType>::local_ordinal_type,
00135                      typename RILUK<MatrixType>::global_ordinal_type,
00136                      typename RILUK<MatrixType>::node_type>&
00137 RILUK<MatrixType>::getD () const
00138 {
00139   TEUCHOS_TEST_FOR_EXCEPTION(
00140     D_.is_null (), std::runtime_error, "Ifpack2::RILUK::getD: The D factor "
00141     "(of diagonal entries) is null.  Please call initialize() (and "
00142     "preferably also compute()) before calling this method.  If the input "
00143     "matrix has not yet been set, you must first call setMatrix() with a "
00144     "nonnull input matrix before you may call initialize() or compute().");
00145   return *D_;
00146 }
00147 
00148 
00149 template<class MatrixType>
00150 const typename RILUK<MatrixType>::crs_matrix_type&
00151 RILUK<MatrixType>::getU () const
00152 {
00153   TEUCHOS_TEST_FOR_EXCEPTION(
00154     U_.is_null (), std::runtime_error, "Ifpack2::RILUK::getU: The U factor "
00155     "is null.  Please call initialize() (and preferably also compute()) "
00156     "before calling this method.  If the input matrix has not yet been set, "
00157     "you must first call setMatrix() with a nonnull input matrix before you "
00158     "may call initialize() or compute().");
00159   return *U_;
00160 }
00161 
00162 
00163 template<class MatrixType>
00164 Teuchos::RCP<const Tpetra::Map<typename RILUK<MatrixType>::local_ordinal_type,
00165                                typename RILUK<MatrixType>::global_ordinal_type,
00166                                typename RILUK<MatrixType>::node_type> >
00167 RILUK<MatrixType>::getDomainMap () const {
00168   TEUCHOS_TEST_FOR_EXCEPTION(
00169     A_.is_null (), std::runtime_error, "Ifpack2::RILUK::getDomainMap: "
00170     "The matrix is null.  Please call setMatrix() with a nonnull input "
00171     "before calling this method.");
00172 
00173   // FIXME (mfh 25 Jan 2014) Shouldn't this just come from A_?
00174   TEUCHOS_TEST_FOR_EXCEPTION(
00175     Graph_.is_null (), std::runtime_error, "Ifpack2::RILUK::getDomainMap: "
00176     "The computed graph is null.  Please call initialize() before calling "
00177     "this method.");
00178   return Graph_->getL_Graph ()->getDomainMap ();
00179 }
00180 
00181 
00182 template<class MatrixType>
00183 Teuchos::RCP<const Tpetra::Map<typename RILUK<MatrixType>::local_ordinal_type,
00184                                typename RILUK<MatrixType>::global_ordinal_type,
00185                                typename RILUK<MatrixType>::node_type> >
00186 RILUK<MatrixType>::getRangeMap () const {
00187   TEUCHOS_TEST_FOR_EXCEPTION(
00188     A_.is_null (), std::runtime_error, "Ifpack2::RILUK::getRangeMap: "
00189     "The matrix is null.  Please call setMatrix() with a nonnull input "
00190     "before calling this method.");
00191 
00192   // FIXME (mfh 25 Jan 2014) Shouldn't this just come from A_?
00193   TEUCHOS_TEST_FOR_EXCEPTION(
00194     Graph_.is_null (), std::runtime_error, "Ifpack2::RILUK::getRangeMap: "
00195     "The computed graph is null.  Please call initialize() before calling "
00196     "this method.");
00197   return Graph_->getL_Graph ()->getRangeMap ();
00198 }
00199 
00200 
00201 template<class MatrixType>
00202 void RILUK<MatrixType>::allocate_L_and_U ()
00203 {
00204   using Teuchos::null;
00205   using Teuchos::rcp;
00206 
00207   if (! isAllocated_) {
00208     // Deallocate any existing storage.  This avoids storing 2x
00209     // memory, since RCP op= does not deallocate until after the
00210     // assignment.
00211     L_ = null;
00212     U_ = null;
00213     D_ = null;
00214 
00215     // Allocate Matrix using ILUK graphs
00216     L_ = rcp (new crs_matrix_type (Graph_->getL_Graph ()));
00217     U_ = rcp (new crs_matrix_type (Graph_->getU_Graph ()));
00218     L_->setAllToScalar (STS::zero ()); // Zero out L and U matrices
00219     U_->setAllToScalar (STS::zero ());
00220 
00221     // FIXME (mfh 24 Jan 2014) This assumes domain == range Map for L and U.
00222     L_->fillComplete ();
00223     U_->fillComplete ();
00224     D_ = rcp (new vec_type (Graph_->getL_Graph ()->getRowMap ()));
00225   }
00226   isAllocated_ = true;
00227 }
00228 
00229 
00230 template<class MatrixType>
00231 void
00232 RILUK<MatrixType>::
00233 setParameters (const Teuchos::ParameterList& params)
00234 {
00235   using Teuchos::as;
00236   using Teuchos::Exceptions::InvalidParameterName;
00237   using Teuchos::Exceptions::InvalidParameterType;
00238 
00239   // Default values of the various parameters.
00240   int fillLevel = 0;
00241   magnitude_type absThresh = STM::zero ();
00242   magnitude_type relThresh = STM::one ();
00243   magnitude_type relaxValue = STM::zero ();
00244 
00245   //
00246   // "fact: iluk level-of-fill" parsing is more complicated, because
00247   // we want to allow as many types as make sense.  int is the native
00248   // type, but we also want to accept magnitude_type (for
00249   // compatibility with ILUT) and double (for backwards compatibilty
00250   // with ILUT).
00251   //
00252 
00253   bool gotFillLevel = false;
00254   try {
00255     fillLevel = params.get<int> ("fact: iluk level-of-fill");
00256     gotFillLevel = true;
00257   }
00258   catch (InvalidParameterType&) {
00259     // Throwing again in the catch block would just unwind the stack.
00260     // Instead, we do nothing here, and check the Boolean outside to
00261     // see if we got the value.
00262   }
00263   catch (InvalidParameterName&) {
00264     gotFillLevel = true; // Accept the default value.
00265   }
00266 
00267   if (! gotFillLevel) {
00268     try {
00269       // Try magnitude_type, for compatibility with ILUT.
00270       // The cast from magnitude_type to int must succeed.
00271       fillLevel = as<int> (params.get<magnitude_type> ("fact: iluk level-of-fill"));
00272       gotFillLevel = true;
00273     }
00274     catch (InvalidParameterType&) {
00275       // Try double next.
00276     }
00277     // Don't catch InvalidParameterName here; we've already done that above.
00278   }
00279 
00280   if (! gotFillLevel) {
00281     try {
00282       // Try double, for compatibility with ILUT.
00283       // The cast from double to int must succeed.
00284       fillLevel = as<int> (params.get<double> ("fact: iluk level-of-fill"));
00285       gotFillLevel = true;
00286     }
00287     catch (InvalidParameterType& e) {
00288       // We're out of options.  The user gave us the parameter, but it
00289       // doesn't have the right type.  The best thing for us to do in
00290       // that case is to throw, telling the user to use the right
00291       // type.
00292       throw e;
00293     }
00294     // Don't catch InvalidParameterName here; we've already done that above.
00295   }
00296 
00297   TEUCHOS_TEST_FOR_EXCEPTION(
00298     ! gotFillLevel,
00299     std::logic_error,
00300     "Ifpack2::RILUK::setParameters: We should never get here!  "
00301     "The method should either have read the \"fact: iluk level-of-fill\"  "
00302     "parameter by this point, or have thrown an exception.  "
00303     "Please let the Ifpack2 developers know about this bug.");
00304 
00305   //
00306   // For the other parameters, we prefer magnitude_type, but allow
00307   // double for backwards compatibility.
00308   //
00309 
00310   try {
00311     absThresh = params.get<magnitude_type> ("fact: absolute threshold");
00312   }
00313   catch (InvalidParameterType&) {
00314     // Try double, for backwards compatibility.
00315     // The cast from double to magnitude_type must succeed.
00316     absThresh = as<magnitude_type> (params.get<double> ("fact: absolute threshold"));
00317   }
00318   catch (InvalidParameterName&) {
00319     // Accept the default value.
00320   }
00321 
00322   try {
00323     relThresh = params.get<magnitude_type> ("fact: relative threshold");
00324   }
00325   catch (InvalidParameterType&) {
00326     // Try double, for backwards compatibility.
00327     // The cast from double to magnitude_type must succeed.
00328     relThresh = as<magnitude_type> (params.get<double> ("fact: relative threshold"));
00329   }
00330   catch (InvalidParameterName&) {
00331     // Accept the default value.
00332   }
00333 
00334   try {
00335     relaxValue = params.get<magnitude_type> ("fact: relax value");
00336   }
00337   catch (InvalidParameterType&) {
00338     // Try double, for backwards compatibility.
00339     // The cast from double to magnitude_type must succeed.
00340     relaxValue = as<magnitude_type> (params.get<double> ("fact: relax value"));
00341   }
00342   catch (InvalidParameterName&) {
00343     // Accept the default value.
00344   }
00345 
00346   // "Commit" the values only after validating all of them.  This
00347   // ensures that there are no side effects if this routine throws an
00348   // exception.
00349 
00350   LevelOfFill_ = fillLevel;
00351 
00352   // mfh 28 Nov 2012: The previous code would not assign Athresh_,
00353   // Rthresh_, or RelaxValue_, if the read-in value was -1.  I don't
00354   // know if keeping this behavior is correct, but I'll keep it just
00355   // so as not to change previous behavior.
00356 
00357   if (absThresh != -STM::one ()) {
00358     Athresh_ = absThresh;
00359   }
00360   if (relThresh != -STM::one ()) {
00361     Rthresh_ = relThresh;
00362   }
00363   if (relaxValue != -STM::one ()) {
00364     RelaxValue_ = relaxValue;
00365   }
00366 }
00367 
00368 
00369 template<class MatrixType>
00370 Teuchos::RCP<const typename RILUK<MatrixType>::row_matrix_type>
00371 RILUK<MatrixType>::getMatrix () const {
00372   return Teuchos::rcp_implicit_cast<const row_matrix_type> (A_);
00373 }
00374 
00375 
00376 template<class MatrixType>
00377 Teuchos::RCP<const typename RILUK<MatrixType>::crs_matrix_type>
00378 RILUK<MatrixType>::getCrsMatrix () const {
00379   return Teuchos::rcp_dynamic_cast<const crs_matrix_type> (A_, true);
00380 }
00381 
00382 
00383 template<class MatrixType>
00384 Teuchos::RCP<const typename RILUK<MatrixType>::row_matrix_type>
00385 RILUK<MatrixType>::makeLocalFilter (const Teuchos::RCP<const row_matrix_type>& A)
00386 {
00387   using Teuchos::RCP;
00388   using Teuchos::rcp;
00389   using Teuchos::rcp_dynamic_cast;
00390   using Teuchos::rcp_implicit_cast;
00391 
00392   // If A_'s communicator only has one process, or if its column and
00393   // row Maps are the same, then it is already local, so use it
00394   // directly.
00395   if (A->getRowMap ()->getComm ()->getSize () == 1 ||
00396       A->getRowMap ()->isSameAs (* (A->getColMap ()))) {
00397     return A;
00398   }
00399 
00400   // If A_ is already a LocalFilter, then use it directly.  This
00401   // should be the case if RILUK is being used through
00402   // AdditiveSchwarz, for example.  There are (unfortunately) two
00403   // kinds of LocalFilter, depending on the template parameter, so we
00404   // have to test for both.
00405   RCP<const LocalFilter<row_matrix_type> > A_lf_r =
00406     rcp_dynamic_cast<const LocalFilter<row_matrix_type> > (A);
00407   if (! A_lf_r.is_null ()) {
00408     return rcp_implicit_cast<const row_matrix_type> (A_lf_r);
00409   }
00410   RCP<const LocalFilter<crs_matrix_type> > A_lf_c =
00411     rcp_dynamic_cast<const LocalFilter<crs_matrix_type> > (A);
00412   if (! A_lf_c.is_null ()) {
00413     return rcp_implicit_cast<const row_matrix_type> (A_lf_c);
00414   }
00415 
00416   // A_'s communicator has more than one process, its row Map and
00417   // its column Map differ, and A_ is not a LocalFilter.  Thus, we
00418   // have to wrap it in a LocalFilter.
00419   return rcp (new LocalFilter<row_matrix_type> (A));
00420 }
00421 
00422 
00423 template<class MatrixType>
00424 void RILUK<MatrixType>::initialize ()
00425 {
00426   using Teuchos::RCP;
00427   using Teuchos::rcp;
00428   using Teuchos::rcp_const_cast;
00429   using Teuchos::rcp_dynamic_cast;
00430   using Teuchos::rcp_implicit_cast;
00431   typedef Tpetra::CrsGraph<local_ordinal_type,
00432                            global_ordinal_type,
00433                            node_type> crs_graph_type;
00434   TEUCHOS_TEST_FOR_EXCEPTION(
00435     A_.is_null (), std::runtime_error, "Ifpack2::RILUK::initialize: "
00436     "The matrix is null.  Please call setMatrix() with a nonnull input "
00437     "before calling this method.");
00438 
00439   Teuchos::Time timer ("RILUK::initialize");
00440   { // Start timing
00441     Teuchos::TimeMonitor timeMon (timer);
00442 
00443     // Calling initialize() means that the user asserts that the graph
00444     // of the sparse matrix may have changed.  We must not just reuse
00445     // the previous graph in that case.
00446     //
00447     // Regarding setting isAllocated_ to false: Eventually, we may want
00448     // some kind of clever memory reuse strategy, but it's always
00449     // correct just to blow everything away and start over.
00450     isInitialized_ = false;
00451     isAllocated_ = false;
00452     isComputed_ = false;
00453     Graph_ = Teuchos::null;
00454 
00455     RCP<const row_matrix_type> A_local = makeLocalFilter (A_);
00456     TEUCHOS_TEST_FOR_EXCEPTION(
00457       A_local.is_null (), std::logic_error, "Ifpack2::RILUK::initialize: "
00458       "makeLocalFilter returned null; it failed to compute A_local.  "
00459       "Please report this bug to the Ifpack2 developers.");
00460 
00461     // FIXME (mfh 24 Jan 2014, 26 Mar 2014) It would be more efficient
00462     // to rewrite RILUK so that it works with any RowMatrix input, not
00463     // just CrsMatrix.  (That would require rewriting IlukGraph to
00464     // handle a Tpetra::RowGraph.)  However, to make it work for now,
00465     // we just copy the input matrix if it's not a CrsMatrix.
00466     {
00467       RCP<const crs_matrix_type> A_local_crs =
00468         rcp_dynamic_cast<const crs_matrix_type> (A_local);
00469       if (A_local_crs.is_null ()) {
00470         // FIXME (mfh 24 Jan 2014) It would be smarter to count up the
00471         // number of elements in each row of A_local, so that we can
00472         // create A_local_crs_nc using static profile.  The code below is
00473         // correct but potentially slow.
00474         RCP<crs_matrix_type> A_local_crs_nc =
00475           rcp (new crs_matrix_type (A_local->getRowMap (),
00476                                     A_local->getColMap (), 0));
00477         // FIXME (mfh 24 Jan 2014) This Import approach will only work
00478         // if A_ has a one-to-one row Map.  This is generally the case
00479         // with matrices given to Ifpack2.
00480         //
00481         // Source and destination Maps are the same in this case.
00482         // That way, the Import just implements a copy.
00483         typedef Tpetra::Import<local_ordinal_type, global_ordinal_type,
00484           node_type> import_type;
00485         import_type import (A_local->getRowMap (), A_local->getRowMap ());
00486         A_local_crs_nc->doImport (*A_local, import, Tpetra::REPLACE);
00487         A_local_crs_nc->fillComplete (A_local->getDomainMap (), A_local->getRangeMap ());
00488         A_local_crs = rcp_const_cast<const crs_matrix_type> (A_local_crs_nc);
00489       }
00490       A_local_crs_ = A_local_crs;
00491       Graph_ = rcp (new Ifpack2::IlukGraph<crs_graph_type> (A_local_crs->getCrsGraph (),
00492                                                             LevelOfFill_, 0));
00493     }
00494 
00495     Graph_->initialize ();
00496     allocate_L_and_U ();
00497     initAllValues (*A_local_crs_);
00498   } // Stop timing
00499 
00500   isInitialized_ = true;
00501   ++numInitialize_;
00502   initializeTime_ += timer.totalElapsedTime ();
00503 }
00504 
00505 
00506 template<class MatrixType>
00507 void
00508 RILUK<MatrixType>::
00509 initAllValues (const row_matrix_type& A)
00510 {
00511   using Teuchos::ArrayRCP;
00512   using Teuchos::Comm;
00513   using Teuchos::ptr;
00514   using Teuchos::RCP;
00515   using Teuchos::REDUCE_SUM;
00516   using Teuchos::reduceAll;
00517   typedef Tpetra::Map<local_ordinal_type,global_ordinal_type,node_type> map_type;
00518 
00519   size_t NumIn = 0, NumL = 0, NumU = 0;
00520   bool DiagFound = false;
00521   size_t NumNonzeroDiags = 0;
00522   size_t MaxNumEntries = A.getGlobalMaxNumRowEntries();
00523 
00524   Teuchos::Array<global_ordinal_type> InI(MaxNumEntries); // Allocate temp space
00525   Teuchos::Array<global_ordinal_type> LI(MaxNumEntries);
00526   Teuchos::Array<global_ordinal_type> UI(MaxNumEntries);
00527   Teuchos::Array<scalar_type> InV(MaxNumEntries);
00528   Teuchos::Array<scalar_type> LV(MaxNumEntries);
00529   Teuchos::Array<scalar_type> UV(MaxNumEntries);
00530 
00531   // Check if values should be inserted or replaced
00532   const bool ReplaceValues = L_->isStaticGraph () || L_->isLocallyIndexed ();
00533 
00534   L_->resumeFill ();
00535   U_->resumeFill ();
00536   if (ReplaceValues) {
00537     L_->setAllToScalar (STS::zero ()); // Zero out L and U matrices
00538     U_->setAllToScalar (STS::zero ());
00539   }
00540 
00541   D_->putScalar (STS::zero ()); // Set diagonal values to zero
00542   ArrayRCP<scalar_type> DV = D_->get1dViewNonConst (); // Get view of diagonal
00543 
00544   RCP<const map_type> rowMap = L_->getRowMap ();
00545 
00546   // First we copy the user's matrix into L and U, regardless of fill level
00547 
00548   Teuchos::ArrayView<const global_ordinal_type> nodeGIDs = rowMap->getNodeElementList();
00549   for (typename Teuchos::ArrayView<const global_ordinal_type>::const_iterator avi = nodeGIDs.begin(); avi != nodeGIDs.end(); avi++)
00550   {
00551     global_ordinal_type global_row = *avi;
00552     local_ordinal_type local_row = rowMap->getLocalElement (global_row);
00553 
00554     A.getGlobalRowCopy (global_row, InI(), InV(), NumIn); // Get Values and Indices
00555 
00556     // Split into L and U (we don't assume that indices are ordered).
00557 
00558     NumL = 0;
00559     NumU = 0;
00560     DiagFound = false;
00561 
00562     for (size_t j = 0; j < NumIn; ++j) {
00563       const global_ordinal_type k = InI[j];
00564 
00565       if (k == global_row) {
00566         DiagFound = true;
00567         // Store perturbed diagonal in Tpetra::Vector D_
00568         DV[local_row] += Rthresh_ * InV[j] + IFPACK2_SGN(InV[j]) * Athresh_;
00569       }
00570       else if (k < 0) { // Out of range
00571         TEUCHOS_TEST_FOR_EXCEPTION(
00572           true, std::runtime_error, "Ifpack2::RILUK::initAllValues: current "
00573           "GID k = " << k << " < 0.  I'm not sure why this is an error; it is "
00574           "probably an artifact of the undocumented assumptions of the "
00575           "original implementation (likely copied and pasted from Ifpack).  "
00576           "Nevertheless, the code I found here insisted on this being an error "
00577           "state, so I will throw an exception here.");
00578       }
00579       else if (k < global_row) {
00580         LI[NumL] = k;
00581         LV[NumL] = InV[j];
00582         NumL++;
00583       }
00584       else if (k <= rowMap->getMaxGlobalIndex()) {
00585         UI[NumU] = k;
00586         UV[NumU] = InV[j];
00587         NumU++;
00588       }
00589 //      else {
00590 //        throw std::runtime_error("out of range in Ifpack2::RILUK::initAllValues");
00591 //      }
00592     }
00593 
00594     // Check in things for this row of L and U
00595 
00596     if (DiagFound) {
00597       ++NumNonzeroDiags;
00598     } else {
00599       DV[local_row] = Athresh_;
00600     }
00601 
00602     if (NumL) {
00603       if (ReplaceValues) {
00604         L_->replaceGlobalValues(global_row, LI(0, NumL), LV(0,NumL));
00605       } else {
00606         L_->insertGlobalValues(global_row, LI(0,NumL), LV(0,NumL));
00607       }
00608     }
00609 
00610     if (NumU) {
00611       if (ReplaceValues) {
00612         U_->replaceGlobalValues(global_row, UI(0,NumU), UV(0,NumU));
00613       } else {
00614         U_->insertGlobalValues(global_row, UI(0,NumU), UV(0,NumU));
00615       }
00616     }
00617   }
00618 
00619   // The domain of L and the range of U are exactly their own row maps
00620   // (there is no communication).  The domain of U and the range of L
00621   // must be the same as those of the original matrix, However if the
00622   // original matrix is a VbrMatrix, these two latter maps are
00623   // translation from a block map to a point map.
00624   L_->fillComplete (L_->getColMap (), A_local_crs_->getRangeMap ());
00625   U_->fillComplete (A_local_crs_->getDomainMap (), U_->getRowMap ());
00626 
00627   // At this point L and U have the values of A in the structure of L
00628   // and U, and diagonal vector D
00629 
00630   isInitialized_ = true;
00631 }
00632 
00633 
00634 template<class MatrixType>
00635 void RILUK<MatrixType>::compute ()
00636 {
00637   // initialize() checks this too, but it's easier for users if the
00638   // error shows them the name of the method that they actually
00639   // called, rather than the name of some internally called method.
00640   TEUCHOS_TEST_FOR_EXCEPTION(
00641     A_.is_null (), std::runtime_error, "Ifpack2::RILUK::compute: "
00642     "The matrix is null.  Please call setMatrix() with a nonnull input "
00643     "before calling this method.");
00644 
00645   if (! isInitialized ()) {
00646     initialize (); // Don't count this in the compute() time
00647   }
00648 
00649   Teuchos::Time timer ("RILUK::compute");
00650   { // Start timing
00651     isComputed_ = false;
00652 
00653     L_->resumeFill ();
00654     U_->resumeFill ();
00655 
00656     // MinMachNum should be officially defined, for now pick something a little
00657     // bigger than IEEE underflow value
00658 
00659     const scalar_type MinDiagonalValue = STS::rmin ();
00660     const scalar_type MaxDiagonalValue = STS::one () / MinDiagonalValue;
00661 
00662     size_t NumIn, NumL, NumU;
00663 
00664     // Get Maximum Row length
00665     const size_t MaxNumEntries =
00666       L_->getNodeMaxNumRowEntries () + U_->getNodeMaxNumRowEntries () + 1;
00667 
00668     Teuchos::Array<local_ordinal_type> InI(MaxNumEntries); // Allocate temp space
00669     Teuchos::Array<scalar_type> InV(MaxNumEntries);
00670     size_t num_cols = U_->getColMap()->getNodeNumElements();
00671     Teuchos::Array<int> colflag(num_cols);
00672 
00673     Teuchos::ArrayRCP<scalar_type> DV = D_->get1dViewNonConst(); // Get view of diagonal
00674 
00675     // Now start the factorization.
00676 
00677     // Need some integer workspace and pointers
00678     size_t NumUU;
00679     Teuchos::ArrayView<const local_ordinal_type> UUI;
00680     Teuchos::ArrayView<const scalar_type> UUV;
00681     for (size_t j = 0; j < num_cols; ++j) {
00682       colflag[j] = -1;
00683     }
00684 
00685     for (size_t i = 0; i < L_->getNodeNumRows (); ++i) {
00686       local_ordinal_type local_row = i;
00687 
00688       // Fill InV, InI with current row of L, D and U combined
00689 
00690       NumIn = MaxNumEntries;
00691       L_->getLocalRowCopy (local_row, InI (), InV (), NumL);
00692 
00693       InV[NumL] = DV[i]; // Put in diagonal
00694       InI[NumL] = local_row;
00695 
00696       U_->getLocalRowCopy (local_row, InI (NumL+1, MaxNumEntries-NumL-1),
00697                            InV (NumL+1, MaxNumEntries-NumL-1), NumU);
00698       NumIn = NumL+NumU+1;
00699 
00700       // Set column flags
00701       for (size_t j = 0; j < NumIn; ++j) {
00702         colflag[InI[j]] = j;
00703       }
00704 
00705       scalar_type diagmod = STS::zero (); // Off-diagonal accumulator
00706 
00707       for (size_t jj = 0; jj < NumL; ++jj) {
00708         local_ordinal_type j = InI[jj];
00709         scalar_type multiplier = InV[jj]; // current_mults++;
00710 
00711         InV[jj] *= DV[j];
00712 
00713         U_->getLocalRowView(j, UUI, UUV); // View of row above
00714         NumUU = UUI.size();
00715 
00716         if (RelaxValue_ == STM::zero ()) {
00717           for (size_t k = 0; k < NumUU; ++k) {
00718             const int kk = colflag[UUI[k]];
00719             // FIXME (mfh 23 Dec 2013) Wait a second, we just set
00720             // colflag above using size_t (which is generally unsigned),
00721             // but now we're querying it using int (which is signed).
00722             if (kk > -1) {
00723               InV[kk] -= multiplier * UUV[k];
00724             }
00725           }
00726         }
00727         else {
00728           for (size_t k = 0; k < NumUU; ++k) {
00729             // FIXME (mfh 23 Dec 2013) Wait a second, we just set
00730             // colflag above using size_t (which is generally unsigned),
00731             // but now we're querying it using int (which is signed).
00732             const int kk = colflag[UUI[k]];
00733             if (kk > -1) {
00734               InV[kk] -= multiplier*UUV[k];
00735             }
00736             else {
00737               diagmod -= multiplier*UUV[k];
00738             }
00739           }
00740         }
00741       }
00742       if (NumL) {
00743         // Replace current row of L
00744         L_->replaceLocalValues (local_row, InI (0, NumL), InV (0, NumL));
00745       }
00746 
00747       DV[i] = InV[NumL]; // Extract Diagonal value
00748 
00749       if (RelaxValue_ != STM::zero ()) {
00750         DV[i] += RelaxValue_*diagmod; // Add off diagonal modifications
00751       }
00752 
00753       if (STS::magnitude (DV[i]) > STS::magnitude (MaxDiagonalValue)) {
00754         if (STS::real (DV[i]) < STM::zero ()) {
00755           DV[i] = -MinDiagonalValue;
00756         }
00757         else {
00758           DV[i] = MinDiagonalValue;
00759         }
00760       }
00761       else {
00762         DV[i] = STS::one () / DV[i]; // Invert diagonal value
00763       }
00764 
00765       for (size_t j = 0; j < NumU; ++j) {
00766         InV[NumL+1+j] *= DV[i]; // Scale U by inverse of diagonal
00767       }
00768 
00769       if (NumU) {
00770         // Replace current row of L and U
00771         U_->replaceLocalValues (local_row, InI (NumL+1, NumU), InV (NumL+1, NumU));
00772       }
00773 
00774       // Reset column flags
00775       for (size_t j = 0; j < NumIn; ++j) {
00776         colflag[InI[j]] = -1;
00777       }
00778     }
00779 
00780     // FIXME (mfh 23 Dec 2013) Do we know that the column Map of L_ is
00781     // always one-to-one?
00782     L_->fillComplete (L_->getColMap (), A_local_crs_->getRangeMap ());
00783     U_->fillComplete (A_local_crs_->getDomainMap (), U_->getRowMap ());
00784 
00785     // Validate that the L and U factors are actually lower and upper triangular
00786     TEUCHOS_TEST_FOR_EXCEPTION(
00787       ! L_->isLowerTriangular (), std::runtime_error,
00788       "Ifpack2::RILUK::compute: L isn't lower triangular.");
00789     TEUCHOS_TEST_FOR_EXCEPTION(
00790       ! U_->isUpperTriangular (), std::runtime_error,
00791       "Ifpack2::RILUK::compute: U isn't lower triangular.");
00792   } // Stop timing
00793 
00794   isComputed_ = true;
00795   ++numCompute_;
00796   computeTime_ += timer.totalElapsedTime ();
00797 }
00798 
00799 
00800 template<class MatrixType>
00801 void
00802 RILUK<MatrixType>::
00803 apply (const Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& X,
00804        Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& Y,
00805        Teuchos::ETransp mode,
00806        scalar_type alpha,
00807        scalar_type beta) const
00808 {
00809   using Teuchos::RCP;
00810   using Teuchos::rcpFromRef;
00811 
00812   TEUCHOS_TEST_FOR_EXCEPTION(
00813     A_.is_null (), std::runtime_error, "Ifpack2::RILUK::apply: The matrix is "
00814     "null.  Please call setMatrix() with a nonnull input, then initialize() "
00815     "and compute(), before calling this method.");
00816   TEUCHOS_TEST_FOR_EXCEPTION(
00817     ! isComputed (), std::runtime_error,
00818     "Ifpack2::RILUK::apply: If you have not yet called compute(), "
00819     "you must call compute() before calling this method.");
00820   TEUCHOS_TEST_FOR_EXCEPTION(
00821     X.getNumVectors () != Y.getNumVectors (), std::invalid_argument,
00822     "Ifpack2::RILUK::apply: X and Y do not have the same number of columns.  "
00823     "X.getNumVectors() = " << X.getNumVectors ()
00824     << " != Y.getNumVectors() = " << Y.getNumVectors () << ".");
00825   TEUCHOS_TEST_FOR_EXCEPTION(
00826     STS::isComplex && mode == Teuchos::CONJ_TRANS, std::logic_error,
00827     "Ifpack2::RILUK::apply: mode = Teuchos::CONJ_TRANS is not implemented for "
00828     "complex Scalar type.  Please talk to the Ifpack2 developers to get this "
00829     "fixed.  There is a FIXME in this file about this very issue.");
00830 
00831   const scalar_type one = STS::one ();
00832   const scalar_type zero = STS::zero ();
00833 
00834   Teuchos::Time timer ("RILUK::apply");
00835   { // Start timing
00836     Teuchos::TimeMonitor timeMon (timer);
00837     if (alpha == one && beta == zero) {
00838       if (mode == Teuchos::NO_TRANS) { // Solve L (D (U Y)) = X for Y.
00839         // Start by solving L C = X for C.  C must have the same Map
00840         // as D.  We have to use a temp multivector, since
00841         // localSolve() does not allow its input and output to alias
00842         // one another.
00843         //
00844         // FIXME (mfh 24 Jan 2014) Cache this temp multivector.
00845         MV C (D_->getMap (), X.getNumVectors ());
00846         L_->localSolve (X, C, mode);
00847 
00848         // Solve D Y_tmp = C.  Y_tmp must have the same Map as C, and
00849         // the operation lets us do this in place in C, so we can
00850         // write "solve D C = C for C."
00851         C.elementWiseMultiply (one, *D_, C, zero);
00852 
00853         U_->localSolve (C, Y, mode); // Solve U Y = C.
00854       }
00855       else { // Solve U^P (D^P (U^P Y)) = X for Y (where P is * or T).
00856 
00857         // Start by solving U^P C = X for C.  C must have the same Map
00858         // as D.  We have to use a temp multivector, since
00859         // localSolve() does not allow its input and output to alias
00860         // one another.
00861         //
00862         // FIXME (mfh 24 Jan 2014) Cache this temp multivector.
00863         MV C (D_->getMap (), X.getNumVectors ());
00864         U_->localSolve (X, C, mode);
00865 
00866         // Solve D^P Y_tmp = C.  Y_tmp must have the same Map as C,
00867         // and the operation lets us do this in place in C, so we can
00868         // write "solve D^P C = C for C."
00869         //
00870         // FIXME (mfh 24 Jan 2014) If mode = Teuchos::CONJ_TRANS, we
00871         // need to do an elementwise multiply with the conjugate of
00872         // D_, not just with D_ itself.
00873         C.elementWiseMultiply (one, *D_, C, zero);
00874 
00875         L_->localSolve (C, Y, mode); // Solve L^P Y = C.
00876       }
00877     }
00878     else { // alpha != 1 or beta != 0
00879       if (alpha == zero) {
00880         if (beta == zero) {
00881           Y.putScalar (zero);
00882         } else {
00883           Y.scale (beta);
00884         }
00885       } else { // alpha != zero
00886         MV Y_tmp (Y.getMap (), Y.getNumVectors ());
00887         apply (X, Y_tmp, mode);
00888         Y.update (alpha, Y_tmp, beta);
00889       }
00890     }
00891   } // Stop timing
00892 
00893   ++numApply_;
00894   applyTime_ = timer.totalElapsedTime ();
00895 }
00896 
00897 
00898 template<class MatrixType>
00899 void RILUK<MatrixType>::
00900 multiply (const Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& X,
00901           Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& Y,
00902           const Teuchos::ETransp mode) const
00903 {
00904   const scalar_type zero = STS::zero ();
00905   const scalar_type one = STS::one ();
00906 
00907   if (mode != Teuchos::NO_TRANS) {
00908     U_->apply (X, Y, mode); //
00909     Y.update (one, X, one); // Y = Y + X (account for implicit unit diagonal)
00910 
00911     // FIXME (mfh 24 Jan 2014) If mode = Teuchos::CONJ_TRANS, we need
00912     // to do an elementwise multiply with the conjugate of D_, not
00913     // just with D_ itself.
00914     Y.elementWiseMultiply (one, *D_, Y, zero); // y = D*y (D_ has inverse of diagonal)
00915 
00916     MV Y_tmp (Y); // Need a temp copy of Y
00917     L_->apply (Y_tmp, Y, mode);
00918     Y.update (one, Y_tmp, one); // (account for implicit unit diagonal)
00919   }
00920   else {
00921     L_->apply (X, Y, mode);
00922     Y.update (one, X, one); // Y = Y + X (account for implicit unit diagonal)
00923     Y.elementWiseMultiply (one, *D_, Y, zero); // y = D*y (D_ has inverse of diagonal)
00924     MV Y_tmp (Y); // Need a temp copy of Y1
00925     U_->apply (Y_tmp, Y, mode);
00926     Y.update (one, Y_tmp, one); // (account for implicit unit diagonal)
00927   }
00928 }
00929 
00930 
00931 template<class MatrixType>
00932 typename Teuchos::ScalarTraits<typename MatrixType::scalar_type>::magnitudeType
00933 RILUK<MatrixType>::computeCondEst (Teuchos::ETransp mode) const
00934 {
00935   if (Condest_ != -Teuchos::ScalarTraits<magnitude_type>::one ()) {
00936     return Condest_;
00937   }
00938   // Create a vector with all values equal to one
00939   vec_type ones (U_->getDomainMap ());
00940   vec_type onesResult (L_->getRangeMap ());
00941   ones.putScalar (Teuchos::ScalarTraits<scalar_type>::one ());
00942 
00943   apply (ones, onesResult, mode); // Compute the effect of the solve on the vector of ones
00944   onesResult.abs (onesResult); // Make all values non-negative
00945   Teuchos::Array<magnitude_type> norms (1);
00946   onesResult.normInf (norms ());
00947   Condest_ = norms[0];
00948   return Condest_;
00949 }
00950 
00951 
00952 template<class MatrixType>
00953 typename Teuchos::ScalarTraits<typename MatrixType::scalar_type>::magnitudeType
00954 RILUK<MatrixType>::
00955 computeCondEst (CondestType CT,
00956                 local_ordinal_type MaxIters,
00957                 magnitude_type Tol,
00958                 const Teuchos::Ptr<const Tpetra::RowMatrix<scalar_type,local_ordinal_type,global_ordinal_type,node_type> >& Matrix)
00959 {
00960   // Forestall "unused variable" compiler warnings.
00961   (void) CT;
00962   (void) MaxIters;
00963   (void) Tol;
00964   (void) Matrix;
00965 
00966   if (Condest_ != -Teuchos::ScalarTraits<magnitude_type>::one() ) {
00967     return Condest_;
00968   }
00969   // Create a vector with all values equal to one
00970   vec_type ones (U_->getDomainMap ());
00971   vec_type onesResult (L_->getRangeMap ());
00972   ones.putScalar (Teuchos::ScalarTraits<scalar_type>::one ());
00973 
00974   // Compute the effect of the solve on the vector of ones
00975   apply (ones, onesResult, Teuchos::NO_TRANS);
00976   onesResult.abs (onesResult); // Make all values non-negative
00977   Teuchos::Array<magnitude_type> norms (1);
00978   onesResult.normInf (norms ());
00979   Condest_ = norms[0];
00980   return Condest_;
00981 }
00982 
00983 
00984 template<class MatrixType>
00985 std::string RILUK<MatrixType>::description () const
00986 {
00987   std::ostringstream os;
00988 
00989   // Output is a valid YAML dictionary in flow style.  If you don't
00990   // like everything on a single line, you should call describe()
00991   // instead.
00992   os << "\"Ifpack2::RILUK\": {";
00993   os << "Initialized: " << (isInitialized () ? "true" : "false") << ", "
00994      << "Computed: " << (isComputed () ? "true" : "false") << ", ";
00995 
00996   os << "Level-of-fill: " << getLevelOfFill() << ", ";
00997 
00998   if (A_.is_null ()) {
00999     os << "Matrix: null";
01000   }
01001   else {
01002     os << "Global matrix dimensions: ["
01003        << A_->getGlobalNumRows () << ", " << A_->getGlobalNumCols () << "]"
01004        << ", Global nnz: " << A_->getGlobalNumEntries();
01005   }
01006 
01007   os << "}";
01008   return os.str ();
01009 }
01010 
01011 
01012 } // namespace Ifpack2
01013 
01014 #define IFPACK2_RILUK_INSTANT(S,LO,GO,N)                            \
01015   template class Ifpack2::RILUK< Tpetra::CrsMatrix<S, LO, GO, N> >; \
01016   template class Ifpack2::RILUK< Tpetra::RowMatrix<S, LO, GO, N> >;
01017 
01018 #endif
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends