Tpetra Matrix/Vector Services: Tpetra_Experimental_BlockCrsMatrix

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_EXPERIMENTAL_BLOCKCRSMATRIX_DECL_HPP
00043 #define TPETRA_EXPERIMENTAL_BLOCKCRSMATRIX_DECL_HPP
00044 
00047 
00048 #include <Tpetra_CrsGraph.hpp>
00049 #include <Tpetra_RowMatrix.hpp>
00050 #include <Tpetra_Experimental_BlockMultiVector.hpp>
00051 #include "Tpetra_ConfigDefs.hpp"
00052 
00053 namespace Tpetra {
00054 namespace Experimental {
00055 
00122 template<class Scalar = BlockMultiVector<>::scalar_type,
00123          class LO = typename BlockMultiVector<Scalar>::local_ordinal_type,
00124          class GO = typename BlockMultiVector<Scalar, LO>::global_ordinal_type,
00125          class Node = typename BlockMultiVector<Scalar, LO, GO>::node_type>
00126 class BlockCrsMatrix :
00127   virtual public Tpetra::RowMatrix<Scalar, LO, GO, Node>,
00128   virtual public Tpetra::DistObject<char, LO, GO, Node>
00129 {
00130 private:
00131   typedef Tpetra::DistObject<char, LO, GO, Node> dist_object_type;
00132   typedef BlockMultiVector<Scalar, LO, GO, Node> BMV;
00133   typedef Teuchos::ScalarTraits<Scalar> STS;
00134 
00135 protected:
00136   typedef char packet_type;
00137 
00138 public:
00140 
00141 
00143   typedef Scalar scalar_type;
00145   typedef LO local_ordinal_type;
00147   typedef GO global_ordinal_type;
00149   typedef Node node_type;
00150 
00151   typedef ::Tpetra::Map<LO, GO, node_type> map_type;
00152   typedef Tpetra::MultiVector<scalar_type, LO, GO, node_type> mv_type;
00153   typedef Tpetra::CrsGraph<LO, GO, node_type> crs_graph_type;
00154 
00155   typedef LittleBlock<Scalar, LO> little_block_type;
00156   typedef LittleBlock<const Scalar, LO> const_little_block_type;
00157   typedef LittleVector<Scalar, LO> little_vec_type;
00158   typedef LittleVector<const Scalar, LO> const_little_vec_type;
00159 
00161 
00162 
00163 
00165   BlockCrsMatrix ();
00166 
00176   BlockCrsMatrix (const crs_graph_type& graph, const LO blockSize);
00177 
00185   BlockCrsMatrix (const crs_graph_type& graph,
00186                   const map_type& domainPointMap,
00187                   const map_type& rangePointMap,
00188                   const LO blockSize);
00189 
00191   virtual ~BlockCrsMatrix () {}
00192 
00194 
00195 
00196 
00198   Teuchos::RCP<const map_type> getDomainMap () const;
00199 
00201   Teuchos::RCP<const map_type> getRangeMap () const;
00202 
00204   Teuchos::RCP<const map_type> getRowMap () const;
00205 
00207   Teuchos::RCP<const map_type> getColMap () const;
00208 
00210   global_size_t getGlobalNumRows() const;
00211 
00213   size_t getNodeNumRows() const;
00214 
00215   size_t getNodeMaxNumRowEntries() const;
00216 
00226   void
00227   apply (const mv_type& X,
00228          mv_type& Y,
00229          Teuchos::ETransp mode = Teuchos::NO_TRANS,
00230          scalar_type alpha = Teuchos::ScalarTraits<scalar_type>::one (),
00231          scalar_type beta = Teuchos::ScalarTraits<scalar_type>::zero ()) const;
00232 
00235   bool hasTransposeApply () const {
00236     // FIXME (mfh 04 May 2014) Transpose and conjugate transpose modes
00237     // are not implemented yet.  Fill in applyBlockTrans() to fix this.
00238     return false;
00239   }
00240 
00242   void setAllToScalar (const Scalar &alpha);
00243 
00245 
00246 
00247 
00249   std::string description () const;
00250 
00274   void
00275   describe (Teuchos::FancyOStream& out,
00276             const Teuchos::EVerbosityLevel verbLevel) const;
00277 
00279 
00280 
00281 
00283   LO getBlockSize () const { return blockSize_; }
00284 
00286   virtual Teuchos::RCP<const Tpetra::RowGraph<LO,GO,Node> > getGraph () const;
00287 
00288   crs_graph_type getCrsGraph () const { return graph_; }
00289 
00294   void
00295   applyBlock (const BlockMultiVector<Scalar, LO, GO, Node>& X,
00296               BlockMultiVector<Scalar, LO, GO, Node>& Y,
00297               Teuchos::ETransp mode = Teuchos::NO_TRANS,
00298               const Scalar alpha = Teuchos::ScalarTraits<Scalar>::one (),
00299               const Scalar beta = Teuchos::ScalarTraits<Scalar>::zero ());
00300 
00304   void
00305   gaussSeidelCopy (MultiVector<Scalar,LO,GO,Node> &X,
00306                    const MultiVector<Scalar,LO,GO,Node> &B,
00307                    const MultiVector<Scalar,LO,GO,Node> &D,
00308                    const Scalar& dampingFactor,
00309                    const ESweepDirection direction,
00310                    const int numSweeps,
00311                    const bool zeroInitialGuess) const;
00312 
00316   void
00317   reorderedGaussSeidelCopy (MultiVector<Scalar,LO,GO,Node>& X,
00318                             const MultiVector<Scalar,LO,GO,Node>& B,
00319                             const MultiVector<Scalar,LO,GO,Node>& D,
00320                             const ArrayView<LO>& rowIndices,
00321                             const Scalar& dampingFactor,
00322                             const ESweepDirection direction,
00323                             const int numSweeps,
00324                             const bool zeroInitialGuess) const;
00325 
00326 
00334   void
00335   localGaussSeidel (const BlockMultiVector<Scalar, LO, GO, Node>& Residual,
00336                           BlockMultiVector<Scalar, LO, GO, Node>& Solution,
00337                           BlockCrsMatrix<Scalar, LO, GO, Node> & factorizedDiagonal,
00338                           const int * factorizationPivots,
00339                           const Scalar omega,
00340                           const ESweepDirection direction) const;
00341 
00360   LO
00361   replaceLocalValues (const LO localRowInd,
00362                       const LO colInds[],
00363                       const Scalar vals[],
00364                       const LO numColInds) const;
00365 
00384   LO
00385   sumIntoLocalValues (const LO localRowInd,
00386                       const LO colInds[],
00387                       const Scalar vals[],
00388                       const LO numColInds) const;
00389 
00399   LO
00400   getLocalRowView (const LO localRowInd,
00401                    const LO*& colInds,
00402                    Scalar*& vals,
00403                    LO& numInds) const;
00404 
00405   void
00406   getLocalRowView (LO LocalRow,
00407                    ArrayView<const LO> &indices,
00408                    ArrayView<const Scalar> &values) const;
00409 
00410   void
00411   getLocalRowCopy (LO LocalRow,
00412                    const ArrayView<LO> &Indices,
00413                    const ArrayView<Scalar> &Values,
00414                    size_t &NumEntries) const;
00415 
00416 
00417   little_block_type getLocalBlock (const LO localRowInd, const LO localColInd) const;
00418 
00441   LO
00442   getLocalRowOffsets (const LO localRowInd,
00443                       ptrdiff_t offsets[],
00444                       const LO colInds[],
00445                       const LO numColInds) const;
00446 
00452   LO
00453   replaceLocalValuesByOffsets (const LO localRowInd,
00454                                const ptrdiff_t offsets[],
00455                                const Scalar vals[],
00456                                const LO numOffsets) const;
00457 
00463   LO
00464   sumIntoLocalValuesByOffsets (const LO localRowInd,
00465                                const ptrdiff_t offsets[],
00466                                const Scalar vals[],
00467                                const LO numOffsets) const;
00468 
00475   size_t getNumEntriesInLocalRow (const LO localRowInd) const;
00476 
00493   bool localError () const {
00494     return *localError_;
00495   }
00496 
00511   std::string errorMessages () const {
00512     return (*errs_).is_null () ? std::string ("") : (*errs_)->str ();
00513   }
00514 
00546   void getLocalDiagOffsets (Teuchos::ArrayRCP<size_t>& offsets) const;
00547 
00557   void
00558   getLocalDiagCopy (BlockCrsMatrix<Scalar,LO,GO,Node>& diag,
00559                     const Teuchos::ArrayView<const size_t>& offsets) const;
00560 
00561 
00563   void computeDiagonalGraph ();
00564 
00566   bool isComputedDiagonalGraph() const { return computedDiagonalGraph_;}
00567 
00568   Teuchos::RCP<crs_graph_type> getDiagonalGraph () const;
00569 
00570 
00571 protected:
00573   LO
00574   absMaxLocalValues (const LO localRowInd,
00575                      const LO colInds[],
00576                      const Scalar vals[],
00577                      const LO numColInds) const;
00578 
00580   LO
00581   absMaxLocalValuesByOffsets (const LO localRowInd,
00582                               const ptrdiff_t offsets[],
00583                               const Scalar vals[],
00584                               const LO numOffsets) const;
00585 
00592 
00593 
00594   virtual bool checkSizes (const Tpetra::SrcDistObject& source);
00595 
00596   virtual void
00597   copyAndPermute (const Tpetra::SrcDistObject& source,
00598                   size_t numSameIDs,
00599                   const Teuchos::ArrayView<const LO>& permuteToLIDs,
00600                   const Teuchos::ArrayView<const LO>& permuteFromLIDs);
00601 
00602   virtual void
00603   packAndPrepare (const Tpetra::SrcDistObject& source,
00604                   const Teuchos::ArrayView<const LO>& exportLIDs,
00605                   Teuchos::Array<packet_type>& exports,
00606                   const Teuchos::ArrayView<size_t>& numPacketsPerLID,
00607                   size_t& constantNumPackets,
00608                   Tpetra::Distributor& distor);
00609 
00610   virtual void
00611   unpackAndCombine (const Teuchos::ArrayView<const LO> &importLIDs,
00612                     const Teuchos::ArrayView<const packet_type> &imports,
00613                     const Teuchos::ArrayView<size_t> &numPacketsPerLID,
00614                     size_t constantNumPackets,
00615                     Tpetra::Distributor& distor,
00616                     Tpetra::CombineMode CM);
00618 
00619 private:
00621   crs_graph_type graph_;
00622   Teuchos::RCP<crs_graph_type> graphRCP_;
00631   map_type rowMeshMap_;
00638   map_type domainPointMap_;
00645   map_type rangePointMap_;
00647   LO blockSize_;
00649   const size_t* ptr_;
00651   const LO* ind_;
00657   Teuchos::ArrayRCP<Scalar> valView_;
00662   Scalar* val_;
00663 
00685   Teuchos::RCP<Teuchos::RCP<BMV> > X_colMap_;
00689   Teuchos::RCP<Teuchos::RCP<BMV> > Y_rowMap_;
00690 
00698   LO columnPadding_;
00699 
00701   bool rowMajor_;
00702 
00714   Teuchos::RCP<bool> localError_;
00715 
00723   Teuchos::RCP<Teuchos::RCP<std::ostringstream> > errs_;
00724 
00734   void
00735   applyBlockTrans (const BlockMultiVector<Scalar, LO, GO, Node>& X,
00736                    BlockMultiVector<Scalar, LO, GO, Node>& Y,
00737                    const Teuchos::ETransp mode,
00738                    const Scalar alpha,
00739                    const Scalar beta);
00740 
00748   void
00749   applyBlockNoTrans (const BlockMultiVector<Scalar, LO, GO, Node>& X,
00750                      BlockMultiVector<Scalar, LO, GO, Node>& Y,
00751                      const Scalar alpha,
00752                      const Scalar beta);
00753 
00761   void
00762   localApplyBlockNoTrans (const BlockMultiVector<Scalar, LO, GO, Node>& X,
00763                           BlockMultiVector<Scalar, LO, GO, Node>& Y,
00764                           const Scalar alpha,
00765                           const Scalar beta);
00766 
00767   Teuchos::RCP<crs_graph_type> diagonalGraph_;
00768   bool computedDiagonalGraph_;
00769 
00770 
00803   size_t
00804   findRelOffsetOfColumnIndex (const LO localRowIndex,
00805                               const LO colIndexToFind,
00806                               const size_t hint) const;
00807 
00810   LO offsetPerBlock () const;
00811 
00812   const_little_block_type
00813   getConstLocalBlockFromInput (const Scalar* val, const size_t pointOffset) const;
00814 
00815   little_block_type
00816   getNonConstLocalBlockFromInput (Scalar* val, const size_t pointOffset) const;
00817 
00818   const_little_block_type
00819   getConstLocalBlockFromAbsOffset (const size_t absBlockOffset) const;
00820 
00821   little_block_type
00822   getNonConstLocalBlockFromAbsOffset (const size_t absBlockOffset) const;
00823 
00824 public:
00826   virtual Teuchos::RCP<const Teuchos::Comm<int> > getComm() const;
00827 
00829   virtual Teuchos::RCP<Node> getNode() const;
00830 
00832   virtual global_size_t getGlobalNumCols() const;
00833 
00834   virtual size_t getNodeNumCols() const;
00835 
00836   virtual GO getIndexBase() const;
00837 
00839   virtual global_size_t getGlobalNumEntries() const;
00840 
00842   virtual size_t getNodeNumEntries() const;
00843 
00853   virtual size_t getNumEntriesInGlobalRow (GO globalRow) const;
00854 
00856   virtual global_size_t getGlobalNumDiags() const;
00857 
00859   virtual size_t getNodeNumDiags() const;
00860 
00862   virtual size_t getGlobalMaxNumRowEntries() const;
00863 
00865   virtual bool hasColMap() const;
00866 
00868   virtual bool isLowerTriangular() const;
00869 
00871   virtual bool isUpperTriangular() const;
00872 
00882   virtual bool isLocallyIndexed() const;
00883 
00893   virtual bool isGloballyIndexed() const;
00894 
00896   virtual bool isFillComplete() const;
00897 
00899   virtual bool supportsRowViews() const;
00900 
00902 
00903 
00904 
00925   virtual void
00926   getGlobalRowCopy (GO GlobalRow,
00927                     const Teuchos::ArrayView<GO> &Indices,
00928                     const Teuchos::ArrayView<Scalar> &Values,
00929                     size_t &NumEntries) const;
00930 
00955   virtual void
00956   getGlobalRowView (GO GlobalRow,
00957                     ArrayView<const GO> &indices,
00958                     ArrayView<const Scalar> &values) const;
00959 
00971   virtual void getLocalDiagCopy (Vector<Scalar,LO,GO,Node> &diag) const;
00972 
00974 
00975 
00976 
00982   virtual void leftScale (const Vector<Scalar, LO, GO, Node>& x);
00983 
00989   virtual void rightScale (const Vector<Scalar, LO, GO, Node>& x);
00990 
00999   virtual typename ScalarTraits<Scalar>::magnitudeType getFrobeniusNorm() const;
01000 
01001 };
01002 
01003 } // namespace Experimental
01004 } // namespace Tpetra
01005 
01006 #endif // TPETRA_EXPERIMENTAL_BLOCKCRSMATRIX_DECL_HPP