Tpetra Matrix/Vector Services Version of the Day
Tpetra_VbrMatrix_decl.hpp
Go to the documentation of this file.
00001 // @HEADER
00002 // ***********************************************************************
00003 //
00004 //          Tpetra: Templated Linear Algebra Services Package
00005 //                 Copyright (2008) Sandia Corporation
00006 //
00007 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
00008 // the U.S. Government retains certain rights in this software.
00009 //
00010 // Redistribution and use in source and binary forms, with or without
00011 // modification, are permitted provided that the following conditions are
00012 // met:
00013 //
00014 // 1. Redistributions of source code must retain the above copyright
00015 // notice, this list of conditions and the following disclaimer.
00016 //
00017 // 2. Redistributions in binary form must reproduce the above copyright
00018 // notice, this list of conditions and the following disclaimer in the
00019 // documentation and/or other materials provided with the distribution.
00020 //
00021 // 3. Neither the name of the Corporation nor the names of the
00022 // contributors may be used to endorse or promote products derived from
00023 // this software without specific prior written permission.
00024 //
00025 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00026 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00027 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00028 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00029 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00030 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00031 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00032 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00033 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00034 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00035 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00036 //
00037 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
00038 //
00039 // ************************************************************************
00040 // @HEADER
00041 
00042 #ifndef TPETRA_VBRMATRIX_DECL_HPP
00043 #define TPETRA_VBRMATRIX_DECL_HPP
00044 
00047 
00048 #include <Tpetra_ConfigDefs.hpp>
00049 #include <Tpetra_DistObject_decl.hpp>
00050 #include <Tpetra_Operator.hpp>
00051 #include <Tpetra_VbrUtils.hpp>
00052 #include <Kokkos_DefaultNode.hpp>
00053 #include <Kokkos_DefaultKernels.hpp>
00054 
00055 // Forward declarations, to avoid including files that we don't really
00056 // need just for declaring VbrMatrix.  The #ifndef ... #endif prevents
00057 // Doxygen from generating spurious documentation for the forward
00058 // declarations.
00059 #ifndef DOXYGEN_SHOULD_SKIP_THIS
00060 namespace Tpetra {
00061   template<class LocalOrdinal, class GlobalOrdinal, class Node>
00062   class BlockMap;
00063 
00064   template<class LocalOrdinal, class GlobalOrdinal, class Node>
00065   class BlockCrsGraph;
00066 
00067   template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
00068   class Vector;
00069 } // namespace Tpetra
00070 
00071 namespace KokkosClassic {
00072   template<class Scalar, class LocalOrdinal, class Node>
00073   class VbrMatrix;
00074 } // namespace KokkosClassic
00075 
00076 namespace Teuchos {
00077   template<class OrdinalType, class ScalarType>
00078   class SerialDenseMatrix;
00079 } // namespace Teuchos
00080 #endif // DOXYGEN_SHOULD_SKIP_THIS
00081 
00082 namespace Tpetra {
00083 
00085 
00113 template <class Scalar,
00114           class LocalOrdinal  = int,
00115           class GlobalOrdinal = LocalOrdinal,
00116           class Node          = KokkosClassic::DefaultNode::DefaultNodeType,
00117           class LocalMatOps   = typename KokkosClassic::DefaultKernels<Scalar,LocalOrdinal,Node>::BlockSparseOps >
00118 class VbrMatrix :
00119     public Tpetra::DistObject<char, LocalOrdinal, GlobalOrdinal, Node>,
00120     public Tpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node> {
00121 public:
00122   typedef Scalar        scalar_type;
00123   typedef LocalOrdinal  local_ordinal_type;
00124   typedef GlobalOrdinal global_ordinal_type;
00125   typedef Node          node_type;
00126   typedef LocalMatOps   mat_vec_type;
00127 
00129 
00130 
00132 
00137   VbrMatrix (const Teuchos::RCP<const BlockMap<LocalOrdinal,GlobalOrdinal,Node> >& blkRowMap, 
00138        size_t maxNumEntriesPerRow, 
00139        ProfileType pftype = DynamicProfile);
00140 
00142 
00152   VbrMatrix (const Teuchos::RCP<const BlockCrsGraph<LocalOrdinal,GlobalOrdinal,Node> >& blkGraph);
00153 
00155   virtual ~VbrMatrix();
00156 
00158 
00159 
00160 
00162 
00167   template <class DomainScalar, class RangeScalar>
00168   void 
00169   multiply (const MultiVector<DomainScalar,LocalOrdinal,GlobalOrdinal,Node>& X, 
00170       MultiVector<RangeScalar,LocalOrdinal,GlobalOrdinal,Node>& Y, 
00171       Teuchos::ETransp trans, 
00172       RangeScalar alpha, 
00173       RangeScalar beta) const;
00174 
00176 
00188   template <class DomainScalar, class RangeScalar>
00189   void 
00190   solve (const MultiVector<RangeScalar,LocalOrdinal,GlobalOrdinal,Node>& Y, 
00191    MultiVector<DomainScalar,LocalOrdinal,GlobalOrdinal,Node>& X, 
00192    Teuchos::ETransp trans) const;
00193 
00195 
00196 
00197 
00199 
00201   Teuchos::RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > getDomainMap() const;
00202 
00204 
00206   Teuchos::RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > getRangeMap() const;
00207 
00209 
00214   void
00215   apply (const MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> &X,
00216    MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> &Y,
00217    Teuchos::ETransp trans = Teuchos::NO_TRANS,
00218    Scalar alpha = Teuchos::ScalarTraits<Scalar>::one(),
00219    Scalar beta = Teuchos::ScalarTraits<Scalar>::zero()) const;
00220 
00222 
00225   void
00226   applyInverse (const MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> & Y,
00227     MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> &X,
00228     Teuchos::ETransp trans) const;
00229 
00231 
00234   bool hasTransposeApply() const;
00235 
00237 
00239 
00240 
00242   Teuchos::RCP<const BlockMap<LocalOrdinal,GlobalOrdinal,Node> > getBlockRowMap() const;
00243 
00245   Teuchos::RCP<const BlockMap<LocalOrdinal,GlobalOrdinal,Node> > getBlockColMap() const;
00246 
00248   Teuchos::RCP<const BlockMap<LocalOrdinal,GlobalOrdinal,Node> > getBlockDomainMap() const;
00249 
00251   Teuchos::RCP<const BlockMap<LocalOrdinal,GlobalOrdinal,Node> > getBlockRangeMap() const;
00252 
00254   Teuchos::RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > getPointRowMap() const;
00255 
00257   Teuchos::RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > getPointColMap() const;
00258 
00260   bool isFillComplete() const;
00261 
00263 
00264 
00265 
00267 
00270   void putScalar(Scalar s);
00271 
00273 
00288   void setGlobalBlockEntry(GlobalOrdinal globalBlockRow, GlobalOrdinal globalBlockCol, const Teuchos::SerialDenseMatrix<int,Scalar>& blockEntry);
00289 
00291 
00298   void setLocalBlockEntry(LocalOrdinal localBlockRow, LocalOrdinal localBlockCol, const Teuchos::SerialDenseMatrix<int,Scalar>& blockEntry);
00299 
00301 
00316   void sumIntoGlobalBlockEntry(GlobalOrdinal globalBlockRow, GlobalOrdinal globalBlockCol, const Teuchos::SerialDenseMatrix<int,Scalar>& blockEntry);
00317 
00319 
00326   void sumIntoLocalBlockEntry(LocalOrdinal localBlockRow, LocalOrdinal localBlockCol, const Teuchos::SerialDenseMatrix<int,Scalar>& blockEntry);
00327 
00329 
00344   void setGlobalBlockEntry(GlobalOrdinal globalBlockRow, GlobalOrdinal globalBlockCol, LocalOrdinal blkRowSize, LocalOrdinal blkColSize, LocalOrdinal LDA, const Teuchos::ArrayView<const Scalar>& blockEntry);
00345 
00347 
00354   void setLocalBlockEntry(LocalOrdinal localBlockRow, LocalOrdinal localBlockCol, LocalOrdinal blkRowSize, LocalOrdinal blkColSize, LocalOrdinal LDA, const Teuchos::ArrayView<const Scalar>& blockEntry);
00355 
00357 
00372   void sumIntoGlobalBlockEntry(GlobalOrdinal globalBlockRow, GlobalOrdinal globalBlockCol, LocalOrdinal blkRowSize, LocalOrdinal blkColSize, LocalOrdinal LDA, const Teuchos::ArrayView<const Scalar>& blockEntry);
00373 
00375 
00382   void sumIntoLocalBlockEntry(LocalOrdinal localBlockRow, LocalOrdinal localBlockCol, LocalOrdinal blkRowSize, LocalOrdinal blkColSize, LocalOrdinal LDA, const Teuchos::ArrayView<const Scalar>& blockEntry);
00383 
00385 
00387 
00388 
00390 
00394   void fillComplete(const Teuchos::RCP<const BlockMap<LocalOrdinal,GlobalOrdinal,Node> >& blockDomainMap, const Teuchos::RCP<const BlockMap<LocalOrdinal,GlobalOrdinal,Node> >& blockRangeMap);
00395 
00397 
00400   void fillComplete();
00401 
00403   void globalAssemble();
00404 
00406 
00407 
00408 
00410 
00412   void getGlobalBlockRowView(GlobalOrdinal globalBlockRow,
00413                              LocalOrdinal& numPtRows,
00414                              Teuchos::ArrayView<const GlobalOrdinal>& blockCols,
00415                              Teuchos::Array<LocalOrdinal>& ptColsPerBlockCol,
00416                              Teuchos::Array<Teuchos::ArrayRCP<const Scalar> >& blockEntries) const;
00417 
00419 
00421   void getLocalBlockRowView(LocalOrdinal localBlockRow,
00422                             LocalOrdinal& numPtRows,
00423                             Teuchos::ArrayView<const LocalOrdinal>& blockCols,
00424                             Teuchos::Array<LocalOrdinal>& ptColsPerBlockCol,
00425                             Teuchos::ArrayRCP<const Scalar>& blockEntries) const;
00426 
00428 
00436   void getGlobalBlockEntryView(GlobalOrdinal globalBlockRow,
00437                                GlobalOrdinal globalBlockCol,
00438                                LocalOrdinal& numPtRows,
00439                                LocalOrdinal& numPtCols,
00440                                Teuchos::ArrayRCP<const Scalar>& blockEntry) const;
00441 
00443 
00453   void getGlobalBlockEntryViewNonConst(GlobalOrdinal globalBlockRow,
00454                                        GlobalOrdinal globalBlockCol,
00455                                        LocalOrdinal& numPtRows,
00456                                        LocalOrdinal& numPtCols,
00457                                        Teuchos::ArrayRCP<Scalar>& blockEntry);
00458 
00460 
00470   void getLocalBlockEntryView(LocalOrdinal localBlockRow,
00471                               LocalOrdinal localBlockCol,
00472                               LocalOrdinal& numPtRows,
00473                               LocalOrdinal& numPtCols,
00474                               Teuchos::ArrayRCP<const Scalar>& blockEntry) const;
00475 
00477 
00493   void getLocalBlockEntryViewNonConst(LocalOrdinal localBlockRow,
00494                                       LocalOrdinal localBlockCol,
00495                                       LocalOrdinal& numPtRows,
00496                                       LocalOrdinal& numPtCols,
00497                                       Teuchos::ArrayRCP<Scalar>& blockEntry);
00498 
00500 
00504   void getLocalDiagCopy (Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node>& diag) const;
00505 
00506   Teuchos::RCP<const BlockCrsGraph<LocalOrdinal, GlobalOrdinal, Node> >
00507   getBlockCrsGraph () {
00508     return constBlkGraph_;
00509   }
00510 
00512 
00513 
00514 
00515   virtual bool checkSizes (const SrcDistObject& source);
00516 
00517   virtual void 
00518   copyAndPermute (const SrcDistObject& source, 
00519       size_t numSameIDs, 
00520       const Teuchos::ArrayView<const LocalOrdinal>& permuteToLIDs, 
00521       const Teuchos::ArrayView<const LocalOrdinal>& permuteFromLIDs);
00522 
00523   virtual void 
00524   packAndPrepare (const SrcDistObject& source, 
00525       const Teuchos::ArrayView<const LocalOrdinal>& exportLIDs, 
00526       Teuchos::Array<char>& exports, 
00527       const Teuchos::ArrayView<size_t>& numPacketsPerLID, 
00528       size_t& constantNumPackets, 
00529       Distributor& distor);
00530 
00531   virtual void
00532   unpackAndCombine (const Teuchos::ArrayView<const LocalOrdinal>& importLIDs, 
00533         const Teuchos::ArrayView<const char>& imports, 
00534         const Teuchos::ArrayView<size_t>& numPacketsPerLID, 
00535         size_t constantNumPackets, 
00536         Distributor& distor, 
00537         CombineMode CM);
00538 
00540 
00541 
00542 
00543   std::string description() const;
00544 
00546   void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel=Teuchos::Describable::verbLevel_default) const;
00548 
00549  private:
00550   Teuchos::RCP<Node> getNode() const;
00551 
00552   void updateImport(const MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>& X) const;
00553   void updateExport(const MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>& Y) const;
00554 
00555   void createImporterExporter();
00556   void optimizeStorage();
00557   void fillLocalMatrix();
00558   void fillLocalMatVec();
00559 
00560   //private data members:
00561 
00562   //We hold two graph pointers, one const and the other non-const.
00563   //If a BlockCrsGraph is provided at construction, it is const and VbrMatrix
00564   //never changes it.
00565   //If a BlockCrsGraph is not provided at construction, VbrMatrix creates one
00566   //internally and fills it as the matrix is filled, up until fillComplete()
00567   //is called.
00568   //
00569   //blkGraph_ is either the internally created graph, or is null.
00570   //constBlkGraph_ is either a pointer to blkGraph_, or a pointer to the
00571   //graph provided at construction time.
00572   //
00573   Teuchos::RCP<BlockCrsGraph<LocalOrdinal,GlobalOrdinal,Node> > blkGraph_;
00574   Teuchos::RCP<const BlockCrsGraph<LocalOrdinal,GlobalOrdinal,Node> > constBlkGraph_;
00575 
00576   typedef KokkosClassic::VbrMatrix<Scalar, LocalOrdinal, Node> local_matrix_type;
00577 
00578   // We use a pointer so that the _decl header file only needs a
00579   // forward declaration, not an include.  We don't really need an RCP
00580   // here, since the local matrix object never gets shared outside the
00581   // class.  std::unique_ptr (C++11) would be more appropriate.
00582   Teuchos::RCP<local_matrix_type> lclMatrix_;
00583 
00584   //A variable-block-row matrix is represented by 6 arrays
00585   //in packed (contiguous storage) form. For a description of these
00586   //arrays, see the text at the bottom of this file.
00587   //(2 of those arrays, rptr and cptr, are represented by arrays in the
00588   //getBlockRowMap() and getBlockColMap() objects, and
00589   //another two of those arrays, bptr and bindx, are represented by arrays in
00590   //the BlockCrsGraph object.)
00591   //This is noted in the comments for rptr,cptr,bptr,bindx below.
00592   //
00593   //These arrays are handled as if they may point to memory that resides on
00594   //a separate device (e.g., a GPU). In other words, when the contents of these
00595   //arrays are manipulated, we use views or buffers obtained from the Node
00596   //object.
00597   Teuchos::ArrayRCP<Scalar> pbuf_values1D_;
00598   Teuchos::ArrayRCP<LocalOrdinal> pbuf_indx_;
00599 
00600   LocalMatOps lclMatOps_;
00601   Teuchos::RCP<Tpetra::Import<LocalOrdinal,GlobalOrdinal,Node> > importer_;
00602   Teuchos::RCP<Tpetra::Export<LocalOrdinal,GlobalOrdinal,Node> > exporter_;
00603   mutable Teuchos::RCP<Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > importedVec_;
00604   mutable Teuchos::RCP<Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > exportedVec_;
00605 
00606   typedef typename std::map<GlobalOrdinal,Teuchos::ArrayRCP<Scalar> > RowGlobalCols;
00607 
00608   //We use an array-of-maps to represent the variable-block-row matrix in
00609   //un-packed '2D' form.
00610   //
00611   //This unpacked data is assumed to be resident in CPU (host) memory.
00612   //It doesn't make sense to copy this data back and forth to a separate
00613   //compute device (e.g., a GPU), since we don't support doing matrix-vector
00614   //products until after fillComplete is called, at which time contiguous
00615   //arrays are allocated on the device and matrix data is copied into them.
00616   Teuchos::RCP<Teuchos::Array<RowGlobalCols> > data_2D_;
00617 
00618   VbrUtils::VbrData<LocalOrdinal,GlobalOrdinal,Scalar> nonlocal_data_;
00619 
00620   bool is_fill_completed_;
00621   bool is_storage_optimized_;
00622 };//class VbrMatrix
00623 
00624 }//namespace Tpetra
00625 
00626 //----------------------------------------------------------------------------
00627 // Description of arrays representing the VBR format:
00628 //
00629 // (For more description, see this URL (valid as of 5/26/2010):
00630 // http://docs.sun.com/source/817-0086-10/prog-sparse-support.html)
00631 // ...and of course more can be found using google...
00632 // The old Aztec manual was a great resource for this but I can't
00633 // find a copy of that these days...
00634 //
00635 //
00636 // Here is a brief description of the 6 arrays that are required to
00637 // represent a VBR matrix in packed (contiguous-memory-storage) format:
00638 //
00639 // rptr: length num_block_rows + 1
00640 //       rptr[i]: the pt-row corresponding to the i-th block-row
00641 //       Note: rptr is getBlockRowMap()->getNodeFirstPointInBlocks().
00642 //
00643 // cptr: length num_distinct_block_cols + 1
00644 //       cptr[j]: the pt-col corresponding to the j-th block-col
00645 //       Note: cptr is getBlockColMap()->getNodeFirstPointInBlocks().
00646 //
00647 // bptr: length num_block_rows + 1
00648 //       bptr[i]: location in bindx of the first nonzero block-entry
00649 //                of the i-th block-row
00650 //       Note: bptr is blkGraph_->getNodeRowOffsets();
00651 //
00652 // bindx: length num-nonzero-block-entries
00653 //        bindx[j]: block-col-index of j-th block-entry
00654 //        Note: bindx is blkGraph_->getNodePackedIndices();
00655 //
00656 // indx: length num-nonzero-block-entries + 1
00657 //       indx[j]: location in vals of the beginning of the j-th
00658 //       block-entry
00659 //
00660 // vals: length num-nonzero-scalar-entries
00661 //
00662 //
00663 // Some example look-ups:
00664 //
00665 // int nbr = num_block_rows;
00666 // int total_num_block_nonzeros = bptr[nbr];
00667 // int total_num_scalar_nonzeros = indx[num_block_nonzeros];
00668 //
00669 // //get arrays for i-th block-row:
00670 // int* bindx_i = &bindx[bptr[i]];
00671 // double* vals_i = &val[indx[bptr[i]]];
00672 // int num_block_nonzeros_in_row_i = bptr[i+1]-bptr[i];
00673 //
00674 //----------------------------------------------------------------------------
00675 
00676 #endif //TPETRA_VBRMATRIX_DECL_HPP
00677 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines