FEI Package Browser (Single Doxygen Collection)  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
fei_impl_utils.cpp
Go to the documentation of this file.
1 /*--------------------------------------------------------------------*/
2 /* Copyright 2008 Sandia Corporation. */
3 /* Under the terms of Contract DE-AC04-94AL85000, there is a */
4 /* non-exclusive license for use of this work by or on behalf */
5 /* of the U.S. Government. Export of this program may require */
6 /* a license from the United States Government. */
7 /*--------------------------------------------------------------------*/
8 
9 #include <fei_CommUtils.hpp>
10 #include <fei_iostream.hpp>
11 #include <fei_impl_utils.hpp>
12 #include <fei_FillableMat.hpp>
13 #include <fei_CSRMat.hpp>
14 #include <fei_CSVec.hpp>
15 #include <fei_Graph.hpp>
16 #include <fei_Matrix.hpp>
17 #include <fei_Reducer.hpp>
18 
19 namespace fei {
20 namespace impl_utils {
21 
22 //----------------------------------------------------------------------------
23 void find_offsets(const std::vector<int>& sources,
24  const std::vector<int>& targets,
25  std::vector<int>& offsets)
26 {
27  offsets.assign(sources.size(), -1);
28 
29  size_t ssize = sources.size(), tsize = targets.size();
30  size_t si = 0, ti = 0;
31 
32  const int* sourcesptr = &sources[0];
33  const int* targetsptr = &targets[0];
34 
35  while(si<ssize && ti<tsize) {
36  if (sourcesptr[si] == targetsptr[ti]) {
37  offsets[si++] = ti++;
38  continue;
39  }
40 
41  while(sourcesptr[si] < targetsptr[ti] && si<ssize) {
42  ++si;
43  }
44 
45  while(sourcesptr[si] > targetsptr[ti] && ti<tsize) {
46  ++ti;
47  }
48  }
49 }
50 
51 //----------------------------------------------------------------------------
53 {
54  int nrows = mat.getNumRows();
55  int nnz = fei::count_nnz(mat);
56 
57  int num_chars_int = (2 + nrows*2 + nnz)*sizeof(int);
58  int num_chars_double = nnz*sizeof(double);
59 
60  return num_chars_int + num_chars_double;
61 }
62 
63 //----------------------------------------------------------------------------
65  char* buffer)
66 {
67  int nrows = mat.getNumRows();
68  int nnz = fei::count_nnz(mat);
69 
70  int num_chars_int = (2 + nrows*2 + nnz)*sizeof(int);
71 
72  int* intdata = reinterpret_cast<int*>(buffer);
73  double* doubledata = reinterpret_cast<double*>(buffer+num_chars_int);
74 
75  int ioffset = 0;
76  int doffset = 0;
77 
78  intdata[ioffset++] = nrows;
79  intdata[ioffset++] = nnz;
80 
81  int ioffsetcols = 2+nrows*2;
82 
84  r_iter = mat.begin(),
85  r_end = mat.end();
86 
87  for(; r_iter!=r_end; ++r_iter) {
88  int rowNumber = r_iter->first;
89  const fei::CSVec* row = r_iter->second;
90 
91  intdata[ioffset++] = rowNumber;
92  const int rowlen = row->size();
93  intdata[ioffset++] = rowlen;
94 
95  const std::vector<int>& rowindices = row->indices();
96  const std::vector<double>& rowcoefs = row->coefs();
97  for(int i=0; i<rowlen; ++i) {
98  intdata[ioffsetcols++] = rowindices[i];
99  doubledata[doffset++] = rowcoefs[i];
100  }
101  }
102 }
103 
104 //----------------------------------------------------------------------------
105 void unpack_FillableMat(const char* buffer_begin, const char* buffer_end,
106  fei::FillableMat& mat,
107  bool clear_mat_on_entry,
108  bool overwrite_entries)
109 {
110  if (clear_mat_on_entry) {
111  mat.clear();
112  }
113 
114  if (buffer_end == buffer_begin) {
115  return;
116  }
117 
118  const int* intdata = reinterpret_cast<const int*>(buffer_begin);
119  int ioffset = 0;
120  int nrows = intdata[ioffset++];
121  int nnz = intdata[ioffset++];
122 
123  int ioffsetcols = 2+nrows*2;
124 
125  int num_chars_int = (2+nrows*2 + nnz)*sizeof(int);
126  const double* doubledata = reinterpret_cast<const double*>(buffer_begin+num_chars_int);
127 
128  int doffset = 0;
129 
130  for(int i=0; i<nrows; ++i) {
131  int row = intdata[ioffset++];
132  int rowlen = intdata[ioffset++];
133 
134  for(int j=0; j<rowlen; ++j) {
135  int col = intdata[ioffsetcols++];
136  double coef = doubledata[doffset++];
137 
138  if (overwrite_entries) {
139  mat.putCoef(row, col, coef);
140  }
141  else {
142  mat.sumInCoef(row, col, coef);
143  }
144  }
145  }
146 
147  if (doffset != nnz) {
148  throw std::runtime_error("fei::impl_utils::unpack_FillableMat: failed, sizes don't agree.");
149  }
150 }
151 
152 //----------------------------------------------------------------------------
153 bool unpack_CSRMat(const char* buffer_begin, const char* buffer_end, fei::CSRMat& mat)
154 {
155  bool all_zeros = true;
156  if (buffer_end == buffer_begin) {
157  return all_zeros;
158  }
159 
160  const int* intdata = reinterpret_cast<const int*>(buffer_begin);
161  int ioffset = 0;
162  int nrows = intdata[ioffset++];
163  int nnz = intdata[ioffset++];
164 
165  fei::SparseRowGraph& srg = mat.getGraph();
166  srg.rowNumbers.resize(nrows);
167  srg.rowOffsets.resize(nrows+1);
168  srg.packedColumnIndices.resize(nnz);
169  std::vector<double>& packed_coefs = mat.getPackedCoefs();
170  packed_coefs.resize(nnz);
171 
172  int ioffsetcols = 2+nrows*2;
173 
174  int num_chars_int = (2+nrows*2 + nnz)*sizeof(int);
175  const double* doubledata = reinterpret_cast<const double*>(buffer_begin+num_chars_int);
176 
177  int doffset = 0;
178 
179  for(int i=0; i<nrows; ++i) {
180  int row = intdata[ioffset++];
181  int rowlen = intdata[ioffset++];
182 
183  srg.rowNumbers[i] = row;
184  srg.rowOffsets[i] = doffset;
185 
186  for(int j=0; j<rowlen; ++j) {
187  int col = intdata[ioffsetcols++];
188  double coef = doubledata[doffset];
189  if (coef != 0.0) all_zeros = false;
190  srg.packedColumnIndices[doffset] = col;
191  packed_coefs[doffset++] = coef;
192  }
193  }
194  srg.rowOffsets[nrows] = nnz;
195  return all_zeros;
196 }
197 
198 size_t num_bytes_indices_coefs(const std::vector<int>& indices,
199  const std::vector<double>& coefs)
200 {
201  int num = indices.size();
202  int num_chars_int = (1+num)*sizeof(int);
203  int num_chars = num_chars_int + num*sizeof(double);
204  return num_chars;
205 }
206 
207 void pack_indices_coefs(const std::vector<int>& indices,
208  const std::vector<double>& coefs,
209  std::vector<char>& buffer,
210  bool resize_buffer)
211 {
212  if (indices.size() != coefs.size()) {
213  throw std::runtime_error("fei::impl_utils::pack_indices_coefs failed, sizes don't match.");
214  }
215 
216  int num = indices.size();
217  int num_chars_int = (1+num)*sizeof(int);
218  int num_chars = num_chars_int + num*sizeof(double);
219  if (resize_buffer) {
220  buffer.resize(num_chars);
221  }
222 
223  int* intdata = reinterpret_cast<int*>(&buffer[0]);
224  double* doubledata = reinterpret_cast<double*>(&buffer[0]+num_chars_int);
225 
226  int ioffset = 0;
227  int doffset = 0;
228  intdata[ioffset++] = num;
229  for(int i=0; i<num; ++i) {
230  intdata[ioffset++] = indices[i];
231  doubledata[doffset++] = coefs[i];
232  }
233 }
234 
235 void unpack_indices_coefs(const std::vector<char>& buffer,
236  std::vector<int>& indices,
237  std::vector<double>& coefs)
238 {
239  if (buffer.size() == 0) return;
240 
241  const int* intdata = reinterpret_cast<const int*>(&buffer[0]);
242  int ioffset = 0;
243  int num = intdata[ioffset++];
244  int num_chars_int = (1+num)*sizeof(int);
245  const double* doubledata = reinterpret_cast<const double*>(&buffer[0]+num_chars_int);
246 
247  indices.resize(num);
248  coefs.resize(num);
249 
250  int doffset = 0;
251  for(int i=0; i<num; ++i) {
252  indices[i] = intdata[ioffset++];
253  coefs[i] = doubledata[doffset++];
254  }
255 }
256 
257 //----------------------------------------------------------------------------
259  std::vector<int>& bcEqns,
260  std::vector<double>& bcVals)
261 {
263  m_iter = mat.begin(),
264  m_end = mat.end();
265 
266  for(; m_iter != m_end; ++m_iter) {
267  int eqn = m_iter->first;
268  const fei::CSVec* row = m_iter->second;
269 
270  std::vector<int>::iterator
271  iter = std::lower_bound(bcEqns.begin(), bcEqns.end(), eqn);
272  if (iter == bcEqns.end() || *iter != eqn) {
273  size_t offset = iter - bcEqns.begin();
274  bcEqns.insert(iter, eqn);
275  std::vector<double>::iterator viter = bcVals.begin();
276  viter += offset;
277  try {
278  double val = get_entry(*row, eqn);
279  bcVals.insert(viter, val);
280  }
281  catch(...) {
282  fei::console_out() << "separate_BC_eqns ERROR, failed to find coef for eqn " << eqn;
283  throw;
284  }
285  }
286  }
287 }
288 
289 //----------------------------------------------------------------------------
291  std::multimap<int,int>& crmap)
292 {
293  crmap.clear();
294 
295  if (mat.getNumRows() == 0) return;
296 
297  std::vector<int> rowNumbers;
298  fei::get_row_numbers(mat, rowNumbers);
299 
301  m_iter = mat.begin(),
302  m_end = mat.end();
303 
304  for(; m_iter != m_end; ++m_iter) {
305  int rowNum = m_iter->first;
306  const fei::CSVec* rowvec = m_iter->second;
307 
308  const std::vector<int>& rowindices = rowvec->indices();
309 
310  for(size_t i=0; i<rowindices.size(); ++i) {
311  int colNum = rowindices[i];
312 
313  crmap.insert(std::make_pair(colNum, rowNum));
314  }
315  }
316 }
317 
318 //----------------------------------------------------------------------------
320 {
321  int levelsOfCoupling = 0;
322 
323  std::vector<int> rowNumbers;
324  fei::get_row_numbers(mat, rowNumbers);
325 
326  bool finished = false;
327  while(!finished) {
328  std::multimap<int,int> crmap;
329  create_col_to_row_map(mat, crmap);
330 
331  typedef std::multimap<int,int>::iterator MM_Iter;
332 
334  m_iter = mat.begin(),
335  m_end = mat.end();
336 
337  bool foundCoupling = false;
338  for(; m_iter != m_end; ++m_iter) {
339  int rownum = m_iter->first;
340  fei::CSVec* mrow = m_iter->second;
341 
342  //now find which rows contain 'rownum' as a column-index:
343  std::pair<MM_Iter,MM_Iter> mmi = crmap.equal_range(rownum);
344 
345  //loop over the rows which contain 'rownum' as a column-index:
346  for(MM_Iter cri = mmi.first; cri != mmi.second; ++cri) {
347  int cri_row = cri->second;
348 
349  fei::CSVec* frow = mat.create_or_getRow(cri_row);
350 
351  double coef = get_entry(*frow,rownum);
352 
353  remove_entry(*frow, rownum);
354 
355  std::vector<int>& indices = mrow->indices();
356  std::vector<double>& coefs = mrow->coefs();
357 
358  size_t rowlen = mrow->size();
359  for(size_t ii=0; ii<rowlen; ++ii) {
360  coefs[ii] *= coef;
361  }
362 
363  int* indPtr = indices.empty() ? NULL : &indices[0];
364  double* coefPtr = coefs.empty() ? NULL : &coefs[0];
365  add_entries(*frow, rowlen, indPtr, coefPtr);
366  foundCoupling = true;
367  }
368  }
369 
370  if (foundCoupling == true) ++levelsOfCoupling;
371  else finished = true;
372  }
373 
374  if (levelsOfCoupling > 1) {
375  fei::console_out() << "fei::removeCouplings WARNING, levelsOfCoupling="
376  << levelsOfCoupling << " (Each level of coupling means that a slave DOF "
377  << "depends on a master which is itself a slave.) Behavior is not well "
378  << "understood for levelsOfCoupling greater than 1..."<<FEI_ENDL;
379  }
380 
381  return levelsOfCoupling;
382 }
383 
384 //----------------------------------------------------------------------------
386  const fei::FillableMat& localMatrix,
387  fei::FillableMat& globalUnionMatrix)
388 {
389  globalUnionMatrix = localMatrix;
390 
391  int localProc = fei::localProc(comm);
392  int numProcs = fei::numProcs(comm);
393 
394  if (numProcs < 2) {
395  return;
396  }
397 
398  //first pack the local matrix into a pair of std::vector objects
399 
400  size_t num_bytes = num_bytes_FillableMat(localMatrix);
401  std::vector<char> localchardata(num_bytes);
402 
403  pack_FillableMat(localMatrix, &localchardata[0]);
404 
405  //next use Allgatherv to place every processor's packed arrays onto every
406  //other processor.
407 
408  std::vector<int> recvdatalengths;
409  std::vector<char> recvdata;
410  int err = fei::Allgatherv(comm, localchardata, recvdatalengths, recvdata);
411  if (err != 0) {
412  throw std::runtime_error("fei::impl_utils::global_union: Allgatherv-int failed.");
413  }
414 
415  //finally unpack the received arrays into matrix objects and combine them
416  //into the globalUnionMatrix object.
417  bool clearMatOnEntry = false;
418  bool overwriteEntries = true;
419 
420  int ioffset = 0;
421  for(size_t p=0; p<recvdatalengths.size(); ++p) {
422  int len = recvdatalengths[p];
423 
424  if (len > 1 && (int)p != localProc) {
425  unpack_FillableMat(&recvdata[ioffset], &recvdata[ioffset]+len,
426  globalUnionMatrix, clearMatOnEntry, overwriteEntries);
427  }
428 
429  ioffset += len;
430  }
431 }
432 
433 //----------------------------------------------------------------------------
435  const fei::CSVec& localVec, fei::CSVec& globalUnionVec)
436 {
437  globalUnionVec.indices().clear();
438  globalUnionVec.coefs().clear();
439 
440  const std::vector<int>& localindices = localVec.indices();
441  const std::vector<double>& localcoefs = localVec.coefs();
442 
443  std::vector<int> localintdata;
444  if (localindices.size() > 0) {
445  localintdata.assign(&localindices[0], &localindices[0]+localindices.size());
446  }
447 
448  std::vector<double> localdoubledata;
449  if (localcoefs.size() > 0) {
450  localdoubledata.assign(&localcoefs[0], &localcoefs[0]+localcoefs.size());
451  }
452 
453  //use Allgatherv to place every processor's arrays onto every
454  //other processor.
455 
456  std::vector<int> recvintdatalengths;
457  std::vector<int> recvintdata;
458  int err = fei::Allgatherv(comm, localintdata, recvintdatalengths, recvintdata);
459  if (err != 0) {
460  throw std::runtime_error("snl_fei::globalUnion csvec: Allgatherv-int failed.");
461  }
462 
463  std::vector<int> recvdoubledatalengths;
464  std::vector<double> recvdoubledata;
465  err = fei::Allgatherv(comm, localdoubledata, recvdoubledatalengths, recvdoubledata);
466  if (err != 0) {
467  throw std::runtime_error("snl_fei::globalUnion csvec: Allgatherv-double failed.");
468  }
469 
470  if (recvintdatalengths.size() != recvdoubledatalengths.size()) {
471  throw std::runtime_error("snl_fei::globalUnion csvec: inconsistent lengths from Allgatherv");
472  }
473 
474  //now unpack the received arrays into the globalUnionVec object.
475 
476  unsigned len = recvintdata.size();
477  if (len > 0) {
478  int* recvintPtr = &recvintdata[0];
479  double* recvdoublePtr = &recvdoubledata[0];
480 
481  for(unsigned i=0; i<len; ++i) {
482  fei::put_entry(globalUnionVec, recvintPtr[i], recvdoublePtr[i]);
483  }
484  }
485 }
486 
487 //----------------------------------------------------------------------------
489 {
490  fei::SparseRowGraph& srg = mat.getGraph();
491 
492  std::vector<int>& rowNumbers = srg.rowNumbers;
493  for(size_t i=0; i<rowNumbers.size(); ++i) {
494  rowNumbers[i] = reducer.translateToReducedEqn(rowNumbers[i]);
495  }
496 
497  std::vector<int>& colIndices = srg.packedColumnIndices;
498  for(size_t i=0; i<colIndices.size(); ++i) {
499  colIndices[i] = reducer.translateToReducedEqn(colIndices[i]);
500  }
501 }
502 
503 //----------------------------------------------------------------------------
505 {
506  std::vector<int>& indices = vec.indices();
507  for(size_t i=0; i<indices.size(); ++i) {
508  indices[i] = reducer.translateToReducedEqn(indices[i]);
509  }
510 }
511 
512 //----------------------------------------------------------------------------
513 void add_to_graph(const fei::CSRMat& inmat, fei::Graph& graph)
514 {
515  const std::vector<int>& rowNumbers = inmat.getGraph().rowNumbers;
516  const std::vector<int>& rowOffsets = inmat.getGraph().rowOffsets;
517  const std::vector<int>& pckColInds = inmat.getGraph().packedColumnIndices;
518 
519  for(size_t i=0; i<rowNumbers.size(); ++i) {
520  int row = rowNumbers[i];
521  int offset = rowOffsets[i];
522  int rowlen = rowOffsets[i+1]-offset;
523  const int* indices = pckColInds.empty() ? NULL : &pckColInds[offset];
524 
525  if (graph.addIndices(row, rowlen, indices) != 0) {
526  throw std::runtime_error("fei::impl_utils::add_to_graph ERROR in graph.addIndices.");
527  }
528  }
529 }
530 
531 //----------------------------------------------------------------------------
532 void add_to_matrix(const fei::CSRMat& inmat, bool sum_into, fei::Matrix& matrix)
533 {
534  const std::vector<int>& rowNumbers = inmat.getGraph().rowNumbers;
535  const std::vector<int>& rowOffsets = inmat.getGraph().rowOffsets;
536  const std::vector<int>& pckColInds = inmat.getGraph().packedColumnIndices;
537  const std::vector<double>& pckCoefs = inmat.getPackedCoefs();
538 
539  for(size_t i=0; i<rowNumbers.size(); ++i) {
540  int row = rowNumbers[i];
541  int offset = rowOffsets[i];
542  int rowlen = rowOffsets[i+1]-offset;
543  const int* indices = &pckColInds[offset];
544  const double* coefs = &pckCoefs[offset];
545 
546  if (sum_into) {
547  if (matrix.sumIn(1, &row, rowlen, indices, &coefs, FEI_DENSE_ROW) != 0) {
548  throw std::runtime_error("fei::impl_utils::add_to_matrix ERROR in matrix.sumIn.");
549  }
550  }
551  else {
552  if (matrix.copyIn(1, &row, rowlen, indices, &coefs, FEI_DENSE_ROW) != 0) {
553  throw std::runtime_error("fei::impl_utils::add_to_matrix ERROR in matrix.copyIn.");
554  }
555  }
556  }
557 }
558 
559 }//namespace impl_utils
560 }//namespace fei
561 
void unpack_indices_coefs(const std::vector< char > &buffer, std::vector< int > &indices, std::vector< double > &coefs)
void remove_entry(CSVec &vec, int eqn)
Definition: fei_CSVec.cpp:77
int Allgatherv(MPI_Comm comm, std::vector< T > &sendbuf, std::vector< int > &recvLengths, std::vector< T > &recvbuf)
feipoolmat::iterator iterator
void translate_to_reduced_eqns(const fei::Reducer &reducer, fei::CSRMat &mat)
void add_entries(CSVec &vec, int num, const int *eqns, const double *coefs)
Definition: fei_CSVec.cpp:34
CSVec * create_or_getRow(int row)
bool unpack_CSRMat(const char *buffer_begin, const char *buffer_end, fei::CSRMat &mat)
std::vector< int > & indices()
Definition: fei_CSVec.hpp:31
int translateToReducedEqn(int unreducedEqn) const
std::vector< int > rowNumbers
size_t size() const
Definition: fei_CSVec.hpp:36
void add_to_matrix(const fei::CSRMat &inmat, bool sum_into, fei::Matrix &matrix)
void get_row_numbers(const FillableMat &mat, std::vector< int > &rows)
#define MPI_Comm
Definition: fei_mpi.h:56
std::vector< int > packedColumnIndices
void pack_indices_coefs(const std::vector< int > &indices, const std::vector< double > &coefs, std::vector< char > &buffer, bool resize_buffer)
std::vector< int > rowOffsets
void put_entry(CSVec &vec, int eqn, double coef)
Definition: fei_CSVec.cpp:39
void find_offsets(const std::vector< int > &sources, const std::vector< int > &targets, std::vector< int > &offsets)
virtual int addIndices(int row, int len, const int *indices)=0
void global_union(MPI_Comm comm, const fei::FillableMat &localMatrix, fei::FillableMat &globalUnionMatrix)
void unpack_FillableMat(const char *buffer_begin, const char *buffer_end, fei::FillableMat &mat, bool clear_mat_on_entry, bool overwrite_entries)
virtual int sumIn(int numRows, const int *rows, int numCols, const int *cols, const double *const *values, int format=0)=0
void pack_FillableMat(const fei::FillableMat &mat, char *buffer)
void putCoef(int row, int col, double coef)
size_t num_bytes_indices_coefs(const std::vector< int > &indices, const std::vector< double > &coefs)
void add_to_graph(const fei::CSRMat &inmat, fei::Graph &graph)
feipoolmat::const_iterator const_iterator
double get_entry(const CSVec &vec, int eqn)
Definition: fei_CSVec.cpp:58
#define FEI_ENDL
SparseRowGraph & getGraph()
Definition: fei_CSRMat.hpp:27
std::ostream & console_out()
size_t num_bytes_FillableMat(const fei::FillableMat &mat)
void separate_BC_eqns(const fei::FillableMat &mat, std::vector< int > &bcEqns, std::vector< double > &bcVals)
int localProc(MPI_Comm comm)
void sumInCoef(int row, int col, double coef)
virtual int copyIn(int numRows, const int *rows, int numCols, const int *cols, const double *const *values, int format=0)=0
unsigned getNumRows() const
std::vector< double > & getPackedCoefs()
Definition: fei_CSRMat.hpp:30
#define FEI_DENSE_ROW
Definition: fei_defs.h:77
void create_col_to_row_map(const fei::FillableMat &mat, std::multimap< int, int > &crmap)
int count_nnz(const FillableMat &mat)
int numProcs(MPI_Comm comm)
std::vector< double > & coefs()
Definition: fei_CSVec.hpp:33
int remove_couplings(fei::FillableMat &mat)