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

Generated on Tue Jul 13 09:39:06 2010 for Tpetra Matrix/Vector Services by  doxygen 1.4.7