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