FEI Version of the Day
Loading...
Searching...
No Matches
fei_CSRMat.cpp
1
2/*--------------------------------------------------------------------*/
3/* Copyright 2005 Sandia Corporation. */
4/* Under the terms of Contract DE-AC04-94AL85000, there is a */
5/* non-exclusive license for use of this work by or on behalf */
6/* of the U.S. Government. Export of this program may require */
7/* a license from the United States Government. */
8/*--------------------------------------------------------------------*/
9
10#include "fei_CSRMat.hpp"
11#include <fei_impl_utils.hpp>
12#include "fei_ArrayUtils.hpp"
13#include <limits>
14#include <cmath>
15
16namespace fei {
17
18CSRMat::CSRMat()
19 : srg_(),
20 packedcoefs_()
21{
22}
23
24CSRMat::CSRMat(const FillableMat& fmat)
25 : srg_(),
26 packedcoefs_()
27{
28 *this = fmat;
29}
30
31CSRMat::~CSRMat()
32{
33}
34
35CSRMat&
36CSRMat::operator=(const fei::FillableMat& src)
37{
38 FillableMat::const_iterator iter = src.begin(), iter_end = src.end();
39
40 unsigned nrows = src.getNumRows();
41
42 srg_.rowNumbers.resize(nrows);
43 srg_.rowOffsets.resize(nrows+1);
44
45 unsigned nnz = 0;
46 unsigned i = 0;
47 for(; iter != iter_end; ++iter, ++i) {
48 srg_.rowNumbers[i] = iter->first;
49 srg_.rowOffsets[i] = nnz;
50 nnz += iter->second->size();
51 }
52
53 srg_.rowOffsets[nrows] = nnz;
54
55 srg_.packedColumnIndices.resize(nnz);
56 packedcoefs_.resize(nnz);
57
58 int* colind_ptr = (srg_.packedColumnIndices.size()
59 ? &(srg_.packedColumnIndices[0]) : 0);
60 double* coef_ptr = (packedcoefs_.size()
61 ? &(packedcoefs_[0]) : 0);
62
63 iter = src.begin();
64
65 unsigned offset = 0;
66 for(; iter != iter_end; ++iter) {
67 const CSVec* v = iter->second;
68 const std::vector<int>& v_ind = v->indices();
69 const std::vector<double>& v_coef = v->coefs();
70 for(size_t i=0; i<v_ind.size(); ++i) {
71 colind_ptr[offset] = v_ind[i];
72 coef_ptr[offset++] = v_coef[i];
73 }
74 }
75
76 return *this;
77}
78
79CSRMat&
80CSRMat::operator+=(const CSRMat& src)
81{
82 FillableMat tmp;
83 add_CSRMat_to_FillableMat(*this, tmp);
84 add_CSRMat_to_FillableMat(src, tmp);
85 *this = tmp;
86 return *this;
87}
88
89bool
90CSRMat::operator==(const CSRMat& rhs) const
91{
92 if (getGraph() != rhs.getGraph()) return false;
93 return getPackedCoefs() == rhs.getPackedCoefs();
94}
95
96bool
97CSRMat::operator!=(const CSRMat& rhs) const
98{
99 return !(*this == rhs);
100}
101
102void multiply_CSRMat_CSVec(const CSRMat& A, const CSVec& x, CSVec& y)
103{
104 //This function is unit-tested in fei/utest_cases/fei_unit_CSRMat_CSVec.cpp
105
106 const std::vector<int>& rows = A.getGraph().rowNumbers;
107 const int* rowoffs = &(A.getGraph().rowOffsets[0]);
108 const std::vector<int>& colinds = A.getGraph().packedColumnIndices;
109 const double* Acoef = A.getPackedCoefs().size() > 0 ? &(A.getPackedCoefs()[0]): NULL;
110
111 const std::vector<int>& xind = x.indices();
112 const std::vector<double>& xcoef = x.coefs();
113
114 const double* xcoef_ptr = xcoef.empty() ? NULL : &xcoef[0];
115 const int* xind_ptr = xind.empty() ? NULL : &xind[0];
116 int xlen = xcoef.size();
117
118 std::vector<int>& yind = y.indices();
119 std::vector<double>& ycoef = y.coefs();
120
121 unsigned nrows = A.getNumRows();
122
123 yind.resize(nrows);
124 ycoef.resize(nrows);
125
126 int* yind_ptr = yind.size() > 0 ? &yind[0] : NULL;
127 double* ycoef_ptr = ycoef.size() > 0 ? &ycoef[0] : NULL;
128
129 int jbeg = *rowoffs++;
130 for(unsigned i=0; i<nrows; ++i) {
131 int jend = *rowoffs++;
132
133 double sum = 0.0;
134 while(jbeg<jend) {
135 int xoff = fei::binarySearch(colinds[jbeg], xind_ptr, xlen);
136
137 if (xoff > -1) {
138 sum += Acoef[jbeg]*xcoef_ptr[xoff];
139 }
140 ++jbeg;
141 }
142
143 yind_ptr[i] = rows[i];
144 ycoef_ptr[i] = sum;
145 }
146}
147
148void multiply_trans_CSRMat_CSVec(const CSRMat& A, const CSVec& x, CSVec& y)
149{
150 const std::vector<int>& rows = A.getGraph().rowNumbers;
151 const int* rowoffs = &(A.getGraph().rowOffsets[0]);
152 const int* colinds = A.getGraph().packedColumnIndices.empty() ? NULL : &(A.getGraph().packedColumnIndices[0]);
153 const double* Acoef = A.getPackedCoefs().empty() ? NULL : &(A.getPackedCoefs()[0]);
154
155 const std::vector<int>& xind = x.indices();
156 const std::vector<double>& xcoef = x.coefs();
157
158 const double* xcoef_ptr = xcoef.empty() ? NULL : &xcoef[0];
159
160 unsigned nrows = A.getNumRows();
161
162 std::vector<int> offsets;
163 fei::impl_utils::find_offsets(rows, xind, offsets);
164 const int* offsetsptr = &offsets[0];
165
166 fei::CSVec fy;
167
168 int jbeg = *rowoffs++;
169 for(unsigned i=0; i<nrows; ++i) {
170 int jend = *rowoffs++;
171
172 int xoff = offsetsptr[i];
173 if (xoff < 0) {
174 jbeg = jend;
175 continue;
176 }
177
178 double xcoeff = xcoef_ptr[xoff];
179
180 while(jbeg<jend) {
181 add_entry(fy, colinds[jbeg],Acoef[jbeg]*xcoeff);
182 ++jbeg;
183 }
184 }
185
186 y = fy;
187}
188
189void multiply_CSRMat_CSRMat(const CSRMat& A, const CSRMat& B, CSRMat& C,
190 bool storeResultZeros)
191{
192 //This function is unit-tested in fei/utest_cases/fei_unit_CSRMat_CSVec.cpp
193
194 fei::FillableMat fc;
195
196 const std::vector<int>& Arows = A.getGraph().rowNumbers;
197 const std::vector<int>& Brows = B.getGraph().rowNumbers;
198 if (Arows.size() < 1 || Brows.size() < 1) {
199 C = fc;
200 return;
201 }
202 const int* Arowoffs = &(A.getGraph().rowOffsets[0]);
203 const int* Acols = &(A.getGraph().packedColumnIndices[0]);
204 const double* Acoefs = &(A.getPackedCoefs()[0]);
205
206 const int* Browoffs = &(B.getGraph().rowOffsets[0]);
207 const std::vector<int>& Bcols = B.getGraph().packedColumnIndices;
208 const double* Bcoefs = B.getPackedCoefs().empty() ? NULL : &(B.getPackedCoefs()[0]);
209
210 static double fei_min = std::numeric_limits<double>::min();
211
212 int jbeg = *Arowoffs++;
213 for(size_t i=0; i<Arows.size(); ++i) {
214 int row = Arows[i];
215 int jend = *Arowoffs++;
216
217 fei::CSVec* fc_row = NULL;
218 if (storeResultZeros) {
219 fc_row = fc.create_or_getRow(row);
220 }
221 else {
222 fc_row = fc.hasRow(row) ? fc.create_or_getRow(row) : NULL;
223 }
224
225 while(jbeg<jend) {
226 ++jbeg;
227 int Acol = *Acols++;
228 double Acoef = *Acoefs++;
229
230 int Brow_offset = fei::binarySearch(Acol, &Brows[0], Brows.size());
231
232 if (Brow_offset < 0) {
233 continue;
234 }
235
236 if (!storeResultZeros) {
237 if (std::abs(Acoef) < fei_min) {
238 continue;
239 }
240 }
241
242 const int* Brow_cols = Bcols.empty() ? NULL : &(Bcols[Browoffs[Brow_offset]]);
243 const double* Brow_coefs = Bcoefs==NULL ? NULL : &(Bcoefs[Browoffs[Brow_offset]]);
244 int Brow_len = Browoffs[Brow_offset+1]-Browoffs[Brow_offset];
245
246 for(int k=0; k<Brow_len; ++k) {
247 double resultCoef = Acoef*Brow_coefs[k];
248 int resultCol = Brow_cols[k];
249
250 if (!storeResultZeros) {
251 if (std::abs(resultCoef) < fei_min) {
252 continue;
253 }
254 }
255
256 if (fc_row == NULL) {
257 fc_row = fc.create_or_getRow(row);
258 }
259
260 add_entry(*fc_row, resultCol, resultCoef);
261 }
262 }
263 }
264
265 C = fc;
266}
267
269 bool storeResultZeros)
270{
271 //This function is unit-tested in fei/utest_cases/fei_unit_CSRMat_CSVec.cpp
272
273 fei::FillableMat fc;
274
275 const std::vector<int>& Arows = A.getGraph().rowNumbers;
276 const std::vector<int>& Brows = B.getGraph().rowNumbers;
277 if (Arows.size() < 1 || Brows.size() < 1) {
278 C = fc;
279 return;
280 }
281
282 const size_t numArows = Arows.size();
283 const int* Arowoffs = &(A.getGraph().rowOffsets[0]);
284 const int* Acols = A.getGraph().packedColumnIndices.empty() ? NULL : &(A.getGraph().packedColumnIndices[0]);
285 const double* Acoefs = A.getPackedCoefs().empty() ? NULL : &(A.getPackedCoefs()[0]);
286
287 const int* Browoffs = &(B.getGraph().rowOffsets[0]);
288 const std::vector<int>& Bcols = B.getGraph().packedColumnIndices;
289 const double* Bcoefs = B.getPackedCoefs().empty() ? NULL : &(B.getPackedCoefs()[0]);
290
291 std::vector<double> row_coefs;
292
293 static double fei_min = std::numeric_limits<double>::min();
294
295 std::vector<int> offsets;
296 fei::impl_utils::find_offsets(Arows, Brows, offsets);
297
298 int jbeg = *Arowoffs++;
299 for(size_t i=0; i<numArows; ++i) {
300 int jend = *Arowoffs++;
301
302 int Brow_offset = offsets[i];
303 if (Brow_offset < 0) {
304 jbeg = jend;
305 continue;
306 }
307
308 const int* Brow_cols = Bcols.empty() ? NULL : &(Bcols[Browoffs[Brow_offset]]);
309 const double* Brow_coefs = Bcoefs==NULL ? NULL : &(Bcoefs[Browoffs[Brow_offset]]);
310 int Brow_len = Browoffs[Brow_offset+1]-Browoffs[Brow_offset];
311
312 if ((int)row_coefs.size() < Brow_len) row_coefs.resize(Brow_len*2);
313 double* row_coefs_ptr = &row_coefs[0];
314
315 while(jbeg<jend) {
316 int Acol = Acols[jbeg];
317 double Acoef = Acoefs[jbeg++];
318
319 if (std::abs(Acoef) < fei_min && !storeResultZeros) {
320 continue;
321 }
322
323 for(int k=0; k<Brow_len; ++k) {
324 row_coefs_ptr[k] = Acoef*Brow_coefs[k];
325 }
326
327 fc.sumInRow(Acol, Brow_cols, row_coefs_ptr, Brow_len);
328 }
329 }
330
331 C = fc;
332}
333
334void add_CSRMat_to_FillableMat(const CSRMat& csrm, FillableMat& fm)
335{
336 const std::vector<int>& rows = csrm.getGraph().rowNumbers;
337 const int* rowoffs = &(csrm.getGraph().rowOffsets[0]);
338 const std::vector<int>& cols = csrm.getGraph().packedColumnIndices;
339 const double* coefs = &(csrm.getPackedCoefs()[0]);
340
341 for(size_t i=0; i<rows.size(); ++i) {
342 int row = rows[i];
343
344 for(int j=rowoffs[i]; j<rowoffs[i+1]; ++j) {
345 fm.sumInCoef(row, cols[j], coefs[j]);
346 }
347 }
348}
349
350}//namespace fei
351
std::vector< int > rowNumbers
std::vector< int > packedColumnIndices
std::vector< int > rowOffsets
void find_offsets(const std::vector< int > &sources, const std::vector< int > &targets, std::vector< int > &offsets)
void multiply_CSRMat_CSRMat(const CSRMat &A, const CSRMat &B, CSRMat &C, bool storeResultZeros)
void multiply_trans_CSRMat_CSRMat(const CSRMat &A, const CSRMat &B, CSRMat &C, bool storeResultZeros)
void multiply_trans_CSRMat_CSVec(const CSRMat &A, const CSVec &x, CSVec &y)
int binarySearch(const T &item, const T *list, int len)
void multiply_CSRMat_CSVec(const CSRMat &A, const CSVec &x, CSVec &y)