FEI  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
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 
16 namespace fei {
17 
18 CSRMat::CSRMat()
19  : srg_(),
20  packedcoefs_()
21 {
22 }
23 
24 CSRMat::CSRMat(const FillableMat& fmat)
25  : srg_(),
26  packedcoefs_()
27 {
28  *this = fmat;
29 }
30 
31 CSRMat::~CSRMat()
32 {
33 }
34 
35 CSRMat&
36 CSRMat::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 
79 CSRMat&
80 CSRMat::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 
89 bool
90 CSRMat::operator==(const CSRMat& rhs) const
91 {
92  if (getGraph() != rhs.getGraph()) return false;
93  return getPackedCoefs() == rhs.getPackedCoefs();
94 }
95 
96 bool
97 CSRMat::operator!=(const CSRMat& rhs) const
98 {
99  return !(*this == rhs);
100 }
101 
102 void 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 
148 void 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 
189 void 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 
268 void multiply_trans_CSRMat_CSRMat(const CSRMat& A, const CSRMat& B, CSRMat& C,
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 
334 void 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 
void multiply_CSRMat_CSRMat(const CSRMat &A, const CSRMat &B, CSRMat &C, bool storeResultZeros)
Definition: fei_CSRMat.cpp:189
std::vector< int > rowNumbers
void multiply_trans_CSRMat_CSVec(const CSRMat &A, const CSVec &x, CSVec &y)
Definition: fei_CSRMat.cpp:148
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)
int binarySearch(const T &item, const T *list, int len)
void multiply_trans_CSRMat_CSRMat(const CSRMat &A, const CSRMat &B, CSRMat &C, bool storeResultZeros)
Definition: fei_CSRMat.cpp:268
void multiply_CSRMat_CSVec(const CSRMat &A, const CSVec &x, CSVec &y)
Definition: fei_CSRMat.cpp:102