FEI Version of the Day
fei_EqnBuffer.cpp
00001 /*
00002 // @HEADER
00003 // ************************************************************************
00004 //             FEI: Finite Element Interface to Linear Solvers
00005 //                  Copyright (2005) Sandia Corporation.
00006 //
00007 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation, the
00008 // U.S. Government retains certain rights in this software.
00009 //
00010 // Redistribution and use in source and binary forms, with or without
00011 // modification, are permitted provided that the following conditions are
00012 // met:
00013 //
00014 // 1. Redistributions of source code must retain the above copyright
00015 // notice, this list of conditions and the following disclaimer.
00016 //
00017 // 2. Redistributions in binary form must reproduce the above copyright
00018 // notice, this list of conditions and the following disclaimer in the
00019 // documentation and/or other materials provided with the distribution.
00020 //
00021 // 3. Neither the name of the Corporation nor the names of the
00022 // contributors may be used to endorse or promote products derived from
00023 // this software without specific prior written permission.
00024 //
00025 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00026 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00027 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00028 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00029 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00030 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00031 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00032 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00033 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00034 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00035 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00036 //
00037 // Questions? Contact Alan Williams (william@sandia.gov) 
00038 //
00039 // ************************************************************************
00040 // @HEADER
00041 */
00042 
00043 
00044 #include <fei_macros.hpp>
00045 
00046 #include <fei_EqnBuffer.hpp>
00047 #include <fei_CSVec.hpp>
00048 
00049 #include <fei_TemplateUtils.hpp>
00050 
00051 //==============================================================================
00052 EqnBuffer::EqnBuffer()
00053  : newCoefData_(0),
00054    newRHSData_(0),
00055    eqnNumbers_(0),
00056    eqns_(),
00057    indices_union_(0),
00058    numRHSs_(1),
00059    rhsCoefs_(),
00060    setNumRHSsCalled_(false),
00061    rhsCoefsAllocated_(false),
00062    dummyCoefs_()
00063 {
00064 }
00065 
00066 //==============================================================================
00067 EqnBuffer::EqnBuffer(const EqnBuffer& src)
00068  : newCoefData_(0),
00069    newRHSData_(0),
00070    eqnNumbers_(0),
00071    eqns_(),
00072    indices_union_(0),
00073    numRHSs_(1),
00074    rhsCoefs_(),
00075    setNumRHSsCalled_(false),
00076    rhsCoefsAllocated_(false),
00077    dummyCoefs_()
00078 {
00079   *this = src;
00080 }
00081 
00082 //==============================================================================
00083 EqnBuffer& EqnBuffer::operator=(const EqnBuffer& src)
00084 {
00085    int i, len = src.eqnNumbers_.size();
00086 
00087    eqnNumbers_ = src.eqnNumbers_;
00088    eqns_.resize(src.eqns_.size());
00089 
00090    numRHSs_ = src.numRHSs_;
00091 
00092    for(i=0; i<len; i++) {
00093      //eqns_ is a table. Each row of the table needs to be allocated and
00094      //copied here. We'll use the fei::CSVec copy constructor to copy the
00095      //contents of each existing row into the 'dest' rows.
00096 
00097       //first get a pointer to the row,
00098       fei::CSVec* row = src.eqns_[i];
00099 
00100       //now allocate the eqns_ row and the coefs row
00101       eqns_[i] = new fei::CSVec(*row);
00102 
00103       //since we allow for multiple rhs's, rhsCoefs_ is a table too...
00104       std::vector<double>* rhsCoefs = src.rhsCoefs_[i];
00105 
00106       rhsCoefs_.push_back( new std::vector<double>(*rhsCoefs) );
00107    }
00108 
00109    return(*this);
00110 }
00111 
00112 //==============================================================================
00113 EqnBuffer::~EqnBuffer() {
00114    deleteMemory();
00115 }
00116 
00117 //==============================================================================
00118 EqnBuffer* EqnBuffer::deepCopy()
00119 {
00120    EqnBuffer* dest = new EqnBuffer;
00121 
00122    *dest = *this;
00123  
00124    return(dest);
00125 }
00126 
00127 //==============================================================================
00128 void EqnBuffer::deleteMemory() {
00129    for(int i=0; i<getNumEqns(); i++) {
00130       delete eqns_[i];
00131 
00132       delete rhsCoefs_[i];
00133    }
00134 
00135    eqns_.clear();
00136    rhsCoefs_.clear();
00137    numRHSs_ = 0;
00138 }
00139 
00140 //==============================================================================
00141 int EqnBuffer::getEqnIndex(int eqn) {
00142 
00143    return(fei::binarySearch(eqn, eqnNumbers_));
00144 }
00145 
00146 //==============================================================================
00147 void EqnBuffer::setNumRHSs(int n) {
00148 
00149    if (n <= 0) { return;}
00150 
00151    numRHSs_ = n;
00152 
00153    if (getNumEqns() <= 0) return;
00154 
00155    for(int i=0; i<getNumEqns(); i++) {
00156      std::vector<double>* rhsCoefs = rhsCoefs_[i];
00157      rhsCoefs->assign(numRHSs_, 0.0);
00158    }
00159 }
00160 
00161 //==============================================================================
00162 int EqnBuffer::addRHS(int eqnNumber, int rhsIndex, double value,
00163           bool accumulate)
00164 {
00165    int index = fei::binarySearch(eqnNumber, eqnNumbers_);
00166 
00167    if (index < 0) {
00168       fei::console_out() << "(deep in FEI) EqnBuffer::addRHS: ERROR, eqnNumber " << eqnNumber
00169            << " not found in send eqns." << FEI_ENDL;
00170       return(-1);
00171    }
00172 
00173    std::vector<double>* rhsCoefs = rhsCoefs_[index];
00174 
00175    if ( (int)rhsCoefs->size() <= rhsIndex) setNumRHSs(rhsIndex+1);
00176 
00177    if (accumulate==true) (*rhsCoefs)[rhsIndex] += value;
00178    else (*rhsCoefs)[rhsIndex] = value;
00179 
00180    return(0);
00181 }
00182 //==============================================================================
00183 int EqnBuffer::isInIndices(int eqn)
00184 {
00185   //
00186   //This function checks the indices_ table to see if 'eqn' is present.
00187   //If it is, the appropriate row index into the table is returned.
00188   //-1 is return otherwise.
00189   //
00190   if (indices_union_.size() > 0) {
00191     int index = fei::binarySearch(eqn, &indices_union_[0], indices_union_.size());
00192     if (index < 0) return(-1);
00193   }
00194 
00195   int numEqns = getNumEqns(), index;
00196   fei::CSVec** eqnsPtr = &eqns_[0];
00197   for(int i=0; i<numEqns; i++) {
00198     std::vector<int>& indices = eqnsPtr[i]->indices();
00199     index = fei::binarySearch(eqn, &indices[0], indices.size());
00200     if (index > -1) return(i);
00201   }
00202 
00203   return(-1);
00204 }
00205 
00206 //==============================================================================
00207 int EqnBuffer::addEqn(int eqnNumber, const double* coefs, const int* indices,
00208                        int len, bool accumulate, bool create_indices_union) 
00209 {
00210   if (len <= 0) return(0);
00211 
00212   int err, insertPoint = -1;
00213   int index = fei::binarySearch(eqnNumber, eqnNumbers_, insertPoint);
00214 
00215   if (index < 0) {
00216     //if eqnNumber isn't already present, insert a new entry into the
00217     //appropriate data structures.
00218     err = insertNewEqn(eqnNumber, insertPoint);
00219     if (err) {return(err);}
00220     index = insertPoint;
00221   }
00222 
00223   //Now add the coef/index values.
00224   err = internalAddEqn(index, coefs, indices, len, accumulate);
00225 
00226   if (create_indices_union) {
00227     for(int i=0; i<len; ++i) {
00228       fei::sortedListInsert(indices[i], indices_union_);
00229     }
00230   }
00231 
00232   return(err);
00233 }
00234 
00235 //==============================================================================
00236 int EqnBuffer::getCoef(int eqnNumber, int colIndex, double& coef)
00237 {
00238   int eqnLoc = fei::binarySearch(eqnNumber, eqnNumbers_);
00239   if (eqnLoc < 0) return(-1);
00240 
00241   int colLoc = fei::binarySearch(colIndex, eqns_[eqnLoc]->indices());
00242   if (colLoc < 0) return(-1);
00243 
00244   coef = eqns_[eqnLoc]->coefs()[colLoc];
00245   return(0);
00246 }
00247 
00248 //==============================================================================
00249 int EqnBuffer::removeIndex(int eqnNumber, int colIndex)
00250 {
00251   int eqnLoc = fei::binarySearch(eqnNumber, eqnNumbers_);
00252   if (eqnLoc < 0) return(-1);
00253 
00254   int colLoc = fei::binarySearch(colIndex, eqns_[eqnLoc]->indices());
00255   if (colLoc < 0) return(0);
00256 
00257   std::vector<int>& indices = eqns_[eqnLoc]->indices();
00258   std::vector<double>& coefs= eqns_[eqnLoc]->coefs();
00259 
00260   int len = indices.size();
00261 
00262   int* indPtr = &indices[0];
00263   double* coefPtr = &coefs[0];
00264 
00265   for(int i=len-1; i>colLoc; --i) {
00266     indPtr[i-1] = indPtr[i];
00267     coefPtr[i-1] = coefPtr[i];
00268   }
00269 
00270   indices.resize(len-1);
00271   coefs.resize(len-1);
00272 
00273   return(0);
00274 }
00275 
00276 //==============================================================================
00277 int EqnBuffer::getCoefAndRemoveIndex(int eqnNumber, int colIndex, double& coef)
00278 {
00279   int eqnLoc = fei::binarySearch(eqnNumber, eqnNumbers_);
00280   if (eqnLoc < 0) return(-1);
00281 
00282   int colLoc = fei::binarySearch(colIndex, eqns_[eqnLoc]->indices());
00283   if (colLoc < 0) return(-1);
00284 
00285   std::vector<int>& indices = eqns_[eqnLoc]->indices();
00286   std::vector<double>& coefs= eqns_[eqnLoc]->coefs();
00287 
00288   coef = coefs[colLoc];
00289   int len = indices.size();
00290 
00291   int* indPtr = &indices[0];
00292   double* coefPtr = &coefs[0];
00293 
00294   for(int i=len-1; i>colLoc; --i) {
00295     indPtr[i-1] = indPtr[i];
00296     coefPtr[i-1] = coefPtr[i];
00297   }
00298 
00299   indices.resize(len-1);
00300   coefs.resize(len-1);
00301 
00302   return(0);
00303 }
00304 
00305 //==============================================================================
00306 int EqnBuffer::addEqns(EqnBuffer& inputEqns, bool accumulate)
00307 {
00308   if (inputEqns.eqnNumbers().size() < 1) {
00309     return(0);
00310   }
00311 
00312   int* eqnNums = &(inputEqns.eqnNumbers()[0]);
00313   fei::CSVec** eqs = &(inputEqns.eqns()[0]);
00314 
00315   int numRHSs = inputEqns.getNumRHSs();
00316   std::vector<double>** rhsCoefs = &((*(inputEqns.rhsCoefsPtr()))[0]);
00317 
00318   for(int i=0; i<inputEqns.getNumEqns(); i++) {
00319     std::vector<int>& indices_i  = eqs[i]->indices();
00320     std::vector<double>& coefs_i = eqs[i]->coefs();
00321 
00322     int err = addEqn(eqnNums[i], &coefs_i[0], &indices_i[0],
00323         eqs[i]->size(), accumulate);
00324     if (err) return(err);
00325 
00326     if (numRHSs > 0) {
00327       for(int j=0; j<numRHSs; ++j) {
00328   addRHS(eqnNums[i], j, (*(rhsCoefs[i]))[j], accumulate);
00329       }
00330     }
00331   }
00332 
00333   return(0);
00334 }
00335 
00336 //==============================================================================
00337 int EqnBuffer::insertNewEqn(int eqn, int insertPoint)
00338 {
00339   //private function. We can safely assume that insertPoint is the correct
00340   //offset at which to insert the new equation.
00341   try {
00342     eqnNumbers_.insert(eqnNumbers_.begin()+insertPoint, eqn);
00343 
00344     fei::CSVec* newEqn = new fei::CSVec;
00345     eqns_.insert(eqns_.begin()+insertPoint, newEqn);
00346 
00347     if (numRHSs_ <= 0) return(-1);
00348 
00349     std::vector<double>* newRhsCoefRow = new std::vector<double>(numRHSs_, 0.0);
00350     rhsCoefs_.insert(rhsCoefs_.begin()+insertPoint, newRhsCoefRow);
00351   }
00352   catch (std::runtime_error& exc) {
00353     fei::console_out() << exc.what() << FEI_ENDL;
00354     return(-1);
00355   }
00356 
00357   return(0);
00358 }
00359 
00360 //==============================================================================
00361 int EqnBuffer::internalAddEqn(int index, const double* coefs,
00362             const int* indices, int len, bool accumulate) 
00363 {
00364   //
00365   //Private EqnBuffer function. We can safely assume that this function is only
00366   //called if indices_ and coefs_ already contain an 'index'th row.
00367   //
00368 
00369   fei::CSVec& eqn = *(eqns_[index]);
00370 
00371   if (accumulate) {
00372     for(int i=0; i<len; ++i) {
00373       fei::add_entry(eqn, indices[i], coefs[i]);
00374     }
00375   }
00376   else {
00377     for(int i=0; i<len; ++i) {
00378       fei::put_entry(eqn, indices[i], coefs[i]);
00379     }
00380   }
00381 
00382   return(0);
00383 }
00384 
00385 //==============================================================================
00386 void EqnBuffer::resetCoefs() {
00387 
00388   for(int i=0; i<getNumEqns(); i++) {
00389     fei::set_values(*eqns_[i], 0.0);
00390   }
00391 }
00392 
00393 //==============================================================================
00394 int EqnBuffer::addIndices(int eqnNumber, const int* indices, int len)
00395 {
00396    int err = 0, insertPoint = -1;
00397    int index = fei::binarySearch(eqnNumber, eqnNumbers_, insertPoint);
00398 
00399    //(we're adding dummy coefs as well, even though there are no
00400    //incoming coefs at this point).
00401 
00402    if ((int)dummyCoefs_.size() < len) {
00403      dummyCoefs_.assign(len, 0.0);
00404    }
00405 
00406    if (index < 0) {
00407      //if eqnNumber was not already present, insert new equation
00408 
00409      err = insertNewEqn(eqnNumber, insertPoint);
00410      if (err) {return(err);}
00411      index = insertPoint;
00412    }
00413 
00414    if (len > 0) {
00415      err = internalAddEqn(index, &dummyCoefs_[0], indices, len, true);
00416    }
00417    return(err);
00418 }
00419 
00420 FEI_OSTREAM& operator<<(FEI_OSTREAM& os, EqnBuffer& eq)
00421 {
00422   std::vector<int>& eqnNums = eq.eqnNumbers();
00423   std::vector<std::vector<double>*>& rhsCoefs = *(eq.rhsCoefsPtr());
00424 
00425   os << "#ereb num-eqns: " << eqnNums.size() << FEI_ENDL;
00426   for(size_t i=0; i<eqnNums.size(); i++) {
00427     os << "#ereb eqn " << eqnNums[i] << ": ";
00428 
00429     std::vector<int>& inds = eq.eqns()[i]->indices();
00430     std::vector<double>& cfs = eq.eqns()[i]->coefs();
00431 
00432     for(size_t j=0; j<inds.size(); j++) {
00433       os << "("<<inds[j] << "," << cfs[j] << ") ";
00434     }
00435 
00436     os << " rhs: ";
00437     std::vector<double>& rhs = *(rhsCoefs[i]);
00438     for(size_t k=0; k<rhs.size(); k++) {
00439       os << rhs[k] << ", ";
00440     }
00441 
00442     os << FEI_ENDL;
00443   }
00444 
00445   return(os);
00446 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Friends