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                 :            :              2022-2024 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                 :            :     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                 :            :     operator()( std::size_t row, std::size_t col, std::size_t pos=0 ) {
      45                 :            :       return const_cast< real& >(
      46 [ +  - ][ +  - ]:    4778210 :                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 [ +  - ][ +  - ]:        721 :     std::size_t rsize() const { return rnz.size()*ncomp; }
         [ +  - ][ +  - ]
         [ +  - ][ +  - ]
         [ -  - ][ -  - ]
         [ -  - ][ -  - ]
         [ -  - ][ -  - ]
      67                 :            : 
      68                 :            :     //! Access the number of scalar components per non-zero matrix entry
      69                 :     690960 :     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                 :      12040 :     void pup( PUP::er &p ) {
      87                 :      12040 :       p | ncomp;
      88                 :      12040 :       p | rnz;
      89                 :      12040 :       p | ia;
      90                 :      12040 :       p | ja;
      91                 :      12040 :       p | a;
      92                 :      12040 :     }
      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 [ -  - ][ -  - ]:       7067 :     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::