arrays.h
Go to the documentation of this file.
1 /*
2  @copyright Russell Standish 2000-2013
3  @author Russell Standish
4  This file is part of EcoLab
5 
6  Open source licensed under the MIT license. See LICENSE for details.
7 */
8 
9 
10 
11 
12 
13 
14 
15 
16 
17 
18 
19 
20 
21 
22 
23 
24 
25 
26 
27 
28 
29 
30 
31 
35 #ifndef ARRAYS_H
36 #define ARRAYS_H
37 
38 #include "random.h"
39 #include "classdesc.h"
40 #include <assert.h>
41 #include <cmath>
42 #include <limits>
43 
44 #include <set>
45 #include <vector>
46 #include <iostream>
47 #include <iterator>
48 
49 #ifdef MKL
50 #ifdef MKL_STUB
51 #include "mkl_stub.h"
52 #else
53 #include <mkl.h>
54 #endif
55 #endif
56 
57 #include <pack_base.h>
58 #include <TCL_obj_base.h>
59 
60 namespace ecolab
61 {
62 
63  namespace array_ns
64  {
65 
79  template <class T> class array;
80  template <class E, class Op> class unop;
81  template <class E1, class E2, class Op> class binop;
82  template <class E, class I> class RVindex;
83  template <class E, class I> class LVindex;
84 
85  /*
86  type traits classes for testing whether a template argument is an
87  expression - an array, a unary operator on an expression, or a
88  binary operator with at least one argument being an expression
89  */
90 
91  template <class T> struct is_expression
92  {static const bool value=false;};
93 
94  template <class T> struct is_expression<array<T> >
95  {static const bool value=true;};
96 
97  template <class T, class Op> struct is_expression<unop<T,Op> >
98  {static const bool value=true;};
99 
100  template <class E1, class E2, class Op> struct is_expression<binop<E1,E2,Op> >
101  {static const bool value=true;};
102 
103  /* is true if at least one argument is an expression */
104  template <class E1, class E2> struct one_is_expression
105  {static const bool value=is_expression<E1>::value || is_expression<E2>::value;};
106 
107  /* is true if both arguments are expressions */
108  template <class E1, class E2> struct both_are_expressions
109  {static const bool value=is_expression<E1>::value && is_expression<E2>::value;};
110 
111  /* type traits class indicating if type is simple data type */
112  template <class T> struct is_scalar {static const bool value=false;};
113 
114  template <> struct is_scalar<bool> {static const bool value=true;};
115  template <> struct is_scalar<char> {static const bool value=true;};
116  template <> struct is_scalar<short> {static const bool value=true;};
117  template <> struct is_scalar<int> {static const bool value=true;};
118  template <> struct is_scalar<long> {static const bool value=true;};
119  template <> struct is_scalar<unsigned char> {static const bool value=true;};
120  template <> struct is_scalar<unsigned short> {static const bool value=true;};
121  template <> struct is_scalar<unsigned int> {static const bool value=true;};
122  template <> struct is_scalar<unsigned long> {static const bool value=true;};
123  template <> struct is_scalar<float> {static const bool value=true;};
124  template <> struct is_scalar<double> {static const bool value=true;};
125 
126  /* type traits class indicating if type is integral */
127  template <class T> struct is_integer {static const bool value=false;};
128  template <> struct is_integer<bool> {static const bool value=true;};
129  template <> struct is_integer<char> {static const bool value=true;};
130  template <> struct is_integer<short> {static const bool value=true;};
131  template <> struct is_integer<int> {static const bool value=true;};
132  template <> struct is_integer<long> {static const bool value=true;};
133  template <> struct is_integer<unsigned char> {static const bool value=true;};
134  template <> struct is_integer<unsigned short> {static const bool value=true;};
135  template <> struct is_integer<unsigned int> {static const bool value=true;};
136  template <> struct is_integer<unsigned long> {static const bool value=true;};
137 
138 #ifdef HAVE_LONGLONG
139  template <> struct is_scalar<long long> {static const bool value=true;};
140  template <> struct is_scalar<unsigned long long> {static const bool value=true;};
141  template <> struct is_integer<long long> {static const bool value=true;};
142  template <> struct is_integer<unsigned long long> {static const bool value=true;};
143 #endif
144 
145  /* type traits class indicating if type is floating point */
146  template <class T> struct is_float {static const bool value=false;};
147  template <> struct is_float<float> {static const bool value=true;};
148  template <> struct is_float<double> {static const bool value=true;};
149 
150  template <class E1, class E2> struct is_expression_and_scalar
151  {static const bool value=is_expression<E1>::value && is_scalar<E2>::value;};
152 
153 
154  template <bool, class type=void> struct enable_if_c {typedef type T;};
155  template <class T> struct enable_if_c<false,T> {};
156  /* used for controlled template specialisation: stolen from boost::enable_if */
157  template <class Cond, class T=void> struct enable_if:
158  public enable_if_c<Cond::value,T> {};
159 
160  // to help distinguish functions templates on old compilers (eg gcc 3.2)
161  template <int> struct dummy {dummy(int) {} };
162 
163 
164  /*
165  \c type_traits<T>::value_type return \c T::value_type is \a T is an
166  \em expression, or \c T if \a T is scalar
167  */
168 
169  /* traits template for non-scalar, non-expressions */
170  template <class T, bool=false, bool=false>
171  struct tt
172  {
173  typedef void value_type;
174  };
175 
176  /* traits template for expressions */
177  template <class E>
178  struct tt<E,true,false>
179  {
180  typedef typename E::value_type value_type;
181  };
182 
183  /* traits template for scalars */
184  template <class S>
185  struct tt<S,false,true>
186  {
187  typedef S value_type;
188  };
189 
190  /* combine the two */
191  template <class T>
192  struct traits: public tt<T,is_expression<T>::value,is_scalar<T>::value> {};
193 
194  /*
195  result type of a binary expression T + U
196  */
197  template <class E1, class E2>
198  struct result
199  {
200  typedef typename result<
201  typename traits<E1>::value_type,
202  typename traits<E2>::value_type
203  >::value_type
204  value_type;
205  };
206 
207  template <class T> struct result<void,T> {typedef void value_type;};
208  template <class T> struct result<T,void> {typedef void value_type;};
209  template <> struct result<void,void> {typedef void value_type;};
210 
211  /*
212  result type of integer/integer binary ops
213  */
214 #define RESULTIII(T1,T2,R) \
215  template <> struct result<T1,T2> {typedef R value_type;}; \
216  template <> struct result<T1,unsigned T2> {typedef R value_type;}; \
217  template <> struct result<unsigned T1, T2> {typedef R value_type;}; \
218  template <> struct result<unsigned T1,unsigned T2>{typedef unsigned R value_type;};
219 
220  /*
221  result type of bool/integer binary ops
222  */
223 #define RESULTBI(T1) \
224  template <> struct result<T1,bool> {typedef T1 value_type;}; \
225  template <> struct result<bool,T1> {typedef T1 value_type;}; \
226  template <> struct result<bool,unsigned T1> {typedef unsigned T1 value_type;}; \
227  template <> struct result<unsigned T1, bool> {typedef unsigned T1 value_type;}; \
228 
229  /*
230  result type of bool/float binary ops
231  */
232 #define RESULTBF(T1) \
233  template <> struct result<T1,bool> {typedef T1 value_type;}; \
234  template <> struct result<bool,T1> {typedef T1 value_type;}; \
235 
236  /*
237  result type of integer/float binary ops
238  */
239 #define RESULTIFF(T1,T2,R) \
240  template <> struct result<T1,T2> {typedef R value_type;}; \
241  template <> struct result<unsigned T1, T2> {typedef R value_type;}; \
242 
243  /*
244  result type of float/integer binary ops
245  */
246 #define RESULTFIF(T1,T2,R) \
247  template <> struct result<T1,T2> {typedef R value_type;}; \
248  template <> struct result<T1,unsigned T2> {typedef R value_type;}; \
249 
250  /*
251  result type of float/float binary ops
252  */
253 #define RESULTFFF(T1,T2,R) \
254  template <> struct result<T1,T2> {typedef R value_type;};
255 
256  template <> struct result<bool,bool> {typedef int value_type;};
257 
258  RESULTBI(char);
259  RESULTBI(short);
260  RESULTBI(int);
261  RESULTBI(long);
262 
263  RESULTBF(float);
264  RESULTBF(double);
265 
266  RESULTIII(char,char,char);
267 
268  RESULTIII(short,short,short);
269  RESULTIII(char,short,short);
270  RESULTIII(short,char,short);
271 
272  RESULTIII(int,int,int);
273  RESULTIII(short,int,int);
274  RESULTIII(int,short,int);
275  RESULTIII(char,int,int);
276  RESULTIII(int,char,int);
277 
278  RESULTIII(long,long,long);
279  RESULTIII(int,long,long);
280  RESULTIII(long,int,long);
281  RESULTIII(short,long,long);
282  RESULTIII(long,short,long);
283  RESULTIII(char,long,long);
284  RESULTIII(long,char,long);
285 
286 #ifdef HAVE_LONGLONG
287  RESULTIII(long long,long long,long long);
288  RESULTIII(long,long long,long long);
289  RESULTIII(long long,long,long long);
290  RESULTIII(int,long long,long long);
291  RESULTIII(long long,int,long long);
292  RESULTIII(short,long long,long long);
293  RESULTIII(long long,short,long long);
294  RESULTIII(char,long long,long long);
295  RESULTIII(long long,char,long long);
296 #endif
297 
298  RESULTFFF(float,float,float);
299  RESULTIFF(long,float,float);
300  RESULTFIF(float,long,float);
301  RESULTIFF(int,float,float);
302  RESULTFIF(float,int,float);
303  RESULTIFF(short,float,float);
304  RESULTFIF(float,short,float);
305  RESULTIFF(char,float,float);
306  RESULTFIF(float,char,float);
307 
308 #ifdef _MSC_VER
309  RESULTFIF(float,__int64,float);
310  RESULTFIF(double,__int64,double);
311 #endif
312 
313  RESULTFFF(double,double,double);
314  RESULTFFF(float,double,double);
315  RESULTFFF(double,float,double);
316  RESULTIFF(long,double,double);
317  RESULTFIF(double,long,double);
318  RESULTIFF(int,double,double);
319  RESULTFIF(double,int,double);
320  RESULTIFF(short,double,double);
321  RESULTFIF(double,short,double);
322  RESULTIFF(char,double,double);
323  RESULTFIF(double,char,double);
324 
325 #ifdef HAVE_LONGLONG
326  RESULTIFF(long long,float,float);
327  RESULTFIF(float,long long,float);
328  RESULTIFF(long long,double,double);
329  RESULTFIF(double,long long,double);
330 #endif
331 
334  template <class E> typename
335  E::value_type& at(E& x, size_t i,
336  typename enable_if< is_expression<E>, int >::T dum=0)
337  {return x[i];}
338 
339  template <class E> typename
340  E::value_type at(const E& x, size_t i,
341  typename enable_if< is_expression<E>, int >::T dum=0)
342  {return x[i];}
343 
344  template <class S>
345  S& at(S& x, size_t i,
346  typename enable_if< is_scalar<S>, int >::T dum=0)
347  {return x;}
348 
349  template <class S>
350  S at(S x, size_t i,
351  typename enable_if< is_scalar<S>, int >::T dum=0)
352  {return x;}
354 
357  template <class E> typename
358  enable_if< is_expression<E>, size_t >::T
359  len(const E& x, dummy<0> d=0) {return x.size();}
360 
361  template <class S> typename
362  enable_if< is_scalar<S>, size_t >::T
363  len(const S& x, dummy<1> d=0) {return 1;}
364 
365  /*
366  \c conformance_check(e1,e2) is an assertion that e1 & e2 are conformant
367  (same size), or one is scalar
368  */
369  template <class E1, class E2>
370  void conformance_check(const E1& e1, const E2& e2)
371  {
372  assert(
374  len(e1)==len(e2)
375  );
376  }
377 
378  /* support class for elementwise negation */
379  template <class T>
380  struct Neg
381  {
382  typedef T value_type;
383  T operator()(T x) const {return -x;}
384  };
385 
386  /* support class for elementwise boolean not */
387 
388  template <class T>
389  struct Not
390  {
391  typedef bool value_type;
392  bool operator()(T x) const {return !x;}
393  };
394 
395  /* support class for elementwise + */
396 
397  template <class E1, class E2>
398  struct Add
399  {
400  typedef typename traits<E1>::value_type E1v;
401  typedef typename traits<E2>::value_type E2v;
402  typedef typename result<E1v,E2v>::value_type value_type;
403  value_type operator()(E1v x, E2v y) const {return x+y;}
404  };
405 
406  /* support class for elementwise - */
407 
408  template <class E1, class E2>
409  struct Sub
410  {
411  typedef typename traits<E1>::value_type E1v;
412  typedef typename traits<E2>::value_type E2v;
413  typedef typename result<E1v,E2v>::value_type value_type;
414  value_type operator()(E1v x, E2v y) const {return x-y;}
415  };
416 
417  /* support class for elementwise * */
418 
419  template <class E1, class E2>
420  struct Mul
421  {
422  typedef typename traits<E1>::value_type E1v;
423  typedef typename traits<E2>::value_type E2v;
424  typedef typename result<E1v,E2v>::value_type value_type;
425  value_type operator()(E1v x, E2v y) const {return x*y;}
426  };
427 
428 
429  /* support class for elementwise / */
430  template <class E1, class E2>
431  struct Div
432  {
433  typedef typename traits<E1>::value_type E1v;
434  typedef typename traits<E2>::value_type E2v;
435  typedef typename result<E1v,E2v>::value_type value_type;
436  value_type operator()(E1v x, E2v y) const {return x/y;}
437  };
438 
439  /* support class for elementwise % */
440  template <class E1, class E2>
441  struct Mod
442  {
443  typedef typename traits<E1>::value_type E1v;
444  typedef typename traits<E2>::value_type E2v;
445  typedef typename result<E1v,E2v>::value_type value_type;
446  value_type operator()(E1v x, E2v y) const {return x%y;}
447  };
448 
449  /* support class for elementwise == */
450  template <class E1, class E2>
451  struct Eq
452  {
453  typedef typename traits<E1>::value_type E1v;
454  typedef typename traits<E2>::value_type E2v;
455  typedef bool value_type;
456  bool operator()(E1v x, E2v y) const {return x==y;}
457  };
458 
459  /* support class for elementwise != */
460  template <class E1, class E2>
461  struct Neq
462  {
463  typedef typename traits<E1>::value_type E1v;
464  typedef typename traits<E2>::value_type E2v;
465  typedef bool value_type;
466  bool operator()(E1v x, E2v y) const {return x!=y;}
467  };
468 
469  /* support class for elementwise > */
470  template <class E1, class E2>
471  struct Gt
472  {
473  typedef typename traits<E1>::value_type E1v;
474  typedef typename traits<E2>::value_type E2v;
475  typedef bool value_type;
476  bool operator()(E1v x, E2v y) const {return x>y;}
477  };
478 
479  /* support class for elementwise >= */
480  template <class E1, class E2>
481  struct Gte
482  {
483  typedef typename traits<E1>::value_type E1v;
484  typedef typename traits<E2>::value_type E2v;
485  typedef bool value_type;
486  bool operator()(E1v x, E2v y) const {return x>=y;}
487  };
488 
489  /* support class for elementwise < */
490  template <class E1, class E2>
491  struct Lt
492  {
493  typedef typename traits<E1>::value_type E1v;
494  typedef typename traits<E2>::value_type E2v;
495  typedef bool value_type;
496  bool operator()(E1v x, E2v y) const {return x<y;}
497  };
498 
499  /* support class for elementwise <= */
500  template <class E1, class E2>
501  struct Lte
502  {
503  typedef typename traits<E1>::value_type E1v;
504  typedef typename traits<E2>::value_type E2v;
505  typedef bool value_type;
506  bool operator()(E1v x, E2v y) const {return x<=y;}
507  };
508 
509  /* support class for elementwise && */
510  template <class E1, class E2>
511  struct And
512  {
513  typedef typename traits<E1>::value_type E1v;
514  typedef typename traits<E2>::value_type E2v;
515  typedef bool value_type;
516  bool operator()(E1v x, E2v y) const {return x&&y;}
517  };
518 
519  /* support class for elementwise || */
520  template <class E1, class E2>
521  struct Or
522  {
523  typedef typename traits<E1>::value_type E1v;
524  typedef typename traits<E2>::value_type E2v;
525  typedef bool value_type;
526  bool operator()(E1v x, E2v y) const {return x||y;}
527  };
528 
529  template <class E>
531  operator-(const E& e) {return unop<E,Neg<typename E::value_type> >(e);}
532 
533  template <class E>
535  operator!(const E& e) {return unop<E,Not<typename E::value_type> >(e);}
536 
537  template <class E1, class E2>
539  operator+(const E1& e1, const E2& e2)
540  {
541  return binop<E1,E2,Add<E1,E2> >(e1,e2);
542  }
543 
544  template <class E1, class E2>
545  typename enable_if< one_is_expression<E1,E2>, binop<E1,E2,Sub<E1,E2> > >::T
546  operator-(const E1& e1, const E2& e2)
547  {
548  return binop<E1,E2,Sub<E1,E2> >(e1,e2);
549  }
550 
551  template <class E1, class E2>
552  typename enable_if< one_is_expression<E1,E2>, binop<E1,E2,Mul<E1,E2> > >::T
553  operator*(const E1& e1, const E2& e2)
554  {
555  return binop<E1,E2,Mul<E1,E2> >(e1,e2);
556  }
557 
558 #ifdef MKL
559  array<float> operator/(float x, const array<float>& y);
560  array<float> operator/(double x, const array<float>& y);
561  array<double> operator/(double x, const array<double>& y);
562  array<float> operator/(const array<float>& x, const array<float>& y);
563  array<double> operator/(const array<double>& x, const array<double>& y);
564 
565  template <class E>
566  typename enable_if<is_expression<E>, array<typename E::value_type> >::T
567  operator/(typename E::value_type x, const E& y)
568  {
569  array<typename E::value_type> Y(y.size()); Y=y;
570  return operator/(x,Y);
571  }
572 
573  template <class E1, class E2>
576  operator/(const E1& x, const E2& y)
577  {
578  array<typename E1::value_type> X(x.size()); X=x;
579  array<typename E2::value_type> Y(y.size()); Y=y;
580  return operator/(X,Y);
581  }
582 
583  template <class E1, class E2>
585  operator/(const E1& e1, const E2& e2)
586  {
587  return binop<E1,E2,Div<E1,E2> >(e1,e2);
588  }
589 
590 #else
591 
592  template <class E1, class E2>
594  operator/(const E1& e1, const E2& e2)
595  {
596  return binop<E1,E2,Div<E1,E2> >(e1,e2);
597  }
598 
599 #endif
600 
601  template <class E1, class E2>
602  typename enable_if< one_is_expression<E1,E2>, binop<E1,E2,Mod<E1,E2> > >::T
603  operator%(const E1& e1, const E2& e2)
604  {
605  return binop<E1,E2,Mod<E1,E2> >(e1,e2);
606  }
607 
608  namespace both_are_expressions_ns
609  {
610  /* concatenation */
611  template <class E1, class E2>
612  typename enable_if<
614  array<
615  typename result<
616  typename E1::value_type,
617  typename E2::value_type
618  >::value_type
619  >
620  >::T
621  operator<<(const E1& e1, const E2& e2)
622  {
623  array<
624  typename result<
625  typename E1::value_type,
626  typename E2::value_type
627  >::value_type
628  > ret(e1.size()+e2.size());
629  for (size_t i=0; i<e1.size(); i++)
630  ret[i]=e1[i];
631  for (size_t i=0; i<e2.size(); i++)
632  ret[i+e1.size()]=e2[i];
633  return ret;
634  }
635  }
636 
637  namespace expression_scalar_ns
638  {
639  /* concatenation */
640  template <class E1, class E2>
641  typename enable_if<
643  typename enable_if< is_scalar<E2>,
644  array<
645  typename result<typename E1::value_type,E2>::value_type
646  >
647  >::T
648  >::T
649  operator<<(const E1& e1, const E2& e2)
650  {
651  array<
652  typename result<typename E1::value_type,E2>::value_type
653  > ret(e1.size()+1);
654  for (size_t i=0; i<e1.size(); i++)
655  ret[i]=e1[i];
656  ret[e1.size()]=e2;
657  return ret;
658  }
659  }
660 
661  namespace scalar_expression_ns
662  {
663  /* concatenation */
664  template <class E1, class E2>
665  typename enable_if<
666  is_scalar<E1>,
669  >::T
670  >::T
671  operator<<(const E1& e1, const E2& e2)
672  {
674  ret(e2.size()+1);
675  ret[0]=e1;
676  for (size_t i=0; i<e2.size(); i++)
677  ret[i+1]=e2[i];
678  return ret;
679  }
680  }
681 
682  using both_are_expressions_ns::operator<<;
683  using expression_scalar_ns::operator<<;
684  using scalar_expression_ns::operator<<;
685 
686  template <class E1, class E2>
687  typename enable_if< one_is_expression<E1,E2>, binop<E1,E2,Eq<E1,E2> > >::T
688  operator==(const E1& e1, const E2& e2)
689  {
690  return binop<E1,E2,Eq<E1,E2> >(e1,e2);
691  }
692 
693  template <class E1, class E2>
694  typename enable_if< one_is_expression<E1,E2>, binop<E1,E2,Neq<E1,E2> > >::T
695  operator!=(const E1& e1, const E2& e2)
696  {
697  return binop<E1,E2,Neq<E1,E2> >(e1,e2);
698  }
699 
700  template <class E1, class E2>
701  typename enable_if< one_is_expression<E1,E2>, binop<E1,E2,Gt<E1,E2> > >::T
702  operator>(const E1& e1, const E2& e2)
703  {
704  return binop<E1,E2,Gt<E1,E2> >(e1,e2);
705  }
706 
707  template <class E1, class E2>
708  typename enable_if< one_is_expression<E1,E2>, binop<E1,E2,Gte<E1,E2> > >::T
709  operator>=(const E1& e1, const E2& e2)
710  {
711  return binop<E1,E2,Gte<E1,E2> >(e1,e2);
712  }
713 
714  template <class E1, class E2>
715  typename enable_if< one_is_expression<E1,E2>, binop<E1,E2,Lt<E1,E2> > >::T
716  operator<(const E1& e1, const E2& e2)
717  {
718  return binop<E1,E2,Lt<E1,E2> >(e1,e2);
719  }
720 
721  template <class E1, class E2>
722  typename enable_if< one_is_expression<E1,E2>, binop<E1,E2,Lte<E1,E2> > >::T
723  operator<=(const E1& e1, const E2& e2)
724  {
725  return binop<E1,E2,Lte<E1,E2> >(e1,e2);
726  }
727 
728  template <class E1, class E2>
729  typename enable_if< one_is_expression<E1,E2>, binop<E1,E2,And<E1,E2> > >::T
730  operator&&(const E1& e1, const E2& e2)
731  {
732  return binop<E1,E2,And<E1,E2> >(e1,e2);
733  }
734 
735  template <class E1, class E2>
736  typename enable_if< one_is_expression<E1,E2>, binop<E1,E2,Or<E1,E2> > >::T
737  operator||(const E1& e1, const E2& e2)
738  {
739  return binop<E1,E2,Or<E1,E2> >(e1,e2);
740  }
741 
746  template <class E1, class E2>
747  bool eq(const E1& e1, const E2& e2)
748  {
749  if(
750  array_ns::len(e1)!=1 && array_ns::len(e2)!=1 &&
751  array_ns::len(e1)!=array_ns::len(e2)
752  )
753  return false;
754  for (size_t i=0; i<std::max(array_ns::len(e1), array_ns::len(e2)); i++)
755  if(array_ns::at(e1,i)!=array_ns::at(e2,i))
756  return false;
757  return true;
758  }
759 
760  /* Functions */
761 #ifndef MKL
762 
763  /* support class for \c pow function */
764  template <class E1, class E2>
765  struct Pow
766  {
767  typedef typename traits<E1>::value_type E1v;
768  typedef typename traits<E2>::value_type E2v;
769  typedef typename result<E1v,E2v>::value_type value_type;
770  value_type operator()(E1v x, E2v y) const {
771  return std::pow(static_cast<value_type>(x),static_cast<value_type>(y));
772  }
773  };
774 
776  template <class E1, class E2>
777  binop<E1,E2,Pow<E1,E2> > pow(const E1& e1, const E2& e2)
778  {
779  return binop<E1,E2,Pow<E1,E2> >(e1,e2);
780  }
781 
782  /* support class for \c exp function */
783  template <class T>
784  struct Exp
785  {
786  typedef T value_type;
787  T operator()(T x) const {return std::exp(x);}
788  };
789 
791  template <class E>
794 
795  /* support class for \c log function */
796  template <class T>
797  struct Log
798  {
799  typedef T value_type;
800  T operator()(T x) const {return std::log(x);}
801  };
802 
804  template <class E>
807 
808  /* support class for \c log10 function */
809  template <class T>
810  struct Log10
811  {
812  typedef T value_type;
813  T operator()(T x) const {return std::log10(x);}
814  };
815 
817  template <class E>
820 
821  /* support class for \c sin function */
822  template <class T>
823  struct Sin
824  {
825  typedef T value_type;
826  T operator()(T x) const {return std::sin(x);}
827  };
828 
830  template <class E>
833 
834  /* support class for \c cos function */
835  template <class T>
836  struct Cos
837  {
838  typedef T value_type;
839  T operator()(T x) const {return std::cos(x);}
840  };
841 
843  template <class E>
846 
847  /* support class for \c tan function */
848  template <class T>
849  struct Tan
850  {
851  typedef T value_type;
852  T operator()(T x) const {return std::tan(x);}
853  };
854 
856  template <class E>
859 
860  /* support class for \c asin function */
861  template <class T>
862  struct Asin
863  {
864  typedef T value_type;
865  T operator()(T x) const {return std::asin(x);}
866  };
867 
869  template <class E>
872 
873  /* support class for \c acos function */
874  template <class T>
875  struct Acos
876  {
877  typedef T value_type;
878  T operator()(T x) const {return std::acos(x);}
879  };
880 
882  template <class E>
885 
886  /* support class for \c atan function */
887  template <class T>
888  struct Atan
889  {
890  typedef T value_type;
891  T operator()(T x) const {return std::atan(x);}
892  };
893 
895  template <class E>
898 
899  /* support class for \c atan2 function */
900  template <class E1, class E2>
901  struct Atan2
902  {
903  typedef typename traits<E1>::value_type E1v;
904  typedef typename traits<E2>::value_type E2v;
905  typedef typename result<E1v,E2v>::value_type value_type;
906  value_type operator()(E1v x, E2v y) const {return std::atan2(x,y);}
907  };
908 
910  template <class E1, class E2>
911  binop<E1,E2,Atan2<E1,E2> > atan2(const E1& e1, const E2& e2)
912  {
913  return binop<E1,E2,Atan2<E1,E2> >(e1,e2);
914  }
915 
916  /* support class for \c sinh function */
917  template <class T>
918  struct Sinh
919  {
920  typedef T value_type;
921  T operator()(T x) const {return std::sinh(x);}
922  };
923 
925  template <class E>
928 
929  /* support class for \c cosh function */
930  template <class T>
931  struct Cosh
932  {
933  typedef T value_type;
934  T operator()(T x) const {return std::cosh(x);}
935  };
936 
938  template <class E>
941 
942  /* support class for \c tanh function */
943  template <class T>
944  struct Tanh
945  {
946  typedef T value_type;
947  T operator()(T x) const {return std::tanh(x);}
948  };
949 
951  template <class E>
954 
955  /* support class for \c sqrt function */
956  template <class T>
957  struct Sqrt
958  {
959  typedef T value_type;
960  T operator()(T x) const {return std::sqrt(x);}
961  };
962 
964  template <class E>
967 #endif // ! MKL
968 
969  /* support class for \c abs function */
970  template <class T>
971  struct Abs
972  {
973  typedef T value_type;
974  T operator()(T x) const {return std::abs(x);}
975  };
976 
977  /* support class for \c SGN function */
978  template <class T>
979  struct Sgn
980  {
981  typedef T value_type;
982  T operator()(T x) const {return x>=0? 1: -1;}
983  };
984 
986  template <class E>
989 
991  template <class E>
994 
1000  template <class E>
1002  {return unop<E,Sgn<typename E::value_type> >(e);}
1003 
1004  /* support class for \c ceil function */
1005  template <class T>
1006  struct Ceil
1007  {
1008  typedef T value_type;
1009  T operator()(T x) const {return std::ceil(x);}
1010  };
1011 
1013  template <class E>
1016 
1017  /* support class for \c floor function */
1018  template <class T>
1019  struct Floor
1020  {
1021  typedef T value_type;
1022  T operator()(T x) const {return std::floor(x);}
1023  };
1024 
1026  template <class E>
1029 
1030  /* support class for \c ldexp function */
1031  template <class E1, class E2>
1032  struct Ldexp
1033  {
1034  typedef typename traits<E1>::value_type E1v;
1035  typedef typename traits<E2>::value_type E2v;
1036  typedef typename result<E1v,E2v>::value_type value_type;
1037  value_type operator()(E1v x, E2v y) const {return std::ldexp(x,y);}
1038  };
1039 
1041  template <class E1, class E2>
1042  binop<E1,E2,Ldexp<E1,E2> > ldexp(const E1& e1, const E2& e2)
1043  {
1044  return binop<E1,E2,Ldexp<E1,E2> >(e1,e2);
1045  }
1046 
1047  /* support class for \c fmod function */
1048  template <class E1, class E2>
1049  struct Fmod
1050  {
1051  typedef typename traits<E1>::value_type E1v;
1052  typedef typename traits<E2>::value_type E2v;
1053  typedef typename result<E1v,E2v>::value_type value_type;
1054  value_type operator()(E1v x, E2v y) const {return std::fmod(x,y);}
1055  };
1056 
1058  template <class E1, class E2>
1059  binop<E1,E2,Fmod<E1,E2> > fmod(const E1& e1, const E2& e2)
1060  {
1061  return binop<E1,E2,Fmod<E1,E2> >(e1,e2);
1062  }
1063 
1064  /* support class for \c max function */
1065  template <class E1, class E2>
1066  struct Max
1067  {
1068  typedef typename traits<E1>::value_type E1v;
1069  typedef typename traits<E2>::value_type E2v;
1070  typedef typename result<E1v,E2v>::value_type value_type;
1071  value_type operator()(E1v x, E2v y) const {
1072  return std::max(static_cast<value_type>(x),static_cast<value_type>(y));
1073  }
1074  };
1075 
1077  template <class E1, class E2>
1078  typename enable_if< one_is_expression<E1,E2>, binop<E1,E2,Max<E1,E2> > >::T
1079  max(const E1& e1, const E2& e2)
1080  {
1081  return binop<E1,E2,Max<E1,E2> >(e1,e2);
1082  }
1083 
1084  /* support class for \c min function */
1085  template <class E1, class E2>
1086  struct Min
1087  {
1088  typedef typename traits<E1>::value_type E1v;
1089  typedef typename traits<E2>::value_type E2v;
1090  typedef typename result<E1v,E2v>::value_type value_type;
1091  value_type operator()(E1v x, E2v y) const {
1092  return std::min(static_cast<value_type>(x),static_cast<value_type>(y));
1093  }
1094  };
1095 
1097  template <class E1, class E2>
1098  typename enable_if< one_is_expression<E1,E2>, binop<E1,E2,Min<E1,E2> > >::T
1099  min(const E1& e1, const E2& e2)
1100  {
1101  return binop<E1,E2,Min<E1,E2> >(e1,e2);
1102  }
1103 
1104 
1105 
1106  /* support class for unassignable vector index expressions */
1107  template <class E, class I>
1108  class RVindex
1109  {
1110  const E& expr;
1111  const I& idx;
1112  void operator=(const RVindex&);
1113  public:
1114  typedef typename E::value_type value_type;
1115  size_t size() const {return idx.size();}
1116  value_type operator[](size_t i) const {return expr[idx[i]];}
1117  RVindex(const E& e,const I& i): expr(e), idx(i) {}
1118  };
1119 
1120  /* support class for assignable vector index expressions (lvalues) */
1121  template <class E, class I>
1122  class LVindex
1123  {
1124  E& expr;
1125  const I& idx;
1126  public:
1127  typedef typename E::value_type value_type;
1128  size_t size() const {return idx.size();}
1129  value_type& operator[](size_t i) {return expr[idx[i]];}
1130  value_type operator[](size_t i) const {return expr[idx[i]];}
1131  LVindex(E& e,const I& i): expr(e), idx(i) {}
1132  template <class E1> typename
1134  operator=(const E1& x) {
1135  conformance_check(idx,x);
1136  for (size_t i=0; i<size(); i++) expr[idx[i]]=x[i];
1137  return RVindex<E,I>(expr,idx);
1138  }
1139  template <class E1> typename
1141  operator+=(const E1& x) {
1142  conformance_check(idx,x);
1143  for (size_t i=0; i<size(); i++)
1144  expr[idx[i]]+=x[i];
1145  // return RVindex<E,I>(expr,idx);
1146  return *this;
1147  }
1148  template <class E1> typename
1149  enable_if <is_expression<E1>, RVindex<E,I> >::T
1150  operator*=(const E1& x) {
1151  conformance_check(idx,x);
1152  for (size_t i=0; i<size(); i++) expr[idx[i]]*=x[i];
1153  return RVindex<E,I>(expr,idx);
1154  }
1155  template <class E1> typename
1156  enable_if <is_expression<E1>, RVindex<E,I> >::T
1157  operator&=(const E1& x) {
1158  conformance_check(idx,x);
1159  for (size_t i=0; i<size(); i++) expr[idx[i]]&=x[i];
1160  return RVindex<E,I>(expr,idx);
1161  }
1162  template <class E1> typename
1163  enable_if <is_expression<E1>, RVindex<E,I> >::T
1164  operator|=(const E1& x) {
1165  conformance_check(idx,x);
1166  for (size_t i=0; i<size(); i++) expr[idx[i]]|=x[i];
1167  return RVindex<E,I>(expr,idx);
1168  }
1169  value_type operator=(value_type x) {
1170  for (size_t i=0; i<size(); i++) expr[idx[i]]=x;
1171  return x;
1172  }
1173  value_type operator+=(value_type x) {
1174  for (size_t i=0; i<size(); i++) expr[idx[i]]+=x;
1175  return x;
1176  }
1177  value_type operator-=(value_type x) {
1178  for (size_t i=0; i<size(); i++) expr[idx[i]]-=x;
1179  return x;
1180  }
1181  value_type operator*=(value_type x) {
1182  for (size_t i=0; i<size(); i++) expr[idx[i]]*=x;
1183  return x;
1184  }
1185  value_type operator/=(value_type x) {
1186  for (size_t i=0; i<size(); i++) expr[idx[i]]/=x;
1187  return x;
1188  }
1189  value_type operator%=(value_type x) {
1190  for (size_t i=0; i<size(); i++) expr[idx[i]]%=x;
1191  return x;
1192  }
1193  value_type operator&=(value_type x) {
1194  for (size_t i=0; i<size(); i++) expr[idx[i]]&=x;
1195  return x;
1196  }
1197  value_type operator|=(value_type x) {
1198  for (size_t i=0; i<size(); i++) expr[idx[i]]|=x;
1199  return x;
1200  }
1201  };
1202 
1203  // Vector indices are also expressions
1204  template <class E, class I> struct is_expression<RVindex<E,I> >
1205  {static const bool value=true;};
1206  // Vector indices are also expressions
1207  template <class E, class I> struct is_expression<LVindex<E,I> >
1208  {static const bool value=true;};
1209 
1210  /* unary operator adaptor */
1211  template <class E, class Op>
1212  class unop
1213  {
1214  public:
1215  const E& e;
1216  Op op;
1217 
1218  unop(const E& expr): e(expr) {}
1220  typedef typename E::value_type value_type;
1221  size_t size() const {return e.size();}
1222 
1223  value_type operator[](size_t i) const {return op(e[i]);}
1224 
1226  template <class I> typename
1228  operator[](const I& i) const {return RVindex<unop,I>(*this,i);}
1229 
1230  };
1231 
1232  /* binary operator adaptor */
1233  template <class E1, class E2, class Op>
1234  class binop
1235  {
1236  public:
1237  const E1& e1;
1238  const E2& e2;
1239  Op op;
1240 
1241  binop(const E1& ex1, const E2& ex2): e1(ex1), e2(ex2) {conformance_check(e1,e2);}
1243  typedef typename Op::value_type value_type;
1244  size_t size() const
1245  {return (array_ns::len(e1)==1)? array_ns::len(e2): array_ns::len(e1);}
1246  typename Op::value_type operator[](size_t i) const {
1247  /* use at function here, as E1 or E2 may be scalar */
1248  return op(array_ns::at(e1,i),array_ns::at(e2,i));
1249  }
1250 
1251  // vector indexing
1252  template <class I> typename
1254  operator[](const I& i) const {return RVindex<binop,I>(*this,i);}
1255 
1256  };
1257 
1258  /*
1259  specialisation that converts divisions into multiplications by a reciprocal
1260  */
1261  template <class E1>
1262  class binop<E1,double,Div<E1,double> >
1263  {
1264  public:
1265  const E1& e1;
1266  double e2;
1267 
1268  binop(const E1& ex1, double ex2): e1(ex1), e2(1/ex2) {conformance_check(e1,e2);}
1269  typedef typename result<typename E1::value_type,double>::value_type value_type;
1270  size_t size() const
1271  {return (array_ns::len(e1)==1)? array_ns::len(e2): array_ns::len(e1);}
1272  /* use at function here, as E1 or E2 may be scalar */
1273  value_type operator[](size_t i) const {return array_ns::at(e1,i)*e2;}
1274 
1276  template <class I> typename
1278  operator[](const I& i) const {return RVindex<binop,I>(*this,i);}
1279  };
1280 
1281  /*
1282  specialisation that converts divisions into multiplications by a reciprocal
1283  */
1284  template <class E1>
1285  class binop<E1,float,Div<E1,float> >
1286  {
1287  public:
1288  const E1& e1;
1289  float e2;
1290 
1291  binop(const E1& ex1, float ex2): e1(ex1), e2(1/ex2) {conformance_check(e1,e2);}
1292  typedef typename result<typename E1::value_type,float>::value_type value_type;
1293  size_t size() const
1294  {return (array_ns::len(e1)==1)? array_ns::len(e2): array_ns::len(e1);}
1295  /* use at function here, as E1 or E2 may be scalar */
1296  value_type operator[](size_t i) const {return array_ns::at(e1,i)*e2;}
1297 
1298  // vector indexing
1299  template <class I> typename
1301  operator[](const I& i) const {return RVindex<binop,I>(*this,i);}
1302  };
1303 
1304  /* internal support */
1305  template <class T>
1306  void asg(T* dt, size_t sz, T x) {
1307 #ifdef __ICC
1308 #pragma loop count(1000)
1309 #pragma ivdep
1310 #pragma vector aligned
1311 #endif
1312  for (size_t i=0; i<sz; i++) dt[i]=x;
1313  }
1314 
1315  template <class T>
1316  void asg_plus(T* dt, size_t sz, T x) {
1317 #ifdef __ICC
1318 #pragma loop count(1000)
1319 #pragma ivdep
1320 #pragma vector aligned
1321 #endif
1322  for (size_t i=0; i<sz; i++) dt[i]+=x;
1323  }
1324 
1325  /* internal support */
1326  template <class T>
1327  void asg_minus(T* dt, size_t sz, T x) {
1328 #ifdef __ICC
1329 #pragma loop count(1000)
1330 #pragma ivdep
1331 #pragma vector aligned
1332 #endif
1333  for (size_t i=0; i<sz; i++) dt[i]-=x;
1334  }
1335  /* internal support */
1336  template <class T>
1337  void asg_mul(T* dt, size_t sz, T x) {
1338 #ifdef __ICC
1339 #pragma loop count(1000)
1340 #pragma ivdep
1341 #pragma vector aligned
1342 #endif
1343  for (size_t i=0; i<sz; i++) dt[i]*=x;
1344  }
1345 
1346  /* internal support */
1347  template <class T>
1348  void asg_div(T* dt, size_t sz, T x) {
1349 #ifdef __ICC
1350 #pragma loop count(1000)
1351 #pragma ivdep
1352 #pragma vector aligned
1353 #endif
1354  for (size_t i=0; i<sz; i++) dt[i]/=x;
1355  }
1356 
1357 
1358  /* internal support */
1359  template <class T> typename
1360  enable_if< is_integer<T>, void >::T
1361  asg_mod(T* dt, size_t sz, T x) {
1362 #ifdef __ICC
1363 #pragma loop count(1000)
1364 #pragma ivdep
1365 #pragma vector aligned
1366 #endif
1367  for (size_t i=0; i<sz; i++) dt[i]%=x;
1368  }
1369 
1370 
1371  /* internal support */
1372  template <class T> typename
1373  enable_if< is_integer<T>, void >::T
1374  asg_and(T* dt, size_t sz, T x) {
1375 #ifdef __ICC
1376 #pragma loop count(1000)
1377 #pragma ivdep
1378 #pragma vector aligned
1379 #endif
1380  for (size_t i=0; i<sz; i++) dt[i]&=x;
1381  }
1382 
1383 
1384  /* internal support */
1385  template <class T> typename
1386  enable_if< is_integer<T>, void >::T
1387  asg_or(T* dt, size_t sz, T x) {
1388 #ifdef __ICC
1389 #pragma loop count(1000)
1390 #pragma ivdep
1391 #pragma vector aligned
1392 #endif
1393  for (size_t i=0; i<sz; i++) dt[i]|=x;
1394  }
1395 
1396 
1397  /* internal support */
1398  template <class T, class U>
1399  void asg_v(T * dt, size_t sz, const U& x) {
1400 #ifdef __ICC
1401 #pragma loop count(1000)
1402 #pragma ivdep
1403 #pragma vector aligned
1404 #endif
1405  for (size_t i=0; i<sz; i++) dt[i]=x[i];
1406  }
1407 
1408  /* internal support */
1409  template <class T, class U>
1410  void asg_plus_v(T * dt, size_t sz, const U& x) {
1411 #ifdef __ICC
1412 #pragma loop count(1000)
1413 #pragma ivdep
1414 #pragma vector aligned
1415 #endif
1416  for (size_t i=0; i<sz; i++) dt[i]+=x[i];
1417  }
1418 
1419  /* internal support */
1420  template <class T, class U>
1421  void asg_minus_v(T * dt, size_t sz, const U& x) {
1422 #ifdef __ICC
1423 #pragma loop count(1000)
1424 #pragma ivdep
1425 #pragma vector aligned
1426 #endif
1427  for (size_t i=0; i<sz; i++) dt[i]-=x[i];
1428  }
1429 
1430  /* internal support */
1431  template <class T, class U>
1432  void asg_mul_v(T * dt, size_t sz, const U& x) {
1433 #ifdef __ICC
1434 #pragma loop count(1000)
1435 #pragma ivdep
1436 #pragma vector aligned
1437 #endif
1438  for (size_t i=0; i<sz; i++) dt[i]*=x[i];
1439  }
1440 
1441  /* internal support */
1442  template <class T, class U>
1443  void asg_div_v(T * dt, size_t sz, const U& x) {
1444 #ifdef __ICC
1445 #pragma loop count(1000)
1446 #pragma ivdep
1447 #pragma vector aligned
1448 #endif
1449  for (size_t i=0; i<sz; i++) dt[i]/=x[i];
1450  }
1451 
1452  /* internal support */
1453  template <class T, class U>
1454  void asg_mod_v(T * dt, size_t sz, const U& x) {
1455 #ifdef __ICC
1456 #pragma loop count(1000)
1457 #pragma ivdep
1458 #pragma vector aligned
1459 #endif
1460  for (size_t i=0; i<sz; i++) dt[i]%=x[i];
1461  }
1462 
1463  /* internal support */
1464  template <class T, class U>
1465  void asg_and_v(T * dt, size_t sz, const U& x) {
1466 #ifdef __ICC
1467 #pragma loop count(1000)
1468 #pragma ivdep
1469 #pragma vector aligned
1470 #endif
1471  for (size_t i=0; i<sz; i++) dt[i]&=x[i];
1472  }
1473 
1474  /* internal support */
1475  template <class T, class U>
1476  void asg_or_v(T * dt, size_t sz, const U& x) {
1477 #ifdef __ICC
1478 #pragma loop count(1000)
1479 #pragma ivdep
1480 #pragma vector aligned
1481 #endif
1482  for (size_t i=0; i<sz; i++) dt[i]|=x[i];
1483  }
1484 
1485  /* layout of array data, including fixed fields */
1486  template <class T>
1488  {
1489  void *allocated_pointer; //actual pointer allocated (for alignment purposes)
1490  std::size_t sz; //array size
1491  unsigned cnt; //reference count
1492  static const unsigned debug_display=10; //allows debuggers to see first debug_display elem
1493  T dt[debug_display]; //data
1494  friend class array<T>;
1495  };
1496 
1501  template <class T>
1502  class array
1503  {
1504  array_data<T> *dt;
1505 
1506  friend class WhereContext;
1507 
1509  array_data<T> *alloc(std::size_t n) const
1510  {
1511  char *p;
1512  array_data<T> *r;
1513  p = (char*)std::malloc((n-array_data<T>::debug_display) * sizeof(T) + sizeof(array_data<T>) + 16);
1514 #ifdef __ICC
1515  // we need to align data onto 16 byte boundaries
1516  size_t d = (size_t)(reinterpret_cast<array_data<T>*>(p)->dt);
1517  size_t offs = (16 - (d&15)) & 15;
1518  r=reinterpret_cast<array_data<T>*>(p+offs);
1519 #else
1520  r=reinterpret_cast<array_data<T>*>(p);
1521 #endif
1522  r->allocated_pointer=p;
1523  r->sz=n;
1524  r->cnt=1;
1525  return r;
1526  }
1527 
1529  void free(array_data<T> *p) const
1530  {
1531  assert(p);
1532  std::free(p->allocated_pointer);
1533  }
1534 
1535  void set_size(size_t s) {dt = alloc(s);}
1536 
1537  unsigned& ref() // access reference counter
1538  {
1539  assert(dt);
1540  return dt->cnt;
1541  }
1542 
1543  void release()
1544  {
1545  if (dt)
1546  {
1547  if (ref()==1)
1548  free(dt);
1549  else
1550  ref()--;
1551  }
1552  }
1553 
1554  protected:
1555 
1556  void copy() //any nonconst method needs to call this
1557  { // to implement copy-on-write semantics
1558  if (dt && ref()>1)
1559  {
1560  array_data<T>* oldData=dt;
1561  bool freeMem = ref()-- == 0;
1562  dt=alloc(size());
1563  memcpy(data(),oldData->dt,size()*sizeof(T));
1564  if (freeMem) free(oldData);
1565  }
1566  }
1567 
1568  public:
1569  typedef T value_type;
1570 
1571  explicit array(size_t s=0)
1572  {
1573  set_size(s);
1574  }
1575 
1576  array(size_t s, T val)
1577  {
1578  set_size(s);
1579  array_ns::asg(data(),size(),val);
1580  }
1581 
1582  array(const array& x)
1583  {
1584  dt=x.dt;
1585  ref()++;
1586  }
1587 
1588  template <class expr>
1589  array(const expr& e,
1590  typename enable_if< is_expression<expr>, void*>::T dummy=0)
1591  {
1592  set_size(e.size());
1593  operator=(e);
1594  }
1595 
1596  ~array() {release();}
1597 
1599  void resize(size_t s) {
1600  if (!dt || s>dt->sz || ref()>1)
1601  {
1602  release();
1603  dt = alloc(s);
1604  }
1605  dt->sz=s;
1606  }
1607 
1609  template <class V>
1610  void resize(size_t s, const V& val) {resize(s); operator=(val);}
1611 
1612  void swap(array& x) {
1613  std::swap(dt, x.dt);
1614  }
1615 
1616  T& operator[](size_t i) {assert(i<size()); copy(); return data()[i];}
1617  T operator[](size_t i) const {assert(i<size()); return data()[i];}
1618 
1620  template <class I> typename
1622  operator[](const I& i) const {return RVindex<array,I>(*this,i);}
1623  template <class I> typename
1625  operator[](const I& i) {copy(); return LVindex<array,I>(*this,i);}
1626 
1627  array& operator=(const array& x) {
1628  if (x.dt==dt) return *this;
1629  release();
1630  dt=x.dt;
1631  ref()++;
1632  return *this;
1633  }
1634 
1635  template <class expr> typename
1637  operator=(const expr& x) {
1638  if ((void*)(&x)==(void*)(this)) return *this;
1639  // since expression x may contain a reference to this, assign to a temporary
1640  array tmp(x.size());
1641  array_ns::asg_v(tmp.data(),tmp.size(),x);
1642  return (*this)=tmp;
1643  }
1644  template <class expr> typename
1646  operator+=(const expr& x) {
1647  assert(size()==x.size());
1648  copy();
1649  array_ns::asg_plus_v(data(),size(),x);
1650  return *this;
1651  }
1652  template <class expr> typename
1654  operator-=(const expr& x) {
1655  assert(size()==x.size());
1656  copy();
1657  array_ns::asg_minus_v(data(),size(),x);
1658  return *this;
1659  }
1660  template <class expr> typename
1662  operator*=(const expr& x) {
1663  assert(size()==x.size());
1664  copy();
1665  array_ns::asg_mul_v(data(),size(),x);
1666  return *this;
1667  }
1668  template <class expr> typename
1670  operator/=(const expr& x) {
1671  assert(size()==x.size());
1672  copy();
1673  array_ns::asg_div_v(data(),size(),x);
1674  return *this;
1675  }
1676  template <class expr> typename
1678  operator%=(const expr& x) {
1679  assert(size()==x.size());
1680  copy();
1681  array_ns::asg_mod_v(data(),size(),x);
1682  return *this;
1683  }
1684  template <class expr> typename
1686  operator&=(const expr& x) {
1687  assert(size()==x.size());
1688  copy();
1689  array_ns::asg_and_v(data(),size(),x);
1690  return *this;
1691  }
1692  template <class expr> typename
1694  operator|=(const expr& x) {
1695  assert(size()==x.size());
1696  copy();
1697  array_ns::asg_or_v(data(),size(),x);
1698  return *this;
1699  }
1700 
1702  template <class E>
1703  typename enable_if<is_expression<E>,array&>::T
1704  operator<<=(const E& x) {
1705  array orig(*this);
1706  resize(orig.size()+x.size());
1707  asg_v(data(),orig.size(),orig);
1708  asg_v(data()+orig.size(),x.size(),x);
1709  return *this;
1710  }
1711 
1712  template <class S>
1713  typename enable_if<is_scalar<S>,array&>::T
1714  operator<<=(S x) {
1715  array orig(*this);
1716  resize(orig.size()+1);
1717  asg_v(data(),orig.size(),orig);
1718  data()[orig.size()]=x;
1719  return *this;
1720  }
1721 
1723  template <class scalar> typename
1725  operator=(scalar x) {copy(); asg(data(),size(),T(x)); return *this;}
1726 
1728  template <class scalar> typename
1730  operator+=(scalar x) {copy(); asg_plus(data(),size(),T(x)); return *this;}
1731 
1733  template <class scalar> typename
1735  operator-=(scalar x) {copy(); asg_minus(data(),size(),T(x)); return *this;}
1736 
1738  template <class scalar> typename
1740  operator*=(scalar x) {copy(); asg_mul(data(),size(),T(x)); return *this;}
1741 
1743  template <class scalar> typename
1745  operator/=(scalar x) {copy(); asg_div(data(),size(),T(x)); return *this;}
1746 
1748  template <class scalar> typename
1750  operator%=(scalar x) {copy(); asg_mod(data(),size(),T(x)); return *this;}
1751 
1753  template <class scalar> typename
1755  operator&=(scalar x) {copy(); asg_and(data(),size(),T(x)); return *this;}
1756 
1758  template <class scalar> typename
1760  operator|=(scalar x) {copy(); asg_or(data(),size(),T(x)); return *this;}
1761 
1762 // binop<array,double,Mul<array,double> > >::T
1763 // operator/(double x) const {return array_ns::operator*(*this,1/x);}
1764 
1766  size_t size() const {return dt? dt->sz: 0;}
1768  T* data() {copy(); return dt? dt->dt: 0;}
1770  const T* data() const {return dt? dt->dt: 0;}
1771 
1772  typedef T *iterator;
1773  typedef const T *const_iterator;
1774 
1775  iterator begin(void) {copy(); return iterator(data()); }
1776  iterator end(void) {copy(); return iterator(data()+size()); }
1777 
1778  const_iterator begin(void) const { return const_iterator(data()); }
1779  const_iterator end(void) const { return const_iterator(data() + size()); }
1780  };
1781 
1782 #ifdef MKL
1783 
1784  inline array<float> operator/(float x, const array<float>& y)
1785  {
1786  array<float> r(y.size());
1787  vsInv(static_cast<int>(y.size()),y.data(),r.data());
1788  if (x!=1) r*=x;
1789  return r;
1790  }
1791 
1792  inline array<float> operator/(double x, const array<float>& y)
1793  {return operator/(static_cast<float>(x),y);}
1794 
1795  inline array<double> operator/(double x, const array<double>& y)
1796  {
1797  array<double> r(y.size());
1798  vdInv(static_cast<int>(y.size()),y.data(),r.data());
1799  if (x!=1) r*=x;
1800  return r;
1801  }
1802 
1803  inline array<float> operator/(const array<float>& x, const array<float>& y)
1804  {
1805  conformance_check(x,y);
1806  array<float> r(x.size());
1807  vsDiv(static_cast<int>(x.size()),x.data(),y.data(),r.data());
1808  return r;
1809  }
1810 
1811  inline array<double> operator/(const array<double>& x, const array<double>& y)
1812  {
1813  conformance_check(x,y);
1814  array<double> r(x.size());
1815  vdDiv(static_cast<int>(x.size()),x.data(),y.data(),r.data());
1816  return r;
1817  }
1818 
1819  /* MKL can perform inplace operations */
1820  inline void asg_div_v(float *dt, size_t sz, array<float>& x)
1821  {vsDiv(static_cast<int>(sz),dt,x.data(),dt);}
1822 
1823  inline void asg_div_v(double *dt, size_t sz, array<double>& x)
1824  {vdDiv(static_cast<int>(sz),dt,x.data(),dt);}
1825 
1826  inline array<float> sqrt(const array<float>& x)
1827  {
1828  array<float> r(x.size());
1829  vsSqrt(static_cast<int>(x.size()),x.data(),r.data());
1830  return r;
1831  }
1832 
1833  inline array<double> sqrt(const array<double>& x)
1834  {
1835  array<double> r(x.size());
1836  vdSqrt(static_cast<int>(x.size()),x.data(),r.data());
1837  return r;
1838  }
1839 
1840  inline array<float> pow(const array<float>& x, float y)
1841  {
1842  array<float> r(x.size());
1843  vsPowx(static_cast<int>(x.size()),x.data(),y,r.data());
1844  return r;
1845  }
1846 
1847  inline array<float> pow(const array<float>& x, double y)
1848  {return pow(x,static_cast<float>(y));}
1849 
1850  inline array<double> pow(const array<double>& x, double y)
1851  {
1852  array<double> r(x.size());
1853  vdPowx(static_cast<int>(x.size()),x.data(),y,r.data());
1854  return r;
1855  }
1856 
1857  inline array<float> pow(const array<float>& x, const array<float>& y)
1858  {
1859  array<float> r(x.size());
1860  conformance_check(x,y);
1861  vsPow(static_cast<int>(x.size()),x.data(),y.data(),r.data());
1862  return r;
1863  }
1864 
1865  inline array<double> pow(const array<double>& x, const array<double>& y)
1866  {
1867  array<double> r(x.size());
1868  conformance_check(x,y);
1869  vdPow(static_cast<int>(x.size()),x.data(),y.data(),r.data());
1870  return r;
1871  }
1872 
1873  inline array<float> exp(const array<float>& x)
1874  {
1875  array<float> r(x.size());
1876  vsExp(static_cast<int>(x.size()),x.data(),r.data());
1877  return r;
1878  }
1879 
1880  inline array<double> exp(const array<double>& x)
1881  {
1882  array<double> r(x.size());
1883  vdExp(static_cast<int>(x.size()),x.data(),r.data());
1884  return r;
1885  }
1886 
1887  inline array<float> log(const array<float>& x)
1888  {
1889  array<float> r(x.size());
1890  vsLn(static_cast<int>(x.size()),x.data(),r.data());
1891  return r;
1892  }
1893 
1894  inline array<double> log(const array<double>& x)
1895  {
1896  array<double> r(x.size());
1897  vdLn(static_cast<int>(x.size()),x.data(),r.data());
1898  return r;
1899  }
1900 
1901  inline array<float> log10(const array<float>& x)
1902  {
1903  array<float> r(x.size());
1904  vsLog10(static_cast<int>(x.size()),x.data(),r.data());
1905  return r;
1906  }
1907 
1908  inline array<double> log10(const array<double>& x)
1909  {
1910  array<double> r(x.size());
1911  vdLog10(static_cast<int>(x.size()),x.data(),r.data());
1912  return r;
1913  }
1914 
1915  inline array<float> cos(const array<float>& x)
1916  {
1917  array<float> r(x.size());
1918  vsCos(static_cast<int>(x.size()),x.data(),r.data());
1919  return r;
1920  }
1921 
1922  inline array<double> cos(const array<double>& x)
1923  {
1924  array<double> r(x.size());
1925  vdCos(static_cast<int>(x.size()),x.data(),r.data());
1926  return r;
1927  }
1928 
1929  inline array<float> sin(const array<float>& x)
1930  {
1931  array<float> r(x.size());
1932  vsSin(static_cast<int>(x.size()),x.data(),r.data());
1933  return r;
1934  }
1935 
1936  inline array<double> sin(const array<double>& x)
1937  {
1938  array<double> r(x.size());
1939  vdSin(static_cast<int>(x.size()),x.data(),r.data());
1940  return r;
1941  }
1942 
1943  inline array<float> tan(const array<float>& x)
1944  {
1945  array<float> r(x.size());
1946  vsTan(static_cast<int>(x.size()),x.data(),r.data());
1947  return r;
1948  }
1949 
1950  inline array<double> tan(const array<double>& x)
1951  {
1952  array<double> r(x.size());
1953  vdTan(static_cast<int>(x.size()),x.data(),r.data());
1954  return r;
1955  }
1956 
1957  inline array<float> asin(const array<float>& x)
1958  {
1959  array<float> r(x.size());
1960  vsAsin(static_cast<int>(x.size()),x.data(),r.data());
1961  return r;
1962  }
1963 
1964  inline array<double> asin(const array<double>& x)
1965  {
1966  array<double> r(x.size());
1967  vdAsin(static_cast<int>(x.size()),x.data(),r.data());
1968  return r;
1969  }
1970 
1971  inline array<float> acos(const array<float>& x)
1972  {
1973  array<float> r(x.size());
1974  vsAcos(static_cast<int>(x.size()),x.data(),r.data());
1975  return r;
1976  }
1977 
1978  inline array<double> acos(const array<double>& x)
1979  {
1980  array<double> r(x.size());
1981  vdAcos(static_cast<int>(x.size()),x.data(),r.data());
1982  return r;
1983  }
1984 
1985  inline array<float> atan(const array<float>& x)
1986  {
1987  array<float> r(x.size());
1988  vsAtan(static_cast<int>(x.size()),x.data(),r.data());
1989  return r;
1990  }
1991 
1992  inline array<double> atan(const array<double>& x)
1993  {
1994  array<double> r(x.size());
1995  vdAtan(static_cast<int>(x.size()),x.data(),r.data());
1996  return r;
1997  }
1998 
1999  inline array<float> atan2(const array<float>& x, const array<float>& y)
2000  {
2001  array<float> r(x.size());
2002  conformance_check(x,y);
2003  vsAtan2(static_cast<int>(x.size()),x.data(),y.data(),r.data());
2004  return r;
2005  }
2006 
2007  inline array<double> atan2(const array<double>& x, const array<double>& y)
2008  {
2009  array<double> r(x.size());
2010  conformance_check(x,y);
2011  vdAtan2(static_cast<int>(x.size()),x.data(),y.data(),r.data());
2012  return r;
2013  }
2014 
2015  inline array<float> cosh(const array<float>& x)
2016  {
2017  array<float> r(x.size());
2018  vsCosh(static_cast<int>(x.size()),x.data(),r.data());
2019  return r;
2020  }
2021 
2022  inline array<double> cosh(const array<double>& x)
2023  {
2024  array<double> r(x.size());
2025  vdCosh(static_cast<int>(x.size()),x.data(),r.data());
2026  return r;
2027  }
2028 
2029  inline array<float> sinh(const array<float>& x)
2030  {
2031  array<float> r(x.size());
2032  vsSinh(static_cast<int>(x.size()),x.data(),r.data());
2033  return r;
2034  }
2035 
2036  inline array<double> sinh(const array<double>& x)
2037  {
2038  array<double> r(x.size());
2039  vdSinh(static_cast<int>(x.size()),x.data(),r.data());
2040  return r;
2041  }
2042 
2043  inline array<float> tanh(const array<float>& x)
2044  {
2045  array<float> r(x.size());
2046  vsTanh(static_cast<int>(x.size()),x.data(),r.data());
2047  return r;
2048  }
2049 
2050  inline array<double> tanh(const array<double>& x)
2051  {
2052  array<double> r(x.size());
2053  vdTanh(static_cast<int>(x.size()),x.data(),r.data());
2054  return r;
2055  }
2056 
2057  /* apply the above functions to expressions */
2058 
2059  template <class E> typename enable_if<is_expression<E>, void >::T
2060  asg_div_v(typename E::value_type *dt, size_t sz, const E& y)
2061  {
2062  array<typename E::value_type> Y(y.size()); Y=y;
2063  asg_div_v(dt,sz,Y);
2064  }
2065 
2066  template <class E> typename enable_if<is_expression<E>, array<typename E::value_type> >::T
2067  pow(const E& x, typename E::value_type y)
2068  {
2070  return pow(X,y);
2071  }
2072 
2073  template <class E1, class E2> typename
2074  enable_if<one_is_expression<E1,E2>,
2076  pow(const E1& x, const E2& y)
2077  {
2079  X(x.size()), Y(y.size()); X=x; Y=y;
2080  return pow(X,Y);
2081  }
2082 
2083  template <class E> typename
2084  enable_if<is_expression<E>, array<typename E::value_type> >::T
2085  sqrt(const E& x) {array<typename E::value_type> X(x.size()); X=x; return sqrt(X);}
2086  template <class E> typename
2087  enable_if<is_expression<E>, array<typename E::value_type> >::T
2088  exp(const E& x) {array<typename E::value_type> X(x.size()); X=x; return exp(X);}
2089  template <class E> typename
2090  enable_if<is_expression<E>, array<typename E::value_type> >::T
2091  log(const E& x) {array<typename E::value_type> X(x.size()); X=x; return log(X);}
2092  template <class E> typename
2093  enable_if<is_expression<E>, array<typename E::value_type> >::T
2094  log10(const E& x) {array<typename E::value_type> X(x.size()); X=x; return log10(X);}
2095  template <class E> typename
2096  enable_if<is_expression<E>, array<typename E::value_type> >::T
2097  cos(const E& x) {array<typename E::value_type> X(x.size()); X=x; return cos(X);}
2098  template <class E> typename
2099  enable_if<is_expression<E>, array<typename E::value_type> >::T
2100  sin(const E& x) {array<typename E::value_type> X(x.size()); X=x; return sin(X);}
2101  template <class E> typename
2102  enable_if<is_expression<E>, array<typename E::value_type> >::T
2103  tan(const E& x) {array<typename E::value_type> X(x.size()); X=x; return tan(X);}
2104  template <class E> typename
2105  enable_if<is_expression<E>, array<typename E::value_type> >::T
2106  acos(const E& x) {array<typename E::value_type> X(x.size()); X=x; return acos(X);}
2107  template <class E> typename
2108  enable_if<is_expression<E>, array<typename E::value_type> >::T
2109  asin(const E& x) {array<typename E::value_type> X(x.size()); X=x; return asin(X);}
2110  template <class E> typename
2111  enable_if<is_expression<E>, array<typename E::value_type> >::T
2112  atan(const E& x) {array<typename E::value_type> X(x.size()); X=x; return atan(X);}
2113 
2114  template <class E1, class E2> typename
2115  enable_if<one_is_expression<E1,E2>, array<typename result<E1,E2>::value_type> >::T
2116  atan2(const E1& x, const E2& y)
2117  {
2118  array<typename result<E1,E2>::value_type> X(x.size()), Y(y.size()); X=x; Y=y;
2119  return atan2(X,Y);
2120  }
2121 
2122  template <class E> typename
2123  enable_if<is_expression<E>, array<typename E::value_type> >::T
2124  cosh(const E& x) {array<typename E::value_type> X(x.size()); X=x; return cosh(X);}
2125  template <class E> typename
2126  enable_if<is_expression<E>, array<typename E::value_type> >::T
2127  sinh(const E& x) {array<typename E::value_type> X(x.size()); X=x; return sinh(X);}
2128  template <class E> typename
2129  enable_if<is_expression<E>, array<typename E::value_type> >::T
2130  tanh(const E& x) {array<typename E::value_type> X(x.size()); X=x; return tanh(X);}
2131 
2132 #endif // MKL
2133 
2134  /* pos takes a boolean array expression x, and returns an integer array giving the
2135  indices where x[i] is true */
2136 
2137  template <class E> typename
2139  pos(const E& x, dummy<0> d=0)
2140  {
2141  array<bool> cond=x;
2142  size_t ntrue=0;
2143  for (size_t i=0; i<cond.size(); i++) ntrue+=cond[i];
2144  array<size_t> r(ntrue);
2145  for (size_t i=0, j=0; i<cond.size(); i++)
2146  if (cond[i]) r[j++]=i;
2147  return r;
2148  }
2149 
2150  template <class S> typename
2151  enable_if< is_scalar<S>, size_t >::T
2152  pos(const S& x, dummy<1> d=1) {return 0;}
2153 
2154  /*
2155  Boolean short cuts
2156  */
2157 
2158  /*
2159  We don't want to evaluate the second argument if first argument is false (&&-case) or
2160  true (||-case). In such a case, return a varArray<bool> sized as a fixedArray<bool>
2161  (since we're expecting this to be used in fixedArray expressions)
2162  */
2163 
2164  template <class E>
2166  operator&&(bool x, const E& e)
2167  {
2168  array<bool> r;
2169  if (x)
2170  r=e;
2171  else
2172  r=array<bool>(1,false);
2173  return r;
2174  }
2175 
2176  template <class T, class E>
2178  operator&&(T *x, const E& e)
2179  {return x!=NULL && e;}
2180 
2181  template <class E>
2183  operator||(bool x, const E& e)
2184  {
2185  array<bool> r;
2186  if (!x)
2187  r=e;
2188  else
2189  r=array<bool>(1,true);
2190  return r;
2191  }
2192 
2193  /* reductions */
2194 
2196  template <class E>
2197  typename enable_if< is_expression<E>, typename E::value_type >::T
2198  max(const E& x)
2199  {
2200  typename E::value_type r;
2201  //strangely numeric_limits<float>::min() differs semantically from numeric_limits<int>::min()
2202  if (std::numeric_limits<typename E::value_type>::is_integer)
2203  r=std::numeric_limits<typename E::value_type>::min();
2204  else
2205  r=-std::numeric_limits<typename E::value_type>::max();
2206  for (size_t i=0; i<x.size(); i++)
2207  r = std::max(x[i],r);
2208  return r;
2209  }
2210 
2212  template <class E>
2213  typename enable_if< is_expression<E>, typename E::value_type >::T
2214  min(const E& x)
2215  {
2216  typename E::value_type r=std::numeric_limits<typename E::value_type>::max();
2217  for (size_t i=0; i<x.size(); i++)
2218  r = std::min(x[i],r);
2219  return r;
2220  }
2221 
2223  template <class E>
2224  typename enable_if< is_expression<E>,
2225  typename result<typename E::value_type, typename E::value_type>::value_type
2226  >::T
2227  sum(const E& x)
2228  {
2229  typename result<typename E::value_type, typename E::value_type>::value_type
2230  r=0;
2231  for (size_t i=0; i<x.size(); i++)
2232  r += x[i];
2233  return r;
2234  }
2235 
2236  template <class E, class M>
2237  typename enable_if< is_expression<E>,
2238  typename result<typename E::value_type, typename E::value_type>::value_type
2239  >::T
2240  sum(const E& x, const M& mask)
2241  {
2242  typename result<typename E::value_type, typename E::value_type>::value_type
2243  r=0;
2244  for (size_t i=0; i<x.size(); i++)
2245  if (mask[i]) r += x[i];
2246  return r;
2247  }
2248 
2250  template <class E>
2251  typename enable_if< is_expression<E>, typename E::value_type >::T
2252  prod(const E& x)
2253  {
2254  typename E::value_type r=1;
2255  for (size_t i=0; i<x.size(); i++)
2256  r *= x[i];
2257  return r;
2258  }
2259 
2261  template <class E, class M>
2262  typename enable_if< is_expression<E>, typename E::value_type >::T
2263  prod(const E& x, const M& mask)
2264  {
2265  typename E::value_type r=1;
2266  for (size_t i=0; i<x.size(); i++)
2267  if (mask[i]) r *= x[i];
2268  return r;
2269  }
2270 
2272  template <class E>
2273  bool any(const E& x)
2274  {
2275  for (size_t i=0; i<x.size(); i++)
2276  if(x[i])
2277  return true;
2278 
2279  return false;
2280  }
2281 
2283  template <class E>
2284  bool all(const E& x)
2285  {
2286  for (size_t i=0; i<x.size(); i++)
2287  if(!x[i])
2288  return false;
2289 
2290  return true;
2291  }
2292 
2294  template <class E1, class E2, class E3> typename
2297  >::T
2298  merge(const E1& mask, const E2& x, const E3& y)
2299  {
2300  conformance_check(mask,x);
2301  conformance_check(mask,y);
2302  array<typename result<E1,E2>::value_type > ret(mask.size());
2303  for (size_t i=0; i<mask.size(); i++)
2304  ret[i]= mask[i]? at(x,i): at(y,i);
2305  return ret;
2306  }
2307 
2309  template <class E1, class E2> typename
2312  >::T
2313  pack(const E1& e, const E2& mask, long ntrue=-1)
2314  {
2315  assert(mask.size()==e.size());
2316  array<typename E1::value_type> r( (ntrue==-1)? sum(mask): ntrue);
2317  for (size_t i=0, j=0; i<mask.size(); i++)
2318  if (mask[i]) r[j++]=e[i];
2319  return r;
2320  }
2321 
2323  template <class E> typename
2324  enable_if<is_expression<E>, array<size_t> >::T
2325  enumerate(const E& x)
2326  {
2327  array<size_t> r(x.size());
2328  for (size_t i=0, j=0; i<x.size(); i++)
2329  {
2330  r[i]=j;
2331  if (x[i]) j++;
2332  }
2333  return r;
2334  }
2335 
2336  template <class E>
2337  std::ostream& put(std::ostream& o, const E& x)
2338  {
2339  if (x.size())
2340  {
2341  for (size_t i=0; i<x.size()-1; i++)
2342  o << x[i] << " ";
2343  o<<x[x.size()-1];
2344  }
2345  return o;
2346  }
2347 
2348  template <class T>
2349  std::ostream& operator<<(std::ostream& o, const array<T>& x)
2350  {return put(o,x);}
2351 
2353  template <class T, class Op>
2354  std::ostream& operator<<(std::ostream& o, const unop<T, Op>& x)
2355  {return put(o,x);}
2356 
2357  template <class E1, class E2, class Op>
2358  std::ostream& operator<<(std::ostream& o, const binop<E1,E2, Op>& x)
2359  {return put(o,x);}
2360 
2362  template <class T>
2363  std::istream& get(std::istream& s, T& x)
2364  {
2365  typename T::value_type v; x.resize(0);
2366  while (s>>v) {x<<=v;}
2367  return s;
2368  }
2369 
2370  template <class T>
2371  std::istream& operator>>(std::istream& i, array<T>& x)
2372  {return get(i,x);}
2373 
2375  inline array<int> pcoord(int n)
2376  {
2377  array<int> r(n);
2378 #ifdef __ICC
2379 #pragma loop count(1000)
2380 #pragma ivdep
2381 #pragma vector aligned
2382 #endif
2383  for (int i=0; i<n; i++)
2384  r[i]=i;
2385  return r;
2386  }
2387 
2398  {
2399  array<int> r(sum(x));
2400  for (size_t i=0,p=0; i<x.size(); i++)
2401  for (int j=0; j<x[i]; j++,p++)
2402  r[p]=i;
2403  return r;
2404  }
2405 
2406 
2407 
2409  void fillrand(array<double>& x);
2411  void fillgrand(array<double>& x);
2413  void fillprand(array<double>& x);
2415  void fill_unique_rand(array<int>& x, unsigned max);
2416 
2418  void lgspread( array<double>& a, const array<double>& s );
2420  void gspread( array<double>& a, const array<double>& s );
2422  template <class E>
2423  void lgspread( array<double>& a, const E& s ) {lgspread(a,array<double>(s));}
2425  template <class E>
2426  void gspread( array<double>& a, const E& s ) {gspread(a,array<double>(s));}
2427 
2428  /* ranking (sort) function */
2429  enum array_dir_t {upwards, downwards};
2430 
2431  template <class T>
2432  class Cmp
2433  {
2434  const array<T>& data;
2435  bool fwd;
2436  public:
2437  array<int> ranks;
2438  Cmp(const array<T>& data, bool fwd): data(data), fwd(fwd) {
2439  ranks=pcoord(data.size());
2440  }
2441  bool operator()(const int& i, const int& j)
2442  {
2443  if (fwd)
2444  return data[i] < data[j];
2445  else
2446  return data[j] < data[i];
2447  }
2448  };
2449 
2451  template <class T> typename
2453  rank(const T& x, enum array_dir_t dir=upwards)
2454  {
2455  Cmp<typename T::value_type> cmp(x,dir==upwards);
2456  std::sort(cmp.ranks.begin(),cmp.ranks.end(),cmp);
2457  return cmp.ranks;
2458  }
2459 
2460 
2461  } // namespace array_ns
2462 
2463  using array_ns::array;
2464  using array_ns::pcoord;
2465 
2466  extern urand array_urand;
2467  extern gaussrand array_grand;
2468 } // namespace ecolab
2469 
2470 using ecolab::TCL_obj;
2471 
2472 namespace classdesc_access
2473 {
2474  /* Definitions of pack, unpack and TCL_obj for array classes */
2475  template <class T>
2476  struct access_pack<ecolab::array_ns::array<T> >
2477  {
2478  template <class U>
2479  void operator()(classdesc::pack_t& targ, const classdesc::string& desc,
2480  U& arg)
2481  {
2482  pack(targ,desc,arg.size());
2483  pack(targ,desc,classdesc::is_array(),*arg.data(),1,arg.size());
2484  }
2485  };
2486 
2487  template <class T>
2488  struct access_unpack<ecolab::array_ns::array<T> >
2489  {
2490  template <class U>
2491  void operator()(classdesc::unpack_t& targ, const classdesc::string& desc,
2492  U& arg)
2493  {
2494  size_t size;
2495  unpack(targ,desc,size);
2496  arg.resize(size);
2497  unpack(targ,desc,classdesc::is_array(),*arg.data(),1,size);
2498  }
2499  };
2500 
2501  template <class T>
2502  struct access_TCL_obj<ecolab::array_ns::array<T> >
2503  {
2504  template <class U>
2505  void operator()(ecolab::TCL_obj_t& targ, const classdesc::string& desc,
2506  U& arg)
2507  {
2508  ecolab::TCL_obj_register(targ,desc,arg);
2509  TCL_obj(targ,desc+".size", arg, &ecolab::array_ns::array<T>::size);
2510  }
2511  };
2512 }
2513 
2514 namespace classdesc
2515 {
2516  template <class T>
2517  struct tn<ecolab::array<T> >
2518  {
2519  static std::string name()
2520  {return "ecolab::array<"+typeName<T>()+">";}
2521  };
2522 }
2523 
2524 #endif
Definition: arrays.h:1086
unop< E, Tanh< typename E::value_type > > tanh(const E &e)
Definition: arrays.h:952
binop< E1, E2, Pow< E1, E2 > > pow(const E1 &e1, const E2 &e2)
Definition: arrays.h:777
Definition: arrays.h:481
unop< E, Floor< typename E::value_type > > floor(const E &e)
Definition: arrays.h:1027
Definition: arrays.h:1066
Definition: arrays.h:91
Definition: arrays.h:192
enable_if< both_are_expressions< E1, E2 >, array< typename E1::value_type > >::T pack(const E1 &e, const E2 &mask, long ntrue=-1)
pack vector (remove elements where mask[i] is false
Definition: arrays.h:2313
Definition: arrays.h:957
unop< E, Acos< typename E::value_type > > acos(const E &e)
Definition: arrays.h:883
Definition: arrays.h:501
enable_if< one_is_expression< E1, E2 >, binop< E1, E2, Max< E1, E2 > > >::T max(const E1 &e1, const E2 &e2)
Definition: arrays.h:1079
Definition: arrays.h:931
Definition: arrays.h:79
Definition: arrays.h:471
enable_if< is_expression< E >, typename E::value_type >::T prod(const E &x)
product of a vector expression
Definition: arrays.h:2252
enable_if< is_scalar< scalar >, array & >::T operator-=(scalar x)
subtracting broadcast
Definition: arrays.h:1735
unop< E, Exp< typename E::value_type > > exp(const E &e)
Definition: arrays.h:792
enable_if< one_is_expression< E1, E2 >, binop< E1, E2, Min< E1, E2 > > >::T min(const E1 &e1, const E2 &e2)
Definition: arrays.h:1099
Definition: arrays.h:441
unop< E, Cosh< typename E::value_type > > cosh(const E &e)
Definition: arrays.h:939
unop< E, Log< typename E::value_type > > log(const E &e)
Definition: arrays.h:805
Namespace for array functionality.
size_t size() const
size of buffer
Definition: pack_base.h:154
enable_if< is_scalar< scalar >, array & >::T operator|=(scalar x)
bitwise or broadcast
Definition: arrays.h:1760
Definition: arrays.h:797
Definition: arrays.h:2432
E::value_type & at(E &x, size_t i, typename enable_if< is_expression< E >, int >::T dum=0)
Definition: arrays.h:335
void fillgrand(array< double > &x)
fill array with gaussian random numbers from
Definition: arrays.h:398
enable_if< is_scalar< scalar >, array & >::T operator/=(scalar x)
dividing broadcast
Definition: arrays.h:1745
Definition: classdesc.h:375
Definition: arrays.h:80
binop< E1, E2, Ldexp< E1, E2 > > ldexp(const E1 &e1, const E2 &e2)
Definition: arrays.h:1042
enable_if< is_expression< E1 >, array< typename result< E2, E3 >::value_type > >::T merge(const E1 &mask, const E2 &x, const E3 &y)
merge two vectors: r[i] = mask[i]? x[i]: y[i]
Definition: arrays.h:2298
Definition: arrays.h:146
Definition: arrays.h:82
serialisation descriptor
void gspread(array< double > &a, const array< double > &s)
Additive process .
unop< E, Abs< typename E::value_type > > abs(const E &e)
Definition: arrays.h:992
Definition: arrays.h:491
unop< E, Asin< typename E::value_type > > asin(const E &e)
Definition: arrays.h:870
Definition: arrays.h:431
Definition: arrays.h:888
unop< E, Abs< typename E::value_type > > fabs(const E &e)
Definition: arrays.h:987
Definition: arrays.h:1006
void fillrand(array< double > &x)
fill array with uniformly random numbers from [0,1)
Definition: arrays.h:127
Definition: classdesc.h:623
void lgspread(array< double > &a, const array< double > &s)
Multiplicative process .
binop< E1, E2, Atan2< E1, E2 > > atan2(const E1 &e1, const E2 &e2)
Definition: arrays.h:911
Definition: arrays.h:104
enable_if< is_scalar< scalar >, array & >::T operator+=(scalar x)
summing broadcast
Definition: arrays.h:1730
bool eq(const E1 &e1, const E2 &e2)
Definition: arrays.h:747
unop< E, Sgn< typename E::value_type > > sign(const E &e)
Definition: arrays.h:1001
Definition: arrays.h:1049
Gaussian (normal) random generator.
Definition: random_basic.h:19
class to allow access to private members
Definition: classdesc_access.h:21
Definition: arrays.h:157
Definition: arrays.h:918
enable_if< is_expression< E >, size_t >::T len(const E &x, dummy< 0 > d=0)
Definition: arrays.h:359
Definition: arrays.h:409
Definition: arrays.h:83
enable_if< is_scalar< scalar >, array & >::T operator*=(scalar x)
multiplying broadcast
Definition: arrays.h:1740
enable_if< is_scalar< scalar >, array & >::T operator=(scalar x)
broadcast
Definition: arrays.h:1725
unop< E, Ceil< typename E::value_type > > ceil(const E &e)
Definition: arrays.h:1014
binop< E1, E2, Fmod< E1, E2 > > fmod(const E1 &e1, const E2 &e2)
Definition: arrays.h:1059
array< int > gen_index(const array< int > &x)
Definition: arrays.h:2397
Definition: arrays.h:944
E::value_type value_type
type of result
Definition: arrays.h:1220
void resize(size_t s, const V &val)
resize array to s elements, and initialise to val
Definition: arrays.h:1610
Definition: arrays.h:451
Definition: TCL_obj_stl.h:185
Definition: arrays.h:511
Definition: arrays.h:1019
class to allow access to private members
Definition: classdesc_access.h:22
Definition: arrays.h:849
bool all(const E &x)
true iff all elements of x are true
Definition: arrays.h:2284
Definition: arrays.h:161
Definition: arrays.h:154
enable_if< is_expression< T >, array< int > >::T rank(const T &x, enum array_dir_t dir=upwards)
rank elements
Definition: arrays.h:2453
void fill_unique_rand(array< int > &x, unsigned max)
fill with uniform numbers drawn from [0...max] without replacement
unop< E, Log10< typename E::value_type > > log10(const E &e)
Definition: arrays.h:818
Definition: arrays.h:901
unop< E, Tan< typename E::value_type > > tan(const E &e)
Definition: arrays.h:857
void fillprand(array< double > &x)
fill array with exponential random numbers
TCL access descriptor.
Definition: arrays.h:389
enable_if< is_scalar< scalar >, array & >::T operator &=(scalar x)
bitwise and broadcast
Definition: arrays.h:1755
Definition: arrays.h:765
Definition: arrays.h:171
enable_if< is_expression< I >, RVindex< unop, I > >::T operator[](const I &i) const
vector indexing
Definition: arrays.h:1228
Definition: arrays.h:784
enable_if< is_expression< I >, RVindex< binop, I > >::T operator[](const I &i) const
vector indexing
Definition: arrays.h:1278
unop< E, Sin< typename E::value_type > > sin(const E &e)
Definition: arrays.h:831
Definition: arrays.h:810
enable_if< is_expression< I >, RVindex< array, I > >::T operator[](const I &i) const
vector indexing
Definition: arrays.h:1622
Definition: TCL_obj_base.h:364
array< int > pcoord(int n)
[0...n-1]
Definition: arrays.h:2375
unop< E, Cos< typename E::value_type > > cos(const E &e)
Definition: arrays.h:844
unop< E, Sqrt< typename E::value_type > > sqrt(const E &e)
Definition: arrays.h:965
size_t size() const
number of elements
Definition: arrays.h:1766
Contains definitions related to classdesc functionality.
Definition: arrays.h:2514
Definition: pack_base.h:124
enable_if< is_expression< E >, typename result< typename E::value_type, typename E::value_type >::value_type >::T sum(const E &x)
sum of a vector expression
Definition: arrays.h:2227
Definition: arrays.h:521
unop< E, Sinh< typename E::value_type > > sinh(const E &e)
Definition: arrays.h:926
Definition: arrays.h:81
_OPENMP
Definition: accessor.h:16
Definition: arrays.h:875
Definition: arrays.h:420
Definition: arrays.h:979
enable_if< is_scalar< scalar >, array & >::T operator%=(scalar x)
remaindering broadcast
Definition: arrays.h:1750
Definition: arrays.h:380
TCL_obj descriptor object.
Definition: TCL_obj_base.h:327
Definition: arrays.h:836
Definition: arrays.h:461
Definition: arrays.h:1487
Definition: TCL_obj_base.h:27
Definition: arrays.h:971
void resize(size_t s)
resize array to s elements
Definition: arrays.h:1599
Definition: arrays.h:1032
Contains access_* structs, and nothing else. These structs are used to gain access to private members...
Definition: accessor.h:55
unop< E, Atan< typename E::value_type > > atan(const E &e)
Definition: arrays.h:896
Definition: random_basic.h:11
bool any(const E &x)
true iff any element x is true
Definition: arrays.h:2273
const T * data() const
obtain raw pointer to data
Definition: arrays.h:1770
enable_if< is_expression< E >, array< size_t > >::T enumerate(const E &x)
running sum of x considered as a boolean mask
Definition: arrays.h:2325
Op::value_type value_type
type of result
Definition: arrays.h:1243
Definition: arrays.h:862
Definition: arrays.h:112
Definition: arrays.h:823
Definition: arrays.h:198
T * data()
obtain raw pointer to data
Definition: arrays.h:1768