Ifpack2 Templated Preconditioning Package Version 1.0
Ifpack2_SparseContainer_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_SPARSECONTAINER_DEF_HPP
00044 #define IFPACK2_SPARSECONTAINER_DEF_HPP
00045 
00046 #include "Ifpack2_SparseContainer_decl.hpp"
00047 #ifdef HAVE_MPI
00048 #include <mpi.h>
00049 #include "Teuchos_DefaultMpiComm.hpp"
00050 #else
00051 #include "Teuchos_DefaultSerialComm.hpp"
00052 #endif
00053 
00054 namespace Ifpack2 {
00055 
00056 //==============================================================================
00057 template<class MatrixType, class InverseType>
00058 SparseContainer<MatrixType,InverseType>::
00059 SparseContainer (const Teuchos::RCP<const row_matrix_type>& matrix,
00060                  const Teuchos::ArrayView<const local_ordinal_type>& localRows) :
00061   Container<MatrixType> (matrix, localRows),
00062   numRows_ (localRows.size ()),
00063   IsInitialized_ (false),
00064   IsComputed_ (false),
00065 #ifdef HAVE_MPI
00066   localComm_ (Teuchos::rcp (new Teuchos::MpiComm<int> (MPI_COMM_SELF)))
00067 #else
00068   localComm_ (Teuchos::rcp (new Teuchos::SerialComm<int> ()))
00069 #endif // HAVE_MPI
00070 {}
00071 
00072 //==============================================================================
00073 template<class MatrixType, class InverseType>
00074 SparseContainer<MatrixType,InverseType>::~SparseContainer()
00075 {}
00076 
00077 //==============================================================================
00078 template<class MatrixType, class InverseType>
00079 size_t SparseContainer<MatrixType,InverseType>::getNumRows() const
00080 {
00081   if (isInitialized ()) {
00082     return numRows_;
00083   } else {
00084     return 0;
00085   }
00086 }
00087 
00088 //==============================================================================
00089 template<class MatrixType, class InverseType>
00090 bool SparseContainer<MatrixType,InverseType>::isInitialized() const
00091 {
00092   return IsInitialized_;
00093 }
00094 
00095 //==============================================================================
00096 template<class MatrixType, class InverseType>
00097 bool SparseContainer<MatrixType,InverseType>::isComputed() const
00098 {
00099   return IsComputed_;
00100 }
00101 
00102 //==============================================================================
00103 template<class MatrixType, class InverseType>
00104 void SparseContainer<MatrixType,InverseType>::setParameters(const Teuchos::ParameterList& List)
00105 {
00106   List_ = List;
00107 }
00108 
00109 //==============================================================================
00110 template<class MatrixType, class InverseType>
00111 void SparseContainer<MatrixType,InverseType>::initialize ()
00112 {
00113   using Teuchos::null;
00114   using Teuchos::rcp;
00115   typedef Tpetra::Map<InverseLocalOrdinal, InverseGlobalOrdinal,
00116                       InverseNode> map_type;
00117   typedef Tpetra::CrsMatrix<InverseScalar, InverseLocalOrdinal,
00118                             InverseGlobalOrdinal, InverseNode> crs_matrix_type;
00119   // We assume that if you called this method, you intend to recompute
00120   // everything.  Thus, we release references to all the internal
00121   // objects.  We do this first to save memory.  (In an RCP
00122   // assignment, the right-hand side and left-hand side coexist before
00123   // the left-hand side's reference count gets updated.)
00124   IsInitialized_ = false;
00125   IsComputed_ = false;
00126   localMap_ = null;
00127   diagBlock_ = null;
00128 
00129   // (Re)create the local Map, and the CrsMatrix that will contain the
00130   // local matrix to use for solves.
00131   localMap_ = rcp (new map_type (numRows_, 0, localComm_));
00132   diagBlock_ = rcp (new crs_matrix_type (localMap_, 0));
00133 
00134   // Create the inverse operator using the local matrix.  We give it
00135   // the matrix here, but don't call its initialize() or compute()
00136   // methods yet, since we haven't initialized the matrix yet.
00137   Inverse_ = rcp (new InverseType (diagBlock_));
00138   Inverse_->setParameters (List_);
00139 
00140   IsInitialized_ = true;
00141 }
00142 
00143 //==============================================================================
00144 template<class MatrixType, class InverseType>
00145 void SparseContainer<MatrixType,InverseType>::compute ()
00146 {
00147   IsComputed_ = false;
00148   if (! this->isInitialized ()) {
00149     this->initialize ();
00150   }
00151 
00152   // Extract the submatrix.
00153   this->extract (this->getMatrix ());
00154 
00155   // The inverse operator already has a pointer to the submatrix.  Now
00156   // that the submatrix is filled in, we can initialize and compute
00157   // the inverse operator.
00158   Inverse_->initialize ();
00159   Inverse_->compute ();
00160 
00161   IsComputed_ = true;
00162 }
00163 
00164 //==============================================================================
00165 template<class MatrixType, class InverseType>
00166 void SparseContainer<MatrixType,InverseType>::
00167 applyImpl (const Tpetra::MultiVector<InverseScalar,InverseLocalOrdinal,InverseGlobalOrdinal,InverseNode>& X,
00168            Tpetra::MultiVector<InverseScalar,InverseLocalOrdinal,InverseGlobalOrdinal,InverseNode>& Y,
00169            Teuchos::ETransp mode,
00170            InverseScalar alpha,
00171            InverseScalar beta) const
00172 {
00173   TEUCHOS_TEST_FOR_EXCEPTION(
00174     Inverse_->getDomainMap ()->getNodeNumElements () != X.getLocalLength (),
00175     std::logic_error, "Ifpack2::SparseContainer::apply: Inverse_ "
00176     "operator and X have incompatible dimensions (" <<
00177     Inverse_->getDomainMap ()->getNodeNumElements () << " resp. "
00178     << X.getLocalLength () << ").  Please report this bug to "
00179     "the Ifpack2 developers.");
00180   TEUCHOS_TEST_FOR_EXCEPTION(
00181     Inverse_->getRangeMap ()->getNodeNumElements () != Y.getLocalLength (),
00182     std::logic_error, "Ifpack2::SparseContainer::apply: Inverse_ "
00183     "operator and Y have incompatible dimensions (" <<
00184     Inverse_->getRangeMap ()->getNodeNumElements () << " resp. "
00185     << Y.getLocalLength () << ").  Please report this bug to "
00186     "the Ifpack2 developers.");
00187 
00188   Inverse_->apply (X, Y, mode, alpha, beta);
00189 }
00190 
00191 //==============================================================================
00192 template<class MatrixType, class InverseType>
00193 void SparseContainer<MatrixType,InverseType>::
00194 apply (const Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& X,
00195        Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& Y,
00196        Teuchos::ETransp mode,
00197        scalar_type alpha,
00198        scalar_type beta) const
00199 {
00200   using Teuchos::ArrayView;
00201   using Teuchos::as;
00202   using Teuchos::RCP;
00203   using Teuchos::rcp;
00204 
00205   // The InverseType template parameter might have different template
00206   // parameters (Scalar, LO, GO, and/or Node) than MatrixType.  For
00207   // example, MatrixType (a global object) might use a bigger GO
00208   // (global ordinal) type than InverseType (which deals with the
00209   // diagonal block, a local object).  This means that we might have
00210   // to convert X and Y to the Tpetra::MultiVector specialization that
00211   // InverseType wants.  This class' X_ and Y_ internal fields are of
00212   // the right type for InverseType, so we can use those as targets.
00213 
00214   // Tpetra::MultiVector specialization corresponding to MatrixType.
00215   typedef Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type> MV_mat;
00216   // Tpetra::MultiVector specialization corresponding to InverseType.
00217   typedef Tpetra::MultiVector<InverseScalar,InverseLocalOrdinal,InverseGlobalOrdinal,InverseNode> MV_inv;
00218   Details::MultiVectorLocalGatherScatter<MV_mat, MV_inv> mvgs;
00219   const size_t numVecs = X.getNumVectors ();
00220 
00221   TEUCHOS_TEST_FOR_EXCEPTION(
00222     ! IsComputed_, std::runtime_error, "Ifpack2::SparseContainer::apply: "
00223     "You must have called the compute() method before you may call apply().  "
00224     "You may call the apply() method as many times as you want after calling "
00225     "compute() once, but you must have called compute() at least once.");
00226   TEUCHOS_TEST_FOR_EXCEPTION(
00227     numVecs != Y.getNumVectors (), std::runtime_error,
00228     "Ifpack2::SparseContainer::apply: X and Y have different numbers of "
00229     "vectors.  X has " << X.getNumVectors ()
00230     << ", but Y has " << X.getNumVectors () << ".");
00231 
00232   if (numVecs == 0) {
00233     return; // done! nothing to do
00234   }
00235 
00236   // The operator Inverse_ works on a permuted subset of the local
00237   // parts of X and Y.  The subset and permutation are defined by the
00238   // index array returned by getLocalRows().  If the permutation is
00239   // trivial and the subset is exactly equal to the local indices,
00240   // then we could use the local parts of X and Y exactly, without
00241   // needing to permute.  Otherwise, we have to use temporary storage
00242   // to permute X and Y.  For now, we always use temporary storage.
00243   //
00244   // Create temporary permuted versions of the input and output.
00245   // (Re)allocate X_ and/or Y_ only if necessary.  We'll use them to
00246   // store the permuted versions of X resp. Y.  Note that X_local has
00247   // the domain Map of the operator, which may be a permuted subset of
00248   // the local Map corresponding to X.getMap().  Similarly, Y_local
00249   // has the range Map of the operator, which may be a permuted subset
00250   // of the local Map corresponding to Y.getMap().  numRows_ here
00251   // gives the number of rows in the row Map of the local Inverse_
00252   // operator.
00253   //
00254   // FIXME (mfh 20 Aug 2013) There might be an implicit assumption
00255   // here that the row Map and the range Map of that operator are
00256   // the same.
00257   //
00258   // FIXME (mfh 20 Aug 2013) This "local permutation" functionality
00259   // really belongs in Tpetra.
00260 
00261   if (X_.is_null ()) {
00262     X_ = rcp (new MV_inv (Inverse_->getDomainMap (), numVecs));
00263   }
00264   RCP<MV_inv> X_local = X_;
00265   TEUCHOS_TEST_FOR_EXCEPTION(
00266     X_local->getLocalLength () != numRows_, std::logic_error,
00267     "Ifpack2::SparseContainer::apply: "
00268     "X_local has length " << X_local->getLocalLength () << ", which does "
00269     "not match numRows_ = " << numRows_ << ".  Please report this bug to "
00270     "the Ifpack2 developers.");
00271   ArrayView<const local_ordinal_type> localRows = this->getLocalRows ();
00272   mvgs.gather (*X_local, X, localRows);
00273 
00274   // We must gather the output multivector Y even on input to
00275   // Inverse_->apply(), since the Inverse_ operator might use it as an
00276   // initial guess for a linear solve.  We have no way of knowing
00277   // whether it does or does not.
00278 
00279   if (Y_.is_null ()) {
00280     Y_ = rcp (new MV_inv (Inverse_->getRangeMap (), numVecs));
00281   }
00282   RCP<MV_inv> Y_local = Y_;
00283   TEUCHOS_TEST_FOR_EXCEPTION(
00284     Y_local->getLocalLength () != numRows_, std::logic_error,
00285     "Ifpack2::SparseContainer::apply: "
00286     "Y_local has length " << X_local->getLocalLength () << ", which does "
00287     "not match numRows_ = " << numRows_ << ".  Please report this bug to "
00288     "the Ifpack2 developers.");
00289   mvgs.gather (*Y_local, Y, localRows);
00290 
00291   // Apply the local operator:
00292   // Y_local := beta*Y_local + alpha*M^{-1}*X_local
00293   this->applyImpl (*X_local, *Y_local, mode, as<InverseScalar> (alpha),
00294                    as<InverseScalar> (beta));
00295 
00296   // Scatter the permuted subset output vector Y_local back into the
00297   // original output multivector Y.
00298   mvgs.scatter (Y, *Y_local, localRows);
00299 }
00300 
00301 //==============================================================================
00302 template<class MatrixType, class InverseType>
00303 void SparseContainer<MatrixType,InverseType>::
00304 weightedApply (const Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& X,
00305                Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& Y,
00306                const Tpetra::Vector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& D,
00307                Teuchos::ETransp mode,
00308                scalar_type alpha,
00309                scalar_type beta) const
00310 {
00311   using Teuchos::ArrayRCP;
00312   using Teuchos::ArrayView;
00313   using Teuchos::Range1D;
00314   using Teuchos::RCP;
00315   using Teuchos::rcp;
00316   using Teuchos::rcp_const_cast;
00317   using std::cerr;
00318   using std::endl;
00319   typedef Teuchos::ScalarTraits<scalar_type> STS;
00320 
00321   // The InverseType template parameter might have different template
00322   // parameters (Scalar, LO, GO, and/or Node) than MatrixType.  For
00323   // example, MatrixType (a global object) might use a bigger GO
00324   // (global ordinal) type than InverseType (which deals with the
00325   // diagonal block, a local object).  This means that we might have
00326   // to convert X and Y to the Tpetra::MultiVector specialization that
00327   // InverseType wants.  This class' X_ and Y_ internal fields are of
00328   // the right type for InverseType, so we can use those as targets.
00329 
00330   // Tpetra::MultiVector specialization corresponding to MatrixType.
00331   typedef Tpetra::MultiVector<scalar_type, local_ordinal_type,
00332                               global_ordinal_type, node_type> MV_mat;
00333   // Tpetra::MultiVector specialization corresponding to InverseType.
00334   typedef Tpetra::MultiVector<InverseScalar, InverseLocalOrdinal,
00335                               InverseGlobalOrdinal, InverseNode> MV_inv;
00336   // Tpetra::Vector specialization corresponding to InverseType.
00337   typedef Tpetra::Vector<InverseScalar, InverseLocalOrdinal,
00338                          InverseGlobalOrdinal, InverseNode> V_inv;
00339   Details::MultiVectorLocalGatherScatter<MV_mat, MV_inv> mvgs;
00340   const size_t numVecs = X.getNumVectors ();
00341 
00342   TEUCHOS_TEST_FOR_EXCEPTION(
00343     ! IsComputed_, std::runtime_error, "Ifpack2::SparseContainer::"
00344     "weightedApply: You must have called the compute() method before you may "
00345     "call apply().  You may call the apply() method as many times as you want "
00346     "after calling compute() once, but you must have called compute() at least "
00347     "once.");
00348   TEUCHOS_TEST_FOR_EXCEPTION(
00349     Inverse_.is_null (), std::logic_error, "Ifpack2::SparseContainer::"
00350     "weightedApply: Inverse_ is null.  Please report this bug to the Ifpack2 "
00351     "developers.");
00352   TEUCHOS_TEST_FOR_EXCEPTION(
00353     numVecs != Y.getNumVectors (), std::runtime_error,
00354     "Ifpack2::SparseContainer::weightedApply: X and Y have different numbers "
00355     "of vectors.  X has " << X.getNumVectors () << ", but Y has "
00356     << X.getNumVectors () << ".");
00357 
00358   if (numVecs == 0) {
00359     return; // done! nothing to do
00360   }
00361 
00362   // The operator Inverse_ works on a permuted subset of the local
00363   // parts of X and Y.  The subset and permutation are defined by the
00364   // index array returned by getLocalRows().  If the permutation is
00365   // trivial and the subset is exactly equal to the local indices,
00366   // then we could use the local parts of X and Y exactly, without
00367   // needing to permute.  Otherwise, we have to use temporary storage
00368   // to permute X and Y.  For now, we always use temporary storage.
00369   //
00370   // Create temporary permuted versions of the input and output.
00371   // (Re)allocate X_ and/or Y_ only if necessary.  We'll use them to
00372   // store the permuted versions of X resp. Y.  Note that X_local has
00373   // the domain Map of the operator, which may be a permuted subset of
00374   // the local Map corresponding to X.getMap().  Similarly, Y_local
00375   // has the range Map of the operator, which may be a permuted subset
00376   // of the local Map corresponding to Y.getMap().  numRows_ here
00377   // gives the number of rows in the row Map of the local Inverse_
00378   // operator.
00379   //
00380   // FIXME (mfh 20 Aug 2013) There might be an implicit assumption
00381   // here that the row Map and the range Map of that operator are
00382   // the same.
00383   //
00384   // FIXME (mfh 20 Aug 2013) This "local permutation" functionality
00385   // really belongs in Tpetra.
00386 
00387   if (X_.is_null ()) {
00388     X_ = rcp (new MV_inv (Inverse_->getDomainMap (), numVecs));
00389   }
00390   RCP<MV_inv> X_local = X_;
00391   TEUCHOS_TEST_FOR_EXCEPTION(
00392     X_local->getLocalLength () != numRows_, std::logic_error,
00393     "Ifpack2::SparseContainer::weightedApply: "
00394     "X_local has length " << X_local->getLocalLength () << ", which does "
00395     "not match numRows_ = " << numRows_ << ".  Please report this bug to "
00396     "the Ifpack2 developers.");
00397   ArrayView<const local_ordinal_type> localRows = this->getLocalRows ();
00398   mvgs.gather (*X_local, X, localRows);
00399 
00400   // We must gather the output multivector Y even on input to
00401   // Inverse_->apply(), since the Inverse_ operator might use it as an
00402   // initial guess for a linear solve.  We have no way of knowing
00403   // whether it does or does not.
00404 
00405   if (Y_.is_null ()) {
00406     Y_ = rcp (new MV_inv (Inverse_->getRangeMap (), numVecs));
00407   }
00408   RCP<MV_inv> Y_local = Y_;
00409   TEUCHOS_TEST_FOR_EXCEPTION(
00410     Y_local->getLocalLength () != numRows_, std::logic_error,
00411     "Ifpack2::SparseContainer::weightedApply: "
00412     "Y_local has length " << X_local->getLocalLength () << ", which does "
00413     "not match numRows_ = " << numRows_ << ".  Please report this bug to "
00414     "the Ifpack2 developers.");
00415   mvgs.gather (*Y_local, Y, localRows);
00416 
00417   // Apply the diagonal scaling D to the input X.  It's our choice
00418   // whether the result has the original input Map of X, or the
00419   // permuted subset Map of X_local.  If the latter, we also need to
00420   // gather D into the permuted subset Map.  We choose the latter, to
00421   // save memory and computation.  Thus, we do the following:
00422   //
00423   // 1. Gather D into a temporary vector D_local.
00424   // 2. Create a temporary X_scaled to hold diag(D_local) * X_local.
00425   // 3. Compute X_scaled := diag(D_loca) * X_local.
00426 
00427   RCP<V_inv> D_local = rcp (new V_inv (Inverse_->getDomainMap ()));
00428   TEUCHOS_TEST_FOR_EXCEPTION(
00429     D_local->getLocalLength () != numRows_, std::logic_error,
00430     "Ifpack2::SparseContainer::weightedApply: "
00431     "D_local has length " << X_local->getLocalLength () << ", which does "
00432     "not match numRows_ = " << numRows_ << ".  Please report this bug to "
00433     "the Ifpack2 developers.");
00434   mvgs.gather (*D_local, D, localRows);
00435   RCP<MV_inv> X_scaled = rcp (new MV_inv (Inverse_->getDomainMap (), numVecs));
00436   X_scaled->elementWiseMultiply (STS::one (), *D_local, *X_local, STS::zero ());
00437 
00438   // Y_temp will hold the result of M^{-1}*X_scaled.  If beta == 0, we
00439   // can write the result of Inverse_->apply() directly to Y_local, so
00440   // Y_temp may alias Y_local.  Otherwise, if beta != 0, we need
00441   // temporary storage for M^{-1}*X_scaled, so Y_temp must be
00442   // different than Y_local.
00443   RCP<MV_inv> Y_temp;
00444   if (beta == STS::zero ()) {
00445     Y_temp = Y_local;
00446   } else {
00447     Y_temp = rcp (new MV_inv (Inverse_->getRangeMap (), numVecs));
00448   }
00449 
00450   // Apply the local operator: Y_temp := M^{-1} * X_scaled
00451   Inverse_->apply (*X_scaled, *Y_temp, mode);
00452   // Y_local := beta * Y_local + alpha * diag(D_local) * Y_temp.
00453   //
00454   // Note that we still use the permuted subset scaling D_local here,
00455   // because Y_temp has the same permuted subset Map.  That's good, in
00456   // fact, because it's a subset: less data to read and multiply.
00457   Y_local->elementWiseMultiply (alpha, *D_local, *Y_temp, beta);
00458 
00459   // Copy the permuted subset output vector Y_local into the original
00460   // output multivector Y.
00461   mvgs.scatter (Y, *Y_local, localRows);
00462 }
00463 
00464 //==============================================================================
00465 template<class MatrixType, class InverseType>
00466 std::ostream& SparseContainer<MatrixType,InverseType>::print(std::ostream& os) const
00467 {
00468   Teuchos::FancyOStream fos(Teuchos::rcp(&os,false));
00469   fos.setOutputToRootOnly(0);
00470   describe(fos);
00471   return(os);
00472 }
00473 
00474 //==============================================================================
00475 template<class MatrixType, class InverseType>
00476 std::string SparseContainer<MatrixType,InverseType>::description() const
00477 {
00478   std::ostringstream oss;
00479   oss << Teuchos::Describable::description();
00480   if (isInitialized()) {
00481     if (isComputed()) {
00482       oss << "{status = initialized, computed";
00483     }
00484     else {
00485       oss << "{status = initialized, not computed";
00486     }
00487   }
00488   else {
00489     oss << "{status = not initialized, not computed";
00490   }
00491 
00492   oss << "}";
00493   return oss.str();
00494 }
00495 
00496 //==============================================================================
00497 template<class MatrixType, class InverseType>
00498 void SparseContainer<MatrixType,InverseType>::describe(Teuchos::FancyOStream &os, const Teuchos::EVerbosityLevel verbLevel) const
00499 {
00500   using std::endl;
00501   if(verbLevel==Teuchos::VERB_NONE) return;
00502   os << "================================================================================" << endl;
00503   os << "Ifpack2::SparseContainer" << endl;
00504   os << "Number of rows          = " << numRows_ << endl;
00505   os << "isInitialized()         = " << IsInitialized_ << endl;
00506   os << "isComputed()            = " << IsComputed_ << endl;
00507   os << "================================================================================" << endl;
00508   os << endl;
00509 }
00510 
00511 //==============================================================================
00512 template<class MatrixType, class InverseType>
00513 void SparseContainer<MatrixType,InverseType>::
00514 extract (const Teuchos::RCP<const row_matrix_type>& globalMatrix)
00515 {
00516   using Teuchos::Array;
00517   using Teuchos::ArrayView;
00518 
00519   const size_t MatrixInNumRows = globalMatrix->getNodeNumRows ();
00520 
00521   // Sanity checking
00522   ArrayView<const local_ordinal_type> localRows = this->getLocalRows ();
00523   for (size_t j = 0; j < numRows_; ++j) {
00524     TEUCHOS_TEST_FOR_EXCEPTION(
00525       localRows[j] < 0 ||
00526       static_cast<size_t> (localRows[j]) >= MatrixInNumRows,
00527       std::runtime_error, "Ifpack2::SparseContainer::extract: localRows[j="
00528       << j << "] = " << localRows[j] << ", which is out of the valid range.  "
00529       "This probably means that compute() has not yet been called.");
00530   }
00531 
00532   const size_t maxNumEntriesInRow = globalMatrix->getNodeMaxNumRowEntries();
00533   Array<scalar_type>         Values;
00534   Array<local_ordinal_type>   Indices;
00535   Array<InverseScalar>        Values_insert;
00536   Array<InverseGlobalOrdinal> Indices_insert;
00537 
00538   Values.resize (maxNumEntriesInRow);
00539   Indices.resize (maxNumEntriesInRow);
00540   Values_insert.resize (maxNumEntriesInRow);
00541   Indices_insert.resize (maxNumEntriesInRow);
00542 
00543   const InverseLocalOrdinal INVALID =
00544     Teuchos::OrdinalTraits<InverseLocalOrdinal>::invalid ();
00545   for (size_t j = 0; j < numRows_; ++j) {
00546     const local_ordinal_type localRow = localRows[j];
00547     size_t numEntries;
00548     globalMatrix->getLocalRowCopy (localRow, Indices (), Values (), numEntries);
00549 
00550     size_t num_entries_found = 0;
00551     for (size_t k = 0; k < numEntries; ++k) {
00552       const local_ordinal_type localCol = Indices[k];
00553 
00554       // Skip off-process elements
00555       //
00556       // FIXME (mfh 24 Aug 2013) This assumes the following:
00557       //
00558       // 1. The column and row Maps begin with the same set of
00559       //    on-process entries, in the same order.  That is,
00560       //    on-process row and column indices are the same.
00561       // 2. All off-process indices in the column Map of the input
00562       //    matrix occur after that initial set.
00563       if (static_cast<size_t> (localCol) >= MatrixInNumRows) {
00564         continue;
00565       }
00566       // for local column IDs, look for each ID in the list
00567       // of columns hosted by this object
00568       InverseLocalOrdinal jj = INVALID;
00569       for (size_t kk = 0; kk < numRows_; ++kk) {
00570         if (localRows[kk] == localCol) {
00571           jj = kk;
00572         }
00573       }
00574 
00575       if (jj != INVALID) {
00576         Indices_insert[num_entries_found] = localMap_->getGlobalElement (jj);
00577         Values_insert[num_entries_found] = Values[k];
00578         num_entries_found++;
00579       }
00580     }
00581     diagBlock_->insertGlobalValues (j, Indices_insert (0, num_entries_found),
00582                                     Values_insert (0, num_entries_found));
00583   }
00584 
00585   // FIXME (mfh 24 Aug 2013) If we generalize the current set of
00586   // assumptions on the column and row Maps (see note above), we may
00587   // need to supply custom domain and range Maps to fillComplete().
00588   diagBlock_->fillComplete ();
00589 }
00590 
00591 } // namespace Ifpack2
00592 
00593 // For ETI
00594 #include "Ifpack2_ILUT.hpp"
00595 
00596 #define IFPACK2_SPARSECONTAINER_INSTANT(S,LO,GO,N) \
00597   template class Ifpack2::SparseContainer< Tpetra::CrsMatrix<S, LO, GO, N>, \
00598                                   Ifpack2::ILUT<Tpetra::CrsMatrix<S,LO,GO,N> > \
00599                                   >;
00600 
00601 #endif // IFPACK2_SPARSECONTAINER_HPP
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends