Sacado Package Browser (Single Doxygen Collection) Version of the Day
Sacado_ELRCacheFad_GeneralFadImp.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 //
00031 // The forward-mode AD classes in Sacado are a derivative work of the
00032 // expression template classes in the Fad package by Nicolas Di Cesare.  
00033 // The following banner is included in the original Fad source code:
00034 //
00035 // ************ DO NOT REMOVE THIS BANNER ****************
00036 //
00037 //  Nicolas Di Cesare <Nicolas.Dicesare@ann.jussieu.fr>
00038 //  http://www.ann.jussieu.fr/~dicesare
00039 //
00040 //            CEMRACS 98 : C++ courses, 
00041 //         templates : new C++ techniques 
00042 //            for scientific computing 
00043 // 
00044 //********************************************************
00045 //
00046 //  A short implementation ( not all operators and 
00047 //  functions are overloaded ) of 1st order Automatic
00048 //  Differentiation in forward mode (FAD) using
00049 //  EXPRESSION TEMPLATES.
00050 //
00051 //********************************************************
00052 // @HEADER
00053 
00054 #include "Sacado_ConfigDefs.h"
00055 #include "Sacado_mpl_range_c.hpp"
00056 #include "Sacado_mpl_for_each.hpp"
00057 
00058 // template <typename T, typename Storage> 
00059 // template <typename S> 
00060 // inline Sacado::ELRCacheFad::GeneralFad<T,Storage>::GeneralFad(const Expr<S>& x) :
00061 //   Storage(T(0.))
00062 // {
00063 //   int sz = x.size();
00064 
00065 //   if (sz != this->size()) 
00066 //     this->resize(sz);
00067 
00068 //   if (sz) {
00069 
00070 //     // Number of arguments
00071 //     const int N = Expr<S>::num_args;
00072 
00073 //     // Array to store partials
00074 //     T partials[N];
00075 
00076 //     // Array to store tangents
00077 //     T dots[N];
00078 
00079 //     // Compute partials
00080 //     x.computePartials(T(1.), partials);
00081   
00082 //     // Compute each tangent direction
00083 //     for(int i=0; i<sz; ++i) {
00084 //       T t = T(0.);
00085 //       x.getTangents(i,dots);
00086 //       for (int j=0; j<N; j++) {
00087 //  t += partials[j]*dots[j];
00088 //       }
00089 //       this->fastAccessDx(i) = t;
00090 //     }
00091 
00092 //   }
00093   
00094 //   // Compute value
00095 //   this->val() = x.val();
00096 // }
00097 
00098 template <typename T, typename Storage> 
00099 template <typename S> 
00100 inline Sacado::ELRCacheFad::GeneralFad<T,Storage>::
00101 GeneralFad(const Expr<S>& x) :
00102   Storage(T(0.)),
00103   update_val_(x.updateValue())
00104 {
00105   x.cache();
00106 
00107   int sz = x.size();
00108 
00109   // Compute value
00110   if (update_val_)
00111     this->val() = x.val();
00112 
00113   if (sz != this->size()) 
00114     this->resize(sz);
00115 
00116   if (sz) {
00117 
00118     if (Expr<S>::is_linear) {
00119       if (x.hasFastAccess())
00120         for(int i=0; i<sz; ++i)
00121           this->fastAccessDx(i) = x.fastAccessDx(i);
00122       else
00123         for(int i=0; i<sz; ++i)
00124           this->fastAccessDx(i) = x.dx(i);
00125     }
00126     else {
00127 
00128       // Number of arguments
00129       const int N = Expr<S>::num_args;
00130 
00131       if (x.hasFastAccess()) {
00132   // Compute partials
00133   FastLocalAccumOp< Expr<S> > op(x);
00134   
00135   // Compute each tangent direction
00136   for(op.i=0; op.i<sz; ++op.i) {
00137     op.t = T(0.);
00138 
00139     // Automatically unrolled loop that computes
00140     // for (int j=0; j<N; j++)
00141     //   op.t += op.partials[j] * x.getTangent<j>(i);
00142     Sacado::mpl::for_each< mpl::range_c< int, 0, N > > f(op);
00143     
00144     this->fastAccessDx(op.i) = op.t;
00145   }
00146       }
00147       else {
00148   // Compute partials
00149   SlowLocalAccumOp< Expr<S> > op(x);
00150   
00151   // Compute each tangent direction
00152   for(op.i=0; op.i<sz; ++op.i) {
00153     op.t = T(0.);
00154 
00155     // Automatically unrolled loop that computes
00156     // for (int j=0; j<N; j++)
00157     //   op.t += op.partials[j] * x.getTangent<j>(i);
00158     Sacado::mpl::for_each< mpl::range_c< int, 0, N > > f(op);
00159     
00160     this->fastAccessDx(op.i) = op.t;
00161   }
00162       }
00163 
00164     }
00165   }
00166 }
00167 
00168 template <typename T, typename Storage> 
00169 inline void 
00170 Sacado::ELRCacheFad::GeneralFad<T,Storage>::
00171 diff(const int ith, const int n) 
00172 { 
00173   if (this->size() != n) 
00174     this->resize(n);
00175 
00176   this->zero();
00177   this->fastAccessDx(ith) = T(1.);
00178 
00179 }
00180 
00181 template <typename T, typename Storage> 
00182 inline Sacado::ELRCacheFad::GeneralFad<T,Storage>& 
00183 Sacado::ELRCacheFad::GeneralFad<T,Storage>::
00184 operator=(const T& v) 
00185 {
00186   this->val() = v;
00187 
00188   if (this->size()) {
00189     this->zero();    // We need to zero out the array for future resizes
00190     this->resize(0);
00191   }
00192 
00193   return *this;
00194 }
00195 
00196 template <typename T, typename Storage> 
00197 inline Sacado::ELRCacheFad::GeneralFad<T,Storage>& 
00198 Sacado::ELRCacheFad::GeneralFad<T,Storage>::
00199 operator=(const Sacado::ELRCacheFad::GeneralFad<T,Storage>& x) 
00200 {
00201   // Copy value & dx_
00202   Storage::operator=(x);
00203   update_val_ = x.update_val_;
00204   
00205   return *this;
00206 }
00207 
00208 template <typename T, typename Storage> 
00209 template <typename S> 
00210 inline Sacado::ELRCacheFad::GeneralFad<T,Storage>& 
00211 Sacado::ELRCacheFad::GeneralFad<T,Storage>::
00212 operator=(const Expr<S>& x) 
00213 {
00214   x.cache();
00215 
00216   int sz = x.size();
00217 
00218   if (sz != this->size()) 
00219     this->resize(sz);
00220 
00221   if (sz) {
00222 
00223     if (Expr<S>::is_linear) {
00224       if (x.hasFastAccess())
00225         for(int i=0; i<sz; ++i)
00226           this->fastAccessDx(i) = x.fastAccessDx(i);
00227       else
00228         for(int i=0; i<sz; ++i)
00229           this->fastAccessDx(i) = x.dx(i);
00230     }
00231     else {
00232 
00233       // Number of arguments
00234       const int N = Expr<S>::num_args;
00235 
00236       if (x.hasFastAccess()) {
00237   // Compute partials
00238   FastLocalAccumOp< Expr<S> > op(x);
00239   
00240   // Compute each tangent direction
00241   for(op.i=0; op.i<sz; ++op.i) {
00242     op.t = T(0.);
00243 
00244     // Automatically unrolled loop that computes
00245     // for (int j=0; j<N; j++)
00246     //   op.t += op.partials[j] * x.getTangent<j>(i);
00247     Sacado::mpl::for_each< mpl::range_c< int, 0, N > > f(op);
00248     
00249     this->fastAccessDx(op.i) = op.t;
00250   }
00251       }
00252       else {
00253   // Compute partials
00254   SlowLocalAccumOp< Expr<S> > op(x);
00255   
00256   // Compute each tangent direction
00257   for(op.i=0; op.i<sz; ++op.i) {
00258     op.t = T(0.);
00259 
00260     // Automatically unrolled loop that computes
00261     // for (int j=0; j<N; j++)
00262     //   op.t += op.partials[j] * x.getTangent<j>(i);
00263     Sacado::mpl::for_each< mpl::range_c< int, 0, N > > f(op);
00264     
00265     this->fastAccessDx(op.i) = op.t;
00266   }
00267       }
00268 
00269     }
00270 
00271   }
00272   
00273   update_val_ = x.updateValue();
00274   if (update_val_)
00275     this->val() = x.val();
00276   
00277   return *this;
00278 }
00279 
00280 template <typename T, typename Storage> 
00281 inline  Sacado::ELRCacheFad::GeneralFad<T,Storage>& 
00282 Sacado::ELRCacheFad::GeneralFad<T,Storage>::
00283 operator += (const T& v)
00284 {
00285   if (update_val_)
00286     this->val() += v;
00287 
00288   return *this;
00289 }
00290 
00291 template <typename T, typename Storage> 
00292 inline Sacado::ELRCacheFad::GeneralFad<T,Storage>& 
00293 Sacado::ELRCacheFad::GeneralFad<T,Storage>::
00294 operator -= (const T& v)
00295 {
00296   if (update_val_)
00297     this->val() -= v;
00298 
00299   return *this;
00300 }
00301 
00302 template <typename T, typename Storage> 
00303 inline Sacado::ELRCacheFad::GeneralFad<T,Storage>& 
00304 Sacado::ELRCacheFad::GeneralFad<T,Storage>::
00305 operator *= (const T& v)
00306 {
00307   int sz = this->size();
00308 
00309   if (update_val_)
00310     this->val() *= v;
00311   for (int i=0; i<sz; ++i)
00312     this->fastAccessDx(i) *= v;
00313 
00314   return *this;
00315 }
00316 
00317 template <typename T, typename Storage> 
00318 inline Sacado::ELRCacheFad::GeneralFad<T,Storage>& 
00319 Sacado::ELRCacheFad::GeneralFad<T,Storage>::
00320 operator /= (const T& v)
00321 {
00322   int sz = this->size();
00323 
00324   if (update_val_)
00325     this->val() /= v;
00326   for (int i=0; i<sz; ++i)
00327     this->fastAccessDx(i) /= v;
00328 
00329   return *this;
00330 }
00331 
00332 template <typename T, typename Storage> 
00333 inline  Sacado::ELRCacheFad::GeneralFad<T,Storage>& 
00334 Sacado::ELRCacheFad::GeneralFad<T,Storage>::
00335 operator += (const typename Sacado::dummy<value_type,scalar_type>::type& v)
00336 {
00337   if (update_val_)
00338     this->val() += v;
00339 
00340   return *this;
00341 }
00342 
00343 template <typename T, typename Storage> 
00344 inline Sacado::ELRCacheFad::GeneralFad<T,Storage>& 
00345 Sacado::ELRCacheFad::GeneralFad<T,Storage>::
00346 operator -= (const typename Sacado::dummy<value_type,scalar_type>::type& v)
00347 {
00348   if (update_val_)
00349     this->val() -= v;
00350 
00351   return *this;
00352 }
00353 
00354 template <typename T, typename Storage> 
00355 inline Sacado::ELRCacheFad::GeneralFad<T,Storage>& 
00356 Sacado::ELRCacheFad::GeneralFad<T,Storage>::
00357 operator *= (const typename Sacado::dummy<value_type,scalar_type>::type& v)
00358 {
00359   int sz = this->size();
00360 
00361   if (update_val_)
00362     this->val() *= v;
00363   for (int i=0; i<sz; ++i)
00364     this->fastAccessDx(i) *= v;
00365 
00366   return *this;
00367 }
00368 
00369 template <typename T, typename Storage> 
00370 inline Sacado::ELRCacheFad::GeneralFad<T,Storage>& 
00371 Sacado::ELRCacheFad::GeneralFad<T,Storage>::
00372 operator /= (const typename Sacado::dummy<value_type,scalar_type>::type& v)
00373 {
00374   int sz = this->size();
00375 
00376   if (update_val_)
00377     this->val() /= v;
00378   for (int i=0; i<sz; ++i)
00379     this->fastAccessDx(i) /= v;
00380 
00381   return *this;
00382 }
00383 
00384 template <typename T, typename Storage> 
00385 template <typename S> 
00386 inline Sacado::ELRCacheFad::GeneralFad<T,Storage>& 
00387 Sacado::ELRCacheFad::GeneralFad<T,Storage>::
00388 operator += (const Sacado::ELRCacheFad::Expr<S>& x)
00389 {
00390   x.cache();
00391 
00392   int xsz = x.size(), sz = this->size();
00393 
00394 #ifdef SACADO_DEBUG
00395   if ((xsz != sz) && (xsz != 0) && (sz != 0))
00396     throw "Fad Error:  Attempt to assign with incompatible sizes";
00397 #endif
00398 
00399   if (Expr<S>::is_linear) {
00400     if (xsz) {
00401       if (sz) {
00402   if (x.hasFastAccess())
00403     for (int i=0; i<sz; ++i)
00404       this->fastAccessDx(i) += x.fastAccessDx(i);
00405   else
00406     for (int i=0; i<sz; ++i)
00407       this->fastAccessDx(i) += x.dx(i);
00408       }
00409       else {
00410   this->resize(xsz);
00411   if (x.hasFastAccess())
00412     for (int i=0; i<xsz; ++i)
00413       this->fastAccessDx(i) = x.fastAccessDx(i);
00414   else
00415     for (int i=0; i<xsz; ++i)
00416       this->fastAccessDx(i) = x.dx(i);
00417       }
00418     }
00419   }
00420   else {
00421 
00422     // Number of arguments
00423     const int N = Expr<S>::num_args;
00424 
00425     if (xsz) {
00426 
00427       if (sz != xsz) {
00428   this->resize(xsz);
00429   this->zero();
00430       }
00431 
00432       if (x.hasFastAccess()) {
00433   // Compute partials
00434   FastLocalAccumOp< Expr<S> > op(x);
00435   
00436   // Compute each tangent direction
00437   for(op.i=0; op.i<xsz; ++op.i) {
00438     op.t = T(0.);
00439     
00440     // Automatically unrolled loop that computes
00441     // for (int j=0; j<N; j++)
00442     //   op.t += op.partials[j] * x.getTangent<j>(i);
00443     Sacado::mpl::for_each< mpl::range_c< int, 0, N > > f(op);
00444     
00445     this->fastAccessDx(op.i) += op.t;
00446   }
00447       }
00448       else {
00449   // Compute partials
00450   SlowLocalAccumOp< Expr<S> > op(x);
00451   
00452   // Compute each tangent direction
00453   for(op.i=0; op.i<xsz; ++op.i) {
00454     op.t = T(0.);
00455 
00456     // Automatically unrolled loop that computes
00457     // for (int j=0; j<N; j++)
00458     //   op.t += op.partials[j] * x.getTangent<j>(i);
00459     Sacado::mpl::for_each< mpl::range_c< int, 0, N > > f(op);
00460     
00461     this->fastAccessDx(op.i) += op.t;
00462   }
00463       }
00464 
00465     }
00466 
00467   }
00468   
00469   update_val_ = x.updateValue();
00470   if (update_val_)
00471     this->val() += x.val();
00472 
00473   return *this;
00474 }
00475 
00476 template <typename T, typename Storage> 
00477 template <typename S> 
00478 inline Sacado::ELRCacheFad::GeneralFad<T,Storage>& 
00479 Sacado::ELRCacheFad::GeneralFad<T,Storage>::
00480 operator -= (const Sacado::ELRCacheFad::Expr<S>& x)
00481 {
00482   x.cache();
00483 
00484   int xsz = x.size(), sz = this->size();
00485 
00486 #ifdef SACADO_DEBUG
00487   if ((xsz != sz) && (xsz != 0) && (sz != 0))
00488     throw "Fad Error:  Attempt to assign with incompatible sizes";
00489 #endif
00490 
00491   if (Expr<S>::is_linear) {
00492     if (xsz) {
00493       if (sz) {
00494   if (x.hasFastAccess())
00495     for(int i=0; i<sz; ++i)
00496       this->fastAccessDx(i) -= x.fastAccessDx(i);
00497   else
00498     for (int i=0; i<sz; ++i)
00499       this->fastAccessDx(i) -= x.dx(i);
00500       }
00501       else {
00502   this->resize(xsz);
00503   if (x.hasFastAccess())
00504     for(int i=0; i<xsz; ++i)
00505       this->fastAccessDx(i) = -x.fastAccessDx(i);
00506   else
00507     for (int i=0; i<xsz; ++i)
00508       this->fastAccessDx(i) = -x.dx(i);
00509       }
00510     }
00511   }
00512   else {
00513 
00514     // Number of arguments
00515     const int N = Expr<S>::num_args;
00516 
00517     if (xsz) {
00518 
00519       if (sz != xsz) {
00520   this->resize(xsz);
00521   this->zero();
00522       }
00523 
00524       if (x.hasFastAccess()) {
00525   // Compute partials
00526   FastLocalAccumOp< Expr<S> > op(x);
00527   
00528   // Compute each tangent direction
00529   for(op.i=0; op.i<xsz; ++op.i) {
00530     op.t = T(0.);
00531     
00532     // Automatically unrolled loop that computes
00533     // for (int j=0; j<N; j++)
00534     //   op.t += op.partials[j] * x.getTangent<j>(i);
00535     Sacado::mpl::for_each< mpl::range_c< int, 0, N > > f(op);
00536     
00537     this->fastAccessDx(op.i) -= op.t;
00538   }
00539       }
00540       else {
00541   // Compute partials
00542   SlowLocalAccumOp< Expr<S> > op(x);
00543   
00544   // Compute each tangent direction
00545   for(op.i=0; op.i<xsz; ++op.i) {
00546     op.t = T(0.);
00547 
00548     // Automatically unrolled loop that computes
00549     // for (int j=0; j<N; j++)
00550     //   op.t += op.partials[j] * x.getTangent<j>(i);
00551     Sacado::mpl::for_each< mpl::range_c< int, 0, N > > f(op);
00552     
00553     this->fastAccessDx(op.i) -= op.t;
00554   } 
00555       }
00556 
00557     }
00558 
00559   }
00560 
00561   update_val_ = x.updateValue();
00562   if (update_val_)
00563     this->val() -= x.val();
00564 
00565   return *this;
00566 }
00567 
00568 template <typename T, typename Storage> 
00569 template <typename S> 
00570 inline Sacado::ELRCacheFad::GeneralFad<T,Storage>& 
00571 Sacado::ELRCacheFad::GeneralFad<T,Storage>::
00572 operator *= (const Sacado::ELRCacheFad::Expr<S>& x)
00573 {
00574   x.cache();
00575 
00576   int xsz = x.size(), sz = this->size();
00577   update_val_ = x.updateValue();
00578   T xval = x.val();
00579   T v = this->val();
00580 
00581 #ifdef SACADO_DEBUG
00582   if ((xsz != sz) && (xsz != 0) && (sz != 0))
00583     throw "Fad Error:  Attempt to assign with incompatible sizes";
00584 #endif
00585 
00586   if (Expr<S>::is_linear) {
00587     if (xsz) {
00588       if (sz) {
00589   if (x.hasFastAccess())
00590     for(int i=0; i<sz; ++i)
00591       this->fastAccessDx(i) = v*x.fastAccessDx(i) + this->fastAccessDx(i)*xval;
00592   else
00593     for (int i=0; i<sz; ++i)
00594       this->fastAccessDx(i) = v*x.dx(i) + this->fastAccessDx(i)*xval;
00595       }
00596       else {
00597   this->resize(xsz);
00598   if (x.hasFastAccess())
00599     for(int i=0; i<xsz; ++i)
00600       this->fastAccessDx(i) = v*x.fastAccessDx(i);
00601   else
00602     for (int i=0; i<xsz; ++i)
00603       this->fastAccessDx(i) = v*x.dx(i);
00604       }
00605     }
00606     else {
00607       if (sz) {
00608   for (int i=0; i<sz; ++i)
00609     this->fastAccessDx(i) *= xval;
00610       }
00611     }
00612   }
00613   else {
00614 
00615     // Number of arguments
00616     const int N = Expr<S>::num_args;
00617 
00618     if (xsz) {
00619 
00620       if (sz) {
00621 
00622   if (x.hasFastAccess()) {
00623     // Compute partials
00624     FastLocalAccumOp< Expr<S> > op(x);
00625     
00626     // Compute each tangent direction
00627     for(op.i=0; op.i<xsz; ++op.i) {
00628       op.t = T(0.);
00629       
00630       // Automatically unrolled loop that computes
00631       // for (int j=0; j<N; j++)
00632       //   op.t += op.partials[j] * x.getTangent<j>(i);
00633       Sacado::mpl::for_each< mpl::range_c< int, 0, N > > f(op);
00634     
00635       this->fastAccessDx(op.i) = 
00636         v * op.t + this->fastAccessDx(op.i) * xval;
00637     }
00638   }
00639   else {
00640     // Compute partials
00641     SlowLocalAccumOp< Expr<S> > op(x);
00642     
00643     // Compute each tangent direction
00644     for(op.i=0; op.i<xsz; ++op.i) {
00645       op.t = T(0.);
00646       
00647       // Automatically unrolled loop that computes
00648       // for (int j=0; j<N; j++)
00649       //   op.t += op.partials[j] * x.getTangent<j>(i);
00650       Sacado::mpl::for_each< mpl::range_c< int, 0, N > > f(op);
00651     
00652       this->fastAccessDx(op.i) = 
00653         v * op.t + this->fastAccessDx(op.i) * xval;
00654     }
00655   }
00656 
00657       }
00658       
00659       else {
00660   
00661   this->resize(xsz);
00662   
00663   if (x.hasFastAccess()) {
00664     // Compute partials
00665     FastLocalAccumOp< Expr<S> > op(x);
00666     
00667     // Compute each tangent direction
00668     for(op.i=0; op.i<xsz; ++op.i) {
00669       op.t = T(0.);
00670       
00671       // Automatically unrolled loop that computes
00672       // for (int j=0; j<N; j++)
00673       //   op.t += op.partials[j] * x.getTangent<j>(i);
00674       Sacado::mpl::for_each< mpl::range_c< int, 0, N > > f(op);
00675     
00676       this->fastAccessDx(op.i) = v * op.t;
00677     }
00678   }
00679   else {
00680     // Compute partials
00681     SlowLocalAccumOp< Expr<S> > op(x);
00682     
00683     // Compute each tangent direction
00684     for(op.i=0; op.i<xsz; ++op.i) {
00685       op.t = T(0.);
00686       
00687       // Automatically unrolled loop that computes
00688       // for (int j=0; j<N; j++)
00689       //   op.t += op.partials[j] * x.getTangent<j>(i);
00690       Sacado::mpl::for_each< mpl::range_c< int, 0, N > > f(op);
00691     
00692       this->fastAccessDx(op.i) = v * op.t;
00693     }
00694   }
00695   
00696       }
00697 
00698     }
00699 
00700     else {
00701 
00702       if (sz) {
00703   for (int i=0; i<sz; ++i)
00704     this->fastAccessDx(i) *= xval;
00705       }
00706 
00707     }
00708 
00709   }
00710 
00711   if (update_val_)
00712     this->val() *= xval;
00713 
00714   return *this;
00715 }
00716 
00717 template <typename T, typename Storage>
00718 template <typename S> 
00719 inline Sacado::ELRCacheFad::GeneralFad<T,Storage>& 
00720 Sacado::ELRCacheFad::GeneralFad<T,Storage>::
00721 operator /= (const Sacado::ELRCacheFad::Expr<S>& x)
00722 {
00723   x.cache();
00724 
00725   int xsz = x.size(), sz = this->size();
00726   update_val_ = x.updateValue();
00727   T xval = x.val();
00728   T v = this->val();
00729 
00730 #ifdef SACADO_DEBUG
00731   if ((xsz != sz) && (xsz != 0) && (sz != 0))
00732     throw "Fad Error:  Attempt to assign with incompatible sizes";
00733 #endif
00734 
00735   if (Expr<S>::is_linear) {
00736     if (xsz) {
00737       if (sz) {
00738   if (x.hasFastAccess())
00739     for(int i=0; i<sz; ++i)
00740       this->fastAccessDx(i) = ( this->fastAccessDx(i)*xval - v*x.fastAccessDx(i) )/ (xval*xval);
00741   else
00742     for (int i=0; i<sz; ++i)
00743       this->fastAccessDx(i) = ( this->fastAccessDx(i)*xval - v*x.dx(i) )/ (xval*xval);
00744       }
00745       else {
00746   this->resize(xsz);
00747   if (x.hasFastAccess())
00748     for(int i=0; i<xsz; ++i)
00749       this->fastAccessDx(i) = - v*x.fastAccessDx(i) / (xval*xval);
00750   else
00751     for (int i=0; i<xsz; ++i)
00752       this->fastAccessDx(i) = -v*x.dx(i) / (xval*xval);
00753       }
00754     }
00755     else {
00756       if (sz) {
00757   for (int i=0; i<sz; ++i)
00758     this->fastAccessDx(i) /= xval;
00759       }
00760     }
00761   }
00762   else {
00763 
00764     // Number of arguments
00765     const int N = Expr<S>::num_args;
00766 
00767     if (xsz) {
00768 
00769       T xval2 = xval*xval;
00770       
00771       if (sz) {
00772 
00773   if (x.hasFastAccess()) {
00774     // Compute partials
00775     FastLocalAccumOp< Expr<S> > op(x);
00776     
00777     // Compute each tangent direction
00778     for(op.i=0; op.i<xsz; ++op.i) {
00779       op.t = T(0.);
00780       
00781       // Automatically unrolled loop that computes
00782       // for (int j=0; j<N; j++)
00783       //   op.t += op.partials[j] * x.getTangent<j>(i);
00784       Sacado::mpl::for_each< mpl::range_c< int, 0, N > > f(op);
00785     
00786       this->fastAccessDx(op.i) = 
00787         (this->fastAccessDx(op.i) * xval - v * op.t) / xval2;
00788     }
00789   }
00790   else {
00791     // Compute partials
00792     SlowLocalAccumOp< Expr<S> > op(x);
00793     
00794     // Compute each tangent direction
00795     for(op.i=0; op.i<xsz; ++op.i) {
00796       op.t = T(0.);
00797       
00798       // Automatically unrolled loop that computes
00799       // for (int j=0; j<N; j++)
00800       //   op.t += op.partials[j] * x.getTangent<j>(i);
00801       Sacado::mpl::for_each< mpl::range_c< int, 0, N > > f(op);
00802     
00803       this->fastAccessDx(op.i) = 
00804         (this->fastAccessDx(op.i) * xval - v * op.t) / xval2;
00805     }
00806   }
00807   
00808       }
00809       
00810       else {
00811   
00812   this->resize(xsz);
00813   
00814   if (x.hasFastAccess()) {
00815     // Compute partials
00816     FastLocalAccumOp< Expr<S> > op(x);
00817     
00818     // Compute each tangent direction
00819     for(op.i=0; op.i<xsz; ++op.i) {
00820       op.t = T(0.);
00821       
00822       // Automatically unrolled loop that computes
00823       // for (int j=0; j<N; j++)
00824       //   op.t += op.partials[j] * x.getTangent<j>(i);
00825       Sacado::mpl::for_each< mpl::range_c< int, 0, N > > f(op);
00826     
00827       this->fastAccessDx(op.i) = (-v * op.t) / xval2;
00828     }
00829   }
00830   else {
00831     // Compute partials
00832     SlowLocalAccumOp< Expr<S> > op(x);
00833     
00834     // Compute each tangent direction
00835     for(op.i=0; op.i<xsz; ++op.i) {
00836       op.t = T(0.);
00837       
00838       // Automatically unrolled loop that computes
00839       // for (int j=0; j<N; j++)
00840       //   op.t += op.partials[j] * x.getTangent<j>(i);
00841       Sacado::mpl::for_each< mpl::range_c< int, 0, N > > f(op);
00842     
00843       this->fastAccessDx(op.i) = (-v * op.t) / xval2;
00844     }
00845   }
00846   
00847       }
00848 
00849     }
00850 
00851     else {
00852 
00853       if (sz) {
00854   for (int i=0; i<sz; ++i)
00855     this->fastAccessDx(i) /= xval;
00856       }
00857 
00858     }
00859 
00860   }
00861 
00862   if (update_val_)
00863     this->val() /= xval;
00864 
00865   return *this;
00866 }
00867 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines