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 <Tpetra_Distributor.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 
00065   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00066   CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::
00067   CrsGraph (const RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > &rowMap,
00068             size_t maxNumEntriesPerRow,
00069             ProfileType pftype,
00070             const RCP<ParameterList>& params)
00071   : DistObject<GlobalOrdinal,LocalOrdinal,GlobalOrdinal,Node>(rowMap)
00072   , rowMap_(rowMap)
00073   , nodeNumEntries_(0)
00074   , nodeNumAllocated_(OrdinalTraits<size_t>::invalid())
00075   , pftype_(pftype)
00076   , numAllocForAllRows_(maxNumEntriesPerRow)
00077   , indicesAreAllocated_(false)
00078   , indicesAreLocal_(false)
00079   , indicesAreGlobal_(false)
00080   , fillComplete_(false)
00081   , indicesAreSorted_ (true)
00082   , noRedundancies_ (true)
00083   , haveLocalConstants_ (false)
00084   , haveGlobalConstants_ (false)
00085   , haveRowInfo_(true)
00086   {
00087     typedef Teuchos::OrdinalTraits<size_t> OTST;
00088     staticAssertions();
00089     TEUCHOS_TEST_FOR_EXCEPTION(maxNumEntriesPerRow == OTST::invalid(),
00090       std::invalid_argument, "The allocation hint must be a valid size_t value, "
00091       "which in this case means it must not be Teuchos::OrdinalTraits<size_t>::"
00092       "invalid().");
00093     resumeFill(params);
00094     checkInternalState();
00095   }
00096 
00099   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00100   CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::
00101   CrsGraph (const RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > &rowMap,
00102             const RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > &colMap,
00103             size_t maxNumEntriesPerRow,
00104             ProfileType pftype,
00105             const RCP<ParameterList>& params)
00106   : DistObject<GlobalOrdinal,LocalOrdinal,GlobalOrdinal,Node>(rowMap)
00107   , rowMap_(rowMap)
00108   , colMap_(colMap)
00109   , nodeNumEntries_(0)
00110   , nodeNumAllocated_(OrdinalTraits<size_t>::invalid())
00111   , pftype_(pftype)
00112   , numAllocForAllRows_(maxNumEntriesPerRow)
00113   , indicesAreAllocated_(false)
00114   , indicesAreLocal_(false)
00115   , indicesAreGlobal_(false)
00116   , fillComplete_(false)
00117   , indicesAreSorted_(true)
00118   , noRedundancies_(true)
00119   , haveLocalConstants_ (false)
00120   , haveGlobalConstants_ (false)
00121   , haveRowInfo_(true)
00122   {
00123     typedef Teuchos::OrdinalTraits<size_t> OTST;
00124     staticAssertions();
00125     TEUCHOS_TEST_FOR_EXCEPTION(maxNumEntriesPerRow == OTST::invalid(),
00126       std::invalid_argument, "The allocation hint must be a valid size_t value, "
00127       "which in this case means it must not be Teuchos::OrdinalTraits<size_t>::"
00128       "invalid().");
00129     resumeFill(params);
00130     checkInternalState();
00131   }
00132 
00135   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00136   CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::
00137   CrsGraph (const RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > &rowMap,
00138             const ArrayRCP<const size_t> &NumEntriesPerRowToAlloc,
00139             ProfileType pftype,
00140             const RCP<ParameterList>& params)
00141   : DistObject<GlobalOrdinal,LocalOrdinal,GlobalOrdinal,Node>(rowMap)
00142   , rowMap_(rowMap)
00143   , nodeNumEntries_(0)
00144   , nodeNumAllocated_(OrdinalTraits<size_t>::invalid())
00145   , pftype_(pftype)
00146   , numAllocPerRow_(NumEntriesPerRowToAlloc)
00147   , numAllocForAllRows_(0)
00148   , indicesAreAllocated_(false)
00149   , indicesAreLocal_(false)
00150   , indicesAreGlobal_(false)
00151   , fillComplete_(false)
00152   , indicesAreSorted_(true)
00153   , noRedundancies_(true)
00154   , haveLocalConstants_ (false)
00155   , haveGlobalConstants_ (false)
00156   , haveRowInfo_(true)
00157   {
00158     typedef Teuchos::OrdinalTraits<size_t> OTST;
00159     const char tfecfFuncName[] = "CrsGraph(rowMap,NumEntriesPerRowToAlloc)";
00160     staticAssertions();
00161     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC((size_t)NumEntriesPerRowToAlloc.size() != getNodeNumRows(), std::invalid_argument,
00162         ": NumEntriesPerRowToAlloc must have as many entries as specified by rowMap for this node.");
00163     for (size_t r=0; r < getNodeNumRows(); ++r) {
00164       const size_t curRowCount = NumEntriesPerRowToAlloc[r];
00165       TEUCHOS_TEST_FOR_EXCEPTION(curRowCount == OTST::invalid(),
00166         std::invalid_argument, "NumEntriesPerRowToAlloc[" << r << "] specifies "
00167         "an invalid number of entries (Teuchos::OrdinalTraits<size_t>::"
00168         "invalid()).");
00169     }
00170     resumeFill(params);
00171     checkInternalState();
00172   }
00173 
00174 
00177   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00178   CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::
00179   CrsGraph (const RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > &rowMap,
00180             const RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > &colMap,
00181             const ArrayRCP<const size_t> &NumEntriesPerRowToAlloc,
00182             ProfileType pftype,
00183             const RCP<ParameterList>& params)
00184   : DistObject<GlobalOrdinal,LocalOrdinal,GlobalOrdinal,Node>(rowMap)
00185   , rowMap_(rowMap)
00186   , colMap_(colMap)
00187   , nodeNumEntries_(0)
00188   , nodeNumAllocated_(OrdinalTraits<size_t>::invalid())
00189   , pftype_(pftype)
00190   , numAllocPerRow_(NumEntriesPerRowToAlloc)
00191   , numAllocForAllRows_(0)
00192   , indicesAreAllocated_(false)
00193   , indicesAreLocal_(false)
00194   , indicesAreGlobal_(false)
00195   , fillComplete_(false)
00196   , indicesAreSorted_(true)
00197   , noRedundancies_(true)
00198   , haveLocalConstants_ (false)
00199   , haveGlobalConstants_ (false)
00200   , haveRowInfo_(true)
00201   {
00202     typedef Teuchos::OrdinalTraits<size_t> OTST;
00203     const char tfecfFuncName[] = "CrsGraph(rowMap,colMap,NumEntriesPerRowToAlloc)";
00204     staticAssertions();
00205     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC((size_t)NumEntriesPerRowToAlloc.size() != getNodeNumRows(), std::invalid_argument,
00206         ": NumEntriesPerRowToAlloc must have as many entries as specified by rowMap for this node.");
00207     for (size_t r=0; r < getNodeNumRows(); ++r) {
00208       const size_t curRowCount = NumEntriesPerRowToAlloc[r];
00209       TEUCHOS_TEST_FOR_EXCEPTION(curRowCount == OTST::invalid(),
00210         std::invalid_argument, "NumEntriesPerRowToAlloc[" << r << "] specifies "
00211         "an invalid number of entries (Teuchos::OrdinalTraits<size_t>::"
00212         "invalid()).");
00213     }
00214     resumeFill(params);
00215     checkInternalState();
00216   }
00217 
00218 
00221   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00222   CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::
00223   CrsGraph (const RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > &rowMap,
00224             const RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > &colMap,
00225             const ArrayRCP<size_t> & rowPointers,
00226             const ArrayRCP<LocalOrdinal> & columnIndices,
00227             const RCP<ParameterList>& params)
00228   : DistObject<GlobalOrdinal,LocalOrdinal,GlobalOrdinal,Node>(rowMap)
00229   , rowMap_(rowMap)
00230   , colMap_(colMap)
00231   , nodeNumEntries_(0)
00232   , nodeNumAllocated_(OrdinalTraits<size_t>::invalid())
00233   , pftype_(StaticProfile)
00234   , numAllocForAllRows_(0)
00235   , indicesAreAllocated_(true)
00236   , indicesAreLocal_(true)
00237   , indicesAreGlobal_(false)
00238   , fillComplete_(false)
00239   , indicesAreSorted_(true)
00240   , noRedundancies_(true)
00241   , haveLocalConstants_ (false)
00242   , haveGlobalConstants_ (false)
00243   , haveRowInfo_(true)
00244   {
00245     staticAssertions();
00246     globalNumEntries_ = globalNumDiags_ = globalMaxNumRowEntries_ = OrdinalTraits<global_size_t>::invalid();
00247     setAllIndices(rowPointers,columnIndices);
00248     checkInternalState();
00249   }
00250 
00253   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00254   CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::~CrsGraph()
00255   {}
00256 
00257   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00258   RCP<const ParameterList>
00259   CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::
00260   getValidParameters () const
00261   {
00262     RCP<ParameterList> params = parameterList ("Tpetra::CrsGraph");
00263 
00264     // Make a sublist for the Import.
00265     RCP<ParameterList> importSublist = parameterList ("Import");
00266 
00267     // FIXME (mfh 02 Apr 2012) We should really have the Import and
00268     // Export objects fill in these lists.  However, we don't want to
00269     // create an Import or Export unless we need them.  For now, we
00270     // know that the Import and Export just pass the list directly to
00271     // their Distributor, so we can create a Distributor here
00272     // (Distributor's constructor is a lightweight operation) and have
00273     // it fill in the list.
00274 
00275     // Fill in Distributor default parameters by creating a
00276     // Distributor and asking it to do the work.
00277     Distributor distributor (rowMap_->getComm(), importSublist);
00278     params->set ("Import", *importSublist, "How the Import performs communication.");
00279 
00280     // Make a sublist for the Export.  For now, it's a clone of the
00281     // Import sublist.  It's not a shallow copy, though, since we
00282     // might like the Import to do communication differently than the
00283     // Export.
00284     params->set ("Export", *importSublist, "How the Export performs communication.");
00285 
00286     return params;
00287   }
00288 
00289   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00290   void
00291   CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::
00292   setParameterList (const RCP<ParameterList>& params)
00293   {
00294     RCP<const ParameterList> validParams = getValidParameters ();
00295     params->validateParametersAndSetDefaults (*validParams);
00296     this->setMyParamList (params);
00297   }
00298 
00299   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00300   global_size_t CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getGlobalNumRows() const
00301   {
00302     return rowMap_->getGlobalNumElements();
00303   }
00304 
00305   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00306   global_size_t CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getGlobalNumCols() const
00307   {
00308     const char tfecfFuncName[] = "getGlobalNumCols()";
00309     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
00310       isFillComplete() == false, std::runtime_error,
00311       ": requires domain map, which requires that fillComplete() has been "
00312       "called.");
00313     return getDomainMap()->getGlobalNumElements();
00314   }
00315 
00316   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00317   size_t CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getNodeNumRows() const
00318   {
00319     return rowMap_->getNodeNumElements();
00320   }
00321 
00322   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00323   size_t CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getNodeNumCols() const
00324   {
00325     const char tfecfFuncName[] = "getNodeNumCols()";
00326     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
00327       hasColMap() == false, std::runtime_error,
00328       ": requires column map.  This requires either that a custom column Map "
00329       "was given to the constructor, or that fillComplete() has been called.");
00330     return colMap_->getNodeNumElements();
00331   }
00332 
00333   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00334   size_t CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getNodeNumDiags() const
00335   {
00336     return nodeNumDiags_;
00337   }
00338 
00339   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00340   global_size_t CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getGlobalNumDiags() const
00341   {
00342     return globalNumDiags_;
00343   }
00344 
00345   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00346   RCP<Node>
00347   CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getNode() const
00348   {
00349     return rowMap_->getNode();
00350   }
00351 
00352   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00353   RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> >
00354   CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getRowMap() const
00355   {
00356     return rowMap_;
00357   }
00358 
00359   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00360   RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> >
00361   CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getColMap() const
00362   {
00363     return colMap_;
00364   }
00365 
00366   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00367   RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> >
00368   CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getDomainMap() const
00369   {
00370     return domainMap_;
00371   }
00372 
00373   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00374   RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> >
00375   CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getRangeMap() const
00376   {
00377     return rangeMap_;
00378   }
00379 
00380   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00381   RCP<const Import<LocalOrdinal,GlobalOrdinal,Node> >
00382   CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getImporter() const
00383   {
00384     return importer_;
00385   }
00386 
00387   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00388   RCP<const Export<LocalOrdinal,GlobalOrdinal,Node> >
00389   CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getExporter() const
00390   {
00391     return exporter_;
00392   }
00393 
00394   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00395   bool CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::hasColMap() const
00396   {
00397     return (colMap_ != null);
00398   }
00399 
00400   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00401   bool CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::isStorageOptimized() const
00402   {
00403     bool isOpt = indicesAreAllocated_ && (numRowEntries_ == null) && (getNodeNumRows() > 0);
00404 #ifdef HAVE_TPETRA_DEBUG
00405     const char tfecfFuncName[] = "isStorageOptimized()";
00406     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC( (isOpt == true) && (getProfileType() == DynamicProfile), std::logic_error,
00407         ": Violated stated post-conditions. Please contact Tpetra team.");
00408 #endif
00409     return isOpt;
00410   }
00411 
00412   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00413   ProfileType CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getProfileType() const
00414   {
00415     return pftype_;
00416   }
00417 
00418   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00419   global_size_t CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getGlobalNumEntries() const
00420   {
00421     return globalNumEntries_;
00422   }
00423 
00424   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00425   size_t CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getNodeNumEntries() const
00426   {
00427     return nodeNumEntries_;
00428   }
00429 
00430   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00431   global_size_t CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getGlobalMaxNumRowEntries() const
00432   {
00433     return globalMaxNumRowEntries_;
00434   }
00435 
00436   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00437   size_t CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getNodeMaxNumRowEntries() const
00438   {
00439     return nodeMaxNumRowEntries_;
00440   }
00441 
00442   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00443   bool CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::isFillComplete() const
00444   {
00445     return fillComplete_;
00446   }
00447 
00448   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00449   bool CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::isFillActive() const
00450   {
00451     return !fillComplete_;
00452   }
00453 
00454   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00455   bool CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::isUpperTriangular() const
00456   {
00457     return upperTriangular_;
00458   }
00459 
00460   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00461   bool CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::isLowerTriangular() const
00462   {
00463     return lowerTriangular_;
00464   }
00465 
00466   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00467   bool CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::isLocallyIndexed() const
00468   {
00469     return indicesAreLocal_;
00470   }
00471 
00472   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00473   bool CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::isGloballyIndexed() const
00474   {
00475     return indicesAreGlobal_;
00476   }
00477 
00478   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00479   size_t CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getNodeAllocationSize() const
00480   {
00481     return nodeNumAllocated_;
00482   }
00483 
00484   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00485   RCP<const Comm<int> >
00486   CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getComm() const
00487   {
00488     return rowMap_->getComm();
00489   }
00490 
00491   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00492   GlobalOrdinal CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getIndexBase() const
00493   {
00494     return rowMap_->getIndexBase();
00495   }
00496 
00497   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00498   bool CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::indicesAreAllocated() const
00499   {
00500     return indicesAreAllocated_;
00501   }
00502 
00503 
00506   //                                                                         //
00507   //                    Internal utility methods                             //
00508   //                                                                         //
00511 
00512 
00515   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00516   bool CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::isSorted() const
00517   {
00518     return indicesAreSorted_;
00519   }
00520 
00521 
00524   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00525   bool CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::isMerged() const
00526   {
00527     return noRedundancies_;
00528   }
00529 
00532   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00533   void CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::setLocallyModified ()
00534   {
00535     // FIXME (mfh 07 May 2013) How do we know that the change
00536     // introduced a redundancy, or even that it invalidated the sorted
00537     // order of indices?  CrsGraph has always made this conservative
00538     // guess.  It could be a bit costly to check at insertion time,
00539     // though.
00540     indicesAreSorted_ = false;
00541     noRedundancies_ = false;
00542 
00543     // We've modified the graph, so we'll have to recompute local
00544     // constants like the number of diagonal entries on this process.
00545     haveLocalConstants_ = false;
00546   }
00547 
00550   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00551   void
00552   CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::
00553   allocateIndices (ELocalGlobal lg)
00554   {
00555     // This is a protected function, only callable by us.  If it was
00556     // called incorrectly, it is our fault.  That's why the tests
00557     // below throw std::logic_error instead of std::invalid_argument.
00558     const char tfecfFuncName[] = "allocateIndices()";
00559     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
00560       isLocallyIndexed() && lg==GlobalIndices, std::logic_error,
00561       ": The graph is locally indexed, but Tpetra code is calling this method "
00562       "with lg=GlobalIndices.  Please report this bug to the Tpetra developers.");
00563     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
00564       isGloballyIndexed() && lg==LocalIndices, std::logic_error,
00565       ": The graph is globally indexed, but Tpetra code is calling this method "
00566       "with lg=LocalIndices.  Please report this bug to the Tpetra developers.");
00567     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
00568       indicesAreAllocated(), std::logic_error, ": The graph's indices are "
00569       "already allocated, but Tpetra code is calling allocateIndices() again.  "
00570       "Please report this bug to the Tpetra developers.");
00571 
00572     const size_t numRows = getNodeNumRows();
00573     indicesAreLocal_  = (lg == LocalIndices);
00574     indicesAreGlobal_ = (lg == GlobalIndices);
00575     nodeNumAllocated_ = 0;
00576     if (numAllocPerRow_ == null && getNodeNumRows() > 0) {
00577       // this wastes memory, temporarily, but it simplifies the code
00578       // and interfaces to follow
00579       ArrayRCP<size_t> tmpnumallocperrow = arcp<size_t>(numRows);
00580       std::fill(tmpnumallocperrow.begin(), tmpnumallocperrow.end(), numAllocForAllRows_);
00581       numAllocPerRow_ = tmpnumallocperrow;
00582     }
00583     //
00584     if (getProfileType() == StaticProfile) {
00585       //
00586       //  STATIC ALLOCATION PROFILE
00587       //
00588       // Have the local sparse kernels object allocate row offsets for
00589       // us, with first-touch allocation if applicable.  This is not
00590       // as important for global indices, because we never use global
00591       // indices in sparse kernels, but we might as well use the code
00592       // that we have for both the local and global indices cases.
00593       // Furthermore, first-touch allocation ensures that we don't
00594       // take up too much memory in any one NUMA affinity region.
00595       rowPtrs_ = LocalMatOps::allocRowPtrs( getRowMap()->getNode(), numAllocPerRow_() );
00596       if (lg == LocalIndices) {
00597         lclInds1D_ = LocalMatOps::template allocStorage<LocalOrdinal>( getRowMap()->getNode(), rowPtrs_() );
00598       }
00599       else {
00600         gblInds1D_ = LocalMatOps::template allocStorage<GlobalOrdinal>( getRowMap()->getNode(), rowPtrs_() );
00601       }
00602       nodeNumAllocated_ = rowPtrs_[numRows];
00603     }
00604     else {
00605       //
00606       //  DYNAMIC ALLOCATION PROFILE
00607       //
00608       typename ArrayRCP<const size_t>::iterator numalloc = numAllocPerRow_.begin();
00609       size_t howmany = numAllocForAllRows_;
00610       if (lg == LocalIndices) {
00611         lclInds2D_ = arcp< Array<LocalOrdinal> >(numRows);
00612         nodeNumAllocated_ = 0;
00613         for (size_t i=0; i < numRows; ++i) {
00614           if (numAllocPerRow_ != null) howmany = *numalloc++;
00615           nodeNumAllocated_ += howmany;
00616           if (howmany > 0) lclInds2D_[i].resize(howmany);
00617         }
00618       }
00619       else { // allocate global indices
00620         gblInds2D_ = arcp< Array<GlobalOrdinal> >(numRows);
00621         nodeNumAllocated_ = 0;
00622         for (size_t i=0; i < numRows; ++i) {
00623           if (numAllocPerRow_ != null) howmany = *numalloc++;
00624           nodeNumAllocated_ += howmany;
00625           if (howmany > 0) gblInds2D_[i].resize(howmany);
00626         }
00627       }
00628     }
00629     if (numRows > 0) {
00630       numRowEntries_ = arcp<size_t>(numRows);
00631       std::fill( numRowEntries_.begin(), numRowEntries_.end(), 0 );
00632     }
00633     // done with these
00634     numAllocForAllRows_ = 0;
00635     numAllocPerRow_     = null;
00636     indicesAreAllocated_ = true;
00637     checkInternalState();
00638   }
00639 
00640 
00643   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00644   template <class T>
00645   ArrayRCP<T> CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::allocateValues1D() const
00646   {
00647     const char tfecfFuncName[] = "allocateValues()";
00648     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
00649         indicesAreAllocated() == false,
00650         std::runtime_error, ": graph indices must be allocated before values."
00651     );
00652     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
00653         getProfileType() != StaticProfile,
00654         std::runtime_error, ": graph indices must be allocated in a static profile."
00655     );
00656     ArrayRCP<T> values1D;
00657     values1D = LocalMatOps::template allocStorage<T>( getRowMap()->getNode(), rowPtrs_() );
00658     return values1D;
00659   }
00660 
00661 
00664   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00665   template <class T>
00666   ArrayRCP<Array<T> >
00667   CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::allocateValues2D() const
00668   {
00669     const char tfecfFuncName[] = "allocateValues2D()";
00670     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
00671         indicesAreAllocated() == false,
00672         std::runtime_error, ": graph indices must be allocated before values."
00673     );
00674     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
00675         getProfileType() != DynamicProfile,
00676         std::runtime_error, ": graph indices must be allocated in a dynamic profile."
00677     );
00678     ArrayRCP<Array<T> > values2D;
00679     values2D = arcp<Array<T> > (getNodeNumRows ());
00680     if (lclInds2D_ != null) {
00681       const size_t numRows = lclInds2D_.size ();
00682       for (size_t r = 0; r < numRows; ++r) {
00683         values2D[r].resize (lclInds2D_[r].size ());
00684       }
00685     }
00686     else if (gblInds2D_ != null) {
00687       const size_t numRows = gblInds2D_.size ();
00688       for (size_t r = 0; r < numRows; ++r) {
00689         values2D[r].resize (gblInds2D_[r].size ());
00690       }
00691     }
00692     return values2D;
00693   }
00694 
00695 
00696 //   /////////////////////////////////////////////////////////////////////////////
00697 //   /////////////////////////////////////////////////////////////////////////////
00698 //   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00699 //   template <ELocalGlobal lg, class T>
00700 //   RowInfo
00701 //   CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::
00702 //   updateAllocAndValues (RowInfo rowinfo,
00703 //                         size_t newAllocSize,
00704 //                         Array<T>& rowVals)
00705 //   {
00706 // #ifdef HAVE_TPETRA_DEBUG
00707 //     TEUCHOS_TEST_FOR_EXCEPT( ! rowMap_->isNodeLocalElement(rowinfo.localRow) );
00708 //     TEUCHOS_TEST_FOR_EXCEPT( newAllocSize < rowinfo.allocSize );
00709 //     TEUCHOS_TEST_FOR_EXCEPT( (lg == LocalIndices && ! isLocallyIndexed()) ||
00710 //                              (lg == GlobalIndices && ! isGloballyIndexed()) );
00711 //     TEUCHOS_TEST_FOR_EXCEPT( newAllocSize == 0 );
00712 //     TEUCHOS_TEST_FOR_EXCEPT( ! indicesAreAllocated() );
00713 // #endif
00714 //     // ArrayRCP::resize automatically copies over values on reallocation.
00715 //     if (lg == LocalIndices) {
00716 //       lclInds2D_[rowinfo.localRow].resize (newAllocSize);
00717 //     }
00718 //     else { // lg == GlobalIndices
00719 //       gblInds2D_[rowinfo.localRow].resize (newAllocSize);
00720 //     }
00721 //     rowVals.resize (newAllocSize);
00722 //     nodeNumAllocated_ += (newAllocSize - rowinfo.allocSize);
00723 //     rowinfo.allocSize = newAllocSize;
00724 //     return rowinfo;
00725 //   }
00726 
00727 
00730   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00731   ArrayView<const LocalOrdinal>
00732   CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::
00733   getLocalView (RowInfo rowinfo) const
00734   {
00735     ArrayView<const LocalOrdinal> view;
00736     if (rowinfo.allocSize > 0) {
00737       if (lclInds1D_ != null) {
00738         view = lclInds1D_ (rowinfo.offset1D, rowinfo.allocSize);
00739       }
00740       else if (! lclInds2D_[rowinfo.localRow].empty ()) {
00741         view = lclInds2D_[rowinfo.localRow] ();
00742       }
00743     }
00744     return view;
00745   }
00746 
00747 
00750   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00751   ArrayView<LocalOrdinal>
00752   CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::
00753   getLocalViewNonConst (RowInfo rowinfo)
00754   {
00755     ArrayView<LocalOrdinal> view;
00756     if (rowinfo.allocSize > 0) {
00757       if (lclInds1D_ != null) {
00758         view = lclInds1D_ (rowinfo.offset1D, rowinfo.allocSize);
00759       }
00760       else if (! lclInds2D_[rowinfo.localRow].empty()) {
00761         view = lclInds2D_[rowinfo.localRow] ();
00762       }
00763     }
00764     return view;
00765   }
00766 
00767 
00770   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00771   ArrayView<const GlobalOrdinal>
00772   CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::
00773   getGlobalView (RowInfo rowinfo) const
00774   {
00775     ArrayView<const GlobalOrdinal> view;
00776     if (rowinfo.allocSize > 0) {
00777       if (gblInds1D_ != null) {
00778         view = gblInds1D_ (rowinfo.offset1D, rowinfo.allocSize);
00779       }
00780       else if (! gblInds2D_[rowinfo.localRow].empty()) {
00781         view = gblInds2D_[rowinfo.localRow] ();
00782       }
00783     }
00784     return view;
00785   }
00786 
00787 
00790   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00791   ArrayView<GlobalOrdinal>
00792   CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::
00793   getGlobalViewNonConst (RowInfo rowinfo)
00794   {
00795     ArrayView<GlobalOrdinal> view;
00796     if (rowinfo.allocSize > 0) {
00797       if (gblInds1D_ != null) {
00798         view = gblInds1D_ (rowinfo.offset1D, rowinfo.allocSize);
00799       }
00800       else if (!gblInds2D_[rowinfo.localRow].empty()) {
00801         view = gblInds2D_[rowinfo.localRow] ();
00802       }
00803     }
00804     return view;
00805   }
00806 
00807 
00810   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00811   RowInfo CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getRowInfo(size_t myRow) const
00812   {
00813 #ifdef HAVE_TPETRA_DEBUG
00814     const char tfecfFuncName[] = "getRowInfo";
00815     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
00816         rowMap_->isNodeLocalElement (myRow) == false,
00817         std::logic_error,
00818         ": The given (local) row index myRow = " << myRow
00819         << " does not belong to the graph's row Map.  "
00820         "This probably indicates a bug in Tpetra::CrsGraph or Tpetra::CrsMatrix.  "
00821         "Please report this to the Tpetra developers.");
00822     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
00823       hasRowInfo() == false, std::logic_error,
00824       ": Late catch! Graph does not have row info anymore.  "
00825       "Error should have been caught earlier.  Please contact Tpetra team.");
00826 #endif // HAVE_TPETRA_DEBUG
00827     const size_t STINV = Teuchos::OrdinalTraits<size_t>::invalid ();
00828     RowInfo ret;
00829     ret.localRow = myRow;
00830     if (nodeNumAllocated_ != 0 && nodeNumAllocated_ != STINV) {
00831       // graph data structures have the info that we need
00832       //
00833       // if static graph, offsets tell us the allocation size
00834       if (getProfileType() == StaticProfile) {
00835         ret.offset1D  = rowPtrs_[myRow];
00836         ret.allocSize = rowPtrs_[myRow+1] - rowPtrs_[myRow];
00837         if (numRowEntries_ == null) {
00838           ret.numEntries = ret.allocSize;
00839         }
00840         else {
00841           ret.numEntries = numRowEntries_[myRow];
00842         }
00843       }
00844       else {
00845         ret.offset1D = STINV;
00846         if (isLocallyIndexed ()) {
00847           ret.allocSize = lclInds2D_[myRow].size();
00848         }
00849         else {
00850           ret.allocSize = gblInds2D_[myRow].size();
00851         }
00852         ret.numEntries = numRowEntries_[myRow];
00853       }
00854     }
00855     else if (nodeNumAllocated_ == 0) {
00856       // have performed allocation, but the graph has no allocation or entries
00857       ret.allocSize = 0;
00858       ret.numEntries = 0;
00859       ret.offset1D = STINV;
00860     }
00861     else if (indicesAreAllocated () == false) {
00862       // haven't performed allocation yet; probably won't hit this code
00863       if (numAllocPerRow_ == null) {
00864         ret.allocSize = numAllocForAllRows_;
00865       }
00866       else {
00867         ret.allocSize = numAllocPerRow_[myRow];
00868       }
00869       ret.numEntries = 0;
00870       ret.offset1D = STINV;
00871     }
00872     else {
00873       // don't know how we ended up here...
00874       TEUCHOS_TEST_FOR_EXCEPT(true);
00875     }
00876     return ret;
00877   }
00878 
00879 
00882   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00883   void CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::staticAssertions() const
00884   {
00885     // Assumption: sizeof(GlobalOrdinal) >= sizeof(LocalOrdinal):
00886     //     This is so that we can store local indices in the memory
00887     //     formerly occupied by global indices.
00888     //
00889     // Assumption: max(GlobalOrdinal) >= max(LocalOrdinal) and
00890     //   max(size_t) >= max(LocalOrdinal)
00891     //     This is so that we can represent any LocalOrdinal as a
00892     //     size_t, and any LocalOrdinal as a GlobalOrdinal
00893     Teuchos::CompileTimeAssert<sizeof(GlobalOrdinal) < sizeof(LocalOrdinal)> cta_size1; (void)cta_size1;
00894     Teuchos::CompileTimeAssert<sizeof(global_size_t) < sizeof(size_t)      > cta_size2; (void)cta_size2;
00895     // can't call max() with CompileTimeAssert, because it isn't a constant expression; will need to make this a runtime check
00896     std::string msg = typeName(*this) + ": Object cannot be allocated with stated template arguments: size assumptions are not valid.";
00897     TEUCHOS_TEST_FOR_EXCEPTION( (size_t)OrdinalTraits<LocalOrdinal>::max() > OrdinalTraits<size_t>::max(),          std::runtime_error, msg);
00898     TEUCHOS_TEST_FOR_EXCEPTION( (global_size_t)OrdinalTraits<LocalOrdinal>::max() > (global_size_t)OrdinalTraits<GlobalOrdinal>::max(),           std::runtime_error, msg);
00899     TEUCHOS_TEST_FOR_EXCEPTION( (size_t)OrdinalTraits<GlobalOrdinal>::max() > OrdinalTraits<global_size_t>::max(),  std::runtime_error, msg);
00900     TEUCHOS_TEST_FOR_EXCEPTION( OrdinalTraits<size_t>::max() > OrdinalTraits<global_size_t>::max(),                 std::runtime_error, msg);
00901   }
00902 
00903 
00906   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00907   template <ELocalGlobal lg>
00908   size_t CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::filterIndices(const SLocalGlobalNCViews &inds) const
00909   {
00910     const Map<LocalOrdinal,GlobalOrdinal,Node> &cmap = *colMap_;
00911     Teuchos::CompileTimeAssert<lg != GlobalIndices && lg != LocalIndices> cta_lg;
00912     (void)cta_lg;
00913     size_t numFiltered = 0;
00914 #ifdef HAVE_TPETRA_DEBUG
00915     size_t numFiltered_debug = 0;
00916 #endif
00917     if (lg == GlobalIndices) {
00918       ArrayView<GlobalOrdinal> ginds = inds.ginds;
00919       typename ArrayView<GlobalOrdinal>::iterator fend = ginds.begin(),
00920                                                   cptr = ginds.begin();
00921       while (cptr != ginds.end()) {
00922         if (cmap.isNodeGlobalElement(*cptr)) {
00923           *fend++ = *cptr;
00924 #ifdef HAVE_TPETRA_DEBUG
00925           ++numFiltered_debug;
00926 #endif
00927         }
00928         ++cptr;
00929       }
00930       numFiltered = fend - ginds.begin();
00931     }
00932     else if (lg == LocalIndices) {
00933       ArrayView<LocalOrdinal> linds = inds.linds;
00934       typename ArrayView<LocalOrdinal>::iterator fend = linds.begin(),
00935                                                  cptr = linds.begin();
00936       while (cptr != linds.end()) {
00937         if (cmap.isNodeLocalElement(*cptr)) {
00938           *fend++ = *cptr;
00939 #ifdef HAVE_TPETRA_DEBUG
00940           ++numFiltered_debug;
00941 #endif
00942         }
00943         ++cptr;
00944       }
00945       numFiltered = fend - linds.begin();
00946     }
00947 #ifdef HAVE_TPETRA_DEBUG
00948     TEUCHOS_TEST_FOR_EXCEPT( numFiltered != numFiltered_debug );
00949 #endif
00950     return numFiltered;
00951   }
00952 
00953 
00954 //   /////////////////////////////////////////////////////////////////////////////
00955 //   /////////////////////////////////////////////////////////////////////////////
00956 //   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00957 //   template <class T>
00958 //   size_t CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::
00959 //   filterGlobalIndicesAndValues (const ArrayView<GlobalOrdinal>& ginds,
00960 //                                 const ArrayView<T>& vals) const
00961 //   {
00962 //     const Map<LocalOrdinal,GlobalOrdinal,Node>& cmap = *colMap_;
00963 //     size_t numFiltered = 0;
00964 //     typename ArrayView<T>::iterator fvalsend = vals.begin();
00965 //     typename ArrayView<T>::iterator valscptr = vals.begin();
00966 // #ifdef HAVE_TPETRA_DEBUG
00967 //     size_t numFiltered_debug = 0;
00968 // #endif
00969 //     typename ArrayView<GlobalOrdinal>::iterator fend = ginds.begin();
00970 //     typename ArrayView<GlobalOrdinal>::iterator cptr = ginds.begin();
00971 //     while (cptr != ginds.end()) {
00972 //       if (cmap.isNodeGlobalElement (*cptr)) {
00973 //         *fend++ = *cptr;
00974 //         *fvalsend++ = *valscptr;
00975 // #ifdef HAVE_TPETRA_DEBUG
00976 //         ++numFiltered_debug;
00977 // #endif
00978 //       }
00979 //       ++cptr;
00980 //       ++valscptr;
00981 //     }
00982 //     numFiltered = fend - ginds.begin();
00983 // #ifdef HAVE_TPETRA_DEBUG
00984 //     TEUCHOS_TEST_FOR_EXCEPT( numFiltered != numFiltered_debug );
00985 //     TEUCHOS_TEST_FOR_EXCEPT( valscptr != vals.end() );
00986 //     const size_t numFilteredActual =
00987 //       Teuchos::as<size_t> (fvalsend - vals.begin ());
00988 //     TEUCHOS_TEST_FOR_EXCEPT( numFiltered != numFilteredActual );
00989 // #endif
00990 //     return numFiltered;
00991 //   }
00992 
00993 
00994 //   /////////////////////////////////////////////////////////////////////////////
00995 //   /////////////////////////////////////////////////////////////////////////////
00996 //   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00997 //   template <class T>
00998 //   size_t
00999 //   CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::
01000 //   filterLocalIndicesAndValues (const ArrayView<LocalOrdinal>& linds,
01001 //                                const ArrayView<T>& vals) const
01002 //   {
01003 //     const Map<LocalOrdinal,GlobalOrdinal,Node>& cmap = *colMap_;
01004 //     size_t numFiltered = 0;
01005 //     typename ArrayView<T>::iterator fvalsend = vals.begin();
01006 //     typename ArrayView<T>::iterator valscptr = vals.begin();
01007 // #ifdef HAVE_TPETRA_DEBUG
01008 //     size_t numFiltered_debug = 0;
01009 // #endif
01010 //     typename ArrayView<LocalOrdinal>::iterator fend = linds.begin();
01011 //     typename ArrayView<LocalOrdinal>::iterator cptr = linds.begin();
01012 //     while (cptr != linds.end()) {
01013 //       if (cmap.isNodeLocalElement (*cptr)) {
01014 //         *fend++ = *cptr;
01015 //         *fvalsend++ = *valscptr;
01016 // #ifdef HAVE_TPETRA_DEBUG
01017 //         ++numFiltered_debug;
01018 // #endif
01019 //       }
01020 //       ++cptr;
01021 //       ++valscptr;
01022 //     }
01023 //     numFiltered = fend - linds.begin();
01024 // #ifdef HAVE_TPETRA_DEBUG
01025 //     TEUCHOS_TEST_FOR_EXCEPT( numFiltered != numFiltered_debug );
01026 //     TEUCHOS_TEST_FOR_EXCEPT( valscptr != vals.end() );
01027 //     const size_t numFilteredActual =
01028 //       Teuchos::as<size_t> (fvalsend - vals.begin ());
01029 //     TEUCHOS_TEST_FOR_EXCEPT( numFiltered != numFilteredActual );
01030 // #endif
01031 //     return numFiltered;
01032 //   }
01033 
01034 
01037   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
01038   size_t
01039   CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::
01040   insertIndices (const RowInfo& rowinfo,
01041                  const SLocalGlobalViews &newInds,
01042                  const ELocalGlobal lg,
01043                  const ELocalGlobal I)
01044   {
01045 #ifdef HAVE_TPETRA_DEBUG
01046     TEUCHOS_TEST_FOR_EXCEPTION(
01047       lg != GlobalIndices && lg != LocalIndices, std::invalid_argument,
01048       "Tpetra::CrsGraph::insertIndices: lg must be either GlobalIndices or "
01049       "LocalIndices.");
01050 #endif // HAVE_TPETRA_DEBUG
01051     size_t numNewInds = 0;
01052     if (lg == GlobalIndices) { // input indices are global
01053       ArrayView<const GlobalOrdinal> new_ginds = newInds.ginds;
01054       numNewInds = new_ginds.size();
01055       if (I == GlobalIndices) { // store global indices
01056         ArrayView<GlobalOrdinal> gind_view = getGlobalViewNonConst(rowinfo);
01057         std::copy(new_ginds.begin(), new_ginds.end(), gind_view.begin()+rowinfo.numEntries);
01058       }
01059       else if (I == LocalIndices) { // store local indices
01060         ArrayView<LocalOrdinal> lind_view = getLocalViewNonConst(rowinfo);
01061         typename ArrayView<const GlobalOrdinal>::iterator         in = new_ginds.begin();
01062         const typename ArrayView<const GlobalOrdinal>::iterator stop = new_ginds.end();
01063         typename ArrayView<LocalOrdinal>::iterator out = lind_view.begin()+rowinfo.numEntries;
01064         while (in != stop) {
01065           *out++ = colMap_->getLocalElement (*in++);
01066         }
01067       }
01068     }
01069     else if (lg == LocalIndices) { // input indices are local
01070       ArrayView<const LocalOrdinal> new_linds = newInds.linds;
01071       numNewInds = new_linds.size();
01072       if (I == LocalIndices) { // store local indices
01073         ArrayView<LocalOrdinal> lind_view = getLocalViewNonConst(rowinfo);
01074         std::copy(new_linds.begin(), new_linds.end(), lind_view.begin()+rowinfo.numEntries);
01075       }
01076       else if (I == GlobalIndices) {
01077         TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error, "Tpetra::CrsGraph::"
01078           "insertIndices: the case where the input indices are local and the "
01079           "indices to write are global (lg=LocalIndices, I=GlobalIndices) is "
01080           "not implemented, because it does not make sense." << std::endl <<
01081           "If you have correct local column indices, that means the graph has "
01082           "a column Map.  In that case, you should be storing local indices.");
01083       }
01084     }
01085     numRowEntries_[rowinfo.localRow] += numNewInds;
01086     nodeNumEntries_ += numNewInds;
01087     setLocallyModified ();
01088     return numNewInds;
01089   }
01090 
01093   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
01094   void CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::
01095   insertGlobalIndicesImpl (const LocalOrdinal myRow,
01096                            const ArrayView<const GlobalOrdinal> &indices)
01097   {
01098     const char* tfecfFuncName ("insertGlobalIndicesImpl");
01099 
01100     RowInfo rowInfo = getRowInfo(myRow);
01101     const size_t numNewInds = indices.size();
01102     const size_t newNumEntries = rowInfo.numEntries + numNewInds;
01103     if (newNumEntries > rowInfo.allocSize) {
01104       TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
01105         getProfileType() == StaticProfile, std::runtime_error,
01106         ": new indices exceed statically allocated graph structure.");
01107 
01108       // update allocation, doubling size to reduce number of reallocations
01109       size_t newAllocSize = 2*rowInfo.allocSize;
01110       if (newAllocSize < newNumEntries)
01111         newAllocSize = newNumEntries;
01112       gblInds2D_[myRow].resize(newAllocSize);
01113       nodeNumAllocated_ += (newAllocSize - rowInfo.allocSize);
01114     }
01115 
01116     // Copy new indices at end of global index array
01117     if (gblInds1D_ != null)
01118       std::copy(indices.begin(), indices.end(),
01119                 gblInds1D_.begin()+rowInfo.offset1D+rowInfo.numEntries);
01120     else
01121       std::copy(indices.begin(), indices.end(),
01122                 gblInds2D_[myRow].begin()+rowInfo.numEntries);
01123     numRowEntries_[myRow] += numNewInds;
01124     nodeNumEntries_ += numNewInds;
01125     setLocallyModified ();
01126 
01127 #ifdef HAVE_TPETRA_DEBUG
01128     {
01129       const size_t chkNewNumEntries = getNumEntriesInLocalRow (myRow);
01130       TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
01131         chkNewNumEntries != newNumEntries, std::logic_error,
01132         ": Internal logic error. Please contact Tpetra team.");
01133     }
01134 #endif
01135   }
01136 
01139   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
01140   void CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::
01141   insertLocalIndicesImpl (const LocalOrdinal myRow,
01142                           const ArrayView<const LocalOrdinal> &indices)
01143   {
01144     const char* tfecfFuncName ("insertLocallIndicesImpl");
01145 
01146     RowInfo rowInfo = getRowInfo(myRow);
01147     const size_t numNewInds = indices.size();
01148     const size_t newNumEntries = rowInfo.numEntries + numNewInds;
01149     if (newNumEntries > rowInfo.allocSize) {
01150       TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
01151         getProfileType() == StaticProfile, std::runtime_error,
01152         ": new indices exceed statically allocated graph structure.");
01153 
01154       // update allocation, doubling size to reduce number of reallocations
01155       size_t newAllocSize = 2*rowInfo.allocSize;
01156       if (newAllocSize < newNumEntries)
01157         newAllocSize = newNumEntries;
01158       lclInds2D_[myRow].resize(newAllocSize);
01159       nodeNumAllocated_ += (newAllocSize - rowInfo.allocSize);
01160     }
01161 
01162     // Insert new indices at end of lclInds array
01163     if (lclInds1D_ != null)
01164       std::copy(indices.begin(), indices.end(),
01165                 lclInds1D_.begin()+rowInfo.offset1D+rowInfo.numEntries);
01166     else
01167       std::copy(indices.begin(), indices.end(),
01168                 lclInds2D_[myRow].begin()+rowInfo.numEntries);
01169     numRowEntries_[myRow] += numNewInds;
01170     nodeNumEntries_ += numNewInds;
01171     setLocallyModified ();
01172 #ifdef HAVE_TPETRA_DEBUG
01173     {
01174       const size_t chkNewNumEntries = getNumEntriesInLocalRow (myRow);
01175       TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
01176         chkNewNumEntries != newNumEntries, std::logic_error,
01177         ": Internal logic error. Please contact Tpetra team.");
01178     }
01179 #endif
01180   }
01181 
01182   // /////////////////////////////////////////////////////////////////////////////
01183   // /////////////////////////////////////////////////////////////////////////////
01184   // template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
01185   // template <class Scalar>
01186   // void CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::
01187   // insertIndicesAndValues (const RowInfo& rowInfo,
01188   //                         const SLocalGlobalViews& newInds,
01189   //                         const ArrayView<Scalar>& oldRowVals,
01190   //                         const ArrayView<const Scalar>& newRowVals,
01191   //                         const ELocalGlobal lg,
01192   //                         const ELocalGlobal I)
01193   // {
01194   //   const size_t numNewInds = insertIndices (rowInfo, newInds, lg, I);
01195   //   typename ArrayView<const Scalar>::const_iterator newRowValsBegin =
01196   //     newRowVals.begin ();
01197   //   std::copy (newRowValsBegin, newRowValsBegin + numNewInds,
01198   //              oldRowVals.begin () + rowInfo.numEntries);
01199   // }
01200 
01201   // /////////////////////////////////////////////////////////////////////////////
01202   // /////////////////////////////////////////////////////////////////////////////
01203   // template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
01204   // template <class Scalar, class BinaryFunction>
01205   // void
01206   // CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::
01207   // transformLocalValues (RowInfo rowInfo,
01208   //                       const Teuchos::ArrayView<Scalar>& rowVals,
01209   //                       const Teuchos::ArrayView<const LocalOrdinal>& inds,
01210   //                       const Teuchos::ArrayView<const Scalar>& newVals,
01211   //                       BinaryFunction f) const
01212   // {
01213   //   const size_t STINV = Teuchos::OrdinalTraits<size_t>::invalid();
01214   //   const size_t numElts = Teuchos::as<size_t> (inds.size ());
01215   //   size_t hint = 0; // Guess for the current index k into rowVals
01216 
01217   //   // Get a view of the column indices in the row.  This amortizes
01218   //   // the cost of getting the view over all the entries of inds.
01219   //   ArrayView<const LocalOrdinal> colInds = getLocalView (rowInfo);
01220 
01221   //   for (size_t j = 0; j < numElts; ++j) {
01222   //     const size_t k = findLocalIndex (rowInfo, inds[j], colInds, hint);
01223   //     if (k != STINV) {
01224   //       rowVals[k] = f( rowVals[k], newVals[j] );
01225   //       hint = k+1;
01226   //     }
01227   //   }
01228   // }
01229 
01230 
01231   // /////////////////////////////////////////////////////////////////////////////
01232   // /////////////////////////////////////////////////////////////////////////////
01233   // template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
01234   // template <class Scalar, class BinaryFunction>
01235   // void
01236   // CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::
01237   // transformGlobalValues (RowInfo rowInfo,
01238   //                        const Teuchos::ArrayView<Scalar>& rowVals,
01239   //                        const Teuchos::ArrayView<const GlobalOrdinal>& inds,
01240   //                        const Teuchos::ArrayView<const Scalar>& newVals,
01241   //                        BinaryFunction f) const
01242   // {
01243   //   const size_t STINV = Teuchos::OrdinalTraits<size_t>::invalid();
01244   //   const size_t numElts = Teuchos::as<size_t> (inds.size ());
01245   //   size_t hint = 0; // hint is a guess as to wheter the index is
01246 
01247   //   for (size_t j = 0; j < numElts; ++j) {
01248   //     const size_t k = findGlobalIndex (rowInfo, inds[j], hint);
01249   //     if (k != STINV) {
01250   //       rowVals[k] = f( rowVals[k], newVals[j] );
01251   //       hint = k+1;
01252   //     }
01253   //   }
01254   // }
01255 
01256 
01257 
01260   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
01261   void CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::sortRowIndices(RowInfo rowinfo)
01262   {
01263     if (rowinfo.numEntries > 0) {
01264       ArrayView<LocalOrdinal> inds_view = getLocalViewNonConst(rowinfo);
01265       std::sort(inds_view.begin(), inds_view.begin() + rowinfo.numEntries);
01266     }
01267   }
01268 
01269 
01272   // in the future, perhaps this could use std::sort with a boost::zip_iterator
01273   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
01274   template <class Scalar>
01275   void CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::sortRowIndicesAndValues(RowInfo rowinfo, ArrayView<Scalar> values)
01276   {
01277     if (rowinfo.numEntries > 0) {
01278       ArrayView<LocalOrdinal> inds_view = getLocalViewNonConst(rowinfo);
01279       sort2(inds_view.begin(), inds_view.begin()+rowinfo.numEntries, values.begin());
01280     }
01281   }
01282 
01283 
01286   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
01287   void CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::mergeRowIndices(RowInfo rowinfo)
01288   {
01289     const char tfecfFuncName[] = "mergRowIndices()";
01290     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
01291       isStorageOptimized() == true, std::logic_error,
01292       ": The graph is already storage optimized, so we shouldn't be merging any indices."
01293       " Please report this bug to the Tpetra developers.");
01294     ArrayView<LocalOrdinal> inds_view = getLocalViewNonConst(rowinfo);
01295     typename ArrayView<LocalOrdinal>::iterator beg, end, newend;
01296     beg = inds_view.begin();
01297     end = inds_view.begin() + rowinfo.numEntries;
01298     newend = std::unique(beg,end);
01299     const size_t mergedEntries = newend - beg;
01300 #ifdef HAVE_TPETRA_DEBUG
01301     // merge should not have eliminated any entries; if so, the assignment below will destory the packed structure
01302     TEUCHOS_TEST_FOR_EXCEPT( isStorageOptimized() && mergedEntries != rowinfo.numEntries );
01303 #endif
01304     numRowEntries_[rowinfo.localRow] = mergedEntries;
01305     nodeNumEntries_ -= (rowinfo.numEntries - mergedEntries);
01306   }
01307 
01308 
01311   // in the future, this could use std::unique with a boost::zip_iterator
01312   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
01313   template<class Scalar>
01314   void CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::
01315   mergeRowIndicesAndValues (RowInfo rowinfo, const ArrayView<Scalar>& rowValues)
01316   {
01317     const char tfecfFuncName[] = "mergeRowIndicesAndValues";
01318     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
01319       isStorageOptimized(), std::logic_error, ": It is invalid to call this "
01320       "method if the graph's storage has already been optimized." << std::endl
01321       << "Please report this bug to the Tpetra developers.");
01322 
01323     typedef typename ArrayView<Scalar>::iterator Iter;
01324     Iter rowValueIter = rowValues.begin ();
01325     ArrayView<LocalOrdinal> inds_view = getLocalViewNonConst (rowinfo);
01326     typename ArrayView<LocalOrdinal>::iterator beg, end, newend;
01327 
01328     // beg,end define a half-exclusive interval over which to iterate.
01329     beg = inds_view.begin();
01330     end = inds_view.begin() + rowinfo.numEntries;
01331     newend = beg;
01332     if (beg != end) {
01333       typename ArrayView<LocalOrdinal>::iterator cur = beg + 1;
01334       Iter vcur = rowValueIter + 1;
01335       Iter vend = rowValueIter;
01336       cur = beg+1;
01337       while (cur != end) {
01338         if (*cur != *newend) {
01339           // new entry; save it
01340           ++newend;
01341           ++vend;
01342           (*newend) = (*cur);
01343           (*vend) = (*vcur);
01344         }
01345         else {
01346           // old entry; merge it
01347           //(*vend) = f (*vend, *vcur);
01348           (*vend) += *vcur;
01349         }
01350         ++cur;
01351         ++vcur;
01352       }
01353       ++newend; // one past the last entry, per typical [beg,end) semantics
01354     }
01355     const size_t mergedEntries = newend - beg;
01356 #ifdef HAVE_TPETRA_DEBUG
01357     // merge should not have eliminated any entries; if so, the
01358     // assignment below will destroy the packed structure
01359     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
01360       isStorageOptimized() && mergedEntries != rowinfo.numEntries,
01361       std::logic_error,
01362       ": Merge was incorrect; it eliminated entries from the graph.  "
01363       << std::endl << "Please report this bug to the Tpetra developers.");
01364 #endif // HAVE_TPETRA_DEBUG
01365     numRowEntries_[rowinfo.localRow] = mergedEntries;
01366     nodeNumEntries_ -= (rowinfo.numEntries - mergedEntries);
01367   }
01368 
01369 
01372   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
01373   void CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::
01374   setDomainRangeMaps (const RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> >& domainMap,
01375                       const RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> >& rangeMap)
01376   {
01377     // simple pointer comparison for equality
01378     if (domainMap_ != domainMap) {
01379       domainMap_ = domainMap;
01380       importer_ = null;
01381     }
01382     if (rangeMap_ != rangeMap) {
01383       rangeMap_  = rangeMap;
01384       exporter_ = null;
01385     }
01386   }
01387 
01388 
01391   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
01392   size_t
01393   CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::
01394   findLocalIndex (RowInfo rowinfo, LocalOrdinal ind, size_t hint) const
01395   {
01396     ArrayView<const LocalOrdinal> colInds = getLocalView (rowinfo);
01397     return this->findLocalIndex (rowinfo, ind, colInds, hint);
01398   }
01399 
01400 
01403   template <class LocalOrdinal,
01404             class GlobalOrdinal,
01405             class Node,
01406             class LocalMatOps>
01407   size_t
01408   CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::
01409   findLocalIndex (RowInfo rowinfo,
01410                   LocalOrdinal ind,
01411                   ArrayView<const LocalOrdinal> colInds,
01412                   size_t hint) const
01413   {
01414     typedef typename ArrayView<const LocalOrdinal>::iterator IT;
01415 
01416     // If the hint was correct, then the hint is the offset to return.
01417     if (hint < rowinfo.numEntries && colInds[hint] == ind) {
01418       return hint;
01419     }
01420 
01421     // The hint was wrong, so we must search for the given column
01422     // index in the column indices for the given row.  How we do the
01423     // search depends on whether the graph's column indices are
01424     // sorted.
01425     IT beg = colInds.begin ();
01426     IT end = beg + rowinfo.numEntries;
01427     IT ptr = beg + rowinfo.numEntries; // "null"
01428     bool found = true;
01429 
01430     if (isSorted ()) {
01431       // binary search
01432       std::pair<IT,IT> p = std::equal_range (beg, end, ind);
01433       if (p.first == p.second) {
01434         found = false;
01435       }
01436       else {
01437         ptr = p.first;
01438       }
01439     }
01440     else {
01441       // direct search
01442       ptr = std::find (beg, end, ind);
01443       if (ptr == end) {
01444         found = false;
01445       }
01446     }
01447 
01448     if (found) {
01449       return Teuchos::as<size_t> (ptr - beg);
01450     }
01451     else {
01452       return Teuchos::OrdinalTraits<size_t>::invalid ();
01453     }
01454   }
01455 
01456 
01459   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
01460   size_t
01461   CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::
01462   findGlobalIndex (RowInfo rowinfo, GlobalOrdinal ind, size_t hint) const
01463   {
01464     typedef typename ArrayView<const GlobalOrdinal>::iterator IT;
01465     ArrayView<const GlobalOrdinal> indices = getGlobalView (rowinfo);
01466 
01467     // We don't actually require that the hint be a valid index.
01468     // If it is not in range, we just ignore it.
01469     if (hint < rowinfo.numEntries && indices[hint] == ind) {
01470       return hint;
01471     }
01472 
01473     IT beg = indices.begin ();
01474     IT end = indices.begin () + rowinfo.numEntries; // not indices.end()
01475     if (isSorted ()) { // use binary search
01476       const std::pair<IT,IT> p = std::equal_range (beg, end, ind);
01477       if (p.first == p.second) { // range of matching entries is empty
01478         return Teuchos::OrdinalTraits<size_t>::invalid ();
01479       } else {
01480         return p.first - beg;
01481       }
01482     }
01483     else { // not sorted; must use linear search
01484       const IT loc = std::find (beg, end, ind);
01485       if (loc == end) {
01486         return Teuchos::OrdinalTraits<size_t>::invalid ();
01487       } else {
01488         return loc - beg;
01489       }
01490     }
01491   }
01492 
01493 
01496   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
01497   void CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::clearGlobalConstants()
01498   {
01499     globalNumEntries_       = OrdinalTraits<global_size_t>::invalid();
01500     globalNumDiags_         = OrdinalTraits<global_size_t>::invalid();
01501     globalMaxNumRowEntries_ = OrdinalTraits<global_size_t>::invalid();
01502     haveGlobalConstants_    = false;
01503   }
01504 
01505 
01508   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
01509   void CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::
01510   checkInternalState () const
01511   {
01512 #ifdef HAVE_TPETRA_DEBUG
01513     const global_size_t GSTI = OrdinalTraits<global_size_t>::invalid();
01514     const size_t         STI = OrdinalTraits<size_t>::invalid();
01515     std::string err = typeName(*this) + "::checkInternalState(): Likely internal logic error. Please contact Tpetra team.";
01516     // check the internal state of this data structure
01517     // this is called by numerous state-changing methods, in a debug build, to ensure that the object
01518     // always remains in a valid state
01519     // the graph should have been allocated with a row map
01520     TEUCHOS_TEST_FOR_EXCEPTION( rowMap_ == null,                     std::logic_error, err );
01521     // am either complete or active
01522     TEUCHOS_TEST_FOR_EXCEPTION( isFillActive() == isFillComplete(),  std::logic_error, err );
01523     // if active, i have no local graph
01524     TEUCHOS_TEST_FOR_EXCEPTION( isFillActive() && lclGraph_ != null, std::logic_error, err );
01525     // if the graph has been fill completed, then all maps should be present
01526     TEUCHOS_TEST_FOR_EXCEPTION( isFillComplete() == true && (colMap_ == null || rangeMap_ == null || domainMap_ == null), std::logic_error, err );
01527     // if storage has been optimized, then indices should have been allocated (even if trivially so)
01528     TEUCHOS_TEST_FOR_EXCEPTION( isStorageOptimized() == true && indicesAreAllocated() == false, std::logic_error, err );
01529     // if storage has been optimized, then number of allocated is now the number of entries
01530     TEUCHOS_TEST_FOR_EXCEPTION( isStorageOptimized() == true && getNodeAllocationSize() != getNodeNumEntries(), std::logic_error, err );
01531     // if graph doesn't have the global constants, then they should all be marked as invalid
01532     TEUCHOS_TEST_FOR_EXCEPTION( haveGlobalConstants_ == false && ( globalNumEntries_ != GSTI || globalNumDiags_ != GSTI || globalMaxNumRowEntries_ != GSTI ), std::logic_error, err );
01533     // if the graph has global cosntants, then they should be valid.
01534     TEUCHOS_TEST_FOR_EXCEPTION( haveGlobalConstants_ == true && ( globalNumEntries_ == GSTI || globalNumDiags_ == GSTI || globalMaxNumRowEntries_ == GSTI ), std::logic_error, err );
01535     TEUCHOS_TEST_FOR_EXCEPTION( haveGlobalConstants_ == true && ( globalNumEntries_ < nodeNumEntries_ || globalNumDiags_ < nodeNumDiags_ || globalMaxNumRowEntries_ < nodeMaxNumRowEntries_ ),
01536                         std::logic_error, err );
01537     // if indices are allocated, then the allocation specifications should have been released
01538     TEUCHOS_TEST_FOR_EXCEPTION( indicesAreAllocated() == true  && (numAllocForAllRows_ != 0 || numAllocPerRow_ != null),                        std::logic_error, err );
01539     // if indices are not allocated, then information dictating allocation quantities should be present
01540     TEUCHOS_TEST_FOR_EXCEPTION( indicesAreAllocated() == false && (nodeNumAllocated_ != STI || nodeNumEntries_ != 0),                           std::logic_error, err );
01541     // if storage is optimized, then profile should be static
01542     TEUCHOS_TEST_FOR_EXCEPTION( isStorageOptimized() && pftype_ != StaticProfile,                                                               std::logic_error, err );
01543     // if rowPtrs_ exists, it is required to have N+1 rows, and rowPtrs_[N] == gblInds1D_.size()/lclInds1D_.size()
01544     TEUCHOS_TEST_FOR_EXCEPTION( isGloballyIndexed() && rowPtrs_ != null && ((size_t)rowPtrs_.size() != getNodeNumRows()+1 || rowPtrs_[getNodeNumRows()] != (size_t)gblInds1D_.size()), std::logic_error, err );
01545     TEUCHOS_TEST_FOR_EXCEPTION(  isLocallyIndexed() && rowPtrs_ != null && ((size_t)rowPtrs_.size() != getNodeNumRows()+1 || rowPtrs_[getNodeNumRows()] != (size_t)lclInds1D_.size()), std::logic_error, err );
01546     // if profile is dynamic and we have allocated, then 2D allocations should be present
01547     TEUCHOS_TEST_FOR_EXCEPTION( pftype_ == DynamicProfile && indicesAreAllocated() && getNodeNumRows() > 0 && lclInds2D_ == null && gblInds2D_ == null,
01548                                                                                                                                         std::logic_error, err );
01549     // if profile is dynamic, then numentries and 2D indices are needed and should be present
01550     TEUCHOS_TEST_FOR_EXCEPTION( pftype_ == DynamicProfile && indicesAreAllocated() && getNodeNumRows() > 0 && (numRowEntries_ == null || (lclInds2D_ == null && gblInds2D_ == null)),
01551                                                                                                                                         std::logic_error, err );
01552     // if profile is dynamic, then 1D allocations should not be present
01553     TEUCHOS_TEST_FOR_EXCEPTION( pftype_ == DynamicProfile && (lclInds1D_ != null || gblInds1D_ != null),                                std::logic_error, err );
01554     // if profile is dynamic, then row offsets should not be present
01555     TEUCHOS_TEST_FOR_EXCEPTION( pftype_ == DynamicProfile && rowPtrs_ != null,                                                          std::logic_error, err );
01556     // if profile is static and we have allocated non-trivially, then 1D allocations should be present
01557     TEUCHOS_TEST_FOR_EXCEPTION( pftype_ == StaticProfile && indicesAreAllocated() && getNodeAllocationSize() > 0 && lclInds1D_ == null && gblInds1D_ == null,
01558                                                                                                                                         std::logic_error, err );
01559     // if profile is static, then 2D allocations should not be present
01560     TEUCHOS_TEST_FOR_EXCEPTION( pftype_ == StaticProfile && (lclInds2D_ != null || gblInds2D_ != null),                                 std::logic_error, err );
01561     // if indices are not allocated, then none of the buffers should be.
01562     TEUCHOS_TEST_FOR_EXCEPTION( indicesAreAllocated() == false && (rowPtrs_ != null || numRowEntries_ != null ||
01563                                                                  lclInds1D_ != null || lclInds2D_ != null ||
01564                                                                  gblInds1D_ != null || gblInds2D_ != null),                             std::logic_error, err );
01565     // indices may be local or global only if they are allocated (numAllocated is redundant; could simply be indicesAreLocal_ || indicesAreGlobal_)
01566     TEUCHOS_TEST_FOR_EXCEPTION( (indicesAreLocal_ == true || indicesAreGlobal_ == true) && indicesAreAllocated_ == false,               std::logic_error, err );
01567     // indices may be local or global, but not both
01568     TEUCHOS_TEST_FOR_EXCEPTION( indicesAreLocal_ == true && indicesAreGlobal_ == true,                                                  std::logic_error, err );
01569     // if indices are local, then global allocations should not be present
01570     TEUCHOS_TEST_FOR_EXCEPTION( indicesAreLocal_ == true && (gblInds1D_ != null || gblInds2D_ != null),                                 std::logic_error, err );
01571     // if indices are global, then local allocations should not be present
01572     TEUCHOS_TEST_FOR_EXCEPTION( indicesAreGlobal_ == true && (lclInds1D_ != null || lclInds2D_ != null),                                std::logic_error, err );
01573     // if indices are local, then local allocations should be present
01574     TEUCHOS_TEST_FOR_EXCEPTION( indicesAreLocal_ == true && getNodeAllocationSize() > 0 && lclInds1D_ == null && getNodeNumRows() > 0 && lclInds2D_ == null,
01575                                                                                                                           std::logic_error, err );
01576     // if indices are global, then global allocations should be present
01577     TEUCHOS_TEST_FOR_EXCEPTION( indicesAreGlobal_ == true && getNodeAllocationSize() > 0 && gblInds1D_ == null && getNodeNumRows() > 0 && gblInds2D_ == null,
01578                                                                                                                           std::logic_error, err );
01579     // if indices are allocated, then we should have recorded how many were allocated
01580     TEUCHOS_TEST_FOR_EXCEPTION( indicesAreAllocated() == true  && nodeNumAllocated_ == STI,                                       std::logic_error, err );
01581     // if indices are not allocated, then the allocation size should be marked invalid
01582     TEUCHOS_TEST_FOR_EXCEPTION( indicesAreAllocated() == false && nodeNumAllocated_ != STI,                                       std::logic_error, err );
01583     // check the actual allocations
01584     if (indicesAreAllocated()) {
01585       size_t actualNumAllocated = 0;
01586       if (pftype_ == DynamicProfile) {
01587         if (isGloballyIndexed() && gblInds2D_ != null) {
01588           for (size_t r = 0; r < getNodeNumRows(); ++r) {
01589             actualNumAllocated += gblInds2D_[r].size();
01590           }
01591         }
01592         else if (isLocallyIndexed() && lclInds2D_ != null) {
01593           for (size_t r = 0; r < getNodeNumRows(); ++r) {
01594             actualNumAllocated += lclInds2D_[r].size();
01595           }
01596         }
01597         TEUCHOS_TEST_FOR_EXCEPTION(actualNumAllocated != nodeNumAllocated_, std::logic_error, err );
01598       }
01599       else if (rowPtrs_ != null) { // pftype_ == StaticProfile
01600         actualNumAllocated = rowPtrs_[getNodeNumRows()];
01601         TEUCHOS_TEST_FOR_EXCEPTION(  isLocallyIndexed() == true && (size_t)lclInds1D_.size() != actualNumAllocated, std::logic_error, err );
01602         TEUCHOS_TEST_FOR_EXCEPTION( isGloballyIndexed() == true && (size_t)gblInds1D_.size() != actualNumAllocated, std::logic_error, err );
01603         TEUCHOS_TEST_FOR_EXCEPTION(actualNumAllocated != nodeNumAllocated_, std::logic_error, err );
01604       }
01605     }
01606 #endif
01607   }
01608 
01609 
01612   //                                                                         //
01613   //                  User-visible class methods                             //
01614   //                                                                         //
01617 
01618 
01621   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
01622   size_t CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getNumEntriesInGlobalRow(GlobalOrdinal globalRow) const
01623   {
01624     const LocalOrdinal lrow = rowMap_->getLocalElement(globalRow);
01625     size_t ret = OrdinalTraits<size_t>::invalid();
01626     if (hasRowInfo() && lrow != OrdinalTraits<LocalOrdinal>::invalid())
01627     {
01628       RowInfo rowinfo = getRowInfo(lrow);
01629       ret = rowinfo.numEntries;
01630     }
01631     return ret;
01632   }
01633 
01634 
01637   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
01638   size_t CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getNumEntriesInLocalRow(LocalOrdinal localRow) const
01639   {
01640     size_t ret = OrdinalTraits<size_t>::invalid();
01641     if (hasRowInfo() && rowMap_->isNodeLocalElement(localRow)) {
01642       RowInfo rowinfo = getRowInfo(localRow);
01643       ret = rowinfo.numEntries;
01644     }
01645     return ret;
01646   }
01647 
01648 
01651   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
01652   size_t CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getNumAllocatedEntriesInGlobalRow(GlobalOrdinal globalRow) const
01653   {
01654     const LocalOrdinal lrow = rowMap_->getLocalElement(globalRow);
01655     size_t ret = OrdinalTraits<size_t>::invalid();
01656     if (hasRowInfo() && lrow != OrdinalTraits<LocalOrdinal>::invalid())
01657     {
01658       RowInfo rowinfo = getRowInfo(lrow);
01659       ret = rowinfo.allocSize;
01660     }
01661     return ret;
01662   }
01663 
01664 
01667   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
01668   size_t CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getNumAllocatedEntriesInLocalRow(LocalOrdinal localRow) const
01669   {
01670     size_t ret = OrdinalTraits<size_t>::invalid();
01671     if (hasRowInfo() && rowMap_->isNodeLocalElement(localRow)) {
01672       RowInfo rowinfo = getRowInfo(localRow);
01673       ret = rowinfo.allocSize;
01674     }
01675     return ret;
01676   }
01677 
01678 
01681   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
01682   ArrayRCP<const size_t> CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getNodeRowPtrs() const
01683   {
01684     return rowPtrs_;
01685   }
01686 
01687 
01690   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
01691   ArrayRCP<const LocalOrdinal> CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getNodePackedIndices() const
01692   {
01693     return lclInds1D_;
01694   }
01695 
01696 
01699   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
01700   void
01701   CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::
01702   getLocalRowCopy (LocalOrdinal localRow,
01703                    const ArrayView<LocalOrdinal>& indices,
01704                    size_t& NumIndices) const
01705   {
01706     const char tfecfFuncName[] = "getLocalRowCopy";
01707     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
01708       isGloballyIndexed () && ! hasColMap (), std::runtime_error,
01709       ": The graph is globally indexed, and does not have a column Map (yet), "
01710       "so it is impossible to return local indices.");
01711     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
01712       ! rowMap_->isNodeLocalElement (localRow), std::runtime_error,
01713       ": localRow (== " << localRow << ") is not valid on the calling process "
01714       "with rank " << getComm ()->getRank () << ".");
01715     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
01716       ! hasRowInfo (), std::runtime_error,
01717       ": graph row information was deleted at fillComplete().");
01718 
01719     const RowInfo rowinfo = getRowInfo (localRow);
01720     NumIndices = rowinfo.numEntries;
01721     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
01722       static_cast<size_t> (indices.size ()) < NumIndices, std::runtime_error,
01723       ": specified storage (size == " << indices.size () << ") does not suffice "
01724       "to hold all entries for this row (NumIndices == " << NumIndices << ").");
01725 
01726     if (isLocallyIndexed ()) {
01727       ArrayView<const LocalOrdinal> lview = getLocalView (rowinfo);
01728       std::copy (lview.begin (), lview.begin () + NumIndices, indices.begin ());
01729     }
01730     else if (isGloballyIndexed ()) {
01731       ArrayView<const GlobalOrdinal> gview = getGlobalView (rowinfo);
01732       for (size_t j = 0; j < NumIndices; ++j) {
01733         indices[j] = colMap_->getLocalElement (gview[j]);
01734       }
01735     }
01736     else {
01737       // If the graph on the calling process is neither locally nor
01738       // globally indexed, that means it owns no column indices.
01739       //
01740       // FIXME (mfh 21 Oct 2013) It's not entirely clear to me whether
01741       // we can reach this branch, given the checks above.  However,
01742       // if that is the case, it should still be correct to call this
01743       // function if the calling process owns no column indices.
01744       NumIndices = 0;
01745     }
01746   }
01747 
01748 
01751   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
01752   void
01753   CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::
01754   getGlobalRowCopy (GlobalOrdinal globalRow,
01755                     const ArrayView<GlobalOrdinal>& indices,
01756                     size_t& NumIndices) const
01757   {
01758     // One of the following is true:
01759     //
01760     // 1. The calling process owns no column indices.
01761     // 2. The graph stores global column indices.
01762     // 3. The graph has a column Map that we can use to convert any
01763     //    stored local indices into global indices for output.
01764 
01765     const char tfecfFuncName[] = "getGlobalRowCopy";
01766     const LocalOrdinal localRow = rowMap_->getLocalElement (globalRow);
01767     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
01768       localRow == Teuchos::OrdinalTraits<LocalOrdinal>::invalid(),
01769       std::runtime_error, ": globalRow (== " << globalRow << ") is not owned "
01770       "by (i.e., is not in the row Map on) the calling process with rank "
01771       << getComm ()->getRank () << ".");
01772     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
01773       ! hasRowInfo (), std::runtime_error,
01774       ": graph row information was deleted at fillComplete().");
01775 
01776     const RowInfo rowinfo = getRowInfo (static_cast<size_t> (localRow));
01777     NumIndices = rowinfo.numEntries;
01778     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
01779       static_cast<size_t> (indices.size ()) < NumIndices, std::runtime_error,
01780       ": specified storage (size == " << indices.size () << ") does not suffice "
01781       "to hold all entries in this row (NumIndices == " << NumIndices << ").");
01782 
01783     if (isLocallyIndexed ()) {
01784       ArrayView<const LocalOrdinal> lview = getLocalView (rowinfo);
01785       for (size_t j = 0; j < NumIndices; ++j) {
01786         indices[j] = colMap_->getGlobalElement (lview[j]);
01787       }
01788     }
01789     else if (isGloballyIndexed ()) {
01790       ArrayView<const GlobalOrdinal> gview = getGlobalView (rowinfo);
01791       std::copy (gview.begin (), gview.begin () + NumIndices, indices.begin ());
01792     }
01793   }
01794 
01795 
01798   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
01799   void
01800   CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::
01801   getLocalRowView (LocalOrdinal localRow, ArrayView<const LocalOrdinal> &indices) const
01802   {
01803     const char tfecfFuncName[] = "getLocalRowView";
01804     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
01805       isGloballyIndexed(), std::runtime_error,
01806       ": The graph is globally indexed, so we cannot return a view with local "
01807       "column indices.  Use getLocalRowCopy() instead.");
01808     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
01809       ! hasRowInfo (), std::runtime_error,
01810       ": graph row information was deleted at fillComplete().");
01811 
01812     indices = null;
01813     if (rowMap_->isNodeLocalElement (localRow)) {
01814       const RowInfo rowinfo = getRowInfo(localRow);
01815       if (rowinfo.numEntries > 0) {
01816         indices = getLocalView (rowinfo);
01817         indices = indices (0, rowinfo.numEntries);
01818       }
01819 #ifdef HAVE_TPETRA_DEBUG
01820       const size_t numEntriesInView = static_cast<size_t> (indices.size ());
01821       TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
01822         numEntriesInView != getNumEntriesInLocalRow (localRow),
01823         std::logic_error, ": The number of entries " << numEntriesInView
01824         << " in the view to return must be the same as getNumEntriesInLocalRow"
01825         "(localRow=" << localRow << ") = " << getNumEntriesInLocalRow (localRow)
01826         << ", but this is not the case.  "
01827         "Please report this bug to the Tpetra developers.");
01828 #endif // HAVE_TPETRA_DEBUG
01829     }
01830   }
01831 
01832 
01835   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
01836   void
01837   CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::
01838   getGlobalRowView (GlobalOrdinal globalRow,
01839                     ArrayView<const GlobalOrdinal>& indices) const
01840   {
01841     using Teuchos::as;
01842     const char tfecfFuncName[] = "getGlobalRowView";
01843     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
01844       isLocallyIndexed (), std::runtime_error,
01845       ": The graph is locally indexed, so we cannot return a view with global "
01846       "column indices.  Use getGlobalRowCopy() instead.");
01847     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
01848       ! hasRowInfo (), std::runtime_error,
01849       ": graph row information was deleted at fillComplete().");
01850 
01851     // We don't have to call isNodeGlobalElement() to test whether
01852     // globalRow belongs to the calling process.  That method requires
01853     // a global to local lookup anyway, and getLocalElement() returns
01854     // invalid() anyway if the global index wasn't found.
01855     const LocalOrdinal localRow = rowMap_->getLocalElement (globalRow);
01856     indices = null;
01857     if (localRow != Teuchos::OrdinalTraits<LocalOrdinal>::invalid ()) {
01858       const RowInfo rowInfo = getRowInfo (as<size_t> (localRow));
01859       if (rowInfo.numEntries > 0) {
01860         indices = (getGlobalView (rowInfo)) (0, rowInfo.numEntries);
01861       }
01862     }
01863 #ifdef HAVE_TPETRA_DEBUG
01864     using std::endl;
01865     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
01866       as<size_t> (indices.size ()) != getNumEntriesInGlobalRow (globalRow),
01867       std::logic_error,
01868       ": Violated stated post-conditions:"
01869       << "  indices.size () = " << indices.size () << endl
01870       << "  as<size_t> (indices.size ()) = " << as<size_t> (indices.size ())
01871       << endl << "  getNumEntriesInGlobalRow (globalRow = " << globalRow
01872       << ") = " << getNumEntriesInGlobalRow (globalRow) << endl
01873       << "Please report this bug to the Tpetra developers.");
01874 #endif // HAVE_TPETRA_DEBUG
01875   }
01876 
01877 
01880   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
01881   void
01882   CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::
01883   insertLocalIndices (const LocalOrdinal localRow,
01884                       const ArrayView<const LocalOrdinal> &indices)
01885   {
01886     using Teuchos::ArrayView;
01887     const char tfecfFuncName[] = "insertLocalIndices";
01888 
01889     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
01890       isFillActive() == false, std::runtime_error,
01891       ": requires that fill is active.");
01892     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
01893       isGloballyIndexed() == true, std::runtime_error,
01894       ": graph indices are global; use insertGlobalIndices().");
01895     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
01896       hasColMap() == false, std::runtime_error,
01897       ": cannot insert local indices without a column map.");
01898     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
01899       rowMap_->isNodeLocalElement(localRow) == false, std::runtime_error,
01900       ": row does not belong to this node.");
01901     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
01902       hasRowInfo() == false, std::runtime_error,
01903       ": graph row information was deleted at fillComplete().");
01904     if (! indicesAreAllocated ()) {
01905       allocateIndices (LocalIndices);
01906     }
01907 
01908 #ifdef HAVE_TPETRA_DEBUG
01909     // In a debug build, if the graph has a column Map, test whether
01910     // any of the given column indices are not in the column Map.
01911     // Keep track of the invalid column indices so we can tell the
01912     // user about them.
01913     if (hasColMap ()) {
01914       using Teuchos::Array;
01915       using Teuchos::toString;
01916       using std::endl;
01917       typedef LocalOrdinal LO;
01918       typedef typename ArrayView<const LO>::size_type size_type;
01919 
01920       const map_type& colMap = * (getColMap ());
01921       Array<LO> badColInds;
01922       bool allInColMap = true;
01923       for (size_type k = 0; k < indices.size (); ++k) {
01924         if (! colMap.isNodeLocalElement (indices[k])) {
01925           allInColMap = false;
01926           badColInds.push_back (indices[k]);
01927         }
01928       }
01929       if (! allInColMap) {
01930         std::ostringstream os;
01931         os << "Tpetra::CrsMatrix::insertLocalIndices: You attempted to insert "
01932           "entries in owned row " << localRow << ", at the following column "
01933           "indices: " << toString (indices) << "." << endl;
01934         os << "Of those, the following indices are not in the column Map on "
01935           "this process: " << toString (badColInds) << "." << endl << "Since "
01936           "the graph has a column Map already, it is invalid to insert entries "
01937           "at those locations.";
01938         TEUCHOS_TEST_FOR_EXCEPTION(! allInColMap, std::invalid_argument, os.str ());
01939       }
01940     }
01941 #endif // HAVE_TPETRA_DEBUG
01942 
01943     insertLocalIndicesImpl (localRow, indices);
01944 
01945 #ifdef HAVE_TPETRA_DEBUG
01946     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
01947       indicesAreAllocated() == false || isLocallyIndexed() == false,
01948       std::logic_error,
01949       ": Violated stated post-conditions. Please contact Tpetra team.");
01950 #endif
01951   }
01952 
01953 
01956   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
01957   void
01958   CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::
01959   insertLocalIndicesFiltered (const LocalOrdinal localRow,
01960                               const ArrayView<const LocalOrdinal> &indices)
01961   {
01962     typedef LocalOrdinal LO;
01963     const char tfecfFuncName[] = "insertLocalIndicesFiltered";
01964 
01965     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
01966       isFillActive() == false, std::runtime_error,
01967       ": requires that fill is active.");
01968     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
01969       isGloballyIndexed() == true, std::runtime_error,
01970       ": graph indices are global; use insertGlobalIndices().");
01971     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
01972       hasColMap() == false, std::runtime_error,
01973       ": cannot insert local indices without a column map.");
01974     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
01975       rowMap_->isNodeLocalElement(localRow) == false, std::runtime_error,
01976       ": row does not belong to this node.");
01977     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
01978       hasRowInfo() == false, std::runtime_error,
01979       ": graph row information was deleted at fillComplete().");
01980     if (indicesAreAllocated() == false) {
01981       allocateIndices(LocalIndices);
01982     }
01983 
01984      // If we have a column map, use it to filter the entries.
01985     if (hasColMap ()) {
01986       Array<LO> filtered_indices(indices);
01987       SLocalGlobalViews inds_view;
01988       SLocalGlobalNCViews inds_ncview;
01989       inds_ncview.linds = filtered_indices();
01990       const size_t numFilteredEntries =
01991         filterIndices<LocalIndices>(inds_ncview);
01992       inds_view.linds = filtered_indices (0, numFilteredEntries);
01993       insertLocalIndicesImpl(localRow, inds_view.linds);
01994     }
01995     else {
01996       insertLocalIndicesImpl(localRow, indices);
01997     }
01998 #ifdef HAVE_TPETRA_DEBUG
01999     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
02000       indicesAreAllocated() == false || isLocallyIndexed() == false,
02001       std::logic_error,
02002       ": Violated stated post-conditions. Please contact Tpetra team.");
02003 #endif
02004   }
02005 
02006 
02009   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
02010   void
02011   CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::
02012   insertGlobalIndices (const GlobalOrdinal grow,
02013                        const ArrayView<const GlobalOrdinal> &indices)
02014   {
02015     typedef LocalOrdinal LO;
02016     typedef GlobalOrdinal GO;
02017     typedef typename ArrayView<const GO>::size_type size_type;
02018     const char tfecfFuncName[] = "insertGlobalIndices";
02019 
02020     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
02021       isLocallyIndexed() == true, std::runtime_error,
02022       ": graph indices are local; use insertLocalIndices().");
02023     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
02024       hasRowInfo() == false, std::runtime_error,
02025       ": graph row information was deleted at fillComplete().");
02026     // This can't really be satisfied for now, because if we are
02027     // fillComplete(), then we are local.  In the future, this may
02028     // change.  However, the rule that modification require active
02029     // fill will not change.
02030     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
02031       isFillActive() == false, std::runtime_error,
02032       ": You are not allowed to call this method if fill is not active.  "
02033       "If fillComplete has been called, you must first call resumeFill "
02034       "before you may insert indices.");
02035     if (! indicesAreAllocated ()) {
02036       allocateIndices (GlobalIndices);
02037     }
02038     const LO myRow = rowMap_->getLocalElement (grow);
02039     if (myRow != Teuchos::OrdinalTraits<LO>::invalid ()) {
02040 #ifdef HAVE_TPETRA_DEBUG
02041       if (hasColMap ()) {
02042         using std::endl;
02043         const map_type& colMap = * (getColMap ());
02044         // In a debug build, keep track of the nonowned ("bad") column
02045         // indices, so that we can display them in the exception
02046         // message.  In a release build, just ditch the loop early if
02047         // we encounter a nonowned column index.
02048         Array<GO> badColInds;
02049         bool allInColMap = true;
02050         for (size_type k = 0; k < indices.size (); ++k) {
02051           if (! colMap.isNodeGlobalElement (indices[k])) {
02052             allInColMap = false;
02053             badColInds.push_back (indices[k]);
02054           }
02055         }
02056         if (! allInColMap) {
02057           std::ostringstream os;
02058           os << "Tpetra::CrsGraph::insertGlobalIndices: You attempted to insert "
02059             "entries in owned row " << grow << ", at the following column "
02060             "indices: " << toString (indices) << "." << endl;
02061           os << "Of those, the following indices are not in the column Map on "
02062             "this process: " << toString (badColInds) << "." << endl << "Since "
02063             "the matrix has a column Map already, it is invalid to insert "
02064             "entries at those locations.";
02065           TEUCHOS_TEST_FOR_EXCEPTION(! allInColMap, std::invalid_argument, os.str ());
02066         }
02067       }
02068 #endif // HAVE_TPETRA_DEBUG
02069       insertGlobalIndicesImpl (myRow, indices);
02070     }
02071     else { // a nonlocal row
02072       const size_type numIndices = indices.size ();
02073       // This creates the Array if it doesn't exist yet.  std::map's
02074       // operator[] does a lookup each time, so it's better to pull
02075       // nonlocals_[grow] out of the loop.
02076       std::vector<GO>& nonlocalRow = nonlocals_[grow];
02077       for (size_type k = 0; k < numIndices; ++k) {
02078         nonlocalRow.push_back (indices[k]);
02079       }
02080     }
02081 #ifdef HAVE_TPETRA_DEBUG
02082     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
02083       indicesAreAllocated() == false || isGloballyIndexed() == false,
02084       std::logic_error,
02085       ": Violated stated post-conditions. Please contact Tpetra team.");
02086 #endif
02087   }
02088 
02089 
02092   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
02093   void
02094   CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::
02095   insertGlobalIndicesFiltered (const GlobalOrdinal grow,
02096                                const ArrayView<const GlobalOrdinal> &indices)
02097   {
02098     typedef LocalOrdinal LO;
02099     typedef GlobalOrdinal GO;
02100     const char tfecfFuncName[] = "insertGlobalIndicesFiltered";
02101 
02102     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
02103       isLocallyIndexed() == true, std::runtime_error,
02104       ": graph indices are local; use insertLocalIndices().");
02105     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
02106       hasRowInfo() == false, std::runtime_error,
02107       ": graph row information was deleted at fillComplete().");
02108     // This can't really be satisfied for now, because if we are
02109     // fillComplete(), then we are local.  In the future, this may
02110     // change.  However, the rule that modification require active
02111     // fill will not change.
02112     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
02113       isFillActive() == false, std::runtime_error,
02114       ": You are not allowed to call this method if fill is not active.  "
02115       "If fillComplete has been called, you must first call resumeFill "
02116       "before you may insert indices.");
02117     if (indicesAreAllocated() == false) {
02118       allocateIndices (GlobalIndices);
02119     }
02120     const LO myRow = rowMap_->getLocalElement (grow);
02121     if (myRow != Teuchos::OrdinalTraits<LO>::invalid ()) {
02122       // If we have a column map, use it to filter the entries.
02123       if (hasColMap ()) {
02124         Array<GO> filtered_indices(indices);
02125         SLocalGlobalViews inds_view;
02126         SLocalGlobalNCViews inds_ncview;
02127         inds_ncview.ginds = filtered_indices();
02128         const size_t numFilteredEntries =
02129           filterIndices<GlobalIndices> (inds_ncview);
02130         inds_view.ginds = filtered_indices (0, numFilteredEntries);
02131         insertGlobalIndicesImpl(myRow, inds_view.ginds);
02132       }
02133       else {
02134        insertGlobalIndicesImpl(myRow, indices);
02135       }
02136     }
02137     else { // a nonlocal row
02138       typedef typename ArrayView<const GO>::size_type size_type;
02139       const size_type numIndices = indices.size ();
02140       for (size_type k = 0; k < numIndices; ++k) {
02141         nonlocals_[grow].push_back (indices[k]);
02142       }
02143     }
02144 #ifdef HAVE_TPETRA_DEBUG
02145     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
02146       indicesAreAllocated() == false || isGloballyIndexed() == false,
02147       std::logic_error,
02148       ": Violated stated post-conditions. Please contact Tpetra team.");
02149 #endif
02150   }
02151 
02154   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
02155   void CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::removeLocalIndices(LocalOrdinal lrow)
02156   {
02157     const char tfecfFuncName[] = "removeLocalIndices()";
02158     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC( isFillActive() == false,                    std::runtime_error, ": requires that fill is active.");
02159     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC( isStorageOptimized() == true,               std::runtime_error, ": cannot remove indices after optimizeStorage() has been called.");
02160     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC( isGloballyIndexed() == true,                std::runtime_error, ": graph indices are global.");
02161     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC( rowMap_->isNodeLocalElement(lrow) == false, std::runtime_error, ": row does not belong to this node.");
02162     if (indicesAreAllocated() == false) {
02163       allocateIndices(LocalIndices);
02164     }
02165     //
02166     clearGlobalConstants();
02167     //
02168     if (numRowEntries_ != null) {
02169       nodeNumEntries_ -= numRowEntries_[lrow];
02170       numRowEntries_[lrow] = 0;
02171     }
02172 #ifdef HAVE_TPETRA_DEBUG
02173     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(getNumEntriesInLocalRow(lrow) != 0 || indicesAreAllocated() == false || isLocallyIndexed() == false, std::logic_error,
02174         ": Violated stated post-conditions. Please contact Tpetra team.");
02175 #endif
02176   }
02177 
02180   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
02181   void CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::setAllIndices(const ArrayRCP<size_t> & rowPointers,const ArrayRCP<LocalOrdinal> & columnIndices)
02182   {
02183     const char tfecfFuncName[] = "setAllIndices()";
02184     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC( hasColMap() == false, std::runtime_error, ": requires a ColMap.");
02185     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC( (size_t)rowPointers.size() != getNodeNumRows()+1, std::runtime_error, ": requires rowPointers.size() == getNodeNumRows()+1");
02186     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC( lclInds1D_ != Teuchos::null || gblInds1D_ != Teuchos::null, std::runtime_error, ": cannot have 1D data structures allocated.");
02187     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC( lclInds2D_ != Teuchos::null || gblInds2D_ != Teuchos::null, std::runtime_error, ": cannot have 2D data structures allocated.");
02188 
02189     indicesAreAllocated_ = true;
02190     indicesAreLocal_     = true;
02191     pftype_              = StaticProfile; // if the profile wasn't static before, it sure is now.
02192     lclInds1D_           = columnIndices;
02193     rowPtrs_             = rowPointers;
02194     nodeNumEntries_ = nodeNumAllocated_ = rowPtrs_[getNodeNumRows()];
02195     numAllocForAllRows_  = 0;
02196     numAllocPerRow_      = null;
02197     checkInternalState();
02198   }
02199 
02200   // TODO: in the future, globalAssemble() should use import/export functionality
02203   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
02204   void CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::globalAssemble()
02205   {
02206     using Teuchos::Array;
02207     using Teuchos::as;
02208     using Teuchos::Comm;
02209     using Teuchos::gatherAll;
02210     using Teuchos::ireceive;
02211     using Teuchos::isend;
02212     using Teuchos::outArg;
02213     using Teuchos::REDUCE_MAX;
02214     using Teuchos::reduceAll;
02215     using Teuchos::toString;
02216     using Teuchos::waitAll;
02217     using std::endl;
02218     using std::make_pair;
02219     using std::pair;
02220     typedef GlobalOrdinal GO;
02221     typedef typename std::map<GO, std::vector<GO> >::const_iterator NLITER;
02222     typedef typename Array<GO>::size_type size_type;
02223 
02224     const char tfecfFuncName[] = "globalAssemble"; // for exception macro
02225     RCP<const Comm<int> > comm = getComm();
02226 
02227     const int numImages = comm->getSize();
02228     const int myImageID = comm->getRank();
02229 #ifdef HAVE_TPETRA_DEBUG
02230     Teuchos::barrier (*comm);
02231 #endif // HAVE_TPETRA_DEBUG
02232 
02233     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC( ! isFillActive(), std::runtime_error,
02234       ": requires that fill is active.");
02235     // Determine if any nodes have global entries to share
02236     {
02237       size_t MyNonlocals = nonlocals_.size(), MaxGlobalNonlocals;
02238       reduceAll (*comm, REDUCE_MAX, MyNonlocals, outArg (MaxGlobalNonlocals));
02239       if (MaxGlobalNonlocals == 0) {
02240         return;  // no entries to share
02241       }
02242     }
02243 
02244     // compute a list of NLRs from nonlocals_ and use it to compute:
02245     //      IdsAndRows: a vector of (id,row) pairs
02246     //          NLR2Id: a map from NLR to the Id that owns it
02247     // globalNeighbors: a global graph of connectivity between images: globalNeighbors(i,j) indicates that j sends to i
02248     //         sendIDs: a list of all images I send to
02249     //         recvIDs: a list of all images I receive from (constructed later)
02250     Array<pair<int, GO> > IdsAndRows;
02251     std::map<GO, int> NLR2Id;
02252     Teuchos::SerialDenseMatrix<int, char> globalNeighbors;
02253     Array<int> sendIDs, recvIDs;
02254     {
02255       // nonlocals_ contains the entries we are holding for all non-local rows
02256       // we want a list of the rows for which we have data
02257       Array<GO> NLRs;
02258       std::set<GO> setOfRows;
02259       for (NLITER iter = nonlocals_.begin(); iter != nonlocals_.end(); ++iter) {
02260         setOfRows.insert(iter->first);
02261       }
02262       // copy the elements in the set into an Array
02263       NLRs.resize(setOfRows.size());
02264       std::copy(setOfRows.begin(), setOfRows.end(), NLRs.begin());
02265 
02266       // get a list of ImageIDs for the non-local rows (NLRs)
02267       Array<int> NLRIds(NLRs.size());
02268       {
02269         LookupStatus stat = rowMap_->getRemoteIndexList(NLRs(),NLRIds());
02270         int lclerror = ( stat == IDNotPresent ? 1 : 0 );
02271         int gblerror;
02272         reduceAll<int, int> (*comm, REDUCE_MAX, lclerror, outArg (gblerror));
02273         if (gblerror != 0) {
02274           const int myRank = comm->getRank ();
02275           std::ostringstream os;
02276           os << "On one or more processes in the communicator, "
02277              << "there were insertions into rows of the graph that do not "
02278              << "exist in the row Map on any process in the communicator."
02279              << endl << "This process " << myRank << " is "
02280              << (lclerror == 0 ? "not " : "") << "one of those offending "
02281              << "processes." << endl;
02282           if (lclerror != 0) {
02283             // If NLRIds[k] is -1, then NLRs[k] is a row index not in
02284             // the row Map.  Collect this list of invalid row indices
02285             // for display in the exception message.
02286             Array<GO> invalidNonlocalRows;
02287             for (size_type k = 0; k < NLRs.size (); ++k) {
02288               if (NLRIds[k] == -1) {
02289                 invalidNonlocalRows.push_back (NLRs[k]);
02290               }
02291             }
02292             const size_type numInvalid = invalidNonlocalRows.size ();
02293             os << "On this process, " << numInvalid << " nonlocal row"
02294                << (numInvalid != 1 ? "s " : " ") << " were inserted that are "
02295                << "not in the row Map on any process." << endl;
02296             // Don't print _too_ many nonlocal rows.
02297             if (numInvalid <= 100) {
02298               os << "Offending row indices: "
02299                  << toString (invalidNonlocalRows ()) << endl;
02300             }
02301           }
02302           TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
02303             gblerror != 0, std::runtime_error,
02304             ": nonlocal entries correspond to invalid rows."
02305             << endl << os.str ());
02306         }
02307       }
02308 
02309       // build up a list of neighbors, as well as a map between NLRs and Ids
02310       // localNeighbors[i] != 0 iff I have data to send to image i
02311       // put NLRs,Ids into an array of pairs
02312       IdsAndRows.reserve(NLRs.size());
02313       Array<char> localNeighbors(numImages, 0);
02314       typename Array<GO>::const_iterator nlr;
02315       typename Array<int>::const_iterator id;
02316       for (nlr = NLRs.begin(), id = NLRIds.begin();
02317            nlr != NLRs.end(); ++nlr, ++id) {
02318         NLR2Id[*nlr] = *id;
02319         localNeighbors[*id] = 1;
02320         // IdsAndRows.push_back(make_pair<int,GlobalOrdinal>(*id,*nlr));
02321         IdsAndRows.push_back(make_pair(*id,*nlr));
02322       }
02323       for (int j=0; j<numImages; ++j) {
02324         if (localNeighbors[j]) {
02325           sendIDs.push_back(j);
02326         }
02327       }
02328       // sort IdsAndRows, by Ids first, then rows
02329       std::sort(IdsAndRows.begin(),IdsAndRows.end());
02330       // gather from other nodes to form the full graph
02331       globalNeighbors.shapeUninitialized(numImages,numImages);
02332       gatherAll (*getComm(), numImages, localNeighbors.getRawPtr(),
02333                  numImages * numImages, globalNeighbors.values());
02334       // globalNeighbors at this point contains (on all images) the
02335       // connectivity between the images.
02336       // globalNeighbors(i,j) != 0 means that j sends to i/that i receives from j
02337     }
02338 
02340     // FIGURE OUT WHO IS SENDING TO WHOM AND HOW MUCH
02341     // DO THIS IN THE PROCESS OF PACKING ALL OUTGOING DATA ACCORDING TO DESTINATION ID
02343 
02344     // loop over all columns to know from which images I can expect to receive something
02345     for (int j=0; j<numImages; ++j) {
02346       if (globalNeighbors(myImageID,j)) {
02347         recvIDs.push_back(j);
02348       }
02349     }
02350     const size_t numRecvs = recvIDs.size();
02351 
02352     // we know how many we're sending to already
02353     // form a contiguous list of all data to be sent
02354     // track the number of entries for each ID
02355     Array<pair<GO, GO> > IJSendBuffer;
02356     Array<size_t> sendSizes(sendIDs.size(), 0);
02357     size_t numSends = 0;
02358     for (typename Array<pair<int, GO> >::const_iterator IdAndRow = IdsAndRows.begin();
02359          IdAndRow != IdsAndRows.end(); ++IdAndRow) {
02360       int id = IdAndRow->first;
02361       GO row = IdAndRow->second;
02362       // have we advanced to a new send?
02363       if (sendIDs[numSends] != id) {
02364         numSends++;
02365         TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(sendIDs[numSends] != id,
02366           std::logic_error, ": internal logic error. Contact Tpetra team.");
02367       }
02368       // copy data for row into contiguous storage
02369       for (typename std::vector<GO>::const_iterator j = nonlocals_[row].begin(); j != nonlocals_[row].end(); ++j)
02370       {
02371         IJSendBuffer.push_back( pair<GlobalOrdinal,GlobalOrdinal>(row,*j) );
02372         sendSizes[numSends]++;
02373       }
02374     }
02375     if (IdsAndRows.size() > 0) {
02376       numSends++; // one last increment, to make it a count instead of an index
02377     }
02378     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(as<typename Array<int>::size_type>(numSends) != sendIDs.size(), std::logic_error, ": internal logic error. Contact Tpetra team.");
02379 
02380     // don't need this data anymore
02381     nonlocals_.clear();
02382 
02384     // TRANSMIT SIZE INFO BETWEEN SENDERS AND RECEIVERS
02386 
02387     // Array of pending nonblocking communication requests.  It's OK
02388     // to mix nonblocking send and receive requests in the same
02389     // waitAll() call.
02390     Array<RCP<Teuchos::CommRequest<int> > > requests;
02391 
02392     // perform non-blocking sends: send sizes to our recipients
02393     for (size_t s = 0; s < numSends ; ++s) {
02394       // We're using a nonowning RCP because all communication
02395       // will be local to this method and the scope of our data
02396       requests.push_back (isend<int, size_t> (*comm,
02397                                               rcp (&sendSizes[s], false),
02398                                               sendIDs[s]));
02399     }
02400     // perform non-blocking receives: receive sizes from our senders
02401     Array<size_t> recvSizes (numRecvs);
02402     for (size_t r = 0; r < numRecvs; ++r) {
02403       // We're using a nonowning RCP because all communication
02404       // will be local to this method and the scope of our data
02405       requests.push_back (ireceive (*comm, rcp (&recvSizes[r], false), recvIDs[r]));
02406     }
02407     // Wait on all the nonblocking sends and receives.
02408     if (! requests.empty()) {
02409       waitAll (*comm, requests());
02410     }
02411 #ifdef HAVE_TPETRA_DEBUG
02412     Teuchos::barrier (*comm);
02413 #endif // HAVE_TPETRA_DEBUG
02414 
02415     // This doesn't necessarily deallocate the array.
02416     requests.resize (0);
02417 
02419     // NOW SEND/RECEIVE ALL ROW DATA
02421     // from the size info, build the ArrayViews into IJSendBuffer
02422     Array<ArrayView<pair<GO,GO> > > sendBuffers(numSends,null);
02423     {
02424       size_t cur = 0;
02425       for (size_t s=0; s<numSends; ++s) {
02426         sendBuffers[s] = IJSendBuffer(cur,sendSizes[s]);
02427         cur += sendSizes[s];
02428       }
02429     }
02430     // perform non-blocking sends
02431     for (size_t s=0; s < numSends ; ++s)
02432     {
02433       // We're using a nonowning RCP because all communication
02434       // will be local to this method and the scope of our data
02435       ArrayRCP<pair<GO,GO> > tmpSendBuf =
02436         arcp (sendBuffers[s].getRawPtr(), 0, sendBuffers[s].size(), false);
02437       requests.push_back (isend<int, pair<GO,GO> > (*comm, tmpSendBuf, sendIDs[s]));
02438     }
02439     // calculate amount of storage needed for receives
02440     // setup pointers for the receives as well
02441     size_t totalRecvSize = std::accumulate (recvSizes.begin(), recvSizes.end(), 0);
02442     Array<pair<GO,GO> > IJRecvBuffer (totalRecvSize);
02443     // from the size info, build the ArrayViews into IJRecvBuffer
02444     Array<ArrayView<pair<GO,GO> > > recvBuffers (numRecvs, null);
02445     {
02446       size_t cur = 0;
02447       for (size_t r=0; r<numRecvs; ++r) {
02448         recvBuffers[r] = IJRecvBuffer(cur,recvSizes[r]);
02449         cur += recvSizes[r];
02450       }
02451     }
02452     // perform non-blocking recvs
02453     for (size_t r = 0; r < numRecvs; ++r) {
02454       // We're using a nonowning RCP because all communication
02455       // will be local to this method and the scope of our data
02456       ArrayRCP<pair<GO,GO> > tmpRecvBuf =
02457         arcp (recvBuffers[r].getRawPtr(), 0, recvBuffers[r].size(), false);
02458       requests.push_back (ireceive (*comm, tmpRecvBuf, recvIDs[r]));
02459     }
02460     // perform waits
02461     if (! requests.empty()) {
02462       waitAll (*comm, requests());
02463     }
02464 #ifdef HAVE_TPETRA_DEBUG
02465     Teuchos::barrier (*comm);
02466 #endif // HAVE_TPETRA_DEBUG
02467 
02469     // NOW PROCESS THE RECEIVED ROW DATA
02471     // TODO: instead of adding one entry at a time, add one row at a time.
02472     //       this requires resorting; they arrived sorted by sending node,
02473     //       so that entries could be non-contiguous if we received
02474     //       multiple entries for a particular row from different processors.
02475     //       it also requires restoring the data, which may make it not worth the trouble.
02476     for (typename Array<pair<GO,GO> >::const_iterator ij = IJRecvBuffer.begin();
02477          ij != IJRecvBuffer.end(); ++ij)
02478     {
02479       insertGlobalIndicesFiltered (ij->first, tuple<GO> (ij->second));
02480     }
02481     checkInternalState();
02482   }
02483 
02484 
02487   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
02488   void
02489   CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::
02490   resumeFill (const RCP<ParameterList> &params)
02491   {
02492     const char tfecfFuncName[] = "resumeFill";
02493     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(! hasRowInfo(), std::runtime_error,
02494       ": Sorry, you cannot resume fill of the CrsGraph, since the graph's row "
02495       "information was deleted in fillComplete().");
02496 
02497 #ifdef HAVE_TPETRA_DEBUG
02498     Teuchos::barrier( *rowMap_->getComm() );
02499 #endif // HAVE_TPETRA_DEBUG
02500     clearGlobalConstants();
02501     lclGraph_ = null;
02502     if (params != null) this->setParameterList (params);
02503     lowerTriangular_  = false;
02504     upperTriangular_  = false;
02505     // either still sorted/merged or initially sorted/merged
02506     indicesAreSorted_ = true;
02507     noRedundancies_ = true;
02508     fillComplete_ = false;
02509 #ifdef HAVE_TPETRA_DEBUG
02510     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
02511       ! isFillActive() || isFillComplete(), std::logic_error,
02512       "::resumeFill(): At end of method, either fill is not active or fill is "
02513       "complete.  This violates stated post-conditions.  Please report this bug "
02514       "to the Tpetra developers.");
02515 #endif // HAVE_TPETRA_DEBUG
02516   }
02517 
02518 
02521   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
02522   void CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::fillComplete(const RCP<ParameterList> &params)
02523   {
02524     fillComplete(rowMap_,rowMap_,params);
02525   }
02526 
02527 
02530   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
02531   void
02532   CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::
02533   fillComplete (const RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > &domainMap,
02534                 const RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > &rangeMap,
02535                 const RCP<ParameterList> &params)
02536   {
02537     const char tfecfFuncName[] = "fillComplete";
02538 
02539 #ifdef HAVE_TPETRA_DEBUG
02540     rowMap_->getComm ()->barrier ();
02541 #endif // HAVE_TPETRA_DEBUG
02542 
02543     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC( ! isFillActive() || isFillComplete(),
02544       std::runtime_error, ": Graph fill state must be active (isFillActive() "
02545       "must be true) before calling fillComplete().");
02546 
02547     const int numProcs = getComm ()->getSize ();
02548 
02549     // allocate if unallocated
02550     if (! indicesAreAllocated()) {
02551       if (hasColMap ()) {
02552         // We have a column Map, so use local indices.
02553         allocateIndices (LocalIndices);
02554       } else {
02555         // We don't have a column Map, so use global indices.
02556         allocateIndices (GlobalIndices);
02557       }
02558     }
02559 
02560     // If true, the caller promises that no process did nonlocal
02561     // changes since the last call to fillComplete.
02562     bool assertNoNonlocalInserts = false;
02563     if (! params.is_null ()) {
02564       assertNoNonlocalInserts =
02565         params->get<bool> ("No Nonlocal Changes", assertNoNonlocalInserts);
02566     }
02567     // We also don't need to do global assembly if there is only one
02568     // process in the communicator.
02569     const bool mayNeedGlobalAssemble = ! assertNoNonlocalInserts && numProcs > 1;
02570     if (mayNeedGlobalAssemble) {
02571       // This first checks if we need to do global assembly.
02572       // The check costs a single all-reduce.
02573       globalAssemble();
02574     }
02575     else {
02576       TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
02577         numProcs > 1 && nonlocals_.size() > 0, std::runtime_error,
02578         ":" << std::endl << "The graph's communicator contains only one "
02579         "process, but there are nonlocal entries.  " << std::endl <<
02580         "This probably means that invalid entries were added to the graph.");
02581     }
02582     // set domain/range map: may clear the import/export objects
02583     setDomainRangeMaps(domainMap,rangeMap);
02584 
02585     // If the graph does not already have a column Map (either from
02586     // the user constructor calling the version of the constructor
02587     // that takes a column Map, or from a previous fillComplete call),
02588     // then create it.
02589     if (! hasColMap ()) {
02590       makeColMap ();
02591     }
02592 
02593     // Make indices local, if they aren't already.
02594     // The method doesn't do any work if the indices are already local.
02595     makeIndicesLocal ();
02596 
02597     if (! isSorted()) {
02598       sortAllIndices();
02599     }
02600     if (! isMerged()) {
02601       mergeAllIndices();
02602     }
02603     makeImportExport(); // Make Import and Export objects
02604     computeGlobalConstants();
02605     // fill local objects
02606     fillLocalGraph(params);
02607     //
02608     fillComplete_ = true;
02609 #ifdef HAVE_TPETRA_DEBUG
02610     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC( isFillActive() == true || isFillComplete() == false, std::logic_error, ": Violated stated post-conditions. Please contact Tpetra team.");
02611 #endif
02612     //
02613     checkInternalState();
02614   }
02615 
02616 
02619   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
02620   void CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::expertStaticFillComplete(const RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > & domainMap,
02621                                                                                        const RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > & rangeMap,
02622                                                                                        const RCP<const Import<LocalOrdinal,GlobalOrdinal,Node> > &importer,
02623                                                                                        const RCP<const Export<LocalOrdinal,GlobalOrdinal,Node> > &exporter,
02624                                                                                        const RCP<ParameterList> &params)
02625   {
02626 #ifdef HAVE_TPETRA_DEBUG
02627     rowMap_->getComm ()->barrier ();
02628 #endif // HAVE_TPETRA_DEBUG
02629     const char tfecfFuncName[] = "expertStaticFillComplete()";
02630     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC( isFillComplete() == true || hasColMap() == false,
02631                                            std::runtime_error, ": fillComplete cannot have already been called and a ColMap is required.");
02632 
02633     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC( getNodeNumRows() > 0 && rowPtrs_==Teuchos::null,
02634                                            std::runtime_error, ": a matrix will getNodeNumRows()>0 requires rowptr to be set.");
02635 
02636     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC( domainMap == Teuchos::null || rangeMap == Teuchos::null,
02637                                            std::runtime_error, ": requires a non-null domainMap & rangeMap.");
02638 
02639     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC( pftype_ !=StaticProfile,
02640                                            std::runtime_error, ": requires StaticProfile.");
02641 
02642     // Note: We don't need to do the following things which are normally done in fillComplete:
02643     // allocateIndices, globalAssemble, makeColMap, makeIndicesLocal, sortAllIndices, mergeAllIndices
02644 
02645     // Note: Need to do this so computeGlobalConstants & fillLocalGraph work
02646     nodeNumEntries_ = nodeNumAllocated_ = rowPtrs_[getNodeNumRows()];
02647 
02648     // Constants from allocateIndices
02649     numAllocForAllRows_  = 0;
02650     numAllocPerRow_      = null;
02651     indicesAreAllocated_ = true;
02652 
02653     // Constants from makeIndicesLocal
02654     indicesAreLocal_  = true;
02655     indicesAreGlobal_ = false;
02656 
02657     // set domain/range map: may clear the import/export objects
02658     setDomainRangeMaps(domainMap,rangeMap);
02659 
02660     // Presume the user sorted and merged the arrays first
02661     indicesAreSorted_ = true;
02662     noRedundancies_ = true;
02663 
02664     // makeImportExport won't create a new importer/exporter if I set one here first.
02665     importer_=Teuchos::null;
02666     exporter_=Teuchos::null;
02667     if(importer != Teuchos::null) {
02668       TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(!importer->getSourceMap()->isSameAs(*getDomainMap()) || !importer->getTargetMap()->isSameAs(*getColMap()),
02669                                             std::invalid_argument,": importer does not match matrix maps.");
02670       importer_ = importer;
02671 
02672     }
02673     if(exporter != Teuchos::null) {
02674       TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(!exporter->getSourceMap()->isSameAs(*getRowMap()) || !exporter->getTargetMap()->isSameAs(*getRangeMap()),
02675                                             std::invalid_argument,": exporter does not match matrix maps.");
02676       exporter_ = exporter;
02677     }
02678     makeImportExport();
02679 
02680     // Compute the constants
02681     computeGlobalConstants();
02682 
02683     // Since we have a StaticProfile, fillLocalGraph will do the right thing...
02684     fillLocalGraph(params);
02685     fillComplete_ = true;
02686 
02687 #ifdef HAVE_TPETRA_DEBUG
02688     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC( isFillActive() == true || isFillComplete() == false, std::logic_error, ": Violated stated post-conditions. Please contact Tpetra team.");
02689 #endif
02690     checkInternalState();
02691   }
02692 
02693 
02694 
02697   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
02698   void CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::fillLocalGraph(const RCP<ParameterList> &params)
02699   {
02700     const size_t numRows = getNodeNumRows();
02701     ArrayRCP<size_t> ptrs;
02702     ArrayRCP<LocalOrdinal> inds;
02703     bool requestOptimizedStorage = true;
02704     if (params != null && params->get("Optimize Storage",true) == false) requestOptimizedStorage = false;
02705     if (getProfileType() == DynamicProfile) {
02706       // 2d -> 1d packed
02707       ptrs = LocalMatOps::allocRowPtrs( getRowMap()->getNode(), numRowEntries_() );
02708       inds = LocalMatOps::template allocStorage<LocalOrdinal>( getRowMap()->getNode(), ptrs() );
02709       for (size_t row=0; row < numRows; ++row) {
02710         const size_t numentrs = numRowEntries_[row];
02711         std::copy( lclInds2D_[row].begin(), lclInds2D_[row].begin()+numentrs, inds+ptrs[row] );
02712       }
02713     }
02714     else if (getProfileType() == StaticProfile) {
02715       // 1d non-packed -> 1d packed
02716       if (nodeNumEntries_ != nodeNumAllocated_) {
02717         ptrs = LocalMatOps::allocRowPtrs( getRowMap()->getNode(), numRowEntries_() );
02718         inds = LocalMatOps::template allocStorage<LocalOrdinal>( getRowMap()->getNode(), ptrs() );
02719         for (size_t row=0; row < numRows; ++row) {
02720           const size_t numentrs = numRowEntries_[row];
02721           std::copy( lclInds1D_+rowPtrs_[row], lclInds1D_+rowPtrs_[row]+numentrs, inds+ptrs[row] );
02722         }
02723       }
02724       else {
02725         inds = lclInds1D_;
02726         ptrs = rowPtrs_;
02727       }
02728     }
02729     // can we ditch the old allocations for the packed one?
02730     if ( requestOptimizedStorage ) {
02731       lclInds2D_ = null;
02732       numRowEntries_ = null;
02733       // keep the new stuff
02734       lclInds1D_ = inds;
02735       rowPtrs_ = ptrs;
02736       nodeNumAllocated_ = nodeNumEntries_;
02737       pftype_ = StaticProfile;
02738     }
02739     // build the local graph, hand over the indices
02740     RCP<ParameterList> lclparams;
02741     if (params == null) lclparams = parameterList();
02742     else                lclparams = sublist(params,"Local Graph");
02743     lclGraph_ = rcp( new local_graph_type( getRowMap()->getNodeNumElements(), getColMap()->getNodeNumElements(), getRowMap()->getNode(), lclparams ) );
02744     lclGraph_->setStructure(ptrs,inds);
02745     ptrs = null;
02746     inds = null;
02747     // finalize local graph
02748     Teuchos::EDiag diag = ( getNodeNumDiags() < getNodeNumRows() ? Teuchos::UNIT_DIAG : Teuchos::NON_UNIT_DIAG );
02749     Teuchos::EUplo uplo = Teuchos::UNDEF_TRI;
02750     if      (isUpperTriangular()) uplo = Teuchos::UPPER_TRI;
02751     else if (isLowerTriangular()) uplo = Teuchos::LOWER_TRI;
02752     LocalMatOps::finalizeGraph(uplo,diag,*lclGraph_,params);
02753   }
02754 
02755 
02758   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
02759   void
02760   CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::
02761   replaceColMap (const Teuchos::RCP<const map_type>& newColMap)
02762   {
02763     // NOTE: This safety check matches the code, but not the documentation of Crsgraph
02764     const char tfecfFuncName[] = "replaceColMap";
02765     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(isLocallyIndexed() || isGloballyIndexed(),  std::runtime_error, " requires matching maps and non-static graph.");
02766     colMap_ = newColMap;
02767   }
02768 
02769 
02772   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
02773   void
02774   CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::
02775   replaceDomainMapAndImporter (const Teuchos::RCP<const map_type>& newDomainMap,
02776                                const Teuchos::RCP<const import_type>& newImporter)
02777   {
02778     const char prefix[] = "Tpetra::CrsGraph::replaceDomainMapAndImporter: ";
02779     TEUCHOS_TEST_FOR_EXCEPTION(
02780       colMap_.is_null (), std::invalid_argument, prefix << "You may not call "
02781       "this method unless the graph already has a column Map.");
02782     TEUCHOS_TEST_FOR_EXCEPTION(
02783       newDomainMap.is_null (), std::invalid_argument,
02784       prefix << "The new domain Map must be nonnull.");
02785 
02786 #ifdef HAVE_TPETRA_DEBUG
02787     if (newImporter.is_null ()) {
02788       // It's not a good idea to put expensive operations in a macro
02789       // clause, even if they are side effect - free, because macros
02790       // don't promise that they won't evaluate their arguments more
02791       // than once.  It's polite for them to do so, but not required.
02792       const bool colSameAsDom = colMap_->isSameAs (*newDomainMap);
02793       TEUCHOS_TEST_FOR_EXCEPTION(
02794         colSameAsDom, std::invalid_argument, "If the new Import is null, "
02795         "then the new domain Map must be the same as the current column Map.");
02796     }
02797     else {
02798       const bool colSameAsTgt =
02799         colMap_->isSameAs (* (newImporter->getTargetMap ()));
02800       const bool newDomSameAsSrc =
02801         newDomainMap->isSameAs (* (newImporter->getSourceMap ()));
02802       TEUCHOS_TEST_FOR_EXCEPTION(
02803         ! colSameAsTgt || ! newDomSameAsSrc, std::invalid_argument, "If the "
02804         "new Import is nonnull, then the current column Map must be the same "
02805         "as the new Import's target Map, and the new domain Map must be the "
02806         "same as the new Import's source Map.");
02807     }
02808 #endif // HAVE_TPETRA_DEBUG
02809 
02810     domainMap_ = newDomainMap;
02811     importer_ = Teuchos::rcp_const_cast<import_type> (newImporter);
02812   }
02813 
02814 
02817   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
02818   const RCP<const typename LocalMatOps::template graph<LocalOrdinal,Node>::graph_type>
02819   CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getLocalGraph() const
02820   {
02821     return lclGraph_;
02822   }
02823 
02824 
02827   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
02828   const RCP<typename LocalMatOps::template graph<LocalOrdinal,Node>::graph_type>
02829   CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getLocalGraphNonConst()
02830   {
02831     return lclGraph_;
02832   }
02833 
02834 
02837   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
02838   void CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::computeGlobalConstants()
02839   {
02840     using Teuchos::as;
02841     using Teuchos::outArg;
02842     using Teuchos::reduceAll;
02843     typedef LocalOrdinal LO;
02844     typedef GlobalOrdinal GO;
02845     typedef global_size_t GST;
02846 
02847 #ifdef HAVE_TPETRA_DEBUG
02848     TEUCHOS_TEST_FOR_EXCEPTION(! hasColMap(), std::logic_error, "Tpetra::"
02849       "CrsGraph::computeGlobalConstants: At this point, the graph should have "
02850       "a column Map, but it does not.  Please report this bug to the Tpetra "
02851       "developers.");
02852 #endif // HAVE_TPETRA_DEBUG
02853 
02854     // If necessary, (re)compute the local constants: nodeNumDiags_,
02855     // lowerTriangular_, upperTriangular_, and nodeMaxNumRowEntries_.
02856     if (! haveLocalConstants_) {
02857       // We have actually already computed nodeNumEntries_.
02858       // nodeNumEntries_ gets updated by methods that insert or remove
02859       // indices (including setAllIndices and
02860       // expertStaticFillComplete).  Before fillComplete, its count
02861       // may include duplicate column indices in the same row.
02862       // However, mergeRowIndices and mergeRowIndicesAndValues both
02863       // subtract off merged indices in each row from the total count.
02864       // Thus, nodeNumEntries_ _should_ be accurate at this point,
02865       // meaning that we don't have to re-count it here.
02866 
02867       // Reset local properties
02868       upperTriangular_ = true;
02869       lowerTriangular_ = true;
02870       nodeMaxNumRowEntries_ = 0;
02871       nodeNumDiags_         = 0;
02872 
02873       // At this point, we know that we have both a row Map and a column Map.
02874       const Map<LO,GO,Node>& rowMap = *rowMap_;
02875       const Map<LO,GO,Node>& colMap = *colMap_;
02876 
02877       // Go through all the entries of the graph.  Count the number of
02878       // diagonal elements we encounter, and figure out whether the
02879       // graph is lower or upper triangular.  Diagonal elements are
02880       // determined using global indices, with respect to the whole
02881       // graph.  However, lower or upper triangularity is a local
02882       // property, and is determined using local indices.
02883       //
02884       // At this point, indices have already been sorted in each row.
02885       // That makes finding out whether the graph is lower / upper
02886       // triangular easier.
02887       if (indicesAreAllocated () && nodeNumAllocated_ > 0) {
02888         const size_t numLocalRows = getNodeNumRows ();
02889         for (size_t localRow = 0; localRow < numLocalRows; ++localRow) {
02890           const GO globalRow = rowMap.getGlobalElement (localRow);
02891           // Find the local (column) index for the diagonal element.
02892           const LO rlcid = colMap.getLocalElement (globalRow);
02893           RowInfo rowInfo = getRowInfo (localRow);
02894           ArrayView<const LO> rview = getLocalView (rowInfo);
02895           typename ArrayView<const LO>::iterator beg, end, cur;
02896           beg = rview.begin();
02897           end = beg + rowInfo.numEntries;
02898           if (beg != end) {
02899             for (cur = beg; cur != end; ++cur) {
02900               // is this the diagonal?
02901               if (rlcid == *cur) ++nodeNumDiags_;
02902             }
02903             // Local column indices are sorted in each row.  That means
02904             // the smallest column index in this row (on this process)
02905             // is *beg, and the largest column index in this row (on
02906             // this process) is *(end - 1).  We know that end - 1 is
02907             // valid because beg != end.
02908             const size_t smallestCol = as<size_t> (*beg);
02909             const size_t largestCol = as<size_t> (*(end - 1));
02910 
02911             if (smallestCol < localRow) {
02912               upperTriangular_ = false;
02913             }
02914             if (localRow < largestCol) {
02915               lowerTriangular_ = false;
02916             }
02917           }
02918           // Update the max number of entries over all rows.
02919           nodeMaxNumRowEntries_ = std::max (nodeMaxNumRowEntries_, rowInfo.numEntries);
02920         }
02921       }
02922       haveLocalConstants_ = true;
02923     } // if my process doesn't have local constants
02924 
02925     // Compute global constants from local constants.  Processes that
02926     // already have local constants still participate in the
02927     // all-reduces, using their previously computed values.
02928     if (haveGlobalConstants_ == false) {
02929       // Promote all the nodeNum* and nodeMaxNum* quantities from
02930       // size_t to global_size_t, when doing the all-reduces for
02931       // globalNum* / globalMaxNum* results.
02932       //
02933       // FIXME (mfh 07 May 2013) Unfortunately, we either have to do
02934       // this in two all-reduces (one for the sum and the other for
02935       // the max), or use a custom MPI_Op that combines the sum and
02936       // the max.  The latter might even be slower than two
02937       // all-reduces on modern network hardware.  It would also be a
02938       // good idea to use nonblocking all-reduces (MPI 3), so that we
02939       // don't have to wait around for the first one to finish before
02940       // starting the second one.
02941       GST lcl[2], gbl[2];
02942       lcl[0] = as<GST> (nodeNumEntries_);
02943       lcl[1] = as<GST> (nodeNumDiags_);
02944       reduceAll<int,GST> (*getComm (), Teuchos::REDUCE_SUM,
02945                           2, lcl, gbl);
02946       globalNumEntries_ = gbl[0];
02947       globalNumDiags_   = gbl[1];
02948       reduceAll<int,GST> (*getComm (), Teuchos::REDUCE_MAX,
02949                           as<GST> (nodeMaxNumRowEntries_),
02950                           outArg (globalMaxNumRowEntries_));
02951       haveGlobalConstants_ = true;
02952     }
02953   }
02954 
02955 
02958   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
02959   void
02960   CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::makeIndicesLocal ()
02961   {
02962     typedef LocalOrdinal LO;
02963     typedef GlobalOrdinal GO;
02964 #ifdef HAVE_TPETRA_DEBUG
02965     const char tfecfFuncName[] = "makeIndicesLocal";
02966 #endif // HAVE_TPETRA_DEBUG
02967 
02968     // Transform indices to local index space
02969     const size_t nlrs = getNodeNumRows ();
02970     if (isGloballyIndexed () && nlrs > 0) {
02971       // allocate data for local indices
02972       if (getProfileType() == StaticProfile) {
02973         // If GO and LO are the same size, we can reuse the existing
02974         // array of 1-D index storage to convert column indices from
02975         // GO to LO.  Otherwise, we'll just allocate a new buffer.
02976         if (nodeNumAllocated_ && sizeof (LO) == sizeof (GO)) {
02977           lclInds1D_ = arcp_reinterpret_cast<LO> (gblInds1D_).persistingView (0, nodeNumAllocated_);
02978         }
02979         else {
02980           lclInds1D_ = LocalMatOps::template allocStorage<LO> (getRowMap ()->getNode (), rowPtrs_ ());
02981         }
02982         for (size_t r = 0; r < nlrs; ++r) {
02983           const size_t offset   = rowPtrs_[r];
02984           const size_t numentry = numRowEntries_[r];
02985           for (size_t j = 0; j < numentry; ++j) {
02986             const GO gid = gblInds1D_[offset + j];
02987             const LO lid = colMap_->getLocalElement (gid);
02988             lclInds1D_[offset + j] = lid;
02989 #ifdef HAVE_TPETRA_DEBUG
02990             TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
02991               lclInds1D_[offset + j] == Teuchos::OrdinalTraits<LO>::invalid(),
02992               std::logic_error,
02993               ": In local row r=" << r << ", global column " << gid << " is "
02994               "not in the column Map.  This should never happen.  Please "
02995               "report this bug to the Tpetra developers.");
02996 #endif // HAVE_TPETRA_DEBUG
02997           }
02998         }
02999         // We've converted column indices from global to local, so we
03000         // can deallocate the global column indices (which we know are
03001         // in 1-D storage, because the graph has static profile).
03002         gblInds1D_ = null;
03003       }
03004       else {  // the graph has dynamic profile (2-D index storage)
03005         lclInds2D_ = arcp<Array<LO> > (nlrs);
03006         for (size_t r = 0; r < getNodeNumRows (); ++r) {
03007           if (! gblInds2D_[r].empty ()) {
03008             const GO* const ginds = gblInds2D_[r].getRawPtr ();
03009             const size_t rna = gblInds2D_[r].size ();
03010             const size_t numentry = numRowEntries_[r];
03011             lclInds2D_[r].resize (rna);
03012             LO* const linds = lclInds2D_[r].getRawPtr ();
03013             for (size_t j = 0; j < numentry; ++j) {
03014               GO gid = ginds[j];
03015               LO lid = colMap_->getLocalElement (gid);
03016               linds[j] = lid;
03017 #ifdef HAVE_TPETRA_DEBUG
03018               TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
03019                 linds[j] == Teuchos::OrdinalTraits<LO>::invalid(),
03020                 std::logic_error,
03021                 ": Global column ginds[j=" << j << "]=" << ginds[j]
03022                 << " of local row r=" << r << " is not in the column Map.  "
03023                 "This should never happen.  Please report this bug to the "
03024                 "Tpetra developers.");
03025 #endif // HAVE_TPETRA_DEBUG
03026             }
03027           }
03028         }
03029         gblInds2D_ = null;
03030       }
03031     } // globallyIndexed() && nlrs > 0
03032     indicesAreLocal_  = true;
03033     indicesAreGlobal_ = false;
03034     checkInternalState();
03035   }
03036 
03037 
03040   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
03041   void CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::sortAllIndices()
03042   {
03043     TEUCHOS_TEST_FOR_EXCEPT(isGloballyIndexed()==true);   // this should be called only after makeIndicesLocal()
03044     if (isSorted() == false) {
03045       for (size_t row=0; row < getNodeNumRows(); ++row) {
03046         RowInfo rowInfo = getRowInfo(row);
03047         sortRowIndices(rowInfo);
03048       }
03049     }
03050     // we just sorted every row
03051     indicesAreSorted_ = true;
03052   }
03053 
03054 
03057   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
03058   void CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::makeColMap()
03059   {
03060     using std::endl;
03061     using Teuchos::REDUCE_MAX;
03062     using Teuchos::reduceAll;
03063     typedef LocalOrdinal LO;
03064     typedef GlobalOrdinal GO;
03065     const char tfecfFuncName[] = "makeColMap";
03066 
03067     if (hasColMap ()) { // The graph already has a column Map.
03068       // FIXME (mfh 26 Feb 2013): This currently prevents structure
03069       // changes that affect the column Map.
03070       return;
03071     }
03072 
03073     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
03074       isLocallyIndexed (), std::runtime_error,
03075       ": The graph is locally indexed.  Calling makeColMap() to make the "
03076       "column Map requires that the graph be globally indexed.");
03077 
03078     // After the calling process is done going through all of the rows
03079     // it owns, myColumns will contain the list of indices owned by
03080     // this process in the column Map.
03081     Array<GO> myColumns;
03082 
03083     // If we reach this point, we don't have a column Map yet, so the
03084     // graph can't be locally indexed.  Thus, isGloballyIndexed() ==
03085     // false means that the graph is empty on this process, so
03086     // myColumns will be left empty.
03087     if (isGloballyIndexed ()) {
03088       // Go through all the rows, finding the populated column indices.
03089       //
03090       // Our final list of indices for the column Map constructor will
03091       // have the following properties (all of which are with respect
03092       // to the calling process):
03093       //
03094       // 1. Indices in the domain Map go first.
03095       // 2. Indices not in the domain Map follow, ordered first
03096       //    contiguously by their owning process rank (in the domain
03097       //    Map), then in increasing order within that.
03098       // 3. No duplicate indices.
03099       //
03100       // This imitates the ordering used by Aztec(OO) and Epetra.
03101       // Storing indices owned by the same process (in the domain Map)
03102       // contiguously permits the use of contiguous send and receive
03103       // buffers.
03104       //
03105       // We begin by partitioning the column indices into "local" GIDs
03106       // (owned by the domain Map) and "remote" GIDs (not owned by the
03107       // domain Map).  We use the same order for local GIDs as the
03108       // domain Map, so we track them in place in their array.  We use
03109       // an std::set (RemoteGIDSet) to keep track of remote GIDs, so
03110       // that we don't have to merge duplicates later.
03111       const LO LINV = Teuchos::OrdinalTraits<LO>::invalid ();
03112       size_t numLocalColGIDs = 0, numRemoteColGIDs = 0;
03113 
03114       // GIDisLocal[lid] == 0 if and only if local index lid in the
03115       // domain Map is remote (not local).
03116       Array<char> GIDisLocal (domainMap_->getNodeNumElements (), 0);
03117       std::set<GO> RemoteGIDSet;
03118       const size_t myNumRows = getNodeNumRows ();
03119       for (size_t r = 0; r < myNumRows; ++r) {
03120         RowInfo rowinfo = getRowInfo (r);
03121         if (rowinfo.numEntries > 0) {
03122           // FIXME (mfh 03 Mar 2013) It's a bit puzzling to me why the
03123           // ArrayView that getGlobalView() returns doesn't return
03124           // rowinfo.numEntries entries.
03125           ArrayView<const GO> rowGids = getGlobalView (rowinfo);
03126           rowGids = rowGids (0, rowinfo.numEntries);
03127 
03128           for (size_t k = 0; k < rowinfo.numEntries; ++k) {
03129             const GO gid = rowGids[k];
03130             const LO lid = domainMap_->getLocalElement (gid);
03131             if (lid != LINV) {
03132               const char alreadyFound = GIDisLocal[lid];
03133               if (alreadyFound == 0) {
03134                 GIDisLocal[lid] = 1;
03135                 ++numLocalColGIDs;
03136               }
03137             }
03138             else {
03139               const bool notAlreadyFound = RemoteGIDSet.insert (gid).second;
03140               if (notAlreadyFound) { // gid did not exist in the set before
03141                 ++numRemoteColGIDs;
03142               }
03143             }
03144           } // for each entry k in row r
03145         } // if row r contains > 0 entries
03146       } // for each locally owned row r
03147 
03148       // Possible short-circuit for serial scenario:
03149       //
03150       // If all domain GIDs are present as column indices, then set
03151       // ColMap=DomainMap.  By construction, LocalGIDs is a subset of
03152       // DomainGIDs.
03153       //
03154       // If we have
03155       //   * Number of remote GIDs is 0, so that ColGIDs == LocalGIDs,
03156       // and
03157       //   * Number of local GIDs is number of domain GIDs
03158       // then
03159       //   * LocalGIDs \subset DomainGIDs && size(LocalGIDs) ==
03160       //     size(DomainGIDs) => DomainGIDs == LocalGIDs == ColGIDs
03161       // on the calling process.
03162       //
03163       // We will concern ourselves only with the special case of a
03164       // serial DomainMap, obviating the need for communication.
03165       //
03166       // If
03167       //   * DomainMap has a serial communicator
03168       // then we can set the column Map as the domain Map
03169       // return. Benefit: this graph won't need an Import object
03170       // later.
03171       //
03172       // Note, for a serial domain map, there can be no RemoteGIDs,
03173       // because there are no remote processes.  Likely explanations
03174       // for this are:
03175       //  * user submitted erroneous column indices
03176       //  * user submitted erroneous domain Map
03177       if (domainMap_->getComm ()->getSize () == 1) {
03178         TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
03179           numRemoteColGIDs != 0, std::runtime_error,
03180           ": " << numRemoteColGIDs << " column "
03181           << (numRemoteColGIDs != 1 ? "indices are" : "index is")
03182           << " not in the domain Map." << endl
03183           << "Either these indices are invalid or the domain Map is invalid."
03184           << endl << "Remember that nonsquare matrices, or matrices where the "
03185           "row and range Maps are different, require calling the version of "
03186           "fillComplete that takes the domain and range Maps as input.");
03187         if (numLocalColGIDs == domainMap_->getNodeNumElements()) {
03188           colMap_ = domainMap_;
03189           checkInternalState ();
03190           return;
03191         }
03192       }
03193 
03194       // Populate myColumns with a list of all column GIDs.  Put
03195       // locally owned (in the domain Map) GIDs at the front: they
03196       // correspond to "same" and "permuted" entries between the
03197       // column Map and the domain Map.  Put remote GIDs at the back.
03198       myColumns.resize (numLocalColGIDs + numRemoteColGIDs);
03199       // get pointers into myColumns for each part
03200       ArrayView<GO> LocalColGIDs  = myColumns (0, numLocalColGIDs);
03201       ArrayView<GO> RemoteColGIDs = myColumns (numLocalColGIDs, numRemoteColGIDs);
03202 
03203       // Copy the remote GIDs into myColumns
03204       std::copy (RemoteGIDSet.begin(), RemoteGIDSet.end(), RemoteColGIDs.begin());
03205 
03206       // Make a list of process ranks corresponding to the remote GIDs.
03207       Array<int> RemoteImageIDs (numRemoteColGIDs);
03208       // Look up the remote process' ranks in the domain Map.
03209       {
03210         const LookupStatus stat =
03211           domainMap_->getRemoteIndexList (RemoteColGIDs, RemoteImageIDs ());
03212 #ifdef HAVE_TPETRA_DEBUG
03213         // If any process returns IDNotPresent, then at least one of
03214         // the remote indices was not present in the domain Map.  This
03215         // means that the Import object cannot be constructed, because
03216         // of incongruity between the column Map and domain Map.
03217         // This has two likely causes:
03218         //   - The user has made a mistake in the column indices
03219         //   - The user has made a mistake with respect to the domain Map
03220         const int missingID_lcl = (stat == IDNotPresent ? 1 : 0);
03221         int missingID_gbl = 0;
03222         reduceAll<int, int> (*getComm (), REDUCE_MAX, missingID_lcl,
03223                              outArg (missingID_gbl));
03224         TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
03225           missingID_gbl == 1, std::runtime_error,
03226           ": Some column indices are not in the domain Map." << endl
03227           << "Either these column indices are invalid or the domain Map is "
03228           "invalid." << endl << "Likely cause: For a nonsquare matrix, you "
03229           "must give the domain and range Maps as input to fillComplete.");
03230 #else
03231         (void) stat; // forestall compiler warning for unused variable
03232 #endif // HAVE_TPETRA_DEBUG
03233       }
03234       // Sort incoming remote column indices so that all columns
03235       // coming from a given remote process are contiguous.  This
03236       // means the Import's Distributor doesn't need to reorder data.
03237       sort2 (RemoteImageIDs.begin(), RemoteImageIDs.end(), RemoteColGIDs.begin());
03238 
03239       // Copy the local GIDs into myColumns. Two cases:
03240       // 1. If the number of Local column GIDs is the same as the
03241       //    number of Local domain GIDs, we can simply read the domain
03242       //    GIDs into the front part of ColIndices (see logic above
03243       //    from the serial short circuit case)
03244       // 2. We step through the GIDs of the DomainMap, checking to see
03245       //    if each domain GID is a column GID.  We want to do this to
03246       //    maintain a consistent ordering of GIDs between the columns
03247       //    and the domain.
03248 
03249       // FIXME (mfh 03 Mar 2013) It's common that the domain Map is
03250       // contiguous.  It would be more efficient in that case to avoid
03251       // calling getNodeElementList(), since that permanently
03252       // constructs and caches the GID list in the contiguous Map.
03253       ArrayView<const GO> domainElts = domainMap_->getNodeElementList ();
03254       const size_t numDomainElts = domainMap_->getNodeNumElements ();
03255       if (numLocalColGIDs == numDomainElts) {
03256         // If the number of locally owned GIDs are the same as the
03257         // number of local domain Map elements, then the local domain
03258         // Map elements are the same as the locally owned GIDs.
03259         std::copy (domainElts.begin(), domainElts.end(), LocalColGIDs.begin());
03260       }
03261       else {
03262         // Count the number of locally owned GIDs, both to keep track
03263         // of the current array index, and as a sanity check.
03264         size_t numLocalCount = 0;
03265         for (size_t i = 0; i < numDomainElts; ++i) {
03266           if (GIDisLocal[i]) {
03267             LocalColGIDs[numLocalCount++] = domainElts[i];
03268           }
03269         }
03270         TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
03271           numLocalCount != numLocalColGIDs, std::logic_error,
03272           ": numLocalCount = " << numLocalCount << " != numLocalColGIDs = "
03273           << numLocalColGIDs << ".  This should never happen.  Please report "
03274           "this bug to the Tpetra developers.");
03275       }
03276 
03277       // FIXME (mfh 03 Apr 2013) Now would be a good time to use the
03278       // information we collected above to construct the Import.  In
03279       // particular, building an Import requires:
03280       //
03281       // 1. numSameIDs (length of initial contiguous sequence of GIDs
03282       //    on this process that are the same in both Maps; this
03283       //    equals the number of domain Map elements on this process)
03284       //
03285       // 2. permuteToLIDs and permuteFromLIDs (both empty in this
03286       //    case, since there's no permutation going on; the column
03287       //    Map starts with the domain Map's GIDs, and immediately
03288       //    after them come the remote GIDs)
03289       //
03290       // 3. remoteGIDs (exactly those GIDs that we found out above
03291       //    were not in the domain Map) and remoteLIDs (which we could
03292       //    have gotten above by using the three-argument version of
03293       //    getRemoteIndexList() that computes local indices as well
03294       //    as process ranks, instead of the two-argument version that
03295       //    was used above)
03296       //
03297       // 4. remotePIDs (which we have from the getRemoteIndexList()
03298       //    call above)
03299       //
03300       // 5. Sorting remotePIDs, and applying that permutation to
03301       //    remoteGIDs and remoteLIDs (by calling sort3 above instead
03302       //    of sort2)
03303       //
03304       // 6. Everything after the sort3 call in Import::setupExport():
03305       //    a. Create the Distributor via createFromRecvs(), which
03306       //       computes exportGIDs and exportPIDs
03307       //    b. Compute exportLIDs from exportGIDs (by asking the
03308       //       source Map, in this case the domain Map, to convert
03309       //       global to local)
03310       //
03311       // Steps 1-5 come for free, since we must do that work anyway in
03312       // order to compute the column Map.  In particular, Step 3 is
03313       // even more expensive than Step 6a, since it involves both
03314       // creating and using a new Distributor object.
03315 
03316     } // if the graph is globally indexed
03317 
03318     const global_size_t gstInv =
03319       Teuchos::OrdinalTraits<global_size_t>::invalid ();
03320     const GO indexBase = domainMap_->getIndexBase ();
03321     colMap_ = rcp (new map_type (gstInv, myColumns, indexBase,
03322                                  domainMap_->getComm (),
03323                                  domainMap_->getNode ()));
03324     checkInternalState ();
03325     return;
03326   }
03327 
03328 
03329   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
03330   void CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::mergeAllIndices()
03331   {
03332     TEUCHOS_TEST_FOR_EXCEPT( isGloballyIndexed() ); // call only after makeIndicesLocal()
03333     TEUCHOS_TEST_FOR_EXCEPT( ! isSorted() ); // call only after sortIndices()
03334     if (! isMerged ()) {
03335       for (size_t row=0; row < getNodeNumRows(); ++row) {
03336         RowInfo rowInfo = getRowInfo(row);
03337         mergeRowIndices(rowInfo);
03338       }
03339       // we just merged every row
03340       noRedundancies_ = true;
03341     }
03342   }
03343 
03344 
03345   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
03346   void
03347   CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::makeImportExport()
03348   {
03349     TEUCHOS_TEST_FOR_EXCEPTION(! hasColMap (), std::logic_error, "Tpetra::"
03350       "CrsGraph: It's not allowed to call makeImportExport() unless the graph "
03351       "has a column Map.");
03352     RCP<ParameterList> params = this->getNonconstParameterList (); // could be null
03353 
03354     // Don't do any checks to see if we need to create the Import, if
03355     // it exists already.
03356     //
03357     // FIXME (mfh 25 Mar 2013) This will become incorrect if we
03358     // change CrsGraph in the future to allow changing the column
03359     // Map after fillComplete.  For now, the column Map is fixed
03360     // after the first fillComplete call.
03361     if (importer_.is_null ()) {
03362       // Create the Import instance if necessary.
03363       if (domainMap_ != colMap_ && (! domainMap_->isSameAs (*colMap_))) {
03364         if (params.is_null () || ! params->isSublist ("Import")) {
03365           importer_ = rcp (new import_type (domainMap_, colMap_));
03366         } else {
03367           RCP<ParameterList> importSublist = sublist (params, "Import", true);
03368           importer_ = rcp (new import_type (domainMap_, colMap_, importSublist));
03369         }
03370       }
03371     }
03372 
03373     // Don't do any checks to see if we need to create the Export, if
03374     // it exists already.
03375     if (exporter_.is_null ()) {
03376       // Create the Export instance if necessary.
03377       if (rangeMap_ != rowMap_ && ! rangeMap_->isSameAs (*rowMap_)) {
03378         if (params.is_null () || ! params->isSublist ("Export")) {
03379           exporter_ = rcp (new export_type (rowMap_, rangeMap_));
03380         }
03381         else {
03382           RCP<ParameterList> exportSublist = sublist (params, "Export", true);
03383           exporter_ = rcp (new export_type (rowMap_, rangeMap_, exportSublist));
03384         }
03385       }
03386     }
03387   }
03388 
03389 
03390   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
03391   std::string
03392   CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::description() const
03393   {
03394     std::ostringstream oss;
03395     oss << DistObject<GlobalOrdinal,LocalOrdinal,GlobalOrdinal,Node>::description();
03396     if (isFillComplete()) {
03397       oss << "{status = fill complete"
03398           << ", global rows = " << getGlobalNumRows()
03399           << ", global cols = " << getGlobalNumCols()
03400           << ", global num entries = " << getGlobalNumEntries()
03401           << "}";
03402     }
03403     else {
03404       oss << "{status = fill not complete"
03405           << ", global rows = " << getGlobalNumRows()
03406           << "}";
03407     }
03408     return oss.str();
03409   }
03410 
03411 
03412   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
03413   void
03414   CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::
03415   describe (Teuchos::FancyOStream &out,
03416             const Teuchos::EVerbosityLevel verbLevel) const
03417   {
03418     const char tfecfFuncName[] = "describe()";
03419     using std::endl;
03420     using std::setw;
03421     using Teuchos::VERB_DEFAULT;
03422     using Teuchos::VERB_NONE;
03423     using Teuchos::VERB_LOW;
03424     using Teuchos::VERB_MEDIUM;
03425     using Teuchos::VERB_HIGH;
03426     using Teuchos::VERB_EXTREME;
03427     Teuchos::EVerbosityLevel vl = verbLevel;
03428     if (vl == VERB_DEFAULT) vl = VERB_LOW;
03429     RCP<const Comm<int> > comm = this->getComm();
03430     const int myImageID = comm->getRank(),
03431               numImages = comm->getSize();
03432     size_t width = 1;
03433     for (size_t dec=10; dec<getGlobalNumRows(); dec *= 10) {
03434       ++width;
03435     }
03436     width = std::max<size_t>(width,Teuchos::as<size_t>(11)) + 2;
03437     Teuchos::OSTab tab(out);
03438     //    none: print nothing
03439     //     low: print O(1) info from node 0
03440     //  medium: print O(P) info, num entries per node
03441     //    high: print O(N) info, num entries per row
03442     // extreme: print O(NNZ) info: print graph indices
03443     //
03444     // for medium and higher, print constituent objects at specified verbLevel
03445     if (vl != VERB_NONE) {
03446       if (myImageID == 0) out << this->description() << std::endl;
03447       // O(1) globals, minus what was already printed by description()
03448       if (isFillComplete() && myImageID == 0) {
03449         out << "Global number of diagonals = " << globalNumDiags_ << std::endl;
03450         out << "Global max number of row entries = " << globalMaxNumRowEntries_ << std::endl;
03451       }
03452       // constituent objects
03453       if (vl == VERB_MEDIUM || vl == VERB_HIGH || vl == VERB_EXTREME) {
03454         if (myImageID == 0) out << "\nRow map: " << std::endl;
03455         rowMap_->describe(out,vl);
03456         if (colMap_ != null) {
03457           if (myImageID == 0) out << "\nColumn map: " << std::endl;
03458           colMap_->describe(out,vl);
03459         }
03460         if (domainMap_ != null) {
03461           if (myImageID == 0) out << "\nDomain map: " << std::endl;
03462           domainMap_->describe(out,vl);
03463         }
03464         if (rangeMap_ != null) {
03465           if (myImageID == 0) out << "\nRange map: " << std::endl;
03466           rangeMap_->describe(out,vl);
03467         }
03468       }
03469       // O(P) data
03470       if (vl == VERB_MEDIUM || vl == VERB_HIGH || vl == VERB_EXTREME) {
03471         for (int imageCtr = 0; imageCtr < numImages; ++imageCtr) {
03472           if (myImageID == imageCtr) {
03473             out << "Node ID = " << imageCtr << std::endl
03474                 << "Node number of entries = " << nodeNumEntries_ << std::endl
03475                 << "Node number of diagonals = " << nodeNumDiags_ << std::endl
03476                 << "Node max number of entries = " << nodeMaxNumRowEntries_ << std::endl;
03477             if (indicesAreAllocated()) {
03478               out << "Node number of allocated entries = " << nodeNumAllocated_ << std::endl;
03479             }
03480             else {
03481               out << "Indices are not allocated." << std::endl;
03482             }
03483           }
03484           comm->barrier();
03485           comm->barrier();
03486           comm->barrier();
03487         }
03488       }
03489       // O(N) and O(NNZ) data
03490       if (vl == VERB_HIGH || vl == VERB_EXTREME) {
03491         TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(hasRowInfo() == false, std::runtime_error, ": reduce verbosity level; graph row information was deleted at fillComplete().");
03492         for (int imageCtr = 0; imageCtr < numImages; ++imageCtr) {
03493           if (myImageID == imageCtr) {
03494             out << std::setw(width) << "Node ID"
03495                 << std::setw(width) << "Global Row"
03496                 << std::setw(width) << "Num Entries";
03497             if (vl == VERB_EXTREME) {
03498               out << "Entries";
03499             }
03500             out << std::endl;
03501             for (size_t r=0; r < getNodeNumRows(); ++r) {
03502               RowInfo rowinfo = getRowInfo(r);
03503               GlobalOrdinal gid = rowMap_->getGlobalElement(r);
03504               out << std::setw(width) << myImageID
03505                   << std::setw(width) << gid
03506                   << std::setw(width) << rowinfo.numEntries;
03507               if (vl == VERB_EXTREME) {
03508                 if (isGloballyIndexed()) {
03509                   ArrayView<const GlobalOrdinal> rowview = getGlobalView(rowinfo);
03510                   for (size_t j=0; j < rowinfo.numEntries; ++j) out << rowview[j] << " ";
03511                 }
03512                 else if (isLocallyIndexed()) {
03513                   ArrayView<const LocalOrdinal> rowview = getLocalView(rowinfo);
03514                   for (size_t j=0; j < rowinfo.numEntries; ++j) out << colMap_->getGlobalElement(rowview[j]) << " ";
03515                 }
03516               }
03517               out << std::endl;
03518             }
03519           }
03520           comm->barrier();
03521           comm->barrier();
03522           comm->barrier();
03523         }
03524       }
03525     }
03526   }
03527 
03528 
03529   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
03530   bool
03531   CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::
03532   checkSizes (const SrcDistObject& source)
03533   {
03534     (void) source; // forestall "unused variable" compiler warnings
03535 
03536     // It's not clear what kind of compatibility checks on sizes can
03537     // be performed here.  Epetra_CrsGraph doesn't check any sizes for
03538     // compatibility.
03539     return true;
03540   }
03541 
03542 
03543   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
03544   void
03545   CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::
03546   copyAndPermute (const SrcDistObject& source,
03547                   size_t numSameIDs,
03548                   const Teuchos::ArrayView<const LocalOrdinal> &permuteToLIDs,
03549                   const Teuchos::ArrayView<const LocalOrdinal> &permuteFromLIDs)
03550   {
03551     using Teuchos::Array;
03552     using Teuchos::ArrayView;
03553     typedef LocalOrdinal LO;
03554     typedef GlobalOrdinal GO;
03555     const char tfecfFuncName[] = "copyAndPermute";
03556     typedef CrsGraph<LO, GO, Node, LocalMatOps> this_type;
03557     typedef RowGraph<LO, GO, Node> row_graph_type;
03558 
03559     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
03560       permuteToLIDs.size() != permuteFromLIDs.size(), std::runtime_error,
03561       ": permuteToLIDs and permuteFromLIDs must have the same size.");
03562     // Make sure that the source object has the right type.  We only
03563     // actually need it to be a RowGraph, with matching first three
03564     // template parameters.  If it's a CrsGraph, we can use view mode
03565     // instead of copy mode to get each row's data.
03566     //
03567     // FIXME (mfh 07 Jul 2013) It should not be necessary for any of
03568     // the template parameters but GO to match.  GO has to match
03569     // because the graph has to send indices as global ordinals, if
03570     // the source and target graphs do not have the same column Map.
03571     // If LO doesn't match, the graphs could communicate using global
03572     // indices.  It could be possible that Node affects the graph's
03573     // storage format, but packAndPrepare should assume a common
03574     // communication format in any case.
03575     const row_graph_type* srcRowGraph = dynamic_cast<const row_graph_type*> (&source);
03576     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
03577       srcRowGraph == NULL, std::invalid_argument,
03578       ": The source object must be a RowGraph with matching first three "
03579       "template parameters.");
03580 
03581     // If the source object is actually a CrsGraph, we can use view
03582     // mode instead of copy mode to access the entries in each row,
03583     // if the graph is not fill complete.
03584     const this_type* srcCrsGraph = dynamic_cast<const this_type*> (&source);
03585 
03586     const map_type& srcRowMap = * (srcRowGraph->getRowMap ());
03587     const map_type& tgtRowMap = * (this->getRowMap ());
03588     const bool src_filled = srcRowGraph->isFillComplete ();
03589     Array<GO> row_copy;
03590     LO myid = 0;
03591 
03592     //
03593     // "Copy" part of "copy and permute."
03594     //
03595     if (src_filled || srcCrsGraph == NULL) {
03596       // If the source graph is fill complete, we can't use view mode,
03597       // because the data might be stored in a different format not
03598       // compatible with the expectations of view mode.  Also, if the
03599       // source graph is not a CrsGraph, we can't use view mode,
03600       // because RowGraph only provides copy mode access to the data.
03601       for (size_t i = 0; i < numSameIDs; ++i, ++myid) {
03602         const GO gid = srcRowMap.getGlobalElement (myid);
03603         size_t row_length = srcRowGraph->getNumEntriesInGlobalRow (gid);
03604         row_copy.resize (row_length);
03605         size_t check_row_length = 0;
03606         srcRowGraph->getGlobalRowCopy (gid, row_copy (), check_row_length);
03607         this->insertGlobalIndices (gid, row_copy ());
03608       }
03609     } else {
03610       for (size_t i = 0; i < numSameIDs; ++i, ++myid) {
03611         const GO gid = srcRowMap.getGlobalElement (myid);
03612         ArrayView<const GO> row;
03613         srcCrsGraph->getGlobalRowView (gid, row);
03614         this->insertGlobalIndices (gid, row);
03615       }
03616     }
03617 
03618     //
03619     // "Permute" part of "copy and permute."
03620     //
03621     if (src_filled || srcCrsGraph == NULL) {
03622       for (LO i = 0; i < permuteToLIDs.size (); ++i) {
03623         const GO mygid = tgtRowMap.getGlobalElement (permuteToLIDs[i]);
03624         const GO srcgid = srcRowMap.getGlobalElement (permuteFromLIDs[i]);
03625         size_t row_length = srcRowGraph->getNumEntriesInGlobalRow (srcgid);
03626         row_copy.resize (row_length);
03627         size_t check_row_length = 0;
03628         srcRowGraph->getGlobalRowCopy (srcgid, row_copy (), check_row_length);
03629         this->insertGlobalIndices (mygid, row_copy ());
03630       }
03631     } else {
03632       for (LO i = 0; i < permuteToLIDs.size (); ++i) {
03633         const GO mygid = tgtRowMap.getGlobalElement (permuteToLIDs[i]);
03634         const GO srcgid = srcRowMap.getGlobalElement (permuteFromLIDs[i]);
03635         ArrayView<const GO> row;
03636         srcCrsGraph->getGlobalRowView (srcgid, row);
03637         this->insertGlobalIndices (mygid, row);
03638       }
03639     }
03640   }
03641 
03642 
03643   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
03644   void
03645   CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::
03646   packAndPrepare (const SrcDistObject& source,
03647                   const Teuchos::ArrayView<const LocalOrdinal> &exportLIDs,
03648                   Teuchos::Array<GlobalOrdinal> &exports,
03649                   const Teuchos::ArrayView<size_t> & numPacketsPerLID,
03650                   size_t& constantNumPackets,
03651                   Distributor &distor)
03652   {
03653     using Teuchos::Array;
03654     const char tfecfFuncName[] = "packAndPrepare";
03655     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
03656       exportLIDs.size() != numPacketsPerLID.size(), std::runtime_error,
03657       ": exportLIDs and numPacketsPerLID must have the same size.");
03658     typedef RowGraph<LocalOrdinal, GlobalOrdinal, Node> row_graph_type;
03659     const row_graph_type& srcGraph = dynamic_cast<const row_graph_type&> (source);
03660 
03661     // We don't check whether src_graph has had fillComplete called,
03662     // because it doesn't matter whether the *source* graph has been
03663     // fillComplete'd. The target graph can not be fillComplete'd yet.
03664     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
03665       this->isFillComplete (), std::runtime_error,
03666       ": The target graph of an Import or Export must not be fill complete.");
03667 
03668     srcGraph.pack (exportLIDs, exports, numPacketsPerLID, constantNumPackets, distor);
03669   }
03670 
03671 
03672   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
03673   void
03674   CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::
03675   pack (const Teuchos::ArrayView<const LocalOrdinal>& exportLIDs,
03676         Teuchos::Array<GlobalOrdinal>& exports,
03677         const Teuchos::ArrayView<size_t>& numPacketsPerLID,
03678         size_t& constantNumPackets,
03679         Distributor& distor) const
03680   {
03681     using Teuchos::Array;
03682     typedef LocalOrdinal LO;
03683     typedef GlobalOrdinal GO;
03684     const char tfecfFuncName[] = "pack";
03685     (void) distor; // forestall "unused argument" compiler warning
03686 
03687     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
03688       exportLIDs.size() != numPacketsPerLID.size(), std::runtime_error,
03689       ": exportLIDs and numPacketsPerLID must have the same size.");
03690 
03691     const map_type& srcMap = * (this->getMap ());
03692     constantNumPackets = 0;
03693 
03694     // Set numPacketsPerLID[i] to the number of entries owned by the
03695     // calling process in (local) row exportLIDs[i] of the graph, that
03696     // the caller wants us to send out.  Compute the total number of
03697     // packets (that is, entries) owned by this process in all the
03698     // rows that the caller wants us to send out.
03699     size_t totalNumPackets = 0;
03700     Array<GO> row;
03701     for (LO i = 0; i < exportLIDs.size (); ++i) {
03702       const GO GID = srcMap.getGlobalElement (exportLIDs[i]);
03703       size_t row_length = this->getNumEntriesInGlobalRow (GID);
03704       numPacketsPerLID[i] = row_length;
03705       totalNumPackets += row_length;
03706     }
03707 
03708     exports.resize (totalNumPackets);
03709 
03710     // Loop again over the rows to export, and pack rows of indices
03711     // into the output buffer.
03712     size_t exportsOffset = 0;
03713     for (LO i = 0; i < exportLIDs.size (); ++i) {
03714       const GO GID = srcMap.getGlobalElement (exportLIDs[i]);
03715       size_t row_length = this->getNumEntriesInGlobalRow (GID);
03716       row.resize (row_length);
03717       size_t check_row_length = 0;
03718       this->getGlobalRowCopy (GID, row (), check_row_length);
03719       typename Array<GO>::const_iterator row_iter = row.begin();
03720       typename Array<GO>::const_iterator row_end = row.end();
03721       size_t j = 0;
03722       for (; row_iter != row_end; ++row_iter, ++j) {
03723         exports[exportsOffset+j] = *row_iter;
03724       }
03725       exportsOffset += row.size();
03726     }
03727   }
03728 
03729 
03730   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
03731   void
03732   CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::
03733   unpackAndCombine (const Teuchos::ArrayView<const LocalOrdinal> &importLIDs,
03734                     const Teuchos::ArrayView<const GlobalOrdinal> &imports,
03735                     const Teuchos::ArrayView<size_t> &numPacketsPerLID,
03736                     size_t constantNumPackets,
03737                     Distributor& /* distor */,
03738                     CombineMode /* CM */)
03739   {
03740     // FIXME (mfh 02 Apr 2012) REPLACE combine mode has a perfectly
03741     // reasonable meaning, whether or not the matrix is fill complete.
03742     // It's just more work to implement.
03743 
03744     // We are not checking the value of the CombineMode input
03745     // argument.  For CrsGraph, we only support import/export
03746     // operations if fillComplete has not yet been called.  Any
03747     // incoming column-indices are inserted into the target graph. In
03748     // this context, CombineMode values of ADD vs INSERT are
03749     // equivalent. What is the meaning of REPLACE for CrsGraph? If a
03750     // duplicate column-index is inserted, it will be compressed out
03751     // when fillComplete is called.
03752     //
03753     // Note: I think REPLACE means that an existing row is replaced by
03754     // the imported row, i.e., the existing indices are cleared. CGB,
03755     // 6/17/2010
03756 
03757     const char tfecfFuncName[] = "unpackAndCombine";
03758     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
03759       importLIDs.size() != numPacketsPerLID.size(), std::runtime_error,
03760       ": importLIDs and numPacketsPerLID must have the same size.");
03761     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
03762       isFillComplete (), std::runtime_error,
03763       ": Import or Export operations are not allowed on the destination "
03764       "CrsGraph if it is fill complete.");
03765     size_t importsOffset = 0;
03766 
03767     typedef typename ArrayView<const LocalOrdinal>::const_iterator iter_type;
03768     iter_type impLIDiter = importLIDs.begin();
03769     iter_type impLIDend = importLIDs.end();
03770 
03771     for (size_t i = 0; impLIDiter != impLIDend; ++impLIDiter, ++i) {
03772       LocalOrdinal row_length = numPacketsPerLID[i];
03773       ArrayView<const GlobalOrdinal> row (&imports[importsOffset], row_length);
03774       insertGlobalIndicesFiltered (this->getMap ()->getGlobalElement (*impLIDiter), row);
03775       importsOffset += row_length;
03776     }
03777   }
03778 
03779 
03780   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
03781   void
03782   CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::
03783   removeEmptyProcessesInPlace (const Teuchos::RCP<const Map<LocalOrdinal, GlobalOrdinal, Node> >& newMap)
03784   {
03785     using Teuchos::Comm;
03786     using Teuchos::null;
03787     using Teuchos::ParameterList;
03788     using Teuchos::RCP;
03789 
03790     // We'll set all the state "transactionally," so that this method
03791     // satisfies the strong exception guarantee.  This object's state
03792     // won't be modified until the end of this method.
03793     RCP<const map_type> rowMap, domainMap, rangeMap, colMap;
03794     RCP<import_type> importer;
03795     RCP<export_type> exporter;
03796 
03797     rowMap = newMap;
03798     RCP<const Comm<int> > newComm =
03799       (newMap.is_null ()) ? null : newMap->getComm ();
03800 
03801     if (! domainMap_.is_null ()) {
03802       if (domainMap_.getRawPtr () == rowMap_.getRawPtr ()) {
03803         // Common case: original domain and row Maps are identical.
03804         // In that case, we need only replace the original domain Map
03805         // with the new Map.  This ensures that the new domain and row
03806         // Maps _stay_ identical.
03807         domainMap = newMap;
03808       } else {
03809         domainMap = domainMap_->replaceCommWithSubset (newComm);
03810       }
03811     }
03812     if (! rangeMap_.is_null ()) {
03813       if (rangeMap_.getRawPtr () == rowMap_.getRawPtr ()) {
03814         // Common case: original range and row Maps are identical.  In
03815         // that case, we need only replace the original range Map with
03816         // the new Map.  This ensures that the new range and row Maps
03817         // _stay_ identical.
03818         rangeMap = newMap;
03819       } else {
03820         rangeMap = rangeMap_->replaceCommWithSubset (newComm);
03821       }
03822     }
03823     if (! colMap.is_null ()) {
03824       colMap = colMap_->replaceCommWithSubset (newComm);
03825     }
03826 
03827     // (Re)create the Export and / or Import if necessary.
03828     if (! newComm.is_null ()) {
03829       RCP<ParameterList> params = this->getNonconstParameterList (); // could be null
03830       //
03831       // The operations below are collective on the new communicator.
03832       //
03833       // (Re)create the Export object if necessary.  If I haven't
03834       // called fillComplete yet, I don't have a rangeMap, so I must
03835       // first check if the _original_ rangeMap is not null.  Ditto
03836       // for the Import object and the domain Map.
03837       if (! rangeMap_.is_null () &&
03838           rangeMap != rowMap &&
03839           ! rangeMap->isSameAs (*rowMap)) {
03840         if (params.is_null () || ! params->isSublist ("Export")) {
03841           exporter = rcp (new export_type (rowMap, rangeMap));
03842         }
03843         else {
03844           RCP<ParameterList> exportSublist = sublist (params, "Export", true);
03845           exporter = rcp (new export_type (rowMap, rangeMap, exportSublist));
03846         }
03847       }
03848       // (Re)create the Import object if necessary.
03849       if (! domainMap_.is_null () &&
03850           domainMap != colMap &&
03851           ! domainMap->isSameAs (*colMap)) {
03852         if (params.is_null () || ! params->isSublist ("Import")) {
03853           importer = rcp (new import_type (domainMap, colMap));
03854         } else {
03855           RCP<ParameterList> importSublist = sublist (params, "Import", true);
03856           importer = rcp (new import_type (domainMap, colMap, importSublist));
03857         }
03858       }
03859     } // if newComm is not null
03860 
03861     // Defer side effects until the end.  If no destructors throw
03862     // exceptions (they shouldn't anyway), then this method satisfies
03863     // the strong exception guarantee.
03864     exporter_ = exporter;
03865     importer_ = importer;
03866     rowMap_ = rowMap;
03867     // mfh 31 Mar 2013: DistObject's map_ is the row Map of a CrsGraph
03868     // or CrsMatrix.  CrsGraph keeps a redundant pointer (rowMap_) to
03869     // the same object.  We might want to get rid of this redundant
03870     // pointer sometime, but for now, we'll leave it alone and just
03871     // set map_ to the same object.
03872     this->map_ = rowMap;
03873     domainMap_ = domainMap;
03874     rangeMap_ = rangeMap;
03875     colMap_ = colMap;
03876   }
03877 
03878 } // namespace Tpetra
03879 
03880 // Include KokkosRefactor partial specialisation if enabled
03881 #if defined(TPETRA_HAVE_KOKKOS_REFACTOR)
03882 #include "Tpetra_KokkosRefactor_CrsGraph_def.hpp"
03883 #endif
03884 
03885 
03886 //
03887 // Explicit instantiation macro
03888 //
03889 // Must be expanded from within the Tpetra namespace!
03890 //
03891 #define TPETRA_CRSGRAPH_GRAPH_INSTANT(LO,GO,NODE) template class CrsGraph< LO , GO , NODE >;
03892 #define TPETRA_CRSGRAPH_SORTROWINDICESANDVALUES_INSTANT(S,LO,GO,NODE) template void CrsGraph< LO , GO , NODE >::sortRowIndicesAndValues< S >(RowInfo, ArrayView< S > );
03893 #define TPETRA_CRSGRAPH_MERGEROWINDICESANDVALUES_INSTANT(S,LO,GO,NODE) template void CrsGraph< LO , GO , NODE >::mergeRowIndicesAndValues< S >(RowInfo, const ArrayView< S >& );
03894 #define TPETRA_CRSGRAPH_ALLOCATEVALUES1D_INSTANT(S,LO,GO,NODE) template ArrayRCP< S > CrsGraph< LO , GO , NODE >::allocateValues1D< S >() const;
03895 #define TPETRA_CRSGRAPH_ALLOCATEVALUES2D_INSTANT(S,LO,GO,NODE) template ArrayRCP< Array< S > > CrsGraph< LO , GO , NODE >::allocateValues2D< S >() const;
03896 
03897 #define TPETRA_CRSGRAPH_INSTANT(S,LO,GO,NODE)                    \
03898   TPETRA_CRSGRAPH_SORTROWINDICESANDVALUES_INSTANT(S,LO,GO,NODE)  \
03899   TPETRA_CRSGRAPH_MERGEROWINDICESANDVALUES_INSTANT(S,LO,GO,NODE) \
03900   TPETRA_CRSGRAPH_ALLOCATEVALUES1D_INSTANT(S,LO,GO,NODE)         \
03901   TPETRA_CRSGRAPH_ALLOCATEVALUES2D_INSTANT(S,LO,GO,NODE)
03902 
03903 #endif // TPETRA_CRSGRAPH_DEF_HPP
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines