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