Tpetra Matrix/Vector Services Version of the Day
Tpetra_VbrMatrix_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_VBRMATRIX_DEF_HPP
00030 #define TPETRA_VBRMATRIX_DEF_HPP
00031 
00032 #include <algorithm>
00033 #include <sstream>
00034 
00035 #include <Kokkos_NodeHelpers.hpp>
00036 #include <Tpetra_Vector.hpp>
00037 
00038 #ifdef DOXYGEN_USE_ONLY
00039 #include "Tpetra_VbrMatrix_decl.hpp"
00040 #endif
00041 
00042 namespace Tpetra {
00043 
00044 template<class Scalar,class Node>
00045 void fill_device_ArrayRCP(Teuchos::RCP<Node>& node, Teuchos::ArrayRCP<Scalar>& ptr, Scalar value)
00046 {
00047   Kokkos::ReadyBufferHelper<Node> rbh(node);
00048   Kokkos::InitOp<Scalar> wdp;
00049   wdp.alpha = value;
00050   rbh.begin();
00051   wdp.x = rbh.addNonConstBuffer(ptr);
00052   rbh.end();
00053   node->template parallel_for<Kokkos::InitOp<Scalar> >(0, ptr.size(), wdp);
00054 }
00055 
00056 //-------------------------------------------------------------------
00057 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00058 VbrMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::VbrMatrix(const Teuchos::RCP<const BlockMap<LocalOrdinal,GlobalOrdinal,Node> > &blkRowMap, size_t maxNumEntriesPerRow, ProfileType pftype)
00059  : blkGraph_(Teuchos::rcp(new BlockCrsGraph<LocalOrdinal,GlobalOrdinal,Node>(blkRowMap, maxNumEntriesPerRow, pftype))),
00060    constBlkGraph_(blkGraph_),
00061    lclMatrix_(blkRowMap->getNodeNumBlocks(), blkRowMap->getPointMap()->getNode()),
00062    pbuf_values1D_(),
00063    pbuf_indx_(),
00064    lclMatOps_(blkRowMap->getPointMap()->getNode()),
00065    importer_(),
00066    exporter_(),
00067    importedVec_(),
00068    exportedVec_(),
00069    col_ind_2D_global_(Teuchos::rcp(new Teuchos::Array<std::map<GlobalOrdinal,Teuchos::ArrayRCP<Scalar> > >)),
00070    values2D_(Teuchos::rcp(new Teuchos::Array<Teuchos::Array<Teuchos::ArrayRCP<Scalar> > >)),
00071    is_fill_completed_(false),
00072    is_storage_optimized_(false)
00073 {
00074   //The graph of this VBR matrix is a BlockCrsGraph, which is a CrsGraph where
00075   //each entry in the graph corresponds to a block-entry in the matrix.
00076   //That is, you can think of a VBR matrix as a Crs matrix of dense
00077   //submatrices...
00078 }
00079 
00080 //-------------------------------------------------------------------
00081 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00082 VbrMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::VbrMatrix(const Teuchos::RCP<const BlockCrsGraph<LocalOrdinal,GlobalOrdinal,Node> > &blkGraph)
00083  : blkGraph_(),
00084    constBlkGraph_(blkGraph),
00085    lclMatrix_(blkGraph->getBlockRowMap()->getNodeNumBlocks(),
00086               blkGraph->getBlockRowMap()->getPointMap()->getNode()),
00087    pbuf_values1D_(),
00088    pbuf_indx_(),
00089    lclMatOps_(blkGraph->getBlockRowMap()->getPointMap()->getNode()),
00090    importer_(),
00091    exporter_(),
00092    importedVec_(),
00093    exportedVec_(),
00094    col_ind_2D_global_(Teuchos::rcp(new Teuchos::Array<std::map<GlobalOrdinal,Teuchos::ArrayRCP<Scalar> > >)),
00095    values2D_(Teuchos::rcp(new Teuchos::Array<Teuchos::Array<Teuchos::ArrayRCP<Scalar> > >)),
00096    is_fill_completed_(false),
00097    is_storage_optimized_(false)
00098 {
00099   //The graph of this VBR matrix is a BlockCrsGraph, which is a CrsGraph where
00100   //each entry in the graph corresponds to a block-entry in the matrix.
00101   //That is, you can think of a VBR matrix as a Crs matrix of dense
00102   //submatrices...
00103 
00104   TEST_FOR_EXCEPTION(blkGraph->isFillComplete() == false, std::runtime_error,
00105    "Tpetra::VbrMatrix::VbrMatrix(BlockCrsGraph) ERROR, this constructor requires graph.isFillComplete()==true.");
00106 
00107   createImporterExporter();
00108 
00109   is_fill_completed_ = true;
00110 
00111   optimizeStorage();
00112 
00113   fillLocalMatrix();
00114   fillLocalMatVec();
00115 }
00116 
00117 //-------------------------------------------------------------------
00118 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00119 VbrMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::~VbrMatrix()
00120 {
00121 }
00122 
00123 //-------------------------------------------------------------------
00124 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00125 const Teuchos::RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > &
00126 VbrMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getDomainMap() const
00127 {
00128   return getBlockDomainMap()->getPointMap();
00129 }
00130 
00131 //-------------------------------------------------------------------
00132 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00133 const Teuchos::RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > &
00134 VbrMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getRangeMap() const
00135 {
00136   return getBlockRangeMap()->getPointMap();
00137 }
00138 
00139 //-------------------------------------------------------------------
00140 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00141 bool
00142 VbrMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::isFillComplete() const
00143 {
00144   return is_fill_completed_;
00145 }
00146 
00147 //-------------------------------------------------------------------
00148 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00149 template <class DomainScalar, class RangeScalar>
00150 void
00151 VbrMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::multiply(const MultiVector<DomainScalar,LocalOrdinal,GlobalOrdinal,Node> & X, MultiVector<RangeScalar,LocalOrdinal,GlobalOrdinal,Node> &Y, Teuchos::ETransp trans, RangeScalar alpha, RangeScalar beta) const
00152 {
00153   TEST_FOR_EXCEPTION(!isFillComplete(), std::runtime_error,
00154     "Tpetra::VbrMatrix::multiply ERROR, multiply may only be called after fillComplete has been called.");
00155 
00156   const Kokkos::MultiVector<Scalar,Node> *lclX = &X.getLocalMV();
00157   Kokkos::MultiVector<Scalar,Node>        *lclY = &Y.getLocalMVNonConst();
00158 
00159   lclMatOps_.template multiply<Scalar,Scalar>(trans,alpha,*lclX,beta,*lclY);
00160 }
00161 
00162 //-------------------------------------------------------------------
00163 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00164 template<class DomainScalar, class RangeScalar>
00165 void
00166 VbrMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::solve(const MultiVector<RangeScalar,LocalOrdinal,GlobalOrdinal,Node>& Y, MultiVector<DomainScalar,LocalOrdinal,GlobalOrdinal,Node>& X, Teuchos::ETransp trans) const
00167 {
00168   TEST_FOR_EXCEPTION(X.isConstantStride() == false || Y.isConstantStride() == false, std::runtime_error,
00169         "Tpetra::VbrMatrix::solve(X,Y): X and Y must be constant stride.");
00170 
00171   TEST_FOR_EXCEPTION(!isFillComplete(), std::runtime_error,
00172     "Tpetra::VbrMatrix::solve ERROR, solve may only be called after fillComplete has been called.");
00173 
00174   TEST_FOR_EXCEPTION(constBlkGraph_->isUpperTriangular()==false && constBlkGraph_->isLowerTriangular()==false, std::runtime_error,
00175     "Tpetra::VbrMatrix::solve ERROR, matrix must be either upper or lower triangular.");
00176 
00177   const Kokkos::MultiVector<RangeScalar,Node> *lclY = &Y.getLocalMV();
00178   Kokkos::MultiVector<DomainScalar,Node>      *lclX = &X.getLocalMVNonConst();
00179 
00180   Teuchos::EUplo triang = constBlkGraph_->isUpperTriangular() ? Teuchos::UPPER_TRI : Teuchos::LOWER_TRI;
00181   Teuchos::EDiag diag = (constBlkGraph_->getNodeNumDiags() < constBlkGraph_->getNodeNumRows()) ? Teuchos::UNIT_DIAG : Teuchos::NON_UNIT_DIAG;
00182 
00183   lclMatOps_.template solve<DomainScalar,RangeScalar>(trans,triang,diag, *lclY,*lclX);
00184 }
00185 
00186 //-------------------------------------------------------------------
00187 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00188 void
00189 VbrMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::updateImport(const MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>& X) const
00190 {
00191   if (importer_ != Teuchos::null) {
00192     if (importedVec_ == Teuchos::null || importedVec_->getNumVectors() != X.getNumVectors()) {
00193       importedVec_ = Teuchos::rcp(new MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>(getBlockColMap()->getPointMap(), X.getNumVectors()));
00194     }
00195 
00196     importedVec_->doImport(X, *importer_, REPLACE);
00197   }
00198 }
00199 
00200 //-------------------------------------------------------------------
00201 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00202 void
00203 VbrMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::updateExport(const MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>& Y) const
00204 {
00205   if (exporter_ != Teuchos::null) {
00206     if (exportedVec_ == Teuchos::null || exportedVec_->getNumVectors() != Y.getNumVectors()) {
00207       exportedVec_ = Teuchos::rcp(new MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>(getBlockColMap()->getPointMap(), Y.getNumVectors()));
00208     }
00209 
00210     exportedVec_->doExport(Y, *exporter_, REPLACE);
00211   }
00212 }
00213 
00214 //-------------------------------------------------------------------
00215 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00216 void
00217 VbrMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::apply(
00218          const MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> &X,
00219                MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> &Y,
00220                Teuchos::ETransp trans,
00221                Scalar alpha,
00222                Scalar beta) const
00223 {
00224   if (trans == Teuchos::NO_TRANS) {
00225     updateImport(X);
00226     const MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>& Xref = importedVec_ == Teuchos::null ? X : *importedVec_;
00227     this->template multiply<Scalar,Scalar>(Xref, Y, trans, alpha, beta);
00228   }
00229   else if (trans == Teuchos::TRANS || trans == Teuchos::CONJ_TRANS) {
00230     updateImport(Y);
00231     MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>& Yref = importedVec_ == Teuchos::null ? Y : *importedVec_;
00232     this->template multiply<Scalar,Scalar>(X, Yref, trans, alpha, beta);
00233     if (importedVec_ != Teuchos::null) {
00234       Y.doExport(*importedVec_, *importer_, ADD);
00235     }
00236   }
00237 }
00238 
00239 //-------------------------------------------------------------------
00240 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00241 void
00242 VbrMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::applyInverse(
00243          const MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> &Y,
00244                MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> &X,
00245                Teuchos::ETransp trans) const
00246 {
00247   if (trans == Teuchos::NO_TRANS) {
00248     updateImport(Y);
00249     MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>& Xref = importedVec_ == Teuchos::null ? X : *importedVec_;
00250     this->template solve<Scalar,Scalar>(Y, Xref, trans);
00251   }
00252   else if (trans == Teuchos::TRANS || trans == Teuchos::CONJ_TRANS) {
00253     updateImport(Y);
00254     const MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>& Yref = importedVec_ == Teuchos::null ? Y : *importedVec_;
00255     this->template solve<Scalar,Scalar>(Yref, X, trans);
00256     if (importedVec_ != Teuchos::null) {
00257       X.doExport(*importedVec_, *importer_, ADD);
00258     }
00259   }
00260 }
00261 
00262 //-------------------------------------------------------------------
00263 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00264 bool
00265 VbrMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::hasTransposeApply() const
00266 {
00267   return false;
00268 }
00269 
00270 //-------------------------------------------------------------------
00271 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00272 const Teuchos::RCP<const BlockMap<LocalOrdinal,GlobalOrdinal,Node> > &
00273 VbrMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getBlockRowMap() const
00274 {
00275   return constBlkGraph_->getBlockRowMap();
00276 }
00277 
00278 //-------------------------------------------------------------------
00279 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00280 const Teuchos::RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > &
00281 VbrMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getPointRowMap() const
00282 {
00283   return constBlkGraph_->getBlockRowMap()->getPointMap();
00284 }
00285 
00286 //-------------------------------------------------------------------
00287 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00288 const Teuchos::RCP<const BlockMap<LocalOrdinal,GlobalOrdinal,Node> > &
00289 VbrMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getBlockColMap() const
00290 {
00291   return constBlkGraph_->getBlockColMap();
00292 }
00293 
00294 //-------------------------------------------------------------------
00295 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00296 const Teuchos::RCP<const BlockMap<LocalOrdinal,GlobalOrdinal,Node> > &
00297 VbrMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getBlockDomainMap() const
00298 {
00299   return constBlkGraph_->getBlockDomainMap();
00300 }
00301 
00302 //-------------------------------------------------------------------
00303 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00304 const Teuchos::RCP<const BlockMap<LocalOrdinal,GlobalOrdinal,Node> > &
00305 VbrMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getBlockRangeMap() const
00306 {
00307   return constBlkGraph_->getBlockRangeMap();
00308 }
00309 
00310 //-------------------------------------------------------------------
00311 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00312 const Teuchos::RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > &
00313 VbrMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getPointColMap() const
00314 {
00315   return constBlkGraph_->getBlockColMap()->getPointMap();
00316 }
00317 
00318 //-------------------------------------------------------------------
00319 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00320 void
00321 VbrMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getGlobalBlockEntryViewNonConst(
00322     GlobalOrdinal globalBlockRow,
00323     GlobalOrdinal globalBlockCol,
00324     LocalOrdinal& numPtRows,
00325     LocalOrdinal& numPtCols,
00326     Teuchos::ArrayRCP<Scalar>& blockEntry)
00327 {
00328   //Return a non-const, read-write view of a block-entry (as an ArrayRCP),
00329   //creating/allocating the block-entry if it doesn't already exist, (but only
00330   //if fillComplete hasn't been called yet, and if the arguments numPtRows
00331   //and numPtCols have been set on input).
00332 
00333   if (is_storage_optimized_) {
00334     TEST_FOR_EXCEPTION(!isFillComplete(), std::runtime_error,
00335       "Tpetra::VbrMatrix::getGlobalBlockEntryViewNonConst internal ERROR, storage is optimized but isFillComplete() is false.");
00336 
00337     LocalOrdinal localBlockRow = getBlockRowMap()->getLocalBlockID(globalBlockRow);
00338     LocalOrdinal localBlockCol = getBlockColMap()->getLocalBlockID(globalBlockCol);
00339 
00340     getLocalBlockEntryViewNonConst(localBlockRow, localBlockCol,
00341                                    numPtRows, numPtCols, blockEntry);
00342     return;
00343   }
00344 
00345   //If we get to here, fillComplete hasn't been called yet, and the matrix data
00346   //is stored in un-packed '2D' form.
00347 
00348   if (values2D_->size() == 0) {
00349     LocalOrdinal numBlockRows = constBlkGraph_->getNodeNumRows();
00350     values2D_->resize(numBlockRows);
00351     col_ind_2D_global_->resize(numBlockRows);
00352   }
00353 
00354   //this acts as a range-check for globalBlockRow:
00355   LocalOrdinal localBlockRow = getBlockRowMap()->getLocalBlockID(globalBlockRow);
00356   TEST_FOR_EXCEPTION( localBlockRow == Teuchos::OrdinalTraits<LocalOrdinal>::invalid(),
00357      std::runtime_error,
00358      "Tpetra::VbrMatrix::getGlobalBlockEntryViewNonConst, globalBlockRow not on the local processor.");
00359 
00360   MapGlobalArrayRCP& blkrow = (*col_ind_2D_global_)[localBlockRow];
00361   typename MapGlobalArrayRCP::iterator col_iter = blkrow.find(globalBlockCol);
00362 
00363   if (col_iter != blkrow.end()) {
00364     blockEntry = col_iter->second;
00365   }
00366   else {
00367     //blockEntry doesn't already exist, so we will create it.
00368 
00369     //make sure block-size is specified:
00370     TEST_FOR_EXCEPTION(numPtRows==0 || numPtCols==0, std::runtime_error,
00371       "Tpetra::VbrMatrix::getGlobalBlockEntryViewNonConst ERROR: creating block-entry, but numPtRows and/or numPtCols is 0.");
00372 
00373     Teuchos::RCP<Node> node = getNode();
00374     size_t blockSize = numPtRows*numPtCols;
00375     blockEntry = Teuchos::arcp(new Scalar[blockSize], 0, blockSize);
00376     std::fill(blockEntry.begin(), blockEntry.end(), 0);
00377     (*values2D_)[localBlockRow].push_back(blockEntry);
00378     blkrow.insert(std::make_pair(globalBlockCol, blockEntry));
00379     blkGraph_->insertGlobalIndices(globalBlockRow, Teuchos::ArrayView<GlobalOrdinal>(&globalBlockCol, 1));
00380   }
00381 }
00382 
00383 //-------------------------------------------------------------------
00384 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00385 void
00386 VbrMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getGlobalBlockEntryView(
00387       GlobalOrdinal globalBlockRow,
00388       GlobalOrdinal globalBlockCol,
00389       LocalOrdinal& numPtRows,
00390       LocalOrdinal& numPtCols,
00391       Teuchos::ArrayRCP<const Scalar>& blockEntry) const
00392 {
00393   //This method returns a const-view of a block-entry (as an ArrayRCP).
00394   //Throws an exception if the block-entry doesn't already exist.
00395 
00396   if (is_storage_optimized_) {
00397     TEST_FOR_EXCEPTION(!isFillComplete(), std::runtime_error,
00398       "Tpetra::VbrMatrix::getGlobalBlockEntryView internal ERROR, storage is optimized but isFillComplete() is false.");
00399 
00400     LocalOrdinal localBlockRow = getBlockRowMap()->getLocalBlockID(globalBlockRow);
00401     LocalOrdinal localBlockCol = getBlockColMap()->getLocalBlockID(globalBlockCol);
00402     getLocalBlockEntryView(localBlockRow, localBlockCol, numPtRows, numPtCols, blockEntry);
00403     return;
00404   }
00405 
00406   TEST_FOR_EXCEPTION(values2D_->size() == 0, std::runtime_error,
00407     "Tpetra::VbrMatrix::getGlobalBlockEntryView ERROR, matrix storage not yet allocated, can't return a const view.");
00408 
00409   //this acts as a range-check for globalBlockRow:
00410   LocalOrdinal localBlockRow = getBlockRowMap()->getLocalBlockID(globalBlockRow);
00411   TEST_FOR_EXCEPTION( localBlockRow == Teuchos::OrdinalTraits<LocalOrdinal>::invalid(),
00412      std::runtime_error,
00413      "Tpetra::VbrMatrix::getGlobalBlockEntryView, globalBlockRow not on the local processor.");
00414 
00415   MapGlobalArrayRCP& blkrow = (*col_ind_2D_global_)[localBlockRow];
00416   typename MapGlobalArrayRCP::iterator col_iter = blkrow.find(globalBlockCol);
00417 
00418   if (col_iter == blkrow.end()) {
00419     throw std::runtime_error("Tpetra::VbrMatrix::getGlobalBlockEntryView ERROR, specified block-entry doesn't exist.");
00420   }
00421 
00422   numPtRows = getBlockRowMap()->getLocalBlockSize(localBlockRow);
00423 
00424   TEST_FOR_EXCEPTION(numPtRows == 0, std::runtime_error,
00425     "Tpetra::VbrMatrix::getGlobalBlockEntryView ERROR, numPtRows == 0.");
00426 
00427   blockEntry = col_iter->second;
00428   numPtCols = blockEntry.size() / numPtRows;
00429 }
00430 
00431 //-------------------------------------------------------------------
00432 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00433 void
00434 VbrMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getLocalBlockEntryViewNonConst(
00435     LocalOrdinal localBlockRow,
00436     LocalOrdinal localBlockCol,
00437     LocalOrdinal& numPtRows,                     
00438     LocalOrdinal& numPtCols,
00439     Teuchos::ArrayRCP<Scalar>& blockEntry)
00440 {
00441   typedef Teuchos::ArrayRCP<const LocalOrdinal> Host_View_LO;
00442   typedef Teuchos::ArrayRCP<const size_t> Host_View;
00443   typedef typename Host_View_LO::iterator ITER;
00444   //This method returns a non-constant view of a block-entry (as an ArrayRCP).
00445 
00446   TEST_FOR_EXCEPTION(isFillComplete() == false, std::runtime_error,
00447    "Tpetra::VbrMatrix::getLocalBlockEntryViewNonConst ERROR, this method can only be called after fillComplete() has been called.");
00448 
00449   TEST_FOR_EXCEPTION(is_storage_optimized_ == false, std::runtime_error,
00450    "Tpetra::VbrMatrix::getLocalBlockEntryViewNonConst ERROR, this method can only be called if storage is optimized.");
00451 
00452   Teuchos::RCP<Node> node = getNode();
00453 
00454   Host_View bptr = constBlkGraph_->getNodeRowOffsets();
00455   LocalOrdinal bindx_offset = bptr[localBlockRow];
00456   LocalOrdinal length = bptr[localBlockRow+1] - bindx_offset;
00457 
00458   TEST_FOR_EXCEPTION( length < 1, std::runtime_error,
00459     "Tpetra::VbrMatrix::getLocalBlockEntryViewNonConst ERROR, specified localBlockCol not found in localBlockRow.");
00460 
00461   Host_View_LO bindx = constBlkGraph_->getNodePackedIndices();
00462   ITER bindx_beg = bindx.begin() + bindx_offset,  
00463        bindx_end = bindx_beg + length;
00464   ITER it = std::lower_bound(bindx_beg, bindx_end, localBlockCol);
00465 
00466   TEST_FOR_EXCEPTION(it == bindx_end || *it != localBlockCol, std::runtime_error,
00467     "Tpetra::VbrMatrix::getLocalBlockEntryViewNonConst ERROR, specified localBlockCol not found.");
00468 
00469   numPtRows = getBlockRowMap()->getLocalBlockSize(localBlockRow);
00470   numPtCols = getBlockColMap()->getLocalBlockSize(localBlockCol);
00471 
00472   const LocalOrdinal blkSize = numPtRows*numPtCols;
00473   Host_View_LO indx = node->template viewBuffer<LocalOrdinal>(1,pbuf_indx_+bptr[localBlockRow]+(it-bindx_beg));
00474   const LocalOrdinal offset = indx[0];
00475   blockEntry = node->template viewBufferNonConst<Scalar>(Kokkos::ReadWrite, blkSize, pbuf_values1D_ + offset);
00476 }
00477 
00478 //-------------------------------------------------------------------
00479 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00480 void
00481 VbrMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getLocalDiagCopy(
00482   Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node>& diag) const
00483 {
00484   const Teuchos::RCP<const BlockMap<LocalOrdinal,GlobalOrdinal,Node> >& rowmap = getBlockRowMap();
00485   TEST_FOR_EXCEPTION(diag.getMap()->isSameAs(*(rowmap->getPointMap())) != true,
00486     std::runtime_error, "Tpetra::VbrMatrix::getLocalDiagCopy ERROR, vector must be distributed the same as this matrix' row-map.");
00487 
00488   Teuchos::ArrayRCP<Scalar> diag_view = diag.get1dViewNonConst();
00489   Teuchos::ArrayView<const GlobalOrdinal> blockIDs = rowmap->getNodeBlockIDs();
00490   size_t offset = 0;
00491   typedef typename Teuchos::ArrayView<const GlobalOrdinal>::size_type Tsize_t;
00492   for(Tsize_t i=0; i<blockIDs.size(); ++i) {
00493     LocalOrdinal localBlockID = rowmap->getLocalBlockID(blockIDs[i]);
00494     LocalOrdinal blockSize = rowmap->getLocalBlockSize(localBlockID);
00495     Teuchos::ArrayRCP<const Scalar> blockEntry;
00496     getGlobalBlockEntryView(blockIDs[i], blockIDs[i], blockSize, blockSize, blockEntry);
00497 
00498     for(LocalOrdinal j=0; j<blockSize; ++j) {
00499       diag_view[offset++] = blockEntry[j*blockSize+j];
00500     }
00501   }
00502 }
00503 
00504 //-------------------------------------------------------------------
00505 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00506 void
00507 VbrMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getLocalBlockEntryView(
00508       LocalOrdinal localBlockRow,
00509       LocalOrdinal localBlockCol,
00510       LocalOrdinal& numPtRows,
00511       LocalOrdinal& numPtCols,
00512       Teuchos::ArrayRCP<const Scalar>& blockEntry) const
00513 {
00514   typedef Teuchos::ArrayRCP<const size_t>       Host_View;
00515   typedef Teuchos::ArrayRCP<const LocalOrdinal> Host_View_LO;
00516   typedef typename Host_View_LO::iterator ITER;
00517   //This method returns a constant view of a block-entry (as an ArrayRCP).
00518 
00519   TEST_FOR_EXCEPTION(isFillComplete() == false, std::runtime_error,
00520    "Tpetra::VbrMatrix::getLocalBlockEntryView ERROR, this method can only be called after fillComplete() has been called.");
00521 
00522   TEST_FOR_EXCEPTION(is_storage_optimized_ == false, std::runtime_error,
00523    "Tpetra::VbrMatrix::getLocalBlockEntryView ERROR, this method can only be called if storage is optimized.");
00524 
00525   Teuchos::RCP<Node> node = getNode();
00526 
00527   Host_View bpr = constBlkGraph_->getNodeRowOffsets();
00528   const size_t bindx_offset = bpr[localBlockRow];
00529   const size_t length       = bpr[localBlockRow+1] - bindx_offset;
00530 
00531   Host_View_LO bindx = constBlkGraph_->getNodePackedIndices();
00532   ITER bindx_beg = bindx.begin()+bindx_offset, 
00533        bindx_end = bindx_beg + length;
00534   ITER it = std::lower_bound(bindx_beg, bindx_end, localBlockCol);
00535 
00536   TEST_FOR_EXCEPTION(it == bindx_end || *it != localBlockCol, std::runtime_error,
00537     "Tpetra::VbrMatrix::getLocalBlockEntryView ERROR, specified localBlockCol not found.");
00538 
00539   numPtRows = getBlockRowMap()->getLocalBlockSize(localBlockRow);
00540   numPtCols = getBlockColMap()->getLocalBlockSize(localBlockCol);
00541 
00542   const LocalOrdinal blkSize = numPtRows*numPtCols;
00543   Host_View_LO indx = node->template viewBuffer<LocalOrdinal>(1, pbuf_indx_+bpr[localBlockRow]+(it-bindx_beg));
00544   const LocalOrdinal offset = indx[0];
00545   blockEntry = node->template viewBuffer<Scalar>(blkSize, pbuf_values1D_ + offset);
00546 }
00547 
00548 //-------------------------------------------------------------------
00549 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00550 Teuchos::RCP<Node>
00551 VbrMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::getNode() const
00552 {
00553   return getBlockRowMap()->getPointMap()->getNode();
00554 }
00555 
00556 //-------------------------------------------------------------------
00557 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00558 void
00559 VbrMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::putScalar(Scalar s)
00560 {
00561   Teuchos::RCP<Node> node = getNode();
00562 
00563   if (isFillComplete()) {
00564     fill_device_ArrayRCP(node, pbuf_values1D_, s);
00565   }
00566   else {
00567     typedef typename Teuchos::Array<Teuchos::Array<Teuchos::ArrayRCP<Scalar> > >::size_type Tsize_t;
00568     Teuchos::Array<Teuchos::Array<Teuchos::ArrayRCP<Scalar> > >& vals2D = *values2D_;
00569     LocalOrdinal numBlockRows = constBlkGraph_->getNodeNumRows();
00570     for(LocalOrdinal r=0; r<numBlockRows; ++r) {
00571       for(Tsize_t c=0; c<vals2D[r].size(); ++c) {
00572         std::fill(vals2D[r][c].begin(), vals2D[r][c].end(), s);
00573       }
00574     }
00575   }
00576 }
00577 
00578 //-------------------------------------------------------------------
00579 template<class ArrayType1,class ArrayType2,class Ordinal>
00580 void set_array_values(ArrayType1& dest, ArrayType2& src, Ordinal rows, Ordinal cols, Ordinal stride, Tpetra::CombineMode mode)
00581 {
00582   size_t offset = 0;
00583   size_t src_offset = 0;
00584 
00585   if (mode == ADD) {
00586     for(Ordinal col=0; col<cols; ++col) {
00587       for(Ordinal row=0; row<rows; ++row) {
00588         dest[offset++] += src[src_offset + row];
00589       }
00590       src_offset += stride;
00591     }
00592   }
00593   else {
00594     for(Ordinal col=0; col<cols; ++col) {
00595       for(Ordinal row=0; row<rows; ++row) {
00596         dest[offset++] = src[src_offset + row];
00597       }
00598       src_offset += stride;
00599     }
00600   }
00601 }
00602 
00603 //-------------------------------------------------------------------
00604 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00605 void
00606 VbrMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::setGlobalBlockEntry(GlobalOrdinal globalBlockRow, GlobalOrdinal globalBlockCol, const Teuchos::SerialDenseMatrix<GlobalOrdinal,Scalar>& blockEntry)
00607 {
00608   //first get an ArrayRCP for the internal storage for this block-entry:
00609   Teuchos::ArrayRCP<Scalar> internalBlockEntry;
00610   LocalOrdinal blkRowSize = blockEntry.numRows();
00611   LocalOrdinal blkColSize = blockEntry.numCols();
00612   getGlobalBlockEntryViewNonConst(globalBlockRow,globalBlockCol, blkRowSize, blkColSize, internalBlockEntry);
00613 
00614   //now copy the incoming block-entry into internal storage:
00615   const Scalar* inputvalues = blockEntry.values();
00616   set_array_values(internalBlockEntry, inputvalues, blkRowSize, blkColSize, blockEntry.stride(), REPLACE);
00617 }
00618 
00619 //-------------------------------------------------------------------
00620 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00621 void
00622 VbrMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::setLocalBlockEntry(LocalOrdinal localBlockRow, LocalOrdinal localBlockCol, const Teuchos::SerialDenseMatrix<LocalOrdinal,Scalar>& blockEntry)
00623 {
00624   //first get an ArrayRCP for the internal storage for this block-entry:
00625   Teuchos::ArrayRCP<Scalar> internalBlockEntry;
00626   LocalOrdinal blkRowSize = blockEntry.numRows();
00627   LocalOrdinal blkColSize = blockEntry.numCols();
00628   getLocalBlockEntryViewNonConst(localBlockRow,localBlockCol, blkRowSize, blkColSize, internalBlockEntry);
00629 
00630   //now copy the incoming block-entry into internal storage:
00631   const Scalar* inputvalues = blockEntry.values();
00632   set_array_values(internalBlockEntry, inputvalues, blkRowSize, blkColSize, blockEntry.stride(), REPLACE);
00633 }
00634 
00635 //-------------------------------------------------------------------
00636 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00637 void
00638 VbrMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::sumIntoGlobalBlockEntry(GlobalOrdinal globalBlockRow, GlobalOrdinal globalBlockCol, const Teuchos::SerialDenseMatrix<GlobalOrdinal,Scalar>& blockEntry)
00639 {
00640   //first get an ArrayRCP for the internal storage for this block-entry:
00641   Teuchos::ArrayRCP<Scalar> internalBlockEntry;
00642   LocalOrdinal blkRowSize = blockEntry.numRows();
00643   LocalOrdinal blkColSize = blockEntry.numCols();
00644   getGlobalBlockEntryViewNonConst(globalBlockRow,globalBlockCol, blkRowSize, blkColSize, internalBlockEntry);
00645 
00646   //now sum (add) the incoming block-entry into internal storage:
00647   const Scalar* inputvalues = blockEntry.values();
00648   set_array_values(internalBlockEntry, inputvalues, blkRowSize, blkColSize, blockEntry.stride(), ADD);
00649 }
00650 
00651 //-------------------------------------------------------------------
00652 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00653 void
00654 VbrMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::sumIntoLocalBlockEntry(LocalOrdinal localBlockRow, LocalOrdinal localBlockCol, const Teuchos::SerialDenseMatrix<LocalOrdinal,Scalar>& blockEntry)
00655 {
00656   //first get an ArrayRCP for the internal storage for this block-entry:
00657   Teuchos::ArrayRCP<Scalar> internalBlockEntry;
00658   LocalOrdinal blkRowSize = blockEntry.numRows();
00659   LocalOrdinal blkColSize = blockEntry.numCols();
00660   getLocalBlockEntryViewNonConst(localBlockRow,localBlockCol, blkRowSize, blkColSize, internalBlockEntry);
00661 
00662   //now sum (add) the incoming block-entry into internal storage:
00663   const Scalar* inputvalues = blockEntry.values();
00664   set_array_values(internalBlockEntry, inputvalues, blkRowSize, blkColSize, blockEntry.stride(), ADD);
00665 }
00666 
00667 //-------------------------------------------------------------------
00668 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00669 void
00670 VbrMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::setGlobalBlockEntry(GlobalOrdinal globalBlockRow, GlobalOrdinal globalBlockCol, LocalOrdinal blkRowSize, LocalOrdinal blkColSize, LocalOrdinal LDA, const Teuchos::ArrayView<const Scalar>& blockEntry)
00671 {
00672   //first get an ArrayRCP for the internal storage for this block-entry:
00673   Teuchos::ArrayRCP<Scalar> internalBlockEntry;
00674   getGlobalBlockEntryViewNonConst(globalBlockRow,globalBlockCol, blkRowSize, blkColSize, internalBlockEntry);
00675 
00676   LocalOrdinal blk_size = blockEntry.size();
00677   TEST_FOR_EXCEPTION(blkColSize*LDA > blk_size, std::runtime_error,
00678     "Tpetra::VbrMatrix::setGlobalBlockEntry ERROR, inconsistent block-entry sizes.");
00679 
00680   //copy the incoming block-entry into internal storage:
00681   set_array_values(internalBlockEntry, blockEntry, blkRowSize, blkColSize, LDA, REPLACE);
00682 }
00683 
00684 //-------------------------------------------------------------------
00685 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00686 void
00687 VbrMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::setLocalBlockEntry(LocalOrdinal localBlockRow, LocalOrdinal localBlockCol, LocalOrdinal blkRowSize, LocalOrdinal blkColSize, LocalOrdinal LDA, const Teuchos::ArrayView<const Scalar>& blockEntry)
00688 {
00689   //first get an ArrayRCP for the internal storage for this block-entry:
00690   Teuchos::ArrayRCP<Scalar> internalBlockEntry;
00691   getLocalBlockEntryViewNonConst(localBlockRow,localBlockCol, blkRowSize, blkColSize, internalBlockEntry);
00692 
00693   LocalOrdinal blk_size = blockEntry.size();
00694   TEST_FOR_EXCEPTION(blkColSize*LDA > blk_size, std::runtime_error,
00695     "Tpetra::VbrMatrix::setLocalBlockEntry ERROR, inconsistent block-entry sizes.");
00696 
00697   //copy the incoming block-entry into internal storage:
00698   set_array_values(internalBlockEntry, blockEntry, blkRowSize, blkColSize, LDA, REPLACE);
00699 }
00700 
00701 //-------------------------------------------------------------------
00702 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00703 void
00704 VbrMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::sumIntoGlobalBlockEntry(GlobalOrdinal globalBlockRow, GlobalOrdinal globalBlockCol, LocalOrdinal blkRowSize, LocalOrdinal blkColSize, LocalOrdinal LDA, const Teuchos::ArrayView<const Scalar>& blockEntry)
00705 {
00706   //first get an ArrayRCP for the internal storage for this block-entry:
00707   Teuchos::ArrayRCP<Scalar> internalBlockEntry;
00708   getGlobalBlockEntryViewNonConst(globalBlockRow,globalBlockCol, blkRowSize, blkColSize, internalBlockEntry);
00709 
00710   LocalOrdinal blk_size = blockEntry.size();
00711   TEST_FOR_EXCEPTION(blkColSize*LDA > blk_size, std::runtime_error,
00712     "Tpetra::VbrMatrix::setGlobalBlockEntry ERROR, inconsistent block-entry sizes.");
00713 
00714   //copy the incoming block-entry into internal storage:
00715   set_array_values(internalBlockEntry, blockEntry, blkRowSize, blkColSize, LDA, ADD);
00716 }
00717 
00718 //-------------------------------------------------------------------
00719 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00720 void
00721 VbrMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::sumIntoLocalBlockEntry(LocalOrdinal localBlockRow, LocalOrdinal localBlockCol, LocalOrdinal blkRowSize, LocalOrdinal blkColSize, LocalOrdinal LDA, const Teuchos::ArrayView<const Scalar>& blockEntry)
00722 {
00723   //first get an ArrayRCP for the internal storage for this block-entry:
00724   Teuchos::ArrayRCP<Scalar> internalBlockEntry;
00725   getLocalBlockEntryViewNonConst(localBlockRow,localBlockCol, blkRowSize, blkColSize, internalBlockEntry);
00726 
00727   LocalOrdinal blk_size = blockEntry.size();
00728   TEST_FOR_EXCEPTION(blkColSize*LDA > blk_size, std::runtime_error,
00729     "Tpetra::VbrMatrix::setLocalBlockEntry ERROR, inconsistent block-entry sizes.");
00730 
00731   //copy the incoming block-entry into internal storage:
00732   set_array_values(internalBlockEntry, blockEntry, blkRowSize, blkColSize, LDA, ADD);
00733 }
00734 
00735 //-------------------------------------------------------------------
00736 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00737 void
00738 VbrMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::optimizeStorage()
00739 {
00740   typedef Teuchos::ArrayRCP<const size_t> Host_View;
00741   typedef Teuchos::ArrayRCP<const LocalOrdinal> Host_View_LO;
00742   typedef typename Teuchos::ArrayRCP<size_t>::size_type Tsize_t;
00743 
00744   if (is_storage_optimized_ == true) return;
00745 
00746   TEST_FOR_EXCEPTION(constBlkGraph_->isFillComplete() != true, std::runtime_error,
00747     "Tpetra::VbrMatrix::optimizeStorage ERROR, isFillComplete() is false, required to be true before optimizeStorage is called.");
00748 
00749   size_t num_block_nonzeros = constBlkGraph_->getNodeNumEntries();
00750 
00751   const Teuchos::RCP<const BlockMap<LocalOrdinal,GlobalOrdinal,Node> >& rowmap = constBlkGraph_->getBlockRowMap();
00752   const Teuchos::RCP<const BlockMap<LocalOrdinal,GlobalOrdinal,Node> >& colmap = constBlkGraph_->getBlockColMap();
00753 
00754   const Teuchos::RCP<Node>& node = getBlockRowMap()->getPointMap()->getNode();
00755 
00756   //need to count the number of point-entries:
00757   size_t num_point_entries = 0;
00758   Host_View bptr = constBlkGraph_->getNodeRowOffsets();
00759 
00760   if (bptr.size() == 0) {
00761     is_storage_optimized_ = true;
00762     return;
00763   }
00764 
00765   Host_View_LO bindx = constBlkGraph_->getNodePackedIndices();
00766 
00767   for(Tsize_t r=0; r<bptr.size()-1; ++r) {
00768     size_t rbeg = bptr[r];
00769     size_t rend = bptr[r+1];
00770 
00771     LocalOrdinal rsize = rowmap->getLocalBlockSize(r);
00772 
00773     for(size_t j=rbeg; j<rend; ++j) {
00774       LocalOrdinal csize = colmap->getLocalBlockSize(bindx[j]);
00775       num_point_entries += rsize*csize;
00776     }
00777   }
00778 
00779   pbuf_indx_ = node->template allocBuffer<LocalOrdinal>(num_block_nonzeros+1);
00780   pbuf_values1D_ = node->template allocBuffer<Scalar>(num_point_entries);
00781 
00782   Teuchos::ArrayRCP<LocalOrdinal> view_indx = node->template viewBufferNonConst<LocalOrdinal>(Kokkos::WriteOnly, num_block_nonzeros+1, pbuf_indx_);
00783   Teuchos::ArrayRCP<Scalar> view_values1D = node->template viewBufferNonConst<Scalar>(Kokkos::WriteOnly, num_point_entries, pbuf_values1D_);
00784 
00785   bool have_2D_data = col_ind_2D_global_->size() > 0;
00786   Scalar zero = Teuchos::ScalarTraits<Scalar>::zero();
00787 
00788   size_t ioffset = 0;
00789   size_t offset = 0;
00790   for(Tsize_t r=0; r<bptr.size()-1; ++r) {
00791     LocalOrdinal rsize = rowmap->getLocalBlockSize(r);
00792 
00793     MapGlobalArrayRCP* blk_row = have_2D_data ? &((*col_ind_2D_global_)[r]) : NULL;
00794 
00795     for(size_t c=bptr[r]; c<bptr[r+1]; ++c) {
00796       view_indx[ioffset++] = offset;
00797 
00798       LocalOrdinal csize = colmap->getLocalBlockSize(bindx[c]);
00799       Tsize_t blkSize = rsize*csize;
00800 
00801       if (blk_row != NULL) {
00802         GlobalOrdinal global_col = colmap->getGlobalBlockID(bindx[c]);
00803         typename MapGlobalArrayRCP::iterator iter = blk_row->find(global_col);
00804         TEST_FOR_EXCEPTION(iter == blk_row->end(), std::runtime_error,
00805           "Tpetra::VbrMatrix::optimizeStorage ERROR, global_col not found in row.");
00806   
00807         Teuchos::ArrayRCP<Scalar> vals = iter->second;
00808         for(Tsize_t n=0; n<blkSize; ++n) {
00809           view_values1D[offset+n] = vals[n];
00810         }
00811       }
00812       else {
00813         for(Tsize_t n=0; n<blkSize; ++n) view_values1D[offset+n] = zero;
00814       }
00815       offset += blkSize;
00816     }
00817   }
00818   view_indx[ioffset] = offset;
00819 
00820   //Final step: release memory for the "2D" storage:
00821   col_ind_2D_global_ = Teuchos::null;
00822   values2D_ = Teuchos::null;
00823 
00824   is_storage_optimized_ = true;
00825 }
00826 
00827 //-------------------------------------------------------------------
00828 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00829 void VbrMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::fillLocalMatrix()
00830 {
00831   //We insist that optimzeStorage has already been called.
00832   //We don't care whether this function (fillLocalMatrix()) is being
00833   //called for the first time or not.
00834   TEST_FOR_EXCEPTION(is_storage_optimized_ != true, std::runtime_error,
00835     "Tpetra::VbrMatrix::fillLocalMatrix ERROR, optimizeStorage is required to have already been called.");
00836 
00837   Teuchos::ArrayRCP<const size_t      > nodeRowOffsets = constBlkGraph_->getNodeRowOffsets();
00838   Teuchos::ArrayRCP<const LocalOrdinal> nodePackedInds = constBlkGraph_->getNodePackedIndices();
00839   if (Node::isHostNode == false) {
00840     RCP<Node> node = getRangeMap()->getNode();
00841     //
00842     Teuchos::ArrayRCP<size_t> dev_nodeRowOffsets = node->template allocBuffer<size_t>(nodeRowOffsets.size());
00843     node->copyToBuffer(nodeRowOffsets.size(), nodeRowOffsets(), dev_nodeRowOffsets);
00844     nodeRowOffsets = dev_nodeRowOffsets;
00845     //
00846     Teuchos::ArrayRCP<LocalOrdinal> dev_nodePackedInds = node->template allocBuffer<LocalOrdinal>(nodePackedInds.size());
00847     node->copyToBuffer(nodePackedInds.size(), nodePackedInds(), dev_nodePackedInds);
00848     nodePackedInds = dev_nodePackedInds;
00849   }
00850   lclMatrix_.setPackedValues(pbuf_values1D_,
00851                              getBlockRowMap()->getNodeFirstPointInBlocks_Device(),
00852                              getBlockColMap()->getNodeFirstPointInBlocks_Device(),
00853                              nodeRowOffsets,
00854                              nodePackedInds,
00855                              pbuf_indx_);
00856 }
00857 
00858 //-------------------------------------------------------------------
00859 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00860 void VbrMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::fillLocalMatVec()
00861 {
00862   //We insist that optimzeStorage has already been called.
00863   //We don't care whether this function (fillLocalMatVec()) is being
00864   //called for the first time or not.
00865   TEST_FOR_EXCEPTION(is_storage_optimized_ != true, std::runtime_error,
00866     "Tpetra::VbrMatrix::fillLocalMatrix ERROR, optimizeStorage is required to have already been called.");
00867 
00868   lclMatOps_.initializeValues(lclMatrix_);
00869 }
00870 
00871 //-------------------------------------------------------------------
00872 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00873 void
00874 VbrMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::createImporterExporter()
00875 {
00876   typedef typename Teuchos::RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > PtMap;
00877 
00878   const PtMap& ptDomMap = (constBlkGraph_->getBlockDomainMap())->getPointMap();
00879   const PtMap& ptRngMap = (constBlkGraph_->getBlockRangeMap())->getPointMap();
00880   const PtMap& ptRowMap = (constBlkGraph_->getBlockRowMap())->getPointMap();
00881   const PtMap& ptColMap = (constBlkGraph_->getBlockColMap())->getPointMap();
00882 
00883   if (!ptDomMap->isSameAs(*ptColMap)) {
00884     importer_ = Teuchos::rcp(new Tpetra::Import<LocalOrdinal,GlobalOrdinal,Node>(ptDomMap, ptColMap));
00885   }
00886   if (!ptRngMap->isSameAs(*ptRowMap)) {
00887     exporter_ = Teuchos::rcp(new Tpetra::Export<LocalOrdinal,GlobalOrdinal,Node>(ptRowMap, ptRngMap));
00888   }
00889 }
00890 
00891 //-------------------------------------------------------------------
00892 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00893 void
00894 VbrMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::fillComplete(const Teuchos::RCP<const BlockMap<LocalOrdinal,GlobalOrdinal,Node> >& blockDomainMap, const Teuchos::RCP<const BlockMap<LocalOrdinal,GlobalOrdinal,Node> >& blockRangeMap)
00895 {
00896   if (isFillComplete()) return;
00897 
00898   blkGraph_->fillComplete(blockDomainMap, blockRangeMap);
00899 
00900   createImporterExporter();
00901 
00902   is_fill_completed_ = true;
00903 
00904   optimizeStorage();
00905 
00906   fillLocalMatrix();
00907   fillLocalMatVec();
00908 }
00909 
00910 //-------------------------------------------------------------------
00911 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00912 void VbrMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::fillComplete()
00913 {
00914   //In this case, our block-row-map will also be our domain-map and
00915   //range-map.
00916 
00917   fillComplete(getBlockRowMap(), getBlockRowMap());
00918 }
00919 
00920 //-------------------------------------------------------------------
00921 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00922 std::string
00923 VbrMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::description() const
00924 {
00925   std::ostringstream oss;
00926   oss << Teuchos::Describable::description();
00927   if (isFillComplete()) {
00928     oss << "{status = fill complete, global num block rows = " << getBlockRowMap()->getGlobalNumBlocks() << "}";
00929   }
00930   else {
00931     oss << "{status = fill not complete, global num block rows = " << getBlockRowMap()->getGlobalNumBlocks() << "}";
00932   }
00933   return oss.str();
00934 }
00935 
00936 //-------------------------------------------------------------------
00937 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatOps>
00938 void
00939 VbrMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>::describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel) const
00940 {
00941   Teuchos::EVerbosityLevel vl = verbLevel;
00942   if (vl == Teuchos::VERB_DEFAULT) vl = Teuchos::VERB_LOW;
00943 
00944   if (vl == Teuchos::VERB_NONE) return;
00945 
00946   Teuchos::RCP<const Teuchos::Comm<int> > comm = getBlockRowMap()->getPointMap()->getComm();
00947   const int myProc = comm->getRank();
00948 
00949   if (myProc == 0) out << "VbrMatrix::describe is under construction..." << std::endl;
00950 }
00951 
00952 }//namespace Tpetra
00953 
00954 //
00955 // Explicit instantiation macro
00956 //
00957 // Must be expanded from within the Tpetra namespace!
00958 //
00959 
00960 #define TPETRA_VBRMATRIX_INSTANT(SCALAR,LO,GO,NODE) \
00961   \
00962   template class VbrMatrix< SCALAR , LO , GO , NODE >;
00963 
00964 #endif //TPETRA_VBRMATRIX_DEF_HPP
00965 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines