Tpetra_CrsGraph_def.hpp

00001 //@HEADER
00002 // ************************************************************************
00003 // 
00004 //               Tpetra: Templated Linear Algebra Services Package 
00005 //                 Copyright (2008) Sandia Corporation
00006 // 
00007 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00008 // license for use of this work by or on behalf of the U.S. Government.
00009 // 
00010 // This library is free software; you can redistribute it and/or modify
00011 // it under the terms of the GNU Lesser General Public License as
00012 // published by the Free Software Foundation; either version 2.1 of the
00013 // License, or (at your option) any later version.
00014 //  
00015 // This library is distributed in the hope that it will be useful, but
00016 // WITHOUT ANY WARRANTY; without even the implied warranty of
00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00018 // Lesser General Public License for more details.
00019 //  
00020 // You should have received a copy of the GNU Lesser General Public
00021 // License along with this library; if not, write to the Free Software
00022 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00023 // USA
00024 // Questions? Contact Michael A. Heroux (maherou@sandia.gov) 
00025 // 
00026 // ************************************************************************
00027 //@HEADER
00028 
00029 #ifndef TPETRA_CRSGRAPH_DEF_HPP
00030 #define TPETRA_CRSGRAPH_DEF_HPP
00031 
00032 // TODO: filter column indices first in insertLocalIndices()
00033 // TODO: filter column indices first in insertGlobalIndices()
00034 
00035 #include <Kokkos_NodeTrace.hpp>
00036 
00037 #ifdef DOXYGEN_USE_ONLY
00038   #include "Tpetra_CrsGraph_decl.hpp"
00039 #endif
00040 
00041 namespace Tpetra {
00042 
00043   template <class LocalOrdinal, class GlobalOrdinal, class Node>
00044   CrsGraph<LocalOrdinal,GlobalOrdinal,Node>::CrsGraph(const Teuchos::RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > &rowMap, size_t maxNumEntriesPerRow, ProfileType pftype)
00045   : DistObject<GlobalOrdinal,LocalOrdinal,GlobalOrdinal,Node>(rowMap)
00046   , rowMap_(rowMap)
00047   , lclGraph_(rowMap->getNodeNumElements(), rowMap->getNode())
00048   , nodeNumEntries_(0)
00049   , nodeNumDiags_(0)
00050   , nodeMaxNumRowEntries_(0)
00051   , nodeNumAllocated_(Teuchos::OrdinalTraits<size_t>::invalid())
00052   , pftype_(pftype)
00053   , numAllocForAllRows_(maxNumEntriesPerRow)
00054   , indicesAreAllocated_(false)
00055   , indicesAreLocal_(false)
00056   , indicesAreGlobal_(false)
00057   , fillComplete_(false)
00058   , storageOptimized_(false)
00059   , lowerTriangular_(false)
00060   , upperTriangular_(false)
00061   , indicesAreSorted_(false)
00062   , noRedundancies_(false) {
00063     staticAssertions();
00064     TEST_FOR_EXCEPTION(maxNumEntriesPerRow > Teuchos::OrdinalTraits<size_t>::max() || (maxNumEntriesPerRow < 1 && maxNumEntriesPerRow != 0), std::runtime_error,
00065         Teuchos::typeName(*this) << "::CrsGraph(rowMap,maxNumEntriesPerRow): maxNumEntriesPerRow must be non-negative.");
00066     clearGlobalConstants();
00067     checkInternalState();
00068   }
00069 
00070 
00071   template <class LocalOrdinal, class GlobalOrdinal, class Node>
00072   CrsGraph<LocalOrdinal,GlobalOrdinal,Node>::CrsGraph(const Teuchos::RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > &rowMap, 
00073                                                       const Teuchos::RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > &colMap, 
00074                                                       size_t maxNumEntriesPerRow, ProfileType pftype)
00075   : DistObject<GlobalOrdinal,LocalOrdinal,GlobalOrdinal,Node>(rowMap)
00076   , rowMap_(rowMap)
00077   , colMap_(colMap)
00078   , lclGraph_(rowMap->getNodeNumElements(), rowMap->getNode())
00079   , nodeNumEntries_(0)
00080   , nodeNumDiags_(0)
00081   , nodeMaxNumRowEntries_(0)
00082   , nodeNumAllocated_(Teuchos::OrdinalTraits<size_t>::invalid())
00083   , pftype_(pftype)
00084   , numAllocForAllRows_(maxNumEntriesPerRow) 
00085   , indicesAreAllocated_(false)
00086   , indicesAreLocal_(false)
00087   , indicesAreGlobal_(false)
00088   , fillComplete_(false)
00089   , storageOptimized_(false)
00090   , lowerTriangular_(false)
00091   , upperTriangular_(false)
00092   , indicesAreSorted_(false)
00093   , noRedundancies_(false) {
00094     staticAssertions();
00095     TEST_FOR_EXCEPTION(maxNumEntriesPerRow > Teuchos::OrdinalTraits<size_t>::max() || (maxNumEntriesPerRow < 1 && maxNumEntriesPerRow != 0), std::runtime_error,
00096         Teuchos::typeName(*this) << "::CrsGraph(rowMap,colMap,maxNumEntriesPerRow): maxNumEntriesPerRow must be non-negative.");
00097     clearGlobalConstants();
00098     checkInternalState();
00099   }
00100 
00101 
00102   template <class LocalOrdinal, class GlobalOrdinal, class Node>
00103   CrsGraph<LocalOrdinal,GlobalOrdinal,Node>::CrsGraph(const Teuchos::RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > &rowMap, 
00104                                                       const Teuchos::ArrayRCP<const size_t> &NumEntriesPerRowToAlloc, ProfileType pftype)
00105   : DistObject<GlobalOrdinal,LocalOrdinal,GlobalOrdinal,Node>(rowMap)
00106   , rowMap_(rowMap)
00107   , lclGraph_(rowMap->getNodeNumElements(), rowMap->getNode())
00108   , nodeNumEntries_(0)
00109   , nodeNumDiags_(0)
00110   , nodeMaxNumRowEntries_(0)
00111   , nodeNumAllocated_(Teuchos::OrdinalTraits<size_t>::invalid())
00112   , pftype_(pftype)
00113   , numAllocPerRow_(NumEntriesPerRowToAlloc) 
00114   , numAllocForAllRows_(0)
00115   , indicesAreAllocated_(false)
00116   , indicesAreLocal_(false)
00117   , indicesAreGlobal_(false)
00118   , fillComplete_(false)
00119   , storageOptimized_(false)
00120   , lowerTriangular_(false)
00121   , upperTriangular_(false)
00122   , indicesAreSorted_(false)
00123   , noRedundancies_(false) {
00124     staticAssertions();
00125     TEST_FOR_EXCEPTION((size_t)NumEntriesPerRowToAlloc.size() != getNodeNumRows(), std::runtime_error,
00126         Teuchos::typeName(*this) << "::CrsGraph(rowMap,NumEntriesPerRowToAlloc): NumEntriesPerRowToAlloc must have as many entries as specified by rowMap for this node.");
00127     size_t numMin = Teuchos::OrdinalTraits<size_t>::max(),
00128            numMax = Teuchos::OrdinalTraits<size_t>::zero();
00129     for (size_t r=0; r < getNodeNumRows(); ++r) {
00130       numMin = std::min<size_t>( numMin, NumEntriesPerRowToAlloc[r] );
00131       numMax = std::max<size_t>( numMax, NumEntriesPerRowToAlloc[r] );
00132     }
00133     TEST_FOR_EXCEPTION((numMin < Teuchos::OrdinalTraits<size_t>::one() && numMin != Teuchos::OrdinalTraits<size_t>::zero()) || numMax > Teuchos::OrdinalTraits<size_t>::max(),
00134         std::runtime_error, Teuchos::typeName(*this) << "::CrsGraph(): Invalid user-specified number of non-zeros per row.");
00135     clearGlobalConstants();
00136     checkInternalState();
00137   }
00138 
00139 
00140   template <class LocalOrdinal, class GlobalOrdinal, class Node>
00141   CrsGraph<LocalOrdinal,GlobalOrdinal,Node>::CrsGraph(const Teuchos::RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > &rowMap, const Teuchos::RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > &colMap, const Teuchos::ArrayRCP<const size_t> &NumEntriesPerRowToAlloc, ProfileType pftype)
00142   : DistObject<GlobalOrdinal,LocalOrdinal,GlobalOrdinal,Node>(rowMap)
00143   , rowMap_(rowMap)
00144   , colMap_(colMap)
00145   , lclGraph_(rowMap->getNodeNumElements(), rowMap->getNode())
00146   , nodeNumEntries_(0)
00147   , nodeNumDiags_(0)
00148   , nodeMaxNumRowEntries_(0)
00149   , nodeNumAllocated_(Teuchos::OrdinalTraits<size_t>::invalid())
00150   , pftype_(pftype)
00151   , numAllocPerRow_(NumEntriesPerRowToAlloc) 
00152   , numAllocForAllRows_(0)
00153   , indicesAreAllocated_(false)
00154   , indicesAreLocal_(false)
00155   , indicesAreGlobal_(false)
00156   , fillComplete_(false)
00157   , storageOptimized_(false)
00158   , lowerTriangular_(false)
00159   , upperTriangular_(false)
00160   , indicesAreSorted_(false)
00161   , noRedundancies_(false) {
00162     staticAssertions();
00163     TEST_FOR_EXCEPTION((size_t)NumEntriesPerRowToAlloc.size() != getNodeNumRows(), std::runtime_error,
00164         Teuchos::typeName(*this) << "::CrsGraph(rowMap,colMap,NumEntriesPerRowToAlloc): NumEntriesPerRowToAlloc must have as many entries as specified by rowMap for this node.");
00165     clearGlobalConstants();
00166     checkInternalState();
00167   }
00168 
00169 
00170   template <class LocalOrdinal, class GlobalOrdinal, class Node>
00171   void CrsGraph<LocalOrdinal,GlobalOrdinal,Node>::allocateIndices(AllocateLocalGlobal lorg) {
00172     using Teuchos::null;
00173     // this is a protected function, only callable by us. if it was called incorrectly, it is our fault.
00174     TEST_FOR_EXCEPTION(isLocallyIndexed() && lorg==AllocateGlobal, std::logic_error,
00175         Teuchos::typeName(*this) << "::allocateIndices(): Internal logic error. Please contact Tpetra team.");
00176     TEST_FOR_EXCEPTION(isGloballyIndexed() && lorg==AllocateLocal, std::logic_error,
00177         Teuchos::typeName(*this) << "::allocateIndices(): Internal logic error. Please contact Tpetra team.");
00178     TEST_FOR_EXCEPTION( indicesAreAllocated() == true, std::logic_error, 
00179         Teuchos::typeName(*this) << "::allocateIndices(): Internal logic error. Please contact Tpetra team.");
00180     Teuchos::RCP<Node> node = lclGraph_.getNode();
00181     const size_t numRows = getNodeNumRows();
00182     indicesAreLocal_  = (lorg == AllocateLocal);
00183     indicesAreGlobal_ = (lorg == AllocateGlobal);
00184     // if we have no row, then we have nothing to do
00185     nodeNumAllocated_ = 0;
00186     if (numRows > 0) {
00187       if (getProfileType() == StaticProfile) {
00188         //
00189         //  STATIC ALLOCATION PROFILE
00190         //
00191         // determine how many entries to allocate and setup offsets into 1D arrays
00192         pbuf_rowOffsets_ = node->template allocBuffer<size_t>(numRows+1);
00193         KOKKOS_NODE_TRACE("CrsGraph::allocateIndices()")
00194         view_rowOffsets_ = node->template viewBufferNonConst<size_t>(Kokkos::WriteOnly,numRows+1,pbuf_rowOffsets_);
00195         if (numAllocPerRow_ != null) {
00196           // allocate offsets, get host view
00197           nodeNumAllocated_ = 0;
00198           for (size_t i=0; i < numRows; ++i) {
00199             view_rowOffsets_[i] = nodeNumAllocated_;
00200             nodeNumAllocated_ += numAllocPerRow_[i];
00201           }
00202           view_rowOffsets_[numRows] = nodeNumAllocated_;
00203         }
00204         else {
00205           nodeNumAllocated_ = numAllocForAllRows_ * numRows;
00206           view_rowOffsets_[0] = 0;
00207           for (size_t i=1; i <= numRows; ++i) {
00208             view_rowOffsets_[i] = view_rowOffsets_[i-1] + numAllocForAllRows_;
00209           }
00210         }
00211         // in hind-sight, we know nodeNumAllocated_ and whether pbuf_rowOffsets_ is needed at all
00212         if (nodeNumAllocated_ == 0) {
00213           view_rowOffsets_ = null;
00214           pbuf_rowOffsets_ = null;
00215         }
00216         // allocate the indices
00217         if (nodeNumAllocated_ > 0) {
00218           if (lorg == AllocateLocal) {
00219             pbuf_lclInds1D_ = node->template allocBuffer<LocalOrdinal>(nodeNumAllocated_);
00220             KOKKOS_NODE_TRACE("CrsGraph::allocateIndices()")
00221             view_lclInds1D_ = node->template viewBufferNonConst<LocalOrdinal>(Kokkos::WriteOnly, pbuf_lclInds1D_.size(), pbuf_lclInds1D_);
00222           }
00223           else {
00224             pbuf_gblInds1D_ = node->template allocBuffer<GlobalOrdinal>(nodeNumAllocated_);
00225             KOKKOS_NODE_TRACE("CrsGraph::allocateIndices()")
00226             view_gblInds1D_ = node->template viewBufferNonConst<GlobalOrdinal>(Kokkos::WriteOnly, pbuf_gblInds1D_.size(), pbuf_gblInds1D_);
00227           }
00228         }
00229       }
00230       else {
00231         //
00232         //  DYNAMIC ALLOCATION PROFILE
00233         //
00234         Teuchos::ArrayRCP<const size_t> numalloc = numAllocPerRow_;
00235         size_t howmany = numAllocForAllRows_;
00236         if (lorg == AllocateLocal) {
00237           pbuf_lclInds2D_ = Teuchos::arcp< Teuchos::ArrayRCP<LocalOrdinal> >(numRows);
00238           nodeNumAllocated_ = 0;
00239           for (size_t i=0; i < numRows; ++i) {
00240             if (numalloc != null) howmany = *numalloc++;
00241             nodeNumAllocated_ += howmany;
00242             if (howmany > 0) pbuf_lclInds2D_[i] = node->template allocBuffer<LocalOrdinal>(howmany);
00243           }
00244           if (nodeNumAllocated_ > 0) {
00245             view_lclInds2D_ = Kokkos::ArrayOfViewsHelper<Node>::template getArrayOfNonConstViews<LocalOrdinal>(node, Kokkos::WriteOnly, pbuf_lclInds2D_);
00246           }
00247           else {
00248             pbuf_lclInds2D_ = null;
00249           }
00250         }
00251         else { // allocate global indices
00252           pbuf_gblInds2D_ = Teuchos::arcp< Teuchos::ArrayRCP<GlobalOrdinal> >(numRows);
00253           nodeNumAllocated_ = 0;
00254           for (size_t i=0; i < numRows; ++i) {
00255             if (numalloc != null) howmany = *numalloc++;
00256             nodeNumAllocated_ += howmany;
00257             if (howmany > 0) pbuf_gblInds2D_[i] = node->template allocBuffer<GlobalOrdinal>(howmany);
00258           }
00259           if (nodeNumAllocated_ > 0) {
00260             view_gblInds2D_ = Kokkos::ArrayOfViewsHelper<Node>::template getArrayOfNonConstViews<GlobalOrdinal>(node, Kokkos::WriteOnly, pbuf_gblInds2D_);
00261           }
00262           else {
00263             pbuf_gblInds2D_ = null;
00264           }
00265         }
00266       }
00267       //
00268       //
00269       // 
00270       if (nodeNumAllocated_ > 0) {
00271         numEntriesPerRow_ = Teuchos::arcp<size_t>(numRows);
00272         std::fill(numEntriesPerRow_.begin(), numEntriesPerRow_.end(), 0);
00273       }
00274     } // if numRows > 0
00275     // done with these
00276     numAllocForAllRows_ = 0;
00277     numAllocPerRow_     = null;
00278     indicesAreAllocated_ = true;    
00279     checkInternalState();
00280   }
00281 
00282 
00283   template <class LocalOrdinal, class GlobalOrdinal, class Node>
00284   CrsGraph<LocalOrdinal,GlobalOrdinal,Node>::~CrsGraph()
00285   {}
00286 
00287 
00288   template <class LocalOrdinal, class GlobalOrdinal, class Node>
00289   global_size_t CrsGraph<LocalOrdinal,GlobalOrdinal,Node>::getGlobalNumRows() const { 
00290     return rowMap_->getGlobalNumElements(); 
00291   }
00292 
00293 
00294   template <class LocalOrdinal, class GlobalOrdinal, class Node>
00295   global_size_t CrsGraph<LocalOrdinal,GlobalOrdinal,Node>::getGlobalNumCols() const {
00296     TEST_FOR_EXCEPTION(hasColMap() == false, std::runtime_error,
00297       Teuchos::typeName(*this) << ": cannot call getGlobalNumCols() without column map.");
00298     return colMap_->getGlobalNumElements();
00299   }
00300 
00301 
00302   template <class LocalOrdinal, class GlobalOrdinal, class Node>
00303   size_t CrsGraph<LocalOrdinal,GlobalOrdinal,Node>::getNodeNumRows() const { 
00304     return rowMap_->getNodeNumElements(); 
00305   }
00306 
00307 
00308   template <class LocalOrdinal, class GlobalOrdinal, class Node>
00309   size_t CrsGraph<LocalOrdinal,GlobalOrdinal,Node>::getNodeNumCols() const {
00310     TEST_FOR_EXCEPTION(hasColMap() != true, std::runtime_error,
00311       Teuchos::typeName(*this) << ": cannot call getNodeNumCols() without a column map.");
00312     return colMap_->getNodeNumElements();
00313   }
00314 
00315 
00316   template <class LocalOrdinal, class GlobalOrdinal, class Node>
00317   size_t CrsGraph<LocalOrdinal,GlobalOrdinal,Node>::getNodeNumDiags() const { 
00318     TEST_FOR_EXCEPTION(!isFillComplete(), std::runtime_error,
00319       Teuchos::typeName(*this) << ": cannot call getNodeNumDiags() until fillComplete() has been called.");
00320     return nodeNumDiags_;
00321   }
00322 
00323 
00324   template <class LocalOrdinal, class GlobalOrdinal, class Node>
00325   global_size_t CrsGraph<LocalOrdinal,GlobalOrdinal,Node>::getGlobalNumDiags() const {
00326     return globalNumDiags_;
00327   }
00328 
00329 
00330   template <class LocalOrdinal, class GlobalOrdinal, class Node>
00331   Teuchos::RCP<Node>
00332   CrsGraph<LocalOrdinal,GlobalOrdinal,Node>::getNode() const {
00333     return lclGraph_.getNode();
00334   }
00335 
00336 
00337   template <class LocalOrdinal, class GlobalOrdinal, class Node>
00338   const Teuchos::RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > & 
00339   CrsGraph<LocalOrdinal,GlobalOrdinal,Node>::getRowMap() const {
00340     return rowMap_; 
00341   }
00342 
00343 
00344   template <class LocalOrdinal, class GlobalOrdinal, class Node>
00345   const Teuchos::RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > & 
00346   CrsGraph<LocalOrdinal,GlobalOrdinal,Node>::getColMap() const { 
00347     return colMap_; 
00348   }
00349 
00350 
00351   template <class LocalOrdinal, class GlobalOrdinal, class Node>
00352   const Teuchos::RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > & 
00353   CrsGraph<LocalOrdinal,GlobalOrdinal,Node>::getDomainMap() const {
00354     return domainMap_; 
00355   }
00356 
00357 
00358   template <class LocalOrdinal, class GlobalOrdinal, class Node>
00359   const Teuchos::RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > & 
00360   CrsGraph<LocalOrdinal,GlobalOrdinal,Node>::getRangeMap() const {
00361     return rangeMap_; 
00362   }
00363 
00364 
00365   template <class LocalOrdinal, class GlobalOrdinal, class Node>
00366   Teuchos::RCP<const Import<LocalOrdinal,GlobalOrdinal,Node> >
00367   CrsGraph<LocalOrdinal,GlobalOrdinal,Node>::getImporter() const { 
00368     return importer_;
00369   }
00370 
00371 
00372   template <class LocalOrdinal, class GlobalOrdinal, class Node>
00373   Teuchos::RCP<const Export<LocalOrdinal,GlobalOrdinal,Node> >
00374   CrsGraph<LocalOrdinal,GlobalOrdinal,Node>::getExporter() const {
00375     return exporter_;
00376   }
00377 
00378 
00379   template <class LocalOrdinal, class GlobalOrdinal, class Node>
00380   bool CrsGraph<LocalOrdinal,GlobalOrdinal,Node>::hasColMap() const { 
00381     return (colMap_ != Teuchos::null); 
00382   }
00383 
00384 
00385   template <class LocalOrdinal, class GlobalOrdinal, class Node>
00386   bool CrsGraph<LocalOrdinal,GlobalOrdinal,Node>::isStorageOptimized() const { 
00387     return storageOptimized_; 
00388   }
00389 
00390 
00391   template <class LocalOrdinal, class GlobalOrdinal, class Node>
00392   ProfileType CrsGraph<LocalOrdinal,GlobalOrdinal,Node>::getProfileType() const { 
00393     return pftype_; 
00394   }
00395 
00396 
00397   template <class LocalOrdinal, class GlobalOrdinal, class Node>
00398   global_size_t CrsGraph<LocalOrdinal,GlobalOrdinal,Node>::getGlobalNumEntries() const { 
00399     return globalNumEntries_; 
00400   }
00401 
00402 
00403   template <class LocalOrdinal, class GlobalOrdinal, class Node>
00404   size_t CrsGraph<LocalOrdinal,GlobalOrdinal,Node>::getNodeNumEntries() const { 
00405     return nodeNumEntries_; 
00406   }
00407 
00408 
00409   template <class LocalOrdinal, class GlobalOrdinal, class Node>
00410   global_size_t CrsGraph<LocalOrdinal,GlobalOrdinal,Node>::getGlobalMaxNumRowEntries() const {
00411     return globalMaxNumRowEntries_;
00412   }
00413 
00414 
00415   template <class LocalOrdinal, class GlobalOrdinal, class Node>
00416   size_t CrsGraph<LocalOrdinal,GlobalOrdinal,Node>::getNodeMaxNumRowEntries() const { 
00417     return nodeMaxNumRowEntries_; 
00418   }
00419 
00420 
00421   template <class LocalOrdinal, class GlobalOrdinal, class Node>
00422   bool CrsGraph<LocalOrdinal,GlobalOrdinal,Node>::isFillComplete() const { 
00423     return fillComplete_; 
00424   }
00425 
00426 
00427   template <class LocalOrdinal, class GlobalOrdinal, class Node>
00428   bool CrsGraph<LocalOrdinal,GlobalOrdinal,Node>::isUpperTriangular() const { 
00429     return upperTriangular_; 
00430   }
00431 
00432 
00433   template <class LocalOrdinal, class GlobalOrdinal, class Node>
00434   bool CrsGraph<LocalOrdinal,GlobalOrdinal,Node>::isLowerTriangular() const { 
00435     return lowerTriangular_; 
00436   }
00437 
00438 
00439   template <class LocalOrdinal, class GlobalOrdinal, class Node>
00440   bool CrsGraph<LocalOrdinal,GlobalOrdinal,Node>::isLocallyIndexed() const { 
00441     return indicesAreLocal_; 
00442   }
00443 
00444 
00445   template <class LocalOrdinal, class GlobalOrdinal, class Node>
00446   bool CrsGraph<LocalOrdinal,GlobalOrdinal,Node>::isGloballyIndexed() const { 
00447     return indicesAreGlobal_; 
00448   }
00449 
00450 
00451   template <class LocalOrdinal, class GlobalOrdinal, class Node>
00452   size_t CrsGraph<LocalOrdinal,GlobalOrdinal,Node>::getNodeAllocationSize() const {
00453     return nodeNumAllocated_;
00454   }
00455 
00456 
00457   template <class LocalOrdinal, class GlobalOrdinal, class Node>
00458   const Teuchos::RCP<const Teuchos::Comm<int> > &
00459   CrsGraph<LocalOrdinal,GlobalOrdinal,Node>::getComm() const { 
00460     return rowMap_->getComm();
00461   }
00462 
00463   template <class LocalOrdinal, class GlobalOrdinal, class Node>
00464   GlobalOrdinal CrsGraph<LocalOrdinal,GlobalOrdinal,Node>::getIndexBase() const { 
00465     return rowMap_->getIndexBase(); 
00466   }
00467 
00468   template <class LocalOrdinal, class GlobalOrdinal, class Node>
00469   bool CrsGraph<LocalOrdinal,GlobalOrdinal,Node>::isSorted() const { 
00470     return indicesAreSorted_; 
00471   }
00472 
00473   template <class LocalOrdinal, class GlobalOrdinal, class Node>
00474   bool CrsGraph<LocalOrdinal,GlobalOrdinal,Node>::notRedundant() const { 
00475     return noRedundancies_; 
00476   }
00477 
00478   template <class LocalOrdinal, class GlobalOrdinal, class Node>
00479   void CrsGraph<LocalOrdinal,GlobalOrdinal,Node>::setSorted(bool sorted) { 
00480     indicesAreSorted_ = sorted; 
00481   }
00482 
00483 
00484   template <class LocalOrdinal, class GlobalOrdinal, class Node>
00485   bool CrsGraph<LocalOrdinal,GlobalOrdinal,Node>::indicesAreAllocated() const { 
00486     return indicesAreAllocated_; 
00487   }
00488 
00489 
00492   //                                                                         //
00493   //                    Internal utility methods                             //
00494   //                                                                         //
00497 
00498 
00501   template <class LocalOrdinal, class GlobalOrdinal, class Node>
00502   RowInfo CrsGraph<LocalOrdinal,GlobalOrdinal,Node>::getRowInfo(size_t myRow) const {
00503 #ifdef HAVE_TPETRA_DEBUG
00504     TEST_FOR_EXCEPTION(rowMap_->isNodeLocalElement(myRow) == false, std::logic_error, 
00505         Teuchos::typeName(*this) << "::getRowInfo(): Internal logic error. Please contact Tpetra team.");
00506 #endif
00507     const size_t STINV = Teuchos::OrdinalTraits<size_t>::invalid();
00508     RowInfo ret;
00509     if (indicesAreAllocated() == false) {
00510       // haven't performed allocation yet; probably won't hit this code
00511       if (numAllocPerRow_ == Teuchos::null) {
00512         ret.allocSize = numAllocForAllRows_;
00513       }
00514       else {
00515         ret.allocSize = numAllocPerRow_[myRow];
00516       }
00517       ret.numEntries = 0;
00518       ret.offset1D = STINV;
00519     }
00520     else if (nodeNumAllocated_ == 0) {
00521       // have performed allocation, but the graph has no allocation or entries
00522       ret.allocSize = 0;
00523       ret.numEntries = 0;
00524       ret.offset1D = STINV;
00525     }
00526     else {
00527       // graph data structures have the info that we need
00528       //
00529       // if static graph, offsets tell us the allocation size
00530       if (getProfileType() == StaticProfile) {
00531         Teuchos::RCP<Node> node = lclGraph_.getNode();
00532         if (view_rowOffsets_ == Teuchos::null) {
00533           KOKKOS_NODE_TRACE("CrsGraph::getRowInfo()")
00534           Teuchos::ArrayRCP<const size_t> offs = node->template viewBuffer<size_t>(2,pbuf_rowOffsets_+myRow);
00535           ret.offset1D = offs[0];
00536           ret.allocSize = offs[1] - ret.offset1D;
00537           offs = Teuchos::null;
00538         }
00539         else {
00540           ret.offset1D = view_rowOffsets_[myRow];
00541           ret.allocSize = view_rowOffsets_[myRow+1] - ret.offset1D;
00542         }
00543       }
00544       else {
00545         if (isLocallyIndexed()) {
00546           ret.allocSize = pbuf_lclInds2D_[myRow].size();
00547         }
00548         else {
00549           ret.allocSize = pbuf_gblInds2D_[myRow].size();
00550         }
00551         ret.offset1D = STINV;
00552       }
00553       //
00554       // if storage is optimized, then allocSize == numnz
00555       if (isStorageOptimized()) {
00556         ret.numEntries = ret.allocSize;
00557       }
00558       else {
00559         ret.numEntries = numEntriesPerRow_[myRow];
00560       }
00561     }
00562     return ret;
00563   }
00564 
00565 
00568   template <class LocalOrdinal, class GlobalOrdinal, class Node>
00569   RowInfo CrsGraph<LocalOrdinal,GlobalOrdinal,Node>::getFullLocalView(
00570                   size_t myRow, 
00571                   Teuchos::ArrayRCP<const LocalOrdinal> &indices) const {
00572 #ifdef HAVE_TPETRA_DEBUG
00573     std::string err = Teuchos::typeName(*this) + "::getFullLocalView(): Internal logic error. Please contact Tpetra team.";
00574     TEST_FOR_EXCEPTION(isGloballyIndexed() == true, std::logic_error, err);
00575     TEST_FOR_EXCEPTION(rowMap_->isNodeLocalElement(myRow) == false, std::logic_error, err);
00576 #endif
00577     RowInfo sizeInfo = getRowInfo(myRow);
00578     indices = Teuchos::null;
00579     // getRowInfo does not specify whether node is allocated or not
00580     if (sizeInfo.allocSize > 0 && indicesAreAllocated_==true) {
00581       Teuchos::RCP<Node> node = lclGraph_.getNode();
00582       if (getProfileType() == StaticProfile) {
00583         if (view_lclInds1D_ == Teuchos::null) {
00584           KOKKOS_NODE_TRACE("CrsGraph::getFullLocalView()")
00585           indices = node->template viewBuffer<LocalOrdinal>(sizeInfo.allocSize, pbuf_lclInds1D_ + sizeInfo.offset1D);
00586         }
00587         else {
00588           indices = view_lclInds1D_.persistingView(sizeInfo.offset1D,sizeInfo.allocSize);
00589         }
00590       }
00591       else {  // dynamic profile
00592         if (view_lclInds2D_ == Teuchos::null) {
00593         KOKKOS_NODE_TRACE("CrsGraph::getFullLocalView()")
00594           indices = node->template viewBuffer<LocalOrdinal>(sizeInfo.allocSize, pbuf_lclInds2D_[myRow]);
00595         }
00596         else {
00597           indices = view_lclInds2D_[myRow];
00598         }
00599       }
00600     }
00601 #ifdef HAVE_TPETRA_DEBUG
00602     TEST_FOR_EXCEPTION(indicesAreAllocated_ == true && indices != Teuchos::null && static_cast<size_t>(indices.size()) != sizeInfo.allocSize, std::logic_error, err);
00603 #endif
00604     return sizeInfo;
00605   }
00606 
00607 
00610   template <class LocalOrdinal, class GlobalOrdinal, class Node>
00611   RowInfo CrsGraph<LocalOrdinal,GlobalOrdinal,Node>::getFullLocalViewNonConst(
00612                   size_t myRow, 
00613                   Teuchos::ArrayRCP<LocalOrdinal> &indices) {
00614 #ifdef HAVE_TPETRA_DEBUG
00615     std::string err = Teuchos::typeName(*this) + "::getFullLocalViewNonConst(): Internal logic error. Please contact Tpetra team.";
00616     TEST_FOR_EXCEPTION(isGloballyIndexed() == true, std::logic_error, err);
00617     TEST_FOR_EXCEPTION(rowMap_->isNodeLocalElement(myRow) == false, std::logic_error, err);
00618 #endif
00619     RowInfo sizeInfo = getRowInfo(myRow);
00620     indices = Teuchos::null;
00621     // getRowInfo does not specify whether node is allocated or not
00622     if (sizeInfo.allocSize > 0 && indicesAreAllocated_==true) {
00623       Teuchos::RCP<Node> node = lclGraph_.getNode();
00624       // if there are no valid entries, then this view can be constructed WriteOnly
00625       Kokkos::ReadWriteOption rw = (sizeInfo.numEntries == 0 ? Kokkos::WriteOnly : Kokkos::ReadWrite);
00626       if (getProfileType() == StaticProfile) {
00627         if (view_lclInds1D_ == Teuchos::null) {
00628           KOKKOS_NODE_TRACE("CrsGraph::getFullLocalViewNonConst()")
00629           indices = node->template viewBufferNonConst<LocalOrdinal>(rw, sizeInfo.allocSize, pbuf_lclInds1D_ + sizeInfo.offset1D);
00630         }
00631         else {
00632           indices = view_lclInds1D_.persistingView(sizeInfo.offset1D,sizeInfo.allocSize);
00633         }
00634       }
00635       else {  // dynamic profile
00636         if (view_lclInds2D_ == Teuchos::null) {
00637           KOKKOS_NODE_TRACE("CrsGraph::getFullLocalViewNonConst()")
00638           indices = node->template viewBufferNonConst<LocalOrdinal>(rw, sizeInfo.allocSize, pbuf_lclInds2D_[myRow]);
00639         }
00640         else {
00641           indices = view_lclInds2D_[myRow];
00642         }
00643       }
00644     }
00645 #ifdef HAVE_TPETRA_DEBUG
00646     TEST_FOR_EXCEPTION(indicesAreAllocated_ == true && indices != Teuchos::null && static_cast<size_t>(indices.size()) != sizeInfo.allocSize, std::logic_error, err);
00647 #endif
00648     return sizeInfo;
00649   }
00650 
00651 
00654   template <class LocalOrdinal, class GlobalOrdinal, class Node>
00655   RowInfo CrsGraph<LocalOrdinal,GlobalOrdinal,Node>::getFullGlobalView(
00656                   size_t myRow, 
00657                   Teuchos::ArrayRCP<const GlobalOrdinal> &indices) const {
00658 #ifdef HAVE_TPETRA_DEBUG
00659     std::string err = Teuchos::typeName(*this) + "::getFullGlobalView(): Internal logic error. Please contact Tpetra team.";
00660     TEST_FOR_EXCEPTION(isLocallyIndexed() == true, std::logic_error, err);
00661     TEST_FOR_EXCEPTION(rowMap_->isNodeLocalElement(myRow) == false, std::logic_error, err);
00662 #endif
00663     RowInfo sizeInfo = getRowInfo(myRow);
00664     indices = Teuchos::null;
00665     // getRowInfo does not specify whether node is allocated or not
00666     if (sizeInfo.allocSize > 0 && indicesAreAllocated_==true) {
00667       Teuchos::RCP<Node> node = lclGraph_.getNode();
00668       // if there are no valid entries, then this view can be constructed WriteOnly
00669       if (getProfileType() == StaticProfile) {
00670         if (view_gblInds1D_ == Teuchos::null) {
00671           KOKKOS_NODE_TRACE("CrsGraph::getFullGlobalView()")
00672           indices = node->template viewBuffer<GlobalOrdinal>(sizeInfo.allocSize, pbuf_gblInds1D_ + sizeInfo.offset1D);
00673         }
00674         else {
00675           indices = view_gblInds1D_.persistingView(sizeInfo.offset1D,sizeInfo.allocSize);
00676         }
00677       }
00678       else {  // dynamic profile
00679         if (view_gblInds2D_ == Teuchos::null) {
00680           KOKKOS_NODE_TRACE("CrsGraph::getFullGlobalView()")
00681           indices = node->template viewBuffer<GlobalOrdinal>(sizeInfo.allocSize, pbuf_gblInds2D_[myRow]);
00682         } 
00683         else {
00684           indices = view_gblInds2D_[myRow];
00685         }
00686       }
00687     }
00688 #ifdef HAVE_TPETRA_DEBUG
00689     TEST_FOR_EXCEPTION(indicesAreAllocated_ == true && indices != Teuchos::null && static_cast<size_t>(indices.size()) != sizeInfo.allocSize, std::logic_error, err);
00690 #endif
00691     return sizeInfo;
00692   }
00693 
00694 
00697   template <class LocalOrdinal, class GlobalOrdinal, class Node>
00698   RowInfo CrsGraph<LocalOrdinal,GlobalOrdinal,Node>::getFullGlobalViewNonConst(
00699                   size_t myRow, 
00700                   Teuchos::ArrayRCP<GlobalOrdinal> &indices) {
00701 #ifdef HAVE_TPETRA_DEBUG
00702     std::string err = Teuchos::typeName(*this) + "::getFullGlobalViewNonConst(): Internal logic error. Please contact Tpetra team.";
00703     TEST_FOR_EXCEPTION(isLocallyIndexed() == true, std::logic_error, err);
00704     TEST_FOR_EXCEPTION(rowMap_->isNodeLocalElement(myRow) == false, std::logic_error, err);
00705 #endif
00706     RowInfo sizeInfo = getRowInfo(myRow);
00707     indices = Teuchos::null;
00708     // getRowInfo does not specify whether node is allocated or not
00709     if (sizeInfo.allocSize > 0 && indicesAreAllocated_==true) {
00710       Teuchos::RCP<Node> node = lclGraph_.getNode();
00711       // if there are no valid entries, then this view can be constructed WriteOnly
00712       Kokkos::ReadWriteOption rw = (sizeInfo.numEntries == 0 ? Kokkos::WriteOnly : Kokkos::ReadWrite);
00713       if (getProfileType() == StaticProfile) {
00714         if (view_gblInds1D_ == Teuchos::null) {
00715           KOKKOS_NODE_TRACE("CrsGraph::getFullGlobalViewNonConst()")
00716           indices = node->template viewBufferNonConst<GlobalOrdinal>(rw, sizeInfo.allocSize, pbuf_gblInds1D_ + sizeInfo.offset1D);
00717         }
00718         else {
00719           indices = view_gblInds1D_.persistingView(sizeInfo.offset1D,sizeInfo.allocSize);
00720         }
00721       }
00722       else {  // dynamic profile
00723         if (view_gblInds2D_ == Teuchos::null) {
00724           KOKKOS_NODE_TRACE("CrsGraph::getFullGlobalViewNonConst()")
00725           indices = node->template viewBufferNonConst<GlobalOrdinal>(rw, sizeInfo.allocSize, pbuf_gblInds2D_[myRow]);
00726         }
00727         else {
00728           indices = view_gblInds2D_[myRow];
00729         }
00730       }
00731     }
00732 #ifdef HAVE_TPETRA_DEBUG
00733     TEST_FOR_EXCEPTION(indicesAreAllocated_ == true && indices != Teuchos::null && static_cast<size_t>(indices.size()) != sizeInfo.allocSize, std::logic_error, err);
00734 #endif
00735     return sizeInfo;
00736   }
00737 
00738 
00741   template <class LocalOrdinal, class GlobalOrdinal, class Node>
00742   void CrsGraph<LocalOrdinal,GlobalOrdinal,Node>::insertLocalIndicesViaView(
00743                                   size_t myRow, 
00744                                   const Teuchos::ArrayView<const LocalOrdinal> &indices, 
00745                                   const Teuchos::ArrayRCP<LocalOrdinal> &inds_view) {
00746 #ifdef HAVE_TPETRA_DEBUG
00747     std::string err = Teuchos::typeName(*this) + "::insertLocalIndicesViaView(): Internal logic error. Please contact Tpetra team.";
00748     TEST_FOR_EXCEPTION(isGloballyIndexed() == true, std::logic_error, err);
00749     TEST_FOR_EXCEPTION(rowMap_->isNodeLocalElement(myRow) == false, std::logic_error, err);
00750     {
00751       RowInfo sizeInfo = getRowInfo(myRow);
00752       TEST_FOR_EXCEPTION( sizeInfo.allocSize < sizeInfo.numEntries + indices.size(), std::logic_error, err );
00753     }
00754 #endif
00755     std::copy( indices.begin(), indices.end(), inds_view.begin() );
00756     numEntriesPerRow_[myRow] += indices.size();
00757     nodeNumEntries_ += indices.size();
00758     indicesAreSorted_ = false;
00759     noRedundancies_ = false;
00760     clearGlobalConstants();
00761   }
00762 
00763 
00766   template <class LocalOrdinal, class GlobalOrdinal, class Node>
00767   void CrsGraph<LocalOrdinal,GlobalOrdinal,Node>::insertGlobalIndicesViaView(
00768                                   size_t myRow, 
00769                                   const Teuchos::ArrayView<const GlobalOrdinal> &indices, 
00770                                   const Teuchos::ArrayRCP<GlobalOrdinal> &inds_view) {
00771 #ifdef HAVE_TPETRA_DEBUG
00772     std::string err = Teuchos::typeName(*this) + "::insertGlobalIndicesViaView(): Internal logic error. Please contact Tpetra team.";
00773     TEST_FOR_EXCEPTION(isLocallyIndexed() == true, std::logic_error, err);
00774     TEST_FOR_EXCEPTION(rowMap_->isNodeLocalElement(myRow) == false, std::logic_error, err);
00775     {
00776       RowInfo sizeInfo = getRowInfo(myRow);
00777       TEST_FOR_EXCEPTION( sizeInfo.allocSize < sizeInfo.numEntries + indices.size(), std::logic_error, err );
00778     }
00779 #endif
00780     std::copy( indices.begin(), indices.end(), inds_view.begin() );
00781     numEntriesPerRow_[myRow] += indices.size();
00782     nodeNumEntries_ += indices.size();
00783     indicesAreSorted_ = false;
00784     noRedundancies_ = false;
00785     clearGlobalConstants();
00786   }
00787 
00788 
00791   template <class LocalOrdinal, class GlobalOrdinal, class Node>
00792   size_t CrsGraph<LocalOrdinal,GlobalOrdinal,Node>::RNNZ(size_t row) const {
00793     size_t rnnz;
00794     if (isStorageOptimized()) {
00795       // if storage is optimized, then numalloc == numnz for every row
00796       rnnz = RNumAlloc(row);
00797     }
00798     else if (indicesAreAllocated() == false || nodeNumAllocated_ == 0) {
00799       // no allocated entries means no valid entries
00800       rnnz = 0;
00801     }
00802     else {
00803       rnnz = numEntriesPerRow_[row];
00804     }
00805     return rnnz;
00806   }
00807 
00808 
00811   template <class LocalOrdinal, class GlobalOrdinal, class Node>
00812   size_t CrsGraph<LocalOrdinal,GlobalOrdinal,Node>::RNumAlloc(size_t row) const {
00813     RowInfo sizeInfo = getRowInfo(row);
00814     return sizeInfo.allocSize;
00815   }
00816 
00817 
00820   template <class LocalOrdinal, class GlobalOrdinal, class Node>
00821   void CrsGraph<LocalOrdinal,GlobalOrdinal,Node>::staticAssertions() {
00822     using Teuchos::OrdinalTraits;
00823     // Assumption: sizeof(GlobalOrdinal) >= sizeof(LocalOrdinal)
00824     //    This is so that we can store LocalOrdinals in the memory formerly occupied by GlobalOrdinals
00825     // Assumption: max(GlobalOrdinal) >= max(LocalOrdinal)  and  max(size_t) >= max(LocalOrdinal)
00826     //    This is so that we can represent any LocalOrdinal as a size_t, and any LocalOrdinal as a GlobalOrdinal
00827     Teuchos::CompileTimeAssert<sizeof(GlobalOrdinal) < sizeof(LocalOrdinal)> cta_size;
00828     (void)cta_size;
00829     // can't call max() with CompileTimeAssert, because it isn't a constant expression; will need to make this a runtime check
00830     const char * err = ": Object cannot be allocated with stated template arguments: size assumptions are not valid.";
00831     TEST_FOR_EXCEPTION( (size_t)OrdinalTraits<LocalOrdinal>::max() > OrdinalTraits<size_t>::max(),          std::runtime_error, Teuchos::typeName(*this) << err);
00832     TEST_FOR_EXCEPTION( OrdinalTraits<LocalOrdinal>::max() > OrdinalTraits<GlobalOrdinal>::max(),   std::runtime_error, Teuchos::typeName(*this) << err);
00833     TEST_FOR_EXCEPTION( (size_t)OrdinalTraits<GlobalOrdinal>::max() > OrdinalTraits<global_size_t>::max(),  std::runtime_error, Teuchos::typeName(*this) << err);
00834     TEST_FOR_EXCEPTION( OrdinalTraits<size_t>::max() > OrdinalTraits<global_size_t>::max(),         std::runtime_error, Teuchos::typeName(*this) << err);
00835   }
00836 
00837 
00840   template <class LocalOrdinal, class GlobalOrdinal, class Node>
00841   void CrsGraph<LocalOrdinal,GlobalOrdinal,Node>::updateLocalAllocation(size_t lrow, size_t allocSize) {
00842     using Teuchos::ArrayRCP;
00843     Teuchos::RCP<Node> node = lclGraph_.getNode();
00844     const size_t curNA = RNumAlloc(lrow),
00845                   rnnz = RNNZ(lrow);
00846 #ifdef HAVE_TPETRA_DEBUG
00847     TEST_FOR_EXCEPT( rowMap_->isNodeLocalElement(lrow) == false );
00848     TEST_FOR_EXCEPT( allocSize < curNA );
00849     TEST_FOR_EXCEPT( isGloballyIndexed() );
00850     TEST_FOR_EXCEPT( allocSize == 0 );
00851     TEST_FOR_EXCEPT( indicesAreAllocated() == false );
00852 #endif
00853     // allocate a larger space for row "lrow"
00854     // copy any existing data from previous allocation to new allocation
00855     // update sizes
00856     // 
00857     // if we already have views of the data, we will create a new view and do the copy on the host
00858     // otherwise, don't create a view, and do the copy on the device
00859     // 
00860     if (pbuf_lclInds2D_ == Teuchos::null) {
00861       pbuf_lclInds2D_ = Teuchos::arcp< ArrayRCP<LocalOrdinal> >(getNodeNumRows());
00862       // if this is our initial allocation (pbuf_lclInds2D_ == null)
00863       // then go ahead and create views; this is our one chance to do it efficiently (i.e., WriteOnly)
00864       view_lclInds2D_ = Kokkos::ArrayOfViewsHelper<Node>::template getArrayOfNonConstViews<LocalOrdinal>(node, Kokkos::WriteOnly, pbuf_lclInds2D_);
00865     }
00866     // it is possible that pbuf_lclInds2D_ and view_lclInds2D_ point to the same space; in this case, we must get view_lclInds2D_ before modifying pbuf_lclInds2D_
00867     ArrayRCP<LocalOrdinal> old_alloc, old_view;
00868     old_alloc = pbuf_lclInds2D_[lrow];
00869     if (view_lclInds2D_ != Teuchos::null) {
00870       old_view = view_lclInds2D_[lrow];
00871     }
00872     // allocate new space
00873     pbuf_lclInds2D_[lrow] = node->template allocBuffer<LocalOrdinal>(allocSize);
00874     if (old_view == Teuchos::null) {
00875       // no view; copy on the device
00876       if (rnnz) {
00877         KOKKOS_NODE_TRACE("CrsGraph::updateLocalAllocation()")
00878         node->template copyBuffers<LocalOrdinal>(rnnz, old_alloc, pbuf_lclInds2D_[lrow]);
00879       }
00880       old_alloc = Teuchos::null;
00881     }
00882     else {
00883       // copy the old values to the new allocation via the views
00884       KOKKOS_NODE_TRACE("CrsGraph::updateLocalAllocation()")
00885       view_lclInds2D_[lrow] = node->template viewBufferNonConst<LocalOrdinal>(Kokkos::WriteOnly, allocSize, pbuf_lclInds2D_[lrow]);
00886       std::copy( old_view.begin(), old_view.begin() + rnnz, view_lclInds2D_[lrow].begin() );
00887       // delete the buffer early, this will prevent a copy back that we neither need nor want to pay for
00888       old_alloc = Teuchos::null;
00889       old_view  = Teuchos::null;
00890     }
00891     nodeNumAllocated_ += (allocSize - curNA);
00892     if (numEntriesPerRow_ == Teuchos::null) {
00893       numEntriesPerRow_ = Teuchos::arcp<size_t>( getNodeNumRows() );
00894       std::fill(numEntriesPerRow_.begin(), numEntriesPerRow_.end(), 0);
00895     }
00896   }
00897 
00898 
00901   template <class LocalOrdinal, class GlobalOrdinal, class Node>
00902   void CrsGraph<LocalOrdinal,GlobalOrdinal,Node>::updateGlobalAllocation(size_t lrow, size_t allocSize) {
00903     using Teuchos::ArrayRCP;
00904     Teuchos::RCP<Node> node = lclGraph_.getNode();
00905     const size_t curNA = RNumAlloc(lrow),
00906                   rnnz = RNNZ(lrow);
00907 #ifdef HAVE_TPETRA_DEBUG
00908     TEST_FOR_EXCEPT( rowMap_->isNodeLocalElement(lrow) == false );
00909     TEST_FOR_EXCEPT( allocSize < curNA );
00910     TEST_FOR_EXCEPT( isLocallyIndexed() );
00911     TEST_FOR_EXCEPT( allocSize == 0 );
00912     TEST_FOR_EXCEPT( indicesAreAllocated() == false );
00913 #endif
00914     // allocate a larger space for row "lrow"
00915     // copy any existing data from previous allocation to new allocation
00916     // update sizes
00917     // 
00918     // if we already have views of the data, we will create a new view and do the copy on the host
00919     // otherwise, don't create a view, and do the copy on the device
00920     // 
00921     if (pbuf_gblInds2D_ == Teuchos::null) {
00922       pbuf_gblInds2D_ = Teuchos::arcp< ArrayRCP<GlobalOrdinal> >(getNodeNumRows());
00923       // if this is our initial allocation (pbuf_gblInds2D_ == null). 
00924       // go ahead and create views; this is our one chance to do it efficiently (i.e., WriteOnly)
00925       view_gblInds2D_ = Kokkos::ArrayOfViewsHelper<Node>::template getArrayOfNonConstViews<GlobalOrdinal>(node, Kokkos::WriteOnly, pbuf_gblInds2D_);
00926     }
00927     // it is possible that pbuf_gblInds2D_ and view_gblInds2D_ point to the same space; in this case, we must get view_gblInds2D_ before modifying pbuf_gblInds2D_
00928     ArrayRCP<GlobalOrdinal> old_alloc, old_view;
00929     old_alloc = pbuf_gblInds2D_[lrow];
00930     if (view_gblInds2D_ != Teuchos::null) {
00931       old_view = view_gblInds2D_[lrow];
00932     }
00933     // allocate new space
00934     pbuf_gblInds2D_[lrow] = node->template allocBuffer<GlobalOrdinal>(allocSize);
00935     if (old_view == Teuchos::null) {
00936       // no view; copy on the device
00937       if (rnnz) {
00938         KOKKOS_NODE_TRACE("CrsGraph::updateGlobalAllocation()")
00939         node->template copyBuffers<GlobalOrdinal>(rnnz, old_alloc, pbuf_gblInds2D_[lrow]);
00940       }
00941       old_alloc = Teuchos::null;
00942     }
00943     else {
00944       // copy the old values to the new allocation via the views
00945       KOKKOS_NODE_TRACE("CrsGraph::updateGlobalAllocation()")
00946       view_gblInds2D_[lrow] = node->template viewBufferNonConst<GlobalOrdinal>(Kokkos::WriteOnly, allocSize, pbuf_gblInds2D_[lrow]);
00947       std::copy( old_view.begin(), old_view.begin() + rnnz, view_gblInds2D_[lrow].begin() );
00948       // delete the buffer early, this will prevent a copy back that we neither need nor want to pay for
00949       old_alloc = Teuchos::null;
00950       old_view = Teuchos::null;
00951     }
00952     nodeNumAllocated_ += (allocSize - curNA);
00953     if (numEntriesPerRow_ == Teuchos::null) {
00954       numEntriesPerRow_ = Teuchos::arcp<size_t>( getNodeNumRows() );
00955       std::fill(numEntriesPerRow_.begin(), numEntriesPerRow_.end(), 0);
00956     }
00957   }
00958 
00959 
00962   template <class LocalOrdinal, class GlobalOrdinal, class Node>
00963   size_t CrsGraph<LocalOrdinal,GlobalOrdinal,Node>::findLocalIndex(
00964                                 size_t row, 
00965                                 LocalOrdinal ind, 
00966                                 const Teuchos::ArrayRCP<const LocalOrdinal> &alreadyHaveAView) const {
00967     typedef typename Teuchos::ArrayRCP<const LocalOrdinal>::iterator IT;
00968     bool found = true;
00969     const size_t nE = RNNZ(row);
00970     // get a view of the row, if it wasn't passed by the caller
00971     Teuchos::ArrayRCP<const LocalOrdinal> rowview = alreadyHaveAView;
00972     if (rowview == Teuchos::null) {
00973       rowview = getLocalRowView(row);
00974     }
00975     IT rptr, locptr;
00976     rptr = rowview.begin();
00977     if (isSorted()) {
00978       // binary search
00979       std::pair<IT,IT> p = std::equal_range(rptr,rptr+nE,ind);
00980       if (p.first == p.second) found = false;
00981       else locptr = p.first;
00982     }
00983     else {
00984       // direct search
00985       locptr = std::find(rptr,rptr+nE,ind);
00986       if (locptr == rptr+nE) found = false;
00987     }
00988     size_t ret;
00989     if (!found) {
00990       ret = Teuchos::OrdinalTraits<size_t>::invalid();
00991     }
00992     else {
00993       ret = (locptr - rptr);
00994     }
00995     locptr = Teuchos::NullIteratorTraits<IT>::getNull();
00996     rptr = Teuchos::NullIteratorTraits<IT>::getNull();
00997     rowview = Teuchos::null;
00998     return ret;
00999   }
01000 
01001 
01004   template <class LocalOrdinal, class GlobalOrdinal, class Node>
01005   size_t CrsGraph<LocalOrdinal,GlobalOrdinal,Node>::findGlobalIndex(
01006                                 size_t row, 
01007                                 GlobalOrdinal ind, 
01008                                 const Teuchos::ArrayRCP<const GlobalOrdinal> &alreadyHaveAView) const {
01009     typedef typename Teuchos::ArrayRCP<const GlobalOrdinal>::iterator IT;
01010     bool found = true;
01011     const size_t nE = RNNZ(row);
01012     // get a view of the row, if it wasn't passed by the caller
01013     Teuchos::ArrayRCP<const GlobalOrdinal> rowview = alreadyHaveAView;
01014     if (rowview == Teuchos::null) {
01015       rowview = getGlobalRowView(row);
01016     }
01017     IT rptr, locptr;
01018     rptr = rowview.begin();
01019     if (isSorted()) {
01020       // binary search
01021       std::pair<IT,IT> p = std::equal_range(rptr,rptr+nE,ind);
01022       if (p.first == p.second) found = false;
01023       else locptr = p.first;
01024     }
01025     else {
01026       // direct search
01027       locptr = std::find(rptr,rptr+nE,ind);
01028       if (locptr == rptr+nE) found = false;
01029     }
01030     size_t ret;
01031     if (!found) {
01032       ret = Teuchos::OrdinalTraits<size_t>::invalid();
01033     }
01034     else {
01035       ret = (locptr - rptr);
01036     }
01037     locptr = Teuchos::NullIteratorTraits<IT>::getNull();
01038     rptr = Teuchos::NullIteratorTraits<IT>::getNull();
01039     rowview = Teuchos::null;
01040     return ret;
01041   }
01042 
01043 
01046   template <class LocalOrdinal, class GlobalOrdinal, class Node>
01047   void CrsGraph<LocalOrdinal,GlobalOrdinal,Node>::clearGlobalConstants() {
01048     globalNumEntries_ = Teuchos::OrdinalTraits<global_size_t>::invalid();
01049     globalNumDiags_ = Teuchos::OrdinalTraits<global_size_t>::invalid();
01050     globalMaxNumRowEntries_ = Teuchos::OrdinalTraits<global_size_t>::invalid();
01051     haveGlobalConstants_ = false;
01052   }
01053 
01054 
01057   template <class LocalOrdinal, class GlobalOrdinal, class Node>
01058   void CrsGraph<LocalOrdinal,GlobalOrdinal,Node>::checkInternalState() const {
01059 #ifdef HAVE_TPETRA_DEBUG
01060     Teuchos::RCP<Node> node = lclGraph_.getNode();
01061     const global_size_t gsti = Teuchos::OrdinalTraits<global_size_t>::invalid();
01062     const size_t         sti = Teuchos::OrdinalTraits<size_t>::invalid();
01063     using Teuchos::null;
01064     std::string err = Teuchos::typeName(*this) + "::checkInternalState(): Likely internal logic error. Please contact Tpetra team.";
01065     // check the internal state of this data structure
01066     // this is called by numerous state-changing methods, in a debug build, to ensure that the object 
01067     // always remains in a valid state
01068     // the graph should have been allocated with a row map
01069     TEST_FOR_EXCEPTION( rowMap_ == null, std::logic_error, err );
01070     // if the graph thinks it has a column map, then it had better have one
01071     TEST_FOR_EXCEPTION( hasColMap() == (colMap_ == null), std::logic_error, err );
01072     // if the graph has been fill completed, then all maps should be present
01073     TEST_FOR_EXCEPTION( isFillComplete() == true && (colMap_ == null || rangeMap_ == null || domainMap_ == null), std::logic_error, err );
01074     // if storage has been optimized, then indices should have been allocated (even if trivially so)
01075     TEST_FOR_EXCEPTION( isStorageOptimized() == true && indicesAreAllocated() == false, std::logic_error, err );
01076     // if storage has been optimized, then number of allocated is now the number of entries
01077     TEST_FOR_EXCEPTION( isStorageOptimized() == true && nodeNumAllocated_ != nodeNumEntries_, std::logic_error, err );
01078     // if graph doesn't have the global constants, then they should all be marked as invalid
01079     TEST_FOR_EXCEPTION( haveGlobalConstants_ == false && ( globalNumEntries_ != gsti || globalNumDiags_ != gsti || globalMaxNumRowEntries_ != gsti ), std::logic_error, err ); 
01080     // if the graph has global cosntants, then they should be valid.
01081     TEST_FOR_EXCEPTION( haveGlobalConstants_ == true && ( globalNumEntries_ == gsti || globalNumDiags_ == gsti || globalMaxNumRowEntries_ == gsti ), std::logic_error, err ); 
01082     TEST_FOR_EXCEPTION( haveGlobalConstants_ == true && ( globalNumEntries_ < nodeNumEntries_ || globalNumDiags_ < nodeNumDiags_ || globalMaxNumRowEntries_ < nodeMaxNumRowEntries_ ),
01083                         std::logic_error, err );
01084     // no view should be present in the absence of a buffer
01085     TEST_FOR_EXCEPTION( (pbuf_lclInds1D_  == Teuchos::null && view_lclInds1D_  != Teuchos::null) ||
01086                         (pbuf_gblInds1D_  == Teuchos::null && view_gblInds1D_  != Teuchos::null) ||
01087                         (pbuf_lclInds2D_  == Teuchos::null && view_lclInds2D_  != Teuchos::null) ||
01088                         (pbuf_gblInds2D_  == Teuchos::null && view_gblInds2D_  != Teuchos::null) ||
01089                         (pbuf_rowOffsets_ == Teuchos::null && view_rowOffsets_ != Teuchos::null), std::logic_error, err );
01090     // if indices are allocated, then the information dictating the allocation quantities should be freed
01091     TEST_FOR_EXCEPTION( indicesAreAllocated() == true  && (numAllocForAllRows_ != 0 || numAllocPerRow_ != null),  std::logic_error, err );
01092     // if indices are not allocated, then information dictating allocation quantities should be present
01093     TEST_FOR_EXCEPTION( indicesAreAllocated() == false && (nodeNumAllocated_ != sti || nodeNumEntries_ != 0),     std::logic_error, err );
01094     // if indices are not allocated, then views should be null
01095     TEST_FOR_EXCEPTION( indicesAreAllocated() == false && (view_lclInds1D_ != null || view_gblInds1D_ != null || view_lclInds2D_ != null || view_gblInds2D_ != null || view_rowOffsets_ != null), std::logic_error, err );
01096     // if storage is optimized, then profile should be static
01097     TEST_FOR_EXCEPTION( isStorageOptimized() && pftype_ != StaticProfile, std::logic_error, err );
01098     // if profile is dynamic and we have a non-trivial allocation, then 2D allocations should be present
01099     TEST_FOR_EXCEPTION( pftype_ == DynamicProfile && indicesAreAllocated() && nodeNumAllocated_ > 0 && pbuf_lclInds2D_ == null && pbuf_gblInds2D_ == null, std::logic_error, err );
01100     // if profile is dynamic, then 1D allocations should not be present
01101     TEST_FOR_EXCEPTION( pftype_ == DynamicProfile && (pbuf_lclInds1D_ != null || pbuf_gblInds1D_ != null), std::logic_error, err );
01102     // if profile is static and we have a non-trivial allocation, then 1D allocations should be present
01103     TEST_FOR_EXCEPTION( pftype_ == StaticProfile && indicesAreAllocated() && nodeNumAllocated_ > 0 && pbuf_lclInds1D_ == null && pbuf_gblInds1D_ == null, std::logic_error, err );
01104     // if profile is static, then 2D allocations should not be present
01105     TEST_FOR_EXCEPTION( pftype_ == StaticProfile && (pbuf_lclInds2D_ != null || pbuf_gblInds2D_ != null), std::logic_error, err );
01106     // if profile is dynamic, then row offsets should not be allocated (they are used only for 1D indexing)
01107     TEST_FOR_EXCEPTION( pftype_ == DynamicProfile && pbuf_rowOffsets_ != null, std::logic_error, err );
01108     // if profile is static and we have a non-trivial application, then roww offsets should be allocated
01109     TEST_FOR_EXCEPTION( pftype_ == StaticProfile && indicesAreAllocated() && nodeNumAllocated_ > 0 && pbuf_rowOffsets_ == null, std::logic_error, err );
01110     // if indices are not allocated, then non of the buffers should be.
01111     TEST_FOR_EXCEPTION( indicesAreAllocated() == false && (pbuf_rowOffsets_ != null || numEntriesPerRow_ != null||
01112                                                           pbuf_lclInds1D_ != null || pbuf_lclInds2D_ != null ||
01113                                                           pbuf_gblInds1D_ != null || pbuf_gblInds2D_ != null), std::logic_error, err );
01114     // for a trivial (i.e., zero) allocation, row offsets and num entries should be freed, because they are not needed
01115     TEST_FOR_EXCEPTION( nodeNumAllocated_ == 0 && (pbuf_rowOffsets_ != null || numEntriesPerRow_ != null), std::logic_error, err );
01116     // for a non-trivial allocation with optimal storage, num entries is redundant and therefore should be freed
01117     TEST_FOR_EXCEPTION( indicesAreAllocated() && nodeNumAllocated_ > 0 && storageOptimized_ == true  && numEntriesPerRow_ != null, std::logic_error, err );
01118     // for a non-trivial allocation without optimal storage, num entries is necessary and should be present
01119     TEST_FOR_EXCEPTION( indicesAreAllocated() && nodeNumAllocated_ > 0 && storageOptimized_ == false && numEntriesPerRow_ == null, std::logic_error, err );
01120     // indices may be local or global only if they are allocated (numAllocated is redundant; could simply be indicesAreLocal_ || indicesAreGlobal_)
01121     TEST_FOR_EXCEPTION( (indicesAreLocal_ == true || indicesAreGlobal_ == true) && indicesAreAllocated_ == false, std::logic_error, err );
01122     // indices may be local or global, but not both
01123     TEST_FOR_EXCEPTION( indicesAreLocal_ == true && indicesAreGlobal_ == true, std::logic_error, err );
01124     // if indices are local, then global allocations should not be present
01125     TEST_FOR_EXCEPTION( indicesAreLocal_ == true && (pbuf_gblInds1D_ != null || pbuf_gblInds2D_ != null), std::logic_error, err );
01126     // if indices are global, then local allocations should not be present
01127     TEST_FOR_EXCEPTION( indicesAreGlobal_ == true && (pbuf_lclInds1D_ != null || pbuf_lclInds2D_ != null), std::logic_error, err );
01128     // if indices are local and non-trivial, then local allocations should be present
01129     TEST_FOR_EXCEPTION( indicesAreLocal_ == true && nodeNumAllocated_ > 0 && pbuf_lclInds1D_ == null && pbuf_lclInds2D_ == null, std::logic_error, err );
01130     // if indices are global and non-trivial, then global allocations should be present
01131     TEST_FOR_EXCEPTION( indicesAreGlobal_ == true && nodeNumAllocated_ > 0 && pbuf_gblInds1D_ == null && pbuf_gblInds2D_ == null, std::logic_error, err );
01132     // if indices are allocated, then we should have recorded how many were allocated
01133     TEST_FOR_EXCEPTION( indicesAreAllocated() == true  && nodeNumAllocated_ == sti, std::logic_error, err );
01134     // if indices are not allocated, then the allocation size should be marked invalid
01135     TEST_FOR_EXCEPTION( indicesAreAllocated() == false && nodeNumAllocated_ != sti, std::logic_error, err );
01136     // check the actual allocations
01137     size_t actualNumAllocated = 0;
01138     if (pftype_ == DynamicProfile) {
01139       if (isGloballyIndexed() && pbuf_gblInds2D_ != Teuchos::null) {
01140         for (size_t r = 0; r < getNodeNumRows(); ++r) {
01141           actualNumAllocated += pbuf_gblInds2D_[r].size();
01142         }
01143       }
01144       else if (isLocallyIndexed() && pbuf_lclInds2D_ != Teuchos::null) {
01145         for (size_t r = 0; r < getNodeNumRows(); ++r) {
01146           actualNumAllocated += pbuf_lclInds2D_[r].size();
01147         }
01148       }
01149     }
01150     else { // pftype_ == StaticProfile)
01151       if (pbuf_rowOffsets_ != Teuchos::null) {
01152         if (view_rowOffsets_ == Teuchos::null) {
01153           KOKKOS_NODE_TRACE("CrsGraph::checkInternalState()")
01154           Teuchos::ArrayRCP<const size_t> last_offset = node->template viewBuffer<size_t>(1,pbuf_rowOffsets_+getNodeNumRows());
01155           actualNumAllocated = last_offset[0];
01156         }
01157         else {
01158           actualNumAllocated = view_rowOffsets_[getNodeNumRows()];
01159         }
01160       }
01161       else {
01162         actualNumAllocated = 0;
01163       }
01164       TEST_FOR_EXCEPTION(  isLocallyIndexed() == true && (size_t)pbuf_lclInds1D_.size() != actualNumAllocated, std::logic_error, err );
01165       TEST_FOR_EXCEPTION( isGloballyIndexed() == true && (size_t)pbuf_gblInds1D_.size() != actualNumAllocated, std::logic_error, err );
01166     }
01167     TEST_FOR_EXCEPTION(indicesAreAllocated() == true && actualNumAllocated != nodeNumAllocated_, std::logic_error, err );
01168 #endif
01169   }
01170 
01171 
01174   template <class LocalOrdinal, class GlobalOrdinal, class Node>
01175   void CrsGraph<LocalOrdinal,GlobalOrdinal,Node>::fillLocalGraph() {
01176 #ifdef HAVE_TPETRA_DEBUG
01177     TEST_FOR_EXCEPTION( view_lclInds1D_ != Teuchos::null || view_lclInds2D_ != Teuchos::null || view_rowOffsets_ != Teuchos::null ||
01178                         view_gblInds1D_ != Teuchos::null || view_gblInds2D_ != Teuchos::null, std::logic_error, 
01179         Teuchos::typeName(*this) << "::fillLocalGraph(): Internal logic error. Please contact Tpetra team.");
01180 #endif
01181     lclGraph_.clear();
01182     if (storageOptimized_) {
01183       // fill packed matrix; it is okay for pbuf_lclInds1D_ to be null; the matrix will flag itself as empty
01184       lclGraph_.setPackedStructure(pbuf_rowOffsets_, pbuf_lclInds1D_);
01185     }
01186     else if (getProfileType() == StaticProfile) {
01187       if (pbuf_lclInds1D_ != Teuchos::null) {
01188         const size_t nlrs = getNodeNumRows();
01189         for (size_t r=0; r < nlrs; ++r) {
01190           RowInfo sizeInfo = getRowInfo(r);
01191           Teuchos::ArrayRCP<const LocalOrdinal> rowinds;
01192           if (sizeInfo.numEntries > 0) {
01193             rowinds = pbuf_lclInds1D_.persistingView(sizeInfo.offset1D, sizeInfo.numEntries);
01194             lclGraph_.set2DIndices(r,rowinds);
01195           }
01196         }
01197       }
01198     }
01199     else if (getProfileType() == DynamicProfile) {
01200       if (pbuf_lclInds2D_ != Teuchos::null) {
01201         const size_t nlrs = getNodeNumRows();
01202         for (size_t r=0; r < nlrs; ++r) {
01203           RowInfo sizeInfo = getRowInfo(r);
01204           Teuchos::ArrayRCP<const LocalOrdinal> rowinds = pbuf_lclInds2D_[r];
01205           if (sizeInfo.numEntries > 0) {
01206             rowinds = rowinds.persistingView(0,sizeInfo.numEntries);
01207             lclGraph_.set2DIndices(r,rowinds);
01208           }
01209         }
01210       }
01211     }
01212   }
01213 
01214 
01217   //                                                                         //
01218   //                  User-visible class methods                             //
01219   //                                                                         //
01222 
01223 
01226   template <class LocalOrdinal, class GlobalOrdinal, class Node>
01227   size_t CrsGraph<LocalOrdinal,GlobalOrdinal,Node>::getNumEntriesInGlobalRow(GlobalOrdinal globalRow) const {
01228     using Teuchos::OrdinalTraits;
01229     const LocalOrdinal rlid = rowMap_->getLocalElement(globalRow);
01230     size_t ret;
01231     if (rlid == OrdinalTraits<LocalOrdinal>::invalid()) {
01232       ret = OrdinalTraits<size_t>::invalid();
01233     }
01234     else { 
01235       ret = RNNZ(rlid);
01236     }
01237     return ret;
01238   }
01239 
01240 
01243   template <class LocalOrdinal, class GlobalOrdinal, class Node>
01244   size_t CrsGraph<LocalOrdinal,GlobalOrdinal,Node>::getNumEntriesInLocalRow(LocalOrdinal localRow) const {
01245     using Teuchos::OrdinalTraits;
01246     size_t ret;
01247     if (!rowMap_->isNodeLocalElement(localRow)) {
01248       ret = OrdinalTraits<size_t>::invalid();
01249     }
01250     else {
01251       ret = RNNZ(localRow);
01252     }
01253     return ret;
01254   }
01255 
01256 
01259   template <class LocalOrdinal, class GlobalOrdinal, class Node>
01260   size_t CrsGraph<LocalOrdinal,GlobalOrdinal,Node>::getNumAllocatedEntriesInGlobalRow(GlobalOrdinal globalRow) const {
01261     using Teuchos::OrdinalTraits;
01262     const LocalOrdinal rlid = rowMap_->getLocalElement(globalRow);
01263     size_t ret;
01264     if (rlid == OrdinalTraits<LocalOrdinal>::invalid()) {
01265       ret = OrdinalTraits<size_t>::invalid();
01266     }
01267     else {
01268       ret = RNumAlloc(rlid);
01269     }
01270     return ret;
01271   }
01272 
01273 
01276   template <class LocalOrdinal, class GlobalOrdinal, class Node>
01277   size_t CrsGraph<LocalOrdinal,GlobalOrdinal,Node>::getNumAllocatedEntriesInLocalRow(LocalOrdinal localRow) const {
01278     using Teuchos::OrdinalTraits;
01279     size_t ret;
01280     if (!rowMap_->isNodeLocalElement(localRow)) {
01281       ret = OrdinalTraits<size_t>::invalid();
01282     }
01283     else {
01284       ret = RNumAlloc(localRow);
01285     }
01286     return ret; 
01287   }
01288 
01289 
01292   template <class LocalOrdinal, class GlobalOrdinal, class Node>
01293   Teuchos::ArrayRCP<const LocalOrdinal> CrsGraph<LocalOrdinal,GlobalOrdinal,Node>::getLocalRowView(LocalOrdinal LocalRow) const {
01294     TEST_FOR_EXCEPTION(isGloballyIndexed() == true, std::runtime_error,
01295         Teuchos::typeName(*this) << "::getLocalRowView(): local indices do not exist.");
01296     TEST_FOR_EXCEPTION(rowMap_->isNodeLocalElement(LocalRow) == false, std::runtime_error,
01297         Teuchos::typeName(*this) << "::getLocalRowView(LocalRow): LocalRow (== " << LocalRow << ") is not valid on this node.");
01298     Teuchos::ArrayRCP<const LocalOrdinal> ret;
01299     RowInfo sizeInfo = getFullLocalView(LocalRow, ret);
01300     if (ret != Teuchos::null) {
01301       ret = ret.persistingView(0,sizeInfo.numEntries);
01302     }
01303     return ret;
01304   }
01305 
01306 
01309   template <class LocalOrdinal, class GlobalOrdinal, class Node>
01310   Teuchos::ArrayRCP<const GlobalOrdinal> CrsGraph<LocalOrdinal,GlobalOrdinal,Node>::getGlobalRowView(GlobalOrdinal GlobalRow) const {
01311     TEST_FOR_EXCEPTION(isLocallyIndexed() == true, std::runtime_error,
01312         Teuchos::typeName(*this) << "::getGlobalRowView(): global indices do not exist.");
01313     const LocalOrdinal lrow = rowMap_->getLocalElement(GlobalRow);
01314     TEST_FOR_EXCEPTION(lrow == Teuchos::OrdinalTraits<LocalOrdinal>::invalid(), std::runtime_error,
01315         Teuchos::typeName(*this) << "::getGlobalRowView(GlobalRow): GlobalRow (== " << GlobalRow << ") does not belong to this node.");
01316     Teuchos::ArrayRCP<const GlobalOrdinal> ret;
01317     RowInfo sizeInfo = getFullGlobalView(lrow, ret);
01318     if (ret != Teuchos::null) {
01319       ret = ret.persistingView(0,sizeInfo.numEntries);
01320     }
01321     return ret;
01322   }
01323 
01324 
01327   template <class LocalOrdinal, class GlobalOrdinal, class Node>
01328   void CrsGraph<LocalOrdinal,GlobalOrdinal,Node>::getLocalRowCopy(LocalOrdinal LocalRow, const Teuchos::ArrayView<LocalOrdinal> &indices, size_t &NumIndices) const {
01329     // can only do this if 
01330     // * we have local indices: isLocallyIndexed()
01331     // * we are capable of producing them: isGloballyIndexed() && hasColMap()
01332     // short circuit if we aren't allocated
01333     TEST_FOR_EXCEPTION(isGloballyIndexed() == true && hasColMap() == false, std::runtime_error,
01334         Teuchos::typeName(*this) << "::getLocalRowCopy(): local indices cannot be produced.");
01335     TEST_FOR_EXCEPTION(rowMap_->isNodeLocalElement(LocalRow) == false, std::runtime_error,
01336         Teuchos::typeName(*this) << "::getLocalRowCopy(LocalRow,...): LocalRow (== " << LocalRow << ") is not valid on this node.");
01337     // use one of the view routines to get the proper view, then copy it over
01338     if (isLocallyIndexed()) {
01339       Teuchos::ArrayRCP<const LocalOrdinal> lview;
01340       RowInfo sizeInfo = getFullLocalView(LocalRow, lview);
01341       NumIndices = sizeInfo.numEntries;
01342       TEST_FOR_EXCEPTION((size_t)indices.size() < NumIndices, std::runtime_error,
01343           Teuchos::typeName(*this) << "::getLocalRowCopy(): specified storage (size==" << indices.size() 
01344           << ") is not large enough to hold all entries for this row (NumIndices == " << NumIndices << ").");
01345       if (NumIndices > 0) {
01346         std::copy( lview.begin(), lview.begin() + NumIndices, indices.begin());
01347       }
01348       lview = Teuchos::null;
01349     }
01350     else if (isGloballyIndexed()) {
01351       Teuchos::ArrayRCP<const GlobalOrdinal> gview;
01352       RowInfo sizeInfo = getFullGlobalView(LocalRow, gview);
01353       NumIndices = sizeInfo.numEntries;
01354       TEST_FOR_EXCEPTION((size_t)indices.size() < NumIndices, std::runtime_error,
01355           Teuchos::typeName(*this) << "::getLocalRowCopy(): specified storage (size==" << indices.size() 
01356           << ") is not large enough to hold all entries for this row (NumIndices == " << NumIndices << ").");
01357       for (size_t j=0; j < NumIndices; ++j) {
01358         indices[j] = colMap_->getLocalElement(gview[j]);
01359       }
01360       gview = Teuchos::null;
01361     }
01362     else {
01363 #ifdef HAVE_TPETRA_DEBUG
01364       // should have fallen in one of the above
01365       TEST_FOR_EXCEPTION( indicesAreAllocated() == true, std::logic_error, 
01366           Teuchos::typeName(*this) << "::getLocalRowCopy(): Internal logic error. Please contact Tpetra team.");
01367 #endif
01368       NumIndices = 0;
01369     }
01370     return;
01371   }
01372 
01373 
01376   template <class LocalOrdinal, class GlobalOrdinal, class Node>
01377   void CrsGraph<LocalOrdinal,GlobalOrdinal,Node>::getGlobalRowCopy(GlobalOrdinal GlobalRow, const Teuchos::ArrayView<GlobalOrdinal> &indices, size_t &NumIndices) const {
01378     // we either currently store global indices, or we have a column map with which to transcribe our local indices for the user
01379     const LocalOrdinal lrow = rowMap_->getLocalElement(GlobalRow);
01380     TEST_FOR_EXCEPTION(lrow == Teuchos::OrdinalTraits<LocalOrdinal>::invalid(), std::runtime_error,
01381         Teuchos::typeName(*this) << "::getGlobalRowCopy(GlobalRow,...): GlobalRow (== " << GlobalRow << ") does not belong to this node.");
01382     // use one of the view routines to get the proper view, then copy it over
01383     if (isLocallyIndexed()) {
01384       Teuchos::ArrayRCP<const LocalOrdinal> lview;
01385       RowInfo sizeInfo = getFullLocalView(static_cast<size_t>(lrow), lview);
01386       NumIndices = sizeInfo.numEntries;
01387       TEST_FOR_EXCEPTION((size_t)indices.size() < NumIndices, std::runtime_error,
01388           Teuchos::typeName(*this) << "::getGlobalRowCopy(): specified storage (size==" << indices.size() 
01389           << ") is not large enough to hold all entries for this row (NumIndices == " << NumIndices << ").");
01390       // copy and convert
01391       for (size_t j=0; j < NumIndices; ++j) {
01392         indices[j] = colMap_->getGlobalElement(lview[j]);
01393       }
01394       lview = Teuchos::null;
01395     }
01396     else if (isGloballyIndexed()) {
01397       Teuchos::ArrayRCP<const GlobalOrdinal> gview;
01398       RowInfo sizeInfo = getFullGlobalView(static_cast<size_t>(lrow), gview);
01399       NumIndices = sizeInfo.numEntries;
01400       TEST_FOR_EXCEPTION((size_t)indices.size() < NumIndices, std::runtime_error,
01401           Teuchos::typeName(*this) << "::getGlobalRowCopy(): specified storage (size==" << indices.size() 
01402           << ") is not large enough to hold all entries for this row (NumIndices == " << NumIndices << ").");
01403       std::copy(gview.begin(), gview.begin() + NumIndices, indices.begin());
01404       gview = Teuchos::null;
01405     }
01406     return;
01407   }
01408 
01409 
01412   template <class LocalOrdinal, class GlobalOrdinal, class Node>
01413   void CrsGraph<LocalOrdinal,GlobalOrdinal,Node>::insertLocalIndices(LocalOrdinal lrow, const Teuchos::ArrayView<const LocalOrdinal> &indices) {
01414     using Teuchos::ArrayRCP;
01415     TEST_FOR_EXCEPTION(isStorageOptimized() == true, std::runtime_error,
01416         Teuchos::typeName(*this) << "::insertLocalIndices(): cannot insert new indices after optimizeStorage() has been called.");
01417     TEST_FOR_EXCEPTION(isGloballyIndexed() == true, std::runtime_error,
01418         Teuchos::typeName(*this) << "::insertLocalIndices(): graph indices are global; use insertGlobalIndices().");
01419     TEST_FOR_EXCEPTION(hasColMap() == false, std::runtime_error,
01420         Teuchos::typeName(*this) << "::insertLocalIndices(): cannot insert local indices without a column map.");
01421     TEST_FOR_EXCEPTION(rowMap_->isNodeLocalElement(lrow) == false, std::runtime_error,
01422         Teuchos::typeName(*this) << "::insertLocalIndices(): row does not belong to this node.");
01423     if (indicesAreAllocated() == false) {
01424       allocateIndices(AllocateLocal);
01425 #ifdef HAVE_TPETRA_DEBUG
01426       TEST_FOR_EXCEPTION(indicesAreAllocated() == false, std::logic_error, 
01427           Teuchos::typeName(*this) << "::insertLocalIndices(): Internal logic error. Please contact Tpetra team.");
01428 #endif
01429     }
01430     //
01431     indicesAreSorted_ = false;
01432     noRedundancies_ = false;
01433     clearGlobalConstants();
01434     //
01435     // add to allocated space
01436     ArrayRCP<LocalOrdinal> rowview;
01437     RowInfo sizeInfo = getFullLocalViewNonConst(lrow, rowview);
01438     const size_t newSize = sizeInfo.numEntries + indices.size();
01439     if (newSize > sizeInfo.allocSize) {
01440       TEST_FOR_EXCEPTION(getProfileType() == StaticProfile, std::runtime_error,
01441           Teuchos::typeName(*this) << "::insertLocalIndices(): new indices exceed statically allocated graph structure.");
01442       TPETRA_EFFICIENCY_WARNING(true, std::runtime_error,
01443           "::insertLocalIndices(): Pre-allocated space has been exceeded, requiring new allocation. To improve efficiency, suggest larger allocation.");
01444       // update allocation only as much as necessary
01445       updateLocalAllocation(lrow,newSize);
01446       // get new view; inefficient, but acceptible in this already inefficient case
01447       sizeInfo = getFullLocalViewNonConst(lrow, rowview);
01448     }
01449     // get pointers to row allocation
01450     ArrayRCP<LocalOrdinal> rowptr = rowview + sizeInfo.numEntries;
01451     // check the local indices against the column map; only add ones that are defined
01452     typename Teuchos::ArrayView<const LocalOrdinal>::iterator srcind = indices.begin();
01453     while (srcind != indices.end()) {
01454       if (colMap_->isNodeLocalElement(*srcind)) {
01455         (*rowptr++) = (*srcind);
01456       }
01457       ++srcind;
01458     }
01459     numEntriesPerRow_[lrow] = rowptr - rowview;
01460     rowptr = Teuchos::null;
01461     rowview = Teuchos::null;
01462     // checkInternalState();
01463   }
01464 
01467   template <class LocalOrdinal, class GlobalOrdinal, class Node>
01468   void CrsGraph<LocalOrdinal,GlobalOrdinal,Node>::removeLocalIndices(LocalOrdinal lrow) {
01469     TEST_FOR_EXCEPTION(isStorageOptimized() == true, std::runtime_error,
01470         Teuchos::typeName(*this) << "::removeLocalIndices(): cannot remove indices after optimizeStorage() has been called.");
01471     TEST_FOR_EXCEPTION(isGloballyIndexed() == true, std::runtime_error,
01472         Teuchos::typeName(*this) << "::removeLocalIndices(): graph indices are global; use removeGlobalIndices().");
01473     TEST_FOR_EXCEPTION(rowMap_->isNodeLocalElement(lrow) == false, std::runtime_error,
01474         Teuchos::typeName(*this) << "::removeLocalIndices(): row does not belong to this node.");
01475     if (indicesAreAllocated() == false) {
01476       allocateIndices(AllocateLocal);
01477 #ifdef HAVE_TPETRA_DEBUG
01478       TEST_FOR_EXCEPTION(indicesAreAllocated() == false, std::logic_error, 
01479           Teuchos::typeName(*this) << "::removeLocalIndices(): Internal logic error. Please contact Tpetra team.");
01480 #endif
01481     }
01482     //
01483     clearGlobalConstants();
01484     //
01485     RowInfo sizeInfo = getRowInfo(lrow);
01486     if (sizeInfo.allocSize > 0 && numEntriesPerRow_ != Teuchos::null) {
01487       numEntriesPerRow_[lrow] = 0;
01488     }
01489   }
01490 
01491 
01494   template <class LocalOrdinal, class GlobalOrdinal, class Node>
01495   void CrsGraph<LocalOrdinal,GlobalOrdinal,Node>::insertGlobalIndices(GlobalOrdinal grow, const Teuchos::ArrayView<const GlobalOrdinal> &indices) {
01496     using Teuchos::OrdinalTraits;
01497     using Teuchos::ArrayRCP;
01498     TEST_FOR_EXCEPTION(isLocallyIndexed() == true, std::runtime_error,
01499         Teuchos::typeName(*this) << "::insertGlobalIndices(): graph indices are local; use insertLocalIndices().");
01500     if (indicesAreAllocated() == false) {
01501       allocateIndices(AllocateGlobal);
01502 #ifdef HAVE_TPETRA_DEBUG
01503       TEST_FOR_EXCEPTION(indicesAreAllocated() == false, std::logic_error, 
01504           Teuchos::typeName(*this) << "::insertGlobalIndices(): Internal logic error. Please contact Tpetra team.");
01505 #endif
01506     }
01507     // 
01508     indicesAreSorted_ = false;
01509     noRedundancies_ = false;
01510     clearGlobalConstants();
01511     //
01512     const LocalOrdinal lrow = rowMap_->getLocalElement(grow);
01513     if (lrow != OrdinalTraits<LocalOrdinal>::invalid()) {
01514       //
01515       // add to allocated space
01516       ArrayRCP<GlobalOrdinal> rowview;
01517       RowInfo sizeInfo = getFullGlobalViewNonConst(static_cast<size_t>(lrow), rowview);
01518       const size_t newSize = sizeInfo.numEntries + indices.size();
01519       if (newSize > sizeInfo.allocSize) {
01520         TEST_FOR_EXCEPTION(getProfileType() == StaticProfile, std::runtime_error,
01521             Teuchos::typeName(*this) << "::insertGlobalIndices(): new indices exceed statically allocated graph structure.");
01522         TPETRA_EFFICIENCY_WARNING(true,std::runtime_error,
01523             "::insertGlobalIndices(): Pre-allocated space has been exceeded, requiring new allocation. To improve efficiency, suggest larger allocation.");
01524         // update allocation
01525         updateGlobalAllocation(static_cast<size_t>(lrow),newSize);
01526         // get new view; inefficient, but acceptible in this already inefficient case
01527         sizeInfo = getFullGlobalViewNonConst(static_cast<size_t>(lrow), rowview);
01528       }
01529       ArrayRCP<GlobalOrdinal> rowptr = rowview + sizeInfo.numEntries;
01530       if (hasColMap()) {
01531         // check the global indices against the column map; only add ones that are defined
01532         typename Teuchos::ArrayView<const GlobalOrdinal>::iterator srcind = indices.begin();
01533         while (srcind != indices.end()) {
01534           if (colMap_->isNodeGlobalElement(*srcind)) {
01535             (*rowptr++) = (*srcind);
01536           }
01537           ++srcind;
01538         }
01539         numEntriesPerRow_[lrow] = rowptr - rowview;
01540       }
01541       else {
01542         std::copy( indices.begin(), indices.end(), rowptr );
01543         numEntriesPerRow_[lrow] += indices.size();
01544       }
01545       rowptr = Teuchos::null;
01546       rowview = Teuchos::null;
01547     }
01548     else {
01549       // a nonlocal row
01550       for (typename Teuchos::ArrayView<const GlobalOrdinal>::iterator i=indices.begin(); i != indices.end(); ++i) {
01551         nonlocals_[grow].push_back(*i);
01552       }
01553     }
01554     // checkInternalState();
01555   }
01556 
01557 
01560   template <class LocalOrdinal, class GlobalOrdinal, class Node>
01561   void CrsGraph<LocalOrdinal,GlobalOrdinal,Node>::globalAssemble() {
01562     using Teuchos::arcp;
01563     using Teuchos::Array;
01564     using Teuchos::ArrayView;
01565     using Teuchos::ArrayRCP;
01566     using Teuchos::RCP;
01567     using Teuchos::rcp;
01568     using Teuchos::outArg;
01569     using std::deque;
01570     using std::pair;
01571     using std::make_pair;
01572     using Teuchos::outArg;
01573     typedef typename std::map<GlobalOrdinal,std::deque<GlobalOrdinal> >::const_iterator NLITER;
01574     int numImages = Teuchos::size(*getComm());
01575     int myImageID = Teuchos::rank(*getComm());
01576     // Determine if any nodes have global entries to share
01577     {
01578       size_t MyNonlocals = nonlocals_.size(), MaxGlobalNonlocals;
01579       Teuchos::reduceAll<int,size_t>(*getComm(),Teuchos::REDUCE_MAX,MyNonlocals,
01580         outArg(MaxGlobalNonlocals));
01581       if (MaxGlobalNonlocals == 0) return;  // no entries to share
01582     }
01583 
01584     // compute a list of NLRs from nonlocals_ and use it to compute:
01585     //      IdsAndRows: a vector of (id,row) pairs
01586     //          NLR2Id: a map from NLR to the Id that owns it
01587     // globalNeighbors: a global graph of connectivity between images: globalNeighbors(i,j) indicates that j sends to i
01588     //         sendIDs: a list of all images I send to
01589     //         recvIDs: a list of all images I receive from (constructed later)
01590     Array<pair<int,GlobalOrdinal> > IdsAndRows;
01591     std::map<GlobalOrdinal,int> NLR2Id;
01592     Teuchos::SerialDenseMatrix<int,char> globalNeighbors;
01593     Array<int> sendIDs, recvIDs;
01594     {
01595       // nonlocals_ contains the entries we are holding for all non-local rows
01596       // we want a list of the rows for which we have data
01597       Array<GlobalOrdinal> NLRs;
01598       std::set<GlobalOrdinal> setOfRows;
01599       for (NLITER iter = nonlocals_.begin(); iter != nonlocals_.end(); ++iter) {
01600         setOfRows.insert(iter->first);
01601       }
01602       // copy the elements in the set into an Array
01603       NLRs.resize(setOfRows.size());
01604       std::copy(setOfRows.begin(), setOfRows.end(), NLRs.begin());
01605 
01606       // get a list of ImageIDs for the non-local rows (NLRs)
01607       Array<int> NLRIds(NLRs.size());
01608       {
01609         LookupStatus stat = rowMap_->getRemoteIndexList(NLRs(),NLRIds());
01610         char lclerror = ( stat == IDNotPresent ? 1 : 0 );
01611         char gblerror;
01612         Teuchos::reduceAll(*getComm(),Teuchos::REDUCE_MAX,lclerror,outArg(gblerror));
01613         TEST_FOR_EXCEPTION(gblerror != 0, std::runtime_error,
01614             Teuchos::typeName(*this) << "::globalAssemble(): non-local entries correspond to invalid rows.");
01615       }
01616 
01617       // build up a list of neighbors, as well as a map between NLRs and Ids
01618       // localNeighbors[i] != 0 iff I have data to send to image i
01619       // put NLRs,Ids into an array of pairs
01620       IdsAndRows.reserve(NLRs.size());
01621       Array<char> localNeighbors(numImages,0);
01622       typename Array<GlobalOrdinal>::const_iterator nlr;
01623       typename Array<int>::const_iterator id;
01624       for (nlr = NLRs.begin(), id = NLRIds.begin();
01625            nlr != NLRs.end(); ++nlr, ++id) {
01626         NLR2Id[*nlr] = *id;
01627         localNeighbors[*id] = 1;
01628         IdsAndRows.push_back(make_pair<int,GlobalOrdinal>(*id,*nlr));
01629       }
01630       for (int j=0; j<numImages; ++j) {
01631         if (localNeighbors[j]) {
01632           sendIDs.push_back(j);
01633         }
01634       }
01635       // sort IdsAndRows, by Ids first, then rows
01636       std::sort(IdsAndRows.begin(),IdsAndRows.end());
01637       // gather from other nodes to form the full graph
01638       globalNeighbors.shapeUninitialized(numImages,numImages);
01639       Teuchos::gatherAll(*getComm(),numImages,localNeighbors.getRawPtr(),numImages*numImages,globalNeighbors.values());
01640       // globalNeighbors at this point contains (on all images) the
01641       // connectivity between the images. 
01642       // globalNeighbors(i,j) != 0 means that j sends to i/that i receives from j
01643     }
01644 
01646     // FIGURE OUT WHO IS SENDING TO WHOM AND HOW MUCH
01647     // DO THIS IN THE PROCESS OF PACKING ALL OUTGOING DATA ACCORDING TO DESTINATION ID
01649 
01650     // loop over all columns to know from which images I can expect to receive something
01651     for (int j=0; j<numImages; ++j) {
01652       if (globalNeighbors(myImageID,j)) {
01653         recvIDs.push_back(j);
01654       }
01655     }
01656     const size_t numRecvs = recvIDs.size();
01657 
01658     // we know how many we're sending to already
01659     // form a contiguous list of all data to be sent
01660     // track the number of entries for each ID
01661     Array<pair<GlobalOrdinal,GlobalOrdinal> > IJSendBuffer;
01662     Array<size_t> sendSizes(sendIDs.size(), 0);
01663     size_t numSends = 0;
01664     for (typename Array<pair<int,GlobalOrdinal> >::const_iterator IdAndRow = IdsAndRows.begin();
01665          IdAndRow != IdsAndRows.end(); ++IdAndRow) {
01666       int            id = IdAndRow->first;
01667       GlobalOrdinal row = IdAndRow->second;
01668       // have we advanced to a new send?
01669       if (sendIDs[numSends] != id) {
01670         numSends++;
01671         TEST_FOR_EXCEPTION(sendIDs[numSends] != id, std::logic_error, Teuchos::typeName(*this) << "::globalAssemble(): internal logic error. Contact Tpetra team.");
01672       }
01673       // copy data for row into contiguous storage
01674       for (typename deque<GlobalOrdinal>::const_iterator j = nonlocals_[row].begin(); j != nonlocals_[row].end(); ++j)
01675       {
01676         IJSendBuffer.push_back( pair<GlobalOrdinal,GlobalOrdinal>(row,*j) );
01677         sendSizes[numSends]++;
01678       }
01679     }
01680     if (IdsAndRows.size() > 0) {
01681       numSends++; // one last increment, to make it a count instead of an index
01682     }
01683     TEST_FOR_EXCEPTION(Teuchos::as<typename Array<int>::size_type>(numSends) != sendIDs.size(), std::logic_error, Teuchos::typeName(*this) << "::globalAssemble(): internal logic error. Contact Tpetra team.");
01684 
01685     // don't need this data anymore
01686     nonlocals_.clear();
01687 
01689     // TRANSMIT SIZE INFO BETWEEN SENDERS AND RECEIVERS
01691     // perform non-blocking sends: send sizes to our recipients
01692     Array<RCP<Teuchos::CommRequest> > sendRequests;
01693     for (size_t s=0; s < numSends ; ++s) {
01694       // we'll fake the memory management, because all communication will be local to this method and the scope of our data
01695       sendRequests.push_back( Teuchos::isend<int,size_t>(*getComm(),rcp<size_t>(&sendSizes[s],false),sendIDs[s]) );
01696     }
01697     // perform non-blocking receives: receive sizes from our senders
01698     Array<RCP<Teuchos::CommRequest> > recvRequests;
01699     Array<size_t> recvSizes(numRecvs);
01700     for (size_t r=0; r < numRecvs; ++r) {
01701       // we'll fake the memory management, because all communication will be local to this method and the scope of our data
01702       recvRequests.push_back( Teuchos::ireceive(*getComm(),rcp(&recvSizes[r],false),recvIDs[r]) );
01703     }
01704     // wait on all 
01705     if (!sendRequests.empty()) {
01706       Teuchos::waitAll(*getComm(),sendRequests());
01707     }
01708     if (!recvRequests.empty()) {
01709       Teuchos::waitAll(*getComm(),recvRequests());
01710     }
01711     Teuchos::barrier(*getComm());
01712     sendRequests.clear();
01713     recvRequests.clear();
01714 
01716     // NOW SEND/RECEIVE ALL ROW DATA
01718     // from the size info, build the ArrayViews into IJSendBuffer
01719     Array<ArrayView<pair<GlobalOrdinal,GlobalOrdinal> > > sendBuffers(numSends,Teuchos::null);
01720     {
01721       size_t cur = 0;
01722       for (size_t s=0; s<numSends; ++s) {
01723         sendBuffers[s] = IJSendBuffer(cur,sendSizes[s]);
01724         cur += sendSizes[s];
01725       }
01726     }
01727     // perform non-blocking sends
01728     for (size_t s=0; s < numSends ; ++s)
01729     {
01730       // we'll fake the memory management, because all communication will be local to this method and the scope of our data
01731       ArrayRCP<pair<GlobalOrdinal,GlobalOrdinal> > tmparcp = arcp(sendBuffers[s].getRawPtr(),0,sendBuffers[s].size(),false);
01732       sendRequests.push_back( Teuchos::isend<int,pair<GlobalOrdinal,GlobalOrdinal> >(*getComm(),tmparcp,sendIDs[s]) );
01733     }
01734     // calculate amount of storage needed for receives
01735     // setup pointers for the receives as well
01736     size_t totalRecvSize = std::accumulate(recvSizes.begin(),recvSizes.end(),0);
01737     Array<pair<GlobalOrdinal,GlobalOrdinal> > IJRecvBuffer(totalRecvSize);
01738     // from the size info, build the ArrayViews into IJRecvBuffer
01739     Array<ArrayView<pair<GlobalOrdinal,GlobalOrdinal> > > recvBuffers(numRecvs,Teuchos::null);
01740     {
01741       size_t cur = 0;
01742       for (size_t r=0; r<numRecvs; ++r) {
01743         recvBuffers[r] = IJRecvBuffer(cur,recvSizes[r]);
01744         cur += recvSizes[r];
01745       }
01746     }
01747     // perform non-blocking recvs
01748     for (size_t r=0; r < numRecvs ; ++r) {
01749       // we'll fake the memory management, because all communication will be local to this method and the scope of our data
01750       ArrayRCP<pair<GlobalOrdinal,GlobalOrdinal> > tmparcp = arcp(recvBuffers[r].getRawPtr(),0,recvBuffers[r].size(),false);
01751       recvRequests.push_back( Teuchos::ireceive(*getComm(),tmparcp,recvIDs[r]) );
01752     }
01753     // perform waits
01754     if (!sendRequests.empty()) {
01755       Teuchos::waitAll(*getComm(),sendRequests());
01756     }
01757     if (!recvRequests.empty()) {
01758       Teuchos::waitAll(*getComm(),recvRequests());
01759     }
01760     Teuchos::barrier(*getComm());
01761     sendRequests.clear();
01762     recvRequests.clear();
01763 
01765     // NOW PROCESS THE RECEIVED ROW DATA
01767     // TODO: instead of adding one entry at a time, add one row at a time.
01768     //       this requires resorting; they arrived sorted by sending node, so that entries could be non-contiguous if we received
01769     //       multiple entries for a particular row from different processors.
01770     //       it also requires restoring the data, which may make it not worth the trouble.
01771     for (typename Array<pair<GlobalOrdinal,GlobalOrdinal> >::const_iterator ij = IJRecvBuffer.begin(); ij != IJRecvBuffer.end(); ++ij) {
01772       insertGlobalIndices(ij->first, Teuchos::tuple<GlobalOrdinal>(ij->second));
01773     }
01774     checkInternalState();
01775   }
01776 
01777 
01780   template <class LocalOrdinal, class GlobalOrdinal, class Node>
01781   void CrsGraph<LocalOrdinal,GlobalOrdinal,Node>::fillComplete(OptimizeOption os) {
01782     fillComplete(rowMap_,rowMap_,os);
01783   }
01784 
01785 
01788   template <class LocalOrdinal, class GlobalOrdinal, class Node>
01789   void CrsGraph<LocalOrdinal,GlobalOrdinal,Node>::fillComplete(const Teuchos::RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > &domainMap, 
01790                                                                const Teuchos::RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > &rangeMap, 
01791                                                                OptimizeOption os) {
01792     using Teuchos::Array;
01793     using Teuchos::ArrayRCP;
01794     using Teuchos::outArg;
01795 
01796     domainMap_ = domainMap;
01797     rangeMap_  = rangeMap;
01798 
01799     if (indicesAreAllocated() == false) {
01800       // must allocate global, because we have no column map
01801       allocateIndices( AllocateGlobal );
01802     }
01803 
01804     if (Teuchos::size(*getComm()) > 1) {
01805       globalAssemble();
01806     }
01807     else {
01808       TEST_FOR_EXCEPTION(nonlocals_.size() > 0, std::runtime_error,
01809           Teuchos::typeName(*this) << "::fillComplete(): cannot have non-local entries on a serial run. Invalid entries were submitted to the CrsMatrix.");
01810     }
01811 
01812     makeIndicesLocal(domainMap, rangeMap);  // create column map, allocate/reapportion space for local indices, transform global indices to local indices
01813     sortIndices();                          // sort local indices
01814     removeRedundantIndices();               // remove redundant indices, count non-zero entries, count diagonals, check triangularity
01815     makeImportExport();                     // create import and export object
01816 
01817     // compute global constants using computed local constants
01818     if (haveGlobalConstants_ == false) {
01819       global_size_t lcl[2], gbl[2];
01820       lcl[0] = nodeNumEntries_;
01821       lcl[1] = nodeNumDiags_;
01822       Teuchos::reduceAll<int,global_size_t>(*getComm(),Teuchos::REDUCE_SUM,2,lcl,gbl);
01823       globalNumEntries_ = gbl[0]; 
01824       globalNumDiags_   = gbl[1];
01825       Teuchos::reduceAll<int,global_size_t>(*getComm(),Teuchos::REDUCE_MAX,nodeMaxNumRowEntries_,
01826         outArg(globalMaxNumRowEntries_));
01827       haveGlobalConstants_ = true;
01828     }
01829 
01830     // mark transformation as successfully completed
01831     fillComplete_ = true;
01832     checkInternalState();
01833 
01834     if (os == DoOptimizeStorage) {
01835       // optimizeStorage will call fillLocalGraph()
01836       optimizeStorage();
01837     }
01838     else {
01839       // clear all views, so that changes can be written
01840       view_lclInds1D_ = Teuchos::null;
01841       view_rowOffsets_ = Teuchos::null;
01842       view_lclInds2D_ = Teuchos::null;
01843       fillLocalGraph();
01844     }
01845   }
01846 
01847 
01850   template <class LocalOrdinal, class GlobalOrdinal, class Node>
01851   void CrsGraph<LocalOrdinal,GlobalOrdinal,Node>::makeIndicesLocal(
01852                                       const Teuchos::RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > &domainMap, 
01853                                       const Teuchos::RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > &rangeMap) {
01854     using Teuchos::ArrayRCP;
01855     using Teuchos::NullIteratorTraits;
01856     // All nodes must be in the same index state.
01857     // Update index state by checking isLocallyIndexed/Global on all nodes
01858     computeIndexState(); 
01859     TEST_FOR_EXCEPTION(isLocallyIndexed() && isGloballyIndexed(), std::logic_error,
01860         Teuchos::typeName(*this) << "::makeIndicesLocal(): indices are marked as both global and local.");
01861     // If user has not prescribed column map, create one from indices
01862     makeColMap(domainMap, rangeMap);
01863     // Transform indices to local index space
01864     const size_t nlrs = getNodeNumRows();
01865     Teuchos::RCP<Node> node = lclGraph_.getNode();
01866     // 
01867     if (isGloballyIndexed() && nlrs > 0) {
01868       // allocate data for local indices
01869       if (nodeNumAllocated_ > 0) {
01870         if (getProfileType() == StaticProfile) {
01871           if (view_gblInds1D_ == Teuchos::null) {
01872             KOKKOS_NODE_TRACE("CrsGraph::makeIndicesLocal()")
01873             view_gblInds1D_ = node->template viewBufferNonConst<GlobalOrdinal>(Kokkos::ReadWrite, pbuf_gblInds1D_.size(), pbuf_gblInds1D_);
01874           }
01875           // do the conversion in situ. this must be done from front to back.
01876           view_lclInds1D_ = Teuchos::arcp_reinterpret_cast<LocalOrdinal>(view_gblInds1D_).persistingView(0,nodeNumAllocated_);
01877           for (size_t r=0; r < getNodeNumRows(); ++r) {
01878             const size_t offset   = view_rowOffsets_[r],
01879                          numentry = numEntriesPerRow_[r];
01880             for (size_t j=0; j<numentry; ++j) {
01881               GlobalOrdinal gid = view_gblInds1D_[offset + j];
01882               LocalOrdinal  lid = colMap_->getLocalElement(gid);
01883               view_lclInds1D_[offset + j] = lid;
01884 #ifdef HAVE_TPETRA_DEBUG
01885               TEST_FOR_EXCEPTION(view_lclInds1D_[offset + j] == Teuchos::OrdinalTraits<LocalOrdinal>::invalid(), std::logic_error,
01886                   Teuchos::typeName(*this) << ": Internal error in fillComplete(). Please contact Tpetra team.");
01887 #endif
01888             }
01889           }
01890           // reinterpret the the compute buffer as LocalOrdinal; keep the original size
01891           pbuf_lclInds1D_ = Teuchos::arcp_reinterpret_cast<LocalOrdinal>(pbuf_gblInds1D_).persistingView(0,nodeNumAllocated_);
01892           // order of deletion doesn't matter; view_lclInds1D_ points to view_gblInds2D_, so there will be no copy-back
01893           pbuf_gblInds1D_ = Teuchos::null;
01894           view_gblInds1D_ = Teuchos::null;
01895         }
01896         else {  // getProfileType() == DynamicProfile
01897           pbuf_lclInds2D_ = Teuchos::arcp< Teuchos::ArrayRCP<LocalOrdinal> >(nlrs);
01898           // if we have views, then make views
01899           if (view_gblInds2D_ != Teuchos::null) {
01900             view_lclInds2D_ = Kokkos::ArrayOfViewsHelper<Node>::template getArrayOfNonConstViews<LocalOrdinal>(node, Kokkos::WriteOnly, pbuf_lclInds2D_);
01901           }
01902           for (size_t r=0; r < getNodeNumRows(); ++r) {
01903             if (pbuf_gblInds2D_[r] != Teuchos::null) {
01904               const size_t rna = pbuf_gblInds2D_[r].size();
01905               ArrayRCP<GlobalOrdinal> view_ginds;
01906               ArrayRCP< LocalOrdinal> view_linds;
01907               if (view_gblInds2D_ == Teuchos::null) {
01908                 KOKKOS_NODE_TRACE("CrsGraph::makeIndicesLocal()")
01909                 view_ginds = node->template viewBufferNonConst<GlobalOrdinal>(Kokkos::ReadWrite,rna,pbuf_gblInds2D_[r]);
01910               }
01911               else {
01912                 view_ginds = view_gblInds2D_[r];
01913               }
01914               // do the conversion in situ. this must be done from front to back.
01915               view_linds = Teuchos::arcp_reinterpret_cast<LocalOrdinal>(view_ginds).persistingView(0,rna);
01916               const size_t numentry = numEntriesPerRow_[r];
01917               for (size_t j=0; j < numentry; ++j) {
01918                 GlobalOrdinal gid = view_ginds[j];
01919                 LocalOrdinal  lid = colMap_->getLocalElement(gid);
01920                 view_linds[j] = lid;
01921 #ifdef HAVE_TPETRA_DEBUG
01922                 TEST_FOR_EXCEPTION(view_linds[j] == Teuchos::OrdinalTraits<LocalOrdinal>::invalid(), std::logic_error,
01923                     Teuchos::typeName(*this) << ": Internal error in makeIndicesLocal(). Please contact Tpetra team.");
01924 #endif
01925               }
01926               if (view_lclInds2D_ != Teuchos::null) { 
01927                 view_lclInds2D_[r] = view_linds;
01928               }
01929               view_linds = Teuchos::null;
01930               view_ginds = Teuchos::null;
01931               // reinterpret the the compute buffer as LocalOrdinal; keep the original size
01932               pbuf_lclInds2D_[r] = Teuchos::arcp_reinterpret_cast<LocalOrdinal>(pbuf_gblInds2D_[r]).persistingView(0,rna);
01933               pbuf_gblInds2D_[r] = Teuchos::null;
01934             }
01935           }
01936           pbuf_gblInds2D_ = Teuchos::null;
01937           view_gblInds2D_ = Teuchos::null;
01938         }
01939       }
01940       // don't set this unless we actually did something
01941       indicesAreLocal_  = true;
01942     }
01943     indicesAreGlobal_ = false;
01944     checkInternalState();
01945   }
01946 
01947 
01950   template <class LocalOrdinal, class GlobalOrdinal, class Node>
01951   void CrsGraph<LocalOrdinal,GlobalOrdinal,Node>::computeIndexState() {
01952     char myIndices[2] = {0,0};
01953     if (indicesAreLocal_)  myIndices[0] = 1;
01954     if (indicesAreGlobal_) myIndices[1] = 1;
01955     char allIndices[2];
01956     Teuchos::reduceAll(*getComm(),Teuchos::REDUCE_MAX,2,myIndices,allIndices);
01957     indicesAreLocal_  = (allIndices[0]==1);  // If indices are local on one PE, should be local on all
01958     indicesAreGlobal_ = (allIndices[1]==1);  // If indices are global on one PE should be local on all
01959   }
01960 
01961 
01964   template <class LocalOrdinal, class GlobalOrdinal, class Node>
01965   void CrsGraph<LocalOrdinal,GlobalOrdinal,Node>::sortIndices() {
01966     TEST_FOR_EXCEPT(isGloballyIndexed()==true && indicesAreAllocated_ && nodeNumAllocated_ > 0);   // this should be called only after makeIndicesLocal()
01967     if (isSorted()) return;
01968     // are there any indices to sort?
01969     if (nodeNumAllocated_ > 0) {
01970       const size_t nlrs = getNodeNumRows();
01971       for (size_t r=0; r < nlrs; ++r) {
01972         // TODO: This is slightly inefficient, because it may query pbuf_rowOffsets_ repeatadly. 
01973         //       However, it is very simple code. Consider rewriting it.
01974         Teuchos::ArrayRCP<LocalOrdinal> row_view;
01975         RowInfo info = getFullLocalViewNonConst(r, row_view);
01976         std::sort(row_view, row_view + info.numEntries);
01977       }
01978     }
01979     setSorted(true);
01980   }
01981 
01982 
01985   template <class LocalOrdinal, class GlobalOrdinal, class Node>
01986   void CrsGraph<LocalOrdinal,GlobalOrdinal,Node>::makeColMap(
01987                                       const Teuchos::RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > &domainMap, 
01988                                       const Teuchos::RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > &rangeMap) {
01989     using Teuchos::ArrayRCP;
01990     using Teuchos::Array;
01991     using Teuchos::ArrayView;
01992     using Teuchos::outArg;
01993     typedef Teuchos::OrdinalTraits<GlobalOrdinal> GOT;
01994     // 
01995     if (hasColMap()) return;
01996     const size_t nlrs = getNodeNumRows();
01997     // 
01998     computeIndexState();
01999     TEST_FOR_EXCEPTION(isLocallyIndexed() == true, std::runtime_error,
02000         Teuchos::typeName(*this) << "::makeColMap(): indices must still be global when making the column map.");
02001     // ultimate goal: list of indices for column map
02002     Array<GlobalOrdinal> myColumns;
02003     // if isGloballyIndexed() == false and isLocallyIndexed() == false, then we have nothing to do
02004     if (isGloballyIndexed() == true) {
02005       // Construct two lists of columns for which we have a non-zero
02006       // Local GIDs: Column GIDs which are locally present on the Domain map
02007       // Remote GIDs: Column GIDs which are not locally present present on the Domain map
02008       //
02009       // instead of list of local GIDs, we will use a set<LocalOrdinal> LocalGID.
02010       //
02011       const LocalOrdinal LINV = Teuchos::OrdinalTraits<LocalOrdinal>::invalid();
02012       size_t numLocalColGIDs = 0, numRemoteColGIDs = 0;
02013       //
02014       // intitial: partitioning into local and remote
02015       Teuchos::Array<char>    GIDisLocal(domainMap->getNodeNumElements(),0);
02016       std::set<GlobalOrdinal> RemoteGIDSet;
02017       for (size_t r=0; r < nlrs; ++r) {
02018         ArrayRCP<GlobalOrdinal> rowgids;
02019         RowInfo info = getFullGlobalViewNonConst(r, rowgids);
02020         if (info.numEntries > 0) {
02021           rowgids = rowgids.persistingView(0, info.numEntries);
02022           for (typename ArrayRCP<GlobalOrdinal>::iterator cind = rowgids.begin(); cind != rowgids.end(); ++cind) {
02023             GlobalOrdinal gid = (*cind);
02024             LocalOrdinal lid = domainMap->getLocalElement(gid);
02025             if (lid != LINV) {
02026               char alreadyFound = GIDisLocal[lid];
02027               if (alreadyFound == 0) {
02028                 GIDisLocal[lid] = 1;
02029                 ++numLocalColGIDs;
02030               }
02031             }
02032             else {
02033               std::pair<typename std::set<GlobalOrdinal>::iterator, bool> ip;
02034               ip = RemoteGIDSet.insert(gid);
02035               if (ip.second == true) { // gid did not exist in the set and was actually inserted
02036                 ++numRemoteColGIDs;
02037               }
02038             }
02039           }
02040         }
02041         rowgids = Teuchos::null;
02042       }
02043 
02044       // Possible short-circuit for serial scenario
02045       // If the all domain GIDs are present as column indices, then set ColMap=DomainMap
02046       // By construction, LocalGIDs \subset DomainGIDs
02047       // If we have
02048       //   * Number of remote GIDs is 0, so that ColGIDs == LocalGIDs,
02049       // and
02050       //   * Number of local GIDs is number of domain GIDs
02051       // then 
02052       //   * LocalGIDs \subset DomainGIDs && size(LocalGIDs) == size(DomainGIDs) => DomainGIDs == LocalGIDs == ColGIDs
02053       // on this node. 
02054       // We will concern ourselves only with the special case of a serial DomainMap, obliviating the need for a communication.
02055       // If 
02056       //   * DomainMap has a serial comm
02057       // then we can set Column map as Domain map and return. Benefit: this graph won't need an Import object later.
02058       // 
02059       // Note, for a serial domain map, there can be no RemoteGIDs, because there are no remote nodes. 
02060       // Likely explanations for this are:
02061       //  * user submitted erroneous column indices
02062       //  * user submitted erroneous domain map
02063       if (Teuchos::size(*domainMap->getComm()) == 1) {
02064         TEST_FOR_EXCEPTION(numRemoteColGIDs != 0, std::runtime_error,
02065             Teuchos::typeName(*this) << "::makeColMap(): Some column IDs are not in the domain map." << std::endl 
02066             << "Either these column IDs are invalid or the domain map is invalid." << std::endl
02067             << "Remember, for a rectangular matrix, the domain map must be passed to fillComplete().");
02068         if (numLocalColGIDs == domainMap->getNodeNumElements()) {
02069           colMap_ = domainMap;
02070           checkInternalState();
02071           return;
02072         }
02073       }
02074 
02075       // Now, populate myColumns() with a list of all column GIDs. 
02076       // Put local GIDs at the front: they correspond to "same" and "permuted" entries between the column map and the domain map
02077       // Put remote GIDs at the back
02078       myColumns.resize(numLocalColGIDs + numRemoteColGIDs);
02079       // get pointers into myColumns for each part
02080       ArrayView<GlobalOrdinal> LocalColGIDs, RemoteColGIDs;
02081       LocalColGIDs  = myColumns(0,numLocalColGIDs);
02082       RemoteColGIDs = myColumns(numLocalColGIDs,numRemoteColGIDs);
02083 
02084       // Copy the Remote GIDs into myColumns
02085       std::copy(RemoteGIDSet.begin(), RemoteGIDSet.end(), RemoteColGIDs.begin());
02086       // We will make a list of corresponding node IDs
02087       Array<int> RemoteImageIDs(numRemoteColGIDs);
02088       // Lookup the Remote Nodes IDs in the Domain map
02089       {
02090         LookupStatus stat = domainMap->getRemoteIndexList(RemoteColGIDs, RemoteImageIDs());
02091         // This error check is crucial: this tells us that one of the remote indices was not present in the domain map. 
02092         // This means that the Import object cannot be constructed, because of incongruity between the column map and domain map.
02093         //   * The user has made a mistake in the column indices
02094         //   * The user has made a mistake w.r.t. the domain map 
02095         // Same error message as above for serial case.
02096         char missingID_lcl = (stat == IDNotPresent ? 1 : 0);
02097         char missingID_gbl;
02098         Teuchos::reduceAll<int,char>(*getComm(),Teuchos::REDUCE_MAX,missingID_lcl,
02099           outArg(missingID_gbl));
02100         TEST_FOR_EXCEPTION(missingID_gbl == 1, std::runtime_error,
02101             Teuchos::typeName(*this) << "::makeColMap(): Some column IDs are not in the domain map." << std::endl 
02102             << "Either these column IDs are invalid or the domain map is invalid." << std::endl
02103             << "Common cause: for a rectangular matrix, the domain map must be passed to fillComplete().");
02104       }
02105       // Sort External column indices so that all columns coming from a given remote processor are contiguous
02106       // This obliviates the need for the Distributor associated with the Import from having to reorder data.
02107       sort2(RemoteImageIDs.begin(), RemoteImageIDs.end(), RemoteColGIDs.begin());
02108 
02109       // Copy the Local GIDs into myColumns. Two cases:
02110       // (1) If the number of Local column GIDs is the same as the number of Local domain GIDs, we
02111       //     can simply read the domain GIDs into the front part of ColIndices (see logic above from the serial short circuit case)
02112       // (2) We step through the GIDs of the DomainMap, checking to see if each domain GID is a column GID.
02113       //     we want to do this to maintain a consistent ordering of GIDs between the columns and the domain.
02114       ArrayView<const GlobalOrdinal> mge = domainMap->getNodeElementList();
02115       if (numLocalColGIDs == domainMap->getNodeNumElements()) {
02116         std::copy(mge.begin(), mge.end(), LocalColGIDs.begin());
02117       }
02118       else {
02119         size_t numlocalagain = 0;
02120         for (size_t i=0; i < domainMap->getNodeNumElements(); ++i) {
02121           if (GIDisLocal[i]) {
02122             LocalColGIDs[numlocalagain++] = mge[i];
02123           }
02124         }
02125         TEST_FOR_EXCEPTION(numlocalagain != numLocalColGIDs, std::logic_error,
02126             Teuchos::typeName(*this) << "::makeColMap(): Internal logic error. Please contact Tpetra team.");
02127       }
02128     }
02129     colMap_ = Teuchos::rcp(new Map<LocalOrdinal,GlobalOrdinal,Node>(GOT::invalid(), myColumns, domainMap->getIndexBase(), domainMap->getComm(), domainMap->getNode()) );
02130     checkInternalState();
02131     return;
02132   }
02133 
02134 
02137   template <class LocalOrdinal, class GlobalOrdinal, class Node>
02138   void CrsGraph<LocalOrdinal,GlobalOrdinal,Node>::removeRedundantIndices() {
02139     TEST_FOR_EXCEPT( isGloballyIndexed() != false );      // this should be called only after makeIndicesLocal()
02140     TEST_FOR_EXCEPT( isSorted() != true );                // this should be called only after sortIndices()
02141     TEST_FOR_EXCEPT( isStorageOptimized() == true );      // assumptions about available data structures
02142     if ( notRedundant() ) return;
02143     const size_t nlrs = getNodeNumRows();
02144     Teuchos::ArrayView<const GlobalOrdinal> myGlobalEntries = rowMap_->getNodeElementList();
02145     // reset all local quantities
02146     upperTriangular_ = true;
02147     lowerTriangular_ = true;
02148     nodeMaxNumRowEntries_ = 0;
02149     nodeNumEntries_       = 0;
02150     nodeNumDiags_         = 0;
02151     // indices are already sorted in each row
02152     if (indicesAreAllocated() == true && nodeNumAllocated_ > 0) {
02153       for (size_t r=0; r < nlrs; ++r) {
02154         GlobalOrdinal rgid = myGlobalEntries[r];
02155         // determine the local column index for this row, used for delimiting the diagonal
02156         const LocalOrdinal rlcid = colMap_->getLocalElement(rgid);   
02157         Teuchos::ArrayRCP<LocalOrdinal> rview;
02158         RowInfo info = getFullLocalViewNonConst(r, rview);
02159         typename Teuchos::ArrayRCP<LocalOrdinal>::iterator beg, end;
02160         beg = rview.begin();
02161         end = beg + info.numEntries;
02162         typename Teuchos::ArrayRCP<LocalOrdinal>::iterator newend, cur;
02163         newend = beg;
02164         if (beg != end) {
02165           if (rlcid == (*newend)) ++nodeNumDiags_;
02166           for (cur = beg + 1; cur != end; ++cur) {
02167             // it is unique; add it and test it for the diagonal
02168             // otherwise, it is a dup and should be ignored
02169             if (*cur != *newend) {
02170               ++newend;
02171               (*newend) = (*cur);
02172               // is this the diagonal?
02173               if (rlcid == (*newend)) ++nodeNumDiags_;
02174             }
02175           }
02176           // because of sorting, smallest column index is (*beg); it indicates upper triangularity
02177           if (Teuchos::as<size_t>(*beg) < r) upperTriangular_ = false;
02178           // because of sorting, largest column index is (*newend); it indicates lower triangularity
02179           if (r < Teuchos::as<size_t>(*newend)) lowerTriangular_ = false;
02180           // increment newend so that our range is [beg,newend) instead of [beg,newend]
02181           ++newend;
02182         }
02183         // compute num entries for this row, accumulate into nodeNumEntries_, update nodeMaxNumRowEntries_
02184         numEntriesPerRow_[r] = newend - beg;
02185         nodeMaxNumRowEntries_ = std::max( nodeMaxNumRowEntries_, numEntriesPerRow_[r] );
02186         nodeNumEntries_ += numEntriesPerRow_[r];
02187         rview = Teuchos::null;
02188       }
02189     }
02190     noRedundancies_ = true;
02191   }
02192 
02193 
02196   template <class LocalOrdinal, class GlobalOrdinal, class Node>
02197   void CrsGraph<LocalOrdinal,GlobalOrdinal,Node>::makeImportExport() {
02198     TEST_FOR_EXCEPT(hasColMap()==false); // must have column map
02199     // create import, export
02200     if (!domainMap_->isSameAs(*colMap_)) {
02201       importer_ = Teuchos::rcp( new Import<LocalOrdinal,GlobalOrdinal,Node>(domainMap_,colMap_) );
02202     }
02203     else {
02204       importer_ = Teuchos::null;
02205     }
02206     if (!rangeMap_->isSameAs(*rowMap_)) {
02207       exporter_ = Teuchos::rcp( new Export<LocalOrdinal,GlobalOrdinal,Node>(rowMap_,rangeMap_) );
02208     }
02209     else {
02210       exporter_ = Teuchos::null;
02211     }
02212   }
02213 
02214 
02217   template <class LocalOrdinal, class GlobalOrdinal, class Node>
02218   void CrsGraph<LocalOrdinal,GlobalOrdinal,Node>::optimizeStorage() {
02219     using Teuchos::ArrayRCP;
02220     // optimizeStorage will perform two functions:
02221     // 1) create a single allocation of memory
02222     // 2) pack data in that allocation
02223     // if getProfileType() == StaticProfile, then 1) has already been done
02224     // 
02225     // post-condition:
02226     //   getProfileType() == StaticProfile
02227     //   numEntriesPerRow_ == Teuchos::null; size and allocation are the same, compute by subtracting offsets
02228     if (isStorageOptimized() == true) return;
02229 
02230     if (indicesAreAllocated() == false) {
02231       // won't ever allocate now
02232       indicesAreAllocated_ = true;
02233       nodeNumAllocated_ = 0;
02234       numAllocPerRow_ = Teuchos::null;
02235       numAllocForAllRows_ = 0;
02236     }
02237 
02238     TEST_FOR_EXCEPTION(isFillComplete() == false || isSorted() == false || notRedundant() == false, std::runtime_error,
02239         Teuchos::typeName(*this) << "::optimizeStorage(): fillComplete() must be called before optimizeStorage().");
02240 
02241     Teuchos::RCP<Node> node = lclGraph_.getNode();
02242 
02243     const size_t nlrs = getNodeNumRows();
02244     if (nlrs > 0 && nodeNumAllocated_ > 0) {
02245       if (getProfileType() == DynamicProfile) {
02246         // any changes in the local indices must be committed before the copyBuffers() calls below
02247         view_lclInds2D_ = Teuchos::null;
02248         // allocate single memory block
02249         pbuf_rowOffsets_ = node->template allocBuffer<size_t>(nlrs+1);
02250         KOKKOS_NODE_TRACE("CrsGraph::optimizeStorage()")
02251         view_rowOffsets_ = node->template viewBufferNonConst<size_t>(Kokkos::WriteOnly,nlrs+1,pbuf_rowOffsets_);
02252         if (nodeNumEntries_ > 0) {
02253           pbuf_lclInds1D_ = node->template allocBuffer<LocalOrdinal>(nodeNumEntries_);
02254           ArrayRCP<LocalOrdinal> curptr = pbuf_lclInds1D_;
02255           size_t sofar = 0;
02256           for (size_t r=0; r<nlrs; ++r) {
02257             const size_t rne = numEntriesPerRow_[r];
02258             if (rne > 0) {
02259               KOKKOS_NODE_TRACE("CrsGraph::optimizeStorage()")
02260               node->template copyBuffers<LocalOrdinal>(rne, pbuf_lclInds2D_[r], curptr);
02261             }
02262             view_rowOffsets_[r] = sofar;
02263             curptr += rne;
02264             sofar += rne;
02265           }
02266           view_rowOffsets_[nlrs] = sofar;
02267 #ifdef HAVE_TPETRA_DEBUG
02268           TEST_FOR_EXCEPTION( nodeNumEntries_ != sofar, std::logic_error, 
02269               Teuchos::typeName(*this) << "::optimizeStorage(): Internal Tpetra logic error. Please contact Tpetra team.");
02270 #endif
02271         }
02272         else {
02273           std::fill(view_rowOffsets_.begin(), view_rowOffsets_.end(), 0);
02274         }
02275         // done with this buffer; we are 1D now
02276         pbuf_lclInds2D_ = Teuchos::null;
02277       }
02278       else {  // StaticProfile
02279         // any changes to the local indices must be committed before the copyBuffers() calls below
02280         view_lclInds1D_ = Teuchos::null;
02281         // storage is already allocated; just need to pack
02282         if (nodeNumEntries_ > 0) {
02283           if (view_rowOffsets_ == Teuchos::null) {
02284             KOKKOS_NODE_TRACE("CrsGraph::optimizeStorage()")
02285             view_rowOffsets_ = node->template viewBufferNonConst<size_t>(Kokkos::ReadWrite,nlrs+1,pbuf_rowOffsets_);
02286           }
02287           ArrayRCP<LocalOrdinal> curptr = pbuf_lclInds1D_,
02288                                  oldptr = pbuf_lclInds1D_;
02289           size_t sofar = 0;
02290           for (size_t r=0; r<nlrs; ++r) {
02291             const size_t rne = numEntriesPerRow_[r],
02292                           na = view_rowOffsets_[r+1] - view_rowOffsets_[r];
02293             if (curptr != oldptr) {
02294               KOKKOS_NODE_TRACE("CrsGraph::optimizeStorage()")
02295               node->template copyBuffers<LocalOrdinal>(rne, oldptr, curptr);
02296               view_rowOffsets_[r] = sofar;
02297             }
02298             sofar += rne;
02299             curptr += rne;
02300             oldptr += na;
02301           }
02302           view_rowOffsets_[nlrs] = sofar;
02303 #ifdef HAVE_TPETRA_DEBUG
02304           TEST_FOR_EXCEPTION( nodeNumEntries_ != sofar, std::logic_error, 
02305               Teuchos::typeName(*this) << "::optimizeStorage(): Internal Tpetra logic error. Please contact Tpetra team.");
02306 #endif
02307           // resize to num allocated
02308           pbuf_lclInds1D_ = pbuf_lclInds1D_.persistingView(0,sofar);
02309         }
02310       }
02311       nodeNumAllocated_ = nodeNumEntries_;
02312     }
02313     // clear it; this is duplicated now by the row pointers
02314     numEntriesPerRow_ = Teuchos::null;
02315     // if we have no storage, then we don't need this buffer at all
02316     // delete it before clearing the view below, to short-circuit the copy-back
02317     if (nodeNumAllocated_ == 0) {
02318       pbuf_rowOffsets_ = Teuchos::null;
02319       pbuf_lclInds1D_  = Teuchos::null;
02320     }
02321     // this must be cleared before the call to fillLocalGraph()
02322     view_rowOffsets_ = Teuchos::null;
02323 
02324     storageOptimized_ = true;
02325     pftype_ = StaticProfile;
02326 
02327     checkInternalState();
02328 
02329     fillLocalGraph();
02330   }
02331 
02332 
02335   template <class LocalOrdinal, class GlobalOrdinal, class Node>
02336   std::string CrsGraph<LocalOrdinal,GlobalOrdinal,Node>::description() const {
02337     std::ostringstream oss;
02338     oss << DistObject<GlobalOrdinal,LocalOrdinal,GlobalOrdinal,Node>::description();
02339     if (isFillComplete()) {
02340       oss << "{status = fill complete"
02341           << ", global rows = " << getGlobalNumRows()
02342           << ", global cols = " << getGlobalNumCols()
02343           << ", global num entries = " << getGlobalNumEntries()
02344           << "}";
02345     }
02346     else {
02347       oss << "{status = fill not complete"
02348           << ", global rows = " << getGlobalNumRows()
02349           << "}";
02350     }
02351     return oss.str();
02352   }
02353 
02354 
02357   template <class LocalOrdinal, class GlobalOrdinal, class Node>
02358   void CrsGraph<LocalOrdinal,GlobalOrdinal,Node>::describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel) const {
02359     using std::endl;
02360     using std::setw;
02361     using Teuchos::VERB_DEFAULT;
02362     using Teuchos::VERB_NONE;
02363     using Teuchos::VERB_LOW;
02364     using Teuchos::VERB_MEDIUM;
02365     using Teuchos::VERB_HIGH;
02366     using Teuchos::VERB_EXTREME;
02367     Teuchos::EVerbosityLevel vl = verbLevel;
02368     if (vl == VERB_DEFAULT) vl = VERB_LOW;
02369     Teuchos::RCP<const Teuchos::Comm<int> > comm = this->getComm();
02370     const int myImageID = comm->getRank(),
02371               numImages = comm->getSize();
02372     size_t width = 1;
02373     for (size_t dec=10; dec<getGlobalNumRows(); dec *= 10) {
02374       ++width;
02375     }
02376     width = std::max<size_t>(width,11) + 2;
02377     Teuchos::OSTab tab(out);
02378     //    none: print nothing
02379     //     low: print O(1) info from node 0
02380     //  medium: print O(P) info, num entries per node
02381     //    high: print O(N) info, num entries per row
02382     // extreme: print O(NNZ) info: print graph indices
02383     // 
02384     // for medium and higher, print constituent objects at specified verbLevel
02385     if (vl != VERB_NONE) {
02386       if (myImageID == 0) out << this->description() << std::endl; 
02387       // O(1) globals, minus what was already printed by description()
02388       if (isFillComplete() && myImageID == 0) {
02389         out << "Global number of diagonals = " << globalNumDiags_ << std::endl;
02390         out << "Global max number of entries = " << globalMaxNumRowEntries_ << std::endl;
02391       }
02392       // constituent objects
02393       if (vl == VERB_MEDIUM || vl == VERB_HIGH || vl == VERB_EXTREME) {
02394         if (myImageID == 0) out << "\nRow map: " << std::endl;
02395         rowMap_->describe(out,vl);
02396         if (colMap_ != Teuchos::null) {
02397           if (myImageID == 0) out << "\nColumn map: " << std::endl;
02398           colMap_->describe(out,vl);
02399         }
02400         if (domainMap_ != Teuchos::null) {
02401           if (myImageID == 0) out << "\nDomain map: " << std::endl;
02402           domainMap_->describe(out,vl);
02403         }
02404         if (rangeMap_ != Teuchos::null) {
02405           if (myImageID == 0) out << "\nRange map: " << std::endl;
02406           rangeMap_->describe(out,vl);
02407         }
02408       }
02409       // O(P) data
02410       if (vl == VERB_MEDIUM || vl == VERB_HIGH || vl == VERB_EXTREME) {
02411         for (int imageCtr = 0; imageCtr < numImages; ++imageCtr) {
02412           if (myImageID == imageCtr) {
02413             out << "Node ID = " << imageCtr << std::endl
02414                 << "Node number of entries = " << nodeNumEntries_ << std::endl
02415                 << "Node number of diagonals = " << nodeNumDiags_ << std::endl
02416                 << "Node max number of entries = " << nodeMaxNumRowEntries_ << std::endl;
02417             if (indicesAreAllocated()) {
02418               out << "Node number of allocated entries = " << nodeNumAllocated_ << std::endl;
02419             }
02420             else {
02421               out << "Indices are not allocated." << std::endl;
02422             }
02423           }
02424           comm->barrier();
02425           comm->barrier();
02426           comm->barrier();
02427         }
02428       }
02429       // O(N) and O(NNZ) data
02430       if (vl == VERB_HIGH || vl == VERB_EXTREME) {
02431         for (int imageCtr = 0; imageCtr < numImages; ++imageCtr) {
02432           if (myImageID == imageCtr) {
02433             out << std::setw(width) << "Node ID"
02434                 << std::setw(width) << "Global Row" 
02435                 << std::setw(width) << "Num Entries";
02436             if (vl == VERB_EXTREME) {
02437               out << "Entries";
02438             }
02439             out << std::endl;
02440             for (size_t r=0; r < getNodeNumRows(); ++r) {
02441               const size_t nE = RNNZ(r);
02442               GlobalOrdinal gid = rowMap_->getGlobalElement(r);
02443               out << std::setw(width) << myImageID 
02444                   << std::setw(width) << gid
02445                   << std::setw(width) << nE;
02446               if (vl == VERB_EXTREME) {
02447                 if (isGloballyIndexed()) {
02448                   Teuchos::ArrayRCP<const GlobalOrdinal> rowview = getGlobalRowView(gid);
02449                   for (size_t j=0; j < nE; ++j) out << rowview[j] << " ";
02450                 }
02451                 else if (isLocallyIndexed()) {
02452                   Teuchos::ArrayRCP<const LocalOrdinal> rowview = getLocalRowView(r);
02453                   for (size_t j=0; j < nE; ++j) out << colMap_->getGlobalElement(rowview[j]) << " ";
02454                 }
02455               }
02456               out << std::endl;
02457             }
02458           }
02459           comm->barrier();
02460           comm->barrier();
02461           comm->barrier();
02462         }
02463       }
02464     }
02465   }
02466 
02467 
02470   template <class LocalOrdinal, class GlobalOrdinal, class Node>
02471   bool CrsGraph<LocalOrdinal,GlobalOrdinal,Node>::checkSizes(const DistObject<GlobalOrdinal,LocalOrdinal,GlobalOrdinal,Node>& source)
02472   {
02473     // It's not clear what kind of compatibility checks on sizes can be performed here.
02474     // Epetra_CrsGraph doesn't check any sizes for compatibility.
02475     return true;
02476   }
02477 
02478 
02481   template <class LocalOrdinal, class GlobalOrdinal, class Node>
02482   void CrsGraph<LocalOrdinal,GlobalOrdinal,Node>::copyAndPermute(
02483                           const DistObject<GlobalOrdinal,LocalOrdinal,GlobalOrdinal,Node> & source,
02484                           size_t numSameIDs,
02485                           const Teuchos::ArrayView<const LocalOrdinal> &permuteToLIDs,
02486                           const Teuchos::ArrayView<const LocalOrdinal> &permuteFromLIDs)
02487   {
02488     const CrsGraph<LocalOrdinal,GlobalOrdinal,Node>& src_graph = dynamic_cast<const CrsGraph<LocalOrdinal,GlobalOrdinal,Node>&>(source);
02489     TEST_FOR_EXCEPTION(permuteToLIDs.size() != permuteFromLIDs.size(), std::runtime_error,
02490          Teuchos::typeName(*this) << "::copyAndPermute: permuteToLIDs and permuteFromLIDs must have the same size.");
02491     bool src_filled = src_graph.isFillComplete();
02492 
02493     Teuchos::Array<GlobalOrdinal> row_copy;
02494     LocalOrdinal myid = 0;
02495     for (size_t i=0; i<numSameIDs; ++i, ++myid) {
02496       GlobalOrdinal gid = src_graph.getMap()->getGlobalElement(myid);
02497       if (src_filled) {
02498         size_t row_length = src_graph.getNumEntriesInGlobalRow(gid);
02499         row_copy.resize(row_length);
02500         size_t check_row_length = 0;
02501         src_graph.getGlobalRowCopy(gid, row_copy(), check_row_length);
02502         insertGlobalIndices( gid, row_copy() );
02503       }
02504       else {
02505         Teuchos::ArrayRCP<const GlobalOrdinal> row = src_graph.getGlobalRowView(gid);
02506         insertGlobalIndices( gid, row() );
02507       }
02508     }
02509 
02510     for (LocalOrdinal i=0; i<permuteToLIDs.size(); ++i) {
02511       GlobalOrdinal mygid = this->getMap()->getGlobalElement(permuteToLIDs[i]);
02512       GlobalOrdinal srcgid= src_graph.getMap()->getGlobalElement(permuteFromLIDs[i]);
02513       if (src_filled) {
02514         size_t row_length = src_graph.getNumEntriesInGlobalRow(srcgid);
02515         row_copy.resize(row_length);
02516         size_t check_row_length = 0;
02517         src_graph.getGlobalRowCopy(srcgid, row_copy(), check_row_length);
02518         insertGlobalIndices( mygid, row_copy() );
02519       }
02520       else {
02521         Teuchos::ArrayRCP<const GlobalOrdinal> row = src_graph.getGlobalRowView(srcgid);
02522         insertGlobalIndices( mygid, row() );
02523       }
02524     }
02525   }
02526 
02527 
02530   template <class LocalOrdinal, class GlobalOrdinal, class Node>
02531   void CrsGraph<LocalOrdinal,GlobalOrdinal,Node>::packAndPrepare(
02532                           const DistObject<GlobalOrdinal,LocalOrdinal,GlobalOrdinal,Node> & source,
02533                           const Teuchos::ArrayView<const LocalOrdinal> &exportLIDs,
02534                           Teuchos::Array<GlobalOrdinal> &exports,
02535                           const Teuchos::ArrayView<size_t> & numPacketsPerLID,
02536                           size_t& constantNumPackets,
02537                           Distributor &distor)
02538   {
02539     TEST_FOR_EXCEPTION(exportLIDs.size() != numPacketsPerLID.size(), std::runtime_error,
02540                     Teuchos::typeName(*this) << "::packAndPrepare: exportLIDs and numPacketsPerLID must have the same size.");
02541     const CrsGraph<LocalOrdinal,GlobalOrdinal,Node>& src_graph = dynamic_cast<const CrsGraph<LocalOrdinal,GlobalOrdinal,Node>&>(source);
02542     // We don't check whether src_graph has had fillComplete called, because it doesn't matter whether the
02543     // *source* graph has been fillComplete'd. The target graph can not be fillComplete'd yet.
02544     TEST_FOR_EXCEPTION(this->isFillComplete() == true, std::runtime_error,
02545          Teuchos::typeName(*this) << "::copyAndPermute: import/export operations are not allowed on destination CrsGraph after fillComplete has been called.");
02546     constantNumPackets = 0;
02547     // first set the contents of numPacketsPerLID, and accumulate a total-num-packets:
02548     size_t totalNumPackets = 0;
02549     Teuchos::Array<GlobalOrdinal> row;
02550     for (LocalOrdinal i=0; i<exportLIDs.size(); ++i) {
02551       GlobalOrdinal GID = src_graph.getMap()->getGlobalElement(exportLIDs[i]);
02552       size_t row_length = src_graph.getNumEntriesInGlobalRow(GID);
02553       numPacketsPerLID[i] = row_length;
02554       totalNumPackets += row_length;
02555     }
02556 
02557     exports.resize(totalNumPackets);
02558 
02559     // now loop again and pack rows of indices into exports:
02560     size_t exportsOffset = 0;
02561     for (LocalOrdinal i=0; i<exportLIDs.size(); ++i) {
02562       GlobalOrdinal GID = src_graph.getMap()->getGlobalElement(exportLIDs[i]);
02563       size_t row_length = src_graph.getNumEntriesInGlobalRow(GID);
02564       row.resize(row_length);
02565       size_t check_row_length = 0;
02566       src_graph.getGlobalRowCopy(GID, row(), check_row_length);
02567       typename Teuchos::Array<GlobalOrdinal>::const_iterator
02568         row_iter = row.begin(), row_end = row.end();
02569       size_t j = 0;
02570       for (; row_iter != row_end; ++row_iter, ++j) {
02571         exports[exportsOffset+j] = *row_iter;
02572       }
02573       exportsOffset += row.size();
02574     }
02575   }
02576 
02577 
02580   template <class LocalOrdinal, class GlobalOrdinal, class Node>
02581   void CrsGraph<LocalOrdinal,GlobalOrdinal,Node>::unpackAndCombine(
02582                             const Teuchos::ArrayView<const LocalOrdinal> &importLIDs,
02583                             const Teuchos::ArrayView<const GlobalOrdinal> &imports,
02584                             const Teuchos::ArrayView<size_t> &numPacketsPerLID,
02585                             size_t constantNumPackets,
02586                             Distributor & /* distor */,
02587                             CombineMode /* CM */)
02588   {
02589     // We are not checking the value of the CombineMode input-argument.
02590     // For CrsGraph, we only support import/export operations if fillComplete has not yet been called.
02591     // Any incoming column-indices are inserted into the target graph. In this context, CombineMode values
02592     // of ADD vs INSERT are equivalent. What is the meaning of REPLACE for CrsGraph? If a duplicate column-index
02593     // is inserted, it will be compressed out when fillComplete is called.
02594     // 
02595     // Note: I think REPLACE means that an existing row is replaced by the imported row, i.e., the existing indices are cleared. CGB, 6/17/2010
02596 
02597     TEST_FOR_EXCEPTION(importLIDs.size() != numPacketsPerLID.size(), std::runtime_error,
02598                     Teuchos::typeName(*this) << "::unpackAndCombine: importLIDs and numPacketsPerLID must have the same size.");
02599     TEST_FOR_EXCEPTION(this->isFillComplete() == true, std::runtime_error,
02600          Teuchos::typeName(*this) << "::unpackAndCombine: import/export operations are not allowed on destination CrsGraph after fillComplete has been called.");
02601     size_t importsOffset = 0;
02602     typename Teuchos::ArrayView<const LocalOrdinal>::iterator
02603       impLIDiter = importLIDs.begin(), impLIDend = importLIDs.end();
02604     size_t i = 0;
02605     for (; impLIDiter != impLIDend; ++impLIDiter, ++i) {
02606       LocalOrdinal row_length = numPacketsPerLID[i];
02607       const Teuchos::ArrayView<const GlobalOrdinal> row(&imports[importsOffset], row_length);
02608       insertGlobalIndices(this->getMap()->getGlobalElement(*impLIDiter), row);
02609       importsOffset += row_length;
02610     }
02611   }
02612 
02613 } // namespace Tpetra
02614 
02615 
02616 //
02617 // Explicit instantiation macro
02618 //
02619 // Must be expanded from within the Tpetra namespace!
02620 //
02621 
02622 #define TPETRA_CRSGRAPH_INSTANT(LO,GO,NODE) \
02623   \
02624   template class CrsGraph< LO , GO , NODE >; \
02625 
02626 
02627 #endif // TPETRA_CRSGRAPH_DEF_HPP
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines
Generated on Wed Apr 13 10:21:40 2011 for Tpetra Matrix/Vector Services by  doxygen 1.6.3