OPAL (Object Oriented Parallel Accelerator Library) 2022.1
OPAL
VpsMap.h
Go to the documentation of this file.
1#ifndef CLASSIC_VpsMap_HH
2#define CLASSIC_VpsMap_HH
3
4// ------------------------------------------------------------------------
5// $RCSfile: VpsMap.h,v $
6// ------------------------------------------------------------------------
7// $Revision: 1.1.1.1.2.1 $
8// ------------------------------------------------------------------------
9// Copyright: see Copyright.readme
10// ------------------------------------------------------------------------
11//
12// Class: VpsMap
13//
14// ------------------------------------------------------------------------
15// Class category: Algebra
16// ------------------------------------------------------------------------
17//
18// $Date: 2004/11/18 22:18:06 $
19// $Author: jsberg $
20//
21// ------------------------------------------------------------------------
22
23#include "Algebra/LUMatrix.h"
24#include "Algebra/Matrix.h"
25#include "Algebra/Tps.h"
26#include "Algebra/TpsData.h"
27#include "Algebra/TpsMonomial.h"
29#include "Algebra/Vector.h"
30#include "Algebra/Vps.h"
31#include "Algebra/VpsMap.h"
32#include "Utilities/SizeError.h"
33#include <iosfwd>
34
35
36// Template class VpsMap<T>, declaration.
37// ------------------------------------------------------------------------
39// This class includes substitution operations.
40// The input and output dimensions need not be the same.
41
42template <class T>
43class VpsMap: public Vps<T> {
44
45public:
46
48 // Construct [b]nDim[/b] series in [b]nVar[/b] variables each.
49 VpsMap(int nDim, int nVar);
50
52 VpsMap(const Vps<T> &rhs);
53
55 // Assign the [b]nDim[/b] components from [b]rhs[/b].
56 VpsMap(int nDim, const Tps<T> rhs[]);
57
59 // The constant terms are zero.
60 // The linear terms are taken from [b]M[/b].
61 VpsMap(const Matrix<T> &M);
62
64 // The constant terms are taken from [b]V[/b].
65 // The linear terms are zero.
66 VpsMap(const Vector<T> &V);
67
68 VpsMap();
69 VpsMap(const VpsMap<T> &rhs);
70 ~VpsMap();
71 VpsMap<T> &operator=(const VpsMap<T> &y);
72
74 VpsMap<T> substitute(const VpsMap<T> &vv) const;
75
77 VpsMap<T> substitute(const Matrix<T> &M) const;
78
80 VpsMap<T> substitute(const VpsMap<T> &y, int trunc) const;
81
83 VpsMap<T> substituteInto(const Matrix<T> &x) const;
84
86 // Return the point [b]y[/b] mapped by the map.
87 Vector<T> constantTerm(const Vector<T> &y) const;
88
90 // Extract constant part.
91 Vector<T> constantTerm() const;
92
94 // Return the linear terms of the map around the point [b]y[/b],
95 // i.e. the linear part of the partial differentials at [b]y[/b].
96 Matrix<T> linearTerms(const Vector<T> &y) const;
97
99 // Return the linear part.
100 Matrix<T> linearTerms() const;
101
103 VpsMap<T> derivative(int var) const;
104
106 VpsMap<T> integral(int var) const;
107};
108
109
110// Global Operators on VpsMap<T>
111// ------------------------------------------------------------------------
112
114template <class T> VpsMap<T> operator*(const VpsMap<T> &x, const Tps<T> &y);
115
117template <class T> VpsMap<T> operator*(const Tps<T> &x, const VpsMap<T> &y);
118
120template <class T> VpsMap<T> operator*(const VpsMap<T> &x, const T &y);
121
123template <class T> VpsMap<T> operator*(const T &x, const VpsMap<T> &y);
124
126// Throw DivideError if constant part of [b]y[/b] is zero.
127template <class T> VpsMap<T> operator/(const VpsMap<T> &x, const Tps<T> &y);
128
130// Throw DivideError if [b]y[/b] is zero.
131template <class T> VpsMap<T> operator/(const VpsMap<T> &x, const T &y);
132
134template <class T> VpsMap<T> operator+(const VpsMap<T> &x, const VpsMap<T> &y);
135
137template <class T> VpsMap<T> operator-(const VpsMap<T> &x, const VpsMap<T> &y);
138
140template <class T> VpsMap<T> operator+(const VpsMap<T> &x, const Vector<T> &y);
141
143template <class T> VpsMap<T> operator-(const VpsMap<T> &x, const Vector<T> &y);
144
146template <class T> VpsMap<T> operator+(const Vector<T> &x, const VpsMap<T> &y);
147
149template <class T> VpsMap<T> operator-(const Vector<T> &x, const VpsMap<T> &y);
150
152template <class T> std::istream &operator>>(std::istream &, VpsMap<T> &x);
153
155template <class T> std::ostream &operator<<(std::ostream &, const VpsMap<T> &x);
156
158template <class T> inline
160{ return y.substituteInto(x); }
161
162
163// Template class VpsMap<T>, implementation.
164// ------------------------------------------------------------------------
165
166template <class T> inline
168 Vps<T>()
169{}
170
171
172template <class T> inline
173VpsMap<T>::VpsMap(int nDim, int nVar):
174 Vps<T>(nDim, nVar)
175{}
176
177
178template <class T> inline
180 Vps<T>(rhs)
181{}
182
183
184template <class T> inline
186 Vps<T>(rhs)
187{}
188
189
190template <class T> inline
191VpsMap<T>::VpsMap(int nDim, const Tps<T> rhs[]):
192 Vps<T>(nDim, rhs)
193{}
194
195
196template <class T> inline
198 Vps<T>(x)
199{}
200
201
202template <class T> inline
204 Vps<T>(x)
205{}
206
207
208template <class T> inline
210{}
211
212
213template <class T> inline
216 return *this;
217}
218
219
220template <class T> inline
222 return substitute(y, this->getTruncOrder());
223}
224
225
226template <class T> inline
228 return substitute(VpsMap<T>(y));
229}
230
231
232template <class T>
234 if(this->getVariables() != y.size()) {
235 throw SizeError("VpsMap::constantTerm()", "Inconsistent dimensions.");
236 }
237
238 // This statement makes sure that the substitution table exists.
239 const Array1D<TpsSubstitution> &table =
240 TpsData::getTpsData(this->getTopOrder(),
241 this->getVariables())->getSubTable();
242
243 // Copy the constant term.
244 Vector<T> z = constantTerm();
245 Array1D<T> product(table.size());
246 product[0] = T(1);
247
248 for(int next = 1; next < table.size();) {
249 const TpsSubstitution &s = table[next];
250 product[s.order] = product[s.order-1] * y[s.variable];
251
252 for(int v = 0; v < this->getDimension(); v++) {
253 if(this->data[v][s.index] != T(0)) {
254 z[v] += this->data[v][s.index] * product[s.order];
255 }
256 }
257
258 next = (s.order < this->getTopOrder()) ? next + 1 : s.skip;
259 }
260
261 return z;
262}
263
264
265template <class T>
267 Vector<T> z(this->getDimension());
268 for(int i = 0; i < this->getDimension(); i++) z[i] = this->data[i][0];
269 return z;
270}
271
272
273template <class T>
275 Matrix<T> z(this->getDimension(), this->getVariables());
276
277 for(int i = 0; i < z.nrows(); i++) {
278 for(int j = 0; j < z.ncols(); j++) {
279 z[i][j] = this->data[i].derivative(j).evaluate(y);
280 }
281 }
282
283 return z;
284}
285
286
287template <class T>
289 Matrix<T> M(this->getDimension(), this->getVariables());
290
291 for(int i = 0; i < this->getDimension(); i++) {
292 for(int j = 0; j < this->getVariables(); j++) {
293 M(i, j) = this->data[i][j+1];
294 }
295 }
296
297 return M;
298}
299
300
301template <class T>
303 VpsMap<T> z(this->getDimension(), this->getVariables());
304
305 for(int i = 0; i < this->getDimension(); i++) {
306 z[i] = this->data[i].derivative(var);
307 }
308
309 return z;
310}
311
312
313template <class T>
315 VpsMap<T> z(this->getDimension(), this->getVariables());
316
317 for(int i = 0; i < this->getDimension(); i++) {
318 z[i] = this->data[i].integral(var);
319 }
320
321 return z;
322}
323
324
325template <class T>
326VpsMap<T> VpsMap<T>::substitute(const VpsMap<T> &y, int trunc) const {
327 // Consistency checks.
328 this->check();
329 y.check();
330
331 if(this->variables != y.getDimension()) {
332 throw SizeError("VpsMap::substitute()", "Inconsistent dimensions.");
333 }
334
335 // This statement makes sure that the substitution table exists.
336 const Array1D<TpsSubstitution> &table =
337 TpsData::getTpsData(trunc, this->getVariables())->getSubTable();
338
339
340 // Highest order present in this Vps.
341 int topOrder = this->getTopOrder();
342 VpsMap<T> z = constantTerm();
343
344 if(topOrder) {
345 Array1D< Tps<T> > product(topOrder + 1);
346 product[0] = Tps<T>(T(1));
347
348 for(int next = 1; next < table.size();) {
349 const TpsSubstitution &s = table[next];
350 product[s.order] = product[s.order-1].multiply(y[s.variable], trunc);
351
352 for(int v = 0; v < this->getDimension(); v++) {
353 int maxOrder = this->data[v].getMaxOrder();
354 int cutOrder = std::min(this->data[v].getTruncOrder(), trunc);
355
356 if(s.order <= maxOrder && this->data[v][s.index] != T(0)) {
357 z.data[v] += this->data[v][s.index] *
358 product[s.order].truncate(cutOrder);
359 }
360 }
361
362 next = (s.order < topOrder) ? next + 1 : s.skip;
363 }
364 }
365
366 return z;
367}
368
369
370template <class T>
372 if(x.ncols() != this->getDimension()) {
373 throw SizeError("VpsMap::substituteInto()", "Inconsistent dimensions.");
374 }
375
376 VpsMap<T> z(x.nrows(), this->getVariables());
377
378 for(int i = 0; i < this->getDimension(); i++) {
379 for(int j = 0; j < this->data.size(); j++) {
380 z[i] += this->data[j] * x(i, j);
381 }
382 }
383
384 return z;
385}
386
387
388// Implementation of global functions for class Vps<T>.
389// ------------------------------------------------------------------------
390
391template <class T> inline
393{ VpsMap<T> z(x); return z *= y; }
394
395template <class T> inline
397{ VpsMap<T> z(x); return z *= y; }
398
399template <class T> inline
400VpsMap<T> operator*(const VpsMap<T> &x, const T &y)
401{ VpsMap<T> z(x); return z *= y; }
402
403template <class T> inline
404VpsMap<T> operator*(const T &x, const VpsMap<T> &y)
405{ VpsMap<T> z(x); return z *= y; }
406
407template <class T> inline
409{ VpsMap<T> z(x); return z /= y; }
410
411template <class T> inline
412VpsMap<T> operator/(const VpsMap<T> &x, const T &y)
413{ VpsMap<T> z(x); return z /= y; }
414
415
416template <class T> inline
418{ VpsMap<T> z(x); return z += y; }
419
420template <class T> inline
422{ VpsMap<T> z(x); return z -= y; }
423
424template <class T> inline
426{ VpsMap<T> z(x); return z += y; }
427
428template <class T> inline
430{ VpsMap<T> z(x); return z -= y; }
431
432template <class T> inline
434{ VpsMap<T> z(y); return z += x; }
435
436template <class T> inline
438{ VpsMap<T> z(- y); return z += x; }
439
440template <class T> inline
441std::istream &operator>>(std::istream &is, VpsMap<T> &x)
442{ return x.get(is); }
443
444template <class T> inline
445std::ostream &operator<<(std::ostream &os, const VpsMap<T> &x)
446{ return x.put(os); }
447
448#endif // CLASSIC_VpsMap_HH
VpsMap< T > operator+(const VpsMap< T > &x, const VpsMap< T > &y)
Add.
Definition: VpsMap.h:417
VpsMap< T > operator*(const VpsMap< T > &x, const Tps< T > &y)
Multiply.
Definition: VpsMap.h:392
std::ostream & operator<<(std::ostream &, const VpsMap< T > &x)
Insert to stream.
Definition: VpsMap.h:445
VpsMap< T > operator-(const VpsMap< T > &x, const VpsMap< T > &y)
Subtract.
Definition: VpsMap.h:421
VpsMap< T > operator/(const VpsMap< T > &x, const Tps< T > &y)
Divide.
Definition: VpsMap.h:408
std::istream & operator>>(std::istream &, VpsMap< T > &x)
Extract from stream.
Definition: VpsMap.h:441
T::PETE_Expr_t::PETE_Return_t min(const PETE_Expr< T > &expr, NDIndex< D > &loc)
Definition: ReductionLoc.h:76
int size() const
Get array size.
Definition: Array1D.h:228
int nrows() const
Get number of rows.
Definition: Array2D.h:301
int ncols() const
Get number of columns.
Definition: Array2D.h:307
Matrix.
Definition: Matrix.h:39
Truncated power series.
Definition: Tps.h:46
Vector.
Definition: Vector.h:37
Truncate power series map.
Definition: VpsMap.h:43
~VpsMap()
Definition: VpsMap.h:209
VpsMap< T > derivative(int var) const
Derivative with respect to variable [b]var[/b].
Definition: VpsMap.h:302
Vector< T > constantTerm() const
Evaluate map at origin.
Definition: VpsMap.h:266
VpsMap()
Definition: VpsMap.h:167
VpsMap< T > integral(int var) const
Integral with respect to variable [b]var[/b].
Definition: VpsMap.h:314
VpsMap< T > substituteInto(const Matrix< T > &x) const
Substitute.
Definition: VpsMap.h:371
VpsMap< T > & operator=(const VpsMap< T > &y)
Definition: VpsMap.h:214
Matrix< T > linearTerms() const
Extract linear terms at origin.
Definition: VpsMap.h:288
VpsMap< T > substitute(const VpsMap< T > &vv) const
Substitute.
Definition: VpsMap.h:221
static TpsData * getTpsData(int nOrd, int nVar)
Definition: TpsData.cpp:43
const Array1D< TpsSubstitution > & getSubTable() const
Definition: TpsData.h:128
Substitution for Tps<T>.
Vector truncated power series.
Definition: Vps.h:38
std::ostream & put(std::ostream &os) const
Put a Vps<T> to stream os.
Definition: Vps.hpp:238
std::istream & get(std::istream &is)
Get a Vps<T> from stream is.
Definition: Vps.hpp:213
Array1D< Tps< T > > data
Definition: Vps.h:143
int getDimension() const
Get dimension (number of Tps<T> components).
Definition: Vps.hpp:247
Vps< T > & operator=(const Vps< T > &)
Definition: Vps.hpp:84
void check() const
Check consistency.
Definition: Vps.hpp:108
Size error exception.
Definition: SizeError.h:33