Sacado Package Browser (Single Doxygen Collection) Version of the Day
Sacado_Tay_CacheTaylorImp.hpp
Go to the documentation of this file.
00001 // $Id$ 
00002 // $Source$ 
00003 // @HEADER
00004 // ***********************************************************************
00005 // 
00006 //                           Sacado Package
00007 //                 Copyright (2006) Sandia Corporation
00008 // 
00009 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
00010 // the U.S. Government retains certain rights in this software.
00011 // 
00012 // This library is free software; you can redistribute it and/or modify
00013 // it under the terms of the GNU Lesser General Public License as
00014 // published by the Free Software Foundation; either version 2.1 of the
00015 // License, or (at your option) any later version.
00016 //  
00017 // This library is distributed in the hope that it will be useful, but
00018 // WITHOUT ANY WARRANTY; without even the implied warranty of
00019 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00020 // Lesser General Public License for more details.
00021 //  
00022 // You should have received a copy of the GNU Lesser General Public
00023 // License along with this library; if not, write to the Free Software
00024 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00025 // USA
00026 // Questions? Contact David M. Gay (dmgay@sandia.gov) or Eric T. Phipps
00027 // (etphipp@sandia.gov).
00028 // 
00029 // ***********************************************************************
00030 // @HEADER
00031 
00032 template <typename T> 
00033 template <typename S> 
00034 inline Sacado::Tay::CacheTaylor<T>::CacheTaylor(const Expr<S>& x) : 
00035   Expr< CacheTaylorImplementation<T> >(x.degree(), T(0.))
00036 {
00037   unsigned int d = this->degree();
00038 
00039   x.allocateCache(d);
00040 
00041   // We copy the coefficients from the highest degree to the lowest just
00042   // to be consistent with operator=(), even though it is not strictly 
00043   // necessary since "this" cannot be on the RHS in the copy constructor.
00044   if (x.hasFastAccess(d))
00045     for(int i=d; i>=0; --i) 
00046       this->coeff_[i] = x.fastAccessCoeff(i);
00047   else
00048     for(int i=d; i>=0; --i) 
00049       this->coeff_[i] = x.coeff(i);
00050 
00051 }
00052 
00053 template <typename T> 
00054 inline Sacado::Tay::CacheTaylor<T>& 
00055 Sacado::Tay::CacheTaylor<T>::operator=(const T& v) 
00056 {
00057   this->coeff_[0] = v;
00058 
00059   for (unsigned int i=1; i<this->coeff_.size(); i++)
00060     this->coeff_[i] = T(0.);
00061 
00062   return *this;
00063 }
00064 
00065 template <typename T> 
00066 inline Sacado::Tay::CacheTaylor<T>& 
00067 Sacado::Tay::CacheTaylor<T>::operator=(const Sacado::Tay::CacheTaylor<T>& x) 
00068 {
00069   if (x.coeff_.size() != this->coeff_.size()) 
00070     this->coeff_.resize(x.coeff_.size());
00071   this->coeff_ = x.coeff_;
00072   
00073   return *this;
00074 }
00075 
00076 template <typename T> 
00077 template <typename S> 
00078 inline Sacado::Tay::CacheTaylor<T>& 
00079 Sacado::Tay::CacheTaylor<T>::operator=(const Expr<S>& x) 
00080 {
00081   unsigned int d = this->degree();
00082   unsigned int xd = x.degree();
00083 
00084   // Resize polynomial for "this" if x has greater degree
00085   if (xd > d) {
00086     this->coeff_.resize(xd+1);
00087     d = xd;
00088   }
00089 
00090   x.allocateCache(d);
00091 
00092   // Copy coefficients.  Note:  we copy from the last term down to the first
00093   // to take into account "this" being involved in an expression on the RHS.
00094   // By overwriting the degree k term, we are guarranteed not to affect any
00095   // of the lower degree terms.  However, this in general will force each
00096   // term in the expression to compute all of its coefficients at once instead
00097   // traversing the expression once for each degree.
00098   if (x.hasFastAccess(d))
00099     for(int i=d; i>=0; --i) 
00100       this->coeff_[i] = x.fastAccessCoeff(i);
00101     else
00102       for(int i=d; i>=0; --i) 
00103   this->coeff_[i] = x.coeff(i);
00104   
00105   return *this;
00106 }
00107 
00108 template <typename T> 
00109 inline  Sacado::Tay::CacheTaylor<T>& 
00110 Sacado::Tay::CacheTaylor<T>::operator += (const T& v)
00111 {
00112   this->coeff_[0] += v;
00113 
00114   return *this;
00115 }
00116 
00117 template <typename T> 
00118 inline Sacado::Tay::CacheTaylor<T>& 
00119 Sacado::Tay::CacheTaylor<T>::operator -= (const T& v)
00120 {
00121   this->coeff_[0] -= v;
00122 
00123   return *this;
00124 }
00125 
00126 template <typename T> 
00127 inline Sacado::Tay::CacheTaylor<T>& 
00128 Sacado::Tay::CacheTaylor<T>::operator *= (const T& v)
00129 {
00130   this->coeff_ *= v;
00131 
00132   return *this;
00133 }
00134 
00135 template <typename T> 
00136 inline Sacado::Tay::CacheTaylor<T>& 
00137 Sacado::Tay::CacheTaylor<T>::operator /= (const T& v)
00138 {
00139   this->coeff_ /= v;
00140 
00141   return *this;
00142 }
00143 
00144 template <typename T> 
00145 template <typename S> 
00146 inline Sacado::Tay::CacheTaylor<T>& 
00147 Sacado::Tay::CacheTaylor<T>::operator += (const S& x)
00148 {
00149   unsigned int xd = x.degree();
00150   unsigned int d = this->degree();
00151 
00152    // Resize polynomial for "this" if x has greater degree
00153   if (xd > d) {
00154     this->resizeCoeffs(xd);
00155     d = xd;
00156   }
00157 
00158   x.allocateCache(d);
00159 
00160   if (x.hasFastAccess(d))
00161     for (int i=d; i>=0; --i)
00162       this->coeff_[i] += x.fastAccessCoeff(i);
00163   else
00164     for (int i=xd; i>=0; --i)
00165       this->coeff_[i] += x.coeff(i);
00166 
00167   return *this;
00168 }
00169 
00170 template <typename T> 
00171 template <typename S> 
00172 inline Sacado::Tay::CacheTaylor<T>& 
00173 Sacado::Tay::CacheTaylor<T>::operator -= (const S& x)
00174 {
00175   unsigned int xd = x.degree();
00176   unsigned int d = this->degree();
00177 
00178    // Resize polynomial for "this" if x has greater degree
00179   if (xd > d) {
00180     this->resizeCoeffs(xd);
00181     d = xd;
00182   }
00183 
00184   x.allocateCache(d);
00185 
00186   if (x.hasFastAccess(d))
00187     for (int i=d; i>=0; --i)
00188       this->coeff_[i] -= x.fastAccessCoeff(i);
00189   else
00190     for (int i=xd; i>=0; --i)
00191       this->coeff_[i] -= x.coeff(i);
00192 
00193   return *this;
00194 }
00195 
00196 template <typename T> 
00197 template <typename S> 
00198 inline Sacado::Tay::CacheTaylor<T>& 
00199 Sacado::Tay::CacheTaylor<T>::operator *= (const S& x)
00200 {
00201   unsigned int xd = x.degree();
00202   unsigned int d = this->degree();
00203   unsigned int dfinal = d;
00204 
00205    // Resize polynomial for "this" if x has greater degree
00206   if (xd > d) {
00207     this->resizeCoeffs(xd);
00208     dfinal = xd;
00209   }
00210 
00211   x.allocateCache(dfinal);
00212 
00213   if (xd) {
00214     if (d) {
00215       T tmp;
00216       if (x.hasFastAccess(dfinal))
00217   for(int i=dfinal; i>=0; --i) {
00218     tmp = T(0.);
00219     for (int k=0; k<=i; ++k)
00220       tmp += this->coeff_[k]*x.fastAccessCoeff(i-k);
00221     this->coeff_[i] = tmp;
00222   }
00223       else
00224   for(int i=dfinal; i>=0; --i) {
00225     tmp = T(0.);
00226     for (int k=0; k<=i; ++k)
00227       tmp += this->coeff_[k]*x.coeff(i-k);
00228     this->coeff_[i] = tmp;
00229   }
00230     }
00231     else {
00232       if (x.hasFastAccess(dfinal))
00233   for(int i=dfinal; i>=0; --i)
00234     this->coeff_[i] = this->coeff_[0] * x.fastAccessCoeff(i);
00235       else
00236   for(int i=dfinal; i>=0; --i)
00237     this->coeff_[i] = this->coeff_[0] * x.coeff(i);
00238     }
00239   }
00240   else 
00241     this->coeff_ *= x.coeff(0);
00242 
00243   return *this;
00244 }
00245 
00246 template <typename T>
00247 template <typename S> 
00248 inline Sacado::Tay::CacheTaylor<T>& 
00249 Sacado::Tay::CacheTaylor<T>::operator /= (const S& x)
00250 {
00251   unsigned int xd = x.degree();
00252   unsigned int d = this->degree();
00253   unsigned int dfinal = d;
00254 
00255    // Resize polynomial for "this" if x has greater degree
00256   if (xd > d) {
00257     this->resizeCoeffs(xd);
00258     dfinal = xd;
00259   }
00260 
00261   x.allocateCache(dfinal);
00262 
00263   if (xd) {
00264     std::valarray<T> tmp(this->coeff_);
00265     if (x.hasFastAccess(dfinal))
00266       for(unsigned int i=0; i<=dfinal; i++) {
00267   for (unsigned int k=1; k<=i; k++)
00268     tmp[i] -= x.fastAccessCoeff(k)*tmp[i-k];
00269   tmp[i] /= x.fastAccessCoeff(0);
00270       }
00271     else
00272       for(unsigned int i=0; i<=dfinal; i++) {
00273   for (unsigned int k=1; k<=i; k++)
00274     tmp[i] -= x.coeff(k)*tmp[i-k];
00275   tmp[i] /= x.coeff(0);
00276       }
00277     this->coeff_ = tmp;
00278   }
00279   else
00280     this->coeff_ /= x.coeff(0);
00281 
00282   return *this;
00283 }
00284 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines