Sacado Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
Sacado_Tay_CacheTaylorImp.hpp
Go to the documentation of this file.
1// @HEADER
2// ***********************************************************************
3//
4// Sacado Package
5// Copyright (2006) Sandia Corporation
6//
7// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
8// the U.S. Government retains certain rights in this software.
9//
10// This library is free software; you can redistribute it and/or modify
11// it under the terms of the GNU Lesser General Public License as
12// published by the Free Software Foundation; either version 2.1 of the
13// License, or (at your option) any later version.
14//
15// This library is distributed in the hope that it will be useful, but
16// WITHOUT ANY WARRANTY; without even the implied warranty of
17// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
18// Lesser General Public License for more details.
19//
20// You should have received a copy of the GNU Lesser General Public
21// License along with this library; if not, write to the Free Software
22// Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
23// USA
24// Questions? Contact David M. Gay (dmgay@sandia.gov) or Eric T. Phipps
25// (etphipp@sandia.gov).
26//
27// ***********************************************************************
28// @HEADER
29
30template <typename T>
31template <typename S>
33 Expr< CacheTaylorImplementation<T> >(x.degree(), T(0.))
34{
35 int d = this->degree();
36
37 x.allocateCache(d);
38
39 // We copy the coefficients from the highest degree to the lowest just
40 // to be consistent with operator=(), even though it is not strictly
41 // necessary since "this" cannot be on the RHS in the copy constructor.
42 if (x.hasFastAccess(d))
43 for(int i=d; i>=0; --i)
44 this->coeff_[i] = x.fastAccessCoeff(i);
45 else
46 for(int i=d; i>=0; --i)
47 this->coeff_[i] = x.coeff(i);
48
49}
50
51template <typename T>
54{
55 this->coeff_[0] = v;
56
57 for (int i=1; i<this->coeff_size(); i++)
58 this->coeff_[i] = T(0.);
59
60 return *this;
61}
62
63template <typename T>
66{
67 if (x.coeff_size() != this->coeff_size())
68 this->coeff_.resize(x.coeff_size());
69 this->coeff_ = x.coeff_;
70
71 return *this;
72}
73
74template <typename T>
75template <typename S>
78{
79 int d = this->degree();
80 int xd = x.degree();
81
82 // Resize polynomial for "this" if x has greater degree
83 if (xd > d) {
84 this->coeff_.resize(xd+1);
85 d = xd;
86 }
87
88 x.allocateCache(d);
89
90 // Copy coefficients. Note: we copy from the last term down to the first
91 // to take into account "this" being involved in an expression on the RHS.
92 // By overwriting the degree k term, we are guarranteed not to affect any
93 // of the lower degree terms. However, this in general will force each
94 // term in the expression to compute all of its coefficients at once instead
95 // traversing the expression once for each degree.
96 if (x.hasFastAccess(d))
97 for(int i=d; i>=0; --i)
98 this->coeff_[i] = x.fastAccessCoeff(i);
99 else
100 for(int i=d; i>=0; --i)
101 this->coeff_[i] = x.coeff(i);
102
103 return *this;
104}
105
106template <typename T>
109{
110 this->coeff_[0] += v;
111
112 return *this;
113}
114
115template <typename T>
118{
119 this->coeff_[0] -= v;
120
121 return *this;
122}
123
124template <typename T>
127{
128 this->coeff_ *= v;
129
130 return *this;
131}
132
133template <typename T>
136{
137 this->coeff_ /= v;
138
139 return *this;
140}
141
142template <typename T>
143template <typename S>
146{
147 int xd = x.degree();
148 int d = this->degree();
149
150 // Resize polynomial for "this" if x has greater degree
151 if (xd > d) {
152 this->resizeCoeffs(xd);
153 d = xd;
154 }
155
156 x.allocateCache(d);
157
158 if (x.hasFastAccess(d))
159 for (int i=d; i>=0; --i)
160 this->coeff_[i] += x.fastAccessCoeff(i);
161 else
162 for (int i=xd; i>=0; --i)
163 this->coeff_[i] += x.coeff(i);
164
165 return *this;
166}
167
168template <typename T>
169template <typename S>
172{
173 int xd = x.degree();
174 int d = this->degree();
175
176 // Resize polynomial for "this" if x has greater degree
177 if (xd > d) {
178 this->resizeCoeffs(xd);
179 d = xd;
180 }
181
182 x.allocateCache(d);
183
184 if (x.hasFastAccess(d))
185 for (int i=d; i>=0; --i)
186 this->coeff_[i] -= x.fastAccessCoeff(i);
187 else
188 for (int i=xd; i>=0; --i)
189 this->coeff_[i] -= x.coeff(i);
190
191 return *this;
192}
193
194template <typename T>
195template <typename S>
198{
199 int xd = x.degree();
200 int d = this->degree();
201 int dfinal = d;
202
203 // Resize polynomial for "this" if x has greater degree
204 if (xd > d) {
205 this->resizeCoeffs(xd);
206 dfinal = xd;
207 }
208
209 x.allocateCache(dfinal);
210
211 if (xd) {
212 if (d) {
213 T tmp;
214 if (x.hasFastAccess(dfinal))
215 for(int i=dfinal; i>=0; --i) {
216 tmp = T(0.);
217 for (int k=0; k<=i; ++k)
218 tmp += this->coeff_[k]*x.fastAccessCoeff(i-k);
219 this->coeff_[i] = tmp;
220 }
221 else
222 for(int i=dfinal; i>=0; --i) {
223 tmp = T(0.);
224 for (int k=0; k<=i; ++k)
225 tmp += this->coeff_[k]*x.coeff(i-k);
226 this->coeff_[i] = tmp;
227 }
228 }
229 else {
230 if (x.hasFastAccess(dfinal))
231 for(int i=dfinal; i>=0; --i)
232 this->coeff_[i] = this->coeff_[0] * x.fastAccessCoeff(i);
233 else
234 for(int i=dfinal; i>=0; --i)
235 this->coeff_[i] = this->coeff_[0] * x.coeff(i);
236 }
237 }
238 else
239 this->coeff_ *= x.coeff(0);
240
241 return *this;
242}
243
244template <typename T>
245template <typename S>
248{
249 int xd = x.degree();
250 int d = this->degree();
251 int dfinal = d;
252
253 // Resize polynomial for "this" if x has greater degree
254 if (xd > d) {
255 this->resizeCoeffs(xd);
256 dfinal = xd;
257 }
258
259 x.allocateCache(dfinal);
260
261 if (xd) {
262 std::valarray<T> tmp(this->coeff_);
263 if (x.hasFastAccess(dfinal))
264 for(int i=0; i<=dfinal; i++) {
265 for (int k=1; k<=i; k++)
266 tmp[i] -= x.fastAccessCoeff(k)*tmp[i-k];
267 tmp[i] /= x.fastAccessCoeff(0);
268 }
269 else
270 for(int i=0; i<=dfinal; i++) {
271 for (int k=1; k<=i; k++)
272 tmp[i] -= x.coeff(k)*tmp[i-k];
273 tmp[i] /= x.coeff(0);
274 }
275 this->coeff_ = tmp;
276 }
277 else
278 this->coeff_ /= x.coeff(0);
279
280 return *this;
281}
282
#define T
Definition: Sacado_rad.hpp:573
Taylor polynomial class using caching expression templates.
int degree() const
Returns degree of polynomial.
std::valarray< T > coeff_
Taylor polynomial coefficients.
Forward-mode AD class using dynamic memory allocation.
CacheTaylor< T > & operator+=(const T &x)
Addition-assignment operator with constant right-hand-side.
CacheTaylor< T > & operator=(const T &v)
Assignment operator with constant right-hand-side.
CacheTaylor< T > & operator/=(const T &x)
Division-assignment operator with constant right-hand-side.
CacheTaylor< T > & operator-=(const T &x)
Subtraction-assignment operator with constant right-hand-side.
CacheTaylor()
Default constructor.
CacheTaylor< T > & operator*=(const T &x)
Multiplication-assignment operator with constant right-hand-side.
Wrapper for a generic expression template.
void allocateCache(unsigned int d) const
Allocate coefficient cache.
value_type fastAccessCoeff(unsigned int i) const
Return degree i term of expression.
value_type coeff(unsigned int i) const
Return degree i term of expression.