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