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 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_DECL_HPP
00030 #define TPETRA_VBRMATRIX_DECL_HPP
00031 
00032 #include <Kokkos_DefaultNode.hpp>
00033 #include <Kokkos_DefaultKernels.hpp>
00034 #include <Kokkos_VbrMatrix.hpp>
00035 
00036 #include "Tpetra_ConfigDefs.hpp"
00037 #include "Tpetra_Operator.hpp"
00038 #include "Tpetra_BlockMap.hpp"
00039 #include "Tpetra_BlockCrsGraph.hpp"
00040 #include "Tpetra_VbrUtils.hpp"
00041 
00046 namespace Tpetra {
00047 
00049 
00077 template <class Scalar, 
00078           class LocalOrdinal  = int, 
00079           class GlobalOrdinal = LocalOrdinal, 
00080           class Node          = Kokkos::DefaultNode::DefaultNodeType, 
00081           class LocalMatOps   = typename Kokkos::DefaultKernels<Scalar,LocalOrdinal,Node>::BlockSparseOps >
00082 class VbrMatrix : public Tpetra::DistObject<char, LocalOrdinal, GlobalOrdinal, Node>, public Tpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node> {
00083  public:
00084   typedef Scalar        scalar_type;
00085   typedef LocalOrdinal  local_ordinal_type;
00086   typedef GlobalOrdinal global_ordinal_type;
00087   typedef Node          node_type;
00088   typedef LocalMatOps   mat_vec_type;
00089 
00091 
00092 
00094 
00099   VbrMatrix(const Teuchos::RCP<const BlockMap<LocalOrdinal,GlobalOrdinal,Node> > &blkRowMap, size_t maxNumEntriesPerRow, ProfileType pftype = DynamicProfile);
00100 
00102 
00112   VbrMatrix(const Teuchos::RCP<const BlockCrsGraph<LocalOrdinal,GlobalOrdinal,Node> >& blkGraph);
00113 
00115   virtual ~VbrMatrix();
00116 
00118 
00120 
00121 
00123 
00128   template <class DomainScalar, class RangeScalar>
00129       void multiply(const MultiVector<DomainScalar,LocalOrdinal,GlobalOrdinal,Node> & X, MultiVector<RangeScalar,LocalOrdinal,GlobalOrdinal,Node> &Y, Teuchos::ETransp trans, RangeScalar alpha, RangeScalar beta) const;
00130 
00132 
00144   template <class DomainScalar, class RangeScalar>
00145       void solve(const MultiVector<RangeScalar,LocalOrdinal,GlobalOrdinal,Node> & Y, MultiVector<DomainScalar,LocalOrdinal,GlobalOrdinal,Node> &X, Teuchos::ETransp trans) const;
00146 
00148 
00150 
00151 
00153 
00155   const Teuchos::RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > & getDomainMap() const;
00156 
00158 
00160   const Teuchos::RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > & getRangeMap() const;
00161 
00163 
00168   void apply(const MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> &X,
00169              MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> &Y,
00170              Teuchos::ETransp trans = Teuchos::NO_TRANS,
00171              Scalar alpha = Teuchos::ScalarTraits<Scalar>::one(),
00172              Scalar beta = Teuchos::ScalarTraits<Scalar>::zero()) const;
00173 
00175 
00178   void applyInverse(const MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> & Y,
00179                     MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> &X,
00180                     Teuchos::ETransp trans) const;
00181 
00183 
00186   bool hasTransposeApply() const;
00187 
00189 
00191 
00192 
00194   const Teuchos::RCP<const BlockMap<LocalOrdinal,GlobalOrdinal,Node> > & getBlockRowMap() const;
00195 
00197   const Teuchos::RCP<const BlockMap<LocalOrdinal,GlobalOrdinal,Node> > & getBlockColMap() const;
00198 
00200   const Teuchos::RCP<const BlockMap<LocalOrdinal,GlobalOrdinal,Node> > & getBlockDomainMap() const;
00201 
00203   const Teuchos::RCP<const BlockMap<LocalOrdinal,GlobalOrdinal,Node> > & getBlockRangeMap() const;
00204 
00206   const Teuchos::RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > & getPointRowMap() const;
00207 
00209   const Teuchos::RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > & getPointColMap() const;
00210 
00212   bool isFillComplete() const;
00214 
00216 
00217 
00219 
00222   void putScalar(Scalar s);
00223 
00225 
00240   void setGlobalBlockEntry(GlobalOrdinal globalBlockRow, GlobalOrdinal globalBlockCol, const Teuchos::SerialDenseMatrix<GlobalOrdinal,Scalar>& blockEntry);
00241 
00243 
00250   void setLocalBlockEntry(LocalOrdinal localBlockRow, LocalOrdinal localBlockCol, const Teuchos::SerialDenseMatrix<LocalOrdinal,Scalar>& blockEntry);
00251 
00253 
00268   void sumIntoGlobalBlockEntry(GlobalOrdinal globalBlockRow, GlobalOrdinal globalBlockCol, const Teuchos::SerialDenseMatrix<GlobalOrdinal,Scalar>& blockEntry);
00269 
00271 
00278   void sumIntoLocalBlockEntry(LocalOrdinal localBlockRow, LocalOrdinal localBlockCol, const Teuchos::SerialDenseMatrix<LocalOrdinal,Scalar>& blockEntry);
00279 
00281 
00296   void setGlobalBlockEntry(GlobalOrdinal globalBlockRow, GlobalOrdinal globalBlockCol, LocalOrdinal blkRowSize, LocalOrdinal blkColSize, LocalOrdinal LDA, const Teuchos::ArrayView<const Scalar>& blockEntry);
00297 
00299 
00306   void setLocalBlockEntry(LocalOrdinal localBlockRow, LocalOrdinal localBlockCol, LocalOrdinal blkRowSize, LocalOrdinal blkColSize, LocalOrdinal LDA, const Teuchos::ArrayView<const Scalar>& blockEntry);
00307 
00309 
00324   void sumIntoGlobalBlockEntry(GlobalOrdinal globalBlockRow, GlobalOrdinal globalBlockCol, LocalOrdinal blkRowSize, LocalOrdinal blkColSize, LocalOrdinal LDA, const Teuchos::ArrayView<const Scalar>& blockEntry);
00325 
00327 
00334   void sumIntoLocalBlockEntry(LocalOrdinal localBlockRow, LocalOrdinal localBlockCol, LocalOrdinal blkRowSize, LocalOrdinal blkColSize, LocalOrdinal LDA, const Teuchos::ArrayView<const Scalar>& blockEntry);
00335 
00337 
00339 
00340 
00342 
00346   void fillComplete(const Teuchos::RCP<const BlockMap<LocalOrdinal,GlobalOrdinal,Node> >& blockDomainMap, const Teuchos::RCP<const BlockMap<LocalOrdinal,GlobalOrdinal,Node> >& blockRangeMap);
00347 
00349 
00352   void fillComplete();
00353 
00355   void globalAssemble();
00357 
00359 
00360 
00362 
00364   void getGlobalBlockRowView(GlobalOrdinal globalBlockRow,
00365                              LocalOrdinal& numPtRows,
00366                              Teuchos::ArrayView<const GlobalOrdinal>& blockCols,
00367                              Teuchos::Array<LocalOrdinal>& ptColsPerBlockCol,
00368                              Teuchos::Array<Teuchos::ArrayRCP<const Scalar> >& blockEntries) const;
00369 
00371 
00373   void getLocalBlockRowView(LocalOrdinal localBlockRow,
00374                             LocalOrdinal& numPtRows,
00375                             Teuchos::ArrayView<const LocalOrdinal>& blockCols,
00376                             Teuchos::Array<LocalOrdinal>& ptColsPerBlockCol,
00377                             Teuchos::ArrayRCP<const Scalar>& blockEntries) const;
00378 
00380 
00388   void getGlobalBlockEntryView(GlobalOrdinal globalBlockRow,
00389                                GlobalOrdinal globalBlockCol,
00390                                LocalOrdinal& numPtRows,
00391                                LocalOrdinal& numPtCols,
00392                                Teuchos::ArrayRCP<const Scalar>& blockEntry) const;
00393 
00395 
00405   void getGlobalBlockEntryViewNonConst(GlobalOrdinal globalBlockRow,
00406                                        GlobalOrdinal globalBlockCol,
00407                                        LocalOrdinal& numPtRows,
00408                                        LocalOrdinal& numPtCols,
00409                                        Teuchos::ArrayRCP<Scalar>& blockEntry);
00410 
00412 
00422   void getLocalBlockEntryView(LocalOrdinal localBlockRow,
00423                               LocalOrdinal localBlockCol,
00424                               LocalOrdinal& numPtRows, 
00425                               LocalOrdinal& numPtCols,
00426                               Teuchos::ArrayRCP<const Scalar>& blockEntry) const;
00427 
00429 
00445   void getLocalBlockEntryViewNonConst(LocalOrdinal localBlockRow,
00446                                       LocalOrdinal localBlockCol,
00447                                       LocalOrdinal& numPtRows,
00448                                       LocalOrdinal& numPtCols,
00449                                       Teuchos::ArrayRCP<Scalar>& blockEntry);
00450 
00452 
00456   void getLocalDiagCopy(Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node>& diag) const;
00457 
00458   const Teuchos::RCP<const BlockCrsGraph<LocalOrdinal, GlobalOrdinal, Node> >& getBlockCrsGraph() {return constBlkGraph_;}
00460 
00462 
00463 
00464   bool checkSizes(const DistObject<char, LocalOrdinal, GlobalOrdinal, Node>& source);
00465 
00466   void copyAndPermute(const DistObject<char, LocalOrdinal, GlobalOrdinal, Node>& source, size_t numSameIDs, const Teuchos::ArrayView<const LocalOrdinal>& permuteToLIDs, const Teuchos::ArrayView<const LocalOrdinal>& permuteFromLIDs);
00467 
00468   void packAndPrepare(const DistObject<char, LocalOrdinal, GlobalOrdinal, Node>& source, const Teuchos::ArrayView<const LocalOrdinal>& exportLIDs, Teuchos::Array<char>& exports, const Teuchos::ArrayView<size_t>& numPacketsPerLID, size_t& constantNumPackets, Distributor& distor);
00469 
00470   void unpackAndCombine(const Teuchos::ArrayView<const LocalOrdinal>& importLIDs, const Teuchos::ArrayView<const char>& imports, const Teuchos::ArrayView<size_t>& numPacketsPerLID, size_t constantNumPackets, Distributor& distor, CombineMode CM);
00471 
00473 
00475 
00476   std::string description() const;
00477 
00480   void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel=Teuchos::Describable::verbLevel_default) const;
00482 
00483  private:
00484   //private methods:
00485 
00486   Teuchos::RCP<Node> getNode() const;
00487 
00488   void updateImport(const MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>& X) const;
00489   void updateExport(const MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>& Y) const;
00490 
00491   void createImporterExporter();
00492   void optimizeStorage();
00493   void fillLocalMatrix();
00494   void fillLocalMatVec();
00495 
00496   //private data members:
00497 
00498   //We hold two graph pointers, one const and the other non-const.
00499   //If a BlockCrsGraph is provided at construction, it is const and VbrMatrix
00500   //never changes it.
00501   //If a BlockCrsGraph is not provided at construction, VbrMatrix creates one
00502   //internally and fills it as the matrix is filled, up until fillComplete()
00503   //is called.
00504   //
00505   //blkGraph_ is either the internally created graph, or is null.
00506   //constBlkGraph_ is either a pointer to blkGraph_, or a pointer to the
00507   //graph provided at construction time.
00508   //
00509   Teuchos::RCP<BlockCrsGraph<LocalOrdinal,GlobalOrdinal,Node> > blkGraph_;
00510   Teuchos::RCP<const BlockCrsGraph<LocalOrdinal,GlobalOrdinal,Node> > constBlkGraph_;
00511 
00512   Kokkos::VbrMatrix<Scalar,LocalOrdinal,Node> lclMatrix_;
00513 
00514   //A variable-block-row matrix is represented by 6 arrays 
00515   //in packed (contiguous storage) form. For a description of these
00516   //arrays, see the text at the bottom of this file.
00517   //(2 of those arrays, rptr and cptr, are represented by arrays in the
00518   //getBlockRowMap() and getBlockColMap() objects, and
00519   //another two of those arrays, bptr and bindx, are represented by arrays in
00520   //the BlockCrsGraph object.)
00521   //This is noted in the comments for rptr,cptr,bptr,bindx below.
00522   //
00523   //These arrays are handled as if they may point to memory that resides on
00524   //a separate device (e.g., a GPU). In other words, when the contents of these
00525   //arrays are manipulated, we use views or buffers obtained from the Node
00526   //object.
00527   Teuchos::ArrayRCP<Scalar> pbuf_values1D_;
00528   Teuchos::ArrayRCP<LocalOrdinal> pbuf_indx_;
00529 
00530   LocalMatOps lclMatOps_;
00531   Teuchos::RCP<Tpetra::Import<LocalOrdinal,GlobalOrdinal,Node> > importer_;
00532   Teuchos::RCP<Tpetra::Export<LocalOrdinal,GlobalOrdinal,Node> > exporter_;
00533   mutable Teuchos::RCP<Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > importedVec_;
00534   mutable Teuchos::RCP<Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > exportedVec_;
00535 
00536   typedef typename std::map<GlobalOrdinal,Teuchos::ArrayRCP<Scalar> > RowGlobalCols;
00537 
00538   //We use an array-of-maps to represent the variable-block-row matrix in
00539   //un-packed '2D' form.
00540   //
00541   //This unpacked data is assumed to be resident in CPU (host) memory.
00542   //It doesn't make sense to copy this data back and forth to a separate
00543   //compute device (e.g., a GPU), since we don't support doing matrix-vector
00544   //products until after fillComplete is called, at which time contiguous
00545   //arrays are allocated on the device and matrix data is copied into them.
00546   Teuchos::RCP<Teuchos::Array<RowGlobalCols> > data_2D_;
00547 
00548   VbrUtils::VbrData<LocalOrdinal,GlobalOrdinal,Scalar> nonlocal_data_;
00549 
00550   bool is_fill_completed_;
00551   bool is_storage_optimized_;
00552 };//class VbrMatrix
00553 
00554 }//namespace Tpetra
00555 
00556 //----------------------------------------------------------------------------
00557 // Description of arrays representing the VBR format:
00558 //
00559 // (For more description, see this URL (valid as of 5/26/2010):
00560 // http://docs.sun.com/source/817-0086-10/prog-sparse-support.html)
00561 // ...and of course more can be found using google...
00562 // The old Aztec manual was a great resource for this but I can't
00563 // find a copy of that these days...
00564 //
00565 //
00566 // Here is a brief description of the 6 arrays that are required to
00567 // represent a VBR matrix in packed (contiguous-memory-storage) format:
00568 //
00569 // rptr: length num_block_rows + 1
00570 //       rptr[i]: the pt-row corresponding to the i-th block-row
00571 //       Note: rptr is getBlockRowMap()->getNodeFirstPointInBlocks().
00572 //
00573 // cptr: length num_distinct_block_cols + 1
00574 //       cptr[j]: the pt-col corresponding to the j-th block-col
00575 //       Note: cptr is getBlockColMap()->getNodeFirstPointInBlocks().
00576 //
00577 // bptr: length num_block_rows + 1
00578 //       bptr[i]: location in bindx of the first nonzero block-entry
00579 //                of the i-th block-row
00580 //       Note: bptr is blkGraph_->getNodeRowOffsets();
00581 //
00582 // bindx: length num-nonzero-block-entries
00583 //        bindx[j]: block-col-index of j-th block-entry
00584 //        Note: bindx is blkGraph_->getNodePackedIndices();
00585 //
00586 // indx: length num-nonzero-block-entries + 1
00587 //       indx[j]: location in vals of the beginning of the j-th
00588 //       block-entry
00589 //
00590 // vals: length num-nonzero-scalar-entries
00591 //
00592 //
00593 // Some example look-ups:
00594 //
00595 // int nbr = num_block_rows;
00596 // int total_num_block_nonzeros = bptr[nbr];
00597 // int total_num_scalar_nonzeros = indx[num_block_nonzeros];
00598 // 
00599 // //get arrays for i-th block-row:
00600 // int* bindx_i = &bindx[bptr[i]];
00601 // double* vals_i = &val[indx[bptr[i]]];
00602 // int num_block_nonzeros_in_row_i = bptr[i+1]-bptr[i];
00603 // 
00604 //----------------------------------------------------------------------------
00605 
00606 #endif //TPETRA_VBRMATRIX_DECL_HPP
00607 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines