FEI Version of the Day
Loading...
Searching...
No Matches
fei_EqnBuffer.cpp
1/*--------------------------------------------------------------------*/
2/* Copyright 2005 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_macros.hpp>
10
11#include <fei_EqnBuffer.hpp>
12#include <fei_CSVec.hpp>
13
14#include <fei_TemplateUtils.hpp>
15
16//==============================================================================
18 : newCoefData_(0),
19 newRHSData_(0),
20 eqnNumbers_(0),
21 eqns_(),
22 indices_union_(0),
23 numRHSs_(1),
24 rhsCoefs_(),
25 setNumRHSsCalled_(false),
26 rhsCoefsAllocated_(false),
27 dummyCoefs_()
28{
29}
30
31//==============================================================================
33 : newCoefData_(0),
34 newRHSData_(0),
35 eqnNumbers_(0),
36 eqns_(),
37 indices_union_(0),
38 numRHSs_(1),
39 rhsCoefs_(),
40 setNumRHSsCalled_(false),
41 rhsCoefsAllocated_(false),
42 dummyCoefs_()
43{
44 *this = src;
45}
46
47//==============================================================================
49{
50 int i, len = src.eqnNumbers_.size();
51
52 eqnNumbers_ = src.eqnNumbers_;
53 eqns_.resize(src.eqns_.size());
54
55 numRHSs_ = src.numRHSs_;
56
57 for(i=0; i<len; i++) {
58 //eqns_ is a table. Each row of the table needs to be allocated and
59 //copied here. We'll use the fei::CSVec copy constructor to copy the
60 //contents of each existing row into the 'dest' rows.
61
62 //first get a pointer to the row,
63 fei::CSVec* row = src.eqns_[i];
64
65 //now allocate the eqns_ row and the coefs row
66 eqns_[i] = new fei::CSVec(*row);
67
68 //since we allow for multiple rhs's, rhsCoefs_ is a table too...
69 std::vector<double>* rhsCoefs = src.rhsCoefs_[i];
70
71 rhsCoefs_.push_back( new std::vector<double>(*rhsCoefs) );
72 }
73
74 return(*this);
75}
76
77//==============================================================================
79 deleteMemory();
80}
81
82//==============================================================================
84{
85 EqnBuffer* dest = new EqnBuffer;
86
87 *dest = *this;
88
89 return(dest);
90}
91
92//==============================================================================
93void EqnBuffer::deleteMemory() {
94 for(int i=0; i<getNumEqns(); i++) {
95 delete eqns_[i];
96
97 delete rhsCoefs_[i];
98 }
99
100 eqns_.clear();
101 rhsCoefs_.clear();
102 numRHSs_ = 0;
103}
104
105//==============================================================================
107
108 return(fei::binarySearch(eqn, eqnNumbers_));
109}
110
111//==============================================================================
113
114 if (n <= 0) { return;}
115
116 numRHSs_ = n;
117
118 if (getNumEqns() <= 0) return;
119
120 for(int i=0; i<getNumEqns(); i++) {
121 std::vector<double>* rhsCoefs = rhsCoefs_[i];
122 rhsCoefs->assign(numRHSs_, 0.0);
123 }
124}
125
126//==============================================================================
127int EqnBuffer::addRHS(int eqnNumber, int rhsIndex, double value,
128 bool accumulate)
129{
130 int index = fei::binarySearch(eqnNumber, eqnNumbers_);
131
132 if (index < 0) {
133 fei::console_out() << "(deep in FEI) EqnBuffer::addRHS: ERROR, eqnNumber " << eqnNumber
134 << " not found in send eqns." << FEI_ENDL;
135 return(-1);
136 }
137
138 std::vector<double>* rhsCoefs = rhsCoefs_[index];
139
140 if ( (int)rhsCoefs->size() <= rhsIndex) setNumRHSs(rhsIndex+1);
141
142 if (accumulate==true) (*rhsCoefs)[rhsIndex] += value;
143 else (*rhsCoefs)[rhsIndex] = value;
144
145 return(0);
146}
147//==============================================================================
149{
150 //
151 //This function checks the indices_ table to see if 'eqn' is present.
152 //If it is, the appropriate row index into the table is returned.
153 //-1 is return otherwise.
154 //
155 if (indices_union_.size() > 0) {
156 int index = fei::binarySearch(eqn, &indices_union_[0], indices_union_.size());
157 if (index < 0) return(-1);
158 }
159
160 int numEqns = getNumEqns(), index;
161 fei::CSVec** eqnsPtr = &eqns_[0];
162 for(int i=0; i<numEqns; i++) {
163 std::vector<int>& indices = eqnsPtr[i]->indices();
164 index = fei::binarySearch(eqn, &indices[0], indices.size());
165 if (index > -1) return(i);
166 }
167
168 return(-1);
169}
170
171//==============================================================================
172int EqnBuffer::addEqn(int eqnNumber, const double* coefs, const int* indices,
173 int len, bool accumulate, bool create_indices_union)
174{
175 if (len <= 0) return(0);
176
177 int err, insertPoint = -1;
178 int index = fei::binarySearch(eqnNumber, eqnNumbers_, insertPoint);
179
180 if (index < 0) {
181 //if eqnNumber isn't already present, insert a new entry into the
182 //appropriate data structures.
183 err = insertNewEqn(eqnNumber, insertPoint);
184 if (err) {return(err);}
185 index = insertPoint;
186 }
187
188 //Now add the coef/index values.
189 err = internalAddEqn(index, coefs, indices, len, accumulate);
190
191 if (create_indices_union) {
192 for(int i=0; i<len; ++i) {
193 fei::sortedListInsert(indices[i], indices_union_);
194 }
195 }
196
197 return(err);
198}
199
200//==============================================================================
201int EqnBuffer::getCoef(int eqnNumber, int colIndex, double& coef)
202{
203 int eqnLoc = fei::binarySearch(eqnNumber, eqnNumbers_);
204 if (eqnLoc < 0) return(-1);
205
206 int colLoc = fei::binarySearch(colIndex, eqns_[eqnLoc]->indices());
207 if (colLoc < 0) return(-1);
208
209 coef = eqns_[eqnLoc]->coefs()[colLoc];
210 return(0);
211}
212
213//==============================================================================
214int EqnBuffer::removeIndex(int eqnNumber, int colIndex)
215{
216 int eqnLoc = fei::binarySearch(eqnNumber, eqnNumbers_);
217 if (eqnLoc < 0) return(-1);
218
219 int colLoc = fei::binarySearch(colIndex, eqns_[eqnLoc]->indices());
220 if (colLoc < 0) return(0);
221
222 std::vector<int>& indices = eqns_[eqnLoc]->indices();
223 std::vector<double>& coefs= eqns_[eqnLoc]->coefs();
224
225 int len = indices.size();
226
227 int* indPtr = &indices[0];
228 double* coefPtr = &coefs[0];
229
230 for(int i=len-1; i>colLoc; --i) {
231 indPtr[i-1] = indPtr[i];
232 coefPtr[i-1] = coefPtr[i];
233 }
234
235 indices.resize(len-1);
236 coefs.resize(len-1);
237
238 return(0);
239}
240
241//==============================================================================
242int EqnBuffer::getCoefAndRemoveIndex(int eqnNumber, int colIndex, double& coef)
243{
244 int eqnLoc = fei::binarySearch(eqnNumber, eqnNumbers_);
245 if (eqnLoc < 0) return(-1);
246
247 int colLoc = fei::binarySearch(colIndex, eqns_[eqnLoc]->indices());
248 if (colLoc < 0) return(-1);
249
250 std::vector<int>& indices = eqns_[eqnLoc]->indices();
251 std::vector<double>& coefs= eqns_[eqnLoc]->coefs();
252
253 coef = coefs[colLoc];
254 int len = indices.size();
255
256 int* indPtr = &indices[0];
257 double* coefPtr = &coefs[0];
258
259 for(int i=len-1; i>colLoc; --i) {
260 indPtr[i-1] = indPtr[i];
261 coefPtr[i-1] = coefPtr[i];
262 }
263
264 indices.resize(len-1);
265 coefs.resize(len-1);
266
267 return(0);
268}
269
270//==============================================================================
271int EqnBuffer::addEqns(EqnBuffer& inputEqns, bool accumulate)
272{
273 if (inputEqns.eqnNumbers().size() < 1) {
274 return(0);
275 }
276
277 int* eqnNums = &(inputEqns.eqnNumbers()[0]);
278 fei::CSVec** eqs = &(inputEqns.eqns()[0]);
279
280 int numRHSs = inputEqns.getNumRHSs();
281 std::vector<double>** rhsCoefs = &((*(inputEqns.rhsCoefsPtr()))[0]);
282
283 for(int i=0; i<inputEqns.getNumEqns(); i++) {
284 std::vector<int>& indices_i = eqs[i]->indices();
285 std::vector<double>& coefs_i = eqs[i]->coefs();
286
287 int err = addEqn(eqnNums[i], &coefs_i[0], &indices_i[0],
288 eqs[i]->size(), accumulate);
289 if (err) return(err);
290
291 if (numRHSs > 0) {
292 for(int j=0; j<numRHSs; ++j) {
293 addRHS(eqnNums[i], j, (*(rhsCoefs[i]))[j], accumulate);
294 }
295 }
296 }
297
298 return(0);
299}
300
301//==============================================================================
302int EqnBuffer::insertNewEqn(int eqn, int insertPoint)
303{
304 //private function. We can safely assume that insertPoint is the correct
305 //offset at which to insert the new equation.
306 try {
307 eqnNumbers_.insert(eqnNumbers_.begin()+insertPoint, eqn);
308
309 fei::CSVec* newEqn = new fei::CSVec;
310 eqns_.insert(eqns_.begin()+insertPoint, newEqn);
311
312 if (numRHSs_ <= 0) return(-1);
313
314 std::vector<double>* newRhsCoefRow = new std::vector<double>(numRHSs_, 0.0);
315 rhsCoefs_.insert(rhsCoefs_.begin()+insertPoint, newRhsCoefRow);
316 }
317 catch (std::runtime_error& exc) {
318 fei::console_out() << exc.what() << FEI_ENDL;
319 return(-1);
320 }
321
322 return(0);
323}
324
325//==============================================================================
326int EqnBuffer::internalAddEqn(int index, const double* coefs,
327 const int* indices, int len, bool accumulate)
328{
329 //
330 //Private EqnBuffer function. We can safely assume that this function is only
331 //called if indices_ and coefs_ already contain an 'index'th row.
332 //
333
334 fei::CSVec& eqn = *(eqns_[index]);
335
336 if (accumulate) {
337 for(int i=0; i<len; ++i) {
338 fei::add_entry(eqn, indices[i], coefs[i]);
339 }
340 }
341 else {
342 for(int i=0; i<len; ++i) {
343 fei::put_entry(eqn, indices[i], coefs[i]);
344 }
345 }
346
347 return(0);
348}
349
350//==============================================================================
352
353 for(int i=0; i<getNumEqns(); i++) {
354 fei::set_values(*eqns_[i], 0.0);
355 }
356}
357
358//==============================================================================
359int EqnBuffer::addIndices(int eqnNumber, const int* indices, int len)
360{
361 int err = 0, insertPoint = -1;
362 int index = fei::binarySearch(eqnNumber, eqnNumbers_, insertPoint);
363
364 //(we're adding dummy coefs as well, even though there are no
365 //incoming coefs at this point).
366
367 if ((int)dummyCoefs_.size() < len) {
368 dummyCoefs_.assign(len, 0.0);
369 }
370
371 if (index < 0) {
372 //if eqnNumber was not already present, insert new equation
373
374 err = insertNewEqn(eqnNumber, insertPoint);
375 if (err) {return(err);}
376 index = insertPoint;
377 }
378
379 if (len > 0) {
380 err = internalAddEqn(index, &dummyCoefs_[0], indices, len, true);
381 }
382 return(err);
383}
384
385FEI_OSTREAM& operator<<(FEI_OSTREAM& os, EqnBuffer& eq)
386{
387 std::vector<int>& eqnNums = eq.eqnNumbers();
388 std::vector<std::vector<double>*>& rhsCoefs = *(eq.rhsCoefsPtr());
389
390 os << "#ereb num-eqns: " << eqnNums.size() << FEI_ENDL;
391 for(size_t i=0; i<eqnNums.size(); i++) {
392 os << "#ereb eqn " << eqnNums[i] << ": ";
393
394 std::vector<int>& inds = eq.eqns()[i]->indices();
395 std::vector<double>& cfs = eq.eqns()[i]->coefs();
396
397 for(size_t j=0; j<inds.size(); j++) {
398 os << "("<<inds[j] << "," << cfs[j] << ") ";
399 }
400
401 os << " rhs: ";
402 std::vector<double>& rhs = *(rhsCoefs[i]);
403 for(size_t k=0; k<rhs.size(); k++) {
404 os << rhs[k] << ", ";
405 }
406
407 os << FEI_ENDL;
408 }
409
410 return(os);
411}
int getNumEqns()
EqnBuffer * deepCopy()
void resetCoefs()
int isInIndices(int eqn)
void setNumRHSs(int n)
int getCoefAndRemoveIndex(int eqnNumber, int colIndex, double &coef)
int getNumRHSs()
virtual ~EqnBuffer()
std::vector< std::vector< double > * > * rhsCoefsPtr()
int addEqns(EqnBuffer &inputEqns, bool accumulate)
int removeIndex(int eqnNumber, int colIndex)
std::vector< int > & eqnNumbers()
int addEqn(int eqnNumber, const double *coefs, const int *indices, int len, bool accumulate, bool create_indices_union=false)
int addRHS(int eqnNumber, int rhsIndex, double value, bool accumulate=true)
EqnBuffer & operator=(const EqnBuffer &src)
int getCoef(int eqnNumber, int colIndex, double &coef)
int getEqnIndex(int eqn)
std::vector< fei::CSVec * > & eqns()
int addIndices(int eqnNumber, const int *indices, int len)
int binarySearch(const T &item, const T *list, int len)
std::ostream & console_out()
int sortedListInsert(const T &item, std::vector< T > &list)