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 : 16 : CSR::CSR( std::size_t nc,
19 : : const std::pair< std::vector< std::size_t >,
20 : 16 : std::vector< std::size_t > >& psup )
21 : : try :
22 : : ncomp( nc ),
23 : 16 : rnz( psup.second.size()-1 ),
24 [ + + ][ + - ]: 37 : 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 [ + + ][ + - ]: 15 : Assert( ncomp > 0, "Sparse matrix ncomp must be positive" );
[ + - ][ + - ]
33 [ - + ][ - - ]: 14 : Assert( rnz.size() > 0, "Sparse matrix size must be positive" );
[ - - ][ - - ]
34 : :
35 : 14 : const auto& psup1 = psup.first;
36 : 14 : 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 [ + + ]: 200 : 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 [ + + ]: 1426 : for (rnz[i]=1, j=psup2[i]+1; j<=psup2[i+1]; ++j)
45 : 1240 : ++rnz[i];
46 : :
47 : : // add up total number of nonzeros
48 : 186 : nnz += rnz[i] * ncomp;
49 : :
50 : : // fill up row index
51 [ + + ]: 474 : for (std::size_t k=0; k<ncomp; ++k)
52 : 288 : ia[i*ncomp+k+1] = ia[i*ncomp+k] + rnz[i];
53 : : }
54 : :
55 : : // Allocate storage for matrix values and column indices
56 [ + - ]: 14 : a.resize( nnz, 0.0 );
57 [ + - ]: 14 : ja.resize( nnz );
58 : :
59 : : // fill column indices
60 [ + + ]: 200 : for (i=0; i<rnz.size(); ++i)
61 [ + + ]: 474 : for (std::size_t k=0; k<ncomp; ++k) {
62 : 288 : auto itmp = i*ncomp+k;
63 : 288 : ja[ia[itmp]-1] = itmp+1; // put in column index of diagonal
64 [ + + ]: 2180 : 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 : 1892 : ja[ia[itmp]-1+(n++)] = psup1[j]*ncomp+k+1;
67 : : }
68 : : }
69 : :
70 : : // (bubble-)sort column indices
71 [ + + ]: 200 : for (i=0; i<rnz.size(); ++i)
72 [ + + ]: 474 : for (std::size_t k=0; k<ncomp; ++k)
73 [ + + ]: 2180 : for (std::size_t j=psup2[i]+1; j<=psup2[i+1]; ++j)
74 [ + + ]: 14936 : for (std::size_t l=1; l<rnz[i]; ++l) // sort column indices of row i
75 [ + + ]: 66352 : for (std::size_t e=0; e<rnz[i]-l; ++e)
76 [ + + ]: 53308 : if (ja[ia[i*ncomp+k]-1+e] > ja[ia[i*ncomp+k]+e])
77 : 946 : 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 : 14 : }
84 : :
85 : : const tk::real&
86 : 10560 : 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 : 10560 : auto rncomp = row * ncomp;
97 : :
98 [ + - ]: 51243 : for (std::size_t n=0, j=ia[rncomp+pos]-1; j<ia[rncomp+pos+1]-1; ++j, ++n)
99 [ + + ]: 51243 : if (col*ncomp+pos+1 == ja[j])
100 : 10560 : return a[ia[rncomp+pos]-1+n];
101 : :
102 [ - - ][ - - ]: 0 : Throw("Sparse matrix index not found");
[ - - ]
103 : : }
104 : :
105 : : void
106 : 12 : 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 : 12 : auto incomp = i * ncomp;
129 [ + + ]: 12 : auto diag = nodecommap.empty() ? 1.0 : 1.0/tk::count(nodecommap,gid[i]);
130 : :
131 : : // zero column
132 [ + + ]: 382 : for (std::size_t r=0; r<rnz.size()*ncomp; ++r)
133 [ + + ]: 3020 : for (std::size_t j=ia[r]-1; j<ia[r+1]-1; ++j)
134 [ + + ]: 2650 : if (incomp+pos+1==ja[j]) a[j] = 0.0;
135 : :
136 : : // zero row and put in diagonal
137 [ + + ]: 80 : for (std::size_t j=ia[incomp+pos]-1; j<ia[incomp+pos+1]-1; ++j)
138 [ + + ]: 68 : if (incomp+pos+1==ja[j]) a[j] = diag; else a[j] = 0.0;
139 : 12 : }
140 : :
141 : : void
142 : 43 : 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 [ + - ]: 43 : std::fill( begin(r), end(r), 0.0 );
156 : :
157 [ + + ]: 1093 : for (std::size_t i=0; i<rnz.size()*ncomp; ++i)
158 [ + + ]: 8582 : for (std::size_t j=ia[i]-1; j<ia[i+1]-1; ++j)
159 : 7532 : r[i] += bcmask[i] * a[j] * x[ja[j]-1];
160 : 43 : }
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 : : }
|
