softsusy is hosted by Hepforge, IPPP Durham
SOFTSUSY  4.1
xpr-base.h
Go to the documentation of this file.
1 
14 #ifndef XPR_BASE_H
15 #define XPR_BASE_H
16 
17 #include <iostream>
18 #include <math.h>
19 
20 /****************** VECTOR **************************/
21 
23 template<typename T, typename E>
24 class Xpr // : public Indexable<T,Xpr> // inheriting doesn't work. Pity.
25 {
26  E e;
27 public:
28  Xpr( const E& e_ ) : e(e_) {}
29  T operator() ( int n ) const { return e(n); }
30  int displayStart() const { return e.displayStart(); }
31  int displayEnd() const { return e.displayEnd(); }
32 };
33 
35 template<typename T, typename Container>
36 class ConstRef {
37  const Container& c;
38 public:
39  ConstRef( const Container& c_) : c(c_) {}
40  T operator() ( int n ) const { return c(n); }
41  int displayStart() const { return c.displayStart(); }
42  int displayEnd() const { return c.displayEnd(); }
43 };
44 
46 template<typename T, typename A, typename B, typename Op>
47 class XprBinOp {
48  A a;
49  B b;
50 public:
51  XprBinOp( const A& a_, const B& b_ )
52  : a(a_), b(b_) {}
53  T operator() ( int n ) const
54  {
55  return Op::apply( a(n), b(n) );
56  }
57  int displayStart() const { return a.displayStart(); }
58  int displayEnd() const { return a.displayEnd(); }
59 };
60 
62 template<typename T, typename B, typename Op>
63 class XprScalOp {
64  T a;
65  B b;
66 public:
67  XprScalOp( T a_, const B& b_ )
68  : a(a_), b(b_) {}
69  T operator() ( int n ) const
70  {
71  return Op::apply( a, b(n) );
72  }
73  int displayStart() const { return b.displayStart(); }
74  int displayEnd() const { return b.displayEnd(); }
75 };
76 
78 template<typename T, typename C>
79 class Indexable {
80 public:
81  explicit Indexable() {}
82  T operator() ( int n ) const
83  {
84  return static_cast<const C*>(this)->operator()(n);
85  }
86 
87  int displayStart() const
88  {
89  return static_cast<const C*>(this)->displayStart();
90  }
91 
92  int displayEnd() const
93  {
94  return static_cast<const C*>(this)->displayEnd();
95  }
96 
97  template<class E>
98  C& assign_from(const Xpr<T,E>& x )
99  {
100  C *me = static_cast<C*>(this);
101  for(int i=me->displayStart(); i <= me->displayEnd(); ++i)
102  me->operator()(i) = x(i);
103  return *me;
104  }
105 
106  template<class E>
107  C copy_from(const Xpr<T,E>& x )
108  {
109  C *meptr = static_cast<C*>(this);
110  C me(meptr->displayStart(), meptr->displayEnd());
111  for(int i=me.displayStart(); i <= me.displayEnd(); ++i)
112  me.operator()(i) = x(i);
113  return me;
114  }
115 
116  template<class E>
117  Indexable& operator+=( const Xpr<T,E>& x )
118  {
119  C *me = static_cast<C*>(this);
120  for ( int i=me->displayStart(); i <= me->displayEnd(); ++i)
121  me->operator() (i) += x(i);
122  return *this;
123  }
124 
125  template<class E>
126  Indexable& operator-=( const Xpr<T,E>& x )
127  {
128  C *me = static_cast<C*>(this);
129  for ( int i=me->displayStart(); i <= me->displayEnd(); ++i)
130  me->operator() (i) -= x(i);
131  return *this;
132  }
133 
134  Indexable& operator*=( T x )
135  {
136  C *me = static_cast<C*>(this);
137  for ( int i=me->displayStart(); i <= me->displayEnd(); ++i)
138  me->operator() (i) *= x;
139  return *this;
140  }
141 
142  // ...
143 };
144 
145 /***************** operations ********************/
146 
148 template<typename T>
149 class OpAdd {
150 public:
151  static inline T apply( const T& a, const T& b )
152  {
153  return a + b;
154  }
155 };
156 
158 template<typename T>
159 class OpSubtract {
160 public:
161  static inline T apply( const T& a, const T& b )
162  {
163  return a - b;
164  }
165 };
166 
168 template<typename T>
169 class OpMultiply {
170 public:
171  static inline T apply( const T& a, const T& b ) {
172  return a * b;
173  }
174 };
175 
177 template<typename T>
178 class OpNegate {
179 public:
180  static inline T apply( const T& a )
181  {
182  return -a ;
183  }
184 };
185 
186 
187 /************** MATRIX ********************/
189 template<typename T, typename E>
190 class MatXpr
191 {
192  E e;
193 public:
194  MatXpr( const E& e_ ) : e(e_) {}
195  T operator() ( int m, int n ) const { return e(m,n); }
196  int displayRows() const { return e.displayRows(); }
197  int displayCols() const { return e.displayCols(); }
198  T trace() const {
199  T retval = T();
200  for (int i=1; i<= displayRows(); ++i)
201  retval += this->operator()(i,i);
202  return retval;
203  }
204 };
205 
206 
208 template<typename T, typename Container>
209 class MatConstRef {
210  const Container& c;
211 public:
212  MatConstRef( const Container& c_) : c(c_) {}
213  T operator() ( int m, int n ) const { return c(m,n); }
214  int displayRows() const { return c.displayRows(); }
215  int displayCols() const { return c.displayCols(); }
216 };
217 
219 template<typename T, typename A, typename B, typename Op>
220 class MatXprBinOp {
221  A a;
222  B b;
223 public:
224  MatXprBinOp( const A& a_, const B& b_ )
225  : a(a_), b(b_) {}
226  T operator() ( int m, int n ) const
227  {
228  return Op::apply( a(m,n), b(m,n) );
229  }
230  int displayRows() const { return a.displayRows(); }
231  int displayCols() const { return a.displayCols(); }
232 };
233 
235 template<typename T, typename A, typename B, typename Op>
237  A a;
238  B b;
239 public:
240  MatXprOuterOp( const A& a_, const B& b_ )
241  : a(a_), b(b_) {}
242  T operator() ( int m, int n ) const
243  {
244  return Op::apply( a(m), b(n) );
245  }
246  int displayRows() const { return a.displayEnd(); }
247  int displayCols() const { return b.displayEnd(); }
248 };
249 
251 template<typename T, typename A, typename B>
253  A a;
254  B b;
255 public:
256  MatXprMatMultOp( const A& a_, const B& b_ )
257  : a(a_), b(b_) {}
258  T operator() ( int m, int n ) const
259  {
260  T value = T();
261  for (int i=1; i<= a.displayCols(); ++i) {
262 // #ifdef ARRAY_BOUNDS_CHECKING
263 // int n2, n3;
264 // frexp(a(m,i), &n2); frexp(b(i,n), &n3);
265 // if (n2+n3>=1024 || n2+n3 <= -1024) {
266 // // throw("Multiply * overload; overflow\n");
267 // if (n2+n3<=-1024) return 1e-290;
268 // else return 1e290;
269 // } else
270 // #endif
271  value += a(m,i) * b(i,n);
272  }
273  return value;
274  }
275  int displayRows() const { return a.displayRows(); }
276  int displayCols() const { return b.displayCols(); }
277 };
278 
280 template<typename T, typename A, typename B>
282  A a;
283  B b;
284 public:
285  XprMatMultOp( const A& a_, const B& b_ )
286  : a(a_), b(b_) {}
287  T operator() ( int n ) const
288  {
289  T value = T();
290  for (int i=1; i<= a.displayCols(); ++i) {
291 #ifdef ARRAY_BOUNDS_CHECKING
292  int n2, n3;
293  frexp(a(n,i), &n2); frexp(b(i), &n3);
294  if (n2+n3>=1024 || n2+n3 <= -1024) {
295  // throw("Multiply * overload; overflow\n");
296  if (n2+n3<=-1024) return 1e-290;
297  else return 1e290;
298  } else
299 #endif
300  value += a(n,i) * b(i);
301  }
302  return value;
303  }
304  int displayStart() const { return 1; }
305  int displayEnd() const { return a.displayRows(); }
306 };
307 
309 template<typename T, typename A, typename Op>
311  A a;
312 public:
313  MatXprUnaryOp( const A& a_)
314  : a(a_) {}
315  T operator() ( int m, int n ) const
316  {
317  return Op::apply( a(m,n) );
318  }
319  int displayRows() const { return a.displayRows(); }
320  int displayCols() const { return a.displayCols(); }
321 };
322 
324 template<typename T, typename B, typename Op>
326  T a;
327  B b;
328 public:
329  MatXprScalOp( T a_, const B& b_ )
330  : a(a_), b(b_) {}
331  T operator() ( int m, int n ) const
332  {
333  return Op::apply( a, b(m,n) );
334  }
335  int displayRows() const { return b.displayRows(); }
336  int displayCols() const { return b.displayCols(); }
337 };
338 
340 template<typename T, typename B, typename Op>
342  T a;
343  B b;
344 public:
345  MatXprScalDiagOp1( T a_, const B& b_ )
346  : a(a_), b(b_) {}
347  T operator() ( int m, int n ) const
348  {
349  return (m == n) ? Op::apply( a, b(m,n) ) : b(m,n);
350  }
351  int displayRows() const { return b.displayRows(); }
352  int displayCols() const { return b.displayCols(); }
353 };
354 
356 template<typename T, typename B, typename Op>
358  T a;
359  B b;
360 public:
361  MatXprScalDiagOp2(const B& b_, T a_ )
362  : a(a_), b(b_) {}
363  T operator() ( int m, int n ) const
364  {
365  return (m == n) ? Op::apply( b(m,n) , a ) : b(m,n);
366  }
367  int displayRows() const { return b.displayRows(); }
368  int displayCols() const { return b.displayCols(); }
369 };
370 
372 template<typename T, typename C>
374 public:
375  explicit MatIndexable() {}
376  T operator() ( int m, int n ) const
377  {
378  return static_cast<const C*>(this)->operator()(m,n);
379  }
380 
381  int displayRows() const
382  {
383  return static_cast<const C*>(this)->displayRows();
384  }
385 
386 
387  int displayCols() const
388  {
389  return static_cast<const C*>(this)->displayCols();
390  }
391 
392  template<class E>
393  C& assign_from(const MatXpr<T,E>& x )
394  {
395  C *me = static_cast<C*>(this);
396  for(int m=1; m <= me->displayRows(); ++m)
397  for (int n=1; n <= me->displayCols(); ++n)
398  me->operator()(m,n) = x(m,n);
399  return *me;
400  }
401 
402  template<class E>
403  C copy_from(const MatXpr<T,E>& x )
404  {
405  C * meptr = static_cast<C*>(this);
406  C me(meptr->displayRows(), meptr->displayCols());
407  for(int m=1; m <= me.displayRows(); ++m)
408  for (int n=1; n <= me.displayCols(); ++n)
409  me.operator()(m,n) = x(m,n);
410  return me;
411  }
412 
413  template<class E>
414  MatIndexable& operator+=( const MatXpr<T,E>& x )
415  {
416  C *me = static_cast<C*>(this);
417  for(int m=1; m <= me->displayRows(); ++m)
418  for (int n=1; n <= me->displayCols(); ++n)
419  me->operator()(m,n) += x(m,n);
420  return *this;
421  }
422 
423 
424  MatIndexable& operator*=(T x)
425  {
426  C *me = static_cast<C*>(this);
427  for(int m=1; m <= me->displayRows(); ++m)
428  for (int n=1; n <= me->displayCols(); ++n)
429  me->operator()(m,n) *= x;
430  return *this;
431  }
432 
433 };
434 
435 
436 
437 
438 #endif
const version of under-object for speed upgrade of linear algebra
Definition: xpr-base.h:36
provides indexing of xpr type objects
Definition: xpr-base.h:79
const version of xpr matrix
Definition: xpr-base.h:209
Indexable form of xpr matrix.
Definition: xpr-base.h:373
binary operator between two matrices
Definition: xpr-base.h:220
xpr Matrix multiplication
Definition: xpr-base.h:252
outer form of xpr matrix operator
Definition: xpr-base.h:236
scalar operator on diagonal parts of xpr matrix
Definition: xpr-base.h:341
scalar operator on diagonal parts of xpr matrix - different place for const
Definition: xpr-base.h:357
scalar operator on xpr matrix
Definition: xpr-base.h:325
unary operator on an xpr matrix
Definition: xpr-base.h:310
matrix form of xpr object
Definition: xpr-base.h:191
Add two xpr objects.
Definition: xpr-base.h:149
multiply two xpr objects
Definition: xpr-base.h:169
take the negative of an xpr object
Definition: xpr-base.h:178
subtract two xpr objects
Definition: xpr-base.h:159
binary operator on two xpr objects
Definition: xpr-base.h:47
another xpr matrix multiplication
Definition: xpr-base.h:281
scalar operator between two xpr object
Definition: xpr-base.h:63
Under-object for speed upgrade of linear algebra.
Definition: xpr-base.h:25
double frexp(const Complex &c, int *i)
Gives exponent of the largest number: either imaginary or real part.
Definition: utils.cpp:12