Branch data     Line data    Source code

       1                 :            : // *****************************************************************************
       2                 :            : /*!
       3                 :            :   \file      src/LinearSolver/CSR.cpp
       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                 :            : 
      13                 :            : #include "Exception.hpp"
      14                 :            : #include "CSR.hpp"
      15                 :            : 
      16                 :            : using tk::CSR;
      17                 :            : 
      18                 :         37 : CSR::CSR( std::size_t nc,
      19                 :            :           const std::pair< std::vector< std::size_t >,
      20                 :         37 :                            std::vector< std::size_t > >& psup )
      21                 :            : try :
      22                 :            :   ncomp( nc ),
      23                 :         37 :   rnz( psup.second.size()-1 ),
      24 [ +  + ][ +  - ]:         79 :   ia( rnz.size()*ncomp+1 )
      25                 :            : // *****************************************************************************
      26                 :            : //  Constructor: Create a CSR symmetric matrix with ncomp scalar components per
      27                 :            : //  non-zero matrix entry, storing only the upper triangular part
      28                 :            : //! \param[in] nc Number of scalar components (degrees of freedom)
      29                 :            : //! \param[in] psup Points surrounding points of mesh graph, see tk::genPsup
      30                 :            : // *****************************************************************************
      31                 :            : {
      32 [ +  + ][ +  - ]:         36 :   Assert( ncomp > 0, "Sparse matrix ncomp must be positive" );
         [ +  - ][ +  - ]
      33 [ -  + ][ -  - ]:         35 :   Assert( rnz.size() > 0, "Sparse matrix size must be positive" );
         [ -  - ][ -  - ]
      34                 :            : 
      35                 :         35 :   const auto& psup1 = psup.first;
      36                 :         35 :   const auto& psup2 = psup.second;
      37                 :            : 
      38                 :            :   // Calculate number of nonzeros in each block row (rnz), total number of
      39                 :            :   // nonzeros (nnz), and fill in row indices (ia)
      40                 :            :   std::size_t nnz, i;
      41         [ +  + ]:      26602 :   for (ia[0]=1, nnz=i=0; i<psup2.size()-1; ++i) {
      42                 :            :     // add up and store nonzeros of row i (only upper triangular part)
      43                 :            :     std::size_t j;
      44         [ +  + ]:     306805 :     for (rnz[i]=1, j=psup2[i]+1; j<=psup2[i+1]; ++j)
      45                 :     280238 :       ++rnz[i];
      46                 :            : 
      47                 :            :     // add up total number of nonzeros
      48                 :      26567 :     nnz += rnz[i] * ncomp;
      49                 :            : 
      50                 :            :     // fill up row index
      51         [ +  + ]:     105360 :     for (std::size_t k=0; k<ncomp; ++k)
      52                 :      78793 :       ia[i*ncomp+k+1] = ia[i*ncomp+k] + rnz[i];
      53                 :            :   }
      54                 :            : 
      55                 :            :   // Allocate storage for matrix values and column indices
      56         [ +  - ]:         35 :   a.resize( nnz, 0.0 );
      57         [ +  - ]:         35 :   ja.resize( nnz );
      58                 :            : 
      59                 :            :   // fill column indices
      60         [ +  + ]:      26602 :   for (i=0; i<rnz.size(); ++i)
      61         [ +  + ]:     105360 :     for (std::size_t k=0; k<ncomp; ++k) {
      62                 :      78793 :       auto itmp = i*ncomp+k;
      63                 :      78793 :       ja[ia[itmp]-1] = itmp+1;  // put in column index of diagonal
      64         [ +  + ]:     912287 :       for (std::size_t n=1, j=psup2[i]+1; j<=psup2[i+1]; ++j) {
      65                 :            :         // put in column index of an off-diagonal
      66                 :     833494 :         ja[ia[itmp]-1+(n++)] = psup1[j]*ncomp+k+1;
      67                 :            :       }
      68                 :            :     }
      69                 :            : 
      70                 :            :   // (bubble-)sort column indices
      71         [ +  + ]:      26602 :   for (i=0; i<rnz.size(); ++i)
      72         [ +  + ]:     105360 :     for (std::size_t k=0; k<ncomp; ++k)
      73         [ +  + ]:     912287 :       for (std::size_t j=psup2[i]+1; j<=psup2[i+1]; ++j)
      74         [ +  + ]:   10481610 :          for (std::size_t l=1; l<rnz[i]; ++l)   // sort column indices of row i
      75         [ +  + ]:   74409488 :             for (std::size_t e=0; e<rnz[i]-l; ++e)
      76         [ +  + ]:   64761372 :               if (ja[ia[i*ncomp+k]-1+e] > ja[ia[i*ncomp+k]+e])
      77                 :     416747 :                 std::swap( ja[ia[i*ncomp+k]-1+e], ja[ia[i*ncomp+k]+e] );
      78                 :            : 
      79                 :            : } // Catch std::exception
      80                 :          4 :   catch (std::exception& se) {
      81                 :            :     // (re-)throw tk::Excpetion
      82 [ +  - ][ +  - ]:          2 :     Throw( std::string("RUNTIME ERROR in CSR constructor: ") + se.what() );
         [ +  - ][ +  - ]
      83                 :         35 :   }
      84                 :            : 
      85                 :            : const tk::real&
      86                 :   13371648 : CSR::operator()( std::size_t row, std::size_t col, std::size_t pos ) const
      87                 :            : // *****************************************************************************
      88                 :            : //  Return const reference to sparse matrix entry at a position specified
      89                 :            : //  using relative addressing
      90                 :            : //! \param[in] row Block row
      91                 :            : //! \param[in] col Block column
      92                 :            : //! \param[in] pos Position in block
      93                 :            : //! \return Const reference to matrix entry at position specified
      94                 :            : // *****************************************************************************
      95                 :            : {
      96                 :   13371648 :   auto rncomp = row * ncomp;
      97                 :            : 
      98         [ +  - ]:   94245243 :   for (std::size_t n=0, j=ia[rncomp+pos]-1; j<ia[rncomp+pos+1]-1; ++j, ++n)
      99         [ +  + ]:   94245243 :     if (col*ncomp+pos+1 == ja[j])
     100                 :   13371648 :       return a[ia[rncomp+pos]-1+n];
     101                 :            : 
     102 [ -  - ][ -  - ]:          0 :   Throw("Sparse matrix index not found");
                 [ -  - ]
     103                 :            : }
     104                 :            : 
     105                 :            : void
     106                 :      93136 : CSR::dirichlet( std::size_t i,
     107                 :            :                 const std::vector< std::size_t >& gid,
     108                 :            :                 const NodeCommMap& nodecommap,
     109                 :            :                 std::size_t pos )
     110                 :            : // *****************************************************************************
     111                 :            : //  Set Dirichlet boundary condition at a node
     112                 :            : //! \param[in] i Local id at which to set Dirichlet BC
     113                 :            : //! \param[in] gid Local->global node id map
     114                 :            : //! \param[in] nodecommap Node communication map with global node ids
     115                 :            : //! \param[in] pos Position in block
     116                 :            : //! \details In parallel there can be multiple contributions to a single node
     117                 :            : //!   on the mesh, and correspondingly, a single matrix row can be partially
     118                 :            : //!   represented on multiple partitions. Setting a Dirichlet BC entails
     119                 :            : //!   zeroing out the row of the matrix and putting 1/N into the diagonal entry,
     120                 :            : //!   where N is the number of partitions that contribute to the mesh node
     121                 :            : //!   (matrix row). As a result, when the matrix participates in a matrix-vector
     122                 :            : //!   product, where the partial contributions across all partitions are
     123                 :            : //!   aggregated, the diagonal will contain 1 after the sum across partitions.
     124                 :            : //! \note Both gid and nodecommap are optional - unused in serial. If nodecommap
     125                 :            : //!   is empty, serial is assumed and gid is unused.
     126                 :            : // *****************************************************************************
     127                 :            : {
     128                 :      93136 :   auto incomp = i * ncomp;
     129         [ +  + ]:      93136 :   auto diag = nodecommap.empty() ? 1.0 : 1.0/tk::count(nodecommap,gid[i]);
     130                 :            : 
     131                 :            :   // zero column
     132         [ +  + ]:  107268783 :   for (std::size_t r=0; r<rnz.size()*ncomp; ++r)
     133         [ +  + ]: 1333502368 :     for (std::size_t j=ia[r]-1; j<ia[r+1]-1; ++j)
     134         [ +  + ]: 1226326721 :       if (incomp+pos+1==ja[j]) a[j] = 0.0;
     135                 :            : 
     136                 :            :   // zero row and put in diagonal
     137         [ +  + ]:     998900 :   for (std::size_t j=ia[incomp+pos]-1; j<ia[incomp+pos+1]-1; ++j)
     138         [ +  + ]:     905764 :     if (incomp+pos+1==ja[j]) a[j] = diag; else a[j] = 0.0;
     139                 :      93136 : }
     140                 :            : 
     141                 :            : void
     142                 :       5465 : CSR::mult( const std::vector< real >& x,
     143                 :            :            std::vector< real >& r,
     144                 :            :            const std::vector< tk::real >& bcmask ) const
     145                 :            : // *****************************************************************************
     146                 :            : //  Multiply CSR matrix with vector from the right: r = A * x
     147                 :            : //! \param[in] x Vector to multiply matrix with from the right
     148                 :            : //! \param[in] r Result vector of product r = A * x
     149                 :            : //! \param[in] bcmask Dirichlet BC mask
     150                 :            : //! \note This is only complete in serial. In parallel, this computes the own
     151                 :            : //!   contributions to the product, so it must be followed by communication
     152                 :            : //!   combining the rows stored on multiple partitions.
     153                 :            : // *****************************************************************************
     154                 :            : {
     155         [ +  - ]:       5465 :   std::fill( begin(r), end(r), 0.0 );
     156                 :            : 
     157         [ +  + ]:   26440842 :   for (std::size_t i=0; i<rnz.size()*ncomp; ++i)
     158         [ +  + ]:  348720804 :     for (std::size_t j=ia[i]-1; j<ia[i+1]-1; ++j)
     159                 :  322285427 :       r[i] += bcmask[i] * a[j] * x[ja[j]-1];
     160                 :       5465 : }
     161                 :            : 
     162                 :            : std::ostream&
     163                 :          1 : CSR::write_stored( std::ostream& os ) const
     164                 :            : // *****************************************************************************
     165                 :            : //  Write out CSR as stored
     166                 :            : //! \param[in,out] os Output stream to write to
     167                 :            : //! \return Updated output stream
     168                 :            : // *****************************************************************************
     169                 :            : {
     170                 :          1 :   os << "size (npoin) = " << rnz.size() << '\n';
     171                 :          1 :   os << "ncomp = " << ncomp << '\n';
     172                 :          1 :   os << "rsize (size*ncomp) = " << rnz.size() * ncomp << '\n';
     173                 :          1 :   os << "nnz = " << a.size() << '\n';
     174                 :            : 
     175                 :            :   std::size_t i;
     176                 :            : 
     177                 :          1 :   os << "rnz[npoin=" << rnz.size() << "] = { ";
     178         [ +  + ]:         14 :   for (i=0; i<rnz.size()-1; ++i) os << rnz[i] << ", ";
     179                 :          1 :   os << rnz[i] << " }\n";
     180                 :            : 
     181                 :          1 :   os << "ia[rsize+1=" << rnz.size()*ncomp+1 << "] = { ";
     182         [ +  + ]:         15 :   for (i=0; i<ia.size()-1; ++i) os << ia[i] << ", ";
     183                 :          1 :   os << ia[i] << " }\n";
     184                 :            : 
     185                 :          1 :   os << "ja[nnz=" << ja.size() << "] = { ";
     186         [ +  + ]:        112 :   for (i=0; i<ja.size()-1; ++i) os << ja[i] << ", ";
     187                 :          1 :   os << ja[i] << " }\n";
     188                 :            : 
     189                 :          1 :   os << "a[nnz=" << a.size() << "] = { ";
     190         [ +  + ]:        112 :   for (i=0; i<a.size()-1; ++i) os << a[i] << ", ";
     191                 :          1 :   os << a[i] << " }\n";
     192                 :            : 
     193                 :          1 :   return os;
     194                 :            : }
     195                 :            : 
     196                 :            : std::ostream&
     197                 :          1 : CSR::write_structure( std::ostream& os ) const
     198                 :            : // *****************************************************************************
     199                 :            : //  Write out CSR nonzero structure
     200                 :            : //! \param[in,out] os Output stream to write to
     201                 :            : //! \return Updated output stream
     202                 :            : // *****************************************************************************
     203                 :            : {
     204         [ +  + ]:         15 :   for (std::size_t i=0; i<rnz.size()*ncomp; ++i) {
     205                 :            :     // leading zeros
     206         [ +  + ]:         28 :     for (std::size_t j=1; j<ja[ia[i]-1]; ++j) os << ". ";
     207         [ +  + ]:        126 :     for (std::size_t n=ia[i]-1; n<ia[i+1]-1; ++n) {
     208                 :            :       // zeros between nonzeros
     209 [ +  + ][ +  + ]:        176 :       if (n>ia[i]-1) for (std::size_t j=ja[n-1]; j<ja[n]-1; ++j) os << ". ";
     210                 :            :       // nonzero
     211                 :        112 :       os << "o ";
     212                 :            :     }
     213                 :            :     // trailing zeros
     214         [ +  + ]:         20 :     for (std::size_t j=ja[ia[i+1]-2]; j<rnz.size()*ncomp; ++j) os << ". ";
     215                 :         14 :     os << '\n';
     216                 :            :   }
     217                 :            : 
     218                 :          1 :   return os;
     219                 :            : }
     220                 :            : 
     221                 :            : std::ostream&
     222                 :          1 : CSR::write_matrix( std::ostream& os ) const
     223                 :            : // *****************************************************************************
     224                 :            : //  Write out CSR as a real matrix
     225                 :            : //! \param[in,out] os Output stream to write to
     226                 :            : //! \return Updated output stream
     227                 :            : // *****************************************************************************
     228                 :            : {
     229         [ +  + ]:         15 :   for (std::size_t i=0; i<rnz.size()*ncomp; ++i) {
     230         [ +  + ]:         28 :     for (std::size_t j=1; j<ja[ia[i]-1]; ++j) os << "0\t";
     231         [ +  + ]:        126 :     for (std::size_t n=ia[i]-1; n<ia[i+1]-1; ++n) {
     232 [ +  + ][ +  + ]:        176 :       if (n>ia[i]-1) for (std::size_t j=ja[n-1]; j<ja[n]-1; ++j ) os << "0\t";
     233                 :        112 :       os << a[n] << '\t';
     234                 :            :     }
     235         [ +  + ]:         20 :     for (std::size_t j=ja[ia[i+1]-2]; j<rnz.size()*ncomp; ++j) os << "0\t";
     236                 :         14 :     os << '\n';
     237                 :            :   }
     238                 :            : 
     239                 :          1 :   return os;
     240                 :            : }
     241                 :            : 
     242                 :            : std::ostream&
     243                 :          1 : CSR::write_matlab( std::ostream& os ) const
     244                 :            : // *****************************************************************************
     245                 :            : //  Write out CSR in Matlab/Octave format
     246                 :            : //! \param[in,out] os Output stream to write to
     247                 :            : //! \return Updated output stream
     248                 :            : // *****************************************************************************
     249                 :            : {
     250                 :          1 :   os << "A = [ ";
     251         [ +  + ]:         15 :   for (std::size_t i=0; i<rnz.size()*ncomp; ++i) {
     252         [ +  + ]:         28 :     for (std::size_t j=1; j<ja[ia[i]-1]; ++j) os << "0 ";
     253         [ +  + ]:        126 :     for ( std::size_t n=ia[i]-1; n<ia[i+1]-1; ++n) {
     254         [ +  + ]:        112 :       if (n > ia[i]-1)
     255         [ +  + ]:        162 :         for (std::size_t j=ja[n-1]; j<ja[n]-1; ++j) os << "0 ";
     256                 :        112 :       os << a[n] << ' ';
     257                 :            :     }
     258         [ +  + ]:         20 :     for (std::size_t j=ja[ia[i+1]-2]; j<rnz.size()*ncomp; ++j) os << "0 ";
     259                 :         14 :     os << ";\n";
     260                 :            :   }
     261                 :          1 :   os << "]\n";
     262                 :            : 
     263                 :          1 :   return os;
     264                 :            : }