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