OPAL (Object Oriented Parallel Accelerator Library)  2.2.0
OPAL
TransportMap.hpp
Go to the documentation of this file.
1 #ifndef CLASSIC_TransportMap_CC
2 #define CLASSIC_TransportMap_CC
3 
4 // ------------------------------------------------------------------------
5 // $RCSfile: TransportMap.cpp,v $
6 // ------------------------------------------------------------------------
7 // $Revision: 1.2 $
8 // ------------------------------------------------------------------------
9 // Copyright: see Copyright.readme
10 // ------------------------------------------------------------------------
11 //
12 // Template Class: TransportMap<T,N>
13 // Second-order map in variables of type T, N variables to N variables.
14 //
15 // ------------------------------------------------------------------------
16 // Class category: FixedAlgebra
17 // ------------------------------------------------------------------------
18 //
19 // $Date: 2000/03/28 13:20:56 $
20 // $Author: mad $
21 //
22 // ------------------------------------------------------------------------
23 
24 #include "FixedAlgebra/FLUMatrix.h"
25 #include "FixedAlgebra/FMatrix.h"
26 #include "FixedAlgebra/TransportFun.h"
27 #include "FixedAlgebra/FVector.h"
28 #include "FixedAlgebra/FVps.h"
29 #include "Utilities/FormatError.h"
30 #include "Utilities/CLRangeError.h"
31 #include <iostream>
32 
33 
34 // Template class TransportMap<T,N>
35 // ------------------------------------------------------------------------
36 
37 template <class T, int N>
38 const int TransportMap<T, N>::SIZE;
39 
40 
41 template <class T, int N>
42 TransportMap<T, N>::TransportMap() {
43  for(int i = 0; i < N; ++i) data[i][i+1] = 1.0;
44 }
45 
46 
47 template <class T, int N>
48 TransportMap<T, N>::TransportMap(const FVps<T, N> &rhs) {
49  for(int i = 0; i < N; ++i) {
50  int size = std::min(rhs[i].getSize(), SIZE);
51  for(int j = 0; j < size; ++j) data[i][j] = rhs[i][j];
52  }
53 }
54 
55 
56 template <class T, int N>
57 TransportMap<T, N>::TransportMap(const FMatrix<T, N, N> &x) {
58  for(int i = 0; i < N; ++i) {
59  for(int j = 0; j < N; j++) data[i][j+1] = x[i][j];
60  }
61 }
62 
63 
64 template <class T, int N>
65 TransportMap<T, N>::TransportMap(const FVector<T, N> &x) {
66  for(int i = 0; i < N; i++) data[i][0] = x[i];
67 }
68 
69 
70 template <class T, int N>
71 const TransportFun<T, N> &TransportMap<T, N>::getComponent(int index) const {
72  if(index > N) {
73  throw CLRangeError("TransportMap::getComponent()", "Index out of range.");
74  }
75 
76  return data[index];
77 }
78 
79 
80 template <class T, int N>
81 void TransportMap<T, N>::setComponent(int index, const TransportFun<T, N> &value) {
82  if(index > N) {
83  throw CLRangeError("TransportMap::setComponent()", "Index out of range.");
84  }
85 
86  data[index] = value;
87 }
88 
89 
90 template <class T, int N> inline
91 TransportFun<T, N> &TransportMap<T, N>::operator[](int index) {
92  return data[index];
93 }
94 
95 template <class T, int N> inline
96 const TransportFun<T, N> &TransportMap<T, N>::operator[](int index) const {
97  return data[index];
98 }
99 
100 
101 template <class T, int N>
102 TransportMap<T, N> TransportMap<T, N>::operator+() const {
103  return *this;
104 }
105 
106 
107 template <class T, int N>
108 TransportMap<T, N> TransportMap<T, N>::operator-() const {
109  TransportMap<T, N> z;
110  for(int i = 0; i < N; i++) z[i] = - data[i];
111  return z;
112 }
113 
114 
115 template <class T, int N>
116 TransportMap<T, N> &TransportMap<T, N>::operator*=(const TransportFun<T, N> &rhs) {
117  for(int i = 0; i < N; i++) data[i] *= rhs;
118  return *this;
119 }
120 
121 
122 template <class T, int N>
123 TransportMap<T, N> &TransportMap<T, N>::operator/=(const TransportFun<T, N> &rhs) {
124  TransportFun<T, N> t = rhs.inverse();
125  for(int i = 0; i < N; i++) data[i] *= t;
126  return *this;
127 }
128 
129 
130 template <class T, int N>
131 TransportMap<T, N> &TransportMap<T, N>::operator*=(const T &rhs) {
132  for(int i = 0; i < N; i++) data[i] *= rhs;
133  return *this;
134 }
135 
136 
137 template <class T, int N>
138 TransportMap<T, N> &TransportMap<T, N>::operator/=(const T &rhs) {
139  for(int i = 0; i < N; i++) data[i] /= rhs;
140  return *this;
141 }
142 
143 
144 template <class T, int N>
145 TransportMap<T, N>& TransportMap<T, N>::operator+=(const TransportMap<T, N> &rhs) {
146  for(int i = 0; i < N; i++) data[i] += rhs[i];
147  return *this;
148 }
149 
150 
151 template <class T, int N>
152 TransportMap<T, N>& TransportMap<T, N>::operator-=(const TransportMap<T, N> &rhs) {
153  for(int i = 0; i < N; i++) data[i] -= rhs[i];
154  return *this;
155 }
156 
157 
158 template <class T, int N>
159 TransportMap<T, N>& TransportMap<T, N>::operator+=(const FVector<T, N> &rhs) {
160  for(int i = 0; i < N; i++) data[i] += rhs[i];
161  return *this;
162 }
163 
164 
165 template <class T, int N>
166 TransportMap<T, N>& TransportMap<T, N>::operator-=(const FVector<T, N> &rhs) {
167  for(int i = 0; i < N; i++) data[i] -= rhs[i];
168  return *this;
169 }
170 
171 
172 template <class T, int N>
173 TransportMap<T, N> TransportMap<T, N>::inverse() const {
174  // Separate out const, linear, and second-order terms.
175  // C := A(0), the constant part of map A,
176  // M := A(1), the linear part of map A,
177  // S := - A(2) = - (A - M - C).
178  FVector<T, N> C;
179  FMatrix<T, N, N> M;
180  TransportMap<T, N> S;
181 
182  for(int i = 0; i < N; ++i) {
183  const TransportFun<T, N> &fun = data[i];
184  C[i] = - fun[0];
185  for(int j = 0; j < N; ++j) {
186  M[i][j] = fun[j+1];
187  }
188  for(int j = N + 1; j < SIZE; ++j) {
189  S[i][j] = - fun[j];
190  }
191  }
192 
193  FLUMatrix<T, N> lu(M);
194  FMatrix<T, N, N> mat(lu.inverse());
195 
196  // Initialize map result with inverse linear map.
197  TransportMap<T, N> result(C);
198 
199  // Compute second-order part of inverse.
200  TransportMap<T, N> L;
201  L -= C;
202  result = mat * (S.substitute(mat * L) + L);
203 
204  return result;
205 }
206 
207 
208 template <class T, int N>
209 void TransportMap<T, N>::identity() {
210  for(int i = 0; i < N; ++i) data[i] = TransportFun<T, N>::makeVariable(i);
211 }
212 
213 
214 template <class T, int N>
215 FVector<T, N> TransportMap<T, N>::constantTerm(const FVector<T, N> &rhs) const {
216  FVector<T, N> z;
217 
218  for(int i = 0; i < N; ++i) {
219  z[i] = data[i].evaluate(rhs);
220  }
221 
222  return z;
223 }
224 
225 
226 template <class T, int N>
227 FVector<T, N> TransportMap<T, N>::constantTerm() const {
228  FVector<T, N> z;
229 
230  for(int i = 0; i < N; i++) {
231  z[i] = data[i][0];
232  }
233 
234  return z;
235 }
236 
237 
238 template <class T, int N>
239 FMatrix<T, N, N> TransportMap<T, N>::linearTerms() const {
240  FMatrix<T, N, N> z;
241 
242  for(int i = 0; i < N; i++) {
243  for(int j = 0; j < N; j++) {
244  z(i, j) = data[i][j+1];
245  }
246  }
247 
248  return z;
249 }
250 
251 
252 template <class T, int N>
253 TransportMap<T, N> TransportMap<T, N>::
254 substitute(const TransportMap<T, N> &rhs) const {
255  TransportMap<T, N> z;
256 
257  for(int n = 0; n < N; ++n) {
258  // Constant terms.
259  z.data[n] = data[n][0];
260  }
261 
262  // Linear terms.
263  for(int i = 0; i < N; ++i) {
264  const TransportFun<T, N> &x = rhs[i];
265  for(int n = 0; n < N; ++n) {
266  z[n] += data[n][i+1] * x;
267  }
268  }
269 
270  // Second-order terms.
271  int k = N;
272  for(int i = 0; i < N; ++i) {
273  for(int j = i; j < N; ++j) {
274  const TransportFun<T, N> x = rhs[i] * rhs[j];
275  ++k;
276  for(int n = 0; n < N; ++n) {
277  z[n] += data[n][k] * x;
278  }
279  }
280  }
281 
282  return z;
283 }
284 
285 
286 template <class T, int N>
287 TransportMap<T, N> TransportMap<T, N>::
288 substitute(const FMatrix<T, N, N> &rhs) const {
289  return substitute(TransportMap<T, N>(rhs));
290 }
291 
292 
293 template <class T, int N>
294 TransportMap<T, N> TransportMap<T, N>::
295 substituteInto(const FMatrix<T, N, N> &lhs) const {
296  TransportMap<T, N> z;
297 
298  for(int i = 0; i < N; ++i) {
299  for(int j = 0; j < N; ++j) {
300  z[i] += lhs[i][j] * data[j];
301  }
302  }
303 
304  return z;
305 }
306 
307 
308 // Global Operators on TransportMap<T,N>
309 // ------------------------------------------------------------------------
310 
311 template <class T, int N>
312 TransportMap<T, N>
313 operator*(const TransportMap<T, N> &lhs, const TransportFun<T, N> &rhs) {
314  TransportMap<T, N> z;
315  for(int i = 0; i < N; ++i) z[i] = lhs[i] * rhs;
316  return z;
317 }
318 
319 
320 template <class T, int N>
321 TransportMap<T, N>
322 operator*(const TransportFun<T, N> &lhs, const TransportMap<T, N> &rhs) {
323  TransportMap<T, N> z;
324  for(int i = 0; i < N; ++i) z[i] = lhs * rhs[i];
325  return z;
326 }
327 
328 
329 template <class T, int N>
330 TransportMap<T, N>
331 operator*(const TransportMap<T, N> &lhs, const T &rhs) {
332  TransportMap<T, N> z;
333  for(int i = 0; i < N; ++i) z[i] = lhs[i] * rhs;
334  return z;
335 }
336 
337 
338 template <class T, int N>
339 TransportMap<T, N> operator*(const T &lhs, const TransportMap<T, N> &rhs) {
340  TransportMap<T, N> z;
341  for(int i = 0; i < N; ++i) z[i] = lhs * rhs[i];
342  return z;
343 }
344 
345 template <class T, int N>
346 TransportMap<T, N> operator*(const FMatrix<T, N, N> &lhs, const TransportMap<T, N> &rhs) {
347  return rhs.substituteInto(lhs);
348 }
349 
350 
351 template <class T, int N>
352 TransportMap<T, N>
353 operator/(const TransportMap<T, N> &lhs, const TransportFun<T, N> &rhs) {
354  TransportMap<T, N> z;
355  for(int i = 0; i < N; ++i) z[i] = lhs[i] / rhs;
356  return z;
357 }
358 
359 
360 template <class T, int N>
361 TransportMap<T, N>
362 operator/(const TransportMap<T, N> &lhs, const T &rhs) {
363  TransportMap<T, N> z;
364  for(int i = 0; i < N; ++i) z[i] = lhs[i] / rhs;
365  return z;
366 }
367 
368 
369 template <class T, int N>
370 TransportMap<T, N>
371 operator+(const TransportMap<T, N> &lhs, const TransportMap<T, N> &rhs) {
372  TransportMap<T, N> z;
373  for(int i = 0; i < N; ++i) z[i] = lhs[i] + rhs[i];
374  return z;
375 }
376 
377 
378 template <class T, int N>
379 TransportMap<T, N>
380 operator-(const TransportMap<T, N> &lhs, const TransportMap<T, N> &rhs) {
381  TransportMap<T, N> z;
382  for(int i = 0; i < N; ++i) z[i] = lhs[i] - rhs[i];
383  return z;
384 }
385 
386 
387 template <class T, int N>
388 TransportMap<T, N>
389 operator+(const TransportMap<T, N> &lhs, const FVector<T, N> &rhs) {
390  TransportMap<T, N> z;
391  for(int i = 0; i < N; ++i) z[i] = lhs[i] + rhs[i];
392  return z;
393 }
394 
395 
396 template <class T, int N>
397 TransportMap<T, N>
398 operator-(const TransportMap<T, N> &lhs, const FVector<T, N> &rhs) {
399  TransportMap<T, N> z;
400  for(int i = 0; i < N; ++i) z[i] = lhs[i] - rhs[i];
401  return z;
402 }
403 
404 
405 template <class T, int N>
406 TransportMap<T, N>
407 operator+(const FVector<T, N> &lhs, const TransportMap<T, N> &rhs) {
408  TransportMap<T, N> z;
409  for(int i = 0; i < N; ++i) z[i] = lhs[i] + rhs[i];
410  return z;
411 }
412 
413 
414 template <class T, int N>
415 TransportMap<T, N>
416 operator-(const FVector<T, N> &lhs, const TransportMap<T, N> &rhs) {
417  TransportMap<T, N> z;
418  for(int i = 0; i < N; ++i) z[i] = lhs[i] - rhs[i];
419  return z;
420 }
421 
422 #endif // CLASSIC_TransportMap_CC
operator+
Matrix< T > operator+(const Matrix< T > &, const Matrix< T > &)
Matrix addition.
Definition: Matrix.h:275
TransportMap::operator*=
TransportMap & operator*=(const TransportFun< T, N > &rhs)
Multiply and assign.
Definition: TransportMap.hpp:116
operator/
Matrix< T > operator/(const Matrix< T > &, const T &)
Matrix divided by scalar.
Definition: Matrix.h:329
T
Definition: rbendmap.h:8
TransportMap
Transport map with values of type [b]T[/b] in [b]N[/b] variables.
Definition: Mapper.h:33
TransportMap::operator+=
TransportMap & operator+=(const TransportMap &rhs)
Add.
Definition: TransportMap.hpp:145
TransportMap::setComponent
void setComponent(int, const TransportFun< T, N > &)
Set component.
Definition: TransportMap.hpp:81
FVector
A templated representation for vectors.
Definition: PartBunchBase.h:26
TransportMap::operator+
TransportMap operator+() const
Unary plus.
Definition: TransportMap.hpp:102
operator*
Matrix< T > operator*(const Matrix< T > &, const Matrix< T > &)
Matrix multiply.
Definition: Matrix.h:297
TransportMap::linearTerms
FMatrix< T, N, N > linearTerms() const
Extract Transport terms at origin.
Definition: TransportMap.hpp:239
TransportMap::getComponent
const TransportFun< T, N > & getComponent(int n) const
Get component.
Definition: TransportMap.hpp:71
TransportMap::operator-=
TransportMap & operator-=(const TransportMap &rhs)
Subtract.
Definition: TransportMap.hpp:152
FLUMatrix
A templated representation of a LU-decomposition.
Definition: FLUMatrix.h:42
TransportFun
Transport function in N variables of type T.
Definition: TransportFun.h:40
TransportMap::inverse
TransportMap inverse() const
Inverse.
Definition: TransportMap.hpp:173
TransportMap::data
TransportFun< T, N > data[N]
Definition: TransportMap.h:144
TransportMap::operator-
TransportMap operator-() const
Unary minus.
Definition: TransportMap.hpp:108
TransportMap::substitute
TransportMap substitute(const FMatrix< T, N, N > &rhs) const
Substitute matrix into map.
Definition: TransportMap.hpp:288
TransportMap::operator/=
TransportMap & operator/=(const TransportFun< T, N > &rhs)
Divide and assign.
Definition: TransportMap.hpp:123
CLRangeError
Range error.
Definition: CLRangeError.h:33
TransportMap::operator[]
TransportFun< T, N > & operator[](int)
Get component.
Definition: TransportMap.hpp:91
TransportMap::TransportMap
TransportMap()
Default constructor.
Definition: TransportMap.hpp:42
TransportMap::substituteInto
TransportMap substituteInto(const FMatrix< T, N, N > &lhs) const
Substitute map into matrix.
Definition: TransportMap.hpp:295
FLUMatrix::inverse
FMatrix< T, N, N > inverse() const
Get inverse.
Definition: FLUMatrix.h:236
operator-
Matrix< T > operator-(const Matrix< T > &, const Matrix< T > &)
Matrix subtraction.
Definition: Matrix.h:282
TransportMap::identity
void identity()
Set to identity.
Definition: TransportMap.hpp:209
TransportMap::constantTerm
FVector< T, N > constantTerm() const
Definition: TransportMap.hpp:227
min
T::PETE_Expr_t::PETE_Return_t min(const PETE_Expr< T > &expr, NDIndex< D > &loc)
Definition: ReductionLoc.h:95
FVps
Vector truncated power series in n variables.
Definition: MPSplitIntegrator.h:32
TransportFun::inverse
TransportFun inverse() const
Approximate reciprocal value 1/(*this).
Definition: TransportFun.hpp:239
Generated on Sat Jan 4 2020 14:47:58 for OPAL (Object Oriented Parallel Accelerator Library) by   doxygen 1.8.5