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 #ifndef TPETRA_CRSGRAPH_DEF_HPP
00043 #define TPETRA_CRSGRAPH_DEF_HPP
00044 
00045 #include <Tpetra_Distributor.hpp>
00046 #include <Teuchos_Assert.hpp>
00047 #include <Teuchos_NullIteratorTraits.hpp>
00048 #include <Teuchos_as.hpp>
00049 #include <algorithm>
00050 #include <string>
00051 #include <utility>
00052 #include <Teuchos_SerialDenseMatrix.hpp>
00053 
00054 #ifdef DOXYGEN_USE_ONLY
00055   #include "Tpetra_CrsGraph_decl.hpp"
00056 #endif
00057 
00058 namespace Tpetra {
00059 
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   , fillComplete_(false)
00078   , indicesAreSorted_ (true)
00079   , noRedundancies_ (true)
00080   , haveLocalConstants_ (false)
00081   , haveGlobalConstants_ (false)
00082   , haveRowInfo_(true)
00083   , sortGhostsAssociatedWithEachProcessor_(true)
00084   {
00085     typedef Teuchos::OrdinalTraits<size_t> OTST;
00086     staticAssertions();
00087     TEUCHOS_TEST_FOR_EXCEPTION(maxNumEntriesPerRow == OTST::invalid(),
00088       std::invalid_argument, "The allocation hint must be a valid size_t value, "
00089       "which in this case means it must not be Teuchos::OrdinalTraits<size_t>::"
00090       "invalid().");
00091     resumeFill(params);
00092     checkInternalState();
00093   }
00094 
00097   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00098   CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::
00099   CrsGraph (const RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > &rowMap,
00100             const RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > &colMap,
00101             size_t maxNumEntriesPerRow,
00102             ProfileType pftype,
00103             const RCP<ParameterList>& params)
00104   : DistObject<GlobalOrdinal,LocalOrdinal,GlobalOrdinal,Node>(rowMap)
00105   , rowMap_(rowMap)
00106   , colMap_(colMap)
00107   , nodeNumEntries_(0)
00108   , nodeNumAllocated_(OrdinalTraits<size_t>::invalid())
00109   , pftype_(pftype)
00110   , numAllocForAllRows_(maxNumEntriesPerRow)
00111   , indicesAreAllocated_(false)
00112   , indicesAreLocal_(false)
00113   , indicesAreGlobal_(false)
00114   , fillComplete_(false)
00115   , indicesAreSorted_(true)
00116   , noRedundancies_(true)
00117   , haveLocalConstants_ (false)
00118   , haveGlobalConstants_ (false)
00119   , haveRowInfo_(true)
00120   , sortGhostsAssociatedWithEachProcessor_(true)
00121   {
00122     typedef Teuchos::OrdinalTraits<size_t> OTST;
00123     staticAssertions();
00124     TEUCHOS_TEST_FOR_EXCEPTION(maxNumEntriesPerRow == OTST::invalid(),
00125       std::invalid_argument, "The allocation hint must be a valid size_t value, "
00126       "which in this case means it must not be Teuchos::OrdinalTraits<size_t>::"
00127       "invalid().");
00128     resumeFill(params);
00129     checkInternalState();
00130   }
00131 
00134   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00135   CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::
00136   CrsGraph (const RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > &rowMap,
00137             const ArrayRCP<const size_t> &NumEntriesPerRowToAlloc,
00138             ProfileType pftype,
00139             const RCP<ParameterList>& params)
00140   : DistObject<GlobalOrdinal,LocalOrdinal,GlobalOrdinal,Node>(rowMap)
00141   , rowMap_(rowMap)
00142   , nodeNumEntries_(0)
00143   , nodeNumAllocated_(OrdinalTraits<size_t>::invalid())
00144   , pftype_(pftype)
00145   , numAllocPerRow_(NumEntriesPerRowToAlloc)
00146   , numAllocForAllRows_(0)
00147   , indicesAreAllocated_(false)
00148   , indicesAreLocal_(false)
00149   , indicesAreGlobal_(false)
00150   , fillComplete_(false)
00151   , indicesAreSorted_(true)
00152   , noRedundancies_(true)
00153   , haveLocalConstants_ (false)
00154   , haveGlobalConstants_ (false)
00155   , haveRowInfo_(true)
00156   , sortGhostsAssociatedWithEachProcessor_(true)
00157   {
00158     typedef Teuchos::OrdinalTraits<size_t> OTST;
00159     const char tfecfFuncName[] = "CrsGraph(rowMap,NumEntriesPerRowToAlloc)";
00160     staticAssertions();
00161     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC((size_t)NumEntriesPerRowToAlloc.size() != getNodeNumRows(), std::invalid_argument,
00162         ": NumEntriesPerRowToAlloc must have as many entries as specified by rowMap for this node.");
00163     for (size_t r=0; r < getNodeNumRows(); ++r) {
00164       const size_t curRowCount = NumEntriesPerRowToAlloc[r];
00165       TEUCHOS_TEST_FOR_EXCEPTION(curRowCount == OTST::invalid(),
00166         std::invalid_argument, "NumEntriesPerRowToAlloc[" << r << "] specifies "
00167         "an invalid number of entries (Teuchos::OrdinalTraits<size_t>::"
00168         "invalid()).");
00169     }
00170     resumeFill(params);
00171     checkInternalState();
00172   }
00173 
00174 
00177   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00178   CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::
00179   CrsGraph (const RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > &rowMap,
00180             const RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > &colMap,
00181             const ArrayRCP<const size_t> &NumEntriesPerRowToAlloc,
00182             ProfileType pftype,
00183             const RCP<ParameterList>& params)
00184   : DistObject<GlobalOrdinal,LocalOrdinal,GlobalOrdinal,Node>(rowMap)
00185   , rowMap_(rowMap)
00186   , colMap_(colMap)
00187   , nodeNumEntries_(0)
00188   , nodeNumAllocated_(OrdinalTraits<size_t>::invalid())
00189   , pftype_(pftype)
00190   , numAllocPerRow_(NumEntriesPerRowToAlloc)
00191   , numAllocForAllRows_(0)
00192   , indicesAreAllocated_(false)
00193   , indicesAreLocal_(false)
00194   , indicesAreGlobal_(false)
00195   , fillComplete_(false)
00196   , indicesAreSorted_(true)
00197   , noRedundancies_(true)
00198   , haveLocalConstants_ (false)
00199   , haveGlobalConstants_ (false)
00200   , haveRowInfo_(true)
00201   , sortGhostsAssociatedWithEachProcessor_(true)
00202   {
00203     typedef Teuchos::OrdinalTraits<size_t> OTST;
00204     const char tfecfFuncName[] = "CrsGraph(rowMap,colMap,NumEntriesPerRowToAlloc)";
00205     staticAssertions();
00206     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC((size_t)NumEntriesPerRowToAlloc.size() != getNodeNumRows(), std::invalid_argument,
00207         ": NumEntriesPerRowToAlloc must have as many entries as specified by rowMap for this node.");
00208     for (size_t r=0; r < getNodeNumRows(); ++r) {
00209       const size_t curRowCount = NumEntriesPerRowToAlloc[r];
00210       TEUCHOS_TEST_FOR_EXCEPTION(curRowCount == OTST::invalid(),
00211         std::invalid_argument, "NumEntriesPerRowToAlloc[" << r << "] specifies "
00212         "an invalid number of entries (Teuchos::OrdinalTraits<size_t>::"
00213         "invalid()).");
00214     }
00215     resumeFill(params);
00216     checkInternalState();
00217   }
00218 
00219 
00222   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00223   CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::
00224   CrsGraph (const RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > &rowMap,
00225             const RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > &colMap,
00226             const ArrayRCP<size_t> & rowPointers,
00227             const ArrayRCP<LocalOrdinal> & columnIndices,
00228             const RCP<ParameterList>& params)
00229   : DistObject<GlobalOrdinal,LocalOrdinal,GlobalOrdinal,Node>(rowMap)
00230   , rowMap_(rowMap)
00231   , colMap_(colMap)
00232   , nodeNumEntries_(0)
00233   , nodeNumAllocated_(OrdinalTraits<size_t>::invalid())
00234   , pftype_(StaticProfile)
00235   , numAllocForAllRows_(0)
00236   , indicesAreAllocated_(true)
00237   , indicesAreLocal_(true)
00238   , indicesAreGlobal_(false)
00239   , fillComplete_(false)
00240   , indicesAreSorted_(true)
00241   , noRedundancies_(true)
00242   , haveLocalConstants_ (false)
00243   , haveGlobalConstants_ (false)
00244   , haveRowInfo_(true)
00245   , sortGhostsAssociatedWithEachProcessor_(true)
00246   {
00247     staticAssertions();
00248     globalNumEntries_ = globalNumDiags_ = globalMaxNumRowEntries_ = OrdinalTraits<global_size_t>::invalid();
00249     setAllIndices(rowPointers,columnIndices);
00250     checkInternalState();
00251   }
00252 
00255   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00256   CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::~CrsGraph()
00257   {}
00258 
00259   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00260   RCP<const ParameterList>
00261   CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::
00262   getValidParameters () const
00263   {
00264     RCP<ParameterList> params = parameterList ("Tpetra::CrsGraph");
00265 
00266     // Make a sublist for the Import.
00267     RCP<ParameterList> importSublist = parameterList ("Import");
00268 
00269     // FIXME (mfh 02 Apr 2012) We should really have the Import and
00270     // Export objects fill in these lists.  However, we don't want to
00271     // create an Import or Export unless we need them.  For now, we
00272     // know that the Import and Export just pass the list directly to
00273     // their Distributor, so we can create a Distributor here
00274     // (Distributor's constructor is a lightweight operation) and have
00275     // it fill in the list.
00276 
00277     // Fill in Distributor default parameters by creating a
00278     // Distributor and asking it to do the work.
00279     Distributor distributor (rowMap_->getComm(), importSublist);
00280     params->set ("Import", *importSublist, "How the Import performs communication.");
00281 
00282     // Make a sublist for the Export.  For now, it's a clone of the
00283     // Import sublist.  It's not a shallow copy, though, since we
00284     // might like the Import to do communication differently than the
00285     // Export.
00286     params->set ("Export", *importSublist, "How the Export performs communication.");
00287 
00288     return params;
00289   }
00290 
00291   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00292   void
00293   CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::
00294   setParameterList (const RCP<ParameterList>& params)
00295   {
00296     RCP<const ParameterList> validParams = getValidParameters ();
00297     params->validateParametersAndSetDefaults (*validParams);
00298     this->setMyParamList (params);
00299   }
00300 
00301   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00302   global_size_t CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getGlobalNumRows() const
00303   {
00304     return rowMap_->getGlobalNumElements();
00305   }
00306 
00307   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00308   global_size_t CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getGlobalNumCols() const
00309   {
00310     const char tfecfFuncName[] = "getGlobalNumCols()";
00311     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
00312       isFillComplete() == false, std::runtime_error,
00313       ": requires domain map, which requires that fillComplete() has been "
00314       "called.");
00315     return getDomainMap()->getGlobalNumElements();
00316   }
00317 
00318   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00319   size_t CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getNodeNumRows() const
00320   {
00321     return rowMap_->getNodeNumElements();
00322   }
00323 
00324   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00325   size_t CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getNodeNumCols() const
00326   {
00327     const char tfecfFuncName[] = "getNodeNumCols()";
00328     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
00329       hasColMap() == false, std::runtime_error,
00330       ": requires column map.  This requires either that a custom column Map "
00331       "was given to the constructor, or that fillComplete() has been called.");
00332     return colMap_->getNodeNumElements();
00333   }
00334 
00335   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00336   size_t CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getNodeNumDiags() const
00337   {
00338     return nodeNumDiags_;
00339   }
00340 
00341   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00342   global_size_t CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getGlobalNumDiags() const
00343   {
00344     return globalNumDiags_;
00345   }
00346 
00347   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00348   RCP<Node>
00349   CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getNode() const
00350   {
00351     return rowMap_->getNode();
00352   }
00353 
00354   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00355   RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> >
00356   CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getRowMap() const
00357   {
00358     return rowMap_;
00359   }
00360 
00361   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00362   RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> >
00363   CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getColMap() const
00364   {
00365     return colMap_;
00366   }
00367 
00368   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00369   RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> >
00370   CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getDomainMap() const
00371   {
00372     return domainMap_;
00373   }
00374 
00375   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00376   RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> >
00377   CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getRangeMap() const
00378   {
00379     return rangeMap_;
00380   }
00381 
00382   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00383   RCP<const Import<LocalOrdinal,GlobalOrdinal,Node> >
00384   CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getImporter() const
00385   {
00386     return importer_;
00387   }
00388 
00389   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00390   RCP<const Export<LocalOrdinal,GlobalOrdinal,Node> >
00391   CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getExporter() const
00392   {
00393     return exporter_;
00394   }
00395 
00396   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00397   bool CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::hasColMap() const
00398   {
00399     return (colMap_ != null);
00400   }
00401 
00402   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00403   bool CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::isStorageOptimized() const
00404   {
00405     bool isOpt = indicesAreAllocated_ && (numRowEntries_ == null) && (getNodeNumRows() > 0);
00406 #ifdef HAVE_TPETRA_DEBUG
00407     const char tfecfFuncName[] = "isStorageOptimized()";
00408     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC( (isOpt == true) && (getProfileType() == DynamicProfile), std::logic_error,
00409         ": Violated stated post-conditions. Please contact Tpetra team.");
00410 #endif
00411     return isOpt;
00412   }
00413 
00414   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00415   ProfileType CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getProfileType() const
00416   {
00417     return pftype_;
00418   }
00419 
00420   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00421   global_size_t CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getGlobalNumEntries() const
00422   {
00423     return globalNumEntries_;
00424   }
00425 
00426   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00427   size_t CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getNodeNumEntries() const
00428   {
00429     return nodeNumEntries_;
00430   }
00431 
00432   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00433   global_size_t CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getGlobalMaxNumRowEntries() const
00434   {
00435     return globalMaxNumRowEntries_;
00436   }
00437 
00438   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00439   size_t CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getNodeMaxNumRowEntries() const
00440   {
00441     return nodeMaxNumRowEntries_;
00442   }
00443 
00444   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00445   bool CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::isFillComplete() const
00446   {
00447     return fillComplete_;
00448   }
00449 
00450   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00451   bool CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::isFillActive() const
00452   {
00453     return !fillComplete_;
00454   }
00455 
00456   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00457   bool CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::isUpperTriangular() const
00458   {
00459     return upperTriangular_;
00460   }
00461 
00462   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00463   bool CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::isLowerTriangular() const
00464   {
00465     return lowerTriangular_;
00466   }
00467 
00468   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00469   bool CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::isLocallyIndexed() const
00470   {
00471     return indicesAreLocal_;
00472   }
00473 
00474   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00475   bool CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::isGloballyIndexed() const
00476   {
00477     return indicesAreGlobal_;
00478   }
00479 
00480   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00481   size_t CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getNodeAllocationSize() const
00482   {
00483     return nodeNumAllocated_;
00484   }
00485 
00486   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00487   RCP<const Comm<int> >
00488   CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getComm() const
00489   {
00490     return rowMap_->getComm();
00491   }
00492 
00493   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00494   GlobalOrdinal CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getIndexBase() const
00495   {
00496     return rowMap_->getIndexBase();
00497   }
00498 
00499   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00500   bool CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::indicesAreAllocated() const
00501   {
00502     return indicesAreAllocated_;
00503   }
00504 
00505 
00508   //                                                                         //
00509   //                    Internal utility methods                             //
00510   //                                                                         //
00513 
00514 
00517   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00518   bool CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::isSorted() const
00519   {
00520     return indicesAreSorted_;
00521   }
00522 
00523 
00526   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00527   bool CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::isMerged() const
00528   {
00529     return noRedundancies_;
00530   }
00531 
00534   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00535   void CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::setLocallyModified ()
00536   {
00537     // FIXME (mfh 07 May 2013) How do we know that the change
00538     // introduced a redundancy, or even that it invalidated the sorted
00539     // order of indices?  CrsGraph has always made this conservative
00540     // guess.  It could be a bit costly to check at insertion time,
00541     // though.
00542     indicesAreSorted_ = false;
00543     noRedundancies_ = false;
00544 
00545     // We've modified the graph, so we'll have to recompute local
00546     // constants like the number of diagonal entries on this process.
00547     haveLocalConstants_ = false;
00548   }
00549 
00552   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00553   void
00554   CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::
00555   allocateIndices (ELocalGlobal lg)
00556   {
00557     // This is a protected function, only callable by us.  If it was
00558     // called incorrectly, it is our fault.  That's why the tests
00559     // below throw std::logic_error instead of std::invalid_argument.
00560     const char tfecfFuncName[] = "allocateIndices()";
00561     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
00562       isLocallyIndexed() && lg==GlobalIndices, std::logic_error,
00563       ": The graph is locally indexed, but Tpetra code is calling this method "
00564       "with lg=GlobalIndices.  Please report this bug to the Tpetra developers.");
00565     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
00566       isGloballyIndexed() && lg==LocalIndices, std::logic_error,
00567       ": The graph is globally indexed, but Tpetra code is calling this method "
00568       "with lg=LocalIndices.  Please report this bug to the Tpetra developers.");
00569     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
00570       indicesAreAllocated(), std::logic_error, ": The graph's indices are "
00571       "already allocated, but Tpetra code is calling allocateIndices() again.  "
00572       "Please report this bug to the Tpetra developers.");
00573 
00574     const size_t numRows = getNodeNumRows();
00575     indicesAreLocal_  = (lg == LocalIndices);
00576     indicesAreGlobal_ = (lg == GlobalIndices);
00577     nodeNumAllocated_ = 0;
00578     if (numAllocPerRow_ == null && getNodeNumRows() > 0) {
00579       // this wastes memory, temporarily, but it simplifies the code
00580       // and interfaces to follow
00581       ArrayRCP<size_t> tmpnumallocperrow = arcp<size_t>(numRows);
00582       std::fill(tmpnumallocperrow.begin(), tmpnumallocperrow.end(), numAllocForAllRows_);
00583       numAllocPerRow_ = tmpnumallocperrow;
00584     }
00585     //
00586     if (getProfileType() == StaticProfile) {
00587       //
00588       //  STATIC ALLOCATION PROFILE
00589       //
00590       // Have the local sparse kernels object allocate row offsets for
00591       // us, with first-touch allocation if applicable.  This is not
00592       // as important for global indices, because we never use global
00593       // indices in sparse kernels, but we might as well use the code
00594       // that we have for both the local and global indices cases.
00595       // Furthermore, first-touch allocation ensures that we don't
00596       // take up too much memory in any one NUMA affinity region.
00597       rowPtrs_ = LocalMatOps::allocRowPtrs( getRowMap()->getNode(), numAllocPerRow_() );
00598       if (lg == LocalIndices) {
00599         lclInds1D_ = LocalMatOps::template allocStorage<LocalOrdinal>( getRowMap()->getNode(), rowPtrs_() );
00600       }
00601       else {
00602         gblInds1D_ = LocalMatOps::template allocStorage<GlobalOrdinal>( getRowMap()->getNode(), rowPtrs_() );
00603       }
00604       nodeNumAllocated_ = rowPtrs_[numRows];
00605     }
00606     else {
00607       //
00608       //  DYNAMIC ALLOCATION PROFILE
00609       //
00610       typename ArrayRCP<const size_t>::iterator numalloc = numAllocPerRow_.begin();
00611       size_t howmany = numAllocForAllRows_;
00612       if (lg == LocalIndices) {
00613         lclInds2D_ = arcp< Array<LocalOrdinal> >(numRows);
00614         nodeNumAllocated_ = 0;
00615         for (size_t i=0; i < numRows; ++i) {
00616           if (numAllocPerRow_ != null) howmany = *numalloc++;
00617           nodeNumAllocated_ += howmany;
00618           if (howmany > 0) lclInds2D_[i].resize(howmany);
00619         }
00620       }
00621       else { // allocate global indices
00622         gblInds2D_ = arcp< Array<GlobalOrdinal> >(numRows);
00623         nodeNumAllocated_ = 0;
00624         for (size_t i=0; i < numRows; ++i) {
00625           if (numAllocPerRow_ != null) howmany = *numalloc++;
00626           nodeNumAllocated_ += howmany;
00627           if (howmany > 0) gblInds2D_[i].resize(howmany);
00628         }
00629       }
00630     }
00631     if (numRows > 0) {
00632       numRowEntries_ = arcp<size_t>(numRows);
00633       std::fill( numRowEntries_.begin(), numRowEntries_.end(), 0 );
00634     }
00635     // done with these
00636     numAllocForAllRows_ = 0;
00637     numAllocPerRow_     = null;
00638     indicesAreAllocated_ = true;
00639     checkInternalState();
00640   }
00641 
00642 
00645   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00646   template <class T>
00647   ArrayRCP<T> CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::allocateValues1D() const
00648   {
00649     const char tfecfFuncName[] = "allocateValues()";
00650     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
00651         indicesAreAllocated() == false,
00652         std::runtime_error, ": graph indices must be allocated before values."
00653     );
00654     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
00655         getProfileType() != StaticProfile,
00656         std::runtime_error, ": graph indices must be allocated in a static profile."
00657     );
00658     ArrayRCP<T> values1D;
00659     values1D = LocalMatOps::template allocStorage<T>( getRowMap()->getNode(), rowPtrs_() );
00660     return values1D;
00661   }
00662 
00663 
00666   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00667   template <class T>
00668   ArrayRCP<Array<T> >
00669   CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::allocateValues2D() const
00670   {
00671     const char tfecfFuncName[] = "allocateValues2D()";
00672     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
00673         indicesAreAllocated() == false,
00674         std::runtime_error, ": graph indices must be allocated before values."
00675     );
00676     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
00677         getProfileType() != DynamicProfile,
00678         std::runtime_error, ": graph indices must be allocated in a dynamic profile."
00679     );
00680     ArrayRCP<Array<T> > values2D;
00681     values2D = arcp<Array<T> > (getNodeNumRows ());
00682     if (lclInds2D_ != null) {
00683       const size_t numRows = lclInds2D_.size ();
00684       for (size_t r = 0; r < numRows; ++r) {
00685         values2D[r].resize (lclInds2D_[r].size ());
00686       }
00687     }
00688     else if (gblInds2D_ != null) {
00689       const size_t numRows = gblInds2D_.size ();
00690       for (size_t r = 0; r < numRows; ++r) {
00691         values2D[r].resize (gblInds2D_[r].size ());
00692       }
00693     }
00694     return values2D;
00695   }
00696 
00697 
00698 //   /////////////////////////////////////////////////////////////////////////////
00699 //   /////////////////////////////////////////////////////////////////////////////
00700 //   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00701 //   template <ELocalGlobal lg, class T>
00702 //   RowInfo
00703 //   CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::
00704 //   updateAllocAndValues (RowInfo rowinfo,
00705 //                         size_t newAllocSize,
00706 //                         Array<T>& rowVals)
00707 //   {
00708 // #ifdef HAVE_TPETRA_DEBUG
00709 //     TEUCHOS_TEST_FOR_EXCEPT( ! rowMap_->isNodeLocalElement(rowinfo.localRow) );
00710 //     TEUCHOS_TEST_FOR_EXCEPT( newAllocSize < rowinfo.allocSize );
00711 //     TEUCHOS_TEST_FOR_EXCEPT( (lg == LocalIndices && ! isLocallyIndexed()) ||
00712 //                              (lg == GlobalIndices && ! isGloballyIndexed()) );
00713 //     TEUCHOS_TEST_FOR_EXCEPT( newAllocSize == 0 );
00714 //     TEUCHOS_TEST_FOR_EXCEPT( ! indicesAreAllocated() );
00715 // #endif
00716 //     // ArrayRCP::resize automatically copies over values on reallocation.
00717 //     if (lg == LocalIndices) {
00718 //       lclInds2D_[rowinfo.localRow].resize (newAllocSize);
00719 //     }
00720 //     else { // lg == GlobalIndices
00721 //       gblInds2D_[rowinfo.localRow].resize (newAllocSize);
00722 //     }
00723 //     rowVals.resize (newAllocSize);
00724 //     nodeNumAllocated_ += (newAllocSize - rowinfo.allocSize);
00725 //     rowinfo.allocSize = newAllocSize;
00726 //     return rowinfo;
00727 //   }
00728 
00729 
00732   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00733   ArrayView<const LocalOrdinal>
00734   CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::
00735   getLocalView (RowInfo rowinfo) const
00736   {
00737     ArrayView<const LocalOrdinal> view;
00738     if (rowinfo.allocSize > 0) {
00739       if (lclInds1D_ != null) {
00740         view = lclInds1D_ (rowinfo.offset1D, rowinfo.allocSize);
00741       }
00742       else if (! lclInds2D_[rowinfo.localRow].empty ()) {
00743         view = lclInds2D_[rowinfo.localRow] ();
00744       }
00745     }
00746     return view;
00747   }
00748 
00749 
00752   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00753   ArrayView<LocalOrdinal>
00754   CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::
00755   getLocalViewNonConst (RowInfo rowinfo)
00756   {
00757     ArrayView<LocalOrdinal> view;
00758     if (rowinfo.allocSize > 0) {
00759       if (lclInds1D_ != null) {
00760         view = lclInds1D_ (rowinfo.offset1D, rowinfo.allocSize);
00761       }
00762       else if (! lclInds2D_[rowinfo.localRow].empty()) {
00763         view = lclInds2D_[rowinfo.localRow] ();
00764       }
00765     }
00766     return view;
00767   }
00768 
00769 
00772   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00773   ArrayView<const GlobalOrdinal>
00774   CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::
00775   getGlobalView (RowInfo rowinfo) const
00776   {
00777     ArrayView<const GlobalOrdinal> view;
00778     if (rowinfo.allocSize > 0) {
00779       if (gblInds1D_ != null) {
00780         view = gblInds1D_ (rowinfo.offset1D, rowinfo.allocSize);
00781       }
00782       else if (! gblInds2D_[rowinfo.localRow].empty()) {
00783         view = gblInds2D_[rowinfo.localRow] ();
00784       }
00785     }
00786     return view;
00787   }
00788 
00789 
00792   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00793   ArrayView<GlobalOrdinal>
00794   CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::
00795   getGlobalViewNonConst (RowInfo rowinfo)
00796   {
00797     ArrayView<GlobalOrdinal> view;
00798     if (rowinfo.allocSize > 0) {
00799       if (gblInds1D_ != null) {
00800         view = gblInds1D_ (rowinfo.offset1D, rowinfo.allocSize);
00801       }
00802       else if (!gblInds2D_[rowinfo.localRow].empty()) {
00803         view = gblInds2D_[rowinfo.localRow] ();
00804       }
00805     }
00806     return view;
00807   }
00808 
00809 
00812   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00813   RowInfo CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getRowInfo(size_t myRow) const
00814   {
00815 #ifdef HAVE_TPETRA_DEBUG
00816     const char tfecfFuncName[] = "getRowInfo";
00817     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
00818         rowMap_->isNodeLocalElement (myRow) == false,
00819         std::logic_error,
00820         ": The given (local) row index myRow = " << myRow
00821         << " does not belong to the graph's row Map.  "
00822         "This probably indicates a bug in Tpetra::CrsGraph or Tpetra::CrsMatrix.  "
00823         "Please report this to the Tpetra developers.");
00824     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
00825       hasRowInfo() == false, std::logic_error,
00826       ": Late catch! Graph does not have row info anymore.  "
00827       "Error should have been caught earlier.  Please contact Tpetra team.");
00828 #endif // HAVE_TPETRA_DEBUG
00829     const size_t STINV = Teuchos::OrdinalTraits<size_t>::invalid ();
00830     RowInfo ret;
00831     ret.localRow = myRow;
00832     if (nodeNumAllocated_ != 0 && nodeNumAllocated_ != STINV) {
00833       // graph data structures have the info that we need
00834       //
00835       // if static graph, offsets tell us the allocation size
00836       if (getProfileType() == StaticProfile) {
00837         ret.offset1D  = rowPtrs_[myRow];
00838         ret.allocSize = rowPtrs_[myRow+1] - rowPtrs_[myRow];
00839         if (numRowEntries_ == null) {
00840           ret.numEntries = ret.allocSize;
00841         }
00842         else {
00843           ret.numEntries = numRowEntries_[myRow];
00844         }
00845       }
00846       else {
00847         ret.offset1D = STINV;
00848         if (isLocallyIndexed ()) {
00849           ret.allocSize = lclInds2D_[myRow].size();
00850         }
00851         else {
00852           ret.allocSize = gblInds2D_[myRow].size();
00853         }
00854         ret.numEntries = numRowEntries_[myRow];
00855       }
00856     }
00857     else if (nodeNumAllocated_ == 0) {
00858       // have performed allocation, but the graph has no allocation or entries
00859       ret.allocSize = 0;
00860       ret.numEntries = 0;
00861       ret.offset1D = STINV;
00862     }
00863     else if (indicesAreAllocated () == false) {
00864       // haven't performed allocation yet; probably won't hit this code
00865       if (numAllocPerRow_ == null) {
00866         ret.allocSize = numAllocForAllRows_;
00867       }
00868       else {
00869         ret.allocSize = numAllocPerRow_[myRow];
00870       }
00871       ret.numEntries = 0;
00872       ret.offset1D = STINV;
00873     }
00874     else {
00875       // don't know how we ended up here...
00876       TEUCHOS_TEST_FOR_EXCEPT(true);
00877     }
00878     return ret;
00879   }
00880 
00881 
00884   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00885   void CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::staticAssertions() const
00886   {
00887     // Assumption: sizeof(GlobalOrdinal) >= sizeof(LocalOrdinal):
00888     //     This is so that we can store local indices in the memory
00889     //     formerly occupied by global indices.
00890     //
00891     // Assumption: max(GlobalOrdinal) >= max(LocalOrdinal) and
00892     //   max(size_t) >= max(LocalOrdinal)
00893     //     This is so that we can represent any LocalOrdinal as a
00894     //     size_t, and any LocalOrdinal as a GlobalOrdinal
00895     Teuchos::CompileTimeAssert<sizeof(GlobalOrdinal) < sizeof(LocalOrdinal)> cta_size1; (void)cta_size1;
00896     Teuchos::CompileTimeAssert<sizeof(global_size_t) < sizeof(size_t)      > cta_size2; (void)cta_size2;
00897     // can't call max() with CompileTimeAssert, because it isn't a constant expression; will need to make this a runtime check
00898     std::string msg = typeName(*this) + ": Object cannot be allocated with stated template arguments: size assumptions are not valid.";
00899     TEUCHOS_TEST_FOR_EXCEPTION( (size_t)OrdinalTraits<LocalOrdinal>::max() > OrdinalTraits<size_t>::max(),          std::runtime_error, msg);
00900     TEUCHOS_TEST_FOR_EXCEPTION( (global_size_t)OrdinalTraits<LocalOrdinal>::max() > (global_size_t)OrdinalTraits<GlobalOrdinal>::max(),           std::runtime_error, msg);
00901     TEUCHOS_TEST_FOR_EXCEPTION( (size_t)OrdinalTraits<GlobalOrdinal>::max() > OrdinalTraits<global_size_t>::max(),  std::runtime_error, msg);
00902     TEUCHOS_TEST_FOR_EXCEPTION( OrdinalTraits<size_t>::max() > OrdinalTraits<global_size_t>::max(),                 std::runtime_error, msg);
00903   }
00904 
00905 
00908   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00909   template <ELocalGlobal lg>
00910   size_t CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::filterIndices(const SLocalGlobalNCViews &inds) const
00911   {
00912     const map_type& cmap = *colMap_;
00913     Teuchos::CompileTimeAssert<lg != GlobalIndices && lg != LocalIndices> cta_lg;
00914     (void)cta_lg;
00915     size_t numFiltered = 0;
00916 #ifdef HAVE_TPETRA_DEBUG
00917     size_t numFiltered_debug = 0;
00918 #endif
00919     if (lg == GlobalIndices) {
00920       ArrayView<GlobalOrdinal> ginds = inds.ginds;
00921       typename ArrayView<GlobalOrdinal>::iterator fend = ginds.begin();
00922       typename ArrayView<GlobalOrdinal>::iterator cptr = ginds.begin();
00923       while (cptr != ginds.end()) {
00924         if (cmap.isNodeGlobalElement(*cptr)) {
00925           *fend++ = *cptr;
00926 #ifdef HAVE_TPETRA_DEBUG
00927           ++numFiltered_debug;
00928 #endif
00929         }
00930         ++cptr;
00931       }
00932       numFiltered = fend - ginds.begin();
00933     }
00934     else if (lg == LocalIndices) {
00935       ArrayView<LocalOrdinal> linds = inds.linds;
00936       typename ArrayView<LocalOrdinal>::iterator fend = linds.begin();
00937       typename ArrayView<LocalOrdinal>::iterator cptr = linds.begin();
00938       while (cptr != linds.end()) {
00939         if (cmap.isNodeLocalElement(*cptr)) {
00940           *fend++ = *cptr;
00941 #ifdef HAVE_TPETRA_DEBUG
00942           ++numFiltered_debug;
00943 #endif
00944         }
00945         ++cptr;
00946       }
00947       numFiltered = fend - linds.begin();
00948     }
00949 #ifdef HAVE_TPETRA_DEBUG
00950     TEUCHOS_TEST_FOR_EXCEPT( numFiltered != numFiltered_debug );
00951 #endif
00952     return numFiltered;
00953   }
00954 
00955 
00956 //   /////////////////////////////////////////////////////////////////////////////
00957 //   /////////////////////////////////////////////////////////////////////////////
00958 //   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00959 //   template <class T>
00960 //   size_t CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::
00961 //   filterGlobalIndicesAndValues (const ArrayView<GlobalOrdinal>& ginds,
00962 //                                 const ArrayView<T>& vals) const
00963 //   {
00964 //     const Map<LocalOrdinal,GlobalOrdinal,Node>& cmap = *colMap_;
00965 //     size_t numFiltered = 0;
00966 //     typename ArrayView<T>::iterator fvalsend = vals.begin();
00967 //     typename ArrayView<T>::iterator valscptr = vals.begin();
00968 // #ifdef HAVE_TPETRA_DEBUG
00969 //     size_t numFiltered_debug = 0;
00970 // #endif
00971 //     typename ArrayView<GlobalOrdinal>::iterator fend = ginds.begin();
00972 //     typename ArrayView<GlobalOrdinal>::iterator cptr = ginds.begin();
00973 //     while (cptr != ginds.end()) {
00974 //       if (cmap.isNodeGlobalElement (*cptr)) {
00975 //         *fend++ = *cptr;
00976 //         *fvalsend++ = *valscptr;
00977 // #ifdef HAVE_TPETRA_DEBUG
00978 //         ++numFiltered_debug;
00979 // #endif
00980 //       }
00981 //       ++cptr;
00982 //       ++valscptr;
00983 //     }
00984 //     numFiltered = fend - ginds.begin();
00985 // #ifdef HAVE_TPETRA_DEBUG
00986 //     TEUCHOS_TEST_FOR_EXCEPT( numFiltered != numFiltered_debug );
00987 //     TEUCHOS_TEST_FOR_EXCEPT( valscptr != vals.end() );
00988 //     const size_t numFilteredActual =
00989 //       Teuchos::as<size_t> (fvalsend - vals.begin ());
00990 //     TEUCHOS_TEST_FOR_EXCEPT( numFiltered != numFilteredActual );
00991 // #endif
00992 //     return numFiltered;
00993 //   }
00994 
00995 
00996 //   /////////////////////////////////////////////////////////////////////////////
00997 //   /////////////////////////////////////////////////////////////////////////////
00998 //   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00999 //   template <class T>
01000 //   size_t
01001 //   CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::
01002 //   filterLocalIndicesAndValues (const ArrayView<LocalOrdinal>& linds,
01003 //                                const ArrayView<T>& vals) const
01004 //   {
01005 //     const Map<LocalOrdinal,GlobalOrdinal,Node>& cmap = *colMap_;
01006 //     size_t numFiltered = 0;
01007 //     typename ArrayView<T>::iterator fvalsend = vals.begin();
01008 //     typename ArrayView<T>::iterator valscptr = vals.begin();
01009 // #ifdef HAVE_TPETRA_DEBUG
01010 //     size_t numFiltered_debug = 0;
01011 // #endif
01012 //     typename ArrayView<LocalOrdinal>::iterator fend = linds.begin();
01013 //     typename ArrayView<LocalOrdinal>::iterator cptr = linds.begin();
01014 //     while (cptr != linds.end()) {
01015 //       if (cmap.isNodeLocalElement (*cptr)) {
01016 //         *fend++ = *cptr;
01017 //         *fvalsend++ = *valscptr;
01018 // #ifdef HAVE_TPETRA_DEBUG
01019 //         ++numFiltered_debug;
01020 // #endif
01021 //       }
01022 //       ++cptr;
01023 //       ++valscptr;
01024 //     }
01025 //     numFiltered = fend - linds.begin();
01026 // #ifdef HAVE_TPETRA_DEBUG
01027 //     TEUCHOS_TEST_FOR_EXCEPT( numFiltered != numFiltered_debug );
01028 //     TEUCHOS_TEST_FOR_EXCEPT( valscptr != vals.end() );
01029 //     const size_t numFilteredActual =
01030 //       Teuchos::as<size_t> (fvalsend - vals.begin ());
01031 //     TEUCHOS_TEST_FOR_EXCEPT( numFiltered != numFilteredActual );
01032 // #endif
01033 //     return numFiltered;
01034 //   }
01035 
01036 
01039   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
01040   size_t
01041   CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::
01042   insertIndices (const RowInfo& rowinfo,
01043                  const SLocalGlobalViews &newInds,
01044                  const ELocalGlobal lg,
01045                  const ELocalGlobal I)
01046   {
01047 #ifdef HAVE_TPETRA_DEBUG
01048     TEUCHOS_TEST_FOR_EXCEPTION(
01049       lg != GlobalIndices && lg != LocalIndices, std::invalid_argument,
01050       "Tpetra::CrsGraph::insertIndices: lg must be either GlobalIndices or "
01051       "LocalIndices.");
01052 #endif // HAVE_TPETRA_DEBUG
01053     size_t numNewInds = 0;
01054     if (lg == GlobalIndices) { // input indices are global
01055       ArrayView<const GlobalOrdinal> new_ginds = newInds.ginds;
01056       numNewInds = new_ginds.size();
01057       if (I == GlobalIndices) { // store global indices
01058         ArrayView<GlobalOrdinal> gind_view = getGlobalViewNonConst(rowinfo);
01059         std::copy(new_ginds.begin(), new_ginds.end(), gind_view.begin()+rowinfo.numEntries);
01060       }
01061       else if (I == LocalIndices) { // store local indices
01062         ArrayView<LocalOrdinal> lind_view = getLocalViewNonConst(rowinfo);
01063         typename ArrayView<const GlobalOrdinal>::iterator         in = new_ginds.begin();
01064         const typename ArrayView<const GlobalOrdinal>::iterator stop = new_ginds.end();
01065         typename ArrayView<LocalOrdinal>::iterator out = lind_view.begin()+rowinfo.numEntries;
01066         while (in != stop) {
01067           *out++ = colMap_->getLocalElement (*in++);
01068         }
01069       }
01070     }
01071     else if (lg == LocalIndices) { // input indices are local
01072       ArrayView<const LocalOrdinal> new_linds = newInds.linds;
01073       numNewInds = new_linds.size();
01074       if (I == LocalIndices) { // store local indices
01075         ArrayView<LocalOrdinal> lind_view = getLocalViewNonConst(rowinfo);
01076         std::copy(new_linds.begin(), new_linds.end(), lind_view.begin()+rowinfo.numEntries);
01077       }
01078       else if (I == GlobalIndices) {
01079         TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error, "Tpetra::CrsGraph::"
01080           "insertIndices: the case where the input indices are local and the "
01081           "indices to write are global (lg=LocalIndices, I=GlobalIndices) is "
01082           "not implemented, because it does not make sense." << std::endl <<
01083           "If you have correct local column indices, that means the graph has "
01084           "a column Map.  In that case, you should be storing local indices.");
01085       }
01086     }
01087     numRowEntries_[rowinfo.localRow] += numNewInds;
01088     nodeNumEntries_ += numNewInds;
01089     setLocallyModified ();
01090     return numNewInds;
01091   }
01092 
01095   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
01096   void CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::
01097   insertGlobalIndicesImpl (const LocalOrdinal myRow,
01098                            const ArrayView<const GlobalOrdinal> &indices)
01099   {
01100     const char* tfecfFuncName ("insertGlobalIndicesImpl");
01101 
01102     RowInfo rowInfo = getRowInfo(myRow);
01103     const size_t numNewInds = indices.size();
01104     const size_t newNumEntries = rowInfo.numEntries + numNewInds;
01105     if (newNumEntries > rowInfo.allocSize) {
01106       TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
01107         getProfileType() == StaticProfile, std::runtime_error,
01108         ": new indices exceed statically allocated graph structure.");
01109 
01110       // update allocation, doubling size to reduce number of reallocations
01111       size_t newAllocSize = 2*rowInfo.allocSize;
01112       if (newAllocSize < newNumEntries)
01113         newAllocSize = newNumEntries;
01114       gblInds2D_[myRow].resize(newAllocSize);
01115       nodeNumAllocated_ += (newAllocSize - rowInfo.allocSize);
01116     }
01117 
01118     // Copy new indices at end of global index array
01119     if (gblInds1D_ != null)
01120       std::copy(indices.begin(), indices.end(),
01121                 gblInds1D_.begin()+rowInfo.offset1D+rowInfo.numEntries);
01122     else
01123       std::copy(indices.begin(), indices.end(),
01124                 gblInds2D_[myRow].begin()+rowInfo.numEntries);
01125     numRowEntries_[myRow] += numNewInds;
01126     nodeNumEntries_ += numNewInds;
01127     setLocallyModified ();
01128 
01129 #ifdef HAVE_TPETRA_DEBUG
01130     {
01131       const size_t chkNewNumEntries = getNumEntriesInLocalRow (myRow);
01132       TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
01133         chkNewNumEntries != newNumEntries, std::logic_error,
01134         ": Internal logic error. Please contact Tpetra team.");
01135     }
01136 #endif
01137   }
01138 
01141   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
01142   void CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::
01143   insertLocalIndicesImpl (const LocalOrdinal myRow,
01144                           const ArrayView<const LocalOrdinal> &indices)
01145   {
01146     const char* tfecfFuncName ("insertLocallIndicesImpl");
01147 
01148     RowInfo rowInfo = getRowInfo(myRow);
01149     const size_t numNewInds = indices.size();
01150     const size_t newNumEntries = rowInfo.numEntries + numNewInds;
01151     if (newNumEntries > rowInfo.allocSize) {
01152       TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
01153         getProfileType() == StaticProfile, std::runtime_error,
01154         ": new indices exceed statically allocated graph structure.");
01155 
01156       // update allocation, doubling size to reduce number of reallocations
01157       size_t newAllocSize = 2*rowInfo.allocSize;
01158       if (newAllocSize < newNumEntries)
01159         newAllocSize = newNumEntries;
01160       lclInds2D_[myRow].resize(newAllocSize);
01161       nodeNumAllocated_ += (newAllocSize - rowInfo.allocSize);
01162     }
01163 
01164     // Insert new indices at end of lclInds array
01165     if (lclInds1D_ != null)
01166       std::copy(indices.begin(), indices.end(),
01167                 lclInds1D_.begin()+rowInfo.offset1D+rowInfo.numEntries);
01168     else
01169       std::copy(indices.begin(), indices.end(),
01170                 lclInds2D_[myRow].begin()+rowInfo.numEntries);
01171     numRowEntries_[myRow] += numNewInds;
01172     nodeNumEntries_ += numNewInds;
01173     setLocallyModified ();
01174 #ifdef HAVE_TPETRA_DEBUG
01175     {
01176       const size_t chkNewNumEntries = getNumEntriesInLocalRow (myRow);
01177       TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
01178         chkNewNumEntries != newNumEntries, std::logic_error,
01179         ": Internal logic error. Please contact Tpetra team.");
01180     }
01181 #endif
01182   }
01183 
01184   // /////////////////////////////////////////////////////////////////////////////
01185   // /////////////////////////////////////////////////////////////////////////////
01186   // template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
01187   // template <class Scalar>
01188   // void CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::
01189   // insertIndicesAndValues (const RowInfo& rowInfo,
01190   //                         const SLocalGlobalViews& newInds,
01191   //                         const ArrayView<Scalar>& oldRowVals,
01192   //                         const ArrayView<const Scalar>& newRowVals,
01193   //                         const ELocalGlobal lg,
01194   //                         const ELocalGlobal I)
01195   // {
01196   //   const size_t numNewInds = insertIndices (rowInfo, newInds, lg, I);
01197   //   typename ArrayView<const Scalar>::const_iterator newRowValsBegin =
01198   //     newRowVals.begin ();
01199   //   std::copy (newRowValsBegin, newRowValsBegin + numNewInds,
01200   //              oldRowVals.begin () + rowInfo.numEntries);
01201   // }
01202 
01203   // /////////////////////////////////////////////////////////////////////////////
01204   // /////////////////////////////////////////////////////////////////////////////
01205   // template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
01206   // template <class Scalar, class BinaryFunction>
01207   // void
01208   // CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::
01209   // transformLocalValues (RowInfo rowInfo,
01210   //                       const Teuchos::ArrayView<Scalar>& rowVals,
01211   //                       const Teuchos::ArrayView<const LocalOrdinal>& inds,
01212   //                       const Teuchos::ArrayView<const Scalar>& newVals,
01213   //                       BinaryFunction f) const
01214   // {
01215   //   const size_t STINV = Teuchos::OrdinalTraits<size_t>::invalid();
01216   //   const size_t numElts = Teuchos::as<size_t> (inds.size ());
01217   //   size_t hint = 0; // Guess for the current index k into rowVals
01218 
01219   //   // Get a view of the column indices in the row.  This amortizes
01220   //   // the cost of getting the view over all the entries of inds.
01221   //   ArrayView<const LocalOrdinal> colInds = getLocalView (rowInfo);
01222 
01223   //   for (size_t j = 0; j < numElts; ++j) {
01224   //     const size_t k = findLocalIndex (rowInfo, inds[j], colInds, hint);
01225   //     if (k != STINV) {
01226   //       rowVals[k] = f( rowVals[k], newVals[j] );
01227   //       hint = k+1;
01228   //     }
01229   //   }
01230   // }
01231 
01232 
01233   // /////////////////////////////////////////////////////////////////////////////
01234   // /////////////////////////////////////////////////////////////////////////////
01235   // template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
01236   // template <class Scalar, class BinaryFunction>
01237   // void
01238   // CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::
01239   // transformGlobalValues (RowInfo rowInfo,
01240   //                        const Teuchos::ArrayView<Scalar>& rowVals,
01241   //                        const Teuchos::ArrayView<const GlobalOrdinal>& inds,
01242   //                        const Teuchos::ArrayView<const Scalar>& newVals,
01243   //                        BinaryFunction f) const
01244   // {
01245   //   const size_t STINV = Teuchos::OrdinalTraits<size_t>::invalid();
01246   //   const size_t numElts = Teuchos::as<size_t> (inds.size ());
01247   //   size_t hint = 0; // hint is a guess as to wheter the index is
01248 
01249   //   for (size_t j = 0; j < numElts; ++j) {
01250   //     const size_t k = findGlobalIndex (rowInfo, inds[j], hint);
01251   //     if (k != STINV) {
01252   //       rowVals[k] = f( rowVals[k], newVals[j] );
01253   //       hint = k+1;
01254   //     }
01255   //   }
01256   // }
01257 
01258 
01259 
01262   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
01263   void CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::sortRowIndices(RowInfo rowinfo)
01264   {
01265     if (rowinfo.numEntries > 0) {
01266       ArrayView<LocalOrdinal> inds_view = getLocalViewNonConst(rowinfo);
01267       std::sort(inds_view.begin(), inds_view.begin() + rowinfo.numEntries);
01268     }
01269   }
01270 
01271 
01274   // in the future, perhaps this could use std::sort with a boost::zip_iterator
01275   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
01276   template <class Scalar>
01277   void CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::sortRowIndicesAndValues(RowInfo rowinfo, ArrayView<Scalar> values)
01278   {
01279     if (rowinfo.numEntries > 0) {
01280       ArrayView<LocalOrdinal> inds_view = getLocalViewNonConst(rowinfo);
01281       sort2(inds_view.begin(), inds_view.begin()+rowinfo.numEntries, values.begin());
01282     }
01283   }
01284 
01285 
01288   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
01289   void CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::mergeRowIndices(RowInfo rowinfo)
01290   {
01291     const char tfecfFuncName[] = "mergRowIndices()";
01292     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
01293       isStorageOptimized() == true, std::logic_error,
01294       ": The graph is already storage optimized, so we shouldn't be merging any indices."
01295       " Please report this bug to the Tpetra developers.");
01296     ArrayView<LocalOrdinal> inds_view = getLocalViewNonConst(rowinfo);
01297     typename ArrayView<LocalOrdinal>::iterator beg, end, newend;
01298     beg = inds_view.begin();
01299     end = inds_view.begin() + rowinfo.numEntries;
01300     newend = std::unique(beg,end);
01301     const size_t mergedEntries = newend - beg;
01302 #ifdef HAVE_TPETRA_DEBUG
01303     // merge should not have eliminated any entries; if so, the assignment below will destory the packed structure
01304     TEUCHOS_TEST_FOR_EXCEPT( isStorageOptimized() && mergedEntries != rowinfo.numEntries );
01305 #endif
01306     numRowEntries_[rowinfo.localRow] = mergedEntries;
01307     nodeNumEntries_ -= (rowinfo.numEntries - mergedEntries);
01308   }
01309 
01310 
01313   // in the future, this could use std::unique with a boost::zip_iterator
01314   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
01315   template<class Scalar>
01316   void CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::
01317   mergeRowIndicesAndValues (RowInfo rowinfo, const ArrayView<Scalar>& rowValues)
01318   {
01319     const char tfecfFuncName[] = "mergeRowIndicesAndValues";
01320     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
01321       isStorageOptimized(), std::logic_error, ": It is invalid to call this "
01322       "method if the graph's storage has already been optimized." << std::endl
01323       << "Please report this bug to the Tpetra developers.");
01324 
01325     typedef typename ArrayView<Scalar>::iterator Iter;
01326     Iter rowValueIter = rowValues.begin ();
01327     ArrayView<LocalOrdinal> inds_view = getLocalViewNonConst (rowinfo);
01328     typename ArrayView<LocalOrdinal>::iterator beg, end, newend;
01329 
01330     // beg,end define a half-exclusive interval over which to iterate.
01331     beg = inds_view.begin();
01332     end = inds_view.begin() + rowinfo.numEntries;
01333     newend = beg;
01334     if (beg != end) {
01335       typename ArrayView<LocalOrdinal>::iterator cur = beg + 1;
01336       Iter vcur = rowValueIter + 1;
01337       Iter vend = rowValueIter;
01338       cur = beg+1;
01339       while (cur != end) {
01340         if (*cur != *newend) {
01341           // new entry; save it
01342           ++newend;
01343           ++vend;
01344           (*newend) = (*cur);
01345           (*vend) = (*vcur);
01346         }
01347         else {
01348           // old entry; merge it
01349           //(*vend) = f (*vend, *vcur);
01350           (*vend) += *vcur;
01351         }
01352         ++cur;
01353         ++vcur;
01354       }
01355       ++newend; // one past the last entry, per typical [beg,end) semantics
01356     }
01357     const size_t mergedEntries = newend - beg;
01358 #ifdef HAVE_TPETRA_DEBUG
01359     // merge should not have eliminated any entries; if so, the
01360     // assignment below will destroy the packed structure
01361     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
01362       isStorageOptimized() && mergedEntries != rowinfo.numEntries,
01363       std::logic_error,
01364       ": Merge was incorrect; it eliminated entries from the graph.  "
01365       << std::endl << "Please report this bug to the Tpetra developers.");
01366 #endif // HAVE_TPETRA_DEBUG
01367     numRowEntries_[rowinfo.localRow] = mergedEntries;
01368     nodeNumEntries_ -= (rowinfo.numEntries - mergedEntries);
01369   }
01370 
01371 
01374   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
01375   void CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::
01376   setDomainRangeMaps (const RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> >& domainMap,
01377                       const RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> >& rangeMap)
01378   {
01379     // simple pointer comparison for equality
01380     if (domainMap_ != domainMap) {
01381       domainMap_ = domainMap;
01382       importer_ = null;
01383     }
01384     if (rangeMap_ != rangeMap) {
01385       rangeMap_  = rangeMap;
01386       exporter_ = null;
01387     }
01388   }
01389 
01390 
01393   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
01394   size_t
01395   CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::
01396   findLocalIndex (RowInfo rowinfo, LocalOrdinal ind, size_t hint) const
01397   {
01398     ArrayView<const LocalOrdinal> colInds = getLocalView (rowinfo);
01399     return this->findLocalIndex (rowinfo, ind, colInds, hint);
01400   }
01401 
01402 
01405   template <class LocalOrdinal,
01406             class GlobalOrdinal,
01407             class Node,
01408             class LocalMatOps>
01409   size_t
01410   CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::
01411   findLocalIndex (RowInfo rowinfo,
01412                   LocalOrdinal ind,
01413                   ArrayView<const LocalOrdinal> colInds,
01414                   size_t hint) const
01415   {
01416     typedef typename ArrayView<const LocalOrdinal>::iterator IT;
01417 
01418     // If the hint was correct, then the hint is the offset to return.
01419     if (hint < rowinfo.numEntries && colInds[hint] == ind) {
01420       return hint;
01421     }
01422 
01423     // The hint was wrong, so we must search for the given column
01424     // index in the column indices for the given row.  How we do the
01425     // search depends on whether the graph's column indices are
01426     // sorted.
01427     IT beg = colInds.begin ();
01428     IT end = beg + rowinfo.numEntries;
01429     IT ptr = beg + rowinfo.numEntries; // "null"
01430     bool found = true;
01431 
01432     if (isSorted ()) {
01433       // binary search
01434       std::pair<IT,IT> p = std::equal_range (beg, end, ind);
01435       if (p.first == p.second) {
01436         found = false;
01437       }
01438       else {
01439         ptr = p.first;
01440       }
01441     }
01442     else {
01443       // direct search
01444       ptr = std::find (beg, end, ind);
01445       if (ptr == end) {
01446         found = false;
01447       }
01448     }
01449 
01450     if (found) {
01451       return Teuchos::as<size_t> (ptr - beg);
01452     }
01453     else {
01454       return Teuchos::OrdinalTraits<size_t>::invalid ();
01455     }
01456   }
01457 
01458 
01461   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
01462   size_t
01463   CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::
01464   findGlobalIndex (RowInfo rowinfo, GlobalOrdinal ind, size_t hint) const
01465   {
01466     using Teuchos::ArrayView;
01467     typedef typename ArrayView<const GlobalOrdinal>::iterator IT;
01468 
01469     // Don't let an invalid global column index through.
01470     if (ind == Teuchos::OrdinalTraits<GlobalOrdinal>::invalid ()) {
01471       return Teuchos::OrdinalTraits<size_t>::invalid ();
01472     }
01473 
01474     ArrayView<const GlobalOrdinal> indices = getGlobalView (rowinfo);
01475 
01476     // We don't actually require that the hint be a valid index.
01477     // If it is not in range, we just ignore it.
01478     if (hint < rowinfo.numEntries && indices[hint] == ind) {
01479       return hint;
01480     }
01481 
01482     IT beg = indices.begin ();
01483     IT end = indices.begin () + rowinfo.numEntries; // not indices.end()
01484     if (isSorted ()) { // use binary search
01485       const std::pair<IT,IT> p = std::equal_range (beg, end, ind);
01486       if (p.first == p.second) { // range of matching entries is empty
01487         return Teuchos::OrdinalTraits<size_t>::invalid ();
01488       } else {
01489         return p.first - beg;
01490       }
01491     }
01492     else { // not sorted; must use linear search
01493       const IT loc = std::find (beg, end, ind);
01494       if (loc == end) {
01495         return Teuchos::OrdinalTraits<size_t>::invalid ();
01496       } else {
01497         return loc - beg;
01498       }
01499     }
01500   }
01501 
01502 
01505   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
01506   void CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::clearGlobalConstants()
01507   {
01508     globalNumEntries_       = OrdinalTraits<global_size_t>::invalid();
01509     globalNumDiags_         = OrdinalTraits<global_size_t>::invalid();
01510     globalMaxNumRowEntries_ = OrdinalTraits<global_size_t>::invalid();
01511     haveGlobalConstants_    = false;
01512   }
01513 
01514 
01517   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
01518   void CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::
01519   checkInternalState () const
01520   {
01521 #ifdef HAVE_TPETRA_DEBUG
01522     const global_size_t GSTI = OrdinalTraits<global_size_t>::invalid();
01523     const size_t         STI = OrdinalTraits<size_t>::invalid();
01524     std::string err = typeName(*this) + "::checkInternalState(): Likely internal logic error. Please contact Tpetra team.";
01525     // check the internal state of this data structure
01526     // this is called by numerous state-changing methods, in a debug build, to ensure that the object
01527     // always remains in a valid state
01528     // the graph should have been allocated with a row map
01529     TEUCHOS_TEST_FOR_EXCEPTION( rowMap_ == null,                     std::logic_error, err );
01530     // am either complete or active
01531     TEUCHOS_TEST_FOR_EXCEPTION( isFillActive() == isFillComplete(),  std::logic_error, err );
01532     // if active, i have no local graph
01533     TEUCHOS_TEST_FOR_EXCEPTION( isFillActive() && lclGraph_ != null, std::logic_error, err );
01534     // if the graph has been fill completed, then all maps should be present
01535     TEUCHOS_TEST_FOR_EXCEPTION( isFillComplete() == true && (colMap_ == null || rangeMap_ == null || domainMap_ == null), std::logic_error, err );
01536     // if storage has been optimized, then indices should have been allocated (even if trivially so)
01537     TEUCHOS_TEST_FOR_EXCEPTION( isStorageOptimized() == true && indicesAreAllocated() == false, std::logic_error, err );
01538     // if storage has been optimized, then number of allocated is now the number of entries
01539     TEUCHOS_TEST_FOR_EXCEPTION( isStorageOptimized() == true && getNodeAllocationSize() != getNodeNumEntries(), std::logic_error, err );
01540     // if graph doesn't have the global constants, then they should all be marked as invalid
01541     TEUCHOS_TEST_FOR_EXCEPTION( haveGlobalConstants_ == false && ( globalNumEntries_ != GSTI || globalNumDiags_ != GSTI || globalMaxNumRowEntries_ != GSTI ), std::logic_error, err );
01542     // if the graph has global cosntants, then they should be valid.
01543     TEUCHOS_TEST_FOR_EXCEPTION( haveGlobalConstants_ == true && ( globalNumEntries_ == GSTI || globalNumDiags_ == GSTI || globalMaxNumRowEntries_ == GSTI ), std::logic_error, err );
01544     TEUCHOS_TEST_FOR_EXCEPTION( haveGlobalConstants_ == true && ( globalNumEntries_ < nodeNumEntries_ || globalNumDiags_ < nodeNumDiags_ || globalMaxNumRowEntries_ < nodeMaxNumRowEntries_ ),
01545                         std::logic_error, err );
01546     // if indices are allocated, then the allocation specifications should have been released
01547     TEUCHOS_TEST_FOR_EXCEPTION( indicesAreAllocated() == true  && (numAllocForAllRows_ != 0 || numAllocPerRow_ != null),                        std::logic_error, err );
01548     // if indices are not allocated, then information dictating allocation quantities should be present
01549     TEUCHOS_TEST_FOR_EXCEPTION( indicesAreAllocated() == false && (nodeNumAllocated_ != STI || nodeNumEntries_ != 0),                           std::logic_error, err );
01550     // if storage is optimized, then profile should be static
01551     TEUCHOS_TEST_FOR_EXCEPTION( isStorageOptimized() && pftype_ != StaticProfile,                                                               std::logic_error, err );
01552     // if rowPtrs_ exists, it is required to have N+1 rows, and rowPtrs_[N] == gblInds1D_.size()/lclInds1D_.size()
01553     TEUCHOS_TEST_FOR_EXCEPTION( isGloballyIndexed() && rowPtrs_ != null && ((size_t)rowPtrs_.size() != getNodeNumRows()+1 || rowPtrs_[getNodeNumRows()] != (size_t)gblInds1D_.size()), std::logic_error, err );
01554     TEUCHOS_TEST_FOR_EXCEPTION(  isLocallyIndexed() && rowPtrs_ != null && ((size_t)rowPtrs_.size() != getNodeNumRows()+1 || rowPtrs_[getNodeNumRows()] != (size_t)lclInds1D_.size()), std::logic_error, err );
01555     // if profile is dynamic and we have allocated, then 2D allocations should be present
01556     TEUCHOS_TEST_FOR_EXCEPTION( pftype_ == DynamicProfile && indicesAreAllocated() && getNodeNumRows() > 0 && lclInds2D_ == null && gblInds2D_ == null,
01557                                                                                                                                         std::logic_error, err );
01558     // if profile is dynamic, then numentries and 2D indices are needed and should be present
01559     TEUCHOS_TEST_FOR_EXCEPTION( pftype_ == DynamicProfile && indicesAreAllocated() && getNodeNumRows() > 0 && (numRowEntries_ == null || (lclInds2D_ == null && gblInds2D_ == null)),
01560                                                                                                                                         std::logic_error, err );
01561     // if profile is dynamic, then 1D allocations should not be present
01562     TEUCHOS_TEST_FOR_EXCEPTION( pftype_ == DynamicProfile && (lclInds1D_ != null || gblInds1D_ != null),                                std::logic_error, err );
01563     // if profile is dynamic, then row offsets should not be present
01564     TEUCHOS_TEST_FOR_EXCEPTION( pftype_ == DynamicProfile && rowPtrs_ != null,                                                          std::logic_error, err );
01565     // if profile is static and we have allocated non-trivially, then 1D allocations should be present
01566     TEUCHOS_TEST_FOR_EXCEPTION( pftype_ == StaticProfile && indicesAreAllocated() && getNodeAllocationSize() > 0 && lclInds1D_ == null && gblInds1D_ == null,
01567                                                                                                                                         std::logic_error, err );
01568     // if profile is static, then 2D allocations should not be present
01569     TEUCHOS_TEST_FOR_EXCEPTION( pftype_ == StaticProfile && (lclInds2D_ != null || gblInds2D_ != null),                                 std::logic_error, err );
01570     // if indices are not allocated, then none of the buffers should be.
01571     TEUCHOS_TEST_FOR_EXCEPTION( indicesAreAllocated() == false && (rowPtrs_ != null || numRowEntries_ != null ||
01572                                                                  lclInds1D_ != null || lclInds2D_ != null ||
01573                                                                  gblInds1D_ != null || gblInds2D_ != null),                             std::logic_error, err );
01574     // indices may be local or global only if they are allocated (numAllocated is redundant; could simply be indicesAreLocal_ || indicesAreGlobal_)
01575     TEUCHOS_TEST_FOR_EXCEPTION( (indicesAreLocal_ == true || indicesAreGlobal_ == true) && indicesAreAllocated_ == false,               std::logic_error, err );
01576     // indices may be local or global, but not both
01577     TEUCHOS_TEST_FOR_EXCEPTION( indicesAreLocal_ == true && indicesAreGlobal_ == true,                                                  std::logic_error, err );
01578     // if indices are local, then global allocations should not be present
01579     TEUCHOS_TEST_FOR_EXCEPTION( indicesAreLocal_ == true && (gblInds1D_ != null || gblInds2D_ != null),                                 std::logic_error, err );
01580     // if indices are global, then local allocations should not be present
01581     TEUCHOS_TEST_FOR_EXCEPTION( indicesAreGlobal_ == true && (lclInds1D_ != null || lclInds2D_ != null),                                std::logic_error, err );
01582     // if indices are local, then local allocations should be present
01583     TEUCHOS_TEST_FOR_EXCEPTION( indicesAreLocal_ == true && getNodeAllocationSize() > 0 && lclInds1D_ == null && getNodeNumRows() > 0 && lclInds2D_ == null,
01584                                                                                                                           std::logic_error, err );
01585     // if indices are global, then global allocations should be present
01586     TEUCHOS_TEST_FOR_EXCEPTION( indicesAreGlobal_ == true && getNodeAllocationSize() > 0 && gblInds1D_ == null && getNodeNumRows() > 0 && gblInds2D_ == null,
01587                                                                                                                           std::logic_error, err );
01588     // if indices are allocated, then we should have recorded how many were allocated
01589     TEUCHOS_TEST_FOR_EXCEPTION( indicesAreAllocated() == true  && nodeNumAllocated_ == STI,                                       std::logic_error, err );
01590     // if indices are not allocated, then the allocation size should be marked invalid
01591     TEUCHOS_TEST_FOR_EXCEPTION( indicesAreAllocated() == false && nodeNumAllocated_ != STI,                                       std::logic_error, err );
01592     // check the actual allocations
01593     if (indicesAreAllocated()) {
01594       size_t actualNumAllocated = 0;
01595       if (pftype_ == DynamicProfile) {
01596         if (isGloballyIndexed() && gblInds2D_ != null) {
01597           for (size_t r = 0; r < getNodeNumRows(); ++r) {
01598             actualNumAllocated += gblInds2D_[r].size();
01599           }
01600         }
01601         else if (isLocallyIndexed() && lclInds2D_ != null) {
01602           for (size_t r = 0; r < getNodeNumRows(); ++r) {
01603             actualNumAllocated += lclInds2D_[r].size();
01604           }
01605         }
01606         TEUCHOS_TEST_FOR_EXCEPTION(actualNumAllocated != nodeNumAllocated_, std::logic_error, err );
01607       }
01608       else if (rowPtrs_ != null) { // pftype_ == StaticProfile
01609         actualNumAllocated = rowPtrs_[getNodeNumRows()];
01610         TEUCHOS_TEST_FOR_EXCEPTION(  isLocallyIndexed() == true && (size_t)lclInds1D_.size() != actualNumAllocated, std::logic_error, err );
01611         TEUCHOS_TEST_FOR_EXCEPTION( isGloballyIndexed() == true && (size_t)gblInds1D_.size() != actualNumAllocated, std::logic_error, err );
01612         TEUCHOS_TEST_FOR_EXCEPTION(actualNumAllocated != nodeNumAllocated_, std::logic_error, err );
01613       }
01614     }
01615 #endif
01616   }
01617 
01618 
01621   //                                                                         //
01622   //                  User-visible class methods                             //
01623   //                                                                         //
01626 
01627 
01630   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
01631   size_t CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getNumEntriesInGlobalRow(GlobalOrdinal globalRow) const
01632   {
01633     const LocalOrdinal lrow = rowMap_->getLocalElement(globalRow);
01634     size_t ret = OrdinalTraits<size_t>::invalid();
01635     if (hasRowInfo() && lrow != OrdinalTraits<LocalOrdinal>::invalid())
01636     {
01637       RowInfo rowinfo = getRowInfo(lrow);
01638       ret = rowinfo.numEntries;
01639     }
01640     return ret;
01641   }
01642 
01643 
01646   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
01647   size_t CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getNumEntriesInLocalRow(LocalOrdinal localRow) const
01648   {
01649     size_t ret = OrdinalTraits<size_t>::invalid();
01650     if (hasRowInfo() && rowMap_->isNodeLocalElement(localRow)) {
01651       RowInfo rowinfo = getRowInfo(localRow);
01652       ret = rowinfo.numEntries;
01653     }
01654     return ret;
01655   }
01656 
01657 
01660   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
01661   size_t CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getNumAllocatedEntriesInGlobalRow(GlobalOrdinal globalRow) const
01662   {
01663     const LocalOrdinal lrow = rowMap_->getLocalElement(globalRow);
01664     size_t ret = OrdinalTraits<size_t>::invalid();
01665     if (hasRowInfo() && lrow != OrdinalTraits<LocalOrdinal>::invalid())
01666     {
01667       RowInfo rowinfo = getRowInfo(lrow);
01668       ret = rowinfo.allocSize;
01669     }
01670     return ret;
01671   }
01672 
01673 
01676   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
01677   size_t CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getNumAllocatedEntriesInLocalRow(LocalOrdinal localRow) const
01678   {
01679     size_t ret = OrdinalTraits<size_t>::invalid();
01680     if (hasRowInfo() && rowMap_->isNodeLocalElement(localRow)) {
01681       RowInfo rowinfo = getRowInfo(localRow);
01682       ret = rowinfo.allocSize;
01683     }
01684     return ret;
01685   }
01686 
01687 
01690   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
01691   ArrayRCP<const size_t> CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getNodeRowPtrs() const
01692   {
01693     return rowPtrs_;
01694   }
01695 
01696 
01699   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
01700   ArrayRCP<const LocalOrdinal> CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getNodePackedIndices() const
01701   {
01702     return lclInds1D_;
01703   }
01704 
01705 
01708   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
01709   void
01710   CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::
01711   getLocalRowCopy (LocalOrdinal localRow,
01712                    const ArrayView<LocalOrdinal>& indices,
01713                    size_t& numEntries) const
01714   {
01715     using Teuchos::ArrayView;
01716     typedef LocalOrdinal LO;
01717     typedef GlobalOrdinal GO;
01718 
01719     TEUCHOS_TEST_FOR_EXCEPTION(
01720       isGloballyIndexed () && ! hasColMap (), std::runtime_error,
01721       "Tpetra::CrsGraph::getLocalRowCopy: The graph is globally indexed and "
01722       "does not have a column Map yet.  That means we don't have local indices "
01723       "for columns yet, so it doesn't make sense to call this method.  If the "
01724       "graph doesn't have a column Map yet, you should call fillComplete on "
01725       "it first.");
01726     TEUCHOS_TEST_FOR_EXCEPTION(
01727       ! hasRowInfo(), std::runtime_error,
01728       "Tpetra::CrsGraph::getLocalRowCopy: graph row information was deleted "
01729       "at fillComplete().");
01730 
01731 
01732     if (! getRowMap ()->isNodeLocalElement (localRow)) {
01733       numEntries = 0;
01734       return;
01735     }
01736 
01737     const RowInfo rowinfo = getRowInfo (localRow);
01738     const size_t theNumEntries = rowinfo.numEntries;
01739 
01740     TEUCHOS_TEST_FOR_EXCEPTION(
01741       static_cast<size_t> (indices.size ()) < theNumEntries,
01742       std::runtime_error,
01743       "Tpetra::CrsGraph::getLocalRowCopy: The given row " << localRow << " has "
01744       << theNumEntries << " entries, but indices.size() = " << indices.size ()
01745       << ", which does not suffice to store the row's indices.");
01746 
01747     numEntries = theNumEntries;
01748 
01749     if (isLocallyIndexed ()) {
01750       ArrayView<const LO> lview = getLocalView (rowinfo);
01751       std::copy (lview.begin (), lview.begin () + numEntries, indices.begin ());
01752     }
01753     else if (isGloballyIndexed ()) {
01754       ArrayView<const GO> gview = getGlobalView (rowinfo);
01755       const map_type& colMap = * (this->getColMap ());
01756       for (size_t j = 0; j < numEntries; ++j) {
01757         indices[j] = colMap.getLocalElement (gview[j]);
01758       }
01759     }
01760     else {
01761       // If the graph on the calling process is neither locally nor
01762       // globally indexed, that means it owns no column indices.
01763       //
01764       // FIXME (mfh 21 Oct 2013) It's not entirely clear to me whether
01765       // we can reach this branch, given the checks above.  However,
01766       // if that is the case, it should still be correct to call this
01767       // function if the calling process owns no column indices.
01768       numEntries = 0;
01769     }
01770   }
01771 
01772 
01775   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
01776   void
01777   CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::
01778   getGlobalRowCopy (GlobalOrdinal globalRow,
01779                     const ArrayView<GlobalOrdinal>& indices,
01780                     size_t& NumIndices) const
01781   {
01782     // One of the following is true:
01783     //
01784     // 1. The calling process owns no column indices.
01785     // 2. The graph stores global column indices.
01786     // 3. The graph has a column Map that we can use to convert any
01787     //    stored local indices into global indices for output.
01788 
01789     const char tfecfFuncName[] = "getGlobalRowCopy";
01790     const LocalOrdinal localRow = rowMap_->getLocalElement (globalRow);
01791     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
01792       localRow == Teuchos::OrdinalTraits<LocalOrdinal>::invalid(),
01793       std::runtime_error, ": globalRow (== " << globalRow << ") is not owned "
01794       "by (i.e., is not in the row Map on) the calling process with rank "
01795       << getComm ()->getRank () << ".");
01796     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
01797       ! hasRowInfo (), std::runtime_error,
01798       ": graph row information was deleted at fillComplete().");
01799 
01800     const RowInfo rowinfo = getRowInfo (static_cast<size_t> (localRow));
01801     NumIndices = rowinfo.numEntries;
01802     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
01803       static_cast<size_t> (indices.size ()) < NumIndices, std::runtime_error,
01804       ": specified storage (size == " << indices.size () << ") does not suffice "
01805       "to hold all entries in this row (NumIndices == " << NumIndices << ").");
01806 
01807     if (isLocallyIndexed ()) {
01808       ArrayView<const LocalOrdinal> lview = getLocalView (rowinfo);
01809       for (size_t j = 0; j < NumIndices; ++j) {
01810         indices[j] = colMap_->getGlobalElement (lview[j]);
01811       }
01812     }
01813     else if (isGloballyIndexed ()) {
01814       ArrayView<const GlobalOrdinal> gview = getGlobalView (rowinfo);
01815       std::copy (gview.begin (), gview.begin () + NumIndices, indices.begin ());
01816     }
01817   }
01818 
01819 
01822   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
01823   void
01824   CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::
01825   getLocalRowView (LocalOrdinal localRow, ArrayView<const LocalOrdinal> &indices) const
01826   {
01827     const char tfecfFuncName[] = "getLocalRowView";
01828     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
01829       isGloballyIndexed(), std::runtime_error,
01830       ": The graph is globally indexed, so we cannot return a view with local "
01831       "column indices.  Use getLocalRowCopy() instead.");
01832     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
01833       ! hasRowInfo (), std::runtime_error,
01834       ": graph row information was deleted at fillComplete().");
01835 
01836     indices = null;
01837     if (rowMap_->isNodeLocalElement (localRow)) {
01838       const RowInfo rowinfo = getRowInfo(localRow);
01839       if (rowinfo.numEntries > 0) {
01840         indices = getLocalView (rowinfo);
01841         indices = indices (0, rowinfo.numEntries);
01842       }
01843 #ifdef HAVE_TPETRA_DEBUG
01844       const size_t numEntriesInView = static_cast<size_t> (indices.size ());
01845       TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
01846         numEntriesInView != getNumEntriesInLocalRow (localRow),
01847         std::logic_error, ": The number of entries " << numEntriesInView
01848         << " in the view to return must be the same as getNumEntriesInLocalRow"
01849         "(localRow=" << localRow << ") = " << getNumEntriesInLocalRow (localRow)
01850         << ", but this is not the case.  "
01851         "Please report this bug to the Tpetra developers.");
01852 #endif // HAVE_TPETRA_DEBUG
01853     }
01854   }
01855 
01856 
01859   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
01860   void
01861   CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::
01862   getGlobalRowView (GlobalOrdinal globalRow,
01863                     ArrayView<const GlobalOrdinal>& indices) const
01864   {
01865     using Teuchos::as;
01866     const char tfecfFuncName[] = "getGlobalRowView";
01867     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
01868       isLocallyIndexed (), std::runtime_error,
01869       ": The graph is locally indexed, so we cannot return a view with global "
01870       "column indices.  Use getGlobalRowCopy() instead.");
01871     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
01872       ! hasRowInfo (), std::runtime_error,
01873       ": graph row information was deleted at fillComplete().");
01874 
01875     // We don't have to call isNodeGlobalElement() to test whether
01876     // globalRow belongs to the calling process.  That method requires
01877     // a global to local lookup anyway, and getLocalElement() returns
01878     // invalid() anyway if the global index wasn't found.
01879     const LocalOrdinal localRow = rowMap_->getLocalElement (globalRow);
01880     indices = null;
01881     if (localRow != Teuchos::OrdinalTraits<LocalOrdinal>::invalid ()) {
01882       const RowInfo rowInfo = getRowInfo (as<size_t> (localRow));
01883       if (rowInfo.numEntries > 0) {
01884         indices = (getGlobalView (rowInfo)) (0, rowInfo.numEntries);
01885       }
01886     }
01887 #ifdef HAVE_TPETRA_DEBUG
01888     using std::endl;
01889     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
01890       as<size_t> (indices.size ()) != getNumEntriesInGlobalRow (globalRow),
01891       std::logic_error,
01892       ": Violated stated post-conditions:"
01893       << "  indices.size () = " << indices.size () << endl
01894       << "  as<size_t> (indices.size ()) = " << as<size_t> (indices.size ())
01895       << endl << "  getNumEntriesInGlobalRow (globalRow = " << globalRow
01896       << ") = " << getNumEntriesInGlobalRow (globalRow) << endl
01897       << "Please report this bug to the Tpetra developers.");
01898 #endif // HAVE_TPETRA_DEBUG
01899   }
01900 
01901 
01904   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
01905   void
01906   CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::
01907   insertLocalIndices (const LocalOrdinal localRow,
01908                       const ArrayView<const LocalOrdinal> &indices)
01909   {
01910     using Teuchos::ArrayView;
01911     const char tfecfFuncName[] = "insertLocalIndices";
01912 
01913     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
01914       isFillActive() == false, std::runtime_error,
01915       ": requires that fill is active.");
01916     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
01917       isGloballyIndexed() == true, std::runtime_error,
01918       ": graph indices are global; use insertGlobalIndices().");
01919     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
01920       hasColMap() == false, std::runtime_error,
01921       ": cannot insert local indices without a column map.");
01922     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
01923       rowMap_->isNodeLocalElement(localRow) == false, std::runtime_error,
01924       ": row does not belong to this node.");
01925     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
01926       hasRowInfo() == false, std::runtime_error,
01927       ": graph row information was deleted at fillComplete().");
01928     if (! indicesAreAllocated ()) {
01929       allocateIndices (LocalIndices);
01930     }
01931 
01932 #ifdef HAVE_TPETRA_DEBUG
01933     // In a debug build, if the graph has a column Map, test whether
01934     // any of the given column indices are not in the column Map.
01935     // Keep track of the invalid column indices so we can tell the
01936     // user about them.
01937     if (hasColMap ()) {
01938       using Teuchos::Array;
01939       using Teuchos::toString;
01940       using std::endl;
01941       typedef LocalOrdinal LO;
01942       typedef typename ArrayView<const LO>::size_type size_type;
01943 
01944       const map_type& colMap = * (getColMap ());
01945       Array<LO> badColInds;
01946       bool allInColMap = true;
01947       for (size_type k = 0; k < indices.size (); ++k) {
01948         if (! colMap.isNodeLocalElement (indices[k])) {
01949           allInColMap = false;
01950           badColInds.push_back (indices[k]);
01951         }
01952       }
01953       if (! allInColMap) {
01954         std::ostringstream os;
01955         os << "Tpetra::CrsMatrix::insertLocalIndices: You attempted to insert "
01956           "entries in owned row " << localRow << ", at the following column "
01957           "indices: " << toString (indices) << "." << endl;
01958         os << "Of those, the following indices are not in the column Map on "
01959           "this process: " << toString (badColInds) << "." << endl << "Since "
01960           "the graph has a column Map already, it is invalid to insert entries "
01961           "at those locations.";
01962         TEUCHOS_TEST_FOR_EXCEPTION(! allInColMap, std::invalid_argument, os.str ());
01963       }
01964     }
01965 #endif // HAVE_TPETRA_DEBUG
01966 
01967     insertLocalIndicesImpl (localRow, indices);
01968 
01969 #ifdef HAVE_TPETRA_DEBUG
01970     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
01971       indicesAreAllocated() == false || isLocallyIndexed() == false,
01972       std::logic_error,
01973       ": Violated stated post-conditions. Please contact Tpetra team.");
01974 #endif
01975   }
01976 
01977 
01980   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
01981   void
01982   CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::
01983   insertLocalIndicesFiltered (const LocalOrdinal localRow,
01984                               const ArrayView<const LocalOrdinal> &indices)
01985   {
01986     typedef LocalOrdinal LO;
01987     const char tfecfFuncName[] = "insertLocalIndicesFiltered";
01988 
01989     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
01990       isFillActive() == false, std::runtime_error,
01991       ": requires that fill is active.");
01992     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
01993       isGloballyIndexed() == true, std::runtime_error,
01994       ": graph indices are global; use insertGlobalIndices().");
01995     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
01996       hasColMap() == false, std::runtime_error,
01997       ": cannot insert local indices without a column map.");
01998     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
01999       rowMap_->isNodeLocalElement(localRow) == false, std::runtime_error,
02000       ": row does not belong to this node.");
02001     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
02002       hasRowInfo() == false, std::runtime_error,
02003       ": graph row information was deleted at fillComplete().");
02004     if (indicesAreAllocated() == false) {
02005       allocateIndices(LocalIndices);
02006     }
02007 
02008      // If we have a column map, use it to filter the entries.
02009     if (hasColMap ()) {
02010       Array<LO> filtered_indices(indices);
02011       SLocalGlobalViews inds_view;
02012       SLocalGlobalNCViews inds_ncview;
02013       inds_ncview.linds = filtered_indices();
02014       const size_t numFilteredEntries =
02015         filterIndices<LocalIndices>(inds_ncview);
02016       inds_view.linds = filtered_indices (0, numFilteredEntries);
02017       insertLocalIndicesImpl(localRow, inds_view.linds);
02018     }
02019     else {
02020       insertLocalIndicesImpl(localRow, indices);
02021     }
02022 #ifdef HAVE_TPETRA_DEBUG
02023     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
02024       indicesAreAllocated() == false || isLocallyIndexed() == false,
02025       std::logic_error,
02026       ": Violated stated post-conditions. Please contact Tpetra team.");
02027 #endif
02028   }
02029 
02030 
02033   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
02034   void
02035   CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::
02036   insertGlobalIndices (const GlobalOrdinal grow,
02037                        const ArrayView<const GlobalOrdinal> &indices)
02038   {
02039     typedef LocalOrdinal LO;
02040     typedef GlobalOrdinal GO;
02041     typedef typename ArrayView<const GO>::size_type size_type;
02042     const char tfecfFuncName[] = "insertGlobalIndices";
02043 
02044     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
02045       isLocallyIndexed() == true, std::runtime_error,
02046       ": graph indices are local; use insertLocalIndices().");
02047     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
02048       hasRowInfo() == false, std::runtime_error,
02049       ": graph row information was deleted at fillComplete().");
02050     // This can't really be satisfied for now, because if we are
02051     // fillComplete(), then we are local.  In the future, this may
02052     // change.  However, the rule that modification require active
02053     // fill will not change.
02054     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
02055       isFillActive() == false, std::runtime_error,
02056       ": You are not allowed to call this method if fill is not active.  "
02057       "If fillComplete has been called, you must first call resumeFill "
02058       "before you may insert indices.");
02059     if (! indicesAreAllocated ()) {
02060       allocateIndices (GlobalIndices);
02061     }
02062     const LO myRow = rowMap_->getLocalElement (grow);
02063     if (myRow != Teuchos::OrdinalTraits<LO>::invalid ()) {
02064 #ifdef HAVE_TPETRA_DEBUG
02065       if (hasColMap ()) {
02066         using std::endl;
02067         const map_type& colMap = * (getColMap ());
02068         // In a debug build, keep track of the nonowned ("bad") column
02069         // indices, so that we can display them in the exception
02070         // message.  In a release build, just ditch the loop early if
02071         // we encounter a nonowned column index.
02072         Array<GO> badColInds;
02073         bool allInColMap = true;
02074         for (size_type k = 0; k < indices.size (); ++k) {
02075           if (! colMap.isNodeGlobalElement (indices[k])) {
02076             allInColMap = false;
02077             badColInds.push_back (indices[k]);
02078           }
02079         }
02080         if (! allInColMap) {
02081           std::ostringstream os;
02082           os << "Tpetra::CrsGraph::insertGlobalIndices: You attempted to insert "
02083             "entries in owned row " << grow << ", at the following column "
02084             "indices: " << toString (indices) << "." << endl;
02085           os << "Of those, the following indices are not in the column Map on "
02086             "this process: " << toString (badColInds) << "." << endl << "Since "
02087             "the matrix has a column Map already, it is invalid to insert "
02088             "entries at those locations.";
02089           TEUCHOS_TEST_FOR_EXCEPTION(! allInColMap, std::invalid_argument, os.str ());
02090         }
02091       }
02092 #endif // HAVE_TPETRA_DEBUG
02093       insertGlobalIndicesImpl (myRow, indices);
02094     }
02095     else { // a nonlocal row
02096       const size_type numIndices = indices.size ();
02097       // This creates the Array if it doesn't exist yet.  std::map's
02098       // operator[] does a lookup each time, so it's better to pull
02099       // nonlocals_[grow] out of the loop.
02100       std::vector<GO>& nonlocalRow = nonlocals_[grow];
02101       for (size_type k = 0; k < numIndices; ++k) {
02102         nonlocalRow.push_back (indices[k]);
02103       }
02104     }
02105 #ifdef HAVE_TPETRA_DEBUG
02106     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
02107       indicesAreAllocated() == false || isGloballyIndexed() == false,
02108       std::logic_error,
02109       ": Violated stated post-conditions. Please contact Tpetra team.");
02110 #endif
02111   }
02112 
02113 
02116   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
02117   void
02118   CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::
02119   insertGlobalIndicesFiltered (const GlobalOrdinal grow,
02120                                const ArrayView<const GlobalOrdinal> &indices)
02121   {
02122     typedef LocalOrdinal LO;
02123     typedef GlobalOrdinal GO;
02124     const char tfecfFuncName[] = "insertGlobalIndicesFiltered";
02125 
02126     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
02127       isLocallyIndexed() == true, std::runtime_error,
02128       ": graph indices are local; use insertLocalIndices().");
02129     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
02130       hasRowInfo() == false, std::runtime_error,
02131       ": graph row information was deleted at fillComplete().");
02132     // This can't really be satisfied for now, because if we are
02133     // fillComplete(), then we are local.  In the future, this may
02134     // change.  However, the rule that modification require active
02135     // fill will not change.
02136     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
02137       isFillActive() == false, std::runtime_error,
02138       ": You are not allowed to call this method if fill is not active.  "
02139       "If fillComplete has been called, you must first call resumeFill "
02140       "before you may insert indices.");
02141     if (indicesAreAllocated() == false) {
02142       allocateIndices (GlobalIndices);
02143     }
02144     const LO myRow = rowMap_->getLocalElement (grow);
02145     if (myRow != Teuchos::OrdinalTraits<LO>::invalid ()) {
02146       // If we have a column map, use it to filter the entries.
02147       if (hasColMap ()) {
02148         Array<GO> filtered_indices(indices);
02149         SLocalGlobalViews inds_view;
02150         SLocalGlobalNCViews inds_ncview;
02151         inds_ncview.ginds = filtered_indices();
02152         const size_t numFilteredEntries =
02153           filterIndices<GlobalIndices> (inds_ncview);
02154         inds_view.ginds = filtered_indices (0, numFilteredEntries);
02155         insertGlobalIndicesImpl(myRow, inds_view.ginds);
02156       }
02157       else {
02158        insertGlobalIndicesImpl(myRow, indices);
02159       }
02160     }
02161     else { // a nonlocal row
02162       typedef typename ArrayView<const GO>::size_type size_type;
02163       const size_type numIndices = indices.size ();
02164       for (size_type k = 0; k < numIndices; ++k) {
02165         nonlocals_[grow].push_back (indices[k]);
02166       }
02167     }
02168 #ifdef HAVE_TPETRA_DEBUG
02169     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
02170       indicesAreAllocated() == false || isGloballyIndexed() == false,
02171       std::logic_error,
02172       ": Violated stated post-conditions. Please contact Tpetra team.");
02173 #endif
02174   }
02175 
02178   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
02179   void CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::removeLocalIndices(LocalOrdinal lrow)
02180   {
02181     const char tfecfFuncName[] = "removeLocalIndices()";
02182     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC( isFillActive() == false,                    std::runtime_error, ": requires that fill is active.");
02183     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC( isStorageOptimized() == true,               std::runtime_error, ": cannot remove indices after optimizeStorage() has been called.");
02184     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC( isGloballyIndexed() == true,                std::runtime_error, ": graph indices are global.");
02185     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC( rowMap_->isNodeLocalElement(lrow) == false, std::runtime_error, ": row does not belong to this node.");
02186     if (indicesAreAllocated() == false) {
02187       allocateIndices(LocalIndices);
02188     }
02189     //
02190     clearGlobalConstants();
02191     //
02192     if (numRowEntries_ != null) {
02193       nodeNumEntries_ -= numRowEntries_[lrow];
02194       numRowEntries_[lrow] = 0;
02195     }
02196 #ifdef HAVE_TPETRA_DEBUG
02197     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(getNumEntriesInLocalRow(lrow) != 0 || indicesAreAllocated() == false || isLocallyIndexed() == false, std::logic_error,
02198         ": Violated stated post-conditions. Please contact Tpetra team.");
02199 #endif
02200   }
02201 
02204   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
02205   void CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::setAllIndices(const ArrayRCP<size_t> & rowPointers,const ArrayRCP<LocalOrdinal> & columnIndices)
02206   {
02207     const char tfecfFuncName[] = "setAllIndices()";
02208     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC( hasColMap() == false, std::runtime_error, ": requires a ColMap.");
02209     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC( (size_t)rowPointers.size() != getNodeNumRows()+1, std::runtime_error, ": requires rowPointers.size() == getNodeNumRows()+1");
02210     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC( lclInds1D_ != Teuchos::null || gblInds1D_ != Teuchos::null, std::runtime_error, ": cannot have 1D data structures allocated.");
02211     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC( lclInds2D_ != Teuchos::null || gblInds2D_ != Teuchos::null, std::runtime_error, ": cannot have 2D data structures allocated.");
02212 
02213     indicesAreAllocated_ = true;
02214     indicesAreLocal_     = true;
02215     pftype_              = StaticProfile; // if the profile wasn't static before, it sure is now.
02216     lclInds1D_           = columnIndices;
02217     rowPtrs_             = rowPointers;
02218     nodeNumEntries_ = nodeNumAllocated_ = rowPtrs_[getNodeNumRows()];
02219     numAllocForAllRows_  = 0;
02220     numAllocPerRow_      = null;
02221     checkInternalState();
02222   }
02223 
02226   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
02227   void
02228   CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::
02229   getNumEntriesPerLocalRowUpperBound (Teuchos::ArrayRCP<const size_t>& boundPerLocalRow,
02230                                       size_t& boundForAllLocalRows,
02231                                       bool& boundSameForAllLocalRows) const
02232   {
02233     // The three output arguments.  We assign them to the actual
02234     // output arguments at the end, in order to implement
02235     // transactional semantics.
02236     Teuchos::ArrayRCP<const size_t> numEntriesPerRow;
02237     size_t numEntriesForAll = 0;
02238     bool allRowsSame = true;
02239 
02240     const ptrdiff_t numRows = static_cast<ptrdiff_t> (this->getNodeNumRows ());
02241 
02242     if (! this->indicesAreAllocated ()) {
02243       if (! this->numAllocPerRow_.is_null ()) {
02244         numEntriesPerRow = this->numAllocPerRow_;
02245         allRowsSame = false; // conservatively; we don't check the array
02246       }
02247       else {
02248         numEntriesForAll = this->numAllocForAllRows_;
02249         allRowsSame = true;
02250       }
02251     }
02252     else if (! this->numRowEntries_.is_null ()) {
02253       numEntriesPerRow = numRowEntries_;
02254         allRowsSame = false; // conservatively; we don't check the array
02255     }
02256     else if (this->nodeNumAllocated_ == 0) {
02257       numEntriesForAll = 0;
02258       allRowsSame = true;
02259     }
02260     else {
02261       // left with the case that we have optimized storage. in this
02262       // case, we have to construct a list of row sizes.
02263       TEUCHOS_TEST_FOR_EXCEPTION(
02264         this->getProfileType () != StaticProfile, std::logic_error,
02265         "Tpetra::CrsGraph::getNumEntriesPerRowUpperBound: "
02266         "The graph is not StaticProfile, but storage appears to be optimized.  "
02267         "Please report this bug to the Tpetra developers.");
02268       TEUCHOS_TEST_FOR_EXCEPTION(
02269         numRows != 0 && this->rowPtrs_.size () == 0, std::logic_error,
02270         "Tpetra::CrsGraph::getNumEntriesPerRowUpperBound: "
02271         "The graph has " << numRows << " (> 0) row" << (numRows != 1 ? "s" : "")
02272         << " on the calling process, but the rowPtrs_ array has zero entries.  "
02273         "Please report this bug to the Tpetra developers.");
02274 
02275       Teuchos::ArrayRCP<size_t> numEnt;
02276       if (numRows != 0) {
02277         numEnt = Teuchos::arcp<size_t> (numRows);
02278       }
02279 
02280       // We have to iterate through the row offsets anyway, so we
02281       // might as well check whether all rows' bounds are the same.
02282       bool allRowsReallySame = false;
02283       for (ptrdiff_t i = 0; i < numRows; ++i) {
02284         numEnt[i] = this->rowPtrs_[i+1] - this->rowPtrs_[i];
02285         if (i != 0 && numEnt[i] != numEnt[i-1]) {
02286           allRowsReallySame = false;
02287         }
02288       }
02289       if (allRowsReallySame) {
02290         if (numRows == 0) {
02291           numEntriesForAll = 0;
02292         } else {
02293           numEntriesForAll = numEnt[1] - numEnt[0];
02294         }
02295         allRowsSame = true;
02296       }
02297       else {
02298         numEntriesPerRow = numEnt; // Teuchos::arcp_const_cast<const size_t> (numEnt);
02299         allRowsSame = false; // conservatively; we don't check the array
02300       }
02301     }
02302 
02303     TEUCHOS_TEST_FOR_EXCEPTION(
02304       numEntriesForAll != 0 && numEntriesPerRow.size () != 0, std::logic_error,
02305       "Tpetra::CrsGraph::getNumEntriesPerLocalRowUpperBound: "
02306       "numEntriesForAll and numEntriesPerRow are not consistent.  The former "
02307       "is " << numEntriesForAll << ", but the latter has nonzero size "
02308       << numEntriesPerRow.size () << ".  "
02309       "Please report this bug to the Tpetra developers.");
02310 
02311     boundPerLocalRow = numEntriesPerRow;
02312     boundForAllLocalRows = numEntriesForAll;
02313     boundSameForAllLocalRows = allRowsSame;
02314   }
02315 
02316   // TODO: in the future, globalAssemble() should use import/export functionality
02319   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
02320   void CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::globalAssemble()
02321   {
02322     using Teuchos::Array;
02323     using Teuchos::as;
02324     using Teuchos::Comm;
02325     using Teuchos::gatherAll;
02326     using Teuchos::ireceive;
02327     using Teuchos::isend;
02328     using Teuchos::outArg;
02329     using Teuchos::REDUCE_MAX;
02330     using Teuchos::reduceAll;
02331     using Teuchos::toString;
02332     using Teuchos::waitAll;
02333     using std::endl;
02334     using std::make_pair;
02335     using std::pair;
02336     typedef GlobalOrdinal GO;
02337     typedef typename std::map<GO, std::vector<GO> >::const_iterator NLITER;
02338     typedef typename Array<GO>::size_type size_type;
02339 
02340     const char tfecfFuncName[] = "globalAssemble"; // for exception macro
02341     RCP<const Comm<int> > comm = getComm();
02342 
02343     const int numImages = comm->getSize();
02344     const int myImageID = comm->getRank();
02345 #ifdef HAVE_TPETRA_DEBUG
02346     Teuchos::barrier (*comm);
02347 #endif // HAVE_TPETRA_DEBUG
02348 
02349     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC( ! isFillActive(), std::runtime_error,
02350       ": requires that fill is active.");
02351     // Determine if any nodes have global entries to share
02352     {
02353       size_t MyNonlocals = nonlocals_.size(), MaxGlobalNonlocals;
02354       reduceAll (*comm, REDUCE_MAX, MyNonlocals, outArg (MaxGlobalNonlocals));
02355       if (MaxGlobalNonlocals == 0) {
02356         return;  // no entries to share
02357       }
02358     }
02359 
02360     // compute a list of NLRs from nonlocals_ and use it to compute:
02361     //      IdsAndRows: a vector of (id,row) pairs
02362     //          NLR2Id: a map from NLR to the Id that owns it
02363     // globalNeighbors: a global graph of connectivity between images: globalNeighbors(i,j) indicates that j sends to i
02364     //         sendIDs: a list of all images I send to
02365     //         recvIDs: a list of all images I receive from (constructed later)
02366     Array<pair<int, GO> > IdsAndRows;
02367     std::map<GO, int> NLR2Id;
02368     Teuchos::SerialDenseMatrix<int, char> globalNeighbors;
02369     Array<int> sendIDs, recvIDs;
02370     {
02371       // nonlocals_ contains the entries we are holding for all non-local rows
02372       // we want a list of the rows for which we have data
02373       Array<GO> NLRs;
02374       std::set<GO> setOfRows;
02375       for (NLITER iter = nonlocals_.begin(); iter != nonlocals_.end(); ++iter) {
02376         setOfRows.insert(iter->first);
02377       }
02378       // copy the elements in the set into an Array
02379       NLRs.resize(setOfRows.size());
02380       std::copy(setOfRows.begin(), setOfRows.end(), NLRs.begin());
02381 
02382       // get a list of ImageIDs for the non-local rows (NLRs)
02383       Array<int> NLRIds(NLRs.size());
02384       {
02385         LookupStatus stat = rowMap_->getRemoteIndexList(NLRs(),NLRIds());
02386         int lclerror = ( stat == IDNotPresent ? 1 : 0 );
02387         int gblerror;
02388         reduceAll<int, int> (*comm, REDUCE_MAX, lclerror, outArg (gblerror));
02389         if (gblerror != 0) {
02390           const int myRank = comm->getRank ();
02391           std::ostringstream os;
02392           os << "On one or more processes in the communicator, "
02393              << "there were insertions into rows of the graph that do not "
02394              << "exist in the row Map on any process in the communicator."
02395              << endl << "This process " << myRank << " is "
02396              << (lclerror == 0 ? "not " : "") << "one of those offending "
02397              << "processes." << endl;
02398           if (lclerror != 0) {
02399             // If NLRIds[k] is -1, then NLRs[k] is a row index not in
02400             // the row Map.  Collect this list of invalid row indices
02401             // for display in the exception message.
02402             Array<GO> invalidNonlocalRows;
02403             for (size_type k = 0; k < NLRs.size (); ++k) {
02404               if (NLRIds[k] == -1) {
02405                 invalidNonlocalRows.push_back (NLRs[k]);
02406               }
02407             }
02408             const size_type numInvalid = invalidNonlocalRows.size ();
02409             os << "On this process, " << numInvalid << " nonlocal row"
02410                << (numInvalid != 1 ? "s " : " ") << " were inserted that are "
02411                << "not in the row Map on any process." << endl;
02412             // Don't print _too_ many nonlocal rows.
02413             if (numInvalid <= 100) {
02414               os << "Offending row indices: "
02415                  << toString (invalidNonlocalRows ()) << endl;
02416             }
02417           }
02418           TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
02419             gblerror != 0, std::runtime_error,
02420             ": nonlocal entries correspond to invalid rows."
02421             << endl << os.str ());
02422         }
02423       }
02424 
02425       // build up a list of neighbors, as well as a map between NLRs and Ids
02426       // localNeighbors[i] != 0 iff I have data to send to image i
02427       // put NLRs,Ids into an array of pairs
02428       IdsAndRows.reserve(NLRs.size());
02429       Array<char> localNeighbors(numImages, 0);
02430       typename Array<GO>::const_iterator nlr;
02431       typename Array<int>::const_iterator id;
02432       for (nlr = NLRs.begin(), id = NLRIds.begin();
02433            nlr != NLRs.end(); ++nlr, ++id) {
02434         NLR2Id[*nlr] = *id;
02435         localNeighbors[*id] = 1;
02436         // IdsAndRows.push_back(make_pair<int,GlobalOrdinal>(*id,*nlr));
02437         IdsAndRows.push_back(make_pair(*id,*nlr));
02438       }
02439       for (int j=0; j<numImages; ++j) {
02440         if (localNeighbors[j]) {
02441           sendIDs.push_back(j);
02442         }
02443       }
02444       // sort IdsAndRows, by Ids first, then rows
02445       std::sort(IdsAndRows.begin(),IdsAndRows.end());
02446       // gather from other nodes to form the full graph
02447       globalNeighbors.shapeUninitialized(numImages,numImages);
02448       gatherAll (*getComm(), numImages, localNeighbors.getRawPtr(),
02449                  numImages * numImages, globalNeighbors.values());
02450       // globalNeighbors at this point contains (on all images) the
02451       // connectivity between the images.
02452       // globalNeighbors(i,j) != 0 means that j sends to i/that i receives from j
02453     }
02454 
02456     // FIGURE OUT WHO IS SENDING TO WHOM AND HOW MUCH
02457     // DO THIS IN THE PROCESS OF PACKING ALL OUTGOING DATA ACCORDING TO DESTINATION ID
02459 
02460     // loop over all columns to know from which images I can expect to receive something
02461     for (int j=0; j<numImages; ++j) {
02462       if (globalNeighbors(myImageID,j)) {
02463         recvIDs.push_back(j);
02464       }
02465     }
02466     const size_t numRecvs = recvIDs.size();
02467 
02468     // we know how many we're sending to already
02469     // form a contiguous list of all data to be sent
02470     // track the number of entries for each ID
02471     Array<pair<GO, GO> > IJSendBuffer;
02472     Array<size_t> sendSizes(sendIDs.size(), 0);
02473     size_t numSends = 0;
02474     for (typename Array<pair<int, GO> >::const_iterator IdAndRow = IdsAndRows.begin();
02475          IdAndRow != IdsAndRows.end(); ++IdAndRow) {
02476       int id = IdAndRow->first;
02477       GO row = IdAndRow->second;
02478       // have we advanced to a new send?
02479       if (sendIDs[numSends] != id) {
02480         numSends++;
02481         TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(sendIDs[numSends] != id,
02482           std::logic_error, ": internal logic error. Contact Tpetra team.");
02483       }
02484       // copy data for row into contiguous storage
02485       for (typename std::vector<GO>::const_iterator j = nonlocals_[row].begin(); j != nonlocals_[row].end(); ++j)
02486       {
02487         IJSendBuffer.push_back( pair<GlobalOrdinal,GlobalOrdinal>(row,*j) );
02488         sendSizes[numSends]++;
02489       }
02490     }
02491     if (IdsAndRows.size() > 0) {
02492       numSends++; // one last increment, to make it a count instead of an index
02493     }
02494     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(as<typename Array<int>::size_type>(numSends) != sendIDs.size(), std::logic_error, ": internal logic error. Contact Tpetra team.");
02495 
02496     // don't need this data anymore
02497     nonlocals_.clear();
02498 
02500     // TRANSMIT SIZE INFO BETWEEN SENDERS AND RECEIVERS
02502 
02503     // Array of pending nonblocking communication requests.  It's OK
02504     // to mix nonblocking send and receive requests in the same
02505     // waitAll() call.
02506     Array<RCP<Teuchos::CommRequest<int> > > requests;
02507 
02508     // perform non-blocking sends: send sizes to our recipients
02509     for (size_t s = 0; s < numSends ; ++s) {
02510       // We're using a nonowning RCP because all communication
02511       // will be local to this method and the scope of our data
02512       requests.push_back (isend<int, size_t> (*comm,
02513                                               rcp (&sendSizes[s], false),
02514                                               sendIDs[s]));
02515     }
02516     // perform non-blocking receives: receive sizes from our senders
02517     Array<size_t> recvSizes (numRecvs);
02518     for (size_t r = 0; r < numRecvs; ++r) {
02519       // We're using a nonowning RCP because all communication
02520       // will be local to this method and the scope of our data
02521       requests.push_back (ireceive (*comm, rcp (&recvSizes[r], false), recvIDs[r]));
02522     }
02523     // Wait on all the nonblocking sends and receives.
02524     if (! requests.empty()) {
02525       waitAll (*comm, requests());
02526     }
02527 #ifdef HAVE_TPETRA_DEBUG
02528     Teuchos::barrier (*comm);
02529 #endif // HAVE_TPETRA_DEBUG
02530 
02531     // This doesn't necessarily deallocate the array.
02532     requests.resize (0);
02533 
02535     // NOW SEND/RECEIVE ALL ROW DATA
02537     // from the size info, build the ArrayViews into IJSendBuffer
02538     Array<ArrayView<pair<GO,GO> > > sendBuffers(numSends,null);
02539     {
02540       size_t cur = 0;
02541       for (size_t s=0; s<numSends; ++s) {
02542         sendBuffers[s] = IJSendBuffer(cur,sendSizes[s]);
02543         cur += sendSizes[s];
02544       }
02545     }
02546     // perform non-blocking sends
02547     for (size_t s=0; s < numSends ; ++s)
02548     {
02549       // We're using a nonowning RCP because all communication
02550       // will be local to this method and the scope of our data
02551       ArrayRCP<pair<GO,GO> > tmpSendBuf =
02552         arcp (sendBuffers[s].getRawPtr(), 0, sendBuffers[s].size(), false);
02553       requests.push_back (isend<int, pair<GO,GO> > (*comm, tmpSendBuf, sendIDs[s]));
02554     }
02555     // calculate amount of storage needed for receives
02556     // setup pointers for the receives as well
02557     size_t totalRecvSize = std::accumulate (recvSizes.begin(), recvSizes.end(), 0);
02558     Array<pair<GO,GO> > IJRecvBuffer (totalRecvSize);
02559     // from the size info, build the ArrayViews into IJRecvBuffer
02560     Array<ArrayView<pair<GO,GO> > > recvBuffers (numRecvs, null);
02561     {
02562       size_t cur = 0;
02563       for (size_t r=0; r<numRecvs; ++r) {
02564         recvBuffers[r] = IJRecvBuffer(cur,recvSizes[r]);
02565         cur += recvSizes[r];
02566       }
02567     }
02568     // perform non-blocking recvs
02569     for (size_t r = 0; r < numRecvs; ++r) {
02570       // We're using a nonowning RCP because all communication
02571       // will be local to this method and the scope of our data
02572       ArrayRCP<pair<GO,GO> > tmpRecvBuf =
02573         arcp (recvBuffers[r].getRawPtr(), 0, recvBuffers[r].size(), false);
02574       requests.push_back (ireceive (*comm, tmpRecvBuf, recvIDs[r]));
02575     }
02576     // perform waits
02577     if (! requests.empty()) {
02578       waitAll (*comm, requests());
02579     }
02580 #ifdef HAVE_TPETRA_DEBUG
02581     Teuchos::barrier (*comm);
02582 #endif // HAVE_TPETRA_DEBUG
02583 
02585     // NOW PROCESS THE RECEIVED ROW DATA
02587     // TODO: instead of adding one entry at a time, add one row at a time.
02588     //       this requires resorting; they arrived sorted by sending node,
02589     //       so that entries could be non-contiguous if we received
02590     //       multiple entries for a particular row from different processors.
02591     //       it also requires restoring the data, which may make it not worth the trouble.
02592     for (typename Array<pair<GO,GO> >::const_iterator ij = IJRecvBuffer.begin();
02593          ij != IJRecvBuffer.end(); ++ij)
02594     {
02595       insertGlobalIndicesFiltered (ij->first, tuple<GO> (ij->second));
02596     }
02597     checkInternalState();
02598   }
02599 
02600 
02603   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
02604   void
02605   CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::
02606   resumeFill (const RCP<ParameterList> &params)
02607   {
02608     const char tfecfFuncName[] = "resumeFill";
02609     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(! hasRowInfo(), std::runtime_error,
02610       ": Sorry, you cannot resume fill of the CrsGraph, since the graph's row "
02611       "information was deleted in fillComplete().");
02612 
02613 #ifdef HAVE_TPETRA_DEBUG
02614     Teuchos::barrier( *rowMap_->getComm() );
02615 #endif // HAVE_TPETRA_DEBUG
02616     clearGlobalConstants();
02617     lclGraph_ = null;
02618     if (params != null) this->setParameterList (params);
02619     lowerTriangular_  = false;
02620     upperTriangular_  = false;
02621     // either still sorted/merged or initially sorted/merged
02622     indicesAreSorted_ = true;
02623     noRedundancies_ = true;
02624     fillComplete_ = false;
02625 #ifdef HAVE_TPETRA_DEBUG
02626     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
02627       ! isFillActive() || isFillComplete(), std::logic_error,
02628       "::resumeFill(): At end of method, either fill is not active or fill is "
02629       "complete.  This violates stated post-conditions.  Please report this bug "
02630       "to the Tpetra developers.");
02631 #endif // HAVE_TPETRA_DEBUG
02632   }
02633 
02634 
02637   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
02638   void CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::fillComplete(const RCP<ParameterList> &params)
02639   {
02640     fillComplete(rowMap_,rowMap_,params);
02641   }
02642 
02643 
02646   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
02647   void
02648   CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::
02649   fillComplete (const RCP<const map_type>& domainMap,
02650                 const RCP<const map_type>& rangeMap,
02651                 const RCP<ParameterList> & params)
02652   {
02653     const char tfecfFuncName[] = "fillComplete";
02654 
02655 #ifdef HAVE_TPETRA_DEBUG
02656     rowMap_->getComm ()->barrier ();
02657 #endif // HAVE_TPETRA_DEBUG
02658 
02659     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC( ! isFillActive() || isFillComplete(),
02660       std::runtime_error, ": Graph fill state must be active (isFillActive() "
02661       "must be true) before calling fillComplete().");
02662 
02663     const int numProcs = getComm ()->getSize ();
02664 
02665     //
02666     // Read and set parameters
02667     //
02668 
02669     // Does the caller want to sort remote GIDs (within those owned by
02670     // the same process) in makeColMap()?
02671     if (! params.is_null ()) {
02672       if (params->isParameter ("sort column map ghost gids")) {
02673         sortGhostsAssociatedWithEachProcessor_ =
02674           params->get<bool> ("sort column map ghost gids",
02675                              sortGhostsAssociatedWithEachProcessor_);
02676       }
02677       else if (params->isParameter ("Sort column Map ghost GIDs")) {
02678         sortGhostsAssociatedWithEachProcessor_ =
02679           params->get<bool> ("Sort column Map ghost GIDs",
02680                              sortGhostsAssociatedWithEachProcessor_);
02681       }
02682     }
02683 
02684     // If true, the caller promises that no process did nonlocal
02685     // changes since the last call to fillComplete.
02686     bool assertNoNonlocalInserts = false;
02687     if (! params.is_null ()) {
02688       assertNoNonlocalInserts =
02689         params->get<bool> ("No Nonlocal Changes", assertNoNonlocalInserts);
02690     }
02691 
02692     //
02693     // Allocate indices, if they haven't already been allocated
02694     //
02695     if (! indicesAreAllocated ()) {
02696       if (hasColMap ()) {
02697         // We have a column Map, so use local indices.
02698         allocateIndices (LocalIndices);
02699       } else {
02700         // We don't have a column Map, so use global indices.
02701         allocateIndices (GlobalIndices);
02702       }
02703     }
02704 
02705     //
02706     // Do global assembly, if requested and if the communicator
02707     // contains more than one process.
02708     //
02709     const bool mayNeedGlobalAssemble = ! assertNoNonlocalInserts && numProcs > 1;
02710     if (mayNeedGlobalAssemble) {
02711       // This first checks if we need to do global assembly.
02712       // The check costs a single all-reduce.
02713       globalAssemble ();
02714     }
02715     else {
02716       TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
02717         numProcs > 1 && nonlocals_.size() > 0, std::runtime_error,
02718         ":" << std::endl << "The graph's communicator contains only one "
02719         "process, but there are nonlocal entries.  " << std::endl <<
02720         "This probably means that invalid entries were added to the graph.");
02721     }
02722 
02723     // Set domain and range Map.  This may clear the Import / Export
02724     // objects if the new Maps differ from any old ones.
02725     setDomainRangeMaps (domainMap, rangeMap);
02726 
02727     // If the graph does not already have a column Map (either from
02728     // the user constructor calling the version of the constructor
02729     // that takes a column Map, or from a previous fillComplete call),
02730     // then create it.
02731     if (! hasColMap ()) {
02732       makeColMap ();
02733     }
02734 
02735     // Make indices local, if they aren't already.
02736     // The method doesn't do any work if the indices are already local.
02737     makeIndicesLocal ();
02738 
02739     if (! isSorted ()) {
02740       // If this process has no indices, then CrsGraph considers it
02741       // already trivially sorted.  Thus, this method need not be
02742       // called on all processes in the row Map's communicator.
02743       sortAllIndices ();
02744     }
02745     if (! isMerged ()) {
02746       mergeAllIndices ();
02747     }
02748     makeImportExport(); // Make Import and Export objects, if necessary
02749     computeGlobalConstants ();
02750     fillLocalGraph (params);
02751     fillComplete_ = true;
02752 
02753 #ifdef HAVE_TPETRA_DEBUG
02754     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
02755       isFillActive() == true || isFillComplete() == false, std::logic_error,
02756       ": Violated stated post-conditions. Please contact Tpetra team.");
02757 #endif
02758 
02759     checkInternalState ();
02760   }
02761 
02762 
02765   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
02766   void CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::expertStaticFillComplete(const RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > & domainMap,
02767                                                                                        const RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > & rangeMap,
02768                                                                                        const RCP<const Import<LocalOrdinal,GlobalOrdinal,Node> > &importer,
02769                                                                                        const RCP<const Export<LocalOrdinal,GlobalOrdinal,Node> > &exporter,
02770                                                                                        const RCP<ParameterList> &params)
02771   {
02772 #ifdef HAVE_TPETRA_DEBUG
02773     rowMap_->getComm ()->barrier ();
02774 #endif // HAVE_TPETRA_DEBUG
02775     const char tfecfFuncName[] = "expertStaticFillComplete()";
02776     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC( isFillComplete() == true || hasColMap() == false,
02777                                            std::runtime_error, ": fillComplete cannot have already been called and a ColMap is required.");
02778 
02779     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC( getNodeNumRows() > 0 && rowPtrs_==Teuchos::null,
02780                                            std::runtime_error, ": a matrix will getNodeNumRows()>0 requires rowptr to be set.");
02781 
02782     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC( domainMap == Teuchos::null || rangeMap == Teuchos::null,
02783                                            std::runtime_error, ": requires a non-null domainMap & rangeMap.");
02784 
02785     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC( pftype_ !=StaticProfile,
02786                                            std::runtime_error, ": requires StaticProfile.");
02787 
02788     // Note: We don't need to do the following things which are normally done in fillComplete:
02789     // allocateIndices, globalAssemble, makeColMap, makeIndicesLocal, sortAllIndices, mergeAllIndices
02790 
02791     // Note: Need to do this so computeGlobalConstants & fillLocalGraph work
02792     nodeNumEntries_ = nodeNumAllocated_ = rowPtrs_[getNodeNumRows()];
02793 
02794     // Constants from allocateIndices
02795     numAllocForAllRows_  = 0;
02796     numAllocPerRow_      = null;
02797     indicesAreAllocated_ = true;
02798 
02799     // Constants from makeIndicesLocal
02800     indicesAreLocal_  = true;
02801     indicesAreGlobal_ = false;
02802 
02803     // set domain/range map: may clear the import/export objects
02804     setDomainRangeMaps(domainMap,rangeMap);
02805 
02806     // Presume the user sorted and merged the arrays first
02807     indicesAreSorted_ = true;
02808     noRedundancies_ = true;
02809 
02810     // makeImportExport won't create a new importer/exporter if I set one here first.
02811     importer_=Teuchos::null;
02812     exporter_=Teuchos::null;
02813     if(importer != Teuchos::null) {
02814       TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(!importer->getSourceMap()->isSameAs(*getDomainMap()) || !importer->getTargetMap()->isSameAs(*getColMap()),
02815                                             std::invalid_argument,": importer does not match matrix maps.");
02816       importer_ = importer;
02817 
02818     }
02819     if(exporter != Teuchos::null) {
02820       TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(!exporter->getSourceMap()->isSameAs(*getRowMap()) || !exporter->getTargetMap()->isSameAs(*getRangeMap()),
02821                                             std::invalid_argument,": exporter does not match matrix maps.");
02822       exporter_ = exporter;
02823     }
02824     makeImportExport();
02825 
02826     // Compute the constants
02827     computeGlobalConstants();
02828 
02829     // Since we have a StaticProfile, fillLocalGraph will do the right thing...
02830     fillLocalGraph(params);
02831     fillComplete_ = true;
02832 
02833 #ifdef HAVE_TPETRA_DEBUG
02834     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC( isFillActive() == true || isFillComplete() == false, std::logic_error, ": Violated stated post-conditions. Please contact Tpetra team.");
02835 #endif
02836     checkInternalState();
02837   }
02838 
02839 
02840 
02843   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
02844   void CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::fillLocalGraph(const RCP<ParameterList> &params)
02845   {
02846     const size_t numRows = getNodeNumRows();
02847     ArrayRCP<size_t> ptrs;
02848     ArrayRCP<LocalOrdinal> inds;
02849     bool requestOptimizedStorage = true;
02850     if (params != null && params->get("Optimize Storage",true) == false) requestOptimizedStorage = false;
02851     if (getProfileType() == DynamicProfile) {
02852       // 2d -> 1d packed
02853       ptrs = LocalMatOps::allocRowPtrs( getRowMap()->getNode(), numRowEntries_() );
02854       inds = LocalMatOps::template allocStorage<LocalOrdinal>( getRowMap()->getNode(), ptrs() );
02855       for (size_t row=0; row < numRows; ++row) {
02856         const size_t numentrs = numRowEntries_[row];
02857         std::copy( lclInds2D_[row].begin(), lclInds2D_[row].begin()+numentrs, inds+ptrs[row] );
02858       }
02859     }
02860     else if (getProfileType() == StaticProfile) {
02861       // 1d non-packed -> 1d packed
02862       if (nodeNumEntries_ != nodeNumAllocated_) {
02863         ptrs = LocalMatOps::allocRowPtrs( getRowMap()->getNode(), numRowEntries_() );
02864         inds = LocalMatOps::template allocStorage<LocalOrdinal>( getRowMap()->getNode(), ptrs() );
02865         for (size_t row=0; row < numRows; ++row) {
02866           const size_t numentrs = numRowEntries_[row];
02867           std::copy( lclInds1D_+rowPtrs_[row], lclInds1D_+rowPtrs_[row]+numentrs, inds+ptrs[row] );
02868         }
02869       }
02870       else {
02871         inds = lclInds1D_;
02872         ptrs = rowPtrs_;
02873       }
02874     }
02875     // can we ditch the old allocations for the packed one?
02876     if ( requestOptimizedStorage ) {
02877       lclInds2D_ = null;
02878       numRowEntries_ = null;
02879       // keep the new stuff
02880       lclInds1D_ = inds;
02881       rowPtrs_ = ptrs;
02882       nodeNumAllocated_ = nodeNumEntries_;
02883       pftype_ = StaticProfile;
02884     }
02885     // build the local graph, hand over the indices
02886     RCP<ParameterList> lclparams;
02887     if (params == null) lclparams = parameterList();
02888     else                lclparams = sublist(params,"Local Graph");
02889     lclGraph_ = rcp( new local_graph_type( getRowMap()->getNodeNumElements(), getColMap()->getNodeNumElements(), getRowMap()->getNode(), lclparams ) );
02890     lclGraph_->setStructure(ptrs,inds);
02891     ptrs = null;
02892     inds = null;
02893     // finalize local graph
02894     Teuchos::EDiag diag = ( getNodeNumDiags() < getNodeNumRows() ? Teuchos::UNIT_DIAG : Teuchos::NON_UNIT_DIAG );
02895     Teuchos::EUplo uplo = Teuchos::UNDEF_TRI;
02896     if      (isUpperTriangular()) uplo = Teuchos::UPPER_TRI;
02897     else if (isLowerTriangular()) uplo = Teuchos::LOWER_TRI;
02898     LocalMatOps::finalizeGraph(uplo,diag,*lclGraph_,params);
02899   }
02900 
02901 
02904   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
02905   void
02906   CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::
02907   replaceColMap (const Teuchos::RCP<const map_type>& newColMap)
02908   {
02909     // NOTE: This safety check matches the code, but not the documentation of Crsgraph
02910     const char tfecfFuncName[] = "replaceColMap";
02911     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(isLocallyIndexed() || isGloballyIndexed(),  std::runtime_error, " requires matching maps and non-static graph.");
02912     colMap_ = newColMap;
02913   }
02914 
02915 
02918   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
02919   void
02920   CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::
02921   replaceDomainMapAndImporter (const Teuchos::RCP<const map_type>& newDomainMap,
02922                                const Teuchos::RCP<const import_type>& newImporter)
02923   {
02924     const char prefix[] = "Tpetra::CrsGraph::replaceDomainMapAndImporter: ";
02925     TEUCHOS_TEST_FOR_EXCEPTION(
02926       colMap_.is_null (), std::invalid_argument, prefix << "You may not call "
02927       "this method unless the graph already has a column Map.");
02928     TEUCHOS_TEST_FOR_EXCEPTION(
02929       newDomainMap.is_null (), std::invalid_argument,
02930       prefix << "The new domain Map must be nonnull.");
02931 
02932 #ifdef HAVE_TPETRA_DEBUG
02933     if (newImporter.is_null ()) {
02934       // It's not a good idea to put expensive operations in a macro
02935       // clause, even if they are side effect - free, because macros
02936       // don't promise that they won't evaluate their arguments more
02937       // than once.  It's polite for them to do so, but not required.
02938       const bool colSameAsDom = colMap_->isSameAs (*newDomainMap);
02939       TEUCHOS_TEST_FOR_EXCEPTION(
02940         colSameAsDom, std::invalid_argument, "If the new Import is null, "
02941         "then the new domain Map must be the same as the current column Map.");
02942     }
02943     else {
02944       const bool colSameAsTgt =
02945         colMap_->isSameAs (* (newImporter->getTargetMap ()));
02946       const bool newDomSameAsSrc =
02947         newDomainMap->isSameAs (* (newImporter->getSourceMap ()));
02948       TEUCHOS_TEST_FOR_EXCEPTION(
02949         ! colSameAsTgt || ! newDomSameAsSrc, std::invalid_argument, "If the "
02950         "new Import is nonnull, then the current column Map must be the same "
02951         "as the new Import's target Map, and the new domain Map must be the "
02952         "same as the new Import's source Map.");
02953     }
02954 #endif // HAVE_TPETRA_DEBUG
02955 
02956     domainMap_ = newDomainMap;
02957     importer_ = Teuchos::rcp_const_cast<import_type> (newImporter);
02958   }
02959 
02960 
02963   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
02964   const RCP<const typename LocalMatOps::template graph<LocalOrdinal,Node>::graph_type>
02965   CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getLocalGraph() const
02966   {
02967     return lclGraph_;
02968   }
02969 
02970 
02973   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
02974   const RCP<typename LocalMatOps::template graph<LocalOrdinal,Node>::graph_type>
02975   CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getLocalGraphNonConst()
02976   {
02977     return lclGraph_;
02978   }
02979 
02980 
02983   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
02984   void CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::computeGlobalConstants()
02985   {
02986     using Teuchos::as;
02987     using Teuchos::outArg;
02988     using Teuchos::reduceAll;
02989     typedef LocalOrdinal LO;
02990     typedef GlobalOrdinal GO;
02991     typedef global_size_t GST;
02992 
02993 #ifdef HAVE_TPETRA_DEBUG
02994     TEUCHOS_TEST_FOR_EXCEPTION(! hasColMap(), std::logic_error, "Tpetra::"
02995       "CrsGraph::computeGlobalConstants: At this point, the graph should have "
02996       "a column Map, but it does not.  Please report this bug to the Tpetra "
02997       "developers.");
02998 #endif // HAVE_TPETRA_DEBUG
02999 
03000     // If necessary, (re)compute the local constants: nodeNumDiags_,
03001     // lowerTriangular_, upperTriangular_, and nodeMaxNumRowEntries_.
03002     if (! haveLocalConstants_) {
03003       // We have actually already computed nodeNumEntries_.
03004       // nodeNumEntries_ gets updated by methods that insert or remove
03005       // indices (including setAllIndices and
03006       // expertStaticFillComplete).  Before fillComplete, its count
03007       // may include duplicate column indices in the same row.
03008       // However, mergeRowIndices and mergeRowIndicesAndValues both
03009       // subtract off merged indices in each row from the total count.
03010       // Thus, nodeNumEntries_ _should_ be accurate at this point,
03011       // meaning that we don't have to re-count it here.
03012 
03013       // Reset local properties
03014       upperTriangular_ = true;
03015       lowerTriangular_ = true;
03016       nodeMaxNumRowEntries_ = 0;
03017       nodeNumDiags_         = 0;
03018 
03019       // At this point, we know that we have both a row Map and a column Map.
03020       const Map<LO,GO,Node>& rowMap = *rowMap_;
03021       const Map<LO,GO,Node>& colMap = *colMap_;
03022 
03023       // Go through all the entries of the graph.  Count the number of
03024       // diagonal elements we encounter, and figure out whether the
03025       // graph is lower or upper triangular.  Diagonal elements are
03026       // determined using global indices, with respect to the whole
03027       // graph.  However, lower or upper triangularity is a local
03028       // property, and is determined using local indices.
03029       //
03030       // At this point, indices have already been sorted in each row.
03031       // That makes finding out whether the graph is lower / upper
03032       // triangular easier.
03033       if (indicesAreAllocated () && nodeNumAllocated_ > 0) {
03034         const size_t numLocalRows = getNodeNumRows ();
03035         for (size_t localRow = 0; localRow < numLocalRows; ++localRow) {
03036           const GO globalRow = rowMap.getGlobalElement (localRow);
03037           // Find the local (column) index for the diagonal element.
03038           const LO rlcid = colMap.getLocalElement (globalRow);
03039           RowInfo rowInfo = getRowInfo (localRow);
03040           ArrayView<const LO> rview = getLocalView (rowInfo);
03041           typename ArrayView<const LO>::iterator beg, end, cur;
03042           beg = rview.begin();
03043           end = beg + rowInfo.numEntries;
03044           if (beg != end) {
03045             for (cur = beg; cur != end; ++cur) {
03046               // is this the diagonal?
03047               if (rlcid == *cur) ++nodeNumDiags_;
03048             }
03049             // Local column indices are sorted in each row.  That means
03050             // the smallest column index in this row (on this process)
03051             // is *beg, and the largest column index in this row (on
03052             // this process) is *(end - 1).  We know that end - 1 is
03053             // valid because beg != end.
03054             const size_t smallestCol = as<size_t> (*beg);
03055             const size_t largestCol = as<size_t> (*(end - 1));
03056 
03057             if (smallestCol < localRow) {
03058               upperTriangular_ = false;
03059             }
03060             if (localRow < largestCol) {
03061               lowerTriangular_ = false;
03062             }
03063           }
03064           // Update the max number of entries over all rows.
03065           nodeMaxNumRowEntries_ = std::max (nodeMaxNumRowEntries_, rowInfo.numEntries);
03066         }
03067       }
03068       haveLocalConstants_ = true;
03069     } // if my process doesn't have local constants
03070 
03071     // Compute global constants from local constants.  Processes that
03072     // already have local constants still participate in the
03073     // all-reduces, using their previously computed values.
03074     if (haveGlobalConstants_ == false) {
03075       // Promote all the nodeNum* and nodeMaxNum* quantities from
03076       // size_t to global_size_t, when doing the all-reduces for
03077       // globalNum* / globalMaxNum* results.
03078       //
03079       // FIXME (mfh 07 May 2013) Unfortunately, we either have to do
03080       // this in two all-reduces (one for the sum and the other for
03081       // the max), or use a custom MPI_Op that combines the sum and
03082       // the max.  The latter might even be slower than two
03083       // all-reduces on modern network hardware.  It would also be a
03084       // good idea to use nonblocking all-reduces (MPI 3), so that we
03085       // don't have to wait around for the first one to finish before
03086       // starting the second one.
03087       GST lcl[2], gbl[2];
03088       lcl[0] = as<GST> (nodeNumEntries_);
03089       lcl[1] = as<GST> (nodeNumDiags_);
03090       reduceAll<int,GST> (*getComm (), Teuchos::REDUCE_SUM,
03091                           2, lcl, gbl);
03092       globalNumEntries_ = gbl[0];
03093       globalNumDiags_   = gbl[1];
03094       reduceAll<int,GST> (*getComm (), Teuchos::REDUCE_MAX,
03095                           as<GST> (nodeMaxNumRowEntries_),
03096                           outArg (globalMaxNumRowEntries_));
03097       haveGlobalConstants_ = true;
03098     }
03099   }
03100 
03101 
03104   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
03105   void
03106   CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::makeIndicesLocal ()
03107   {
03108     typedef LocalOrdinal LO;
03109     typedef GlobalOrdinal GO;
03110 #ifdef HAVE_TPETRA_DEBUG
03111     const char tfecfFuncName[] = "makeIndicesLocal";
03112 #endif // HAVE_TPETRA_DEBUG
03113 
03114     // Transform indices to local index space
03115     const size_t nlrs = getNodeNumRows ();
03116     if (isGloballyIndexed () && nlrs > 0) {
03117       // allocate data for local indices
03118       if (getProfileType() == StaticProfile) {
03119         // If GO and LO are the same size, we can reuse the existing
03120         // array of 1-D index storage to convert column indices from
03121         // GO to LO.  Otherwise, we'll just allocate a new buffer.
03122         if (nodeNumAllocated_ && sizeof (LO) == sizeof (GO)) {
03123           lclInds1D_ = arcp_reinterpret_cast<LO> (gblInds1D_).persistingView (0, nodeNumAllocated_);
03124         }
03125         else {
03126           lclInds1D_ = LocalMatOps::template allocStorage<LO> (getRowMap ()->getNode (), rowPtrs_ ());
03127         }
03128         for (size_t r = 0; r < nlrs; ++r) {
03129           const size_t offset   = rowPtrs_[r];
03130           const size_t numentry = numRowEntries_[r];
03131           for (size_t j = 0; j < numentry; ++j) {
03132             const GO gid = gblInds1D_[offset + j];
03133             const LO lid = colMap_->getLocalElement (gid);
03134             lclInds1D_[offset + j] = lid;
03135 #ifdef HAVE_TPETRA_DEBUG
03136             TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
03137               lclInds1D_[offset + j] == Teuchos::OrdinalTraits<LO>::invalid(),
03138               std::logic_error,
03139               ": In local row r=" << r << ", global column " << gid << " is "
03140               "not in the column Map.  This should never happen.  Please "
03141               "report this bug to the Tpetra developers.");
03142 #endif // HAVE_TPETRA_DEBUG
03143           }
03144         }
03145         // We've converted column indices from global to local, so we
03146         // can deallocate the global column indices (which we know are
03147         // in 1-D storage, because the graph has static profile).
03148         gblInds1D_ = null;
03149       }
03150       else {  // the graph has dynamic profile (2-D index storage)
03151         lclInds2D_ = arcp<Array<LO> > (nlrs);
03152         for (size_t r = 0; r < getNodeNumRows (); ++r) {
03153           if (! gblInds2D_[r].empty ()) {
03154             const GO* const ginds = gblInds2D_[r].getRawPtr ();
03155             const size_t rna = gblInds2D_[r].size ();
03156             const size_t numentry = numRowEntries_[r];
03157             lclInds2D_[r].resize (rna);
03158             LO* const linds = lclInds2D_[r].getRawPtr ();
03159             for (size_t j = 0; j < numentry; ++j) {
03160               GO gid = ginds[j];
03161               LO lid = colMap_->getLocalElement (gid);
03162               linds[j] = lid;
03163 #ifdef HAVE_TPETRA_DEBUG
03164               TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
03165                 linds[j] == Teuchos::OrdinalTraits<LO>::invalid(),
03166                 std::logic_error,
03167                 ": Global column ginds[j=" << j << "]=" << ginds[j]
03168                 << " of local row r=" << r << " is not in the column Map.  "
03169                 "This should never happen.  Please report this bug to the "
03170                 "Tpetra developers.");
03171 #endif // HAVE_TPETRA_DEBUG
03172             }
03173           }
03174         }
03175         gblInds2D_ = null;
03176       }
03177     } // globallyIndexed() && nlrs > 0
03178     indicesAreLocal_  = true;
03179     indicesAreGlobal_ = false;
03180     checkInternalState();
03181   }
03182 
03183 
03186   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
03187   void CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::sortAllIndices()
03188   {
03189     TEUCHOS_TEST_FOR_EXCEPT(isGloballyIndexed()==true);   // this should be called only after makeIndicesLocal()
03190     if (isSorted() == false) {
03191       for (size_t row=0; row < getNodeNumRows(); ++row) {
03192         RowInfo rowInfo = getRowInfo(row);
03193         sortRowIndices(rowInfo);
03194       }
03195     }
03196     // we just sorted every row
03197     indicesAreSorted_ = true;
03198   }
03199 
03200 
03203   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
03204   void CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::makeColMap ()
03205   {
03206     using std::endl;
03207     using Teuchos::REDUCE_MAX;
03208     using Teuchos::reduceAll;
03209     typedef LocalOrdinal LO;
03210     typedef GlobalOrdinal GO;
03211     const char tfecfFuncName[] = "makeColMap";
03212 
03213     if (hasColMap ()) { // The graph already has a column Map.
03214       // FIXME (mfh 26 Feb 2013): This currently prevents structure
03215       // changes that affect the column Map.
03216       return;
03217     }
03218 
03219     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
03220       isLocallyIndexed (), std::runtime_error,
03221       ": The graph is locally indexed.  Calling makeColMap() to make the "
03222       "column Map requires that the graph be globally indexed.");
03223 
03224     // After the calling process is done going through all of the rows
03225     // it owns, myColumns will contain the list of indices owned by
03226     // this process in the column Map.
03227     Array<GO> myColumns;
03228 
03229     // If we reach this point, we don't have a column Map yet, so the
03230     // graph can't be locally indexed.  Thus, isGloballyIndexed() ==
03231     // false means that the graph is empty on this process, so
03232     // myColumns will be left empty.
03233     if (isGloballyIndexed ()) {
03234       // Go through all the rows, finding the populated column indices.
03235       //
03236       // Our final list of indices for the column Map constructor will
03237       // have the following properties (all of which are with respect
03238       // to the calling process):
03239       //
03240       // 1. Indices in the domain Map go first.
03241       // 2. Indices not in the domain Map follow, ordered first
03242       //    contiguously by their owning process rank (in the domain
03243       //    Map), then in increasing order within that.
03244       // 3. No duplicate indices.
03245       //
03246       // This imitates the ordering used by Aztec(OO) and Epetra.
03247       // Storing indices owned by the same process (in the domain Map)
03248       // contiguously permits the use of contiguous send and receive
03249       // buffers.
03250       //
03251       // We begin by partitioning the column indices into "local" GIDs
03252       // (owned by the domain Map) and "remote" GIDs (not owned by the
03253       // domain Map).  We use the same order for local GIDs as the
03254       // domain Map, so we track them in place in their array.  We use
03255       // an std::set (RemoteGIDSet) to keep track of remote GIDs, so
03256       // that we don't have to merge duplicates later.
03257       const LO LINV = Teuchos::OrdinalTraits<LO>::invalid ();
03258       size_t numLocalColGIDs = 0, numRemoteColGIDs = 0;
03259 
03260       // GIDisLocal[lid] == 0 if and only if local index lid in the
03261       // domain Map is remote (not local).
03262       Array<char> GIDisLocal (domainMap_->getNodeNumElements (), 0);
03263       std::set<GO> RemoteGIDSet;
03264       // This preserves the not-sorted Epetra order of GIDs.
03265       std::vector<GO> RemoteGIDUnorderedVector;
03266       const size_t myNumRows = getNodeNumRows ();
03267       for (size_t r = 0; r < myNumRows; ++r) {
03268         RowInfo rowinfo = getRowInfo (r);
03269         if (rowinfo.numEntries > 0) {
03270           // FIXME (mfh 03 Mar 2013) It's a bit puzzling to me why the
03271           // ArrayView that getGlobalView() returns doesn't return
03272           // rowinfo.numEntries entries.
03273           ArrayView<const GO> rowGids = getGlobalView (rowinfo);
03274           rowGids = rowGids (0, rowinfo.numEntries);
03275 
03276           for (size_t k = 0; k < rowinfo.numEntries; ++k) {
03277             const GO gid = rowGids[k];
03278             const LO lid = domainMap_->getLocalElement (gid);
03279             if (lid != LINV) {
03280               const char alreadyFound = GIDisLocal[lid];
03281               if (alreadyFound == 0) {
03282                 GIDisLocal[lid] = 1;
03283                 ++numLocalColGIDs;
03284               }
03285             }
03286             else {
03287               const bool notAlreadyFound = RemoteGIDSet.insert (gid).second;
03288               if (notAlreadyFound) { // gid did not exist in the set before
03289                 if (sortGhostsAssociatedWithEachProcessor_) {
03290                   RemoteGIDUnorderedVector.push_back (gid);
03291                 }
03292                 ++numRemoteColGIDs;
03293               }
03294             }
03295           } // for each entry k in row r
03296         } // if row r contains > 0 entries
03297       } // for each locally owned row r
03298 
03299       // Possible short-circuit for serial scenario:
03300       //
03301       // If all domain GIDs are present as column indices, then set
03302       // ColMap=DomainMap.  By construction, LocalGIDs is a subset of
03303       // DomainGIDs.
03304       //
03305       // If we have
03306       //   * Number of remote GIDs is 0, so that ColGIDs == LocalGIDs,
03307       // and
03308       //   * Number of local GIDs is number of domain GIDs
03309       // then
03310       //   * LocalGIDs \subset DomainGIDs && size(LocalGIDs) ==
03311       //     size(DomainGIDs) => DomainGIDs == LocalGIDs == ColGIDs
03312       // on the calling process.
03313       //
03314       // We will concern ourselves only with the special case of a
03315       // serial DomainMap, obviating the need for communication.
03316       //
03317       // If
03318       //   * DomainMap has a serial communicator
03319       // then we can set the column Map as the domain Map
03320       // return. Benefit: this graph won't need an Import object
03321       // later.
03322       //
03323       // Note, for a serial domain map, there can be no RemoteGIDs,
03324       // because there are no remote processes.  Likely explanations
03325       // for this are:
03326       //  * user submitted erroneous column indices
03327       //  * user submitted erroneous domain Map
03328       if (domainMap_->getComm ()->getSize () == 1) {
03329         TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
03330           numRemoteColGIDs != 0, std::runtime_error,
03331           ": " << numRemoteColGIDs << " column "
03332           << (numRemoteColGIDs != 1 ? "indices are" : "index is")
03333           << " not in the domain Map." << endl
03334           << "Either these indices are invalid or the domain Map is invalid."
03335           << endl << "Remember that nonsquare matrices, or matrices where the "
03336           "row and range Maps are different, require calling the version of "
03337           "fillComplete that takes the domain and range Maps as input.");
03338         if (numLocalColGIDs == domainMap_->getNodeNumElements()) {
03339           colMap_ = domainMap_;
03340           checkInternalState ();
03341           return;
03342         }
03343       }
03344 
03345       // Populate myColumns with a list of all column GIDs.  Put
03346       // locally owned (in the domain Map) GIDs at the front: they
03347       // correspond to "same" and "permuted" entries between the
03348       // column Map and the domain Map.  Put remote GIDs at the back.
03349       myColumns.resize (numLocalColGIDs + numRemoteColGIDs);
03350       // get pointers into myColumns for each part
03351       ArrayView<GO> LocalColGIDs  = myColumns (0, numLocalColGIDs);
03352       ArrayView<GO> RemoteColGIDs = myColumns (numLocalColGIDs, numRemoteColGIDs);
03353 
03354       // Copy the remote GIDs into myColumns
03355       if (sortGhostsAssociatedWithEachProcessor_) {
03356         std::copy (RemoteGIDSet.begin(), RemoteGIDSet.end(),
03357                    RemoteColGIDs.begin());
03358       } else {
03359         std::copy (RemoteGIDUnorderedVector.begin(),
03360                    RemoteGIDUnorderedVector.end(), RemoteColGIDs.begin());
03361       }
03362 
03363       // Make a list of process ranks corresponding to the remote GIDs.
03364       Array<int> RemoteImageIDs (numRemoteColGIDs);
03365       // Look up the remote process' ranks in the domain Map.
03366       {
03367         const LookupStatus stat =
03368           domainMap_->getRemoteIndexList (RemoteColGIDs, RemoteImageIDs ());
03369 #ifdef HAVE_TPETRA_DEBUG
03370         // If any process returns IDNotPresent, then at least one of
03371         // the remote indices was not present in the domain Map.  This
03372         // means that the Import object cannot be constructed, because
03373         // of incongruity between the column Map and domain Map.
03374         // This has two likely causes:
03375         //   - The user has made a mistake in the column indices
03376         //   - The user has made a mistake with respect to the domain Map
03377         const int missingID_lcl = (stat == IDNotPresent ? 1 : 0);
03378         int missingID_gbl = 0;
03379         reduceAll<int, int> (*getComm (), REDUCE_MAX, missingID_lcl,
03380                              outArg (missingID_gbl));
03381         TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
03382           missingID_gbl == 1, std::runtime_error,
03383           ": Some column indices are not in the domain Map." << endl
03384           << "Either these column indices are invalid or the domain Map is "
03385           "invalid." << endl << "Likely cause: For a nonsquare matrix, you "
03386           "must give the domain and range Maps as input to fillComplete.");
03387 #else
03388         (void) stat; // forestall compiler warning for unused variable
03389 #endif // HAVE_TPETRA_DEBUG
03390       }
03391       // Sort incoming remote column indices so that all columns
03392       // coming from a given remote process are contiguous.  This
03393       // means the Import's Distributor doesn't need to reorder data.
03394       sort2 (RemoteImageIDs.begin(), RemoteImageIDs.end(), RemoteColGIDs.begin());
03395 
03396       // Copy the local GIDs into myColumns. Two cases:
03397       // 1. If the number of Local column GIDs is the same as the
03398       //    number of Local domain GIDs, we can simply read the domain
03399       //    GIDs into the front part of ColIndices (see logic above
03400       //    from the serial short circuit case)
03401       // 2. We step through the GIDs of the DomainMap, checking to see
03402       //    if each domain GID is a column GID.  We want to do this to
03403       //    maintain a consistent ordering of GIDs between the columns
03404       //    and the domain.
03405 
03406       // FIXME (mfh 03 Mar 2013) It's common that the domain Map is
03407       // contiguous.  It would be more efficient in that case to avoid
03408       // calling getNodeElementList(), since that permanently
03409       // constructs and caches the GID list in the contiguous Map.
03410       ArrayView<const GO> domainElts = domainMap_->getNodeElementList ();
03411       const size_t numDomainElts = domainMap_->getNodeNumElements ();
03412       if (numLocalColGIDs == numDomainElts) {
03413         // If the number of locally owned GIDs are the same as the
03414         // number of local domain Map elements, then the local domain
03415         // Map elements are the same as the locally owned GIDs.
03416         std::copy (domainElts.begin(), domainElts.end(), LocalColGIDs.begin());
03417       }
03418       else {
03419         // Count the number of locally owned GIDs, both to keep track
03420         // of the current array index, and as a sanity check.
03421         size_t numLocalCount = 0;
03422         for (size_t i = 0; i < numDomainElts; ++i) {
03423           if (GIDisLocal[i]) {
03424             LocalColGIDs[numLocalCount++] = domainElts[i];
03425           }
03426         }
03427         TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
03428           numLocalCount != numLocalColGIDs, std::logic_error,
03429           ": numLocalCount = " << numLocalCount << " != numLocalColGIDs = "
03430           << numLocalColGIDs << ".  This should never happen.  Please report "
03431           "this bug to the Tpetra developers.");
03432       }
03433 
03434       // FIXME (mfh 03 Apr 2013) Now would be a good time to use the
03435       // information we collected above to construct the Import.  In
03436       // particular, building an Import requires:
03437       //
03438       // 1. numSameIDs (length of initial contiguous sequence of GIDs
03439       //    on this process that are the same in both Maps; this
03440       //    equals the number of domain Map elements on this process)
03441       //
03442       // 2. permuteToLIDs and permuteFromLIDs (both empty in this
03443       //    case, since there's no permutation going on; the column
03444       //    Map starts with the domain Map's GIDs, and immediately
03445       //    after them come the remote GIDs)
03446       //
03447       // 3. remoteGIDs (exactly those GIDs that we found out above
03448       //    were not in the domain Map) and remoteLIDs (which we could
03449       //    have gotten above by using the three-argument version of
03450       //    getRemoteIndexList() that computes local indices as well
03451       //    as process ranks, instead of the two-argument version that
03452       //    was used above)
03453       //
03454       // 4. remotePIDs (which we have from the getRemoteIndexList()
03455       //    call above)
03456       //
03457       // 5. Sorting remotePIDs, and applying that permutation to
03458       //    remoteGIDs and remoteLIDs (by calling sort3 above instead
03459       //    of sort2)
03460       //
03461       // 6. Everything after the sort3 call in Import::setupExport():
03462       //    a. Create the Distributor via createFromRecvs(), which
03463       //       computes exportGIDs and exportPIDs
03464       //    b. Compute exportLIDs from exportGIDs (by asking the
03465       //       source Map, in this case the domain Map, to convert
03466       //       global to local)
03467       //
03468       // Steps 1-5 come for free, since we must do that work anyway in
03469       // order to compute the column Map.  In particular, Step 3 is
03470       // even more expensive than Step 6a, since it involves both
03471       // creating and using a new Distributor object.
03472 
03473     } // if the graph is globally indexed
03474 
03475     const global_size_t gstInv =
03476       Teuchos::OrdinalTraits<global_size_t>::invalid ();
03477     // FIXME (mfh 05 Mar 2014) Doesn't the index base of a Map have to
03478     // be the same as the Map's min GID? If the first column is empty
03479     // (contains no entries), then the column Map's min GID won't
03480     // necessarily be the same as the domain Map's index base.
03481     const GO indexBase = domainMap_->getIndexBase ();
03482     colMap_ = rcp (new map_type (gstInv, myColumns, indexBase,
03483                                  domainMap_->getComm (),
03484                                  domainMap_->getNode ()));
03485     checkInternalState ();
03486     return;
03487   }
03488 
03489 
03490   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
03491   void CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::mergeAllIndices()
03492   {
03493     TEUCHOS_TEST_FOR_EXCEPT( isGloballyIndexed() ); // call only after makeIndicesLocal()
03494     TEUCHOS_TEST_FOR_EXCEPT( ! isSorted() ); // call only after sortIndices()
03495     if (! isMerged ()) {
03496       for (size_t row=0; row < getNodeNumRows(); ++row) {
03497         RowInfo rowInfo = getRowInfo(row);
03498         mergeRowIndices(rowInfo);
03499       }
03500       // we just merged every row
03501       noRedundancies_ = true;
03502     }
03503   }
03504 
03505 
03506   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
03507   void
03508   CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::makeImportExport()
03509   {
03510     TEUCHOS_TEST_FOR_EXCEPTION(! hasColMap (), std::logic_error, "Tpetra::"
03511       "CrsGraph: It's not allowed to call makeImportExport() unless the graph "
03512       "has a column Map.");
03513     RCP<ParameterList> params = this->getNonconstParameterList (); // could be null
03514 
03515     // Don't do any checks to see if we need to create the Import, if
03516     // it exists already.
03517     //
03518     // FIXME (mfh 25 Mar 2013) This will become incorrect if we
03519     // change CrsGraph in the future to allow changing the column
03520     // Map after fillComplete.  For now, the column Map is fixed
03521     // after the first fillComplete call.
03522     if (importer_.is_null ()) {
03523       // Create the Import instance if necessary.
03524       if (domainMap_ != colMap_ && (! domainMap_->isSameAs (*colMap_))) {
03525         if (params.is_null () || ! params->isSublist ("Import")) {
03526           importer_ = rcp (new import_type (domainMap_, colMap_));
03527         } else {
03528           RCP<ParameterList> importSublist = sublist (params, "Import", true);
03529           importer_ = rcp (new import_type (domainMap_, colMap_, importSublist));
03530         }
03531       }
03532     }
03533 
03534     // Don't do any checks to see if we need to create the Export, if
03535     // it exists already.
03536     if (exporter_.is_null ()) {
03537       // Create the Export instance if necessary.
03538       if (rangeMap_ != rowMap_ && ! rangeMap_->isSameAs (*rowMap_)) {
03539         if (params.is_null () || ! params->isSublist ("Export")) {
03540           exporter_ = rcp (new export_type (rowMap_, rangeMap_));
03541         }
03542         else {
03543           RCP<ParameterList> exportSublist = sublist (params, "Export", true);
03544           exporter_ = rcp (new export_type (rowMap_, rangeMap_, exportSublist));
03545         }
03546       }
03547     }
03548   }
03549 
03550 
03551   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
03552   std::string
03553   CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::description() const
03554   {
03555     std::ostringstream oss;
03556     oss << DistObject<GlobalOrdinal,LocalOrdinal,GlobalOrdinal,Node>::description();
03557     if (isFillComplete()) {
03558       oss << "{status = fill complete"
03559           << ", global rows = " << getGlobalNumRows()
03560           << ", global cols = " << getGlobalNumCols()
03561           << ", global num entries = " << getGlobalNumEntries()
03562           << "}";
03563     }
03564     else {
03565       oss << "{status = fill not complete"
03566           << ", global rows = " << getGlobalNumRows()
03567           << "}";
03568     }
03569     return oss.str();
03570   }
03571 
03572 
03573   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
03574   void
03575   CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::
03576   describe (Teuchos::FancyOStream &out,
03577             const Teuchos::EVerbosityLevel verbLevel) const
03578   {
03579     const char tfecfFuncName[] = "describe()";
03580     using std::endl;
03581     using std::setw;
03582     using Teuchos::VERB_DEFAULT;
03583     using Teuchos::VERB_NONE;
03584     using Teuchos::VERB_LOW;
03585     using Teuchos::VERB_MEDIUM;
03586     using Teuchos::VERB_HIGH;
03587     using Teuchos::VERB_EXTREME;
03588     Teuchos::EVerbosityLevel vl = verbLevel;
03589     if (vl == VERB_DEFAULT) vl = VERB_LOW;
03590     RCP<const Comm<int> > comm = this->getComm();
03591     const int myImageID = comm->getRank(),
03592               numImages = comm->getSize();
03593     size_t width = 1;
03594     for (size_t dec=10; dec<getGlobalNumRows(); dec *= 10) {
03595       ++width;
03596     }
03597     width = std::max<size_t>(width,Teuchos::as<size_t>(11)) + 2;
03598     Teuchos::OSTab tab(out);
03599     //    none: print nothing
03600     //     low: print O(1) info from node 0
03601     //  medium: print O(P) info, num entries per node
03602     //    high: print O(N) info, num entries per row
03603     // extreme: print O(NNZ) info: print graph indices
03604     //
03605     // for medium and higher, print constituent objects at specified verbLevel
03606     if (vl != VERB_NONE) {
03607       if (myImageID == 0) out << this->description() << std::endl;
03608       // O(1) globals, minus what was already printed by description()
03609       if (isFillComplete() && myImageID == 0) {
03610         out << "Global number of diagonals = " << globalNumDiags_ << std::endl;
03611         out << "Global max number of row entries = " << globalMaxNumRowEntries_ << std::endl;
03612       }
03613       // constituent objects
03614       if (vl == VERB_MEDIUM || vl == VERB_HIGH || vl == VERB_EXTREME) {
03615         if (myImageID == 0) out << "\nRow map: " << std::endl;
03616         rowMap_->describe(out,vl);
03617         if (colMap_ != null) {
03618           if (myImageID == 0) out << "\nColumn map: " << std::endl;
03619           colMap_->describe(out,vl);
03620         }
03621         if (domainMap_ != null) {
03622           if (myImageID == 0) out << "\nDomain map: " << std::endl;
03623           domainMap_->describe(out,vl);
03624         }
03625         if (rangeMap_ != null) {
03626           if (myImageID == 0) out << "\nRange map: " << std::endl;
03627           rangeMap_->describe(out,vl);
03628         }
03629       }
03630       // O(P) data
03631       if (vl == VERB_MEDIUM || vl == VERB_HIGH || vl == VERB_EXTREME) {
03632         for (int imageCtr = 0; imageCtr < numImages; ++imageCtr) {
03633           if (myImageID == imageCtr) {
03634             out << "Node ID = " << imageCtr << std::endl
03635                 << "Node number of entries = " << nodeNumEntries_ << std::endl
03636                 << "Node number of diagonals = " << nodeNumDiags_ << std::endl
03637                 << "Node max number of entries = " << nodeMaxNumRowEntries_ << std::endl;
03638             if (isUpperTriangular()) {
03639               out << "Graph is upper triangular" << std::endl;
03640             }
03641             if (isLowerTriangular()) {
03642               out << "Graph is lower triangular" << std::endl;
03643             }
03644             if (indicesAreAllocated()) {
03645               out << "Node number of allocated entries = " << nodeNumAllocated_ << std::endl;
03646             }
03647             else {
03648               out << "Indices are not allocated." << std::endl;
03649             }
03650           }
03651           comm->barrier();
03652           comm->barrier();
03653           comm->barrier();
03654         }
03655       }
03656       // O(N) and O(NNZ) data
03657       if (vl == VERB_HIGH || vl == VERB_EXTREME) {
03658         TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(hasRowInfo() == false, std::runtime_error, ": reduce verbosity level; graph row information was deleted at fillComplete().");
03659         for (int imageCtr = 0; imageCtr < numImages; ++imageCtr) {
03660           if (myImageID == imageCtr) {
03661             out << std::setw(width) << "Node ID"
03662                 << std::setw(width) << "Global Row"
03663                 << std::setw(width) << "Num Entries";
03664             if (vl == VERB_EXTREME) {
03665               out << "Entries";
03666             }
03667             out << std::endl;
03668             for (size_t r=0; r < getNodeNumRows(); ++r) {
03669               RowInfo rowinfo = getRowInfo(r);
03670               GlobalOrdinal gid = rowMap_->getGlobalElement(r);
03671               out << std::setw(width) << myImageID
03672                   << std::setw(width) << gid
03673                   << std::setw(width) << rowinfo.numEntries;
03674               if (vl == VERB_EXTREME) {
03675                 if (isGloballyIndexed()) {
03676                   ArrayView<const GlobalOrdinal> rowview = getGlobalView(rowinfo);
03677                   for (size_t j=0; j < rowinfo.numEntries; ++j) out << rowview[j] << " ";
03678                 }
03679                 else if (isLocallyIndexed()) {
03680                   ArrayView<const LocalOrdinal> rowview = getLocalView(rowinfo);
03681                   for (size_t j=0; j < rowinfo.numEntries; ++j) out << colMap_->getGlobalElement(rowview[j]) << " ";
03682                 }
03683               }
03684               out << std::endl;
03685             }
03686           }
03687           comm->barrier();
03688           comm->barrier();
03689           comm->barrier();
03690         }
03691       }
03692     }
03693   }
03694 
03695 
03696   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
03697   bool
03698   CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::
03699   checkSizes (const SrcDistObject& source)
03700   {
03701     (void) source; // forestall "unused variable" compiler warnings
03702 
03703     // It's not clear what kind of compatibility checks on sizes can
03704     // be performed here.  Epetra_CrsGraph doesn't check any sizes for
03705     // compatibility.
03706     return true;
03707   }
03708 
03709 
03710   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
03711   void
03712   CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::
03713   copyAndPermute (const SrcDistObject& source,
03714                   size_t numSameIDs,
03715                   const Teuchos::ArrayView<const LocalOrdinal> &permuteToLIDs,
03716                   const Teuchos::ArrayView<const LocalOrdinal> &permuteFromLIDs)
03717   {
03718     using Teuchos::Array;
03719     using Teuchos::ArrayView;
03720     typedef LocalOrdinal LO;
03721     typedef GlobalOrdinal GO;
03722     const char tfecfFuncName[] = "copyAndPermute";
03723     typedef CrsGraph<LO, GO, Node, LocalMatOps> this_type;
03724     typedef RowGraph<LO, GO, Node> row_graph_type;
03725 
03726     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
03727       permuteToLIDs.size() != permuteFromLIDs.size(), std::runtime_error,
03728       ": permuteToLIDs and permuteFromLIDs must have the same size.");
03729     // Make sure that the source object has the right type.  We only
03730     // actually need it to be a RowGraph, with matching first three
03731     // template parameters.  If it's a CrsGraph, we can use view mode
03732     // instead of copy mode to get each row's data.
03733     //
03734     // FIXME (mfh 07 Jul 2013) It should not be necessary for any of
03735     // the template parameters but GO to match.  GO has to match
03736     // because the graph has to send indices as global ordinals, if
03737     // the source and target graphs do not have the same column Map.
03738     // If LO doesn't match, the graphs could communicate using global
03739     // indices.  It could be possible that Node affects the graph's
03740     // storage format, but packAndPrepare should assume a common
03741     // communication format in any case.
03742     const row_graph_type* srcRowGraph = dynamic_cast<const row_graph_type*> (&source);
03743     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
03744       srcRowGraph == NULL, std::invalid_argument,
03745       ": The source object must be a RowGraph with matching first three "
03746       "template parameters.");
03747 
03748     // If the source object is actually a CrsGraph, we can use view
03749     // mode instead of copy mode to access the entries in each row,
03750     // if the graph is not fill complete.
03751     const this_type* srcCrsGraph = dynamic_cast<const this_type*> (&source);
03752 
03753     const map_type& srcRowMap = * (srcRowGraph->getRowMap ());
03754     const map_type& tgtRowMap = * (this->getRowMap ());
03755     const bool src_filled = srcRowGraph->isFillComplete ();
03756     Array<GO> row_copy;
03757     LO myid = 0;
03758 
03759     //
03760     // "Copy" part of "copy and permute."
03761     //
03762     if (src_filled || srcCrsGraph == NULL) {
03763       // If the source graph is fill complete, we can't use view mode,
03764       // because the data might be stored in a different format not
03765       // compatible with the expectations of view mode.  Also, if the
03766       // source graph is not a CrsGraph, we can't use view mode,
03767       // because RowGraph only provides copy mode access to the data.
03768       for (size_t i = 0; i < numSameIDs; ++i, ++myid) {
03769         const GO gid = srcRowMap.getGlobalElement (myid);
03770         size_t row_length = srcRowGraph->getNumEntriesInGlobalRow (gid);
03771         row_copy.resize (row_length);
03772         size_t check_row_length = 0;
03773         srcRowGraph->getGlobalRowCopy (gid, row_copy (), check_row_length);
03774         this->insertGlobalIndices (gid, row_copy ());
03775       }
03776     } else {
03777       for (size_t i = 0; i < numSameIDs; ++i, ++myid) {
03778         const GO gid = srcRowMap.getGlobalElement (myid);
03779         ArrayView<const GO> row;
03780         srcCrsGraph->getGlobalRowView (gid, row);
03781         this->insertGlobalIndices (gid, row);
03782       }
03783     }
03784 
03785     //
03786     // "Permute" part of "copy and permute."
03787     //
03788     if (src_filled || srcCrsGraph == NULL) {
03789       for (LO i = 0; i < permuteToLIDs.size (); ++i) {
03790         const GO mygid = tgtRowMap.getGlobalElement (permuteToLIDs[i]);
03791         const GO srcgid = srcRowMap.getGlobalElement (permuteFromLIDs[i]);
03792         size_t row_length = srcRowGraph->getNumEntriesInGlobalRow (srcgid);
03793         row_copy.resize (row_length);
03794         size_t check_row_length = 0;
03795         srcRowGraph->getGlobalRowCopy (srcgid, row_copy (), check_row_length);
03796         this->insertGlobalIndices (mygid, row_copy ());
03797       }
03798     } else {
03799       for (LO i = 0; i < permuteToLIDs.size (); ++i) {
03800         const GO mygid = tgtRowMap.getGlobalElement (permuteToLIDs[i]);
03801         const GO srcgid = srcRowMap.getGlobalElement (permuteFromLIDs[i]);
03802         ArrayView<const GO> row;
03803         srcCrsGraph->getGlobalRowView (srcgid, row);
03804         this->insertGlobalIndices (mygid, row);
03805       }
03806     }
03807   }
03808 
03809 
03810   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
03811   void
03812   CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::
03813   packAndPrepare (const SrcDistObject& source,
03814                   const Teuchos::ArrayView<const LocalOrdinal> &exportLIDs,
03815                   Teuchos::Array<GlobalOrdinal> &exports,
03816                   const Teuchos::ArrayView<size_t> & numPacketsPerLID,
03817                   size_t& constantNumPackets,
03818                   Distributor &distor)
03819   {
03820     using Teuchos::Array;
03821     const char tfecfFuncName[] = "packAndPrepare";
03822     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
03823       exportLIDs.size() != numPacketsPerLID.size(), std::runtime_error,
03824       ": exportLIDs and numPacketsPerLID must have the same size.");
03825     typedef RowGraph<LocalOrdinal, GlobalOrdinal, Node> row_graph_type;
03826     const row_graph_type& srcGraph = dynamic_cast<const row_graph_type&> (source);
03827 
03828     // We don't check whether src_graph has had fillComplete called,
03829     // because it doesn't matter whether the *source* graph has been
03830     // fillComplete'd. The target graph can not be fillComplete'd yet.
03831     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
03832       this->isFillComplete (), std::runtime_error,
03833       ": The target graph of an Import or Export must not be fill complete.");
03834 
03835     srcGraph.pack (exportLIDs, exports, numPacketsPerLID, constantNumPackets, distor);
03836   }
03837 
03838 
03839   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
03840   void
03841   CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::
03842   pack (const Teuchos::ArrayView<const LocalOrdinal>& exportLIDs,
03843         Teuchos::Array<GlobalOrdinal>& exports,
03844         const Teuchos::ArrayView<size_t>& numPacketsPerLID,
03845         size_t& constantNumPackets,
03846         Distributor& distor) const
03847   {
03848     using Teuchos::Array;
03849     typedef LocalOrdinal LO;
03850     typedef GlobalOrdinal GO;
03851     const char tfecfFuncName[] = "pack";
03852     (void) distor; // forestall "unused argument" compiler warning
03853 
03854     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
03855       exportLIDs.size() != numPacketsPerLID.size(), std::runtime_error,
03856       ": exportLIDs and numPacketsPerLID must have the same size.");
03857 
03858     const map_type& srcMap = * (this->getMap ());
03859     constantNumPackets = 0;
03860 
03861     // Set numPacketsPerLID[i] to the number of entries owned by the
03862     // calling process in (local) row exportLIDs[i] of the graph, that
03863     // the caller wants us to send out.  Compute the total number of
03864     // packets (that is, entries) owned by this process in all the
03865     // rows that the caller wants us to send out.
03866     size_t totalNumPackets = 0;
03867     Array<GO> row;
03868     for (LO i = 0; i < exportLIDs.size (); ++i) {
03869       const GO GID = srcMap.getGlobalElement (exportLIDs[i]);
03870       size_t row_length = this->getNumEntriesInGlobalRow (GID);
03871       numPacketsPerLID[i] = row_length;
03872       totalNumPackets += row_length;
03873     }
03874 
03875     exports.resize (totalNumPackets);
03876 
03877     // Loop again over the rows to export, and pack rows of indices
03878     // into the output buffer.
03879     size_t exportsOffset = 0;
03880     for (LO i = 0; i < exportLIDs.size (); ++i) {
03881       const GO GID = srcMap.getGlobalElement (exportLIDs[i]);
03882       size_t row_length = this->getNumEntriesInGlobalRow (GID);
03883       row.resize (row_length);
03884       size_t check_row_length = 0;
03885       this->getGlobalRowCopy (GID, row (), check_row_length);
03886       typename Array<GO>::const_iterator row_iter = row.begin();
03887       typename Array<GO>::const_iterator row_end = row.end();
03888       size_t j = 0;
03889       for (; row_iter != row_end; ++row_iter, ++j) {
03890         exports[exportsOffset+j] = *row_iter;
03891       }
03892       exportsOffset += row.size();
03893     }
03894   }
03895 
03896 
03897   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
03898   void
03899   CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::
03900   unpackAndCombine (const Teuchos::ArrayView<const LocalOrdinal> &importLIDs,
03901                     const Teuchos::ArrayView<const GlobalOrdinal> &imports,
03902                     const Teuchos::ArrayView<size_t> &numPacketsPerLID,
03903                     size_t constantNumPackets,
03904                     Distributor& /* distor */,
03905                     CombineMode /* CM */)
03906   {
03907     // FIXME (mfh 02 Apr 2012) REPLACE combine mode has a perfectly
03908     // reasonable meaning, whether or not the matrix is fill complete.
03909     // It's just more work to implement.
03910 
03911     // We are not checking the value of the CombineMode input
03912     // argument.  For CrsGraph, we only support import/export
03913     // operations if fillComplete has not yet been called.  Any
03914     // incoming column-indices are inserted into the target graph. In
03915     // this context, CombineMode values of ADD vs INSERT are
03916     // equivalent. What is the meaning of REPLACE for CrsGraph? If a
03917     // duplicate column-index is inserted, it will be compressed out
03918     // when fillComplete is called.
03919     //
03920     // Note: I think REPLACE means that an existing row is replaced by
03921     // the imported row, i.e., the existing indices are cleared. CGB,
03922     // 6/17/2010
03923 
03924     const char tfecfFuncName[] = "unpackAndCombine";
03925     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
03926       importLIDs.size() != numPacketsPerLID.size(), std::runtime_error,
03927       ": importLIDs and numPacketsPerLID must have the same size.");
03928     TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
03929       isFillComplete (), std::runtime_error,
03930       ": Import or Export operations are not allowed on the destination "
03931       "CrsGraph if it is fill complete.");
03932     size_t importsOffset = 0;
03933 
03934     typedef typename ArrayView<const LocalOrdinal>::const_iterator iter_type;
03935     iter_type impLIDiter = importLIDs.begin();
03936     iter_type impLIDend = importLIDs.end();
03937 
03938     for (size_t i = 0; impLIDiter != impLIDend; ++impLIDiter, ++i) {
03939       LocalOrdinal row_length = numPacketsPerLID[i];
03940       ArrayView<const GlobalOrdinal> row (&imports[importsOffset], row_length);
03941       insertGlobalIndicesFiltered (this->getMap ()->getGlobalElement (*impLIDiter), row);
03942       importsOffset += row_length;
03943     }
03944   }
03945 
03946 
03947   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
03948   void
03949   CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::
03950   removeEmptyProcessesInPlace (const Teuchos::RCP<const Map<LocalOrdinal, GlobalOrdinal, Node> >& newMap)
03951   {
03952     using Teuchos::Comm;
03953     using Teuchos::null;
03954     using Teuchos::ParameterList;
03955     using Teuchos::RCP;
03956 
03957     // We'll set all the state "transactionally," so that this method
03958     // satisfies the strong exception guarantee.  This object's state
03959     // won't be modified until the end of this method.
03960     RCP<const map_type> rowMap, domainMap, rangeMap, colMap;
03961     RCP<import_type> importer;
03962     RCP<export_type> exporter;
03963 
03964     rowMap = newMap;
03965     RCP<const Comm<int> > newComm =
03966       (newMap.is_null ()) ? null : newMap->getComm ();
03967 
03968     if (! domainMap_.is_null ()) {
03969       if (domainMap_.getRawPtr () == rowMap_.getRawPtr ()) {
03970         // Common case: original domain and row Maps are identical.
03971         // In that case, we need only replace the original domain Map
03972         // with the new Map.  This ensures that the new domain and row
03973         // Maps _stay_ identical.
03974         domainMap = newMap;
03975       } else {
03976         domainMap = domainMap_->replaceCommWithSubset (newComm);
03977       }
03978     }
03979     if (! rangeMap_.is_null ()) {
03980       if (rangeMap_.getRawPtr () == rowMap_.getRawPtr ()) {
03981         // Common case: original range and row Maps are identical.  In
03982         // that case, we need only replace the original range Map with
03983         // the new Map.  This ensures that the new range and row Maps
03984         // _stay_ identical.
03985         rangeMap = newMap;
03986       } else {
03987         rangeMap = rangeMap_->replaceCommWithSubset (newComm);
03988       }
03989     }
03990     if (! colMap.is_null ()) {
03991       colMap = colMap_->replaceCommWithSubset (newComm);
03992     }
03993 
03994     // (Re)create the Export and / or Import if necessary.
03995     if (! newComm.is_null ()) {
03996       RCP<ParameterList> params = this->getNonconstParameterList (); // could be null
03997       //
03998       // The operations below are collective on the new communicator.
03999       //
04000       // (Re)create the Export object if necessary.  If I haven't
04001       // called fillComplete yet, I don't have a rangeMap, so I must
04002       // first check if the _original_ rangeMap is not null.  Ditto
04003       // for the Import object and the domain Map.
04004       if (! rangeMap_.is_null () &&
04005           rangeMap != rowMap &&
04006           ! rangeMap->isSameAs (*rowMap)) {
04007         if (params.is_null () || ! params->isSublist ("Export")) {
04008           exporter = rcp (new export_type (rowMap, rangeMap));
04009         }
04010         else {
04011           RCP<ParameterList> exportSublist = sublist (params, "Export", true);
04012           exporter = rcp (new export_type (rowMap, rangeMap, exportSublist));
04013         }
04014       }
04015       // (Re)create the Import object if necessary.
04016       if (! domainMap_.is_null () &&
04017           domainMap != colMap &&
04018           ! domainMap->isSameAs (*colMap)) {
04019         if (params.is_null () || ! params->isSublist ("Import")) {
04020           importer = rcp (new import_type (domainMap, colMap));
04021         } else {
04022           RCP<ParameterList> importSublist = sublist (params, "Import", true);
04023           importer = rcp (new import_type (domainMap, colMap, importSublist));
04024         }
04025       }
04026     } // if newComm is not null
04027 
04028     // Defer side effects until the end.  If no destructors throw
04029     // exceptions (they shouldn't anyway), then this method satisfies
04030     // the strong exception guarantee.
04031     exporter_ = exporter;
04032     importer_ = importer;
04033     rowMap_ = rowMap;
04034     // mfh 31 Mar 2013: DistObject's map_ is the row Map of a CrsGraph
04035     // or CrsMatrix.  CrsGraph keeps a redundant pointer (rowMap_) to
04036     // the same object.  We might want to get rid of this redundant
04037     // pointer sometime, but for now, we'll leave it alone and just
04038     // set map_ to the same object.
04039     this->map_ = rowMap;
04040     domainMap_ = domainMap;
04041     rangeMap_ = rangeMap;
04042     colMap_ = colMap;
04043   }
04044 
04045 
04046   template <class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
04047   bool
04048   CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::
04049   hasRowInfo () const
04050   {
04051 #ifdef HAVE_TPETRA_DEBUG
04052     bool actuallyHasRowInfo = true;
04053     if (indicesAreAllocated() && getProfileType() == StaticProfile && rowPtrs_ == null) {
04054       actuallyHasRowInfo = false;
04055     }
04056     TEUCHOS_TEST_FOR_EXCEPTION(
04057       actuallyHasRowInfo != haveRowInfo_,
04058       std::logic_error, "Internal logic error. Please contact Tpetra team.");
04059 #endif // HAVE_TPETRA_DEBUG
04060     return haveRowInfo_;
04061   }
04062 
04063 } // namespace Tpetra
04064 
04065 // Include KokkosRefactor partial specialisation if enabled
04066 #if defined(TPETRA_HAVE_KOKKOS_REFACTOR)
04067 #include "Tpetra_KokkosRefactor_CrsGraph_def.hpp"
04068 #endif
04069 
04070 
04071 //
04072 // Explicit instantiation macro
04073 //
04074 // Must be expanded from within the Tpetra namespace!
04075 //
04076 #define TPETRA_CRSGRAPH_GRAPH_INSTANT(LO,GO,NODE) template class CrsGraph< LO , GO , NODE >;
04077 #define TPETRA_CRSGRAPH_SORTROWINDICESANDVALUES_INSTANT(S,LO,GO,NODE) template void CrsGraph< LO , GO , NODE >::sortRowIndicesAndValues< S >(RowInfo, ArrayView< S > );
04078 #define TPETRA_CRSGRAPH_MERGEROWINDICESANDVALUES_INSTANT(S,LO,GO,NODE) template void CrsGraph< LO , GO , NODE >::mergeRowIndicesAndValues< S >(RowInfo, const ArrayView< S >& );
04079 #define TPETRA_CRSGRAPH_ALLOCATEVALUES1D_INSTANT(S,LO,GO,NODE) template ArrayRCP< S > CrsGraph< LO , GO , NODE >::allocateValues1D< S >() const;
04080 #define TPETRA_CRSGRAPH_ALLOCATEVALUES2D_INSTANT(S,LO,GO,NODE) template ArrayRCP< Array< S > > CrsGraph< LO , GO , NODE >::allocateValues2D< S >() const;
04081 
04082 #define TPETRA_CRSGRAPH_INSTANT(S,LO,GO,NODE)                    \
04083   TPETRA_CRSGRAPH_SORTROWINDICESANDVALUES_INSTANT(S,LO,GO,NODE)  \
04084   TPETRA_CRSGRAPH_MERGEROWINDICESANDVALUES_INSTANT(S,LO,GO,NODE) \
04085   TPETRA_CRSGRAPH_ALLOCATEVALUES1D_INSTANT(S,LO,GO,NODE)         \
04086   TPETRA_CRSGRAPH_ALLOCATEVALUES2D_INSTANT(S,LO,GO,NODE)
04087 
04088 #endif // TPETRA_CRSGRAPH_DEF_HPP
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines