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