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   // First check that the local row map ordering is the same as the local portion of the column map.
00525   // The extraction of the strictly lower/upper parts of A, as well as the factorization,
00526   // implicitly assume that this is the case.
00527   Teuchos::ArrayView<const global_ordinal_type> rowGIDs = A.getRowMap()->getNodeElementList();
00528   Teuchos::ArrayView<const global_ordinal_type> colGIDs = A.getColMap()->getNodeElementList();
00529   bool gidsAreConsistentlyOrdered=true;
00530   global_ordinal_type indexOfInconsistentGID=0;
00531   for (global_ordinal_type i=0; i<rowGIDs.size(); ++i) {
00532     if (rowGIDs[i] != colGIDs[i]) {
00533       gidsAreConsistentlyOrdered=false;
00534       indexOfInconsistentGID=i;
00535       break;
00536     }
00537   }
00538   TEUCHOS_TEST_FOR_EXCEPTION(gidsAreConsistentlyOrdered==false, std::runtime_error, 
00539                              "The ordering of the local GIDs in the row and column maps is not the same"
00540                              << std::endl << "at index " << indexOfInconsistentGID
00541                              << ".  Consistency is required, as all calculations are done with"
00542                              << std::endl << "local indexing.");
00543 
00544   // Allocate temporary space for extracting the strictly
00545   // lower and upper parts of the matrix A.
00546   Teuchos::Array<local_ordinal_type> InI(MaxNumEntries);
00547   Teuchos::Array<local_ordinal_type> LI(MaxNumEntries);
00548   Teuchos::Array<local_ordinal_type> UI(MaxNumEntries);
00549   Teuchos::Array<scalar_type> InV(MaxNumEntries);
00550   Teuchos::Array<scalar_type> LV(MaxNumEntries);
00551   Teuchos::Array<scalar_type> UV(MaxNumEntries);
00552 
00553   // Check if values should be inserted or replaced
00554   const bool ReplaceValues = L_->isStaticGraph () || L_->isLocallyIndexed ();
00555 
00556   L_->resumeFill ();
00557   U_->resumeFill ();
00558   if (ReplaceValues) {
00559     L_->setAllToScalar (STS::zero ()); // Zero out L and U matrices
00560     U_->setAllToScalar (STS::zero ());
00561   }
00562 
00563   D_->putScalar (STS::zero ()); // Set diagonal values to zero
00564   ArrayRCP<scalar_type> DV = D_->get1dViewNonConst (); // Get view of diagonal
00565 
00566   RCP<const map_type> rowMap = L_->getRowMap ();
00567 
00568   // First we copy the user's matrix into L and U, regardless of fill level.
00569   // It is important to note that L and U are populated using local indices.
00570   // This means that if the row map GIDs are not monotonically increasing
00571   // (i.e., permuted or gappy), then the strictly lower (upper) part of the
00572   // matrix is not the one that you would get if you based L (U) on GIDs.
00573   // This is ok, as the *order* of the GIDs in the rowmap is a better
00574   // expression of the user's intent than the GIDs themselves.
00575 
00576   Teuchos::ArrayView<const global_ordinal_type> nodeGIDs = rowMap->getNodeElementList();
00577   for (size_t myRow=0; myRow<A.getNodeNumRows(); ++myRow) {
00578     local_ordinal_type local_row = myRow;
00579 
00580     //TODO JJH 4April2014 An optimization is to use getLocalRowView.  Not all matrices support this,
00581     //                    we'd need to check via the Tpetra::RowMatrix method supportsRowViews().
00582     A.getLocalRowCopy (local_row, InI(), InV(), NumIn); // Get Values and Indices
00583 
00584     // Split into L and U (we don't assume that indices are ordered).
00585 
00586     NumL = 0;
00587     NumU = 0;
00588     DiagFound = false;
00589 
00590     for (size_t j = 0; j < NumIn; ++j) {
00591       const local_ordinal_type k = InI[j];
00592 
00593       if (k == local_row) {
00594         DiagFound = true;
00595         // Store perturbed diagonal in Tpetra::Vector D_
00596         DV[local_row] += Rthresh_ * InV[j] + IFPACK2_SGN(InV[j]) * Athresh_;
00597       }
00598       else if (k < 0) { // Out of range
00599         TEUCHOS_TEST_FOR_EXCEPTION(
00600           true, std::runtime_error, "Ifpack2::RILUK::initAllValues: current "
00601           "GID k = " << k << " < 0.  I'm not sure why this is an error; it is "
00602           "probably an artifact of the undocumented assumptions of the "
00603           "original implementation (likely copied and pasted from Ifpack).  "
00604           "Nevertheless, the code I found here insisted on this being an error "
00605           "state, so I will throw an exception here.");
00606       }
00607       else if (k < local_row) {
00608         LI[NumL] = k;
00609         LV[NumL] = InV[j];
00610         NumL++;
00611       }
00612       else if (Teuchos::as<size_t>(k) <= rowMap->getNodeNumElements()) {
00613         UI[NumU] = k;
00614         UV[NumU] = InV[j];
00615         NumU++;
00616       }
00617     }
00618 
00619     // Check in things for this row of L and U
00620 
00621     if (DiagFound) {
00622       ++NumNonzeroDiags;
00623     } else {
00624       DV[local_row] = Athresh_;
00625     }
00626 
00627     if (NumL) {
00628       if (ReplaceValues) {
00629         L_->replaceLocalValues(local_row, LI(0, NumL), LV(0,NumL));
00630       } else {
00631         //FIXME JJH 24April2014 Is this correct?  I believe this case is when there aren't already values
00632         //FIXME in this row in the column locations corresponding to UI.
00633         L_->insertLocalValues(local_row, LI(0,NumL), LV(0,NumL));
00634       }
00635     }
00636 
00637     if (NumU) {
00638       if (ReplaceValues) {
00639         U_->replaceLocalValues(local_row, UI(0,NumU), UV(0,NumU));
00640       } else {
00641         //FIXME JJH 24April2014 Is this correct?  I believe this case is when there aren't already values
00642         //FIXME in this row in the column locations corresponding to UI.
00643         U_->insertLocalValues(local_row, UI(0,NumU), UV(0,NumU));
00644       }
00645     }
00646   }
00647 
00648   // The domain of L and the range of U are exactly their own row maps
00649   // (there is no communication).  The domain of U and the range of L
00650   // must be the same as those of the original matrix, However if the
00651   // original matrix is a VbrMatrix, these two latter maps are
00652   // translation from a block map to a point map.
00653   L_->fillComplete (L_->getColMap (), A_local_crs_->getRangeMap ());
00654   U_->fillComplete (A_local_crs_->getDomainMap (), U_->getRowMap ());
00655 
00656   // At this point L and U have the values of A in the structure of L
00657   // and U, and diagonal vector D
00658 
00659   isInitialized_ = true;
00660 }
00661 
00662 
00663 template<class MatrixType>
00664 void RILUK<MatrixType>::compute ()
00665 {
00666   // initialize() checks this too, but it's easier for users if the
00667   // error shows them the name of the method that they actually
00668   // called, rather than the name of some internally called method.
00669   TEUCHOS_TEST_FOR_EXCEPTION(
00670     A_.is_null (), std::runtime_error, "Ifpack2::RILUK::compute: "
00671     "The matrix is null.  Please call setMatrix() with a nonnull input "
00672     "before calling this method.");
00673 
00674   if (! isInitialized ()) {
00675     initialize (); // Don't count this in the compute() time
00676   }
00677 
00678   Teuchos::Time timer ("RILUK::compute");
00679   { // Start timing
00680     isComputed_ = false;
00681 
00682     L_->resumeFill ();
00683     U_->resumeFill ();
00684 
00685     // MinMachNum should be officially defined, for now pick something a little
00686     // bigger than IEEE underflow value
00687 
00688     const scalar_type MinDiagonalValue = STS::rmin ();
00689     const scalar_type MaxDiagonalValue = STS::one () / MinDiagonalValue;
00690 
00691     size_t NumIn, NumL, NumU;
00692 
00693     // Get Maximum Row length
00694     const size_t MaxNumEntries =
00695       L_->getNodeMaxNumRowEntries () + U_->getNodeMaxNumRowEntries () + 1;
00696 
00697     Teuchos::Array<local_ordinal_type> InI(MaxNumEntries); // Allocate temp space
00698     Teuchos::Array<scalar_type> InV(MaxNumEntries);
00699     size_t num_cols = U_->getColMap()->getNodeNumElements();
00700     Teuchos::Array<int> colflag(num_cols);
00701 
00702     Teuchos::ArrayRCP<scalar_type> DV = D_->get1dViewNonConst(); // Get view of diagonal
00703 
00704     // Now start the factorization.
00705 
00706     // Need some integer workspace and pointers
00707     size_t NumUU;
00708     Teuchos::ArrayView<const local_ordinal_type> UUI;
00709     Teuchos::ArrayView<const scalar_type> UUV;
00710     for (size_t j = 0; j < num_cols; ++j) {
00711       colflag[j] = -1;
00712     }
00713 
00714     for (size_t i = 0; i < L_->getNodeNumRows (); ++i) {
00715       local_ordinal_type local_row = i;
00716 
00717       // Fill InV, InI with current row of L, D and U combined
00718 
00719       NumIn = MaxNumEntries;
00720       L_->getLocalRowCopy (local_row, InI (), InV (), NumL);
00721 
00722       InV[NumL] = DV[i]; // Put in diagonal
00723       InI[NumL] = local_row;
00724 
00725       U_->getLocalRowCopy (local_row, InI (NumL+1, MaxNumEntries-NumL-1),
00726                            InV (NumL+1, MaxNumEntries-NumL-1), NumU);
00727       NumIn = NumL+NumU+1;
00728 
00729       // Set column flags
00730       for (size_t j = 0; j < NumIn; ++j) {
00731         colflag[InI[j]] = j;
00732       }
00733 
00734       scalar_type diagmod = STS::zero (); // Off-diagonal accumulator
00735 
00736       for (size_t jj = 0; jj < NumL; ++jj) {
00737         local_ordinal_type j = InI[jj];
00738         scalar_type multiplier = InV[jj]; // current_mults++;
00739 
00740         InV[jj] *= DV[j];
00741 
00742         U_->getLocalRowView(j, UUI, UUV); // View of row above
00743         NumUU = UUI.size();
00744 
00745         if (RelaxValue_ == STM::zero ()) {
00746           for (size_t k = 0; k < NumUU; ++k) {
00747             const int kk = colflag[UUI[k]];
00748             // FIXME (mfh 23 Dec 2013) Wait a second, we just set
00749             // colflag above using size_t (which is generally unsigned),
00750             // but now we're querying it using int (which is signed).
00751             if (kk > -1) {
00752               InV[kk] -= multiplier * UUV[k];
00753             }
00754           }
00755         }
00756         else {
00757           for (size_t k = 0; k < NumUU; ++k) {
00758             // FIXME (mfh 23 Dec 2013) Wait a second, we just set
00759             // colflag above using size_t (which is generally unsigned),
00760             // but now we're querying it using int (which is signed).
00761             const int kk = colflag[UUI[k]];
00762             if (kk > -1) {
00763               InV[kk] -= multiplier*UUV[k];
00764             }
00765             else {
00766               diagmod -= multiplier*UUV[k];
00767             }
00768           }
00769         }
00770       }
00771       if (NumL) {
00772         // Replace current row of L
00773         L_->replaceLocalValues (local_row, InI (0, NumL), InV (0, NumL));
00774       }
00775 
00776       DV[i] = InV[NumL]; // Extract Diagonal value
00777 
00778       if (RelaxValue_ != STM::zero ()) {
00779         DV[i] += RelaxValue_*diagmod; // Add off diagonal modifications
00780       }
00781 
00782       if (STS::magnitude (DV[i]) > STS::magnitude (MaxDiagonalValue)) {
00783         if (STS::real (DV[i]) < STM::zero ()) {
00784           DV[i] = -MinDiagonalValue;
00785         }
00786         else {
00787           DV[i] = MinDiagonalValue;
00788         }
00789       }
00790       else {
00791         DV[i] = STS::one () / DV[i]; // Invert diagonal value
00792       }
00793 
00794       for (size_t j = 0; j < NumU; ++j) {
00795         InV[NumL+1+j] *= DV[i]; // Scale U by inverse of diagonal
00796       }
00797 
00798       if (NumU) {
00799         // Replace current row of L and U
00800         U_->replaceLocalValues (local_row, InI (NumL+1, NumU), InV (NumL+1, NumU));
00801       }
00802 
00803       // Reset column flags
00804       for (size_t j = 0; j < NumIn; ++j) {
00805         colflag[InI[j]] = -1;
00806       }
00807     }
00808 
00809     // FIXME (mfh 23 Dec 2013) Do we know that the column Map of L_ is
00810     // always one-to-one?
00811     L_->fillComplete (L_->getColMap (), A_local_crs_->getRangeMap ());
00812     U_->fillComplete (A_local_crs_->getDomainMap (), U_->getRowMap ());
00813 
00814     // Validate that the L and U factors are actually lower and upper triangular
00815     TEUCHOS_TEST_FOR_EXCEPTION(
00816       ! L_->isLowerTriangular (), std::runtime_error,
00817       "Ifpack2::RILUK::compute: L isn't lower triangular.");
00818     TEUCHOS_TEST_FOR_EXCEPTION(
00819       ! U_->isUpperTriangular (), std::runtime_error,
00820       "Ifpack2::RILUK::compute: U isn't lower triangular.");
00821   } // Stop timing
00822 
00823   isComputed_ = true;
00824   ++numCompute_;
00825   computeTime_ += timer.totalElapsedTime ();
00826 }
00827 
00828 
00829 template<class MatrixType>
00830 void
00831 RILUK<MatrixType>::
00832 apply (const Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& X,
00833        Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& Y,
00834        Teuchos::ETransp mode,
00835        scalar_type alpha,
00836        scalar_type beta) const
00837 {
00838   using Teuchos::RCP;
00839   using Teuchos::rcpFromRef;
00840 
00841   TEUCHOS_TEST_FOR_EXCEPTION(
00842     A_.is_null (), std::runtime_error, "Ifpack2::RILUK::apply: The matrix is "
00843     "null.  Please call setMatrix() with a nonnull input, then initialize() "
00844     "and compute(), before calling this method.");
00845   TEUCHOS_TEST_FOR_EXCEPTION(
00846     ! isComputed (), std::runtime_error,
00847     "Ifpack2::RILUK::apply: If you have not yet called compute(), "
00848     "you must call compute() before calling this method.");
00849   TEUCHOS_TEST_FOR_EXCEPTION(
00850     X.getNumVectors () != Y.getNumVectors (), std::invalid_argument,
00851     "Ifpack2::RILUK::apply: X and Y do not have the same number of columns.  "
00852     "X.getNumVectors() = " << X.getNumVectors ()
00853     << " != Y.getNumVectors() = " << Y.getNumVectors () << ".");
00854   TEUCHOS_TEST_FOR_EXCEPTION(
00855     STS::isComplex && mode == Teuchos::CONJ_TRANS, std::logic_error,
00856     "Ifpack2::RILUK::apply: mode = Teuchos::CONJ_TRANS is not implemented for "
00857     "complex Scalar type.  Please talk to the Ifpack2 developers to get this "
00858     "fixed.  There is a FIXME in this file about this very issue.");
00859 
00860   const scalar_type one = STS::one ();
00861   const scalar_type zero = STS::zero ();
00862 
00863   Teuchos::Time timer ("RILUK::apply");
00864   { // Start timing
00865     Teuchos::TimeMonitor timeMon (timer);
00866     if (alpha == one && beta == zero) {
00867       if (mode == Teuchos::NO_TRANS) { // Solve L (D (U Y)) = X for Y.
00868         // Start by solving L C = X for C.  C must have the same Map
00869         // as D.  We have to use a temp multivector, since
00870         // localSolve() does not allow its input and output to alias
00871         // one another.
00872         //
00873         // FIXME (mfh 24 Jan 2014) Cache this temp multivector.
00874         MV C (D_->getMap (), X.getNumVectors ());
00875         L_->localSolve (X, C, mode);
00876 
00877         // Solve D Y_tmp = C.  Y_tmp must have the same Map as C, and
00878         // the operation lets us do this in place in C, so we can
00879         // write "solve D C = C for C."
00880         C.elementWiseMultiply (one, *D_, C, zero);
00881 
00882         U_->localSolve (C, Y, mode); // Solve U Y = C.
00883       }
00884       else { // Solve U^P (D^P (U^P Y)) = X for Y (where P is * or T).
00885 
00886         // Start by solving U^P C = X for C.  C must have the same Map
00887         // as D.  We have to use a temp multivector, since
00888         // localSolve() does not allow its input and output to alias
00889         // one another.
00890         //
00891         // FIXME (mfh 24 Jan 2014) Cache this temp multivector.
00892         MV C (D_->getMap (), X.getNumVectors ());
00893         U_->localSolve (X, C, mode);
00894 
00895         // Solve D^P Y_tmp = C.  Y_tmp must have the same Map as C,
00896         // and the operation lets us do this in place in C, so we can
00897         // write "solve D^P C = C for C."
00898         //
00899         // FIXME (mfh 24 Jan 2014) If mode = Teuchos::CONJ_TRANS, we
00900         // need to do an elementwise multiply with the conjugate of
00901         // D_, not just with D_ itself.
00902         C.elementWiseMultiply (one, *D_, C, zero);
00903 
00904         L_->localSolve (C, Y, mode); // Solve L^P Y = C.
00905       }
00906     }
00907     else { // alpha != 1 or beta != 0
00908       if (alpha == zero) {
00909         if (beta == zero) {
00910           Y.putScalar (zero);
00911         } else {
00912           Y.scale (beta);
00913         }
00914       } else { // alpha != zero
00915         MV Y_tmp (Y.getMap (), Y.getNumVectors ());
00916         apply (X, Y_tmp, mode);
00917         Y.update (alpha, Y_tmp, beta);
00918       }
00919     }
00920   } // Stop timing
00921 
00922   ++numApply_;
00923   applyTime_ = timer.totalElapsedTime ();
00924 }
00925 
00926 
00927 template<class MatrixType>
00928 void RILUK<MatrixType>::
00929 multiply (const Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& X,
00930           Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& Y,
00931           const Teuchos::ETransp mode) const
00932 {
00933   const scalar_type zero = STS::zero ();
00934   const scalar_type one = STS::one ();
00935 
00936   if (mode != Teuchos::NO_TRANS) {
00937     U_->apply (X, Y, mode); //
00938     Y.update (one, X, one); // Y = Y + X (account for implicit unit diagonal)
00939 
00940     // FIXME (mfh 24 Jan 2014) If mode = Teuchos::CONJ_TRANS, we need
00941     // to do an elementwise multiply with the conjugate of D_, not
00942     // just with D_ itself.
00943     Y.elementWiseMultiply (one, *D_, Y, zero); // y = D*y (D_ has inverse of diagonal)
00944 
00945     MV Y_tmp (Y); // Need a temp copy of Y
00946     L_->apply (Y_tmp, Y, mode);
00947     Y.update (one, Y_tmp, one); // (account for implicit unit diagonal)
00948   }
00949   else {
00950     L_->apply (X, Y, mode);
00951     Y.update (one, X, one); // Y = Y + X (account for implicit unit diagonal)
00952     Y.elementWiseMultiply (one, *D_, Y, zero); // y = D*y (D_ has inverse of diagonal)
00953     MV Y_tmp (Y); // Need a temp copy of Y1
00954     U_->apply (Y_tmp, Y, mode);
00955     Y.update (one, Y_tmp, one); // (account for implicit unit diagonal)
00956   }
00957 }
00958 
00959 
00960 template<class MatrixType>
00961 typename Teuchos::ScalarTraits<typename MatrixType::scalar_type>::magnitudeType
00962 RILUK<MatrixType>::computeCondEst (Teuchos::ETransp mode) const
00963 {
00964   if (Condest_ != -Teuchos::ScalarTraits<magnitude_type>::one ()) {
00965     return Condest_;
00966   }
00967   // Create a vector with all values equal to one
00968   vec_type ones (U_->getDomainMap ());
00969   vec_type onesResult (L_->getRangeMap ());
00970   ones.putScalar (Teuchos::ScalarTraits<scalar_type>::one ());
00971 
00972   apply (ones, onesResult, mode); // Compute the effect of the solve on the vector of ones
00973   onesResult.abs (onesResult); // Make all values non-negative
00974   Teuchos::Array<magnitude_type> norms (1);
00975   onesResult.normInf (norms ());
00976   Condest_ = norms[0];
00977   return Condest_;
00978 }
00979 
00980 
00981 template<class MatrixType>
00982 typename Teuchos::ScalarTraits<typename MatrixType::scalar_type>::magnitudeType
00983 RILUK<MatrixType>::
00984 computeCondEst (CondestType CT,
00985                 local_ordinal_type MaxIters,
00986                 magnitude_type Tol,
00987                 const Teuchos::Ptr<const Tpetra::RowMatrix<scalar_type,local_ordinal_type,global_ordinal_type,node_type> >& Matrix)
00988 {
00989   // Forestall "unused variable" compiler warnings.
00990   (void) CT;
00991   (void) MaxIters;
00992   (void) Tol;
00993   (void) Matrix;
00994 
00995   if (Condest_ != -Teuchos::ScalarTraits<magnitude_type>::one() ) {
00996     return Condest_;
00997   }
00998   // Create a vector with all values equal to one
00999   vec_type ones (U_->getDomainMap ());
01000   vec_type onesResult (L_->getRangeMap ());
01001   ones.putScalar (Teuchos::ScalarTraits<scalar_type>::one ());
01002 
01003   // Compute the effect of the solve on the vector of ones
01004   apply (ones, onesResult, Teuchos::NO_TRANS);
01005   onesResult.abs (onesResult); // Make all values non-negative
01006   Teuchos::Array<magnitude_type> norms (1);
01007   onesResult.normInf (norms ());
01008   Condest_ = norms[0];
01009   return Condest_;
01010 }
01011 
01012 
01013 template<class MatrixType>
01014 std::string RILUK<MatrixType>::description () const
01015 {
01016   std::ostringstream os;
01017 
01018   // Output is a valid YAML dictionary in flow style.  If you don't
01019   // like everything on a single line, you should call describe()
01020   // instead.
01021   os << "\"Ifpack2::RILUK\": {";
01022   os << "Initialized: " << (isInitialized () ? "true" : "false") << ", "
01023      << "Computed: " << (isComputed () ? "true" : "false") << ", ";
01024 
01025   os << "Level-of-fill: " << getLevelOfFill() << ", ";
01026 
01027   if (A_.is_null ()) {
01028     os << "Matrix: null";
01029   }
01030   else {
01031     os << "Global matrix dimensions: ["
01032        << A_->getGlobalNumRows () << ", " << A_->getGlobalNumCols () << "]"
01033        << ", Global nnz: " << A_->getGlobalNumEntries();
01034   }
01035 
01036   os << "}";
01037   return os.str ();
01038 }
01039 
01040 
01041 } // namespace Ifpack2
01042 
01043 #define IFPACK2_RILUK_INSTANT(S,LO,GO,N)                            \
01044   template class Ifpack2::RILUK< Tpetra::CrsMatrix<S, LO, GO, N> >; \
01045   template class Ifpack2::RILUK< Tpetra::RowMatrix<S, LO, GO, N> >;
01046 
01047 #endif
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends