Branch data     Line data    Source code

       1                 :            : // *****************************************************************************
       2                 :            : /*!
       3                 :            :   \file      src/LinearSolver/CSR.hpp
       4                 :            :   \copyright 2012-2015 J. Bakosi,
       5                 :            :              2016-2018 Los Alamos National Security, LLC.,
       6                 :            :              2019-2021 Triad National Security, LLC.
       7                 :            :              All rights reserved. See the LICENSE file for details.
       8                 :            :   \brief     Compressed sparse row (CSR) storage for a sparse matrix
       9                 :            :   \details   Compressed sparse row (CSR) storage for a sparse matrix.
      10                 :            : */
      11                 :            : // *****************************************************************************
      12                 :            : #ifndef CSR_h
      13                 :            : #define CSR_h
      14                 :            : 
      15                 :            : #include <vector>
      16                 :            : #include <ostream>
      17                 :            : 
      18                 :            : #include "NoWarning/pup_stl.hpp"
      19                 :            : 
      20                 :            : #include "Types.hpp"
      21                 :            : #include "CommMap.hpp"
      22                 :            : 
      23                 :            : namespace tk {
      24                 :            : 
      25                 :            : //! Compressed sparse row (CSR) storage for a sparse matrix
      26                 :            : class CSR {
      27                 :            : 
      28                 :            :   public:
      29                 :            :     //! \brief Constructor: Create a CSR symmetric matrix with ncomp scalar
      30                 :            :     //!   components, storing only the upper triangular part
      31                 :            :     explicit CSR( std::size_t nc,
      32                 :            :                   const std::pair< std::vector< std::size_t >,
      33                 :            :                                    std::vector< std::size_t > >& psup );
      34                 :            : 
      35                 :            :     //! Empty constructor for Charm++
      36                 :         80 :     explicit CSR() = default;
      37                 :            : 
      38                 :            :     //! Return const reference to sparse matrix entry at a position
      39                 :            :     const real&
      40                 :            :     operator()( std::size_t row, std::size_t col, std::size_t pos=0 ) const;
      41                 :            : 
      42                 :            :     //! Return non-const reference to sparse matrix entry at a position
      43                 :            :     //! \see "Avoid Duplication in const and Non-const Member Function," and
      44                 :            :     //!   "Use const whenever possible," Scott Meyers, Effective C++, 3d ed.
      45                 :            :     real&
      46                 :   13371648 :     operator()( std::size_t row, std::size_t col, std::size_t pos=0 ) {
      47                 :            :       return const_cast< real& >(
      48                 :   13371648 :                static_cast< const CSR& >( *this ).operator()( row, col, pos ) );
      49                 :            :     }
      50                 :            : 
      51                 :            :     //! Set Dirichlet boundary condition at a node
      52                 :            :     void dirichlet(
      53                 :            :       std::size_t i,
      54                 :            :       const std::vector< std::size_t >& gid = {},
      55                 :            :       const NodeCommMap& nodecommap = {},
      56                 :            :       std::size_t pos=0 );
      57                 :            : 
      58                 :            :     //! Multiply CSR matrix with vector from the right: r = A * x
      59                 :            :     void mult( const std::vector< real >& x, std::vector< real >& r,
      60                 :            :                const std::vector< tk::real >& bcmask ) const;
      61                 :            : 
      62                 :            :     //! Access real size of matrix
      63                 :        140 :     std::size_t rsize() const { return rnz.size()*ncomp; }
      64                 :            : 
      65                 :            :     //! Access the number of scalar components per non-zero matrix entry
      66                 :      24282 :     std::size_t Ncomp() const { return ncomp; }
      67                 :            : 
      68                 :            :     //! Write out CSR as stored
      69                 :            :     std::ostream& write_stored( std::ostream &os ) const;
      70                 :            :     //! Write out CSR nonzero structure
      71                 :            :     std::ostream& write_structure( std::ostream &os ) const;
      72                 :            :     //! Write out CSR contribute structure
      73                 :            :     std::ostream& write_contribute( std::ostream &os ) const;
      74                 :            :     //! Write out CSR as a real matrix
      75                 :            :     std::ostream& write_matrix( std::ostream &os ) const;
      76                 :            :     //! Write out CSR in Matlab/Octave format
      77                 :            :     std::ostream& write_matlab( std::ostream &os ) const;
      78                 :            : 
      79                 :            :     /** @name Pack/unpack (Charm++ serialization) routines */
      80                 :            :     ///@{
      81                 :            :     //! \brief Pack/Unpack serialize member function
      82                 :            :     //! \param[in,out] p Charm++'s PUP::er serializer object reference
      83                 :        240 :     void pup( PUP::er &p ) {
      84                 :        240 :       p | ncomp;
      85                 :        240 :       p | rnz;
      86                 :        240 :       p | ia;
      87                 :        240 :       p | ja;
      88                 :        240 :       p | a;
      89                 :        240 :     }
      90                 :            :     //! \brief Pack/Unpack serialize operator|
      91                 :            :     //! \param[in,out] p Charm++'s PUP::er serializer object reference
      92                 :            :     //! \param[in,out] c CSR object reference
      93                 :        240 :     friend void operator|( PUP::er& p, CSR& c ) { c.pup(p); }
      94                 :            :     ///@}
      95                 :            : 
      96                 :            :   private:
      97                 :            :     std::size_t ncomp;                  //!< Number of scalars per non-zero
      98                 :            :     std::vector< std::size_t > rnz;     //!< Number of nonzeros in each row
      99                 :            :     std::vector< std::size_t > ia;      //!< Row pointers
     100                 :            :     std::vector< std::size_t > ja;      //!< Column indices
     101                 :            :     std::vector< real > a;              //!< Nonzero matrix values
     102                 :            : };
     103                 :            : 
     104                 :            : } // tk::
     105                 :            : 
     106                 :            : #endif // CSR_h