Tpetra Matrix/Vector Services Version of the Day
Tpetra_CrsGraph_def.hpp
00001 // @HEADER
00002 // ***********************************************************************
00003 //
00004 //          Tpetra: Templated Linear Algebra Services Package
00005 //                 Copyright (2008) Sandia Corporation
00006 //
00007 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
00008 // the U.S. Government retains certain rights in this software.
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 // FINISH: can't use rowPtrs_ without checking that it exists
00043 // FINISH: add code to fillComplete() and CrsMatrix::fillComplete() to delete the Tpetra data
00044 
00045 #ifndef TPETRA_CRSGRAPH_DEF_HPP
00046 #define TPETRA_CRSGRAPH_DEF_HPP
00047 
00048 #include <Kokkos_NodeTrace.hpp>
00049 #include <Teuchos_Assert.hpp>
00050 #include <Teuchos_NullIteratorTraits.hpp>
00051 #include <Teuchos_as.hpp>
00052 #include <algorithm>
00053 #include <string>
00054 #include <utility>
00055 #include <Teuchos_SerialDenseMatrix.hpp>
00056 
00057 #ifdef DOXYGEN_USE_ONLY
00058   #include "Tpetra_CrsGraph_decl.hpp"
00059 #endif
00060 
00061 namespace Tpetra {
00062 
00063   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00064   CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::
00065   CrsGraph (const RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > &rowMap,
00066             size_t maxNumEntriesPerRow,
00067             ProfileType pftype,
00068             const RCP<ParameterList>& params)
00069   : DistObject<GlobalOrdinal,LocalOrdinal,GlobalOrdinal,Node>(rowMap)
00070   , rowMap_(rowMap)
00071   , nodeNumEntries_(0)
00072   , nodeNumAllocated_(OrdinalTraits<size_t>::invalid())
00073   , pftype_(pftype)
00074   , numAllocForAllRows_(maxNumEntriesPerRow)
00075   , indicesAreAllocated_(false)
00076   , indicesAreLocal_(false)
00077   , indicesAreGlobal_(false)
00078   , haveRowInfo_(true)
00079   , insertGlobalIndicesWarnedEfficiency_(false)
00080   , insertLocalIndicesWarnedEfficiency_(false)
00081   {
00082     typedef Teuchos::OrdinalTraits<size_t> OTST;
00083     staticAssertions();
00084     TEUCHOS_TEST_FOR_EXCEPTION(maxNumEntriesPerRow == OTST::invalid(),
00085       std::invalid_argument, "The allocation hint must be a valid size_t value, "
00086       "which in this case means it must not be Teuchos::OrdinalTraits<size_t>::"
00087       "invalid().");
00088     resumeFill(params);
00089     checkInternalState();
00090   }
00091 
00092 
00093   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00094   CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::
00095   CrsGraph (const RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > &rowMap,
00096             const RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > &colMap,
00097             size_t maxNumEntriesPerRow,
00098             ProfileType pftype,
00099             const RCP<ParameterList>& params)
00100   : DistObject<GlobalOrdinal,LocalOrdinal,GlobalOrdinal,Node>(rowMap)
00101   , rowMap_(rowMap)
00102   , colMap_(colMap)
00103   , nodeNumEntries_(0)
00104   , nodeNumAllocated_(OrdinalTraits<size_t>::invalid())
00105   , pftype_(pftype)
00106   , numAllocForAllRows_(maxNumEntriesPerRow)
00107   , indicesAreAllocated_(false)
00108   , indicesAreLocal_(false)
00109   , indicesAreGlobal_(false)
00110   , haveRowInfo_(true)
00111   , insertGlobalIndicesWarnedEfficiency_(false)
00112   , insertLocalIndicesWarnedEfficiency_(false)
00113   {
00114     typedef Teuchos::OrdinalTraits<size_t> OTST;
00115     staticAssertions();
00116     TEUCHOS_TEST_FOR_EXCEPTION(maxNumEntriesPerRow == OTST::invalid(),
00117       std::invalid_argument, "The allocation hint must be a valid size_t value, "
00118       "which in this case means it must not be Teuchos::OrdinalTraits<size_t>::"
00119       "invalid().");
00120     resumeFill(params);
00121     checkInternalState();
00122   }
00123 
00124 
00125   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00126   CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::
00127   CrsGraph (const RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > &rowMap,
00128             const ArrayRCP<const size_t> &NumEntriesPerRowToAlloc,
00129             ProfileType pftype,
00130             const RCP<ParameterList>& params)
00131   : DistObject<GlobalOrdinal,LocalOrdinal,GlobalOrdinal,Node>(rowMap)
00132   , rowMap_(rowMap)
00133   , nodeNumEntries_(0)
00134   , nodeNumAllocated_(OrdinalTraits<size_t>::invalid())
00135   , pftype_(pftype)
00136   , numAllocPerRow_(NumEntriesPerRowToAlloc)
00137   , numAllocForAllRows_(0)
00138   , indicesAreAllocated_(false)
00139   , indicesAreLocal_(false)
00140   , indicesAreGlobal_(false)
00141   , haveRowInfo_(true)
00142   , insertGlobalIndicesWarnedEfficiency_(false)
00143   , insertLocalIndicesWarnedEfficiency_(false)
00144   {
00145     typedef Teuchos::OrdinalTraits<size_t> OTST;
00146     std::string tfecfFuncName("CrsGraph(rowMap,NumEntriesPerRowToAlloc)");
00147     staticAssertions();
00148     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC((size_t)NumEntriesPerRowToAlloc.size() != getNodeNumRows(), std::invalid_argument,
00149         ": NumEntriesPerRowToAlloc must have as many entries as specified by rowMap for this node.");
00150     for (size_t r=0; r < getNodeNumRows(); ++r) {
00151       const size_t curRowCount = NumEntriesPerRowToAlloc[r];
00152       TEUCHOS_TEST_FOR_EXCEPTION(curRowCount == OTST::invalid(),
00153         std::invalid_argument, "NumEntriesPerRowToAlloc[" << r << "] specifies "
00154         "an invalid number of entries (Teuchos::OrdinalTraits<size_t>::"
00155         "invalid()).");
00156     }
00157     resumeFill(params);
00158     checkInternalState();
00159   }
00160 
00161 
00162   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00163   CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::
00164   CrsGraph (const RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > &rowMap,
00165             const RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > &colMap,
00166             const ArrayRCP<const size_t> &NumEntriesPerRowToAlloc,
00167             ProfileType pftype,
00168             const RCP<ParameterList>& params)
00169   : DistObject<GlobalOrdinal,LocalOrdinal,GlobalOrdinal,Node>(rowMap)
00170   , rowMap_(rowMap)
00171   , colMap_(colMap)
00172   , nodeNumEntries_(0)
00173   , nodeNumAllocated_(OrdinalTraits<size_t>::invalid())
00174   , pftype_(pftype)
00175   , numAllocPerRow_(NumEntriesPerRowToAlloc)
00176   , numAllocForAllRows_(0)
00177   , indicesAreAllocated_(false)
00178   , indicesAreLocal_(false)
00179   , indicesAreGlobal_(false)
00180   , haveRowInfo_(true)
00181   , insertGlobalIndicesWarnedEfficiency_(false)
00182   , insertLocalIndicesWarnedEfficiency_(false)
00183   {
00184     typedef Teuchos::OrdinalTraits<size_t> OTST;
00185     std::string tfecfFuncName("CrsGraph(rowMap,colMap,NumEntriesPerRowToAlloc)");
00186     staticAssertions();
00187     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC((size_t)NumEntriesPerRowToAlloc.size() != getNodeNumRows(), std::invalid_argument,
00188         ": NumEntriesPerRowToAlloc must have as many entries as specified by rowMap for this node.");
00189     for (size_t r=0; r < getNodeNumRows(); ++r) {
00190       const size_t curRowCount = NumEntriesPerRowToAlloc[r];
00191       TEUCHOS_TEST_FOR_EXCEPTION(curRowCount == OTST::invalid(),
00192         std::invalid_argument, "NumEntriesPerRowToAlloc[" << r << "] specifies "
00193         "an invalid number of entries (Teuchos::OrdinalTraits<size_t>::"
00194         "invalid()).");
00195     }
00196     resumeFill(params);
00197     checkInternalState();
00198   }
00199 
00200 
00201   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00202   CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::~CrsGraph()
00203   {}
00204 
00205   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00206   RCP<const ParameterList>
00207   CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::
00208   getValidParameters () const
00209   {
00210     RCP<ParameterList> params = parameterList ("Tpetra::CrsGraph");
00211 
00212     // Make a sublist for the Import.
00213     RCP<ParameterList> importSublist = parameterList ("Import");
00214 
00215     // FIXME (mfh 02 Apr 2012) We should really have the Import and
00216     // Export objects fill in these lists.  However, we don't want to
00217     // create an Import or Export unless we need them.  For now, we
00218     // know that the Import and Export just pass the list directly to
00219     // their Distributor, so we can create a Distributor here
00220     // (Distributor's constructor is a lightweight operation) and have
00221     // it fill in the list.
00222 
00223     // Fill in Distributor default parameters by creating a
00224     // Distributor and asking it to do the work.
00225     Distributor distributor (rowMap_->getComm(), importSublist);
00226     params->set ("Import", *importSublist, "How the Import performs communication.");
00227 
00228     // Make a sublist for the Export.  For now, it's a clone of the
00229     // Import sublist.  It's not a shallow copy, though, since we
00230     // might like the Import to do communication differently than the
00231     // Export.
00232     params->set ("Export", *importSublist, "How the Export performs communication.");
00233 
00234     return params;
00235   }
00236 
00237   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00238   void
00239   CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::
00240   setParameterList (const RCP<ParameterList>& params)
00241   {
00242     RCP<const ParameterList> validParams = getValidParameters ();
00243     params->validateParametersAndSetDefaults (*validParams);
00244     this->setMyParamList (params);
00245   }
00246 
00247   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00248   global_size_t CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getGlobalNumRows() const
00249   {
00250     return rowMap_->getGlobalNumElements();
00251   }
00252 
00253   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00254   global_size_t CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getGlobalNumCols() const
00255   {
00256     std::string tfecfFuncName("getGlobalNumCols()");
00257     TEUCHOS_TEST_FOR_EXCEPTION(isFillComplete() == false, std::runtime_error, ": requires domain map (which requires fillComplete()).");
00258     return getDomainMap()->getGlobalNumElements();
00259   }
00260 
00261   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00262   size_t CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getNodeNumRows() const
00263   {
00264     return rowMap_->getNodeNumElements();
00265   }
00266 
00267   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00268   size_t CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getNodeNumCols() const
00269   {
00270     std::string tfecfFuncName("getNodeNumCols()");
00271     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(hasColMap() != true, std::runtime_error, ": requires column map.");
00272     return colMap_->getNodeNumElements();
00273   }
00274 
00275   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00276   size_t CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getNodeNumDiags() const
00277   {
00278     return nodeNumDiags_;
00279   }
00280 
00281   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00282   global_size_t CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getGlobalNumDiags() const
00283   {
00284     return globalNumDiags_;
00285   }
00286 
00287   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00288   RCP<Node>
00289   CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getNode() const
00290   {
00291     return rowMap_->getNode();
00292   }
00293 
00294   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00295   const RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > &
00296   CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getRowMap() const
00297   {
00298     return rowMap_;
00299   }
00300 
00301   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00302   const RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > &
00303   CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getColMap() const
00304   {
00305     return colMap_;
00306   }
00307 
00308   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00309   const RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > &
00310   CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getDomainMap() const
00311   {
00312     return domainMap_;
00313   }
00314 
00315   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00316   const RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > &
00317   CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getRangeMap() const
00318   {
00319     return rangeMap_;
00320   }
00321 
00322   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00323   RCP<const Import<LocalOrdinal,GlobalOrdinal,Node> >
00324   CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getImporter() const
00325   {
00326     return importer_;
00327   }
00328 
00329   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00330   RCP<const Export<LocalOrdinal,GlobalOrdinal,Node> >
00331   CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getExporter() const
00332   {
00333     return exporter_;
00334   }
00335 
00336   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00337   bool CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::hasColMap() const
00338   {
00339     return (colMap_ != null);
00340   }
00341 
00342   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00343   bool CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::isStorageOptimized() const
00344   {
00345     bool isOpt = indicesAreAllocated_ && (numRowEntries_ == null) && (getNodeNumRows() > 0);
00346 #ifdef HAVE_TPETRA_DEBUG
00347     std::string tfecfFuncName("isStorageOptimized()");
00348     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC( (isOpt == true) && (getProfileType() == DynamicProfile), std::logic_error,
00349         ": Violated stated post-conditions. Please contact Tpetra team.");
00350 #endif
00351     return isOpt;
00352   }
00353 
00354   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00355   ProfileType CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getProfileType() const
00356   {
00357     return pftype_;
00358   }
00359 
00360   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00361   global_size_t CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getGlobalNumEntries() const
00362   {
00363     return globalNumEntries_;
00364   }
00365 
00366   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00367   size_t CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getNodeNumEntries() const
00368   {
00369     return nodeNumEntries_;
00370   }
00371 
00372   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00373   global_size_t CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getGlobalMaxNumRowEntries() const
00374   {
00375     return globalMaxNumRowEntries_;
00376   }
00377 
00378   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00379   size_t CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getNodeMaxNumRowEntries() const
00380   {
00381     return nodeMaxNumRowEntries_;
00382   }
00383 
00384   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00385   bool CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::isFillComplete() const
00386   {
00387     return fillComplete_;
00388   }
00389 
00390   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00391   bool CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::isFillActive() const
00392   {
00393     return !fillComplete_;
00394   }
00395 
00396   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00397   bool CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::isUpperTriangular() const
00398   {
00399     return upperTriangular_;
00400   }
00401 
00402   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00403   bool CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::isLowerTriangular() const
00404   {
00405     return lowerTriangular_;
00406   }
00407 
00408   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00409   bool CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::isLocallyIndexed() const
00410   {
00411     return indicesAreLocal_;
00412   }
00413 
00414   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00415   bool CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::isGloballyIndexed() const
00416   {
00417     return indicesAreGlobal_;
00418   }
00419 
00420   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00421   size_t CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getNodeAllocationSize() const
00422   {
00423     return nodeNumAllocated_;
00424   }
00425 
00426   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00427   const RCP<const Comm<int> > &
00428   CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getComm() const
00429   {
00430     return rowMap_->getComm();
00431   }
00432 
00433   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00434   GlobalOrdinal CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getIndexBase() const
00435   {
00436     return rowMap_->getIndexBase();
00437   }
00438 
00439   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00440   bool CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::indicesAreAllocated() const
00441   {
00442     return indicesAreAllocated_;
00443   }
00444 
00445 
00448   //                                                                         //
00449   //                    Internal utility methods                             //
00450   //                                                                         //
00453 
00454 
00457   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00458   bool CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::isSorted() const
00459   {
00460     return indicesAreSorted_;
00461   }
00462 
00463 
00466   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00467   bool CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::isMerged() const
00468   {
00469     return noRedundancies_;
00470   }
00471 
00472 
00475   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00476   void CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::setSorted(bool sorted)
00477   {
00478     indicesAreSorted_ = sorted;
00479   }
00480 
00481 
00484   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00485   void CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::setMerged(bool merged)
00486   {
00487     noRedundancies_ = merged;
00488   }
00489 
00490 
00493   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00494   void CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::allocateIndices(ELocalGlobal lg)
00495   {
00496     // This is a protected function, only callable by us.  If it was
00497     // called incorrectly, it is our fault.  That's why the tests
00498     // below throw std::logic_error instead of std::invalid_argument.
00499     std::string tfecfFuncName("allocateIndices()");
00500     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
00501       isLocallyIndexed() && lg==GlobalIndices, std::logic_error,
00502       ": The graph is locally indexed, but Tpetra code is calling this method "
00503       "with lg=GlobalIndices.  Please report this bug to the Tpetra developers.");
00504     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
00505       isGloballyIndexed() && lg==LocalIndices, std::logic_error,
00506       ": The graph is globally indexed, but Tpetra code is calling this method "
00507       "with lg=LocalIndices.  Please report this bug to the Tpetra developers.");
00508     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
00509       indicesAreAllocated(), std::logic_error, ": The graph's indices are "
00510       "already allocated, but Tpetra code is calling allocateIndices() again.  "
00511       "Please report this bug to the Tpetra developers.");
00512 
00513     const size_t numRows = getNodeNumRows();
00514     indicesAreLocal_  = (lg == LocalIndices);
00515     indicesAreGlobal_ = (lg == GlobalIndices);
00516     nodeNumAllocated_ = 0;
00517     if (numAllocPerRow_ == null && getNodeNumRows() > 0) {
00518       // this wastes memory, temporarily, but it simplifies the code
00519       // and interfaces to follow
00520       ArrayRCP<size_t> tmpnumallocperrow = arcp<size_t>(numRows);
00521       std::fill(tmpnumallocperrow.begin(), tmpnumallocperrow.end(), numAllocForAllRows_);
00522       numAllocPerRow_ = tmpnumallocperrow;
00523     }
00524     //
00525     if (getProfileType() == StaticProfile) {
00526       //
00527       //  STATIC ALLOCATION PROFILE
00528       //
00529       // Have the local sparse kernels object allocate row offsets for
00530       // us, with first-touch allocation if applicable.  This is not
00531       // as important for global indices, because we never use global
00532       // indices in sparse kernels, but we might as well use the code
00533       // that we have for both the local and global indices cases.
00534       // Furthermore, first-touch allocation ensures that we don't
00535       // take up too much memory in any one NUMA affinity region.
00536       rowPtrs_ = LocalMatOps::allocRowPtrs( getRowMap()->getNode(), numAllocPerRow_() );
00537       if (lg == LocalIndices) {
00538         lclInds1D_ = LocalMatOps::template allocStorage<LocalOrdinal>( getRowMap()->getNode(), rowPtrs_() );
00539       }
00540       else {
00541         gblInds1D_ = LocalMatOps::template allocStorage<GlobalOrdinal>( getRowMap()->getNode(), rowPtrs_() );
00542       }
00543       nodeNumAllocated_ = rowPtrs_[numRows];
00544     }
00545     else {
00546       //
00547       //  DYNAMIC ALLOCATION PROFILE
00548       //
00549       typename ArrayRCP<const size_t>::iterator numalloc = numAllocPerRow_.begin();
00550       size_t howmany = numAllocForAllRows_;
00551       if (lg == LocalIndices) {
00552         lclInds2D_ = arcp< ArrayRCP<LocalOrdinal> >(numRows);
00553         nodeNumAllocated_ = 0;
00554         for (size_t i=0; i < numRows; ++i) {
00555           if (numAllocPerRow_ != null) howmany = *numalloc++;
00556           nodeNumAllocated_ += howmany;
00557           if (howmany > 0) lclInds2D_[i] = arcp<LocalOrdinal>(howmany);
00558         }
00559       }
00560       else { // allocate global indices
00561         gblInds2D_ = arcp< ArrayRCP<GlobalOrdinal> >(numRows);
00562         nodeNumAllocated_ = 0;
00563         for (size_t i=0; i < numRows; ++i) {
00564           if (numAllocPerRow_ != null) howmany = *numalloc++;
00565           nodeNumAllocated_ += howmany;
00566           if (howmany > 0) gblInds2D_[i] = arcp<GlobalOrdinal>(howmany);
00567         }
00568       }
00569     }
00570     if (numRows > 0) {
00571       numRowEntries_ = arcp<size_t>(numRows);
00572       std::fill( numRowEntries_.begin(), numRowEntries_.end(), 0 );
00573     }
00574     // done with these
00575     numAllocForAllRows_ = 0;
00576     numAllocPerRow_     = null;
00577     indicesAreAllocated_ = true;
00578     checkInternalState();
00579   }
00580 
00581 
00584   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00585   template <class T>
00586   ArrayRCP<T> CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::allocateValues1D() const
00587   {
00588     std::string tfecfFuncName("allocateValues()");
00589     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
00590         indicesAreAllocated() == false,
00591         std::runtime_error, ": graph indices must be allocated before values."
00592     );
00593     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
00594         getProfileType() != StaticProfile,
00595         std::runtime_error, ": graph indices must be allocated in a static profile."
00596     );
00597     ArrayRCP<T> values1D;
00598     values1D = LocalMatOps::template allocStorage<T>( getRowMap()->getNode(), rowPtrs_() );
00599     return values1D;
00600   }
00601 
00602 
00605   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00606   template <class T>
00607   ArrayRCP<ArrayRCP<T> > CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::allocateValues2D() const
00608   {
00609     std::string tfecfFuncName("allocateValues2D()");
00610     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
00611         indicesAreAllocated() == false,
00612         std::runtime_error, ": graph indices must be allocated before values."
00613     );
00614     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
00615         getProfileType() != DynamicProfile,
00616         std::runtime_error, ": graph indices must be allocated in a dynamic profile."
00617     );
00618     ArrayRCP<ArrayRCP<T> > values2D;
00619     values2D = arcp<ArrayRCP<T> >(getNodeNumRows());
00620     if (lclInds2D_ != null) {
00621       const size_t numRows = lclInds2D_.size();
00622       for (size_t r=0; r != numRows; ++r) {
00623         if (lclInds2D_[r] != null) {
00624           values2D[r] = arcp<T>(lclInds2D_[r].size());
00625         }
00626       }
00627     }
00628     else if (gblInds2D_ != null) {
00629       const size_t numRows = gblInds2D_.size();
00630       for (size_t r=0; r != numRows; ++r) {
00631         if (gblInds2D_[r] != null) {
00632           values2D[r] = arcp<T>(gblInds2D_[r].size());
00633         }
00634       }
00635     }
00636     return values2D;
00637   }
00638 
00639 
00642   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00643   template <ELocalGlobal lg>
00644   RowInfo CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::updateAlloc(RowInfo rowinfo, size_t newAllocSize)
00645   {
00646 #ifdef HAVE_TPETRA_DEBUG
00647     TEUCHOS_TEST_FOR_EXCEPT( rowMap_->isNodeLocalElement(rowinfo.localRow) == false );
00648     TEUCHOS_TEST_FOR_EXCEPT( newAllocSize < rowinfo.allocSize );
00649     TEUCHOS_TEST_FOR_EXCEPT( (lg == LocalIndices && isLocallyIndexed() == false) || (lg == GlobalIndices && isGloballyIndexed() == false) );
00650     TEUCHOS_TEST_FOR_EXCEPT( newAllocSize == 0 );
00651     TEUCHOS_TEST_FOR_EXCEPT( indicesAreAllocated() == false );
00652 #endif
00653     // allocate a larger space for row "lrow"
00654     // copy any existing data from previous allocation to new allocation
00655     // update sizes
00656     if (lg == LocalIndices) {
00657       ArrayRCP<LocalOrdinal> old_alloc = lclInds2D_[rowinfo.localRow];
00658       lclInds2D_[rowinfo.localRow] = arcp<LocalOrdinal>(newAllocSize);
00659       std::copy(old_alloc.begin(), old_alloc.begin() + rowinfo.numEntries, lclInds2D_[rowinfo.localRow].begin());
00660     }
00661     else /* if lg == GlobalIndices */ {
00662       ArrayRCP<GlobalOrdinal> old_alloc = gblInds2D_[rowinfo.localRow];
00663       gblInds2D_[rowinfo.localRow] = arcp<GlobalOrdinal>(newAllocSize);
00664       std::copy(old_alloc.begin(), old_alloc.begin() + rowinfo.numEntries, gblInds2D_[rowinfo.localRow].begin());
00665     }
00666     //
00667     nodeNumAllocated_ += (newAllocSize - rowinfo.allocSize);
00668     rowinfo.allocSize = newAllocSize;
00669     return rowinfo;
00670   }
00671 
00672 
00675   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00676   template <ELocalGlobal lg, class T>
00677   RowInfo CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::updateAllocAndValues(RowInfo rowinfo, size_t newAllocSize, ArrayRCP<T> &rowVals)
00678   {
00679 #ifdef HAVE_TPETRA_DEBUG
00680     TEUCHOS_TEST_FOR_EXCEPT( ! rowMap_->isNodeLocalElement(rowinfo.localRow) );
00681     TEUCHOS_TEST_FOR_EXCEPT( newAllocSize < rowinfo.allocSize );
00682     TEUCHOS_TEST_FOR_EXCEPT( (lg == LocalIndices && ! isLocallyIndexed()) || (lg == GlobalIndices && ! isGloballyIndexed()) );
00683     TEUCHOS_TEST_FOR_EXCEPT( newAllocSize == 0 );
00684     TEUCHOS_TEST_FOR_EXCEPT( ! indicesAreAllocated() );
00685 #endif
00686     // allocate a larger space for row "lrow"
00687     // copy any existing data from previous allocation to new allocation
00688     // update sizes
00689     if (lg == LocalIndices) {
00690       ArrayRCP<LocalOrdinal> old_alloc = lclInds2D_[rowinfo.localRow];
00691       lclInds2D_[rowinfo.localRow] = arcp<LocalOrdinal>(newAllocSize);
00692       std::copy(old_alloc.begin(), old_alloc.begin() + rowinfo.numEntries, lclInds2D_[rowinfo.localRow].begin());
00693     }
00694     else /* if lg == GlobalIndices */ {
00695       ArrayRCP<GlobalOrdinal> old_alloc = gblInds2D_[rowinfo.localRow];
00696       gblInds2D_[rowinfo.localRow] = arcp<GlobalOrdinal>(newAllocSize);
00697       std::copy(old_alloc.begin(), old_alloc.begin() + rowinfo.numEntries, gblInds2D_[rowinfo.localRow].begin());
00698     }
00699     ArrayRCP<const T> oldVals = rowVals;
00700     rowVals = arcp<T>(newAllocSize);
00701     std::copy(oldVals.begin(), oldVals.begin() + rowinfo.numEntries, rowVals.begin());
00702     //
00703     nodeNumAllocated_ += (newAllocSize - rowinfo.allocSize);
00704     rowinfo.allocSize = newAllocSize;
00705     return rowinfo;
00706   }
00707 
00708 
00711   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00712   ArrayView<const LocalOrdinal> CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getLocalView(RowInfo rowinfo) const
00713   {
00714     ArrayView<const LocalOrdinal> view;
00715     if (rowinfo.allocSize > 0) {
00716       if (lclInds1D_ != null) {
00717         view = lclInds1D_(rowinfo.offset1D,rowinfo.allocSize);
00718       }
00719       else if (lclInds2D_[rowinfo.localRow] != null) {
00720         view = lclInds2D_[rowinfo.localRow]();
00721       }
00722     }
00723     return view;
00724   }
00725 
00726 
00729   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00730   ArrayView<LocalOrdinal> CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getLocalViewNonConst(RowInfo rowinfo)
00731   {
00732     ArrayView<LocalOrdinal> view;
00733     if (rowinfo.allocSize > 0) {
00734       if (lclInds1D_ != null) {
00735         view = lclInds1D_(rowinfo.offset1D,rowinfo.allocSize);
00736       }
00737       else if (lclInds2D_[rowinfo.localRow] != null) {
00738         view = lclInds2D_[rowinfo.localRow]();
00739       }
00740     }
00741     return view;
00742   }
00743 
00744 
00747   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00748   ArrayView<const GlobalOrdinal> CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getGlobalView(RowInfo rowinfo) const
00749   {
00750     ArrayView<const GlobalOrdinal> view;
00751     if (rowinfo.allocSize > 0) {
00752       if (gblInds1D_ != null) {
00753         view = gblInds1D_(rowinfo.offset1D,rowinfo.allocSize);
00754       }
00755       else if (gblInds2D_[rowinfo.localRow] != null) {
00756         view = gblInds2D_[rowinfo.localRow]();
00757       }
00758     }
00759     return view;
00760   }
00761 
00762 
00765   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00766   ArrayView<GlobalOrdinal> CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getGlobalViewNonConst(RowInfo rowinfo)
00767   {
00768     ArrayView<GlobalOrdinal> view;
00769     if (rowinfo.allocSize > 0) {
00770       if (gblInds1D_ != null) {
00771         view = gblInds1D_(rowinfo.offset1D,rowinfo.allocSize);
00772       }
00773       else if (gblInds2D_[rowinfo.localRow] != null) {
00774         view = gblInds2D_[rowinfo.localRow]();
00775       }
00776     }
00777     return view;
00778   }
00779 
00780 
00783   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00784   RowInfo CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getRowInfo(size_t myRow) const
00785   {
00786 #ifdef HAVE_TPETRA_DEBUG
00787     std::string tfecfFuncName("getRowInfo()");
00788     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
00789         rowMap_->isNodeLocalElement(myRow) == false,
00790         std::logic_error, ": Internal logic error. Please contact Tpetra team."
00791     )
00792     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
00793         hasRowInfo() == false,
00794         std::logic_error, ": Late catch! Graph does not have row info anymore. Error should have been caught earlier. Please contact Tpetra team."
00795     )
00796 #endif
00797     const size_t STINV = OrdinalTraits<size_t>::invalid();
00798     RowInfo ret;
00799     ret.localRow = myRow;
00800     if (nodeNumAllocated_ != 0 && nodeNumAllocated_ != STINV) {
00801       // graph data structures have the info that we need
00802       //
00803       // if static graph, offsets tell us the allocation size
00804       if (getProfileType() == StaticProfile) {
00805         ret.offset1D   = rowPtrs_[myRow];
00806         ret.allocSize  = rowPtrs_[myRow+1] - rowPtrs_[myRow];
00807         if (numRowEntries_ == null) ret.numEntries = ret.allocSize;
00808         else                        ret.numEntries = numRowEntries_[myRow];
00809       }
00810       else {
00811         ret.offset1D = STINV;
00812         if (isLocallyIndexed()) {
00813           ret.allocSize = lclInds2D_[myRow].size();
00814         }
00815         else {
00816           ret.allocSize = gblInds2D_[myRow].size();
00817         }
00818         ret.numEntries = numRowEntries_[myRow];
00819       }
00820     }
00821     else if (nodeNumAllocated_ == 0) {
00822       // have performed allocation, but the graph has no allocation or entries
00823       ret.allocSize = 0;
00824       ret.numEntries = 0;
00825       ret.offset1D = STINV;
00826     }
00827     else if (indicesAreAllocated() == false) {
00828       // haven't performed allocation yet; probably won't hit this code
00829       if (numAllocPerRow_ == null) {
00830         ret.allocSize = numAllocForAllRows_;
00831       }
00832       else {
00833         ret.allocSize = numAllocPerRow_[myRow];
00834       }
00835       ret.numEntries = 0;
00836       ret.offset1D = STINV;
00837     }
00838     else {
00839       // don't know how we ended up here...
00840       TEUCHOS_TEST_FOR_EXCEPT(true);
00841     }
00842     return ret;
00843   }
00844 
00845 
00848   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00849   void CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::staticAssertions() const
00850   {
00851     // Assumption: sizeof(GlobalOrdinal) >= sizeof(LocalOrdinal)
00852     //    This is so that we can store LocalOrdinals in the memory formerly occupied by GlobalOrdinals
00853     // Assumption: max(GlobalOrdinal) >= max(LocalOrdinal)  and  max(size_t) >= max(LocalOrdinal)
00854     //    This is so that we can represent any LocalOrdinal as a size_t, and any LocalOrdinal as a GlobalOrdinal
00855     Teuchos::CompileTimeAssert<sizeof(GlobalOrdinal) < sizeof(LocalOrdinal)> cta_size1; (void)cta_size1;
00856     Teuchos::CompileTimeAssert<sizeof(global_size_t) < sizeof(size_t)      > cta_size2; (void)cta_size2;
00857     // can't call max() with CompileTimeAssert, because it isn't a constant expression; will need to make this a runtime check
00858     std::string msg = typeName(*this) + ": Object cannot be allocated with stated template arguments: size assumptions are not valid.";
00859     TEUCHOS_TEST_FOR_EXCEPTION( (size_t)OrdinalTraits<LocalOrdinal>::max() > OrdinalTraits<size_t>::max(),          std::runtime_error, msg);
00860     TEUCHOS_TEST_FOR_EXCEPTION( (global_size_t)OrdinalTraits<LocalOrdinal>::max() > (global_size_t)OrdinalTraits<GlobalOrdinal>::max(),           std::runtime_error, msg);
00861     TEUCHOS_TEST_FOR_EXCEPTION( (size_t)OrdinalTraits<GlobalOrdinal>::max() > OrdinalTraits<global_size_t>::max(),  std::runtime_error, msg);
00862     TEUCHOS_TEST_FOR_EXCEPTION( OrdinalTraits<size_t>::max() > OrdinalTraits<global_size_t>::max(),                 std::runtime_error, msg);
00863   }
00864 
00865 
00868   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00869   template <ELocalGlobal lg>
00870   size_t CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::filterIndices(const SLocalGlobalNCViews &inds) const
00871   {
00872     const Map<LocalOrdinal,GlobalOrdinal,Node> &cmap = *colMap_;
00873     Teuchos::CompileTimeAssert<lg != GlobalIndices && lg != LocalIndices> cta_lg;
00874     (void)cta_lg;
00875     size_t numFiltered = 0;
00876 #ifdef HAVE_TPETRA_DEBUG
00877     size_t numFiltered_debug = 0;
00878 #endif
00879     if (lg == GlobalIndices) {
00880       ArrayView<GlobalOrdinal> ginds = inds.ginds;
00881       typename ArrayView<GlobalOrdinal>::iterator fend = ginds.begin(),
00882                                                   cptr = ginds.begin();
00883       while (cptr != ginds.end()) {
00884         if (cmap.isNodeGlobalElement(*cptr)) {
00885           *fend++ = *cptr;
00886 #ifdef HAVE_TPETRA_DEBUG
00887           ++numFiltered_debug;
00888 #endif
00889         }
00890         ++cptr;
00891       }
00892       numFiltered = fend - ginds.begin();
00893     }
00894     else if (lg == LocalIndices) {
00895       ArrayView<LocalOrdinal> linds = inds.linds;
00896       typename ArrayView<LocalOrdinal>::iterator fend = linds.begin(),
00897                                                  cptr = linds.begin();
00898       while (cptr != linds.end()) {
00899         if (cmap.isNodeLocalElement(*cptr)) {
00900           *fend++ = *cptr;
00901 #ifdef HAVE_TPETRA_DEBUG
00902           ++numFiltered_debug;
00903 #endif
00904         }
00905         ++cptr;
00906       }
00907       numFiltered = fend - linds.begin();
00908     }
00909 #ifdef HAVE_TPETRA_DEBUG
00910     TEUCHOS_TEST_FOR_EXCEPT( numFiltered != numFiltered_debug );
00911 #endif
00912     return numFiltered;
00913   }
00914 
00915 
00918   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00919   template <ELocalGlobal lg, class T>
00920   size_t CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::filterIndicesAndValues(const SLocalGlobalNCViews &inds, const ArrayView<T> &vals) const
00921   {
00922     const Map<LocalOrdinal,GlobalOrdinal,Node> &cmap = *colMap_;
00923     Teuchos::CompileTimeAssert<lg != GlobalIndices && lg != LocalIndices> cta_lg;
00924     (void)cta_lg;
00925     size_t numFiltered = 0;
00926     typename ArrayView<T>::iterator fvalsend = vals.begin(),
00927                                     valscptr = vals.begin();
00928 #ifdef HAVE_TPETRA_DEBUG
00929     size_t numFiltered_debug = 0;
00930 #endif
00931     if (lg == GlobalIndices) {
00932       ArrayView<GlobalOrdinal> ginds = inds.ginds;
00933       typename ArrayView<GlobalOrdinal>::iterator fend = ginds.begin(),
00934                                                   cptr = ginds.begin();
00935       while (cptr != ginds.end()) {
00936         if (cmap.isNodeGlobalElement(*cptr)) {
00937           *fend++ = *cptr;
00938           *fvalsend++ = *valscptr;
00939 #ifdef HAVE_TPETRA_DEBUG
00940           ++numFiltered_debug;
00941 #endif
00942         }
00943         ++cptr;
00944         ++valscptr;
00945       }
00946       numFiltered = fend - ginds.begin();
00947     }
00948     else if (lg == LocalIndices) {
00949       ArrayView<LocalOrdinal> linds = inds.linds;
00950       typename ArrayView<LocalOrdinal>::iterator fend = linds.begin(),
00951                                                  cptr = linds.begin();
00952       while (cptr != linds.end()) {
00953         if (cmap.isNodeLocalElement(*cptr)) {
00954           *fend++ = *cptr;
00955           *fvalsend++ = *valscptr;
00956 #ifdef HAVE_TPETRA_DEBUG
00957           ++numFiltered_debug;
00958 #endif
00959         }
00960         ++cptr;
00961         ++valscptr;
00962       }
00963       numFiltered = fend - linds.begin();
00964     }
00965 #ifdef HAVE_TPETRA_DEBUG
00966     TEUCHOS_TEST_FOR_EXCEPT( numFiltered != numFiltered_debug );
00967     TEUCHOS_TEST_FOR_EXCEPT( valscptr != vals.end() );
00968     TEUCHOS_TEST_FOR_EXCEPT( numFiltered != (size_t)(fvalsend - vals.begin()) );
00969 #endif
00970     return numFiltered;
00971   }
00972 
00973 
00976   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00977   template <ELocalGlobal lg, ELocalGlobal I>
00978   size_t CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::insertIndices(RowInfo rowinfo, const SLocalGlobalViews &newInds)
00979   {
00980     Teuchos::CompileTimeAssert<lg != GlobalIndices && lg != LocalIndices> cta_lg; (void)cta_lg;
00981     size_t numNewInds = 0;
00982     if (lg == GlobalIndices) {
00983       ArrayView<const GlobalOrdinal> new_ginds = newInds.ginds;
00984       numNewInds = new_ginds.size();
00985       if (I == GlobalIndices) {
00986         ArrayView<GlobalOrdinal> gind_view = getGlobalViewNonConst(rowinfo);
00987         std::copy(new_ginds.begin(), new_ginds.end(), gind_view.begin()+rowinfo.numEntries);
00988       }
00989       else if (I == LocalIndices) {
00990         ArrayView<LocalOrdinal> lind_view = getLocalViewNonConst(rowinfo);
00991         typename ArrayView<const GlobalOrdinal>::iterator         in = new_ginds.begin();
00992         const typename ArrayView<const GlobalOrdinal>::iterator stop = new_ginds.end();
00993         typename ArrayView<LocalOrdinal>::iterator out = lind_view.begin()+rowinfo.numEntries;
00994         while (in != stop) {
00995           *out++ = colMap_->getLocalElement (*in++);
00996         }
00997       }
00998     }
00999     else if (lg == LocalIndices) {
01000       ArrayView<const LocalOrdinal> new_linds = newInds.linds;
01001       numNewInds = new_linds.size();
01002       if (I == LocalIndices) {
01003         ArrayView<LocalOrdinal> lind_view = getLocalViewNonConst(rowinfo);
01004         std::copy(new_linds.begin(), new_linds.end(), lind_view.begin()+rowinfo.numEntries);
01005       }
01006       else if (I == GlobalIndices) {
01007         // not needed yet
01008         TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error, "Tpetra::CrsGraph::"
01009           "insertIndices: the case where the input indices are local and the "
01010           "indices to write are global (lg=LocalIndices, I=GlobalIndices) has "
01011           "not yet been implemented.");
01012       }
01013     }
01014     numRowEntries_[rowinfo.localRow] += numNewInds;
01015     nodeNumEntries_ += numNewInds;
01016     setSorted(false);
01017     noRedundancies_ = false;
01018     return numNewInds;
01019   }
01020 
01021 
01024   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
01025   template <ELocalGlobal lg, ELocalGlobal I, class IterO, class IterN>
01026   void CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::insertIndicesAndValues(RowInfo rowinfo, const SLocalGlobalViews &newInds, IterO rowVals, IterN newVals)
01027   {
01028     size_t numNewInds = insertIndices<lg,I>(rowinfo,newInds);
01029     std::copy( newVals, newVals + numNewInds, rowVals + rowinfo.numEntries );
01030   }
01031 
01032 
01035   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
01036   template <ELocalGlobal lg, class IterO, class IterN, class BinaryFunction>
01037   void CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::transformValues(RowInfo rowinfo, const SLocalGlobalViews &inds, IterO rowVals, IterN newVals, BinaryFunction f) const
01038   {
01039     Teuchos::CompileTimeAssert<lg != GlobalIndices && lg != LocalIndices> cta_lg; (void)cta_lg;
01040     const size_t STINV = OrdinalTraits<size_t>::invalid();
01041     if (lg == GlobalIndices) {
01042       ArrayView<const GlobalOrdinal> search_ginds = inds.ginds;
01043       for (size_t j=0; j < (size_t)search_ginds.size(); ++j) {
01044         const size_t k = findGlobalIndex(rowinfo, search_ginds[j]);
01045         if (k != STINV) {
01046           rowVals[k] = f( rowVals[k], newVals[j] );
01047         }
01048       }
01049     }
01050     else if (lg == LocalIndices) {
01051       ArrayView<const LocalOrdinal> search_linds = inds.linds;
01052       for (size_t j=0; j < (size_t)search_linds.size(); ++j) {
01053         const size_t k = findLocalIndex(rowinfo, search_linds[j]);
01054         if (k != STINV) {
01055           rowVals[k] = f( rowVals[k], newVals[j] );
01056         }
01057       }
01058     }
01059   }
01060 
01061 
01064   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
01065   void CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::sortRowIndices(RowInfo rowinfo)
01066   {
01067     if (rowinfo.numEntries > 0) {
01068       ArrayView<LocalOrdinal> inds_view = getLocalViewNonConst(rowinfo);
01069       std::sort(inds_view.begin(), inds_view.begin() + rowinfo.numEntries);
01070     }
01071   }
01072 
01073 
01076   // in the future, perhaps this could use std::sort with a boost::zip_iterator
01077   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
01078   template <class Scalar>
01079   void CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::sortRowIndicesAndValues(RowInfo rowinfo, ArrayView<Scalar> values)
01080   {
01081     if (rowinfo.numEntries > 0) {
01082       ArrayView<LocalOrdinal> inds_view = getLocalViewNonConst(rowinfo);
01083       sort2(inds_view.begin(), inds_view.begin()+rowinfo.numEntries, values.begin());
01084     }
01085   }
01086 
01087 
01090   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
01091   void CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::mergeRowIndices(RowInfo rowinfo)
01092   {
01093     ArrayView<LocalOrdinal> inds_view = getLocalViewNonConst(rowinfo);
01094     typename ArrayView<LocalOrdinal>::iterator beg, end, newend;
01095     beg = inds_view.begin();
01096     end = inds_view.begin() + rowinfo.numEntries;
01097     newend = std::unique(beg,end);
01098     const size_t mergedEntries = newend - beg;
01099 #ifdef HAVE_TPETRA_DEBUG
01100     // merge should not have eliminated any entries; if so, the assignment below will destory the packed structure
01101     TEUCHOS_TEST_FOR_EXCEPT( isStorageOptimized() && mergedEntries != rowinfo.numEntries );
01102 #endif
01103     numRowEntries_[rowinfo.localRow] = mergedEntries;
01104     nodeNumEntries_ -= (rowinfo.numEntries - mergedEntries);
01105   }
01106 
01107 
01110   // in the future, this could use std::unique with a boost::zip_iterator
01111   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
01112   template <class Iter, class BinaryFunction>
01113   void CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::
01114   mergeRowIndicesAndValues(RowInfo rowinfo, Iter rowValueIter, BinaryFunction f)
01115   {
01116     ArrayView<LocalOrdinal> inds_view = getLocalViewNonConst(rowinfo);
01117     typename ArrayView<LocalOrdinal>::iterator beg, end, newend;
01118 
01119     // beg,end define a half-exclusive interval over which to iterate.
01120     beg = inds_view.begin();
01121     end = inds_view.begin() + rowinfo.numEntries;
01122     newend = beg;
01123     if (beg != end) {
01124       typename ArrayView<LocalOrdinal>::iterator cur = beg + 1;
01125       Iter vcur = rowValueIter + 1,
01126            vend = rowValueIter;
01127       cur = beg+1;
01128       while (cur != end) {
01129         if (*cur != *newend) {
01130           // new entry; save it
01131           ++newend;
01132           ++vend;
01133           (*newend) = (*cur);
01134           (*vend) = (*vcur);
01135         }
01136         else {
01137           // old entry; merge it
01138           (*vend) = f (*vend, *vcur);
01139         }
01140         ++cur;
01141         ++vcur;
01142       }
01143       ++newend; // one past the last entry, per typical [beg,end) semantics
01144     }
01145     const size_t mergedEntries = newend - beg;
01146 #ifdef HAVE_TPETRA_DEBUG
01147     // merge should not have eliminated any entries; if so, the
01148     // assignment below will destroy the packed structure
01149     TEUCHOS_TEST_FOR_EXCEPT( isStorageOptimized() && mergedEntries != rowinfo.numEntries );
01150 #endif // HAVE_TPETRA_DEBUG
01151     numRowEntries_[rowinfo.localRow] = mergedEntries;
01152     nodeNumEntries_ -= (rowinfo.numEntries - mergedEntries);
01153   }
01154 
01155 
01158   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
01159   void CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::setDomainRangeMaps(
01160                                 const RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > &domainMap,
01161                                 const RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > &rangeMap)
01162   {
01163     // simple pointer comparison for equality
01164     if (domainMap_ != domainMap) {
01165       domainMap_ = domainMap;
01166       importer_ = null;
01167     }
01168     if (rangeMap_ != rangeMap) {
01169       rangeMap_  = rangeMap;
01170       exporter_ = null;
01171     }
01172   }
01173 
01174 
01177   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
01178   size_t CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::findLocalIndex(RowInfo rowinfo, LocalOrdinal ind) const
01179   {
01180     typedef typename ArrayView<const LocalOrdinal>::iterator IT;
01181     bool found = true;
01182     // get a view of the row, if it wasn't passed by the caller
01183     ArrayView<const LocalOrdinal> rowinds = getLocalView(rowinfo);
01184     IT rptr, locptr = Teuchos::NullIteratorTraits<IT>::getNull();
01185     rptr = rowinds.begin();
01186     if (isSorted()) {
01187       // binary search
01188       std::pair<IT,IT> p = std::equal_range(rptr,rptr+rowinfo.numEntries,ind);
01189       if (p.first == p.second) found = false;
01190       else locptr = p.first;
01191     }
01192     else {
01193       // direct search
01194       locptr = std::find(rptr,rptr+rowinfo.numEntries,ind);
01195       if (locptr == rptr+rowinfo.numEntries) found = false;
01196     }
01197     size_t ret = OrdinalTraits<size_t>::invalid();
01198     if (found) {
01199       ret = (locptr - rptr);
01200     }
01201     return ret;
01202   }
01203 
01204 
01207   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
01208   size_t CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::findGlobalIndex(RowInfo rowinfo, GlobalOrdinal ind) const
01209   {
01210     typedef typename ArrayView<const GlobalOrdinal>::iterator IT;
01211     bool found = true;
01212     // get a view of the row, if it wasn't passed by the caller
01213     ArrayView<const GlobalOrdinal> rowinds = getGlobalView(rowinfo);
01214     IT rptr, locptr = Teuchos::NullIteratorTraits<IT>::getNull();
01215     rptr = rowinds.begin();
01216     if (isSorted()) {
01217       // binary search
01218       std::pair<IT,IT> p = std::equal_range(rptr,rptr+rowinfo.numEntries,ind);
01219       if (p.first == p.second) found = false;
01220       else locptr = p.first;
01221     }
01222     else {
01223       // direct search
01224       locptr = std::find(rptr,rptr+rowinfo.numEntries,ind);
01225       if (locptr == rptr+rowinfo.numEntries) found = false;
01226     }
01227     size_t ret = OrdinalTraits<size_t>::invalid();
01228     if (found) {
01229       ret = (locptr - rptr);
01230     }
01231     return ret;
01232   }
01233 
01234 
01237   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
01238   void CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::clearGlobalConstants()
01239   {
01240     globalNumEntries_       = OrdinalTraits<global_size_t>::invalid();
01241     globalNumDiags_         = OrdinalTraits<global_size_t>::invalid();
01242     globalMaxNumRowEntries_ = OrdinalTraits<global_size_t>::invalid();
01243     nodeNumDiags_           = OrdinalTraits<       size_t>::invalid();
01244     nodeMaxNumRowEntries_   = OrdinalTraits<       size_t>::invalid();
01245     haveGlobalConstants_    = false;
01246   }
01247 
01248 
01251   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
01252   void CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::checkInternalState() const
01253   {
01254 #ifdef HAVE_TPETRA_DEBUG
01255     const global_size_t GSTI = OrdinalTraits<global_size_t>::invalid();
01256     const size_t         STI = OrdinalTraits<size_t>::invalid();
01257     std::string err = typeName(*this) + "::checkInternalState(): Likely internal logic error. Please contact Tpetra team.";
01258     // check the internal state of this data structure
01259     // this is called by numerous state-changing methods, in a debug build, to ensure that the object
01260     // always remains in a valid state
01261     // the graph should have been allocated with a row map
01262     TEUCHOS_TEST_FOR_EXCEPTION( rowMap_ == null,                     std::logic_error, err );
01263     // am either complete or active
01264     TEUCHOS_TEST_FOR_EXCEPTION( isFillActive() == isFillComplete(),  std::logic_error, err );
01265     // if active, i have no local graph
01266     TEUCHOS_TEST_FOR_EXCEPTION( isFillActive() && lclGraph_ != null, std::logic_error, err );
01267     // if the graph has been fill completed, then all maps should be present
01268     TEUCHOS_TEST_FOR_EXCEPTION( isFillComplete() == true && (colMap_ == null || rangeMap_ == null || domainMap_ == null), std::logic_error, err );
01269     // if storage has been optimized, then indices should have been allocated (even if trivially so)
01270     TEUCHOS_TEST_FOR_EXCEPTION( isStorageOptimized() == true && indicesAreAllocated() == false, std::logic_error, err );
01271     // if storage has been optimized, then number of allocated is now the number of entries
01272     TEUCHOS_TEST_FOR_EXCEPTION( isStorageOptimized() == true && getNodeAllocationSize() != getNodeNumEntries(), std::logic_error, err );
01273     // if graph doesn't have the global constants, then they should all be marked as invalid
01274     TEUCHOS_TEST_FOR_EXCEPTION( haveGlobalConstants_ == false && ( globalNumEntries_ != GSTI || globalNumDiags_ != GSTI || globalMaxNumRowEntries_ != GSTI ), std::logic_error, err );
01275     // if the graph has global cosntants, then they should be valid.
01276     TEUCHOS_TEST_FOR_EXCEPTION( haveGlobalConstants_ == true && ( globalNumEntries_ == GSTI || globalNumDiags_ == GSTI || globalMaxNumRowEntries_ == GSTI ), std::logic_error, err );
01277     TEUCHOS_TEST_FOR_EXCEPTION( haveGlobalConstants_ == true && ( globalNumEntries_ < nodeNumEntries_ || globalNumDiags_ < nodeNumDiags_ || globalMaxNumRowEntries_ < nodeMaxNumRowEntries_ ),
01278                         std::logic_error, err );
01279     // if indices are allocated, then the allocation specifications should have been released
01280     TEUCHOS_TEST_FOR_EXCEPTION( indicesAreAllocated() == true  && (numAllocForAllRows_ != 0 || numAllocPerRow_ != null),                        std::logic_error, err );
01281     // if indices are not allocated, then information dictating allocation quantities should be present
01282     TEUCHOS_TEST_FOR_EXCEPTION( indicesAreAllocated() == false && (nodeNumAllocated_ != STI || nodeNumEntries_ != 0),                           std::logic_error, err );
01283     // if storage is optimized, then profile should be static
01284     TEUCHOS_TEST_FOR_EXCEPTION( isStorageOptimized() && pftype_ != StaticProfile,                                                               std::logic_error, err );
01285     // if rowPtrs_ exists, it is required to have N+1 rows, and rowPtrs_[N] == gblInds1D_.size()/lclInds1D_.size()
01286     TEUCHOS_TEST_FOR_EXCEPTION( isGloballyIndexed() && rowPtrs_ != null && ((size_t)rowPtrs_.size() != getNodeNumRows()+1 || rowPtrs_[getNodeNumRows()] != (size_t)gblInds1D_.size()), std::logic_error, err );
01287     TEUCHOS_TEST_FOR_EXCEPTION(  isLocallyIndexed() && rowPtrs_ != null && ((size_t)rowPtrs_.size() != getNodeNumRows()+1 || rowPtrs_[getNodeNumRows()] != (size_t)lclInds1D_.size()), std::logic_error, err );
01288     // if profile is dynamic and we have allocated, then 2D allocations should be present
01289     TEUCHOS_TEST_FOR_EXCEPTION( pftype_ == DynamicProfile && indicesAreAllocated() && getNodeNumRows() > 0 && lclInds2D_ == null && gblInds2D_ == null,
01290                                                                                                                                         std::logic_error, err );
01291     // if profile is dynamic, then numentries and 2D indices are needed and should be present
01292     TEUCHOS_TEST_FOR_EXCEPTION( pftype_ == DynamicProfile && indicesAreAllocated() && getNodeNumRows() > 0 && (numRowEntries_ == null || (lclInds2D_ == null && gblInds2D_ == null)),
01293                                                                                                                                         std::logic_error, err );
01294     // if profile is dynamic, then 1D allocations should not be present
01295     TEUCHOS_TEST_FOR_EXCEPTION( pftype_ == DynamicProfile && (lclInds1D_ != null || gblInds1D_ != null),                                std::logic_error, err );
01296     // if profile is dynamic, then row offsets should not be present
01297     TEUCHOS_TEST_FOR_EXCEPTION( pftype_ == DynamicProfile && rowPtrs_ != null,                                                          std::logic_error, err );
01298     // if profile is static and we have allocated non-trivially, then 1D allocations should be present
01299     TEUCHOS_TEST_FOR_EXCEPTION( pftype_ == StaticProfile && indicesAreAllocated() && getNodeAllocationSize() > 0 && lclInds1D_ == null && gblInds1D_ == null,
01300                                                                                                                                         std::logic_error, err );
01301     // if profile is static, then 2D allocations should not be present
01302     TEUCHOS_TEST_FOR_EXCEPTION( pftype_ == StaticProfile && (lclInds2D_ != null || gblInds2D_ != null),                                 std::logic_error, err );
01303     // if indices are not allocated, then none of the buffers should be.
01304     TEUCHOS_TEST_FOR_EXCEPTION( indicesAreAllocated() == false && (rowPtrs_ != null || numRowEntries_ != null ||
01305                                                                  lclInds1D_ != null || lclInds2D_ != null ||
01306                                                                  gblInds1D_ != null || gblInds2D_ != null),                             std::logic_error, err );
01307     // indices may be local or global only if they are allocated (numAllocated is redundant; could simply be indicesAreLocal_ || indicesAreGlobal_)
01308     TEUCHOS_TEST_FOR_EXCEPTION( (indicesAreLocal_ == true || indicesAreGlobal_ == true) && indicesAreAllocated_ == false,               std::logic_error, err );
01309     // indices may be local or global, but not both
01310     TEUCHOS_TEST_FOR_EXCEPTION( indicesAreLocal_ == true && indicesAreGlobal_ == true,                                                  std::logic_error, err );
01311     // if indices are local, then global allocations should not be present
01312     TEUCHOS_TEST_FOR_EXCEPTION( indicesAreLocal_ == true && (gblInds1D_ != null || gblInds2D_ != null),                                 std::logic_error, err );
01313     // if indices are global, then local allocations should not be present
01314     TEUCHOS_TEST_FOR_EXCEPTION( indicesAreGlobal_ == true && (lclInds1D_ != null || lclInds2D_ != null),                                std::logic_error, err );
01315     // if indices are local, then local allocations should be present
01316     TEUCHOS_TEST_FOR_EXCEPTION( indicesAreLocal_ == true && getNodeAllocationSize() > 0 && lclInds1D_ == null && getNodeNumRows() > 0 && lclInds2D_ == null,
01317                                                                                                                           std::logic_error, err );
01318     // if indices are global, then global allocations should be present
01319     TEUCHOS_TEST_FOR_EXCEPTION( indicesAreGlobal_ == true && getNodeAllocationSize() > 0 && gblInds1D_ == null && getNodeNumRows() > 0 && gblInds2D_ == null,
01320                                                                                                                           std::logic_error, err );
01321     // if indices are allocated, then we should have recorded how many were allocated
01322     TEUCHOS_TEST_FOR_EXCEPTION( indicesAreAllocated() == true  && nodeNumAllocated_ == STI,                                       std::logic_error, err );
01323     // if indices are not allocated, then the allocation size should be marked invalid
01324     TEUCHOS_TEST_FOR_EXCEPTION( indicesAreAllocated() == false && nodeNumAllocated_ != STI,                                       std::logic_error, err );
01325     // check the actual allocations
01326     if (indicesAreAllocated()) {
01327       size_t actualNumAllocated = 0;
01328       if (pftype_ == DynamicProfile) {
01329         if (isGloballyIndexed() && gblInds2D_ != null) {
01330           for (size_t r = 0; r < getNodeNumRows(); ++r) {
01331             actualNumAllocated += gblInds2D_[r].size();
01332           }
01333         }
01334         else if (isLocallyIndexed() && lclInds2D_ != null) {
01335           for (size_t r = 0; r < getNodeNumRows(); ++r) {
01336             actualNumAllocated += lclInds2D_[r].size();
01337           }
01338         }
01339         TEUCHOS_TEST_FOR_EXCEPTION(actualNumAllocated != nodeNumAllocated_, std::logic_error, err );
01340       }
01341       else if (rowPtrs_ != null) { // pftype_ == StaticProfile
01342         actualNumAllocated = rowPtrs_[getNodeNumRows()];
01343         TEUCHOS_TEST_FOR_EXCEPTION(  isLocallyIndexed() == true && (size_t)lclInds1D_.size() != actualNumAllocated, std::logic_error, err );
01344         TEUCHOS_TEST_FOR_EXCEPTION( isGloballyIndexed() == true && (size_t)gblInds1D_.size() != actualNumAllocated, std::logic_error, err );
01345         TEUCHOS_TEST_FOR_EXCEPTION(actualNumAllocated != nodeNumAllocated_, std::logic_error, err );
01346       }
01347     }
01348 #endif
01349   }
01350 
01351 
01354   //                                                                         //
01355   //                  User-visible class methods                             //
01356   //                                                                         //
01359 
01360 
01363   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
01364   size_t CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getNumEntriesInGlobalRow(GlobalOrdinal globalRow) const
01365   {
01366     const LocalOrdinal lrow = rowMap_->getLocalElement(globalRow);
01367     size_t ret = OrdinalTraits<size_t>::invalid();
01368     if (hasRowInfo() && lrow != OrdinalTraits<LocalOrdinal>::invalid())
01369     {
01370       RowInfo rowinfo = getRowInfo(lrow);
01371       ret = rowinfo.numEntries;
01372     }
01373     return ret;
01374   }
01375 
01376 
01379   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
01380   size_t CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getNumEntriesInLocalRow(LocalOrdinal localRow) const
01381   {
01382     size_t ret = OrdinalTraits<size_t>::invalid();
01383     if (hasRowInfo() && rowMap_->isNodeLocalElement(localRow)) {
01384       RowInfo rowinfo = getRowInfo(localRow);
01385       ret = rowinfo.numEntries;
01386     }
01387     return ret;
01388   }
01389 
01390 
01393   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
01394   size_t CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getNumAllocatedEntriesInGlobalRow(GlobalOrdinal globalRow) const
01395   {
01396     const LocalOrdinal lrow = rowMap_->getLocalElement(globalRow);
01397     size_t ret = OrdinalTraits<size_t>::invalid();
01398     if (hasRowInfo() && lrow != OrdinalTraits<LocalOrdinal>::invalid())
01399     {
01400       RowInfo rowinfo = getRowInfo(lrow);
01401       ret = rowinfo.allocSize;
01402     }
01403     return ret;
01404   }
01405 
01406 
01409   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
01410   size_t CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getNumAllocatedEntriesInLocalRow(LocalOrdinal localRow) const
01411   {
01412     size_t ret = OrdinalTraits<size_t>::invalid();
01413     if (hasRowInfo() && rowMap_->isNodeLocalElement(localRow)) {
01414       RowInfo rowinfo = getRowInfo(localRow);
01415       ret = rowinfo.allocSize;
01416     }
01417     return ret;
01418   }
01419 
01420 
01423   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
01424   ArrayRCP<const size_t> CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getNodeRowPtrs() const
01425   {
01426     return rowPtrs_;
01427   }
01428 
01429 
01432   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
01433   ArrayRCP<const LocalOrdinal> CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getNodePackedIndices() const
01434   {
01435     return lclInds1D_;
01436   }
01437 
01438 
01441   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
01442   void CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getLocalRowCopy(LocalOrdinal localRow, const ArrayView<LocalOrdinal> &indices, size_t &NumIndices) const
01443   {
01444     // can only do this if
01445     // * we have local indices: isLocallyIndexed()
01446     // or
01447     // * we are capable of producing them: isGloballyIndexed() && hasColMap()
01448     // short circuit if we aren't allocated
01449     std::string tfecfFuncName("getLocalRowCopy(localRow,...)");
01450     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(isGloballyIndexed() == true && hasColMap() == false, std::runtime_error, ": local indices cannot be produced.");
01451     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(rowMap_->isNodeLocalElement(localRow) == false, std::runtime_error,
01452         ": localRow (== " << localRow << ") is not valid on this node.");
01453     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(hasRowInfo() == false, std::runtime_error, ": graph row information was deleted at fillComplete().");
01454     const RowInfo rowinfo = getRowInfo(localRow);
01455     NumIndices = rowinfo.numEntries;
01456     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC((size_t)indices.size() < NumIndices, std::runtime_error,
01457         ": specified storage (size==" << indices.size()
01458         << ") is not large enough to hold all entries for this row (NumIndices == " << NumIndices << ").");
01459     if (isLocallyIndexed()) {
01460       ArrayView<const LocalOrdinal> lview = getLocalView(rowinfo);
01461       std::copy( lview.begin(), lview.begin() + NumIndices, indices.begin());
01462     }
01463     else if (isGloballyIndexed()) {
01464       ArrayView<const GlobalOrdinal> gview = getGlobalView(rowinfo);
01465       for (size_t j=0; j < NumIndices; ++j) {
01466         indices[j] = colMap_->getLocalElement(gview[j]);
01467       }
01468     }
01469     else {
01470 #ifdef HAVE_TPETRA_DEBUG
01471       // should have fallen in one of the above
01472       TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC( indicesAreAllocated() == true, std::logic_error, ": Internal logic error. Please contact Tpetra team.");
01473 #endif
01474       NumIndices = 0;
01475     }
01476     return;
01477   }
01478 
01479 
01482   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
01483   void CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getGlobalRowCopy(GlobalOrdinal globalRow, const ArrayView<GlobalOrdinal> &indices, size_t &NumIndices) const
01484   {
01485     // we either currently store global indices, or we have a column map with which to transcribe our local indices for the user
01486     const LocalOrdinal lrow = rowMap_->getLocalElement(globalRow);
01487     std::string tfecfFuncName("getGlobalRowCopy(globalRow,...)");
01488     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(lrow == OrdinalTraits<LocalOrdinal>::invalid(), std::runtime_error,
01489         ": globalRow (== " << globalRow << ") does not belong to this node.");
01490     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(hasRowInfo() == false, std::runtime_error, ": graph row information was deleted at fillComplete().");
01491     const RowInfo rowinfo = getRowInfo((size_t)lrow);
01492     NumIndices = rowinfo.numEntries;
01493     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC((size_t)indices.size() < NumIndices, std::runtime_error,
01494         ": specified storage (size==" << indices.size()
01495         << ") is not large enough to hold all entries for this row (NumIndices == " << NumIndices << ").");
01496     if (isLocallyIndexed()) {
01497       ArrayView<const LocalOrdinal> lview = getLocalView(rowinfo);
01498       for (size_t j=0; j < NumIndices; ++j) {
01499         indices[j] = colMap_->getGlobalElement(lview[j]);
01500       }
01501     }
01502     else if (isGloballyIndexed()) {
01503       ArrayView<const GlobalOrdinal> gview = getGlobalView(rowinfo);
01504       std::copy(gview.begin(), gview.begin() + NumIndices, indices.begin());
01505     }
01506     return;
01507   }
01508 
01509 
01512   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
01513   void CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getLocalRowView(LocalOrdinal localRow, ArrayView<const LocalOrdinal> &indices) const
01514   {
01515     std::string tfecfFuncName("getLocalRowView()");
01516     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(isGloballyIndexed() == true, std::runtime_error, ": local indices cannot be provided.");
01517     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(hasRowInfo() == false, std::runtime_error, ": graph row information was deleted at fillComplete().");
01518     indices = null;
01519     if (rowMap_->isNodeLocalElement(localRow) == true) {
01520       const RowInfo rowinfo = getRowInfo(localRow);
01521       if (rowinfo.numEntries > 0) {
01522         indices = getLocalView(rowinfo);
01523         indices = indices(0,rowinfo.numEntries);
01524       }
01525     }
01526 #ifdef HAVE_TPETRA_DEBUG
01527     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC( (size_t)indices.size() != getNumEntriesInLocalRow(localRow), std::logic_error, ": Violated stated post-conditions. Please contact Tpetra team.");
01528 #endif
01529     return;
01530   }
01531 
01532 
01535   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
01536   void CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getGlobalRowView(GlobalOrdinal globalRow, ArrayView<const GlobalOrdinal> &indices) const
01537   {
01538     std::string tfecfFuncName("getGlobalRowView()");
01539     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(isLocallyIndexed() == true, std::runtime_error, ": global indices cannot be provided.");
01540     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(hasRowInfo() == false, std::runtime_error, ": graph row information was deleted at fillComplete().");
01541     indices = null;
01542     if (rowMap_->isNodeGlobalElement(globalRow) == true) {
01543       const LocalOrdinal lrow = rowMap_->getLocalElement(globalRow);
01544       const RowInfo rowinfo = getRowInfo(lrow);
01545       if (rowinfo.numEntries > 0) {
01546         indices = getGlobalView(rowinfo);
01547         indices = indices(0,rowinfo.numEntries);
01548       }
01549     }
01550 #ifdef HAVE_TPETRA_DEBUG
01551     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC( (size_t)indices.size() != getNumEntriesInGlobalRow(globalRow), std::logic_error, ": Violated stated post-conditions. Please contact Tpetra team.");
01552 #endif
01553     return;
01554   }
01555 
01556 
01559   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
01560   void CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::insertLocalIndices(
01561             LocalOrdinal localRow, const ArrayView<const LocalOrdinal> &indices)
01562   {
01563     std::string tfecfFuncName("insertLocalIndices()");
01564     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC( isFillActive() == false,                        std::runtime_error, ": requires that fill is active.");
01565     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC( isGloballyIndexed() == true,                    std::runtime_error, ": graph indices are global; use insertGlobalIndices().");
01566     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC( hasColMap() == false,                           std::runtime_error, ": cannot insert local indices without a column map.");
01567     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC( rowMap_->isNodeLocalElement(localRow) == false, std::runtime_error, ": row does not belong to this node.");
01568     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(hasRowInfo() == false, std::runtime_error, ": graph row information was deleted at fillComplete().");
01569     if (indicesAreAllocated() == false) {
01570       allocateIndices(LocalIndices);
01571     }
01572     // use column map to filter the entries
01573     Array<LocalOrdinal> f_inds(indices);
01574     SLocalGlobalNCViews inds_ncview;
01575     inds_ncview.linds = f_inds();
01576     const size_t numFilteredEntries = filterIndices<LocalIndices>(inds_ncview);
01577     if (numFilteredEntries > 0) {
01578       RowInfo rowInfo = getRowInfo(localRow);
01579       const size_t curNumEntries = rowInfo.numEntries;
01580       const size_t newNumEntries = curNumEntries + numFilteredEntries;
01581       if (newNumEntries > rowInfo.allocSize) {
01582         TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(getProfileType() == StaticProfile, std::runtime_error, ": new indices exceed statically allocated graph structure.");
01583         // Only print an efficiency warning once per CrsGraph
01584         // instance, per method name (insertLocalIndices() or
01585         // insertGlobalIndices()).
01586         if (! insertLocalIndicesWarnedEfficiency_) {
01587           TPETRA_EFFICIENCY_WARNING(true, std::runtime_error,
01588             "::insertLocalIndices():" << std::endl
01589             << "Pre-allocated space has been exceeded, requiring new allocation.  "
01590             "This is allowed but not efficient in terms of run time.  "
01591             "To improve efficiency, suggest a larger number of entries per row in the constructor.  "
01592             "You may either specify a maximum number of entries for all the rows, or a per-row maximum.  "
01593             "This CrsGraph instance will not print further messages of this kind, in order not to clutter output.");
01594           insertLocalIndicesWarnedEfficiency_ = true;
01595         }
01596         // update allocation only as much as necessary
01597         rowInfo = updateAlloc<LocalIndices>(rowInfo, newNumEntries);
01598       }
01599       SLocalGlobalViews inds_view;
01600       inds_view.linds = f_inds(0,numFilteredEntries);
01601       insertIndices<LocalIndices,LocalIndices>(rowInfo, inds_view);
01602     }
01603 #ifdef HAVE_TPETRA_DEBUG
01604     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(indicesAreAllocated() == false || isLocallyIndexed() == false, std::logic_error,
01605         ": Violated stated post-conditions. Please contact Tpetra team.");
01606 #endif
01607   }
01608 
01609 
01612   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
01613   void CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::insertGlobalIndices(
01614             GlobalOrdinal grow, const ArrayView<const GlobalOrdinal> &indices)
01615   {
01616     std::string tfecfFuncName("insertGlobalIndices()");
01617     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(isLocallyIndexed() == true, std::runtime_error, ": graph indices are local; use insertLocalIndices().");
01618     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(hasRowInfo() == false, std::runtime_error, ": graph row information was deleted at fillComplete().");
01619     // this can't really be satisfied for now, because if we are fillComplete(), then we are local
01620     // in the future, this may change. however, the rule that modification require active fill will not change.
01621     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC( isFillActive() == false, std::runtime_error, ": requires that fill is active.");
01622     if (indicesAreAllocated() == false) {
01623       allocateIndices(GlobalIndices);
01624     }
01625     const LocalOrdinal myRow = rowMap_->getLocalElement(grow);
01626     if (myRow != OrdinalTraits<LocalOrdinal>::invalid())
01627     {
01628       // if we have a column map, use it to filter the entries.
01629       Array<GlobalOrdinal> filtered_indices;
01630       SLocalGlobalViews inds_view;
01631       if (hasColMap()) {
01632         SLocalGlobalNCViews inds_ncview;
01633         // filter indices and values through the column map
01634         filtered_indices.assign(indices.begin(), indices.end());
01635         inds_ncview.ginds = filtered_indices();
01636         const size_t numFilteredEntries = filterIndices<GlobalIndices>(inds_ncview);
01637         inds_view.ginds = filtered_indices(0,numFilteredEntries);
01638       }
01639       else {
01640         inds_view.ginds = indices;
01641       }
01642       const size_t numFilteredEntries = inds_view.ginds.size();
01643       // add the new indices and values
01644       if (numFilteredEntries > 0) {
01645         RowInfo rowInfo = getRowInfo(myRow);
01646         const size_t curNumEntries = rowInfo.numEntries;
01647         const size_t newNumEntries = curNumEntries + numFilteredEntries;
01648         if (newNumEntries > rowInfo.allocSize) {
01649           TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(getProfileType() == StaticProfile, std::runtime_error, ": new indices exceed statically allocated graph structure.");
01650           // Only print an efficiency warning once per CrsGraph
01651           // instance, per method name (insertLocalIndices() or
01652           // insertGlobalIndices()).
01653           if (! insertGlobalIndicesWarnedEfficiency_) {
01654             TPETRA_EFFICIENCY_WARNING(true, std::runtime_error,
01655               "::insertGlobalIndices():" << std::endl
01656               << "Pre-allocated space has been exceeded, requiring new allocation.  "
01657               "This is allowed but not efficient in terms of run time.  "
01658               "To improve efficiency, suggest a larger number of entries per row in the constructor.  "
01659               "You may either specify a maximum number of entries for all the rows, or a per-row maximum.  "
01660               "This CrsGraph instance will not print further messages of this kind, in order not to clutter output.");
01661             insertGlobalIndicesWarnedEfficiency_ = true;
01662           }
01663           // update allocation only as much as necessary
01664           rowInfo = updateAlloc<GlobalIndices>(rowInfo, newNumEntries);
01665         }
01666         insertIndices<GlobalIndices,GlobalIndices>(rowInfo, inds_view);
01667 #ifdef HAVE_TPETRA_DEBUG
01668         {
01669           const size_t chkNewNumEntries = getNumEntriesInLocalRow(myRow);
01670           TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(chkNewNumEntries != newNumEntries, std::logic_error, ": Internal logic error. Please contact Tpetra team.");
01671         }
01672 #endif
01673       }
01674     }
01675     else {
01676       // a nonlocal row
01677       for (typename ArrayView<const GlobalOrdinal>::iterator i=indices.begin(); i != indices.end(); ++i) {
01678         nonlocals_[grow].push_back(*i);
01679       }
01680     }
01681 #ifdef HAVE_TPETRA_DEBUG
01682     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(indicesAreAllocated() == false || isGloballyIndexed() == false, std::logic_error,
01683         ": Violated stated post-conditions. Please contact Tpetra team.");
01684 #endif
01685   }
01686 
01687 
01690   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
01691   void CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::removeLocalIndices(LocalOrdinal lrow)
01692   {
01693     std::string tfecfFuncName("removeLocalIndices()");
01694     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC( isFillActive() == false,                    std::runtime_error, ": requires that fill is active.");
01695     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC( isStorageOptimized() == true,               std::runtime_error, ": cannot remove indices after optimizeStorage() has been called.");
01696     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC( isGloballyIndexed() == true,                std::runtime_error, ": graph indices are global.");
01697     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC( rowMap_->isNodeLocalElement(lrow) == false, std::runtime_error, ": row does not belong to this node.");
01698     if (indicesAreAllocated() == false) {
01699       allocateIndices(LocalIndices);
01700     }
01701     //
01702     clearGlobalConstants();
01703     //
01704     if (numRowEntries_ != null) {
01705       nodeNumEntries_ -= numRowEntries_[lrow];
01706       numRowEntries_[lrow] = 0;
01707     }
01708 #ifdef HAVE_TPETRA_DEBUG
01709     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(getNumEntriesInLocalRow(lrow) != 0 || indicesAreAllocated() == false || isLocallyIndexed() == false, std::logic_error,
01710         ": Violated stated post-conditions. Please contact Tpetra team.");
01711 #endif
01712   }
01713 
01714 
01715   // TODO: in the future, globalAssemble() should use import/export functionality
01718   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
01719   void CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::globalAssemble()
01720   {
01721     using Teuchos::as;
01722     using Teuchos::gatherAll;
01723     using Teuchos::ireceive;
01724     using Teuchos::isend;
01725     using Teuchos::outArg;
01726     using Teuchos::REDUCE_MAX;
01727     using Teuchos::reduceAll;
01728     using Teuchos::waitAll;
01729     using std::deque;
01730     using std::pair;
01731     using std::make_pair;
01732     typedef GlobalOrdinal GO;
01733     typedef typename std::map<GO, std::deque<GO> >::const_iterator NLITER;
01734 
01735     std::string tfecfFuncName("globalAssemble()"); // for exception macro
01736     RCP<const Comm<int> > comm = getComm();
01737 
01738     const int numImages = comm->getSize();
01739     const int myImageID = comm->getRank();
01740 #ifdef HAVE_TPETRA_DEBUG
01741     Teuchos::barrier (*comm);
01742 #endif // HAVE_TPETRA_DEBUG
01743 
01744     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC( ! isFillActive(), std::runtime_error,
01745       ": requires that fill is active.");
01746     // Determine if any nodes have global entries to share
01747     {
01748       size_t MyNonlocals = nonlocals_.size(), MaxGlobalNonlocals;
01749       reduceAll (*comm, REDUCE_MAX, MyNonlocals, outArg (MaxGlobalNonlocals));
01750       if (MaxGlobalNonlocals == 0) {
01751         return;  // no entries to share
01752       }
01753     }
01754 
01755     // compute a list of NLRs from nonlocals_ and use it to compute:
01756     //      IdsAndRows: a vector of (id,row) pairs
01757     //          NLR2Id: a map from NLR to the Id that owns it
01758     // globalNeighbors: a global graph of connectivity between images: globalNeighbors(i,j) indicates that j sends to i
01759     //         sendIDs: a list of all images I send to
01760     //         recvIDs: a list of all images I receive from (constructed later)
01761     Array<pair<int, GO> > IdsAndRows;
01762     std::map<GO, int> NLR2Id;
01763     Teuchos::SerialDenseMatrix<int, char> globalNeighbors;
01764     Array<int> sendIDs, recvIDs;
01765     {
01766       // nonlocals_ contains the entries we are holding for all non-local rows
01767       // we want a list of the rows for which we have data
01768       Array<GO> NLRs;
01769       std::set<GO> setOfRows;
01770       for (NLITER iter = nonlocals_.begin(); iter != nonlocals_.end(); ++iter) {
01771         setOfRows.insert(iter->first);
01772       }
01773       // copy the elements in the set into an Array
01774       NLRs.resize(setOfRows.size());
01775       std::copy(setOfRows.begin(), setOfRows.end(), NLRs.begin());
01776 
01777       // get a list of ImageIDs for the non-local rows (NLRs)
01778       Array<int> NLRIds(NLRs.size());
01779       {
01780         LookupStatus stat = rowMap_->getRemoteIndexList(NLRs(),NLRIds());
01781         char lclerror = ( stat == IDNotPresent ? 1 : 0 );
01782         char gblerror;
01783         reduceAll (*getComm(), REDUCE_MAX, lclerror, outArg (gblerror));
01784         TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(gblerror != 0, std::runtime_error,
01785           ": nonlocal entries correspond to invalid rows.");
01786       }
01787 
01788       // build up a list of neighbors, as well as a map between NLRs and Ids
01789       // localNeighbors[i] != 0 iff I have data to send to image i
01790       // put NLRs,Ids into an array of pairs
01791       IdsAndRows.reserve(NLRs.size());
01792       Array<char> localNeighbors(numImages, 0);
01793       typename Array<GO>::const_iterator nlr;
01794       typename Array<int>::const_iterator id;
01795       for (nlr = NLRs.begin(), id = NLRIds.begin();
01796            nlr != NLRs.end(); ++nlr, ++id) {
01797         NLR2Id[*nlr] = *id;
01798         localNeighbors[*id] = 1;
01799         // IdsAndRows.push_back(make_pair<int,GlobalOrdinal>(*id,*nlr));
01800         IdsAndRows.push_back(make_pair(*id,*nlr));
01801       }
01802       for (int j=0; j<numImages; ++j) {
01803         if (localNeighbors[j]) {
01804           sendIDs.push_back(j);
01805         }
01806       }
01807       // sort IdsAndRows, by Ids first, then rows
01808       std::sort(IdsAndRows.begin(),IdsAndRows.end());
01809       // gather from other nodes to form the full graph
01810       globalNeighbors.shapeUninitialized(numImages,numImages);
01811       gatherAll (*getComm(), numImages, localNeighbors.getRawPtr(),
01812                  numImages * numImages, globalNeighbors.values());
01813       // globalNeighbors at this point contains (on all images) the
01814       // connectivity between the images.
01815       // globalNeighbors(i,j) != 0 means that j sends to i/that i receives from j
01816     }
01817 
01819     // FIGURE OUT WHO IS SENDING TO WHOM AND HOW MUCH
01820     // DO THIS IN THE PROCESS OF PACKING ALL OUTGOING DATA ACCORDING TO DESTINATION ID
01822 
01823     // loop over all columns to know from which images I can expect to receive something
01824     for (int j=0; j<numImages; ++j) {
01825       if (globalNeighbors(myImageID,j)) {
01826         recvIDs.push_back(j);
01827       }
01828     }
01829     const size_t numRecvs = recvIDs.size();
01830 
01831     // we know how many we're sending to already
01832     // form a contiguous list of all data to be sent
01833     // track the number of entries for each ID
01834     Array<pair<GO, GO> > IJSendBuffer;
01835     Array<size_t> sendSizes(sendIDs.size(), 0);
01836     size_t numSends = 0;
01837     for (typename Array<pair<int, GO> >::const_iterator IdAndRow = IdsAndRows.begin();
01838          IdAndRow != IdsAndRows.end(); ++IdAndRow) {
01839       int id = IdAndRow->first;
01840       GO row = IdAndRow->second;
01841       // have we advanced to a new send?
01842       if (sendIDs[numSends] != id) {
01843         numSends++;
01844         TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(sendIDs[numSends] != id,
01845           std::logic_error, ": internal logic error. Contact Tpetra team.");
01846       }
01847       // copy data for row into contiguous storage
01848       for (typename deque<GO>::const_iterator j = nonlocals_[row].begin(); j != nonlocals_[row].end(); ++j)
01849       {
01850         IJSendBuffer.push_back( pair<GlobalOrdinal,GlobalOrdinal>(row,*j) );
01851         sendSizes[numSends]++;
01852       }
01853     }
01854     if (IdsAndRows.size() > 0) {
01855       numSends++; // one last increment, to make it a count instead of an index
01856     }
01857     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(as<typename Array<int>::size_type>(numSends) != sendIDs.size(), std::logic_error, ": internal logic error. Contact Tpetra team.");
01858 
01859     // don't need this data anymore
01860     nonlocals_.clear();
01861 
01863     // TRANSMIT SIZE INFO BETWEEN SENDERS AND RECEIVERS
01865 
01866     // Array of pending nonblocking communication requests.  It's OK
01867     // to mix nonblocking send and receive requests in the same
01868     // waitAll() call.
01869     Array<RCP<Teuchos::CommRequest> > requests;
01870 
01871     // perform non-blocking sends: send sizes to our recipients
01872     for (size_t s = 0; s < numSends ; ++s) {
01873       // We're using a nonowning RCP because all communication
01874       // will be local to this method and the scope of our data
01875       requests.push_back (isend<int, size_t> (*comm,
01876                                               rcp (&sendSizes[s], false),
01877                                               sendIDs[s]));
01878     }
01879     // perform non-blocking receives: receive sizes from our senders
01880     Array<size_t> recvSizes (numRecvs);
01881     for (size_t r = 0; r < numRecvs; ++r) {
01882       // We're using a nonowning RCP because all communication
01883       // will be local to this method and the scope of our data
01884       requests.push_back (ireceive (*comm, rcp (&recvSizes[r], false), recvIDs[r]));
01885     }
01886     // Wait on all the nonblocking sends and receives.
01887     if (! requests.empty()) {
01888       waitAll (*comm, requests());
01889     }
01890 #ifdef HAVE_TPETRA_DEBUG
01891     Teuchos::barrier (*comm);
01892 #endif // HAVE_TPETRA_DEBUG
01893 
01894     // This doesn't necessarily deallocate the array.
01895     requests.resize (0);
01896 
01898     // NOW SEND/RECEIVE ALL ROW DATA
01900     // from the size info, build the ArrayViews into IJSendBuffer
01901     Array<ArrayView<pair<GO,GO> > > sendBuffers(numSends,null);
01902     {
01903       size_t cur = 0;
01904       for (size_t s=0; s<numSends; ++s) {
01905         sendBuffers[s] = IJSendBuffer(cur,sendSizes[s]);
01906         cur += sendSizes[s];
01907       }
01908     }
01909     // perform non-blocking sends
01910     for (size_t s=0; s < numSends ; ++s)
01911     {
01912       // We're using a nonowning RCP because all communication
01913       // will be local to this method and the scope of our data
01914       ArrayRCP<pair<GO,GO> > tmpSendBuf =
01915         arcp (sendBuffers[s].getRawPtr(), 0, sendBuffers[s].size(), false);
01916       requests.push_back (isend<int, pair<GO,GO> > (*comm, tmpSendBuf, sendIDs[s]));
01917     }
01918     // calculate amount of storage needed for receives
01919     // setup pointers for the receives as well
01920     size_t totalRecvSize = std::accumulate (recvSizes.begin(), recvSizes.end(), 0);
01921     Array<pair<GO,GO> > IJRecvBuffer (totalRecvSize);
01922     // from the size info, build the ArrayViews into IJRecvBuffer
01923     Array<ArrayView<pair<GO,GO> > > recvBuffers (numRecvs, null);
01924     {
01925       size_t cur = 0;
01926       for (size_t r=0; r<numRecvs; ++r) {
01927         recvBuffers[r] = IJRecvBuffer(cur,recvSizes[r]);
01928         cur += recvSizes[r];
01929       }
01930     }
01931     // perform non-blocking recvs
01932     for (size_t r = 0; r < numRecvs; ++r) {
01933       // We're using a nonowning RCP because all communication
01934       // will be local to this method and the scope of our data
01935       ArrayRCP<pair<GO,GO> > tmpRecvBuf =
01936         arcp (recvBuffers[r].getRawPtr(), 0, recvBuffers[r].size(), false);
01937       requests.push_back (ireceive (*comm, tmpRecvBuf, recvIDs[r]));
01938     }
01939     // perform waits
01940     if (! requests.empty()) {
01941       waitAll (*comm, requests());
01942     }
01943 #ifdef HAVE_TPETRA_DEBUG
01944     Teuchos::barrier (*comm);
01945 #endif // HAVE_TPETRA_DEBUG
01946 
01948     // NOW PROCESS THE RECEIVED ROW DATA
01950     // TODO: instead of adding one entry at a time, add one row at a time.
01951     //       this requires resorting; they arrived sorted by sending node, so that entries could be non-contiguous if we received
01952     //       multiple entries for a particular row from different processors.
01953     //       it also requires restoring the data, which may make it not worth the trouble.
01954     for (typename Array<pair<GO,GO> >::const_iterator ij = IJRecvBuffer.begin();
01955          ij != IJRecvBuffer.end(); ++ij)
01956     {
01957       insertGlobalIndices(ij->first, tuple<GO> (ij->second));
01958     }
01959     checkInternalState();
01960   }
01961 
01962 
01965   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
01966   void
01967   CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::
01968   resumeFill (const RCP<ParameterList> &params)
01969   {
01970     std::string tfecfFuncName("resumeFill");
01971     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(! hasRowInfo(), std::runtime_error,
01972       ": Sorry, you cannot resume fill of the CrsGraph, since the graph's row "
01973       "information was deleted in fillComplete().");
01974 
01975 #ifdef HAVE_TPETRA_DEBUG
01976     Teuchos::barrier( *rowMap_->getComm() );
01977 #endif // HAVE_TPETRA_DEBUG
01978     clearGlobalConstants();
01979     lclGraph_ = null;
01980     if (params != null) this->setParameterList (params);
01981     setSorted(true);
01982     lowerTriangular_  = false;
01983     upperTriangular_  = false;
01984     noRedundancies_   = false;
01985     fillComplete_ = false;
01986 #ifdef HAVE_TPETRA_DEBUG
01987     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
01988       ! isFillActive() || isFillComplete(), std::logic_error,
01989       "::resumeFill(): At end of method, either fill is not active or fill is "
01990       "complete.  This violates stated post-conditions.  Please report this bug "
01991       "to the Tpetra developers.");
01992 #endif // HAVE_TPETRA_DEBUG
01993   }
01994 
01995 
01998   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
01999   void CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::fillComplete(const RCP<ParameterList> &params)
02000   {
02001     fillComplete(rowMap_,rowMap_,params);
02002   }
02003 
02004 
02007   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
02008   void CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::fillComplete(
02009                                     const RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > &domainMap,
02010                                     const RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > &rangeMap,
02011                                     const RCP<ParameterList> &params)
02012   {
02013 #ifdef HAVE_TPETRA_DEBUG
02014     Teuchos::barrier( *rowMap_->getComm() );
02015 #endif // HAVE_TPETRA_DEBUG
02016     std::string tfecfFuncName("fillComplete()");
02017     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC( ! isFillActive() || isFillComplete(),
02018       std::runtime_error, ": Graph fill state must be active (isFillActive() "
02019       "must be true) before calling fillComplete().");
02020 
02021     // allocate if unallocated
02022     if (! indicesAreAllocated()) {
02023       // allocate global, in case we do not have a column map
02024       allocateIndices( GlobalIndices );
02025     }
02026     // Global assemble, if we need to (we certainly don't need to if
02027     // there's only one process).  This call only costs a single
02028     // all-reduce if we don't need global assembly.
02029     if (getComm()->getSize() > 1) {
02030       globalAssemble();
02031     }
02032     else {
02033       TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
02034           nonlocals_.size() > 0,
02035           std::runtime_error, ": cannot have non-local entries on a serial run. Invalid entries were submitted to the CrsMatrix.");
02036     }
02037     // set domain/range map: may clear the import/export objects
02038     setDomainRangeMaps(domainMap,rangeMap);
02039     // make column map
02040     if (! hasColMap()) {
02041       makeColMap();
02042     }
02043     if (isGloballyIndexed()) {
02044       makeIndicesLocal();
02045     }
02046     if (! isSorted()) {
02047       sortAllIndices();
02048     }
02049     if (! isMerged()) {
02050       mergeAllIndices();
02051     }
02052     makeImportExport(); // Make Import and Export objects
02053     computeGlobalConstants();
02054     // fill local objects
02055     fillLocalGraph(params);
02056     //
02057     fillComplete_ = true;
02058 #ifdef HAVE_TPETRA_DEBUG
02059     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC( isFillActive() == true || isFillComplete() == false, std::logic_error, ": Violated stated post-conditions. Please contact Tpetra team.");
02060 #endif
02061     //
02062     checkInternalState();
02063   }
02064 
02065 
02068   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
02069   void CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::fillLocalGraph(const RCP<ParameterList> &params)
02070   {
02071     const size_t numRows = getNodeNumRows();
02072     ArrayRCP<size_t> ptrs;
02073     ArrayRCP<LocalOrdinal> inds;
02074     bool requestOptimizedStorage = true;
02075     if (params != null && params->get("Optimize Storage",true) == false) requestOptimizedStorage = false;
02076     if (getProfileType() == DynamicProfile) {
02077       // 2d -> 1d packed
02078       ptrs = LocalMatOps::allocRowPtrs( getRowMap()->getNode(), numRowEntries_() );
02079       inds = LocalMatOps::template allocStorage<LocalOrdinal>( getRowMap()->getNode(), ptrs() );
02080       for (size_t row=0; row < numRows; ++row) {
02081         const size_t numentrs = numRowEntries_[row];
02082         std::copy( lclInds2D_[row].begin(), lclInds2D_[row].begin()+numentrs, inds+ptrs[row] );
02083       }
02084     }
02085     else if (getProfileType() == StaticProfile) {
02086       // 1d non-packed -> 1d packed
02087       if (nodeNumEntries_ != nodeNumAllocated_) {
02088         ptrs = LocalMatOps::allocRowPtrs( getRowMap()->getNode(), numRowEntries_() );
02089         inds = LocalMatOps::template allocStorage<LocalOrdinal>( getRowMap()->getNode(), ptrs() );
02090         for (size_t row=0; row < numRows; ++row) {
02091           const size_t numentrs = numRowEntries_[row];
02092           std::copy( lclInds1D_+rowPtrs_[row], lclInds1D_+rowPtrs_[row]+numentrs, inds+ptrs[row] );
02093         }
02094       }
02095       else {
02096         inds = lclInds1D_;
02097         ptrs = rowPtrs_;
02098       }
02099     }
02100     // can we ditch the old allocations for the packed one?
02101     if ( requestOptimizedStorage ) {
02102       lclInds2D_ = null;
02103       numRowEntries_ = null;
02104       // keep the new stuff
02105       lclInds1D_ = inds;
02106       rowPtrs_ = ptrs;
02107       nodeNumAllocated_ = nodeNumEntries_;
02108       pftype_ = StaticProfile;
02109     }
02110     // build the local graph, hand over the indices
02111     RCP<ParameterList> lclparams;
02112     if (params == null) lclparams = parameterList();
02113     else                lclparams = sublist(params,"Local Graph");
02114     lclGraph_ = rcp( new local_graph_type( getRowMap()->getNodeNumElements(), getColMap()->getNodeNumElements(), getRowMap()->getNode(), lclparams ) );
02115     lclGraph_->setStructure(ptrs,inds);
02116     ptrs = null;
02117     inds = null;
02118     // finalize local graph
02119     Teuchos::EDiag diag = ( getNodeNumDiags() < getNodeNumRows() ? Teuchos::UNIT_DIAG : Teuchos::NON_UNIT_DIAG );
02120     Teuchos::EUplo uplo = Teuchos::UNDEF_TRI;
02121     if      (isUpperTriangular()) uplo = Teuchos::UPPER_TRI;
02122     else if (isLowerTriangular()) uplo = Teuchos::LOWER_TRI;
02123     LocalMatOps::finalizeGraph(uplo,diag,*lclGraph_,params);
02124   }
02125 
02126 
02129   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
02130   const RCP<const typename LocalMatOps::template graph<LocalOrdinal,Node>::graph_type>
02131   CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getLocalGraph() const
02132   {
02133     return lclGraph_;
02134   }
02135 
02136 
02139   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
02140   const RCP<typename LocalMatOps::template graph<LocalOrdinal,Node>::graph_type>
02141   CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getLocalGraphNonConst()
02142   {
02143     return lclGraph_;
02144   }
02145 
02146 
02149   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
02150   void CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::computeGlobalConstants()
02151   {
02152     // compute the local constants first
02153     const size_t nlrs = getNodeNumRows();
02154     // reset all local properties
02155     upperTriangular_ = true;
02156     lowerTriangular_ = true;
02157     nodeMaxNumRowEntries_ = 0;
02158     nodeNumDiags_         = 0;
02159     // indices are already sorted in each row
02160     const Map<LocalOrdinal,GlobalOrdinal,Node> &rowMap = *rowMap_;
02161     if (indicesAreAllocated() == true && nodeNumAllocated_ > 0) {
02162       for (size_t r=0; r < nlrs; ++r) {
02163         GlobalOrdinal rgid = rowMap.getGlobalElement(r);
02164         // determine the local column index for this row, used for delimiting the diagonal
02165         const LocalOrdinal rlcid = colMap_->getLocalElement(rgid);
02166         RowInfo rowinfo = getRowInfo(r);
02167         ArrayView<const LocalOrdinal> rview = getLocalView(rowinfo);
02168         typename ArrayRCP<const LocalOrdinal>::iterator beg, end, cur;
02169         beg = rview.begin();
02170         end = beg + rowinfo.numEntries;
02171         if (beg != end) {
02172           for (cur = beg; cur != end; ++cur) {
02173             // is this the diagonal?
02174             if (rlcid == (*cur)) ++nodeNumDiags_;
02175           }
02176           // because of sorting, smallest column index is (*beg); it indicates upper triangularity
02177           if (Teuchos::as<size_t>(beg[0]) < r) upperTriangular_ = false;
02178           // because of sorting, largest column index is (*newend); it indicates lower triangularity
02179           if (r < Teuchos::as<size_t>(end[-1])) lowerTriangular_ = false;
02180         }
02181         // compute num entries for this row, accumulate into nodeNumEntries_, update nodeMaxNumRowEntries_
02182         nodeMaxNumRowEntries_ = std::max( nodeMaxNumRowEntries_, rowinfo.numEntries );
02183       }
02184     }
02185 
02186     // compute global constants using computed local constants
02187     if (haveGlobalConstants_ == false) {
02188       global_size_t lcl[2], gbl[2];
02189       lcl[0] = nodeNumEntries_;
02190       lcl[1] = nodeNumDiags_;
02191       Teuchos::reduceAll<int,global_size_t>(*getComm(),Teuchos::REDUCE_SUM,2,lcl,gbl);
02192       globalNumEntries_ = gbl[0];
02193       globalNumDiags_   = gbl[1];
02194       Teuchos::reduceAll<int,global_size_t>(*getComm(),Teuchos::REDUCE_MAX,nodeMaxNumRowEntries_,
02195           outArg(globalMaxNumRowEntries_));
02196       haveGlobalConstants_ = true;
02197     }
02198   }
02199 
02200 
02203   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
02204   void CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::makeIndicesLocal()
02205   {
02206     // All nodes must be in the same index state.
02207     // Update index state by checking isLocallyIndexed/Global on all nodes
02208     computeIndexState();
02209     std::string tfecfFuncName("makeIndicesLocal()");
02210     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(isLocallyIndexed() && isGloballyIndexed(), std::logic_error,
02211         ": inconsistent index state. Indices must be local on all nodes or global on all nodes.");
02212     // If user has not prescribed column map, create one from indices
02213     makeColMap();
02214     // Transform indices to local index space
02215     const size_t nlrs = getNodeNumRows();
02216     //
02217     if (isGloballyIndexed() && nlrs > 0) {
02218       // allocate data for local indices
02219       if (getProfileType() == StaticProfile) {
02220         // reinterpret the the compute buffer as LocalOrdinal if they are the same size
02221         // otherwise, just reallocate
02222         if (nodeNumAllocated_ && sizeof(LocalOrdinal) == sizeof(GlobalOrdinal) ) {
02223           lclInds1D_ = arcp_reinterpret_cast<LocalOrdinal>(gblInds1D_).persistingView(0,nodeNumAllocated_);
02224         }
02225         else {
02226           lclInds1D_ = LocalMatOps::template allocStorage<LocalOrdinal>( getRowMap()->getNode(), rowPtrs_() );
02227         }
02228         for (size_t r=0; r < getNodeNumRows(); ++r) {
02229           const size_t offset   = rowPtrs_[r],
02230                        numentry = numRowEntries_[r];
02231           for (size_t j=0; j<numentry; ++j) {
02232             GlobalOrdinal gid = gblInds1D_[offset + j];
02233             LocalOrdinal  lid = colMap_->getLocalElement(gid);
02234             lclInds1D_[offset + j] = lid;
02235 #ifdef HAVE_TPETRA_DEBUG
02236             TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(lclInds1D_[offset + j] == OrdinalTraits<LocalOrdinal>::invalid(), std::logic_error,
02237                 ": Internal error in fillComplete(). Please contact Tpetra team.");
02238 #endif
02239           }
02240         }
02241         // done with this pointer (allocation will persist in lclInds1D_)
02242         gblInds1D_ = null;
02243       }
02244       else {  // getProfileType() == DynamicProfile
02245         lclInds2D_ = arcp< ArrayRCP<LocalOrdinal> >(nlrs);
02246         // if we have views, then make views
02247         for (size_t r=0; r < getNodeNumRows(); ++r) {
02248           if (gblInds2D_[r] != null) {
02249             ArrayRCP<GlobalOrdinal> ginds = gblInds2D_[r];
02250             const size_t rna = gblInds2D_[r].size();
02251             ArrayRCP< LocalOrdinal> linds = arcp_reinterpret_cast<LocalOrdinal>(ginds).persistingView(0,rna);
02252             // do the conversion in situ. this must be done from front to back.
02253             const size_t numentry = numRowEntries_[r];
02254             for (size_t j=0; j < numentry; ++j) {
02255               GlobalOrdinal gid = ginds[j];
02256               LocalOrdinal  lid = colMap_->getLocalElement(gid);
02257               linds[j] = lid;
02258 #ifdef HAVE_TPETRA_DEBUG
02259               TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(linds[j] == OrdinalTraits<LocalOrdinal>::invalid(), std::logic_error,
02260                   ": Internal error in makeIndicesLocal(). Please contact Tpetra team.");
02261 #endif
02262             }
02263             // reinterpret the the compute buffer as LocalOrdinal; keep the original size
02264             lclInds2D_[r] = linds;
02265             gblInds2D_[r] = null;
02266           }
02267         }
02268         gblInds2D_ = null;
02269       }
02270     }
02271     indicesAreLocal_  = true;
02272     indicesAreGlobal_ = false;
02273     checkInternalState();
02274   }
02275 
02276 
02279   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
02280   void CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::computeIndexState()
02281   {
02282     char myIndices[2] = {0,0};
02283     if (indicesAreLocal_)  myIndices[0] = 1;
02284     if (indicesAreGlobal_) myIndices[1] = 1;
02285     char allIndices[2];
02286     Teuchos::reduceAll(*getComm(),Teuchos::REDUCE_MAX,2,myIndices,allIndices);
02287     indicesAreLocal_  = (allIndices[0]==1);  // If indices are local on one PE, should be local on all
02288     indicesAreGlobal_ = (allIndices[1]==1);  // If indices are global on one PE should be local on all
02289   }
02290 
02291 
02294   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
02295   void CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::sortAllIndices()
02296   {
02297     TEUCHOS_TEST_FOR_EXCEPT(isGloballyIndexed()==true);   // this should be called only after makeIndicesLocal()
02298     if (isSorted() == false) {
02299       for (size_t row=0; row < getNodeNumRows(); ++row) {
02300         RowInfo rowInfo = getRowInfo(row);
02301         sortRowIndices(rowInfo);
02302       }
02303     }
02304     // we just sorted every row
02305     setSorted(true);
02306   }
02307 
02308 
02311   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
02312   void CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::makeColMap()
02313   {
02314     std::string tfecfFuncName("makeColMap()");
02315     typedef OrdinalTraits<GlobalOrdinal> GOT;
02316     //
02317     if (hasColMap()) return;
02318     const size_t nlrs = getNodeNumRows();
02319     //
02320     computeIndexState();
02321     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(isLocallyIndexed() == true, std::runtime_error, ": indices must still be global when making the column map.");
02322     // ultimate goal: list of indices for column map
02323     Array<GlobalOrdinal> myColumns;
02324     // if isGloballyIndexed() == false and isLocallyIndexed() == false, then we have nothing to do
02325     if (isGloballyIndexed() == true) {
02326       // Construct two lists of columns for which we have a non-zero
02327       // Local GIDs: Column GIDs which are locally present on the Domain map
02328       // Remote GIDs: Column GIDs which are not locally present present on the Domain map
02329       //
02330       // instead of list of local GIDs, we will use a set<LocalOrdinal> LocalGID.
02331       //
02332       const LocalOrdinal LINV = OrdinalTraits<LocalOrdinal>::invalid();
02333       size_t numLocalColGIDs = 0, numRemoteColGIDs = 0;
02334       //
02335       // initial: partitioning into local and remote
02336       Array<char> GIDisLocal(domainMap_->getNodeNumElements(),0);
02337       std::set<GlobalOrdinal> RemoteGIDSet;
02338       for (size_t r=0; r < nlrs; ++r) {
02339         RowInfo rowinfo = getRowInfo(r);
02340         if (rowinfo.numEntries > 0) {
02341           ArrayView<const GlobalOrdinal> rowgids = getGlobalView(rowinfo);
02342           rowgids = rowgids(0,rowinfo.numEntries);
02343           for (typename ArrayView<const GlobalOrdinal>::iterator cind = rowgids.begin(); cind != rowgids.end(); ++cind)
02344           {
02345             GlobalOrdinal gid = (*cind);
02346             LocalOrdinal lid = domainMap_->getLocalElement(gid);
02347             if (lid != LINV) {
02348               char alreadyFound = GIDisLocal[lid];
02349               if (alreadyFound == 0) {
02350                 GIDisLocal[lid] = 1;
02351                 ++numLocalColGIDs;
02352               }
02353             }
02354             else {
02355               std::pair<typename std::set<GlobalOrdinal>::iterator, bool> ip;
02356               ip = RemoteGIDSet.insert(gid);
02357               if (ip.second == true) { // gid did not exist in the set and was actually inserted
02358                 ++numRemoteColGIDs;
02359               }
02360             }
02361           }
02362         }
02363       }
02364 
02365       // Possible short-circuit for serial scenario
02366       // If the all domain GIDs are present as column indices, then set ColMap=DomainMap
02367       // By construction, LocalGIDs \subset DomainGIDs
02368       // If we have
02369       //   * Number of remote GIDs is 0, so that ColGIDs == LocalGIDs,
02370       // and
02371       //   * Number of local GIDs is number of domain GIDs
02372       // then
02373       //   * LocalGIDs \subset DomainGIDs && size(LocalGIDs) == size(DomainGIDs) => DomainGIDs == LocalGIDs == ColGIDs
02374       // on this node.
02375       // We will concern ourselves only with the special case of a serial DomainMap, obviating the need for a communication.
02376       // If
02377       //   * DomainMap has a serial comm
02378       // then we can set Column map as Domain map and return. Benefit: this graph won't need an Import object later.
02379       //
02380       // Note, for a serial domain map, there can be no RemoteGIDs, because there are no remote nodes.
02381       // Likely explanations for this are:
02382       //  * user submitted erroneous column indices
02383       //  * user submitted erroneous domain map
02384       if (Teuchos::size(*domainMap_->getComm()) == 1)
02385       {
02386         TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(numRemoteColGIDs != 0, std::runtime_error,
02387             ": Some column IDs are not in the domain map." << std::endl
02388             << "Either these column IDs are invalid or the domain map is invalid." << std::endl
02389             << "Remember, for a rectangular matrix, the domain map must be passed to fillComplete().");
02390         if (numLocalColGIDs == domainMap_->getNodeNumElements()) {
02391           colMap_ = domainMap_;
02392           checkInternalState();
02393           return;
02394         }
02395       }
02396 
02397       // Now, populate myColumns() with a list of all column GIDs.
02398       // Put local GIDs at the front: they correspond to "same" and "permuted" entries between the column map and the domain map
02399       // Put remote GIDs at the back
02400       myColumns.resize(numLocalColGIDs + numRemoteColGIDs);
02401       // get pointers into myColumns for each part
02402       ArrayView<GlobalOrdinal> LocalColGIDs, RemoteColGIDs;
02403       LocalColGIDs  = myColumns(0,numLocalColGIDs);
02404       RemoteColGIDs = myColumns(numLocalColGIDs,numRemoteColGIDs);
02405 
02406       // Copy the Remote GIDs into myColumns
02407       std::copy(RemoteGIDSet.begin(), RemoteGIDSet.end(), RemoteColGIDs.begin());
02408       // We will make a list of corresponding node IDs
02409       Array<int> RemoteImageIDs(numRemoteColGIDs);
02410       // Lookup the Remote Nodes IDs in the Domain map
02411       {
02412         LookupStatus stat = domainMap_->getRemoteIndexList(RemoteColGIDs, RemoteImageIDs());
02413         // This error check is crucial: this tells us that one of the remote indices was not present in the domain map.
02414         // This means that the Import object cannot be constructed, because of incongruity between the column map and domain map.
02415         //   * The user has made a mistake in the column indices
02416         //   * The user has made a mistake w.r.t. the domain map
02417         // Same error message as above for serial case.
02418         char missingID_lcl = (stat == IDNotPresent ? 1 : 0);
02419         char missingID_gbl;
02420         Teuchos::reduceAll<int,char>(*getComm(),Teuchos::REDUCE_MAX,missingID_lcl,
02421           outArg(missingID_gbl));
02422         TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(missingID_gbl == 1, std::runtime_error,
02423             ": Some column IDs are not in the domain map." << std::endl
02424             << "Either these column IDs are invalid or the domain map is invalid." << std::endl
02425             << "Common cause: for a rectangular matrix, the domain map must be passed to fillComplete().");
02426       }
02427       // Sort External column indices so that all columns coming from a given remote processor are contiguous
02428       // This obviates the need for the Distributor associated with the Import from having to reorder data.
02429       sort2(RemoteImageIDs.begin(), RemoteImageIDs.end(), RemoteColGIDs.begin());
02430 
02431       // Copy the Local GIDs into myColumns. Two cases:
02432       // (1) If the number of Local column GIDs is the same as the number of Local domain GIDs, we
02433       //     can simply read the domain GIDs into the front part of ColIndices (see logic above from the serial short circuit case)
02434       // (2) We step through the GIDs of the DomainMap, checking to see if each domain GID is a column GID.
02435       //     we want to do this to maintain a consistent ordering of GIDs between the columns and the domain.
02436       ArrayView<const GlobalOrdinal> mge = domainMap_->getNodeElementList();
02437       if (numLocalColGIDs == domainMap_->getNodeNumElements()) {
02438         std::copy(mge.begin(), mge.end(), LocalColGIDs.begin());
02439       }
02440       else {
02441         size_t numlocalagain = 0;
02442         for (size_t i=0; i < domainMap_->getNodeNumElements(); ++i) {
02443           if (GIDisLocal[i]) {
02444             LocalColGIDs[numlocalagain++] = mge[i];
02445           }
02446         }
02447         TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(numlocalagain != numLocalColGIDs, std::logic_error, ": Internal logic error. Please contact Tpetra team.");
02448       }
02449     }
02450     colMap_ = rcp(new Map<LocalOrdinal,GlobalOrdinal,Node>(GOT::invalid(), myColumns, domainMap_->getIndexBase(), domainMap_->getComm(), domainMap_->getNode()) );
02451     checkInternalState();
02452     return;
02453   }
02454 
02455 
02456   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
02457   void CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::mergeAllIndices()
02458   {
02459     TEUCHOS_TEST_FOR_EXCEPT( isGloballyIndexed() ); // call only after makeIndicesLocal()
02460     TEUCHOS_TEST_FOR_EXCEPT( ! isSorted() ); // call only after sortIndices()
02461     if (! isMerged ()) {
02462       for (size_t row=0; row < getNodeNumRows(); ++row) {
02463         RowInfo rowInfo = getRowInfo(row);
02464         mergeRowIndices(rowInfo);
02465       }
02466       // we just merged every row
02467       setMerged(true);
02468     }
02469   }
02470 
02471 
02472   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
02473   void
02474   CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::makeImportExport()
02475   {
02476     typedef Import<LocalOrdinal,GlobalOrdinal,Node> import_type;
02477     typedef Export<LocalOrdinal,GlobalOrdinal,Node> export_type;
02478     TEUCHOS_TEST_FOR_EXCEPTION(! hasColMap (), std::logic_error, "Tpetra::"
02479       "CrsGraph: It's not allowed to call makeImportExport() unless the graph "
02480       "has a column Map.");
02481     RCP<ParameterList> params = this->getNonconstParameterList (); // could be null
02482     // Create the Import instance if necessary.
02483     if (domainMap_ != colMap_ && (! domainMap_->isSameAs (*colMap_))) {
02484       if (params.is_null () || ! params->isSublist ("Import")) {
02485         importer_ = rcp (new import_type (domainMap_, colMap_));
02486       }
02487       else {
02488         RCP<ParameterList> importSublist = sublist (params, "Import", true);
02489         importer_ = rcp (new import_type (domainMap_, colMap_, importSublist));
02490       }
02491     }
02492     else {
02493       importer_ = null;
02494     }
02495     // Create the Export instance if necessary.
02496     if (rangeMap_ != rowMap_ && (!rangeMap_->isSameAs(*rowMap_))) {
02497       if (params.is_null () || ! params->isSublist ("Export")) {
02498         exporter_ = rcp (new export_type (rowMap_, rangeMap_));
02499       }
02500       else {
02501         RCP<ParameterList> exportSublist = sublist (params, "Export", true);
02502         exporter_ = rcp (new export_type (rowMap_, rangeMap_, exportSublist));
02503       }
02504     }
02505     else {
02506       exporter_ = null;
02507     }
02508   }
02509 
02510 
02511   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
02512   std::string
02513   CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::description() const
02514   {
02515     std::ostringstream oss;
02516     oss << DistObject<GlobalOrdinal,LocalOrdinal,GlobalOrdinal,Node>::description();
02517     if (isFillComplete()) {
02518       oss << "{status = fill complete"
02519           << ", global rows = " << getGlobalNumRows()
02520           << ", global cols = " << getGlobalNumCols()
02521           << ", global num entries = " << getGlobalNumEntries()
02522           << "}";
02523     }
02524     else {
02525       oss << "{status = fill not complete"
02526           << ", global rows = " << getGlobalNumRows()
02527           << "}";
02528     }
02529     return oss.str();
02530   }
02531 
02532 
02533   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
02534   void
02535   CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::
02536   describe (Teuchos::FancyOStream &out,
02537             const Teuchos::EVerbosityLevel verbLevel) const
02538   {
02539     std::string tfecfFuncName("describe()");
02540     using std::endl;
02541     using std::setw;
02542     using Teuchos::VERB_DEFAULT;
02543     using Teuchos::VERB_NONE;
02544     using Teuchos::VERB_LOW;
02545     using Teuchos::VERB_MEDIUM;
02546     using Teuchos::VERB_HIGH;
02547     using Teuchos::VERB_EXTREME;
02548     Teuchos::EVerbosityLevel vl = verbLevel;
02549     if (vl == VERB_DEFAULT) vl = VERB_LOW;
02550     RCP<const Comm<int> > comm = this->getComm();
02551     const int myImageID = comm->getRank(),
02552               numImages = comm->getSize();
02553     size_t width = 1;
02554     for (size_t dec=10; dec<getGlobalNumRows(); dec *= 10) {
02555       ++width;
02556     }
02557     width = std::max<size_t>(width,Teuchos::as<size_t>(11)) + 2;
02558     Teuchos::OSTab tab(out);
02559     //    none: print nothing
02560     //     low: print O(1) info from node 0
02561     //  medium: print O(P) info, num entries per node
02562     //    high: print O(N) info, num entries per row
02563     // extreme: print O(NNZ) info: print graph indices
02564     //
02565     // for medium and higher, print constituent objects at specified verbLevel
02566     if (vl != VERB_NONE) {
02567       if (myImageID == 0) out << this->description() << std::endl;
02568       // O(1) globals, minus what was already printed by description()
02569       if (isFillComplete() && myImageID == 0) {
02570         out << "Global number of diagonals = " << globalNumDiags_ << std::endl;
02571         out << "Global max number of row entries = " << globalMaxNumRowEntries_ << std::endl;
02572       }
02573       // constituent objects
02574       if (vl == VERB_MEDIUM || vl == VERB_HIGH || vl == VERB_EXTREME) {
02575         if (myImageID == 0) out << "\nRow map: " << std::endl;
02576         rowMap_->describe(out,vl);
02577         if (colMap_ != null) {
02578           if (myImageID == 0) out << "\nColumn map: " << std::endl;
02579           colMap_->describe(out,vl);
02580         }
02581         if (domainMap_ != null) {
02582           if (myImageID == 0) out << "\nDomain map: " << std::endl;
02583           domainMap_->describe(out,vl);
02584         }
02585         if (rangeMap_ != null) {
02586           if (myImageID == 0) out << "\nRange map: " << std::endl;
02587           rangeMap_->describe(out,vl);
02588         }
02589       }
02590       // O(P) data
02591       if (vl == VERB_MEDIUM || vl == VERB_HIGH || vl == VERB_EXTREME) {
02592         for (int imageCtr = 0; imageCtr < numImages; ++imageCtr) {
02593           if (myImageID == imageCtr) {
02594             out << "Node ID = " << imageCtr << std::endl
02595                 << "Node number of entries = " << nodeNumEntries_ << std::endl
02596                 << "Node number of diagonals = " << nodeNumDiags_ << std::endl
02597                 << "Node max number of entries = " << nodeMaxNumRowEntries_ << std::endl;
02598             if (indicesAreAllocated()) {
02599               out << "Node number of allocated entries = " << nodeNumAllocated_ << std::endl;
02600             }
02601             else {
02602               out << "Indices are not allocated." << std::endl;
02603             }
02604           }
02605           comm->barrier();
02606           comm->barrier();
02607           comm->barrier();
02608         }
02609       }
02610       // O(N) and O(NNZ) data
02611       if (vl == VERB_HIGH || vl == VERB_EXTREME) {
02612         TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(hasRowInfo() == false, std::runtime_error, ": reduce verbosity level; graph row information was deleted at fillComplete().");
02613         for (int imageCtr = 0; imageCtr < numImages; ++imageCtr) {
02614           if (myImageID == imageCtr) {
02615             out << std::setw(width) << "Node ID"
02616                 << std::setw(width) << "Global Row"
02617                 << std::setw(width) << "Num Entries";
02618             if (vl == VERB_EXTREME) {
02619               out << "Entries";
02620             }
02621             out << std::endl;
02622             for (size_t r=0; r < getNodeNumRows(); ++r) {
02623               RowInfo rowinfo = getRowInfo(r);
02624               GlobalOrdinal gid = rowMap_->getGlobalElement(r);
02625               out << std::setw(width) << myImageID
02626                   << std::setw(width) << gid
02627                   << std::setw(width) << rowinfo.numEntries;
02628               if (vl == VERB_EXTREME) {
02629                 if (isGloballyIndexed()) {
02630                   ArrayView<const GlobalOrdinal> rowview = getGlobalView(rowinfo);
02631                   for (size_t j=0; j < rowinfo.numEntries; ++j) out << rowview[j] << " ";
02632                 }
02633                 else if (isLocallyIndexed()) {
02634                   ArrayView<const LocalOrdinal> rowview = getLocalView(rowinfo);
02635                   for (size_t j=0; j < rowinfo.numEntries; ++j) out << colMap_->getGlobalElement(rowview[j]) << " ";
02636                 }
02637               }
02638               out << std::endl;
02639             }
02640           }
02641           comm->barrier();
02642           comm->barrier();
02643           comm->barrier();
02644         }
02645       }
02646     }
02647   }
02648 
02649 
02650   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
02651   bool CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::checkSizes(const DistObject<GlobalOrdinal,LocalOrdinal,GlobalOrdinal,Node>& source)
02652   {
02653     // It's not clear what kind of compatibility checks on sizes can be performed here.
02654     // Epetra_CrsGraph doesn't check any sizes for compatibility.
02655     return true;
02656   }
02657 
02658 
02659   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
02660   void CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::copyAndPermute(
02661                           const DistObject<GlobalOrdinal,LocalOrdinal,GlobalOrdinal,Node> & source,
02662                           size_t numSameIDs,
02663                           const ArrayView<const LocalOrdinal> &permuteToLIDs,
02664                           const ArrayView<const LocalOrdinal> &permuteFromLIDs)
02665   {
02666     const CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>& src_graph = dynamic_cast<const CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>&>(source);
02667     std::string tfecfFuncName("copyAndPermute()");
02668     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(permuteToLIDs.size() != permuteFromLIDs.size(), std::runtime_error, ": permuteToLIDs and permuteFromLIDs must have the same size.");
02669     bool src_filled = src_graph.isFillComplete();
02670 
02671     Array<GlobalOrdinal> row_copy;
02672     LocalOrdinal myid = 0;
02673     for (size_t i=0; i<numSameIDs; ++i, ++myid) {
02674       GlobalOrdinal gid = src_graph.getMap()->getGlobalElement(myid);
02675       if (src_filled) {
02676         size_t row_length = src_graph.getNumEntriesInGlobalRow(gid);
02677         row_copy.resize(row_length);
02678         size_t check_row_length = 0;
02679         src_graph.getGlobalRowCopy(gid, row_copy(), check_row_length);
02680         insertGlobalIndices( gid, row_copy() );
02681       }
02682       else {
02683         ArrayView<const GlobalOrdinal> row;
02684         src_graph.getGlobalRowView( gid,row );
02685         insertGlobalIndices( gid, row );
02686       }
02687     }
02688 
02689     for (LocalOrdinal i=0; i<permuteToLIDs.size(); ++i) {
02690       GlobalOrdinal mygid = this->getMap()->getGlobalElement(permuteToLIDs[i]);
02691       GlobalOrdinal srcgid= src_graph.getMap()->getGlobalElement(permuteFromLIDs[i]);
02692       if (src_filled) {
02693         size_t row_length = src_graph.getNumEntriesInGlobalRow(srcgid);
02694         row_copy.resize(row_length);
02695         size_t check_row_length = 0;
02696         src_graph.getGlobalRowCopy(srcgid, row_copy(), check_row_length);
02697         insertGlobalIndices( mygid, row_copy() );
02698       }
02699       else {
02700         ArrayView<const GlobalOrdinal> row;
02701         src_graph.getGlobalRowView( srcgid, row );
02702         insertGlobalIndices( mygid, row );
02703       }
02704     }
02705   }
02706 
02707 
02708   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
02709   void
02710   CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::
02711   packAndPrepare (const DistObject<GlobalOrdinal,LocalOrdinal,GlobalOrdinal,Node> & source,
02712                   const ArrayView<const LocalOrdinal> &exportLIDs,
02713                   Array<GlobalOrdinal> &exports,
02714                   const ArrayView<size_t> & numPacketsPerLID,
02715                   size_t& constantNumPackets,
02716                   Distributor &distor)
02717   {
02718     std::string tfecfFuncName("packAndPrepare()");
02719     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(exportLIDs.size() != numPacketsPerLID.size(), std::runtime_error,
02720         ": exportLIDs and numPacketsPerLID must have the same size.");
02721     const CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>& src_graph = dynamic_cast<const CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>&>(source);
02722     // We don't check whether src_graph has had fillComplete called, because it doesn't matter whether the
02723     // *source* graph has been fillComplete'd. The target graph can not be fillComplete'd yet.
02724     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(this->isFillComplete() == true, std::runtime_error,
02725          ": import/export operations are not allowed on destination CrsGraph after fillComplete has been called.");
02726     constantNumPackets = 0;
02727     // first set the contents of numPacketsPerLID, and accumulate a total-num-packets:
02728     size_t totalNumPackets = 0;
02729     Array<GlobalOrdinal> row;
02730     for (LocalOrdinal i=0; i<exportLIDs.size(); ++i) {
02731       GlobalOrdinal GID = src_graph.getMap()->getGlobalElement(exportLIDs[i]);
02732       size_t row_length = src_graph.getNumEntriesInGlobalRow(GID);
02733       numPacketsPerLID[i] = row_length;
02734       totalNumPackets += row_length;
02735     }
02736 
02737     exports.resize(totalNumPackets);
02738 
02739     // now loop again and pack rows of indices into exports:
02740     size_t exportsOffset = 0;
02741     for (LocalOrdinal i=0; i<exportLIDs.size(); ++i) {
02742       GlobalOrdinal GID = src_graph.getMap()->getGlobalElement(exportLIDs[i]);
02743       size_t row_length = src_graph.getNumEntriesInGlobalRow(GID);
02744       row.resize(row_length);
02745       size_t check_row_length = 0;
02746       src_graph.getGlobalRowCopy(GID, row(), check_row_length);
02747       typename Array<GlobalOrdinal>::const_iterator
02748         row_iter = row.begin(), row_end = row.end();
02749       size_t j = 0;
02750       for (; row_iter != row_end; ++row_iter, ++j) {
02751         exports[exportsOffset+j] = *row_iter;
02752       }
02753       exportsOffset += row.size();
02754     }
02755   }
02756 
02757 
02758   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
02759   void CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::unpackAndCombine(
02760                             const ArrayView<const LocalOrdinal> &importLIDs,
02761                             const ArrayView<const GlobalOrdinal> &imports,
02762                             const ArrayView<size_t> &numPacketsPerLID,
02763                             size_t constantNumPackets,
02764                             Distributor & /* distor */,
02765                             CombineMode /* CM */)
02766   {
02767     // FIXME (mfh 02 Apr 2012) REPLACE combine mode has a perfectly
02768     // reasonable meaning, whether or not the matrix is fill complete.
02769     // It's just more work to implement.
02770 
02771     // We are not checking the value of the CombineMode input-argument.
02772     // For CrsGraph, we only support import/export operations if fillComplete has not yet been called.
02773     // Any incoming column-indices are inserted into the target graph. In this context, CombineMode values
02774     // of ADD vs INSERT are equivalent. What is the meaning of REPLACE for CrsGraph? If a duplicate column-index
02775     // is inserted, it will be compressed out when fillComplete is called.
02776     //
02777     // Note: I think REPLACE means that an existing row is replaced by the imported row, i.e., the existing indices are cleared. CGB, 6/17/2010
02778 
02779     std::string tfecfFuncName("unpackAndCombine()");
02780     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(importLIDs.size() != numPacketsPerLID.size(), std::runtime_error,
02781         ": importLIDs and numPacketsPerLID must have the same size.");
02782     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(this->isFillComplete() == true, std::runtime_error,
02783         ": import/export operations are not allowed on destination CrsGraph after fillComplete has been called.");
02784     size_t importsOffset = 0;
02785     typename ArrayView<const LocalOrdinal>::iterator
02786       impLIDiter = importLIDs.begin(), impLIDend = importLIDs.end();
02787     size_t i = 0;
02788     for (; impLIDiter != impLIDend; ++impLIDiter, ++i) {
02789       LocalOrdinal row_length = numPacketsPerLID[i];
02790       const ArrayView<const GlobalOrdinal> row(&imports[importsOffset], row_length);
02791       insertGlobalIndices(this->getMap()->getGlobalElement(*impLIDiter), row);
02792       importsOffset += row_length;
02793     }
02794   }
02795 
02796 
02799   //                                                                         //
02800   //                         Deprecated methods                              //
02801   //                                                                         //
02804 
02805 
02806 } // namespace Tpetra
02807 
02808 
02809 //
02810 // Explicit instantiation macro
02811 //
02812 // Must be expanded from within the Tpetra namespace!
02813 //
02814 
02815 #define TPETRA_CRSGRAPH_INSTANT(LO,GO,NODE) \
02816   \
02817   template class CrsGraph< LO , GO , NODE >; \
02818 
02819 
02820 #endif // TPETRA_CRSGRAPH_DEF_HPP
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines