Tpetra_CrsMatrix_def.hpp

00001 //@HEADER
00002 // ************************************************************************
00003 // 
00004 //               Tpetra: Templated Linear Algebra Services Package 
00005 //                 Copyright (2008) Sandia Corporation
00006 // 
00007 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00008 // license for use of this work by or on behalf of the U.S. Government.
00009 // 
00010 // This library is free software; you can redistribute it and/or modify
00011 // it under the terms of the GNU Lesser General Public License as
00012 // published by the Free Software Foundation; either version 2.1 of the
00013 // License, or (at your option) any later version.
00014 //  
00015 // This library is distributed in the hope that it will be useful, but
00016 // WITHOUT ANY WARRANTY; without even the implied warranty of
00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00018 // Lesser General Public License for more details.
00019 //  
00020 // You should have received a copy of the GNU Lesser General Public
00021 // License along with this library; if not, write to the Free Software
00022 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00023 // USA
00024 // Questions? Contact Michael A. Heroux (maherou@sandia.gov) 
00025 // 
00026 // ************************************************************************
00027 //@HEADER
00028 
00029 #ifndef TPETRA_CRSMATRIX_DEF_HPP
00030 #define TPETRA_CRSMATRIX_DEF_HPP
00031 
00032 // TODO: row-wise insertion of entries in globalAssemble() may be more efficient
00033 
00034 #include <Kokkos_NodeHelpers.hpp>
00035 #include <Kokkos_NodeTrace.hpp>
00036 
00037 #include <Teuchos_SerialDenseMatrix.hpp>
00038 #include <Teuchos_CommHelpers.hpp>
00039 #include <Teuchos_ScalarTraits.hpp>
00040 #include <Teuchos_OrdinalTraits.hpp>
00041 #include <Teuchos_TypeNameTraits.hpp>
00042 
00043 #include "Tpetra_CrsMatrixMultiplyOp.hpp" // must include for implicit instantiation to work
00044 #ifdef DOXYGEN_USE_ONLY
00045   #include "Tpetra_CrsMatrix_decl.hpp"
00046 #endif
00047 
00048 #ifndef DOXYGEN_SHOULD_SKIP_THIS
00049 
00050 template <class Ordinal, class Scalar>
00051 bool std::operator<(const Tpetra::CrsIJV<Ordinal,Scalar> &ijv1, const Tpetra::CrsIJV<Ordinal,Scalar> &ijv2) {
00052   return ijv1.i < ijv2.i;
00053 }
00054 #endif
00055 
00056 namespace Tpetra {
00057 
00058   template <class Ordinal, class Scalar>
00059   CrsIJV<Ordinal,Scalar>::CrsIJV() {}
00060 
00061   template <class Ordinal, class Scalar>
00062   CrsIJV<Ordinal,Scalar>::CrsIJV(Ordinal row, Ordinal col, const Scalar &val) {
00063     i = row;
00064     j = col;
00065     v = val;
00066   }
00067 
00070 
00071   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00072   CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::CrsMatrix(
00073                                           const Teuchos::RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > &rowMap, 
00074                                           size_t maxNumEntriesPerRow, 
00075                                           ProfileType pftype)
00076   : DistObject<char, LocalOrdinal,GlobalOrdinal,Node>(rowMap)
00077   , lclMatrix_(rowMap->getNodeNumElements(), rowMap->getNode())
00078   , lclMatOps_(rowMap->getNode())
00079   , constructedWithOptimizedGraph_(false)
00080   , fillComplete_(false)
00081   {
00082     try {
00083       myGraph_ = Teuchos::rcp( new CrsGraph<LocalOrdinal,GlobalOrdinal,Node>(rowMap,maxNumEntriesPerRow,pftype) );
00084     }
00085     catch (std::exception &e) {
00086       TEST_FOR_EXCEPTION(true, std::runtime_error,
00087           Teuchos::typeName(*this) << "::CrsMatrix(): caught exception while allocating CrsGraph object: " 
00088           << std::endl << e.what() << std::endl);
00089     }
00090     // it is okay to create this now; this will prevent us from having to check for it on every call to apply()
00091     // we will use a non-owning rcp to wrap *this; this is safe as long as we do not shared sameScalarMultiplyOp_ with anyone, 
00092     // which would allow it to persist past the destruction of *this
00093     sameScalarMultiplyOp_ = createCrsMatrixMultiplyOp<Scalar>( rcp(this,false).getConst() );
00094 
00095     checkInternalState();
00096   }
00097 
00098   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00099   CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::CrsMatrix(
00100                                           const Teuchos::RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > &rowMap, 
00101                                           const Teuchos::ArrayRCP<const size_t> &NumEntriesPerRowToAlloc, 
00102                                           ProfileType pftype)
00103   : DistObject<char, LocalOrdinal,GlobalOrdinal,Node>(rowMap)
00104   , lclMatrix_(rowMap->getNodeNumElements(), rowMap->getNode())
00105   , lclMatOps_(rowMap->getNode())
00106   , constructedWithOptimizedGraph_(false)
00107   , fillComplete_(false)
00108   {
00109     try {
00110       myGraph_ = Teuchos::rcp( new CrsGraph<LocalOrdinal,GlobalOrdinal,Node>(rowMap,NumEntriesPerRowToAlloc,pftype) );
00111     }
00112     catch (std::exception &e) {
00113       TEST_FOR_EXCEPTION(true, std::runtime_error,
00114           Teuchos::typeName(*this) << "::CrsMatrix(): caught exception while allocating CrsGraph object: " 
00115           << std::endl << e.what() << std::endl);
00116     }
00117     // it is okay to create this now; this will prevent us from having to check for it on every call to apply()
00118     // we will use a non-owning rcp to wrap *this; this is safe as long as we do not shared sameScalarMultiplyOp_ with anyone, 
00119     // which would allow it to persist past the destruction of *this
00120     sameScalarMultiplyOp_ = createCrsMatrixMultiplyOp<Scalar>( rcp(this,false).getConst() );
00121 
00122     checkInternalState();
00123   }
00124 
00125   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00126   CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::CrsMatrix(
00127                                           const Teuchos::RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > &rowMap, 
00128                                           const Teuchos::RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > &colMap, 
00129                                           size_t maxNumEntriesPerRow, 
00130                                           ProfileType pftype)
00131   : DistObject<char, LocalOrdinal,GlobalOrdinal,Node>(rowMap)
00132   , lclMatrix_(rowMap->getNodeNumElements(), rowMap->getNode())
00133   , lclMatOps_(rowMap->getNode())
00134   , constructedWithOptimizedGraph_(false)
00135   , fillComplete_(false)
00136   {
00137     try {
00138       myGraph_ = Teuchos::rcp( new CrsGraph<LocalOrdinal,GlobalOrdinal,Node>(rowMap,colMap,maxNumEntriesPerRow,pftype) );
00139     }
00140     catch (std::exception &e) {
00141       TEST_FOR_EXCEPTION(true, std::runtime_error,
00142           Teuchos::typeName(*this) << "::CrsMatrix(): caught exception while allocating CrsGraph object: " 
00143           << std::endl << e.what() << std::endl);
00144     }
00145     // it is okay to create this now; this will prevent us from having to check for it on every call to apply()
00146     // we will use a non-owning rcp to wrap *this; this is safe as long as we do not shared sameScalarMultiplyOp_ with anyone, 
00147     // which would allow it to persist past the destruction of *this
00148     sameScalarMultiplyOp_ = createCrsMatrixMultiplyOp<Scalar>( rcp(this,false).getConst() );
00149 
00150     checkInternalState();
00151   }
00152 
00153   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00154   CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::CrsMatrix(
00155                                           const Teuchos::RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > &rowMap, 
00156                                           const Teuchos::RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > &colMap, 
00157                                           const Teuchos::ArrayRCP<const size_t> &NumEntriesPerRowToAlloc, 
00158                                           ProfileType pftype)
00159   : DistObject<char, LocalOrdinal,GlobalOrdinal,Node>(rowMap)
00160   , lclMatrix_(rowMap->getNodeNumElements(), rowMap->getNode())
00161   , lclMatOps_(rowMap->getNode())
00162   , constructedWithOptimizedGraph_(false)
00163   , fillComplete_(false)
00164   {
00165     try {
00166       myGraph_ = Teuchos::rcp( new CrsGraph<LocalOrdinal,GlobalOrdinal,Node>(rowMap,colMap,NumEntriesPerRowToAlloc,pftype) );
00167     }
00168     catch (std::exception &e) {
00169       TEST_FOR_EXCEPTION(true, std::runtime_error,
00170           Teuchos::typeName(*this) << "::CrsMatrix(): caught exception while allocating CrsGraph object: " 
00171           << std::endl << e.what() << std::endl);
00172     }
00173     // it is okay to create this now; this will prevent us from having to check for it on every call to apply()
00174     // we will use a non-owning rcp to wrap *this; this is safe as long as we do not shared sameScalarMultiplyOp_ with anyone, 
00175     // which would allow it to persist past the destruction of *this
00176     sameScalarMultiplyOp_ = createCrsMatrixMultiplyOp<Scalar>( rcp(this,false).getConst() );
00177 
00178     checkInternalState();
00179   }
00180 
00181   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00182   CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::CrsMatrix(const Teuchos::RCP<const CrsGraph<LocalOrdinal,GlobalOrdinal,Node> > &graph)
00183   : DistObject<char, LocalOrdinal,GlobalOrdinal,Node>(graph->getRowMap())
00184   , staticGraph_(graph)
00185   , lclMatrix_(graph->getRowMap()->getNodeNumElements(), graph->getRowMap()->getNode())
00186   , lclMatOps_(graph->getNode())
00187   , fillComplete_(false) 
00188   {
00189     TEST_FOR_EXCEPTION(staticGraph_ == Teuchos::null, std::runtime_error,
00190         Teuchos::typeName(*this) << "::CrsMatrix(graph): specified pointer is null.");
00191     // we prohibit the case where the graph is not yet filled
00192     // we do not require optimized storage, but check to ensure that optimizeStorage() isn't called between now and later.
00193     TEST_FOR_EXCEPTION( staticGraph_->isFillComplete() == false, std::runtime_error, 
00194         Teuchos::typeName(*this) << "::CrsMatrix(graph): specified graph is not fill-complete. You must fillComplete() the graph before using it to construct a CrsMatrix.");
00195     constructedWithOptimizedGraph_ = staticGraph_->isStorageOptimized();
00196 
00197     // it is okay to create this now; this will prevent us from having to check for it on every call to apply()
00198     // we will use a non-owning rcp to wrap *this; this is safe as long as we do not shared sameScalarMultiplyOp_ with anyone, 
00199     // which would allow it to persist past the destruction of *this
00200     sameScalarMultiplyOp_ = createCrsMatrixMultiplyOp<Scalar>( rcp(this,false).getConst() );
00201 
00202     // the graph has entries, and the matrix should have entries as well, set to zero. no need or point in lazy allocating in this case.
00203     // first argument doesn't actually matter, due to the third
00204     allocateValues( CrsGraph<LocalOrdinal,GlobalOrdinal,Node>::AllocateLocal, Teuchos::ScalarTraits<Scalar>::zero(), GraphAlreadyAllocated );
00205 
00206     checkInternalState();
00207   }
00208 
00209 
00210   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00211   CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::~CrsMatrix() {
00212   }
00213 
00214 
00215   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00216   const Teuchos::RCP<const Teuchos::Comm<int> > &
00217   CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getComm() const {
00218     return getCrsGraph()->getComm();
00219   }
00220 
00221 
00222   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00223   Teuchos::RCP<Node>
00224   CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getNode() const {
00225     return lclMatOps_.getNode();
00226   }
00227 
00228 
00229   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00230   ProfileType CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getProfileType() const {
00231     return getCrsGraph()->getProfileType();
00232   }
00233 
00234 
00235   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00236   bool CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::isFillComplete() const {
00237     return fillComplete_; 
00238   }
00239 
00240 
00241   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00242   bool CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::isFillResumed() const {
00243     return !fillComplete_; 
00244   }
00245 
00246 
00247   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00248   bool CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::isStorageOptimized() const {
00249     return getCrsGraph()->isStorageOptimized();
00250   }
00251 
00252 
00253   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00254   bool CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::isLocallyIndexed() const {
00255     return getCrsGraph()->isLocallyIndexed();
00256   }
00257 
00258 
00259   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00260   bool CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::isGloballyIndexed() const {
00261     return getCrsGraph()->isGloballyIndexed();
00262   }
00263 
00264 
00265   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00266   bool CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::hasColMap() const {
00267     return getCrsGraph()->hasColMap();
00268   }
00269 
00270 
00271   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00272   global_size_t CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getGlobalNumEntries() const {
00273     return getCrsGraph()->getGlobalNumEntries();
00274   }
00275 
00276 
00277   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00278   size_t CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getNodeNumEntries() const {
00279     return getCrsGraph()->getNodeNumEntries();
00280   }
00281 
00282 
00283   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00284   global_size_t CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getGlobalNumRows() const {
00285     return getCrsGraph()->getGlobalNumRows(); 
00286   }
00287 
00288 
00289   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00290   global_size_t CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getGlobalNumCols() const { 
00291     return getCrsGraph()->getGlobalNumCols(); 
00292   }
00293 
00294 
00295   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00296   size_t CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getNodeNumRows() const { 
00297     return getCrsGraph()->getNodeNumRows(); 
00298   }
00299 
00300 
00301   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00302   size_t CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getNodeNumCols() const { 
00303     return getCrsGraph()->getNodeNumCols(); 
00304   }
00305 
00306 
00307   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00308   global_size_t CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getGlobalNumDiags() const { 
00309     return getCrsGraph()->getGlobalNumDiags(); 
00310   }
00311 
00312 
00313   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00314   size_t CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getNodeNumDiags() const { 
00315     return getCrsGraph()->getNodeNumDiags(); 
00316   }
00317 
00318 
00319   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00320   size_t CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getNumEntriesInGlobalRow(GlobalOrdinal globalRow) const { 
00321     return getCrsGraph()->getNumEntriesInGlobalRow(globalRow); 
00322   }
00323 
00324 
00325   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00326   size_t CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getNumEntriesInLocalRow(LocalOrdinal localRow) const { 
00327     return getCrsGraph()->getNumEntriesInLocalRow(localRow);
00328   }
00329 
00330 
00331   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00332   size_t CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getGlobalMaxNumRowEntries() const { 
00333     return getCrsGraph()->getGlobalMaxNumRowEntries(); 
00334   }
00335 
00336 
00337   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00338   size_t CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getNodeMaxNumRowEntries() const { 
00339     return getCrsGraph()->getNodeMaxNumRowEntries(); 
00340   }
00341 
00342 
00343   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00344   GlobalOrdinal CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getIndexBase() const { 
00345     return getRowMap()->getIndexBase(); 
00346   }
00347 
00348 
00349   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00350   const Teuchos::RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > & 
00351   CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getRowMap() const { 
00352     return getCrsGraph()->getRowMap(); 
00353   }
00354 
00355 
00356   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00357   const Teuchos::RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > & 
00358   CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getColMap() const {
00359     return getCrsGraph()->getColMap(); 
00360   }
00361 
00362 
00363   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00364   const Teuchos::RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > & 
00365   CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getDomainMap() const { 
00366     return getCrsGraph()->getDomainMap(); 
00367   }
00368 
00369 
00370   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00371   const Teuchos::RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > & 
00372   CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getRangeMap() const { 
00373     return getCrsGraph()->getRangeMap(); 
00374   }
00375 
00376 
00377   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00378   Teuchos::RCP<const RowGraph<LocalOrdinal,GlobalOrdinal,Node> >
00379   CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getGraph() const { 
00380     if (staticGraph_ != Teuchos::null) return staticGraph_;
00381     return myGraph_;
00382   }
00383 
00384 
00385   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00386   Teuchos::RCP<const CrsGraph<LocalOrdinal,GlobalOrdinal,Node> >
00387   CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getCrsGraph() const { 
00388     if (staticGraph_ != Teuchos::null) return staticGraph_;
00389     return myGraph_;
00390   }
00391 
00392 
00393   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00394   bool CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::isLowerTriangular() const { 
00395     return getCrsGraph()->isLowerTriangular(); 
00396   }
00397 
00398 
00399   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00400   bool CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::isUpperTriangular() const { 
00401     return getCrsGraph()->isUpperTriangular(); 
00402   }
00403 
00404 
00405   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00406   bool CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::isStaticGraph() const { 
00407     return (staticGraph_ != Teuchos::null);
00408   }
00409 
00410 
00411   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00412   bool CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::hasTransposeApply() const {
00413     return true;
00414   }
00415 
00416 
00419   //                                                                         //
00420   //                    Internal utility methods                             //
00421   //                                                                         //
00424 
00425 
00428   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00429   void CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::allocateValues(typename CrsGraph<LocalOrdinal,GlobalOrdinal,Node>::AllocateLocalGlobal lg, Scalar alpha, GraphAllocationStatus gas) {
00430     Teuchos::RCP< const CrsGraph<LocalOrdinal,GlobalOrdinal,Node> > graph = getCrsGraph();
00431     // allocate values and, optionally, ask graph to allocate indices
00432 #ifdef HAVE_TPETRA_DEBUG
00433     // if the graph is already allocated, then gas should be GraphAlreadyAllocated
00434     // otherwise, gas should be GraphNotYetAllocated
00435     TEST_FOR_EXCEPTION((gas == GraphAlreadyAllocated) != graph->indicesAreAllocated(), std::logic_error,
00436         Teuchos::typeName(*this) << "::insertLocalValues(): Internal logic error. Please contact Tpetra team.");
00437     // if the graph is unallocated, then it better be a matrix-owned graph
00438     TEST_FOR_EXCEPTION(graph->indicesAreAllocated() == false && myGraph_ == Teuchos::null, std::logic_error,
00439         Teuchos::typeName(*this) << "::insertLocalValues(): Internal logic error. Please contact Tpetra team.");
00440 #endif
00441     if (gas == GraphNotYetAllocated) {
00442       myGraph_->allocateIndices(lg);
00443     }
00444     const size_t nlrs = getRowMap()->getNodeNumElements(),
00445                   nta = graph->getNodeAllocationSize(),
00446            numEntries = graph->getNodeNumEntries();
00447     // do we even have anything to allocate?
00448     if (nta > 0) {
00449       Teuchos::RCP<Node> node = lclMatrix_.getNode();
00451       if (getProfileType() == StaticProfile) {
00452         //
00453         //  STATIC ALLOCATION PROFILE
00454         //
00455         // determine how many entries to allocate and setup offsets into 1D arrays
00456         values1D_ = node->template allocBuffer<Scalar>(nta);
00457         // init values if the graph already has valid entries
00458         if (numEntries > 0) {
00459           Kokkos::ReadyBufferHelper<Node> rbh(node);
00460           Kokkos::InitOp<Scalar> wdp;
00461           wdp.alpha = alpha;
00462           rbh.begin();
00463           wdp.x   = rbh.addNonConstBuffer(values1D_);
00464           rbh.end();
00465           node->template parallel_for<Kokkos::InitOp<Scalar> >(0,nta,wdp);
00466           KOKKOS_NODE_TRACE("CrsMatrix::allocateValues()")
00467         }
00468         else {
00469           KOKKOS_NODE_TRACE("CrsMatrix::allocateValues()")
00470           // a simple WriteOnly view will suffice
00471         }
00472       }
00473       else {
00474         //
00475         //  DYNAMIC ALLOCATION PROFILE
00476         //
00477         Kokkos::InitOp<Scalar> wdp;
00478         wdp.alpha = alpha;
00479         // allocate array of buffers
00480         values2D_ = Teuchos::arcp< Teuchos::ArrayRCP<Scalar> >(nlrs);
00481         bool someRowWasInitialized = false;
00482         Kokkos::ReadyBufferHelper<Node> rbh(node);
00483         for (size_t r=0; r<nlrs; ++r) {
00484           // this call to getNumAllocatedEntries() is cheap for the DynamicProfile case
00485           const size_t ntarow = graph->getNumAllocatedEntriesInLocalRow(r),
00486                         nErow = graph->getNumEntriesInLocalRow(r);
00487           if (ntarow > 0) {
00488             // allocate values for this row
00489             values2D_[r] = node->template allocBuffer<Scalar>(ntarow);
00490             // initi values in parallel, if the graph already has valid entries
00491             if (nErow > 0) {
00492               rbh.begin();
00493               wdp.x   = rbh.addNonConstBuffer(values2D_[r]);
00494               rbh.end();
00495               node->template parallel_for<Kokkos::InitOp<Scalar> >(0,ntarow,wdp);
00496               someRowWasInitialized = true;
00497             }
00498           }
00499         }
00500       }
00501     } // num to allocate > 0
00502   }
00503 
00504 
00507   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00508   void CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::checkInternalState() const {
00509 #ifdef HAVE_TPETRA_DEBUG
00510     Teuchos::RCP<Node> node = getNode();
00511     Teuchos::RCP< const CrsGraph<LocalOrdinal,GlobalOrdinal,Node> > graph = getCrsGraph();
00512     using Teuchos::null;
00513     std::string err = Teuchos::typeName(*this) + "::checkInternalState(): Likely internal logic error. Please contact Tpetra team.";
00514     // check the internal state of this data structure
00515     // this is called by numerous state-changing methods, in a debug build, to ensure that the object 
00516     // always remains in a valid state
00517 
00518     // we must have a single graph, static or matrix-owned
00519     TEST_FOR_EXCEPTION( staticGraph_ == null && myGraph_ == null,                                      std::logic_error, err );
00520     TEST_FOR_EXCEPTION( staticGraph_ != null && myGraph_ != null,                                      std::logic_error, err );
00521     // constructedWithOptimizedGraph_ should only be true if matrix was constructed with a graph, in which case staticGraph_ should be non-null
00522     TEST_FOR_EXCEPTION( constructedWithOptimizedGraph_ == true && staticGraph_ == null,                std::logic_error, err ); 
00523     // if matrix is fill complete, then graph must be fill complete
00524     TEST_FOR_EXCEPTION( fillComplete_ == true && graph->isFillComplete() == false,             std::logic_error, err );
00525     // if matrix is storage optimized, it should have a 1D allocation 
00526     TEST_FOR_EXCEPTION( isStorageOptimized() == true && values2D_ != Teuchos::null,               std::logic_error, err );
00527     // if matrix/graph are static profile, then 2D allocation should not be present
00528     TEST_FOR_EXCEPTION( graph->getProfileType() == StaticProfile  && values2D_ != Teuchos::null, std::logic_error, err );
00529     // if matrix/graph are dynamic profile, then 1D allocation should not be present
00530     TEST_FOR_EXCEPTION( graph->getProfileType() == DynamicProfile && values1D_ != Teuchos::null, std::logic_error, err );
00531     // if values are allocated and they are non-zero in number, then one of the allocations should be present
00532     TEST_FOR_EXCEPTION( graph->indicesAreAllocated() && graph->getNodeAllocationSize() > 0 && values2D_ == Teuchos::null && values1D_ == Teuchos::null,
00533                         std::logic_error, err );
00534     // we can nae have both a 1D and 2D allocation
00535     TEST_FOR_EXCEPTION( values1D_ != Teuchos::null && values2D_ != Teuchos::null, std::logic_error, err );
00536     // compare matrix allocations against graph allocations
00537     if (graph->indicesAreAllocated() && graph->getNodeAllocationSize() > 0) {
00538       if (graph->getProfileType() == StaticProfile) {
00539         if (graph->isLocallyIndexed()) {
00540           TEST_FOR_EXCEPTION( values1D_.size() != graph->pbuf_lclInds1D_.size(),  std::logic_error, err );
00541         }
00542         else {
00543           TEST_FOR_EXCEPTION( values1D_.size() != graph->pbuf_gblInds1D_.size(),  std::logic_error, err );
00544         }
00545       } 
00546       else { // graph->getProfileType() == DynamicProfile
00547         if (graph->isLocallyIndexed()) {
00548           for (size_t r=0; r < getNodeNumRows(); ++r) {
00549             TEST_FOR_EXCEPTION( (values2D_[r] == Teuchos::null) != (graph->pbuf_lclInds2D_[r] == Teuchos::null), std::logic_error, err );
00550             if (values2D_[r] != Teuchos::null && graph->pbuf_lclInds2D_[r] != Teuchos::null) {
00551               TEST_FOR_EXCEPTION( values2D_[r].size() != graph->pbuf_lclInds2D_[r].size(), std::logic_error, err );
00552             }
00553           }
00554         }
00555         else {
00556           for (size_t r=0; r < getNodeNumRows(); ++r) {
00557             TEST_FOR_EXCEPTION( (values2D_[r] == Teuchos::null) != (graph->pbuf_gblInds2D_[r] == Teuchos::null), std::logic_error, err );
00558             if (values2D_[r] != Teuchos::null && graph->pbuf_gblInds2D_[r] != Teuchos::null) {
00559               TEST_FOR_EXCEPTION( values2D_[r].size() != graph->pbuf_gblInds2D_[r].size(), std::logic_error, err );
00560             }
00561           }
00562         }
00563       }
00564     }
00565 #endif
00566   }
00567 
00568 
00571   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00572   Teuchos::ArrayRCP<const Scalar> 
00573   CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getFullView(size_t myRow, RowInfo sizeInfo) const {
00574 #ifdef HAVE_TPETRA_DEBUG
00575     std::string err = Teuchos::typeName(*this) + "::getFullView(): Internal logic error. Please contact Tpetra team.";
00576     TEST_FOR_EXCEPTION(getRowMap()->isNodeLocalElement(myRow) == false, std::logic_error, err);
00577 #endif
00578     Teuchos::ArrayRCP<const Scalar> values = Teuchos::null;
00579     // sizeInfo indicates the allocation size for this row, whether it has actually been allocated or not
00580     if (sizeInfo.allocSize > 0 && getCrsGraph()->indicesAreAllocated()) {
00581       Teuchos::RCP<Node> node = getNode();
00582       if (getCrsGraph()->getProfileType() == StaticProfile) {
00583         KOKKOS_NODE_TRACE("CrsMatrix::getFullView()")
00584         values = node->template viewBuffer<Scalar>(sizeInfo.allocSize, values1D_ + sizeInfo.offset1D);
00585       }
00586       else {  // dynamic profile
00587         KOKKOS_NODE_TRACE("CrsMatrix::getFullView()")
00588         values = node->template viewBuffer<Scalar>(sizeInfo.allocSize, values2D_[myRow]);
00589       }
00590     }
00591 #ifdef HAVE_TPETRA_DEBUG
00592     TEST_FOR_EXCEPTION(getCrsGraph()->indicesAreAllocated() && getCrsGraph()->getNodeAllocationSize() > 0 && values != Teuchos::null && static_cast<size_t>(values.size()) != sizeInfo.allocSize, std::logic_error, err);
00593 #endif
00594     return values;
00595   }
00596 
00597 
00600   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00601   Teuchos::ArrayRCP<Scalar> 
00602   CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getFullViewNonConst(size_t myRow, RowInfo sizeInfo) {
00603 #ifdef HAVE_TPETRA_DEBUG
00604     std::string err = Teuchos::typeName(*this) + "::getFullViewNonConst(): Internal logic error. Please contact Tpetra team.";
00605     TEST_FOR_EXCEPTION(getRowMap()->isNodeLocalElement(myRow) == false, std::logic_error, err);
00606 #endif
00607     Teuchos::ArrayRCP<Scalar> values = Teuchos::null;
00608     // sizeInfo indicates the allocation size for this row, whether it has actually been allocated or not
00609     if (sizeInfo.allocSize > 0 && getCrsGraph()->indicesAreAllocated()) {
00610       Teuchos::RCP<Node> node = getNode();
00611       // if there are no valid entries, then this view can be constructed WriteOnly
00612       Kokkos::ReadWriteOption rw = (sizeInfo.numEntries == 0 ? Kokkos::WriteOnly : Kokkos::ReadWrite);
00613       if (getCrsGraph()->getProfileType() == StaticProfile) {
00614         KOKKOS_NODE_TRACE("CrsMatrix::getFullViewNonConst()")
00615         values = node->template viewBufferNonConst<Scalar>(rw, sizeInfo.allocSize, values1D_ + sizeInfo.offset1D);
00616       }
00617       else {  // dynamic profile
00618         KOKKOS_NODE_TRACE("CrsMatrix::getFullViewNonConst()")
00619         values = node->template viewBufferNonConst<Scalar>(rw, sizeInfo.allocSize, values2D_[myRow]);
00620       }
00621     }
00622 #ifdef HAVE_TPETRA_DEBUG
00623     TEST_FOR_EXCEPTION(getCrsGraph()->indicesAreAllocated() && getCrsGraph()->getNodeAllocationSize() > 0 && values != Teuchos::null && static_cast<size_t>(values.size()) != sizeInfo.allocSize, std::logic_error, err);
00624 #endif
00625     return values;
00626   }
00627 
00628 
00631   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00632   void CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::updateAllocation(size_t lrow, size_t allocSize) {
00633     using Teuchos::ArrayRCP;
00634     Teuchos::RCP<Node> node = getNode();
00635     RowInfo sizeInfo = getCrsGraph()->getRowInfo(lrow);
00636     // allocate a larger space for row "lrow"
00637     // copy any existing data from previous allocation to new allocation
00638     // update sizes
00639     // 
00640     // if we already have views of the data, we will create a new view and do the copy on the host
00641     // otherwise, don't create a view, and do the copy on the device
00642     // 
00643     if (values2D_ == Teuchos::null) {
00644       values2D_ = Teuchos::arcp< ArrayRCP<Scalar> >(getNodeNumRows());
00645     }
00646 #ifdef HAVE_TPETRA_DEBUG
00647     TEST_FOR_EXCEPT( getRowMap()->isNodeLocalElement(lrow) == false );
00648     TEST_FOR_EXCEPT( values2D_[lrow] != Teuchos::null && allocSize < static_cast<size_t>(values2D_[lrow].size()) );
00649     TEST_FOR_EXCEPT( allocSize == 0 );
00650 #endif
00651     ArrayRCP<Scalar> old_alloc, new_row;
00652     old_alloc = values2D_[lrow];
00653     values2D_[lrow] = node->template allocBuffer<Scalar>(allocSize);
00654     if (sizeInfo.numEntries) {
00655       KOKKOS_NODE_TRACE("CrsMatrix::updateAllocation()")
00656       node->template copyBuffers<Scalar>(sizeInfo.numEntries,old_alloc,values2D_[lrow]);
00657     }
00658     old_alloc = Teuchos::null;
00659   }
00660 
00661 
00664   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00665   void CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::fillLocalMatrix() {
00666     Teuchos::RCP< const CrsGraph<LocalOrdinal,GlobalOrdinal,Node> > graph = getCrsGraph();
00667     lclMatrix_.clear();
00668     if (isStorageOptimized()) {
00669       // fill packed matrix; it is okay for values1D_ to be null; the matrix will flag itself as empty
00670       lclMatrix_.setPackedValues(values1D_);
00671     }
00672     else if (graph->getProfileType() == StaticProfile) {
00673       if (values1D_ != Teuchos::null) {
00674         const size_t nlrs = getNodeNumRows();
00675         for (size_t r=0; r < nlrs; ++r) {
00676           RowInfo sizeInfo = graph->getRowInfo(r);
00677           Teuchos::ArrayRCP<const Scalar> rowvals;
00678           if (sizeInfo.numEntries > 0) {
00679             rowvals = values1D_.persistingView(sizeInfo.offset1D, sizeInfo.numEntries);
00680             lclMatrix_.set2DValues(r,rowvals);
00681           }
00682         }
00683       }
00684     }
00685     else if (graph->getProfileType() == DynamicProfile) {
00686       if (values2D_ != Teuchos::null) {
00687         const size_t nlrs = getNodeNumRows();
00688         for (size_t r=0; r < nlrs; ++r) {
00689           RowInfo sizeInfo = graph->getRowInfo(r);
00690           Teuchos::ArrayRCP<const Scalar> rowvals = values2D_[r];
00691           if (sizeInfo.numEntries > 0) {
00692             rowvals = rowvals.persistingView(0,sizeInfo.numEntries);
00693             lclMatrix_.set2DValues(r,rowvals);
00694           }
00695         }
00696       }
00697     }
00698 
00699     // submit local matrix and local graph to lclMatOps_
00700     // this is permitted to view, but we don't care whether they do or not
00701     lclMatOps_.clear();
00702     Teuchos::DataAccess ret;
00703     ret = lclMatOps_.initializeStructure( graph->lclGraph_, Teuchos::View );
00704     ret = lclMatOps_.initializeValues( lclMatrix_, Teuchos::View );
00705   }
00706 
00707 
00710   //                                                                         //
00711   //                  User-visible class methods                             //
00712   //                                                                         //
00715 
00716 
00719   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00720   void CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::insertLocalValues(
00721                                                          LocalOrdinal localRow, 
00722                          const Teuchos::ArrayView<const LocalOrdinal> &indices,
00723                          const Teuchos::ArrayView<const Scalar>       &values) {
00724     using Teuchos::ArrayRCP;
00725     TEST_FOR_EXCEPTION(isStaticGraph() == true, std::runtime_error,
00726         Teuchos::typeName(*this) << "::insertLocalValues(): matrix was constructed with static graph; cannot insert new entries.");
00727     TEST_FOR_EXCEPTION(isStorageOptimized() == true, std::runtime_error,
00728         Teuchos::typeName(*this) << "::insertLocalValues(): cannot insert new values after optimizeStorage() has been called.");
00729     TEST_FOR_EXCEPTION(myGraph_->isGloballyIndexed() == true, std::runtime_error,
00730         Teuchos::typeName(*this) << "::insertLocalValues(): graph indices are global; use insertGlobalValues().");
00731     TEST_FOR_EXCEPTION(hasColMap() == false, std::runtime_error,
00732         Teuchos::typeName(*this) << "::insertLocalValues(): cannot insert local indices without a column map; ");
00733     TEST_FOR_EXCEPTION(values.size() != indices.size(), std::runtime_error,
00734         Teuchos::typeName(*this) << "::insertLocalValues(): values.size() must equal indices.size().");
00735     TEST_FOR_EXCEPTION(getRowMap()->isNodeLocalElement(localRow) == false, std::runtime_error,
00736         Teuchos::typeName(*this) << "::insertLocalValues(): row does not belong to this node.");
00737     Teuchos::Array<LocalOrdinal> finds;
00738     Teuchos::Array<Scalar>       fvals;
00739     finds.reserve(indices.size());
00740     fvals.reserve(values.size());
00741     // use column map to filter the entries:
00742     const Map<LocalOrdinal,GlobalOrdinal,Node> &cmap = *getColMap();
00743     for (size_t i=0; i < static_cast<size_t>(indices.size()); ++i) {
00744       if (cmap.isNodeLocalElement(indices[i])) {
00745         finds.push_back(indices[i]);
00746         fvals.push_back(values[i]);
00747       }
00748     }
00749     if (finds.size() > 0) {
00750       if (myGraph_->indicesAreAllocated() == false) {
00751         allocateValues(CrsGraph<LocalOrdinal,GlobalOrdinal,Node>::AllocateLocal, ST::zero(), GraphNotYetAllocated);
00752       }
00753       //
00754       ArrayRCP<LocalOrdinal> rowindsview;
00755       ArrayRCP<Scalar>       rowvalsview;
00756       RowInfo sizeInfo = myGraph_->getFullLocalViewNonConst(localRow, rowindsview);
00757       rowvalsview = getFullViewNonConst(localRow, sizeInfo);
00758       const size_t newSize = sizeInfo.numEntries + finds.size();
00759       if (newSize > sizeInfo.allocSize) {
00760         TEST_FOR_EXCEPTION(myGraph_->getProfileType() == StaticProfile, std::runtime_error,
00761             Teuchos::typeName(*this) << "::insertLocalValues(): new indices exceed statically allocated graph structure.");
00762         TPETRA_EFFICIENCY_WARNING(true,std::runtime_error,
00763             "::insertLocalValues(): Pre-allocated space has been exceeded, requiring new allocation. To improve efficiency, suggest larger allocation.");
00764         // update allocation only as much as necessary
00765         myGraph_->updateLocalAllocation(localRow,newSize);
00766         updateAllocation(localRow,newSize);
00767         // get new views; inefficient, but acceptible in this already inefficient case
00768         sizeInfo = myGraph_->getFullLocalViewNonConst(localRow, rowindsview);
00769         rowvalsview = getFullViewNonConst(localRow, sizeInfo);
00770       }
00771       // add new indices, values to graph and matrix rows
00772       myGraph_->insertLocalIndicesViaView( localRow, finds, rowindsview + sizeInfo.numEntries );
00773       std::copy( fvals.begin(), fvals.end(), rowvalsview + sizeInfo.numEntries);
00774 #ifdef HAVE_TPETRA_DEBUG
00775       {
00776         RowInfo sizeInfoNew = myGraph_->getRowInfo(localRow);
00777         TEST_FOR_EXCEPTION(sizeInfoNew.numEntries != newSize, std::logic_error,
00778             Teuchos::typeName(*this) << "::insertLocalValues(): Internal logic error. Please contact Tpetra team.");
00779       }
00780 #endif
00781       // 
00782       rowindsview = Teuchos::null;
00783       rowvalsview = Teuchos::null;
00784       // checkInternalState();
00785     }
00786   }
00787 
00788 
00791   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00792   void CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::insertGlobalValues(GlobalOrdinal globalRow, 
00793                          const Teuchos::ArrayView<const GlobalOrdinal> &indices,
00794                          const Teuchos::ArrayView<const Scalar>        &values) {
00795     using Teuchos::ArrayRCP;
00796     TEST_FOR_EXCEPTION(isStaticGraph() == true, std::runtime_error,
00797         Teuchos::typeName(*this) << "::insertGlobalValues(): matrix was constructed with static graph. Cannot insert new entries.");
00798     TEST_FOR_EXCEPTION(myGraph_->isLocallyIndexed() == true, std::runtime_error,
00799         Teuchos::typeName(*this) << "::insertGlobalValues(): graph indices are local; use insertLocalValues().");
00800     TEST_FOR_EXCEPTION(values.size() != indices.size(), std::runtime_error,
00801         Teuchos::typeName(*this) << "::insertGlobalValues(): values.size() must equal indices.size().");
00802     const LocalOrdinal myRow = getRowMap()->getLocalElement(globalRow);
00803     if (myRow != LOT::invalid()) {
00804       // if we have a column map, use it to filter the entries.
00805       // only filter if this is our row.
00806       Teuchos::Array<GlobalOrdinal> finds_is_temporary;
00807       Teuchos::Array<Scalar>        fvals_is_temporary;
00808       Teuchos::ArrayView<const GlobalOrdinal> findices = indices;
00809       Teuchos::ArrayView<const Scalar       > fvalues  = values;
00810       if (hasColMap()) {
00811         // filter indices and values through the column map
00812         finds_is_temporary.reserve(indices.size());
00813         fvals_is_temporary.reserve(values.size());
00814         const Map<LocalOrdinal,GlobalOrdinal,Node> &cmap = *getColMap();
00815         for (size_t i=0; i< static_cast<size_t>(indices.size()); ++i) {
00816           if (cmap.isNodeGlobalElement(indices[i])) {
00817             finds_is_temporary.push_back(indices[i]);
00818             fvals_is_temporary.push_back(values[i]);
00819           }
00820         }
00821         findices = finds_is_temporary();
00822         fvalues  = fvals_is_temporary();
00823       }
00824       // add the new indices and values
00825       if (findices.size() > 0) {
00826         if (myGraph_->indicesAreAllocated() == false) {
00827           allocateValues(CrsGraph<LocalOrdinal,GlobalOrdinal,Node>::AllocateGlobal, ST::zero(), GraphNotYetAllocated);
00828         }
00829         // 
00830         ArrayRCP<GlobalOrdinal> rowindsview;
00831         ArrayRCP<Scalar>        rowvalsview;
00832         RowInfo sizeInfo = myGraph_->getFullGlobalViewNonConst(myRow, rowindsview);
00833         rowvalsview = getFullViewNonConst(myRow, sizeInfo);
00834         const size_t newSize = sizeInfo.numEntries + findices.size();
00835         if (newSize > sizeInfo.allocSize) {
00836           TEST_FOR_EXCEPTION(myGraph_->getProfileType() == StaticProfile, std::runtime_error,
00837               Teuchos::typeName(*this) << "::insertGlobalValues(): new indices exceed statically allocated graph structure.");
00838           TPETRA_EFFICIENCY_WARNING(true,std::runtime_error,
00839               "::insertGlobalValues(): Pre-allocated space has been exceeded, requiring new allocation. To improve efficiency, suggest larger allocation.");
00840           // update allocation only as much as necessary
00841           myGraph_->updateGlobalAllocation(myRow,newSize);
00842           updateAllocation(myRow,newSize);
00843           // get new views; inefficient, but acceptible in this already inefficient case
00844           sizeInfo = myGraph_->getFullGlobalViewNonConst(myRow, rowindsview);
00845           rowvalsview = getFullViewNonConst(myRow, sizeInfo);
00846         }
00847         // add new indices, values to graph and matrix rows
00848         myGraph_->insertGlobalIndicesViaView( myRow, findices, rowindsview + sizeInfo.numEntries );
00849         std::copy(  fvalues.begin(),  fvalues.end(), rowvalsview + sizeInfo.numEntries );
00850 #ifdef HAVE_TPETRA_DEBUG
00851         {
00852           RowInfo sizeInfoNew = myGraph_->getRowInfo(myRow);
00853           TEST_FOR_EXCEPTION(sizeInfoNew.numEntries != newSize, std::logic_error,
00854               Teuchos::typeName(*this) << "::insertGlobalValues(): Internal logic error. Please contact Tpetra team.");
00855         }
00856 #endif
00857         // 
00858         rowindsview = Teuchos::null;
00859         rowvalsview = Teuchos::null;
00860         // checkInternalState();
00861       }
00862     }
00863     else {
00864       typename Teuchos::ArrayView<const GlobalOrdinal>::iterator ind = indices.begin();
00865       typename Teuchos::ArrayView<const Scalar       >::iterator val =  values.begin();
00866       nonlocals_[globalRow].reserve( nonlocals_[globalRow].size() + indices.size() );
00867       for (; val != values.end(); ++val, ++ind) {
00868         nonlocals_[globalRow].push_back(std::make_pair(*ind, *val));
00869       }
00870     }
00871   }
00872 
00873 
00876   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00877   void CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::replaceGlobalValues(      
00878                                         GlobalOrdinal globalRow, 
00879                                         const Teuchos::ArrayView<const GlobalOrdinal> &indices,
00880                                         const Teuchos::ArrayView<const Scalar>        &values) {
00881     using Teuchos::ArrayRCP;
00882     // find the values for the specified indices
00883     // if the row is not ours, throw an exception
00884     // ignore values not in the matrix (indices not found)
00885     // operate whether indices are local or global
00886     const size_t STINV = Teuchos::OrdinalTraits<size_t>::invalid();
00887     TEST_FOR_EXCEPTION(values.size() != indices.size(), std::runtime_error,
00888         Teuchos::typeName(*this) << "::replaceGlobalValues(): values.size() must equal indices.size().");
00889     typename Teuchos::ArrayView<const GlobalOrdinal>::iterator ind = indices.begin();
00890     typename Teuchos::ArrayView<const        Scalar>::iterator val = values.begin();
00891     LocalOrdinal lrow = getRowMap()->getLocalElement(globalRow);
00892     TEST_FOR_EXCEPTION(lrow == LOT::invalid(), std::runtime_error,
00893         Teuchos::typeName(*this) << "::replaceGlobalValues(): specified global row does not belong to this processor.");
00894     Teuchos::RCP< const CrsGraph<LocalOrdinal,GlobalOrdinal,Node> > graph = getCrsGraph();
00895     // 
00896     if (isLocallyIndexed() == true) {
00897       ArrayRCP<const LocalOrdinal> lindrowview;
00898       RowInfo sizeInfo = graph->getFullLocalView(lrow, lindrowview);
00899       ArrayRCP<Scalar> valrowview = getFullViewNonConst(lrow, sizeInfo);
00900       while (ind != indices.end()) {
00901         LocalOrdinal lind = getColMap()->getLocalElement(*ind);
00902         size_t loc = graph->findLocalIndex(lrow,lind,lindrowview);
00903         if (loc != STINV) {
00904           valrowview[loc] = (*val);
00905         }
00906         ++ind;
00907         ++val;
00908       }
00909       valrowview = Teuchos::null;
00910       lindrowview = Teuchos::null;
00911     }
00912     else if (isGloballyIndexed() == true) {
00913       ArrayRCP<const GlobalOrdinal> gindrowview;
00914       RowInfo sizeInfo = graph->getFullGlobalView(lrow, gindrowview);
00915       ArrayRCP<Scalar> valrowview = getFullViewNonConst(lrow, sizeInfo);
00916       while (ind != indices.end()) {
00917         size_t loc = graph->findGlobalIndex(lrow,*ind,gindrowview);
00918         if (loc != STINV) {
00919           valrowview[loc] = (*val);
00920         }
00921         ++ind;
00922         ++val;
00923       }
00924       valrowview = Teuchos::null;
00925       gindrowview = Teuchos::null;
00926     }
00927     //else {
00928     // graph indices are not allocated, i.e., are non-existant; nothing to do
00929     //}
00930   }
00931 
00932 
00935   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00936   void CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::replaceLocalValues(      
00937                                         LocalOrdinal localRow,
00938                                         const Teuchos::ArrayView<const LocalOrdinal> &indices,
00939                                         const Teuchos::ArrayView<const Scalar>        &values) {
00940     using Teuchos::ArrayRCP;
00941     // find the values for the specified indices
00942     // if the row is not ours, throw an exception
00943     // ignore values not in the matrix (indices not found)
00944     // operate whether indices are local or global
00945     const size_t STINV = Teuchos::OrdinalTraits<size_t>::invalid();
00946     TEST_FOR_EXCEPTION(values.size() != indices.size(), std::runtime_error,
00947         Teuchos::typeName(*this) << "::replaceLocalValues(): values.size() must equal indices.size().");
00948     typename Teuchos::ArrayView<const LocalOrdinal>::iterator ind = indices.begin();
00949     typename Teuchos::ArrayView<const       Scalar>::iterator val = values.begin();
00950     bool isLocalRow = getRowMap()->isNodeLocalElement(localRow);
00951     TEST_FOR_EXCEPTION(isLocalRow == false, std::runtime_error,
00952         Teuchos::typeName(*this) << "::replaceLocalValues(): specified local row does not belong to this processor.");
00953     Teuchos::RCP< const CrsGraph<LocalOrdinal,GlobalOrdinal,Node> > graph = getCrsGraph();
00954     // 
00955     if (isLocallyIndexed() == true) {
00956       ArrayRCP<const LocalOrdinal> lindrowview;
00957       RowInfo sizeInfo = graph->getFullLocalView(localRow, lindrowview);
00958       ArrayRCP<Scalar> valrowview = getFullViewNonConst(localRow, sizeInfo);
00959       while (ind != indices.end()) {
00960         LocalOrdinal lind = *ind;
00961         size_t loc = graph->findLocalIndex(localRow,lind,lindrowview);
00962         if (loc != STINV) {
00963           valrowview[loc] = (*val);
00964         }
00965         ++ind;
00966         ++val;
00967       }
00968       valrowview = Teuchos::null;
00969       lindrowview = Teuchos::null;
00970     }
00971     else if (isGloballyIndexed() == true) {
00972       ArrayRCP<const GlobalOrdinal> gindrowview;
00973       GlobalOrdinal g_row = graph->getRowMap()->getGlobalElement(localRow);
00974       RowInfo sizeInfo = graph->getFullGlobalView(g_row, gindrowview);
00975       ArrayRCP<Scalar> valrowview = getFullViewNonConst(localRow, sizeInfo);
00976       while (ind != indices.end()) {
00977         GlobalOrdinal gind = graph->getColMap()->getGlobalElement(*ind);
00978         size_t loc = graph->findGlobalIndex(g_row,gind,gindrowview);
00979         if (loc != STINV) {
00980           valrowview[loc] = (*val);
00981         }
00982         ++ind;
00983         ++val;
00984       }
00985       valrowview = Teuchos::null;
00986       gindrowview = Teuchos::null;
00987     }
00988     //else {
00989     // graph indices are not allocated, i.e., are non-existant; nothing to do
00990     //}
00991   }
00992 
00993 
00996   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00997   void CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::sumIntoGlobalValues(GlobalOrdinal globalRow, 
00998                          const Teuchos::ArrayView<const GlobalOrdinal> &indices,
00999                          const Teuchos::ArrayView<const Scalar>        &values) {
01000     using Teuchos::ArrayRCP;
01001     // find the values for the specified indices
01002     // if the row is not ours, throw an exception
01003     // ignore values not in the matrix (indices not found)
01004     // operate whether indices are local or global
01005     const size_t STINV = Teuchos::OrdinalTraits<size_t>::invalid();
01006     TEST_FOR_EXCEPTION(values.size() != indices.size(), std::runtime_error,
01007         Teuchos::typeName(*this) << "::sumIntoGlobalValues(): values.size() must equal indices.size().");
01008     typename Teuchos::ArrayView<const GlobalOrdinal>::iterator ind = indices.begin();
01009     typename Teuchos::ArrayView<const        Scalar>::iterator val = values.begin();
01010     LocalOrdinal lrow = getRowMap()->getLocalElement(globalRow);
01011     TEST_FOR_EXCEPTION(lrow == LOT::invalid(), std::runtime_error,
01012         Teuchos::typeName(*this) << "::sumIntoGlobalValues(): specified global row does not belong to this processor.");
01013     Teuchos::RCP< const CrsGraph<LocalOrdinal,GlobalOrdinal,Node> > graph = getCrsGraph();
01014     //
01015     if (isLocallyIndexed() == true) {
01016       ArrayRCP<const LocalOrdinal> lindrowview;
01017       RowInfo sizeInfo = graph->getFullLocalView(lrow, lindrowview);
01018       ArrayRCP<Scalar> valrowview = getFullViewNonConst(lrow, sizeInfo);
01019       while (ind != indices.end()) {
01020         LocalOrdinal lind = getColMap()->getLocalElement(*ind);
01021         size_t loc = graph->findLocalIndex(lrow,lind,lindrowview);
01022         if (loc != STINV) {
01023           valrowview[loc] += (*val);
01024         }
01025         ++ind;
01026         ++val;
01027       }
01028       valrowview = Teuchos::null;
01029       lindrowview = Teuchos::null;
01030     }
01031     else if (isGloballyIndexed() == true) {
01032       ArrayRCP<const GlobalOrdinal> gindrowview;
01033       RowInfo sizeInfo = graph->getFullGlobalView(lrow, gindrowview);
01034       ArrayRCP<Scalar> valrowview = getFullViewNonConst(lrow, sizeInfo);
01035       while (ind != indices.end()) {
01036         size_t loc = graph->findGlobalIndex(lrow,*ind,gindrowview);
01037         if (loc != STINV) {
01038           valrowview[loc] += (*val);
01039         }
01040         ++ind;
01041         ++val;
01042       }
01043       valrowview = Teuchos::null;
01044       gindrowview = Teuchos::null;
01045     }
01046     //else {
01047       // indices are not allocated; nothing to do
01048     //}
01049   }
01050 
01051 
01054   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
01055   void CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getLocalRowCopy(
01056                                 LocalOrdinal LocalRow, 
01057                                 const Teuchos::ArrayView<LocalOrdinal> &indices, 
01058                                 const Teuchos::ArrayView<Scalar>       &values,
01059                                 size_t &numEntries) const {
01060     using Teuchos::ArrayRCP;
01061     TEST_FOR_EXCEPTION(isGloballyIndexed()==true && hasColMap()==false, std::runtime_error,
01062         Teuchos::typeName(*this) << "::getLocalRowCopy(): local indices cannot be produced.");
01063     TEST_FOR_EXCEPTION(getRowMap()->isNodeLocalElement(LocalRow) == false, std::runtime_error,
01064         Teuchos::typeName(*this) << "::getLocalRowCopy(LocalRow,...): specified row (==" << LocalRow << ") is not valid on this node.");
01065     Teuchos::RCP< const CrsGraph<LocalOrdinal,GlobalOrdinal,Node> > graph = getCrsGraph();
01066     if (graph->isLocallyIndexed()) {
01067       ArrayRCP<const LocalOrdinal> indrowview;
01068       ArrayRCP<const Scalar>       valrowview; 
01069       RowInfo sizeInfo = graph->getFullLocalView(LocalRow, indrowview);
01070       valrowview = getFullView(LocalRow, sizeInfo);
01071       numEntries = sizeInfo.numEntries;
01072       TEST_FOR_EXCEPTION(static_cast<size_t>(indices.size()) < numEntries || static_cast<size_t>(values.size()) < numEntries, std::runtime_error, 
01073           Teuchos::typeName(*this) << "::getLocalRowCopy(LocalRow,indices,values): size of indices,values must be sufficient to store the specified row.");
01074       if (numEntries > 0) {
01075         std::copy( indrowview.begin(), indrowview.begin() + numEntries, indices.begin() );
01076         std::copy( valrowview.begin(), valrowview.begin() + numEntries, values.begin() );
01077       }
01078       valrowview = Teuchos::null;
01079       indrowview = Teuchos::null;
01080     }
01081     else if (graph->isGloballyIndexed()) {
01082       ArrayRCP<const GlobalOrdinal> indrowview;
01083       ArrayRCP<const Scalar>        valrowview; 
01084       RowInfo sizeInfo = graph->getFullGlobalView(LocalRow, indrowview);
01085       valrowview = getFullView(LocalRow, sizeInfo);
01086       numEntries = sizeInfo.numEntries;
01087       TEST_FOR_EXCEPTION(static_cast<size_t>(indices.size()) < numEntries || static_cast<size_t>(values.size()) < numEntries, std::runtime_error, 
01088           Teuchos::typeName(*this) << "::getLocalRowCopy(LocalRow,indices,values): size of indices,values must be sufficient to store the specified row.");
01089       if (numEntries > 0) {
01090         std::copy( valrowview.begin(), valrowview.begin() + numEntries, values.begin() );
01091       }
01092       for (size_t j=0; j < numEntries; ++j) {
01093         indices[j] = getColMap()->getLocalElement(indrowview[j]);
01094       }
01095       valrowview = Teuchos::null;
01096       indrowview = Teuchos::null;
01097     }
01098     else {
01099 #ifdef HAVE_TPETRA_DEBUG
01100       // should have fallen in one of the above if indices are allocated
01101       TEST_FOR_EXCEPTION( graph->indicesAreAllocated() == true, std::logic_error, 
01102           Teuchos::typeName(*this) << "::getLocalRowCopy(): Internal logic error. Please contact Tpetra team.");
01103 #endif
01104       numEntries = 0;
01105     }
01106   }
01107 
01108 
01111   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
01112   void CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getGlobalRowCopy(
01113                                 GlobalOrdinal globalRow, 
01114                                 const Teuchos::ArrayView<GlobalOrdinal> &indices,
01115                                 const Teuchos::ArrayView<Scalar>        &values,
01116                                 size_t &numEntries) const {
01117     using Teuchos::ArrayRCP;
01118     // Only locally owned rows can be queried, otherwise complain
01119     LocalOrdinal myRow = getRowMap()->getLocalElement(globalRow);
01120     TEST_FOR_EXCEPTION(myRow == LOT::invalid(), std::runtime_error,
01121         Teuchos::typeName(*this) << "::getGlobalRowCopy(globalRow,...): globalRow does not belong to this node.");
01122     Teuchos::RCP< const CrsGraph<LocalOrdinal,GlobalOrdinal,Node> > graph = getCrsGraph();
01123     if (graph->isGloballyIndexed()) {
01124       ArrayRCP<const GlobalOrdinal> indrowview;
01125       ArrayRCP<const Scalar>        valrowview; 
01126       RowInfo sizeInfo = graph->getFullGlobalView(myRow, indrowview);
01127       valrowview = getFullView(myRow, sizeInfo);
01128       numEntries = sizeInfo.numEntries;
01129       TEST_FOR_EXCEPTION(static_cast<size_t>(indices.size()) < numEntries || static_cast<size_t>(values.size()) < numEntries, std::runtime_error, 
01130           Teuchos::typeName(*this) << "::getGlobalRowCopy(GlobalRow,indices,values): size of indices,values must be sufficient to store the specified row.");
01131       if (numEntries > 0) {
01132         std::copy( indrowview.begin(), indrowview.begin() + numEntries, indices.begin() );
01133         std::copy( valrowview.begin(), valrowview.begin() + numEntries, values.begin() );
01134       }
01135       valrowview = Teuchos::null;
01136       indrowview = Teuchos::null;
01137     }
01138     else if (graph->isLocallyIndexed()) {
01139       ArrayRCP<const LocalOrdinal> indrowview;
01140       ArrayRCP<const Scalar>       valrowview; 
01141       RowInfo sizeInfo = graph->getFullLocalView(myRow, indrowview);
01142       valrowview = getFullView(myRow, sizeInfo);
01143       numEntries = sizeInfo.numEntries;
01144       TEST_FOR_EXCEPTION(static_cast<size_t>(indices.size()) < numEntries || static_cast<size_t>(values.size()) < numEntries, std::runtime_error, 
01145           Teuchos::typeName(*this) << "::getGlobalRowCopy(GlobalRow,indices,values): size of indices,values must be sufficient to store the specified row.");
01146       if (numEntries > 0) {
01147         std::copy( valrowview.begin(), valrowview.begin() + numEntries, values.begin() );
01148       }
01149       for (size_t j=0; j < numEntries; ++j) {
01150         indices[j] = getColMap()->getGlobalElement(indrowview[j]);
01151       }
01152       valrowview = Teuchos::null;
01153       indrowview = Teuchos::null;
01154     }
01155     else {
01156 #ifdef HAVE_TPETRA_DEBUG
01157       // should have fallen in one of the above if indices are allocated
01158       TEST_FOR_EXCEPTION( graph->indicesAreAllocated() == true, std::logic_error, 
01159           Teuchos::typeName(*this) << "::getGlobalRowCopy(): Internal logic error. Please contact Tpetra team.");
01160 #endif
01161       numEntries = 0;
01162     }
01163   }
01164 
01165 
01168   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
01169   void CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getGlobalRowView(
01170                                 GlobalOrdinal GlobalRow, 
01171                                 Teuchos::ArrayRCP<const GlobalOrdinal> &indices,
01172                                 Teuchos::ArrayRCP<const Scalar>        &values) const {
01173     TEST_FOR_EXCEPTION(isLocallyIndexed() == true, std::runtime_error,
01174         Teuchos::typeName(*this) << "::getGlobalRowView(): global indices do not exist; call getLocalRowView().");
01175     LocalOrdinal lrow = getRowMap()->getLocalElement(GlobalRow);
01176     TEST_FOR_EXCEPTION(lrow == LOT::invalid(), std::runtime_error,
01177         Teuchos::typeName(*this) << "::getGlobalRowView(GlobalRow,...): GlobalRow (== " << GlobalRow << ") does not belong to this node.");
01178     Teuchos::RCP< const CrsGraph<LocalOrdinal,GlobalOrdinal,Node> > graph = getCrsGraph();
01179     RowInfo sizeInfo = graph->getFullGlobalView(lrow,indices);
01180     values = getFullView(lrow,sizeInfo);
01181     if (sizeInfo.numEntries != sizeInfo.allocSize) {
01182 #ifdef HAVE_TPETRA_DEBUG
01183       TEST_FOR_EXCEPTION( indices == Teuchos::null && values != Teuchos::null, std::logic_error, 
01184           Teuchos::typeName(*this) << "::getGlobalRowView(): Internal logic error. Please contact Tpetra team.");
01185 #endif
01186       if (indices != Teuchos::null) {
01187         indices = indices.persistingView(0,sizeInfo.numEntries);
01188         if (values != Teuchos::null) {
01189           values  =  values.persistingView(0,sizeInfo.numEntries);
01190         }
01191       }
01192     }
01193     return;
01194   }
01195 
01196 
01199   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
01200   void CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getLocalRowView(
01201                                 LocalOrdinal LocalRow, 
01202                                 Teuchos::ArrayRCP<const LocalOrdinal> &indices,
01203                                 Teuchos::ArrayRCP<const Scalar>        &values) const {
01204     TEST_FOR_EXCEPTION(isGloballyIndexed() == true, std::runtime_error,
01205         Teuchos::typeName(*this) << "::getLocalRowView(): local indices do not exist; call getGlobalRowView().");
01206     TEST_FOR_EXCEPTION(getRowMap()->isNodeLocalElement(LocalRow) == false, std::runtime_error,
01207         Teuchos::typeName(*this) << "::getLocalRowView(LocalRow,...): LocalRow (== " << LocalRow << ") is not valid on this node.");
01208     Teuchos::RCP< const CrsGraph<LocalOrdinal,GlobalOrdinal,Node> > graph = getCrsGraph();
01209     RowInfo sizeInfo = graph->getFullLocalView(LocalRow,indices);
01210     values = getFullView(LocalRow,sizeInfo);
01211     if (sizeInfo.numEntries != sizeInfo.allocSize) {
01212 #ifdef HAVE_TPETRA_DEBUG
01213       TEST_FOR_EXCEPTION( indices == Teuchos::null && values != Teuchos::null, std::logic_error, 
01214           Teuchos::typeName(*this) << "::getLocalRowView(): Internal logic error. Please contact Tpetra team.");
01215 #endif
01216       if (indices != Teuchos::null) {
01217         indices = indices.persistingView(0,sizeInfo.numEntries);
01218         if (values != Teuchos::null) {
01219           values  =  values.persistingView(0,sizeInfo.numEntries);
01220         }
01221       }
01222     }
01223     return;
01224   }
01225 
01226 
01229   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
01230   void CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::scale(const Scalar &alpha) {
01231     // scale all values in the matrix
01232     // it is easiest to scale all allocated values, instead of scaling only the ones with valid entries
01233     // however, if there are no valid entries, we can short-circuit
01234     // furthermore, if the values aren't allocated, we can short-circuit (unallocated values are zero, scaling to zero)
01235     Teuchos::RCP< const CrsGraph<LocalOrdinal,GlobalOrdinal,Node> > graph = getCrsGraph();
01236     const size_t     nlrs = graph->getNodeNumRows(),
01237                  numAlloc = graph->getNodeAllocationSize(),
01238                numEntries = graph->getNodeNumEntries();
01239     if (graph->indicesAreAllocated() == false || numAlloc == 0 || numEntries == 0) {
01240       // do nothing
01241     }
01242     else {
01243       // do the scale in parallel
01244       Teuchos::RCP<Node> node = lclMatrix_.getNode();
01245       Kokkos::ReadyBufferHelper<Node> rbh(node);
01246       Kokkos::SingleScaleOp<Scalar> wdp;
01247       wdp.alpha = alpha;
01248       if (graph->getProfileType() == StaticProfile) {
01249         rbh.begin();
01250         wdp.x = rbh.addNonConstBuffer(values1D_);
01251         rbh.end();
01252         node->template parallel_for<Kokkos::SingleScaleOp<Scalar> >(0,numAlloc,wdp);
01253       }
01254       else if (graph->getProfileType() == DynamicProfile) {
01255         for (size_t row=0; row < nlrs; ++row) {
01256           if (values2D_[row] != Teuchos::null) {
01257             rbh.begin();
01258             wdp.x = rbh.addNonConstBuffer(values2D_[row]);
01259             rbh.end();
01260             node->template parallel_for<Kokkos::SingleScaleOp<Scalar> >(0,values2D_[row].size(),wdp);
01261           }
01262         }
01263       }
01264     }
01265   }
01266 
01267 
01270   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
01271   void CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::setAllToScalar(const Scalar &alpha) {
01272     // set all values in the matrix
01273     // it is easiest to set all allocated values, instead of setting only the ones with valid entries
01274     // however, if there are no entries, we can short-circuit
01275     // this method is equivalent replacing all valid entries with alpha.
01276     Teuchos::RCP< const CrsGraph<LocalOrdinal,GlobalOrdinal,Node> > graph = getCrsGraph();
01277     const size_t     nlrs = graph->getNodeNumRows(),
01278                  numAlloc = graph->getNodeAllocationSize(),
01279                numEntries = graph->getNodeNumEntries();
01280     if (graph->indicesAreAllocated() == false || numAlloc == 0 || numEntries == 0) {
01281       // do nothing
01282     }
01283     else {
01284 #ifdef HAVE_TPETRA_DEBUG
01285       TEST_FOR_EXCEPTION( graph->indicesAreAllocated() == false, std::logic_error, 
01286           Teuchos::typeName(*this) << "::setAllToScalar(): internal logic error. Please contact Tpetra team.");
01287 #endif
01288       // set the values in parallel
01289       Teuchos::RCP<Node> node = lclMatrix_.getNode();
01290       Kokkos::ReadyBufferHelper<Node> rbh(node);
01291       Kokkos::InitOp<Scalar> wdp;
01292       wdp.alpha = alpha;
01293       if (graph->getProfileType() == StaticProfile) {
01294         rbh.begin();
01295         wdp.x = rbh.addNonConstBuffer(values1D_);
01296         rbh.end();
01297         node->template parallel_for<Kokkos::InitOp<Scalar> >(0,numAlloc,wdp);
01298       }
01299       else if (graph->getProfileType() == DynamicProfile) {
01300         for (size_t row=0; row < nlrs; ++row) {
01301           if (values2D_[row] != Teuchos::null) {
01302             rbh.begin();
01303             wdp.x = rbh.addNonConstBuffer(values2D_[row]);
01304             rbh.end();
01305             node->template parallel_for<Kokkos::InitOp<Scalar> >(0,values2D_[row].size(),wdp);
01306           }
01307         }
01308       }
01309     }
01310   }
01311 
01312 
01315   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
01316   void CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getLocalDiagCopy(Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> &dvec) const {
01317     TEST_FOR_EXCEPTION(isFillComplete() == false, std::runtime_error,
01318         Teuchos::typeName(*this) << ": cannot call getLocalDiagCopy() until fillComplete() has been called.");
01319     TEST_FOR_EXCEPTION(dvec.getMap()->isSameAs(*getRowMap()) == false, std::runtime_error,
01320         Teuchos::typeName(*this) << "::getLocalDiagCopy(dvec): dvec must have the same map as the CrsMatrix.");
01321     Teuchos::RCP< const CrsGraph<LocalOrdinal,GlobalOrdinal,Node> > graph = getCrsGraph();
01322     const size_t STINV = Teuchos::OrdinalTraits<size_t>::invalid();
01323 #ifdef HAVE_TPETRA_DEBUG
01324     size_t numDiagFound = 0;
01325 #endif
01326     const size_t nlrs = getNodeNumRows();
01327     Teuchos::ArrayRCP<Scalar> vecView = dvec.get1dViewNonConst();
01328     for (size_t r=0; r < nlrs; ++r) {
01329       vecView[r] = Teuchos::ScalarTraits<Scalar>::zero();
01330       GlobalOrdinal rgid = getRowMap()->getGlobalElement(r);
01331       if (getColMap()->isNodeGlobalElement(rgid)) {
01332         LocalOrdinal rlid = getColMap()->getLocalElement(rgid);
01333         Teuchos::ArrayRCP<const LocalOrdinal> inds;
01334         RowInfo sizeInfo = graph->getFullLocalView(r,inds);
01335         if (sizeInfo.numEntries > 0) {
01336           const size_t j = graph->findLocalIndex(r, rlid, inds);
01337           if (j != STINV) {
01338             Teuchos::ArrayRCP<const Scalar> vals;
01339             vals = getFullView(r,sizeInfo);
01340             vecView[r] = vals[j];
01341             vals = Teuchos::null;
01342 #ifdef HAVE_TPETRA_DEBUG
01343             ++numDiagFound;
01344 #endif
01345           }
01346         }
01347         inds = Teuchos::null;
01348       }
01349     }
01350     vecView = Teuchos::null;
01351 #ifdef HAVE_TPETRA_DEBUG
01352     TEST_FOR_EXCEPTION(numDiagFound != getNodeNumDiags(), std::logic_error, 
01353         "CrsMatrix::getLocalDiagCopy(): logic error. Please contact Tpetra team.");
01354 #endif
01355   }
01356 
01357 
01360   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
01361   void CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::globalAssemble() {
01362     using Teuchos::arcp;
01363     using Teuchos::OrdinalTraits;
01364     using Teuchos::Array;
01365     using Teuchos::SerialDenseMatrix;
01366     using Teuchos::ArrayView;
01367     using Teuchos::ArrayRCP;
01368     using Teuchos::rcp;
01369     using Teuchos::outArg;
01370     using std::pair;
01371     using std::make_pair;
01372     using Teuchos::tuple;
01373     typedef typename std::map<GlobalOrdinal,Teuchos::Array<pair<GlobalOrdinal,Scalar> > >::const_iterator NLITER;
01374     typedef typename Teuchos::Array<pair<GlobalOrdinal,Scalar> >::const_iterator NLRITER;
01375     const int numImages = getComm()->getSize();
01376     const int myImageID = getComm()->getRank();
01377     // Determine if any nodes have global entries to share
01378     size_t MyNonlocals = nonlocals_.size(), 
01379            MaxGlobalNonlocals;
01380     Teuchos::reduceAll<int,size_t>(*getComm(),Teuchos::REDUCE_MAX,MyNonlocals,
01381       outArg(MaxGlobalNonlocals));
01382     if (MaxGlobalNonlocals == 0) return;  // no entries to share
01383 
01384     // compute a list of NLRs from nonlocals_ and use it to compute:
01385     //      IdsAndRows: a vector of (id,row) pairs
01386     //          NLR2Id: a map from NLR to the Id that owns it
01387     // globalNeighbors: a global graph of connectivity between images: globalNeighbors(i,j) indicates that j sends to i
01388     //         sendIDs: a list of all images I send to
01389     //         recvIDs: a list of all images I receive from (constructed later)
01390     Array<pair<int,GlobalOrdinal> > IdsAndRows;
01391     std::map<GlobalOrdinal,int> NLR2Id;
01392     SerialDenseMatrix<int,char> globalNeighbors;
01393     Array<int> sendIDs, recvIDs;
01394     {
01395       // nonlocals_ contains the entries we are holding for all non-local rows
01396       // we want a list of the rows for which we have data
01397       Array<GlobalOrdinal> NLRs;
01398       std::set<GlobalOrdinal> setOfRows;
01399       for (NLITER iter = nonlocals_.begin(); iter != nonlocals_.end(); ++iter)
01400       {
01401         setOfRows.insert(iter->first);
01402       }
01403       // copy the elements in the set into an Array
01404       NLRs.resize(setOfRows.size());
01405       std::copy(setOfRows.begin(), setOfRows.end(), NLRs.begin());
01406 
01407       // get a list of ImageIDs for the non-local rows (NLRs)
01408       Array<int> NLRIds(NLRs.size());
01409       {
01410         LookupStatus stat = getRowMap()->getRemoteIndexList(NLRs(),NLRIds());
01411         char lclerror = ( stat == IDNotPresent ? 1 : 0 );
01412         char gblerror;
01413         Teuchos::reduceAll(*getComm(),Teuchos::REDUCE_MAX,lclerror,outArg(gblerror));
01414         TEST_FOR_EXCEPTION(gblerror, std::runtime_error,
01415             Teuchos::typeName(*this) << "::globalAssemble(): non-local entries correspond to invalid rows.");
01416       }
01417 
01418       // build up a list of neighbors, as well as a map between NLRs and Ids
01419       // localNeighbors[i] != 0 iff I have data to send to image i
01420       // put NLRs,Ids into an array of pairs
01421       IdsAndRows.reserve(NLRs.size());
01422       Array<char> localNeighbors(numImages,0);
01423       typename Array<GlobalOrdinal>::const_iterator nlr;
01424       typename Array<int>::const_iterator id;
01425       for (nlr = NLRs.begin(), id = NLRIds.begin();
01426            nlr != NLRs.end(); ++nlr, ++id) 
01427       {
01428         NLR2Id[*nlr] = *id;
01429         localNeighbors[*id] = 1;
01430         IdsAndRows.push_back(make_pair<int,GlobalOrdinal>(*id,*nlr));
01431       }
01432       for (int j=0; j<numImages; ++j)
01433       {
01434         if (localNeighbors[j]) {
01435           sendIDs.push_back(j);
01436         }
01437       }
01438       // sort IdsAndRows, by Ids first, then rows
01439       std::sort(IdsAndRows.begin(),IdsAndRows.end());
01440       // gather from other nodes to form the full graph
01441       globalNeighbors.shapeUninitialized(numImages,numImages);
01442       Teuchos::gatherAll(*getComm(),numImages,localNeighbors.getRawPtr(),numImages*numImages,globalNeighbors.values());
01443       // globalNeighbors at this point contains (on all images) the
01444       // connectivity between the images. 
01445       // globalNeighbors(i,j) != 0 means that j sends to i/that i receives from j
01446     }
01447 
01449     // FIGURE OUT WHO IS SENDING TO WHOM AND HOW MUCH
01450     // DO THIS IN THE PROCESS OF PACKING ALL OUTGOING DATA ACCORDING TO DESTINATION ID
01452 
01453     // loop over all columns to know from which images I can expect to receive something
01454     for (int j=0; j<numImages; ++j)
01455     {
01456       if (globalNeighbors(myImageID,j)) {
01457         recvIDs.push_back(j);
01458       }
01459     }
01460     size_t numRecvs = recvIDs.size();
01461 
01462     // we know how many we're sending to already
01463     // form a contiguous list of all data to be sent
01464     // track the number of entries for each ID
01465     Array<CrsIJV<GlobalOrdinal,Scalar> > IJVSendBuffer;
01466     Array<size_t> sendSizes(sendIDs.size(), 0);
01467     size_t numSends = 0;
01468     for (typename Array<pair<int,GlobalOrdinal> >::const_iterator IdAndRow = IdsAndRows.begin();
01469          IdAndRow != IdsAndRows.end(); ++IdAndRow) 
01470     {
01471       int            id = IdAndRow->first;
01472       GlobalOrdinal row = IdAndRow->second;
01473       // have we advanced to a new send?
01474       if (sendIDs[numSends] != id) {
01475         numSends++;
01476         TEST_FOR_EXCEPTION(sendIDs[numSends] != id, std::logic_error, Teuchos::typeName(*this) << "::globalAssemble(): internal logic error. Contact Tpetra team.");
01477       }
01478       // copy data for row into contiguous storage
01479       for (NLRITER jv = nonlocals_[row].begin(); jv != nonlocals_[row].end(); ++jv)
01480       {
01481         IJVSendBuffer.push_back( CrsIJV<GlobalOrdinal,Scalar>(row,jv->first,jv->second) );
01482         sendSizes[numSends]++;
01483       }
01484     }
01485     if (IdsAndRows.size() > 0) {
01486       numSends++; // one last increment, to make it a count instead of an index
01487     }
01488     TEST_FOR_EXCEPTION(Teuchos::as<typename Array<int>::size_type>(numSends) != sendIDs.size(), std::logic_error, Teuchos::typeName(*this) << "::globalAssemble(): internal logic error. Contact Tpetra team.");
01489 
01490     // don't need this data anymore
01491     nonlocals_.clear();
01492 
01494     // TRANSMIT SIZE INFO BETWEEN SENDERS AND RECEIVERS
01496     // perform non-blocking sends: send sizes to our recipients
01497     Array<Teuchos::RCP<Teuchos::CommRequest> > sendRequests;
01498     for (size_t s=0; s < numSends ; ++s) {
01499       // we'll fake the memory management, because all communication will be local to this method and the scope of our data
01500       sendRequests.push_back( Teuchos::isend<int,size_t>(*getComm(),Teuchos::rcpFromRef(sendSizes[s]),sendIDs[s]) );
01501     }
01502     // perform non-blocking receives: receive sizes from our senders
01503     Array<Teuchos::RCP<Teuchos::CommRequest> > recvRequests;
01504     Array<size_t> recvSizes(numRecvs);
01505     for (size_t r=0; r < numRecvs; ++r) {
01506       // we'll fake the memory management, because all communication will be local to this method and the scope of our data
01507       recvRequests.push_back( Teuchos::ireceive(*getComm(),rcp(&recvSizes[r],false),recvIDs[r]) );
01508     }
01509     // wait on all 
01510     if (!sendRequests.empty()) {
01511       Teuchos::waitAll(*getComm(),sendRequests());
01512     }
01513     if (!recvRequests.empty()) {
01514       Teuchos::waitAll(*getComm(),recvRequests());
01515     }
01516     Teuchos::barrier(*getComm());
01517     sendRequests.clear();
01518     recvRequests.clear();
01519 
01521     // NOW SEND/RECEIVE ALL ROW DATA
01523     // from the size info, build the ArrayViews into IJVSendBuffer
01524     Array<ArrayView<CrsIJV<GlobalOrdinal,Scalar> > > sendBuffers(numSends,Teuchos::null);
01525     {
01526       size_t cur = 0;
01527       for (size_t s=0; s<numSends; ++s) {
01528         sendBuffers[s] = IJVSendBuffer(cur,sendSizes[s]);
01529         cur += sendSizes[s];
01530       }
01531     }
01532     // perform non-blocking sends
01533     for (size_t s=0; s < numSends ; ++s)
01534     {
01535       // we'll fake the memory management, because all communication will be local to this method and the scope of our data
01536       ArrayRCP<CrsIJV<GlobalOrdinal,Scalar> > tmparcp = arcp(sendBuffers[s].getRawPtr(),0,sendBuffers[s].size(),false);
01537       sendRequests.push_back( Teuchos::isend<int,CrsIJV<GlobalOrdinal,Scalar> >(*getComm(),tmparcp,sendIDs[s]) );
01538     }
01539     // calculate amount of storage needed for receives
01540     // setup pointers for the receives as well
01541     size_t totalRecvSize = std::accumulate(recvSizes.begin(),recvSizes.end(),0);
01542     Array<CrsIJV<GlobalOrdinal,Scalar> > IJVRecvBuffer(totalRecvSize);
01543     // from the size info, build the ArrayViews into IJVRecvBuffer
01544     Array<ArrayView<CrsIJV<GlobalOrdinal,Scalar> > > recvBuffers(numRecvs,Teuchos::null);
01545     {
01546       size_t cur = 0;
01547       for (size_t r=0; r<numRecvs; ++r) {
01548         recvBuffers[r] = IJVRecvBuffer(cur,recvSizes[r]);
01549         cur += recvSizes[r];
01550       }
01551     }
01552     // perform non-blocking recvs
01553     for (size_t r=0; r < numRecvs ; ++r)
01554     {
01555       // we'll fake the memory management, because all communication will be local to this method and the scope of our data
01556       ArrayRCP<CrsIJV<GlobalOrdinal,Scalar> > tmparcp = arcp(recvBuffers[r].getRawPtr(),0,recvBuffers[r].size(),false);
01557       recvRequests.push_back( Teuchos::ireceive(*getComm(),tmparcp,recvIDs[r]) );
01558     }
01559     // perform waits
01560     if (!sendRequests.empty()) {
01561       Teuchos::waitAll(*getComm(),sendRequests());
01562     }
01563     if (!recvRequests.empty()) {
01564       Teuchos::waitAll(*getComm(),recvRequests());
01565     }
01566     Teuchos::barrier(*getComm());
01567     sendRequests.clear();
01568     recvRequests.clear();
01569 
01570 
01572     // NOW PROCESS THE RECEIVED ROW DATA
01574     // TODO: instead of adding one entry at a time, add one row at a time.
01575     //       this requires resorting; they arrived sorted by sending node, so that entries could be non-contiguous if we received
01576     //       multiple entries for a particular row from different processors.
01577     //       it also requires restoring the data, which may make it not worth the trouble.
01578     for (typename Array<CrsIJV<GlobalOrdinal,Scalar> >::const_iterator ijv = IJVRecvBuffer.begin(); ijv != IJVRecvBuffer.end(); ++ijv)
01579     {
01580       try {
01581         insertGlobalValues(ijv->i, tuple(ijv->j), tuple(ijv->v));
01582       }
01583       catch (std::runtime_error &e) {
01584         std::ostringstream omsg;
01585         omsg << e.what() << std::endl
01586           << "caught in globalAssemble() in " << __FILE__ << ":" << __LINE__ << std::endl ;
01587         throw std::runtime_error(omsg.str());
01588       }
01589     }
01590 
01591     // WHEW! THAT WAS TIRING!
01592   }
01593 
01594 
01597   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
01598   void CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::fillComplete(OptimizeOption os) {
01599     fillComplete(getRowMap(),getRowMap(),os);
01600   }
01601 
01602 
01605   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
01606   void CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::fillComplete(
01607                                             const Teuchos::RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > &domainMap, 
01608                                             const Teuchos::RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > &rangeMap, 
01609                                             OptimizeOption os) {
01610 
01611     if (getCrsGraph()->indicesAreAllocated() == false) {
01612       // must allocate global, because we have no column map
01613       allocateValues( CrsGraph<LocalOrdinal,GlobalOrdinal,Node>::AllocateGlobal, Teuchos::ScalarTraits<Scalar>::zero(), GraphNotYetAllocated );
01614     }
01615 
01616     if (getComm()->getSize() > 1) {
01617       globalAssemble();
01618     }
01619     else {
01620       TEST_FOR_EXCEPTION(nonlocals_.size() > 0, std::runtime_error,
01621           Teuchos::typeName(*this) << "::fillComplete(): cannot have non-local entries on a serial run. Invalid entry was submitted to the CrsMatrix.");
01622     }
01623 
01624     // if the matrix was constructed with an un-optimized graph, it must not have been optimized in the interim.
01625     // this would likely have wrecked the fragile dance between graph and matrix.
01626     if (isStaticGraph()) {
01627       if (constructedWithOptimizedGraph_ == false && staticGraph_->isStorageOptimized() == true) {
01628         TPETRA_ABUSE_WARNING(true,std::runtime_error,
01629             "::fillComplete(): optimizeStorage() has been called on graph since matrix was constructed.");
01630         return;
01631       }
01632     }
01633     
01634     // if we're not allowed to change a static graph, then we can't call optimizeStorage() on it.
01635     // then we can't call fillComplete() on the matrix with DoOptimizeStorage
01636     // throw a warning. this is an unfortunate late evaluation; however, we couldn't know when we received the graph that the user would try to 
01637     // optimize the storage later on.
01638     if (os == DoOptimizeStorage && isStaticGraph() && staticGraph_->isStorageOptimized() == false) {
01639       TPETRA_ABUSE_WARNING(true,std::runtime_error,
01640           "::fillComplete(): requested optimized storage, but static graph does not have optimized storage. Ignoring request to optimize storage.");
01641       os = DoNotOptimizeStorage;
01642     }
01643 
01644     if (isStaticGraph() == false) {
01645       myGraph_->makeIndicesLocal(domainMap,rangeMap);
01646       sortEntries();
01647       mergeRedundantEntries();
01648       // can't optimize graph storage before optimizing our storage; therefore, do not optimize storage yet.
01649       myGraph_->fillComplete(domainMap,rangeMap,DoNotOptimizeStorage);
01650     }
01651 
01652     fillComplete_ = true;
01653 
01654     if (os == DoOptimizeStorage && isStorageOptimized() == false) {
01655       // this will also:
01656       // * call optimizeStorage() on the graph as well, if isStaticGraph() == false
01657       //   this will call fillLocalGraph() on the graph
01658       // * call fillLocalMatrix()
01659       //   this will initialize the local mat-vec and (if appropriate) the local solve
01660       optimizeStorage();
01661     }
01662     else { 
01663       // local graph already filled.
01664       // fill the local matrix.
01665       fillLocalMatrix();
01666     }
01667   
01668     checkInternalState();
01669   }
01670 
01671 
01674   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
01675   void CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::sortEntries() {
01676     using Teuchos::ArrayRCP;
01677     TEST_FOR_EXCEPTION(myGraph_->isGloballyIndexed() == true, std::logic_error,
01678         Teuchos::typeName(*this) << "::sortEntries(): sortEntries() must be called after indices are transformed to local.\n"
01679                                  << "Likely internal logic error. Please contact Tpetra team.");
01680     if (myGraph_->isSorted()) return;
01681     if (myGraph_->indicesAreAllocated() && myGraph_->getNodeAllocationSize() > 0) {
01682       const size_t nlrs = getNodeNumRows();
01683       for (size_t r=0; r < nlrs; ++r) {
01684         // TODO: This is slightly inefficient, because it may query pbuf_rowOffsets_ repeatadly. 
01685         //       However, it is very simple code. Consider rewriting it.
01686         Teuchos::ArrayRCP<LocalOrdinal> inds;
01687         Teuchos::ArrayRCP<Scalar>       vals;
01688         RowInfo sizeInfo = myGraph_->getFullLocalViewNonConst(r, inds);
01689         vals = getFullViewNonConst(r, sizeInfo);
01690         sort2(inds.begin(), inds.begin() + sizeInfo.numEntries, vals);
01691       }
01692     }
01693     myGraph_->setSorted(true);  // we just sorted them
01694     return;
01695   }
01696 
01697 
01700   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
01701   void CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::mergeRedundantEntries() {
01702     using Teuchos::ArrayRCP;
01703     TEST_FOR_EXCEPTION(myGraph_->isSorted() == false, std::runtime_error,
01704         Teuchos::typeName(*this) << "::mergeRedundantEntries() cannot be called before indices are sorted.\n"
01705                                  << "Likely interal logic error. Please contact Tpetra team.");
01706     if (myGraph_->notRedundant()) return;
01707     if (myGraph_->indicesAreAllocated() && myGraph_->getNodeAllocationSize() > 0) {
01708       for (size_t r=0; r<getNodeNumRows(); ++r) {
01709         // TODO: This is slightly inefficient, because it may query pbuf_rowOffsets_ repeatadly. 
01710         //       However, it is very simple code. Consider rewriting it.
01711         Teuchos::ArrayRCP<LocalOrdinal> inds;
01712         Teuchos::ArrayRCP<Scalar>       vals;
01713         RowInfo sizeInfo = myGraph_->getFullLocalViewNonConst(r, inds);
01714         vals = getFullViewNonConst(r, sizeInfo);
01715         if (sizeInfo.numEntries > 0) {
01716           size_t curEntry = 0;
01717           Scalar curValue = vals[curEntry];
01718           for (size_t k=1; k < sizeInfo.numEntries; ++k) {
01719             if (inds[k] == inds[k-1]) {
01720               curValue += vals[k];
01721             }
01722             else {
01723               vals[curEntry++] = curValue;
01724               curValue = vals[k];
01725             }
01726           }
01727           vals[curEntry] = curValue;
01728         }
01729         vals = Teuchos::null;
01730         inds = Teuchos::null;
01731       }
01732       myGraph_->removeRedundantIndices();
01733     }
01734   }
01735 
01736 
01739   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
01740   void CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::optimizeStorage() {
01741     // optimizeStorage will perform two functions:
01742     // 1) create a single allocation of memory
01743     // 2) pack data in that allocation
01744     // if getProfileType() == StaticProfile, then 1) has already been done
01745 #ifdef HAVE_TPETRA_DEBUG
01746     std::string err = Teuchos::typeName(*this) + "::optimizeStorage(): Internal logic error. Please contact Tpetra team.";
01747 #endif
01748     if (isStorageOptimized() == true) return;
01749     // if storage is not optimized and we have a static graph, we will not be able to optimize the storage. throw an exception.
01750     TEST_FOR_EXCEPTION( isStaticGraph(), std::logic_error, 
01751         Teuchos::typeName(*this) << "::optimizeStorage(): Cannot optimize storage with a static graph. Possible internal logic error. Please contact Tpetra team.");
01752     TEST_FOR_EXCEPTION(isFillComplete() == false || myGraph_->isSorted() == false || myGraph_->notRedundant() == false, std::runtime_error,
01753         Teuchos::typeName(*this) << "::optimizeStorage(): fillComplete() must be called before optimizeStorage().");
01754     // 
01755     Teuchos::RCP<Node> node = lclMatrix_.getNode();
01756     const size_t       nlrs = getNodeNumRows(),
01757                  numEntries = myGraph_->getNodeNumEntries();
01758     if (nlrs > 0) {
01759       if (myGraph_->getProfileType() == DynamicProfile) {
01760         // allocate a single memory block
01761         if (numEntries > 0) {
01762           // numEntries > 0  =>  allocSize > 0  =>  values2D_ != Teuchos::null
01763 #ifdef HAVE_TPETRA_DEBUG
01764           TEST_FOR_EXCEPTION( values2D_ == Teuchos::null, std::logic_error, err);
01765 #endif
01766           values1D_ = node->template allocBuffer<Scalar>(numEntries);
01767           size_t sofar = 0;
01768           for (size_t row=0; row<nlrs; ++row) {
01769             RowInfo sizeInfo = myGraph_->getRowInfo(row);
01770             if (sizeInfo.numEntries > 0) {
01771               KOKKOS_NODE_TRACE("CrsMatrix::optimizeStorage()")
01772               node->template copyBuffers<Scalar>( sizeInfo.numEntries, values2D_[row], values1D_ + sofar );
01773               values2D_[row] = Teuchos::null;
01774             }
01775             sofar += sizeInfo.numEntries;
01776           }
01777 #ifdef HAVE_TPETRA_DEBUG
01778           TEST_FOR_EXCEPTION(sofar != numEntries, std::logic_error, err);
01779 #endif
01780         }
01781         // done with 2D allocation
01782         values2D_ = Teuchos::null;
01783       }
01784       else {
01785         // storage is already allocated; just need to pack
01786         if (numEntries > 0) {
01787 #ifdef HAVE_TPETRA_DEBUG
01788           TEST_FOR_EXCEPTION( values1D_ == Teuchos::null, std::logic_error, err);
01789 #endif
01790           size_t sofar = 0;
01791           for (size_t row=0; row<nlrs; ++row) {
01792             RowInfo sizeInfo = myGraph_->getRowInfo(row);
01793             if (sizeInfo.numEntries > 0 && sizeInfo.offset1D != sofar) {
01794               KOKKOS_NODE_TRACE("CrsMatrix::optimizeStorage()")
01795               node->template copyBuffers<Scalar>( sizeInfo.numEntries, 
01796                                                   values1D_ + sizeInfo.offset1D,
01797                                                   values1D_ + sofar );
01798             }
01799             sofar += sizeInfo.numEntries;
01800           }
01801           values1D_ = values1D_.persistingView(0,sofar);
01802 #ifdef HAVE_TPETRA_DEBUG
01803           TEST_FOR_EXCEPTION(sofar != numEntries, std::logic_error, err);
01804 #endif
01805         }
01806       }
01807     }
01808 
01809     // now we can optimize the graph storage; this process is symmetric w.r.t. the process undertaken in the matrix
01810     myGraph_->optimizeStorage();
01811 
01812     // if empty, then release buffers
01813     // this mirrors similar code in CrsGraph::optimizeStorage, and is necessary so that 
01814     // both the local matrix and local graph are marked empty
01815     if (myGraph_->getNodeAllocationSize() == 0) {
01816       values1D_ = Teuchos::null;
01817     }
01818 
01819     // local graph was filled during myGraph_->optimizeStorage()
01820     fillLocalMatrix(); 
01821   }
01822 
01823 
01826   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
01827   void CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::apply(
01828                                         const MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> &X, 
01829                                         MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> &Y, 
01830                                         Teuchos::ETransp mode, Scalar alpha, Scalar beta) const {
01831     TEST_FOR_EXCEPTION( isFillComplete() == false, std::runtime_error, 
01832         Teuchos::typeName(*this) << "::apply(): cannot call apply() until fillComplete() has been called.");
01833     sameScalarMultiplyOp_->apply(X,Y,mode,alpha,beta);
01834   }
01835 
01836 
01839   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
01840   template <class DomainScalar, class RangeScalar>
01841   void CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::multiply(
01842                                         const MultiVector<DomainScalar,LocalOrdinal,GlobalOrdinal,Node> &X, 
01843                                               MultiVector<RangeScalar,LocalOrdinal,GlobalOrdinal,Node> &Y, 
01844                                               Teuchos::ETransp mode, RangeScalar alpha, RangeScalar beta) const {
01845     typedef Teuchos::ScalarTraits<RangeScalar> RST;
01846     const Kokkos::MultiVector<DomainScalar,Node> *lclX = &X.getLocalMV();
01847     Kokkos::MultiVector<RangeScalar,Node>        *lclY = &Y.getLocalMVNonConst();
01848 #ifdef HAVE_TPETRA_DEBUG
01849     TEST_FOR_EXCEPTION(!isFillComplete(), std::runtime_error, 
01850         Teuchos::typeName(*this) << ": cannot call multiply() until fillComplete() has been called.");
01851     TEST_FOR_EXCEPTION(X.getNumVectors() != Y.getNumVectors(), std::runtime_error,
01852         Teuchos::typeName(*this) << "::multiply(X,Y): X and Y must have the same number of vectors.");
01853     TEST_FOR_EXCEPTION(X.isConstantStride() == false || Y.isConstantStride() == false, std::runtime_error,
01854         Teuchos::typeName(*this) << "::multiply(X,Y): X and Y must be constant stride.");
01855     TEST_FOR_EXCEPTION(lclX==lclY, std::runtime_error,
01856         Teuchos::typeName(*this) << "::multiply(X,Y): X and Y cannot share data.");
01857 #endif
01858     //
01859     // Call the matvec
01860     if (beta == RST::zero()) {
01861       // Y = alpha*op(M)*X with overwrite semantics
01862       lclMatOps_.template multiply<DomainScalar,RangeScalar>(mode, alpha, *lclX, *lclY);
01863     }
01864     else {
01865       // Y = alpha*op(M) + beta*Y
01866       lclMatOps_.template multiply<DomainScalar,RangeScalar>(mode, alpha, *lclX, beta, *lclY);
01867     }
01868   }
01869 
01870 
01873   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
01874   template <class DomainScalar, class RangeScalar>
01875   void CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::solve(
01876                                     const MultiVector<RangeScalar,LocalOrdinal,GlobalOrdinal,Node>  &Y, 
01877                                           MultiVector<DomainScalar,LocalOrdinal,GlobalOrdinal,Node> &X,
01878                                           Teuchos::ETransp mode) const {
01879     const Kokkos::MultiVector<RangeScalar,Node> *lclY = &Y.getLocalMV();
01880     Kokkos::MultiVector<DomainScalar,Node>      *lclX = &X.getLocalMVNonConst();
01881 #ifdef HAVE_TPETRA_DEBUG
01882     TEST_FOR_EXCEPTION(!isFillComplete(), std::runtime_error, 
01883         Teuchos::typeName(*this) << ": cannot call solve() until fillComplete() has been called.");
01884     TEST_FOR_EXCEPTION(X.getNumVectors() != Y.getNumVectors(), std::runtime_error,
01885         Teuchos::typeName(*this) << "::solve(X,Y): X and Y must have the same number of vectors.");
01886     TEST_FOR_EXCEPTION(X.isConstantStride() == false || Y.isConstantStride() == false, std::runtime_error,
01887         Teuchos::typeName(*this) << "::solve(X,Y): X and Y must be constant stride.");
01888     TEST_FOR_EXCEPTION(isUpperTriangular() == false && isLowerTriangular() == false, std::runtime_error,
01889         Teuchos::typeName(*this) << "::solve(): can only solve() triangular matrices.");
01890     TEST_FOR_EXCEPTION(Teuchos::ScalarTraits<Scalar>::isComplex && mode == Teuchos::TRANS, std::logic_error,
01891         Teuchos::typeName(*this) << "::solve() does not currently support transposed solve for complex scalar types.");
01892 #endif
01893     //
01894     // Call the solve
01895     Teuchos::EDiag diag = ( getNodeNumDiags() < getNodeNumRows() ? Teuchos::UNIT_DIAG : Teuchos::NON_UNIT_DIAG );
01896     if (mode == Teuchos::NO_TRANS) {
01897       if (isUpperTriangular()) {
01898         lclMatOps_.template solve<DomainScalar,RangeScalar>(Teuchos::NO_TRANS, Teuchos::UPPER_TRI, diag, *lclY, *lclX);
01899       }
01900       else {
01901         lclMatOps_.template solve<DomainScalar,RangeScalar>(Teuchos::NO_TRANS, Teuchos::LOWER_TRI, diag, *lclY, *lclX);
01902       }
01903     }
01904     else {
01905       if (isUpperTriangular()) {
01906         lclMatOps_.template solve<DomainScalar,RangeScalar>(Teuchos::CONJ_TRANS, Teuchos::UPPER_TRI, diag, *lclY, *lclX);
01907       }
01908       else {
01909         lclMatOps_.template solve<DomainScalar,RangeScalar>(Teuchos::CONJ_TRANS, Teuchos::LOWER_TRI, diag, *lclY, *lclX);
01910       }
01911     }
01912   }
01913 
01914 
01917   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
01918   std::string CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::description() const {
01919     std::ostringstream oss;
01920     oss << DistObject<char, LocalOrdinal,GlobalOrdinal,Node>::description();
01921     if (isFillComplete()) {
01922       oss << "{status = fill complete"
01923           << ", global rows = " << getGlobalNumRows()
01924           << ", global cols = " << getGlobalNumCols()
01925           << ", global num entries = " << getGlobalNumEntries()
01926           << "}";
01927     }
01928     else {
01929       oss << "{status = fill not complete"
01930           << ", global rows = " << getGlobalNumRows()
01931           << "}";
01932     }
01933     return oss.str();
01934   }
01935 
01936 
01939   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
01940   void CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel) const {
01941     using std::endl;
01942     using std::setw;
01943     using Teuchos::VERB_DEFAULT;
01944     using Teuchos::VERB_NONE;
01945     using Teuchos::VERB_LOW;
01946     using Teuchos::VERB_MEDIUM;
01947     using Teuchos::VERB_HIGH;
01948     using Teuchos::VERB_EXTREME;
01949     Teuchos::EVerbosityLevel vl = verbLevel;
01950     if (vl == VERB_DEFAULT) vl = VERB_LOW;
01951     Teuchos::RCP<const Teuchos::Comm<int> > comm = this->getComm();
01952     const int myImageID = comm->getRank(),
01953               numImages = comm->getSize();
01954     size_t width = 1;
01955     for (size_t dec=10; dec<getGlobalNumRows(); dec *= 10) {
01956       ++width;
01957     }
01958     width = std::max<size_t>(width,11) + 2;
01959     Teuchos::OSTab tab(out);
01960     //    none: print nothing
01961     //     low: print O(1) info from node 0
01962     //  medium: print O(P) info, num entries per node
01963     //    high: print O(N) info, num entries per row
01964     // extreme: print O(NNZ) info: print indices and values
01965     // 
01966     // for medium and higher, print constituent objects at specified verbLevel
01967     Teuchos::RCP< const CrsGraph<LocalOrdinal,GlobalOrdinal,Node> > graph = getCrsGraph();
01968     if (vl != VERB_NONE) {
01969       if (myImageID == 0) out << this->description() << std::endl; 
01970       // O(1) globals, minus what was already printed by description()
01971       if (isFillComplete() && myImageID == 0) {
01972         out << "Global number of diagonals = " << getGlobalNumDiags() << std::endl;
01973         out << "Global max number of entries = " << getGlobalMaxNumRowEntries() << std::endl;
01974       }
01975       // constituent objects
01976       if (vl == VERB_MEDIUM || vl == VERB_HIGH || vl == VERB_EXTREME) {
01977         if (myImageID == 0) out << "\nRow map: " << std::endl;
01978         getRowMap()->describe(out,vl);
01979         //
01980         if (getColMap() != Teuchos::null) {
01981           if (getColMap() == getRowMap()) {
01982             if (myImageID == 0) out << "\nColumn map is row map.";
01983           }
01984           else {
01985             if (myImageID == 0) out << "\nColumn map: " << std::endl;
01986             getColMap()->describe(out,vl);
01987           }
01988         }
01989         if (getDomainMap() != Teuchos::null) {
01990           if (getDomainMap() == getRowMap()) {
01991             if (myImageID == 0) out << "\nDomain map is row map.";
01992           }
01993           else if (getDomainMap() == getColMap()) {
01994             if (myImageID == 0) out << "\nDomain map is row map.";
01995           }
01996           else {
01997             if (myImageID == 0) out << "\nDomain map: " << std::endl;
01998             getDomainMap()->describe(out,vl);
01999           }
02000         }
02001         if (getRangeMap() != Teuchos::null) {
02002           if (getRangeMap() == getDomainMap()) {
02003             if (myImageID == 0) out << "\nRange map is domain map." << std::endl;
02004           }
02005           else if (getRangeMap() == getRowMap()) {
02006             if (myImageID == 0) out << "\nRange map is row map." << std::endl;
02007           }
02008           else {
02009             if (myImageID == 0) out << "\nRange map: " << std::endl;
02010             getRangeMap()->describe(out,vl);
02011           }
02012         }
02013         if (myImageID == 0) out << std::endl;
02014       }
02015       // O(P) data
02016       if (vl == VERB_MEDIUM || vl == VERB_HIGH || vl == VERB_EXTREME) {
02017         for (int imageCtr = 0; imageCtr < numImages; ++imageCtr) {
02018           if (myImageID == imageCtr) {
02019             out << "Node ID = " << imageCtr << std::endl;
02020             if (graph->indicesAreAllocated() == false) {
02021               out << "Node not allocated" << std::endl;
02022             }
02023             else {
02024               out << "Node number of allocated entries = " << graph->getNodeAllocationSize() << std::endl;
02025             }
02026             out << "Node number of entries = " << getNodeNumEntries() << std::endl;
02027             if (isFillComplete()) {
02028               out << "Node number of diagonals = " << getNodeNumDiags() << std::endl;
02029             }
02030             out << "Node max number of entries = " << getNodeMaxNumRowEntries() << std::endl;
02031           }
02032           comm->barrier();
02033           comm->barrier();
02034           comm->barrier();
02035         }
02036       }
02037       // O(N) and O(NNZ) data
02038       if (vl == VERB_HIGH || vl == VERB_EXTREME) {
02039         for (int imageCtr = 0; imageCtr < numImages; ++imageCtr) {
02040           if (myImageID == imageCtr) {
02041             out << std::setw(width) << "Node ID"
02042                 << std::setw(width) << "Global Row" 
02043                 << std::setw(width) << "Num Entries";
02044             if (vl == VERB_EXTREME) {
02045               out << std::setw(width) << "(Index,Value)";
02046             }
02047             out << std::endl;
02048             for (size_t r=0; r < getNodeNumRows(); ++r) {
02049               const size_t nE = getNumEntriesInLocalRow(r);
02050               GlobalOrdinal gid = getRowMap()->getGlobalElement(r);
02051               out << std::setw(width) << myImageID 
02052                   << std::setw(width) << gid
02053                   << std::setw(width) << nE;
02054               if (vl == VERB_EXTREME) {
02055                 if (isGloballyIndexed()) {
02056                   Teuchos::ArrayRCP<const GlobalOrdinal> rowinds;
02057                   Teuchos::ArrayRCP<const Scalar> rowvals;
02058                   getGlobalRowView(gid,rowinds,rowvals);
02059                   for (size_t j=0; j < nE; ++j) {
02060                     out << " (" << rowinds[j]
02061                         << ", " << rowvals[j]
02062                         << ") ";
02063                   }
02064                 }
02065                 else if (isLocallyIndexed()) {
02066                   Teuchos::ArrayRCP<const LocalOrdinal> rowinds;
02067                   Teuchos::ArrayRCP<const Scalar> rowvals;
02068                   getLocalRowView(r,rowinds,rowvals);
02069                   for (size_t j=0; j < nE; ++j) {
02070                     out << " (" << getColMap()->getGlobalElement(rowinds[j]) 
02071                         << ", " << rowvals[j]
02072                         << ") ";
02073                   }
02074                 }
02075               }
02076               out << std::endl;
02077             }
02078           }
02079           comm->barrier();
02080           comm->barrier();
02081           comm->barrier();
02082         }
02083       }
02084     }
02085   }
02086 
02087 
02090   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
02091   bool CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::checkSizes(const DistObject<char, LocalOrdinal,GlobalOrdinal,Node> & source)
02092   {
02093     // It's not clear what kind of compatibility checks on sizes can be performed here.
02094     // Epetra_CrsGraph doesn't check any sizes for compatibility.
02095 
02096     // right now, we'll only support import/exporting between CrsMatrix<Scalar>
02097     // if the source dist object isn't CrsMatrix or some offspring, flag this operation as incompatible.
02098     try  {
02099       const CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node,LocalMatOps> & A = dynamic_cast<const CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node,LocalMatOps> &>(source);
02100       (void)A;
02101     }
02102     catch (...) {
02103       return false;
02104     }
02105     return true;
02106   }
02107 
02108 
02111   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
02112   void CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::copyAndPermute(
02113                           const DistObject<char, LocalOrdinal,GlobalOrdinal,Node> & source,
02114                           size_t numSameIDs,
02115                           const Teuchos::ArrayView<const LocalOrdinal> &permuteToLIDs,
02116                           const Teuchos::ArrayView<const LocalOrdinal> &permuteFromLIDs)
02117   {
02118     using Teuchos::Array;
02119     using Teuchos::ArrayRCP;
02120     // this should succeed, because we already tested compatibility in checkSizes()
02121     const CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node,LocalMatOps> & src_mat = dynamic_cast<const CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node,LocalMatOps> &>(source);
02122     TEST_FOR_EXCEPTION(permuteToLIDs.size() != permuteFromLIDs.size(), std::runtime_error,
02123         Teuchos::typeName(*this) << "::copyAndPermute: permuteToLIDs and permuteFromLIDs must have the same size.");
02124     const bool src_is_locally_indexed = src_mat.isLocallyIndexed();
02125 
02126     // do numSame: copy the first numSame row from the source to *this
02127     // specifically, copy rows corresponding to Local Elements 0,numSame-1
02128     Array<GlobalOrdinal> row_indices;
02129     Array<Scalar>        row_values;
02130     LocalOrdinal mylid = 0;
02131     for (size_t i=0; i<numSameIDs; ++i, ++mylid) {
02132       // get Global ID for this row
02133       GlobalOrdinal gid = src_mat.getMap()->getGlobalElement(mylid);
02134       if (src_is_locally_indexed) {
02135         const size_t row_length = src_mat.getNumEntriesInGlobalRow(gid);
02136         row_indices.resize( row_length );
02137         row_values.resize( row_length );
02138         size_t check_row_length = 0;
02139         src_mat.getGlobalRowCopy(gid, row_indices(), row_values(), check_row_length);
02140 #ifdef HAVE_TPETRA_DEBUG
02141         TEST_FOR_EXCEPTION(row_length != check_row_length, std::logic_error,
02142             Teuchos::typeName(*this) << "::copyAndPermute(): Internal logic error. Please contact Tpetra team.");
02143 #endif
02144         insertGlobalValues( gid, row_indices(), row_values() );
02145       }
02146       else {
02147         ArrayRCP<const GlobalOrdinal> row_inds; 
02148         ArrayRCP<const Scalar>        row_vals; 
02149         src_mat.getGlobalRowView(gid, row_inds, row_vals);
02150         insertGlobalValues( gid, row_inds(), row_vals() );
02151       }
02152     }
02153 
02154     // handle the permuted rows.
02155     for (size_t p=0; p<(size_t)permuteToLIDs.size(); ++p) {
02156       const GlobalOrdinal  mygid =   this->getMap()->getGlobalElement(permuteToLIDs[p]);
02157       const GlobalOrdinal srcgid = src_mat.getMap()->getGlobalElement(permuteFromLIDs[p]);
02158       if (src_is_locally_indexed) {
02159         const size_t row_length = src_mat.getNumEntriesInGlobalRow(srcgid);
02160         row_indices.resize( row_length );
02161         row_values.resize( row_length );
02162         size_t check_row_length = 0;
02163         src_mat.getGlobalRowCopy(srcgid, row_indices(), row_values(), check_row_length);
02164 #ifdef HAVE_TPETRA_DEBUG
02165         TEST_FOR_EXCEPTION(row_length != check_row_length, std::logic_error,
02166             Teuchos::typeName(*this) << "::copyAndPermute(): Internal logic error. Please contact Tpetra team.");
02167 #endif
02168         insertGlobalValues( mygid, row_indices(), row_values() );
02169       }
02170       else {
02171         ArrayRCP<const GlobalOrdinal> row_inds;
02172         ArrayRCP<const Scalar>        row_vals;
02173         src_mat.getGlobalRowView( srcgid, row_inds, row_vals);
02174         insertGlobalValues( mygid, row_inds(), row_vals());
02175       }
02176     }
02177   }
02178 
02179 
02182   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
02183   void CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::packAndPrepare(
02184                           const DistObject<char, LocalOrdinal,GlobalOrdinal,Node> & source,
02185                           const Teuchos::ArrayView<const LocalOrdinal> &exportLIDs,
02186                           Teuchos::Array<char> &exports,
02187                           const Teuchos::ArrayView<size_t> & numPacketsPerLID,
02188                           size_t& constantNumPackets,
02189                           Distributor &distor)
02190   {
02191     using Teuchos::Array;
02192     using Teuchos::ArrayRCP;
02193     using Teuchos::ArrayView;
02194 
02195     TEST_FOR_EXCEPTION(exportLIDs.size() != numPacketsPerLID.size(), std::runtime_error,
02196         Teuchos::typeName(*this) << "::packAndPrepare: exportLIDs and numPacketsPerLID must have the same size.");
02197     // this should succeed, because we already tested compatibility in checkSizes() and performed this cast in packAndPrepare()
02198     const CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node,LocalMatOps> & src_mat = dynamic_cast<const CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node,LocalMatOps> &>(source);
02199     const bool src_is_locally_indexed = src_mat.isLocallyIndexed();
02200     constantNumPackets = 0;
02201 
02202     // first, set the contents of numPacketsPerLID, and accumulate a total-num-packets:
02203     // grab the max row size, while we're at it. may need it below.
02204     // Subtle: numPacketsPerLID is for byte-packets, so it needs to be multiplied
02205     const size_t SizeOfOrdValPair = sizeof(GlobalOrdinal)+sizeof(Scalar);
02206     size_t totalNumEntries = 0;
02207     size_t maxExpRowLength = 0;
02208     for (size_t i=0; i<(size_t)exportLIDs.size(); ++i) {
02209       GlobalOrdinal expGID = src_mat.getMap()->getGlobalElement(exportLIDs[i]);
02210       const size_t row_length = src_mat.getNumEntriesInGlobalRow(expGID);
02211       numPacketsPerLID[i] = row_length * SizeOfOrdValPair;
02212       totalNumEntries += row_length;
02213       maxExpRowLength = (row_length > maxExpRowLength ? row_length : maxExpRowLength);
02214     }
02215 
02216     // Need to do the following:
02217     // [inds_row0 vals_row0 inds_row1 vals_row1 ... inds_rowN vals_rowN]
02218     if (totalNumEntries > 0) {
02219       // exports is an array of char (bytes). it needs room for all of the indices and values
02220       const size_t totalNumBytes = totalNumEntries * SizeOfOrdValPair;
02221       exports.resize(totalNumBytes);
02222 
02223       ArrayView<char> avIndsC, avValsC;
02224       ArrayView<GlobalOrdinal> avInds;
02225       ArrayView<Scalar>        avVals;
02226 
02227       // now loop again and pack rows of indices into exports:
02228       // if global indices exist in the source, then we can use view semantics
02229       // otherwise, we are forced to use copy semantics (for the indices; for simplicity, we'll use them for values as well)
02230       size_t curOffsetInBytes = 0;
02231       if (src_is_locally_indexed) {
02232         Array<GlobalOrdinal> row_inds(maxExpRowLength);
02233         Array<Scalar>        row_vals(maxExpRowLength);
02234         for (size_t i=0; i<(size_t)exportLIDs.size(); ++i) {
02235           // get copy
02236           const GlobalOrdinal GID = src_mat.getMap()->getGlobalElement(exportLIDs[i]);
02237           size_t rowSize;
02238           src_mat.getGlobalRowCopy(GID, row_inds(), row_vals(), rowSize);
02239           // get export views
02240           avIndsC = exports(curOffsetInBytes,rowSize*sizeof(GlobalOrdinal));
02241           avValsC = exports(curOffsetInBytes+rowSize*sizeof(GlobalOrdinal),rowSize*sizeof(Scalar));
02242           avInds = Teuchos::av_reinterpret_cast<GlobalOrdinal>(avIndsC);
02243           avVals = Teuchos::av_reinterpret_cast<Scalar       >(avValsC);
02244           // copy
02245           std::copy( row_inds.begin(), row_inds.begin()+rowSize, avInds.begin());
02246           std::copy( row_vals.begin(), row_vals.begin()+rowSize, avVals.begin());
02247           curOffsetInBytes += SizeOfOrdValPair * rowSize;
02248         }
02249       }
02250       else {
02251         ArrayRCP<const GlobalOrdinal> row_inds;
02252         ArrayRCP<const Scalar>        row_vals;
02253         for (size_t i=0; i<(size_t)exportLIDs.size(); ++i) {
02254           // get view
02255           const GlobalOrdinal GID = src_mat.getMap()->getGlobalElement(exportLIDs[i]);
02256           src_mat.getGlobalRowView(GID, row_inds, row_vals);
02257           const size_t rowSize = (size_t)row_inds.size();
02258           // get export views
02259           avIndsC = exports(curOffsetInBytes,rowSize*sizeof(GlobalOrdinal));
02260           avValsC = exports(curOffsetInBytes+rowSize*sizeof(GlobalOrdinal),rowSize*sizeof(Scalar));
02261           avInds = Teuchos::av_reinterpret_cast<GlobalOrdinal>(avIndsC);
02262           avVals = Teuchos::av_reinterpret_cast<Scalar       >(avValsC);
02263           // copy
02264           std::copy( row_inds.begin(), row_inds.end(), avInds.begin());
02265           std::copy( row_vals.begin(), row_vals.end(), avVals.begin());
02266           curOffsetInBytes += SizeOfOrdValPair * rowSize;
02267         }
02268       }
02269 #ifdef HAVE_TPETRA_DEBUG
02270       TEST_FOR_EXCEPTION(curOffsetInBytes != totalNumBytes, std::logic_error,
02271           Teuchos::typeName(*this) << "::packAndPrepare(): Internal logic error. Please contact Tpetra team.");
02272 #endif
02273     }
02274   }
02275 
02276 
02279   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
02280   void CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::unpackAndCombine(
02281                             const Teuchos::ArrayView<const LocalOrdinal> &importLIDs,
02282                             const Teuchos::ArrayView<const char> &imports,
02283                             const Teuchos::ArrayView<size_t> &numPacketsPerLID,
02284                             size_t constantNumPackets,
02285                             Distributor & /* distor */,
02286                             CombineMode /* CM */)
02287   {
02288     using Teuchos::ArrayView;
02289     // We are not checking the value of the CombineMode input-argument.
02290     // Any incoming column-indices are inserted into the target graph. In this context, CombineMode values
02291     // of ADD vs INSERT are equivalent. What is the meaning of REPLACE for CrsGraph? If a duplicate column-index
02292     // is inserted, it will be compressed out when fillComplete is called.
02293     // NOTE: I have added a note to the Tpetra todo list to revisit this discussion. CGB, 6/18/2010
02294 
02295     TEST_FOR_EXCEPTION(importLIDs.size() != numPacketsPerLID.size(), std::runtime_error,
02296         Teuchos::typeName(*this) << "::unpackAndCombine: importLIDs and numPacketsPerLID must have the same size.");
02297 
02298     const size_t SizeOfOrdValPair = sizeof(GlobalOrdinal)+sizeof(Scalar);
02299     const size_t totalNumBytes = imports.size(); // * sizeof(char), which is one.
02300     const size_t totalNumEntries = totalNumBytes / SizeOfOrdValPair;
02301 
02302     if (totalNumEntries > 0) {
02303       // data packed as follows:
02304       // [inds_row0 vals_row0 inds_row1 vals_row1 ...]
02305       ArrayView<const char> avIndsC, avValsC;
02306       ArrayView<const GlobalOrdinal> avInds;
02307       ArrayView<const Scalar>        avVals;
02308 
02309       size_t curOffsetInBytes = 0;
02310       for (size_t i=0; i<(size_t)importLIDs.size(); ++i) {
02311         // get row info
02312         const LocalOrdinal LID = importLIDs[i];
02313         const GlobalOrdinal myGID = this->getMap()->getGlobalElement(LID);
02314         const size_t rowSize = numPacketsPerLID[i] / SizeOfOrdValPair;
02315         // get import views
02316         avIndsC = imports(curOffsetInBytes,rowSize*sizeof(GlobalOrdinal));
02317         avValsC = imports(curOffsetInBytes+rowSize*sizeof(GlobalOrdinal),rowSize*sizeof(Scalar));
02318         avInds = Teuchos::av_reinterpret_cast<const GlobalOrdinal>(avIndsC);
02319         avVals = Teuchos::av_reinterpret_cast<const Scalar       >(avValsC);
02320         // do insert
02321         insertGlobalValues(myGID, avInds(), avVals());
02322         curOffsetInBytes += rowSize * SizeOfOrdValPair;
02323       }
02324 #ifdef HAVE_TPETRA_DEBUG
02325       TEST_FOR_EXCEPTION(curOffsetInBytes != totalNumBytes, std::logic_error,
02326           Teuchos::typeName(*this) << "::packAndPrepare(): Internal logic error. Please contact Tpetra team.");
02327 #endif
02328     }
02329   }
02330 
02331 
02334   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
02335   void CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::createViews() const {
02336   }
02337 
02338 
02341   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
02342   void CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::createViewsNonConst(Kokkos::ReadWriteOption rwo) {
02343   }
02344 
02345 
02348   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
02349   void CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::releaseViews() const {
02350   }
02351 
02352 
02353 } // namespace Tpetra
02354 
02355 //
02356 // Explicit instantiation macro
02357 //
02358 // Must be expanded from within the Tpetra namespace!
02359 //
02360 
02361 #define TPETRA_CRSMATRIX_INSTANT(SCALAR,LO,GO,NODE) \
02362   \
02363   template class CrsMatrix< SCALAR , LO , GO , NODE >;
02364 
02365 #endif
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines
Generated on Wed Apr 13 10:21:41 2011 for Tpetra Matrix/Vector Services by  doxygen 1.6.3