FEI Version of the Day
Loading...
Searching...
No Matches
fei_impl_utils.cpp
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
19namespace fei {
20namespace impl_utils {
21
22//----------------------------------------------------------------------------
23void 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//----------------------------------------------------------------------------
52size_t num_bytes_FillableMat(const fei::FillableMat& mat)
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//----------------------------------------------------------------------------
64void pack_FillableMat(const fei::FillableMat& mat,
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
83 fei::FillableMat::const_iterator
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//----------------------------------------------------------------------------
105void 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//----------------------------------------------------------------------------
153bool 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
198size_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
207void 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
235void 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//----------------------------------------------------------------------------
258void separate_BC_eqns(const fei::FillableMat& mat,
259 std::vector<int>& bcEqns,
260 std::vector<double>& bcVals)
261{
262 fei::FillableMat::const_iterator
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//----------------------------------------------------------------------------
290void create_col_to_row_map(const fei::FillableMat& mat,
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
300 fei::FillableMat::const_iterator
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//----------------------------------------------------------------------------
319int remove_couplings(fei::FillableMat& mat)
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
333 fei::FillableMat::iterator
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//----------------------------------------------------------------------------
385void global_union(MPI_Comm comm,
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//----------------------------------------------------------------------------
434void global_union(MPI_Comm comm,
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//----------------------------------------------------------------------------
488void translate_to_reduced_eqns(const fei::Reducer& reducer, fei::CSRMat& mat)
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//----------------------------------------------------------------------------
504void translate_to_reduced_eqns(const fei::Reducer& reducer, fei::CSVec& vec)
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//----------------------------------------------------------------------------
513void 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//----------------------------------------------------------------------------
532void 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
virtual int addIndices(int row, int len, const int *indices)=0
virtual int copyIn(int numRows, const int *rows, int numCols, const int *cols, const double *const *values, int format=0)=0
virtual int sumIn(int numRows, const int *rows, int numCols, const int *cols, const double *const *values, int format=0)=0
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)
bool unpack_CSRMat(const char *buffer_begin, const char *buffer_end, fei::CSRMat &mat)
void unpack_FillableMat(const char *buffer_begin, const char *buffer_end, fei::FillableMat &mat, bool clear_mat_on_entry, bool overwrite_entries)
int count_nnz(const FillableMat &mat)
int localProc(MPI_Comm comm)
int Allgatherv(MPI_Comm comm, std::vector< T > &sendbuf, std::vector< int > &recvLengths, std::vector< T > &recvbuf)
std::ostream & console_out()
int numProcs(MPI_Comm comm)
void get_row_numbers(const FillableMat &mat, std::vector< int > &rows)