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