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

Generated on Wed May 12 21:40:14 2010 for Tpetra Matrix/Vector Services by  doxygen 1.4.7