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