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