fei_EqnBuffer.cpp

00001 /*--------------------------------------------------------------------*/
00002 /*    Copyright 2005 Sandia Corporation.                              */
00003 /*    Under the terms of Contract DE-AC04-94AL85000, there is a       */
00004 /*    non-exclusive license for use of this work by or on behalf      */
00005 /*    of the U.S. Government.  Export of this program may require     */
00006 /*    a license from the United States Government.                    */
00007 /*--------------------------------------------------------------------*/
00008 
00009 #include <fei_macros.hpp>
00010 
00011 #include <fei_EqnBuffer.hpp>
00012 #include <fei_CSVec.hpp>
00013 
00014 #include <fei_TemplateUtils.hpp>
00015 
00016 //==============================================================================
00017 EqnBuffer::EqnBuffer()
00018  : newCoefData_(0),
00019    newRHSData_(0),
00020    eqnNumbers_(0),
00021    eqns_(),
00022    indices_union_(0),
00023    numRHSs_(1),
00024    rhsCoefs_(),
00025    setNumRHSsCalled_(false),
00026    rhsCoefsAllocated_(false),
00027    dummyCoefs_()
00028 {
00029 }
00030 
00031 //==============================================================================
00032 EqnBuffer::EqnBuffer(const EqnBuffer& src)
00033  : newCoefData_(0),
00034    newRHSData_(0),
00035    eqnNumbers_(0),
00036    eqns_(),
00037    indices_union_(0),
00038    numRHSs_(1),
00039    rhsCoefs_(),
00040    setNumRHSsCalled_(false),
00041    rhsCoefsAllocated_(false),
00042    dummyCoefs_()
00043 {
00044   *this = src;
00045 }
00046 
00047 //==============================================================================
00048 EqnBuffer& EqnBuffer::operator=(const EqnBuffer& src)
00049 {
00050    int i, len = src.eqnNumbers_.size();
00051 
00052    eqnNumbers_ = src.eqnNumbers_;
00053    eqns_.resize(src.eqns_.size());
00054 
00055    numRHSs_ = src.numRHSs_;
00056 
00057    for(i=0; i<len; i++) {
00058      //eqns_ is a table. Each row of the table needs to be allocated and
00059      //copied here. We'll use the fei::CSVec copy constructor to copy the
00060      //contents of each existing row into the 'dest' rows.
00061 
00062       //first get a pointer to the row,
00063       fei::CSVec* row = src.eqns_[i];
00064 
00065       //now allocate the eqns_ row and the coefs row
00066       eqns_[i] = new fei::CSVec(*row);
00067 
00068       //since we allow for multiple rhs's, rhsCoefs_ is a table too...
00069       std::vector<double>* rhsCoefs = src.rhsCoefs_[i];
00070 
00071       rhsCoefs_.push_back( new std::vector<double>(*rhsCoefs) );
00072    }
00073 
00074    return(*this);
00075 }
00076 
00077 //==============================================================================
00078 EqnBuffer::~EqnBuffer() {
00079    deleteMemory();
00080 }
00081 
00082 //==============================================================================
00083 EqnBuffer* EqnBuffer::deepCopy()
00084 {
00085    EqnBuffer* dest = new EqnBuffer;
00086 
00087    *dest = *this;
00088  
00089    return(dest);
00090 }
00091 
00092 //==============================================================================
00093 void EqnBuffer::deleteMemory() {
00094    for(int i=0; i<getNumEqns(); i++) {
00095       delete eqns_[i];
00096 
00097       delete rhsCoefs_[i];
00098    }
00099 
00100    eqns_.clear();
00101    rhsCoefs_.clear();
00102    numRHSs_ = 0;
00103 }
00104 
00105 //==============================================================================
00106 int EqnBuffer::getEqnIndex(int eqn) {
00107 
00108    return(fei::binarySearch(eqn, eqnNumbers_));
00109 }
00110 
00111 //==============================================================================
00112 void EqnBuffer::setNumRHSs(int n) {
00113 
00114    if (n <= 0) { return;}
00115 
00116    numRHSs_ = n;
00117 
00118    if (getNumEqns() <= 0) return;
00119 
00120    for(int i=0; i<getNumEqns(); i++) {
00121      std::vector<double>* rhsCoefs = rhsCoefs_[i];
00122      rhsCoefs->assign(numRHSs_, 0.0);
00123    }
00124 }
00125 
00126 //==============================================================================
00127 int EqnBuffer::addRHS(int eqnNumber, int rhsIndex, double value,
00128           bool accumulate)
00129 {
00130    int index = fei::binarySearch(eqnNumber, eqnNumbers_);
00131 
00132    if (index < 0) {
00133       FEI_CERR << "(deep in FEI) EqnBuffer::addRHS: ERROR, eqnNumber " << eqnNumber
00134            << " not found in send eqns." << FEI_ENDL;
00135       return(-1);
00136    }
00137 
00138    std::vector<double>* rhsCoefs = rhsCoefs_[index];
00139 
00140    if ( (int)rhsCoefs->size() <= rhsIndex) setNumRHSs(rhsIndex+1);
00141 
00142    if (accumulate==true) (*rhsCoefs)[rhsIndex] += value;
00143    else (*rhsCoefs)[rhsIndex] = value;
00144 
00145    return(0);
00146 }
00147 //==============================================================================
00148 int EqnBuffer::isInIndices(int eqn)
00149 {
00150   //
00151   //This function checks the indices_ table to see if 'eqn' is present.
00152   //If it is, the appropriate row index into the table is returned.
00153   //-1 is return otherwise.
00154   //
00155   if (indices_union_.size() > 0) {
00156     int index = fei::binarySearch(eqn, &indices_union_[0], indices_union_.size());
00157     if (index < 0) return(-1);
00158   }
00159 
00160   int numEqns = getNumEqns(), index;
00161   fei::CSVec** eqnsPtr = &eqns_[0];
00162   for(int i=0; i<numEqns; i++) {
00163     std::vector<int>& indices = eqnsPtr[i]->indices();
00164     index = fei::binarySearch(eqn, &indices[0], indices.size());
00165     if (index > -1) return(i);
00166   }
00167 
00168   return(-1);
00169 }
00170 
00171 //==============================================================================
00172 int EqnBuffer::addEqn(int eqnNumber, const double* coefs, const int* indices,
00173                        int len, bool accumulate, bool create_indices_union) 
00174 {
00175   if (len <= 0) return(0);
00176 
00177   int err, insertPoint = -1;
00178   int index = fei::binarySearch(eqnNumber, eqnNumbers_, insertPoint);
00179 
00180   if (index < 0) {
00181     //if eqnNumber isn't already present, insert a new entry into the
00182     //appropriate data structures.
00183     err = insertNewEqn(eqnNumber, insertPoint);
00184     if (err) {return(err);}
00185     index = insertPoint;
00186   }
00187 
00188   //Now add the coef/index values.
00189   err = internalAddEqn(index, coefs, indices, len, accumulate);
00190 
00191   if (create_indices_union) {
00192     for(int i=0; i<len; ++i) {
00193       fei::sortedListInsert(indices[i], indices_union_);
00194     }
00195   }
00196 
00197   return(err);
00198 }
00199 
00200 //==============================================================================
00201 int EqnBuffer::getCoef(int eqnNumber, int colIndex, double& coef)
00202 {
00203   int eqnLoc = fei::binarySearch(eqnNumber, eqnNumbers_);
00204   if (eqnLoc < 0) return(-1);
00205 
00206   int colLoc = fei::binarySearch(colIndex, eqns_[eqnLoc]->indices());
00207   if (colLoc < 0) return(-1);
00208 
00209   coef = eqns_[eqnLoc]->coefs()[colLoc];
00210   return(0);
00211 }
00212 
00213 //==============================================================================
00214 int EqnBuffer::removeIndex(int eqnNumber, int colIndex)
00215 {
00216   int eqnLoc = fei::binarySearch(eqnNumber, eqnNumbers_);
00217   if (eqnLoc < 0) return(-1);
00218 
00219   int colLoc = fei::binarySearch(colIndex, eqns_[eqnLoc]->indices());
00220   if (colLoc < 0) return(0);
00221 
00222   std::vector<int>& indices = eqns_[eqnLoc]->indices();
00223   std::vector<double>& coefs= eqns_[eqnLoc]->coefs();
00224 
00225   int len = indices.size();
00226 
00227   int* indPtr = &indices[0];
00228   double* coefPtr = &coefs[0];
00229 
00230   for(int i=len-1; i>colLoc; --i) {
00231     indPtr[i-1] = indPtr[i];
00232     coefPtr[i-1] = coefPtr[i];
00233   }
00234 
00235   indices.resize(len-1);
00236   coefs.resize(len-1);
00237 
00238   return(0);
00239 }
00240 
00241 //==============================================================================
00242 int EqnBuffer::getCoefAndRemoveIndex(int eqnNumber, int colIndex, double& coef)
00243 {
00244   int eqnLoc = fei::binarySearch(eqnNumber, eqnNumbers_);
00245   if (eqnLoc < 0) return(-1);
00246 
00247   int colLoc = fei::binarySearch(colIndex, eqns_[eqnLoc]->indices());
00248   if (colLoc < 0) return(-1);
00249 
00250   std::vector<int>& indices = eqns_[eqnLoc]->indices();
00251   std::vector<double>& coefs= eqns_[eqnLoc]->coefs();
00252 
00253   coef = coefs[colLoc];
00254   int len = indices.size();
00255 
00256   int* indPtr = &indices[0];
00257   double* coefPtr = &coefs[0];
00258 
00259   for(int i=len-1; i>colLoc; --i) {
00260     indPtr[i-1] = indPtr[i];
00261     coefPtr[i-1] = coefPtr[i];
00262   }
00263 
00264   indices.resize(len-1);
00265   coefs.resize(len-1);
00266 
00267   return(0);
00268 }
00269 
00270 //==============================================================================
00271 int EqnBuffer::addEqns(EqnBuffer& inputEqns, bool accumulate)
00272 {
00273   if (inputEqns.eqnNumbers().size() < 1) {
00274     return(0);
00275   }
00276 
00277   int* eqnNums = &(inputEqns.eqnNumbers()[0]);
00278   fei::CSVec** eqs = &(inputEqns.eqns()[0]);
00279 
00280   int numRHSs = inputEqns.getNumRHSs();
00281   std::vector<double>** rhsCoefs = &((*(inputEqns.rhsCoefsPtr()))[0]);
00282 
00283   for(int i=0; i<inputEqns.getNumEqns(); i++) {
00284     std::vector<int>& indices_i  = eqs[i]->indices();
00285     std::vector<double>& coefs_i = eqs[i]->coefs();
00286 
00287     int err = addEqn(eqnNums[i], &coefs_i[0], &indices_i[0],
00288         eqs[i]->size(), accumulate);
00289     if (err) return(err);
00290 
00291     if (numRHSs > 0) {
00292       for(int j=0; j<numRHSs; ++j) {
00293   addRHS(eqnNums[i], j, (*(rhsCoefs[i]))[j], accumulate);
00294       }
00295     }
00296   }
00297 
00298   return(0);
00299 }
00300 
00301 //==============================================================================
00302 int EqnBuffer::insertNewEqn(int eqn, int insertPoint)
00303 {
00304   //private function. We can safely assume that insertPoint is the correct
00305   //offset at which to insert the new equation.
00306   try {
00307     eqnNumbers_.insert(eqnNumbers_.begin()+insertPoint, eqn);
00308 
00309     fei::CSVec* newEqn = new fei::CSVec;
00310     eqns_.insert(eqns_.begin()+insertPoint, newEqn);
00311 
00312     if (numRHSs_ <= 0) return(-1);
00313 
00314     std::vector<double>* newRhsCoefRow = new std::vector<double>(numRHSs_, 0.0);
00315     rhsCoefs_.insert(rhsCoefs_.begin()+insertPoint, newRhsCoefRow);
00316   }
00317   catch (std::runtime_error& exc) {
00318     FEI_CERR << exc.what() << FEI_ENDL;
00319     return(-1);
00320   }
00321 
00322   return(0);
00323 }
00324 
00325 //==============================================================================
00326 int EqnBuffer::internalAddEqn(int index, const double* coefs,
00327             const int* indices, int len, bool accumulate) 
00328 {
00329   //
00330   //Private EqnBuffer function. We can safely assume that this function is only
00331   //called if indices_ and coefs_ already contain an 'index'th row.
00332   //
00333 
00334   fei::CSVec& eqn = *(eqns_[index]);
00335 
00336   if (accumulate) {
00337     for(int i=0; i<len; ++i) {
00338       fei::add_entry(eqn, indices[i], coefs[i]);
00339     }
00340   }
00341   else {
00342     for(int i=0; i<len; ++i) {
00343       fei::put_entry(eqn, indices[i], coefs[i]);
00344     }
00345   }
00346 
00347   return(0);
00348 }
00349 
00350 //==============================================================================
00351 void EqnBuffer::resetCoefs() {
00352 
00353   for(int i=0; i<getNumEqns(); i++) {
00354     fei::set_values(*eqns_[i], 0.0);
00355   }
00356 }
00357 
00358 //==============================================================================
00359 int EqnBuffer::addIndices(int eqnNumber, const int* indices, int len)
00360 {
00361    int err = 0, insertPoint = -1;
00362    int index = fei::binarySearch(eqnNumber, eqnNumbers_, insertPoint);
00363 
00364    //(we're adding dummy coefs as well, even though there are no
00365    //incoming coefs at this point).
00366 
00367    if ((int)dummyCoefs_.size() < len) {
00368      dummyCoefs_.assign(len, 0.0);
00369    }
00370 
00371    if (index < 0) {
00372      //if eqnNumber was not already present, insert new equation
00373 
00374      err = insertNewEqn(eqnNumber, insertPoint);
00375      if (err) {return(err);}
00376      index = insertPoint;
00377    }
00378 
00379    if (len > 0) {
00380      err = internalAddEqn(index, &dummyCoefs_[0], indices, len, true);
00381    }
00382    return(err);
00383 }
00384 
00385 FEI_OSTREAM& operator<<(FEI_OSTREAM& os, EqnBuffer& eq)
00386 {
00387   std::vector<int>& eqnNums = eq.eqnNumbers();
00388   std::vector<std::vector<double>*>& rhsCoefs = *(eq.rhsCoefsPtr());
00389 
00390   os << "#ereb num-eqns: " << eqnNums.size() << FEI_ENDL;
00391   for(size_t i=0; i<eqnNums.size(); i++) {
00392     os << "#ereb eqn " << eqnNums[i] << ": ";
00393 
00394     std::vector<int>& inds = eq.eqns()[i]->indices();
00395     std::vector<double>& cfs = eq.eqns()[i]->coefs();
00396 
00397     for(size_t j=0; j<inds.size(); j++) {
00398       os << "("<<inds[j] << "," << cfs[j] << ") ";
00399     }
00400 
00401     os << " rhs: ";
00402     std::vector<double>& rhs = *(rhsCoefs[i]);
00403     for(size_t k=0; k<rhs.size(); k++) {
00404       os << rhs[k] << ", ";
00405     }
00406 
00407     os << FEI_ENDL;
00408   }
00409 
00410   return(os);
00411 }

Generated on Tue Jul 13 09:27:44 2010 for FEI by  doxygen 1.4.7