FEI  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
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 //==============================================================================
93 void 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 //==============================================================================
127 int 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 //==============================================================================
172 int 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 //==============================================================================
201 int 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 //==============================================================================
214 int 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 //==============================================================================
242 int 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 //==============================================================================
271 int 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 //==============================================================================
302 int 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 //==============================================================================
326 int 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 //==============================================================================
359 int 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 
385 FEI_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 isInIndices(int eqn)
int sortedListInsert(const T &item, std::vector< T > &list)
std::ostream & operator<<(std::ostream &os, const Teuchos::SerialBandDenseMatrix< OrdinalType, ScalarType > &obj)
int addEqn(int eqnNumber, const double *coefs, const int *indices, int len, bool accumulate, bool create_indices_union=false)
std::vector< fei::CSVec * > & eqns()
EqnBuffer & operator=(const EqnBuffer &src)
void setNumRHSs(int n)
std::vector< int > & eqnNumbers()
int addRHS(int eqnNumber, int rhsIndex, double value, bool accumulate=true)
std::vector< std::vector< double > * > * rhsCoefsPtr()
int binarySearch(const T &item, const T *list, int len)
int getNumRHSs()
EqnBuffer * deepCopy()
int getEqnIndex(int eqn)
int getCoefAndRemoveIndex(int eqnNumber, int colIndex, double &coef)
int addEqns(EqnBuffer &inputEqns, bool accumulate)
int addIndices(int eqnNumber, const int *indices, int len)
std::ostream & console_out()
virtual ~EqnBuffer()
int removeIndex(int eqnNumber, int colIndex)
int getCoef(int eqnNumber, int colIndex, double &coef)
void resetCoefs()
int getNumEqns()