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

Generated on Tue Jul 13 09:39:06 2010 for Tpetra Matrix/Vector Services by  doxygen 1.4.7