TBCI Numerical high perf. C++ Library  2.8.0
vector.h
Go to the documentation of this file.
1 
13 #ifndef TBCI_VECTOR_H
14 #define TBCI_VECTOR_H
15 
16 #include "tbci/bvector.h"
17 //#include "cplx.h"
18 
19 // Include vector kernels
20 #include "tbci/vec_kern_unr_pref.h"
21 
26 #ifndef TBCI_NO_SIMD
27 # if defined(__AVX__) && !defined(NO_AVX)
28 # include "tbci/vec_kern_special2.h"
29 # else
30 # include "tbci/vec_kern_special.h"
31 # endif
32 #endif
33 
34 // Avoid -fguiding-decls
35 #if !defined(NO_GD) && !defined(AUTO_DECL)
36 # include "tbci/vector_gd.h"
37 #endif
38 
39 NAMESPACE_TBCI
40 
41 //template <typename T> class cplx;
42 
43 template <typename T> class Vector;
44 template <typename T> class TVector;
45 template <typename T> class TSVector;
46 template <typename T> class Matrix;
47 template <typename T> class TMatrix;
48 template <typename T> class TSMatrix;
49 template <typename T> class F_Matrix;
50 template <typename T> class F_TMatrix;
51 template <typename T> class F_TSMatrix;
52 template <typename T> class Mat_Brack;
53 template <typename T> class CRMatrix;
54 template <typename T> class CSCMatrix;
55 template <typename T> class BdMatrix;
56 template <typename T> class Symm_BdMatrix;
57 
58 #ifdef PRAGMA_I
59 # pragma interface "vector.h"
60 #endif
61 
62 //------------------------------------------------------------------------
63 
71 template <typename T>
72 class TVector : public BVector<T>, public Vector_Sig<T>
73 {
74  public:
75  friend class Vector<T>;
76  friend class TSVector<T>; // shouldn't be necessary ...
77  friend class Matrix<T>;
78  friend class BdMatrix<T>;
79 
80  typedef T value_type;
81  typedef T element_type;
82  typedef T aligned_value_type TALIGN(MIN_ALIGN2);
83 
84  explicit TVector (const unsigned long d = 0) : BVector<T> (d) {}
85  TVector (const T& val, const unsigned long d) : BVector<T> (val, d) {}
86  TVector (const BVector<T>& bv) : BVector<T> (bv) {}
87  TVector (const Vector<T>& v) : BVector<T> ( v) {}
88  //alias
89  TVector (const TVector<T>& tv) HOT
90  { this->dim = tv.dim; this->vec = tv.vec; }
91  explicit TVector (const TSVector<T>& ts)
92  { ts.eval (); this->vec = ts.vec; this->dim = ts.dim; }
93  TVector (const Mat_Brack<T>&);
94  // dtor: Don't destroy!
95  ~TVector () { this->keep = true; }
96  /* TODO: Check whether it would be more efficient to save time by
97  * setting keep in c'tors or even only when aliasing takes place;
98  * we'd be doing refcounting more explicitly this way.
99  * Question is whether we'd not end up with full ref counting, which
100  * is what TBCI tries to avoid. -KG, 2004-04-27.
101  */
102 
103  // avoid ambiguity
104  unsigned long size () const { return this->dim; }
105 
106  TVector<T>& operator = (const T& a)
107  { this->BVector<T>::operator = (a); return *this; }
108  TVector<T>& operator = (const Vector<T>& v)
109  { this->BVector<T>::operator = (v); return *this; }
111  TVector<T>& operator = (const TVector<T>& tv);
112  TVector<T>& operator = (const TSVector<T>& ts);
113 
114  bool operator == (const TVector<T>& tv) const
115  { Vector<T> v(*this); if (LIKELY(this->vec == tv.vec)) return true;
116  /* else { */ Vector<T> v2(tv); return (v == v2); /*}*/ }
117 
118  bool operator != (const TVector<T>& tv) const { return !(*this == tv); }
119  bool operator == (const Vector<T>& v) const
120  { Vector<T> tv(*this); return (tv == v); }
121  bool operator != (const Vector<T>& v) const { return !(*this == v); }
122  bool operator == (const BVector<T>& bv) const
123  { Vector<T> v(*this); return bv.operator == (v); }
124  bool operator != (const BVector<T>& v) const { return !(*this == v); }
125  bool operator == (const TSVector<T>& tsv) const { return (tsv == *this); }
126  bool operator != (const TSVector<T>& tsv) const { return !(*this == tsv); }
127 
128  TVector<T>& operator += (const T&);
129  TVector<T>& operator -= (const T&);
130  TVector<T>& operator *= (const T&);
131  TVector<T>& operator /= (const T&);
132 
133  TVector<T>& operator += (const Vector<T>&) HOT;
134  TVector<T>& operator -= (const Vector<T>&) HOT;
135 
136  TVector<T>& operator += (const TVector<T>& tv)
137  { return this->operator += (Vector<T> (tv)); }
138  TVector<T>& operator -= (const TVector<T>& tv)
139  { return this->operator -= (Vector<T> (tv)); }
140 
141  TVector<T>& operator += (const TSVector<T>& tsv) HOT;
142  TVector<T>& operator -= (const TSVector<T>& tsv) HOT;
143 
144  // addtinional arithmetics (valarray like)
145 
146  T min () { Vector<T> th (*this); return th.min (); }
147  T max () { Vector<T> th (*this); return th.max (); }
148  T sum () { Vector<T> th (*this); return th.sum (); }
149 
150  double fabssqr () { Vector<T> th (*this); return th.fabssqr (); }
151  double fabs () { Vector<T> th (*this); return th.fabs (); }
152  T abs () { Vector<T> th (*this); return th.abs (); }
153  /* NOTE: The const is wrong! Just to make the compiler happy, when compiling
154  * operator / (const T&, const TVector<T>&) */
155 
156  // basic operations
157 
158  TVector<T>& operator + (const Vector<T>&) HOT;
159  TVector<T>& operator - (const Vector<T>&) HOT;
160  TVector<T>& operator + (const TSVector<T>& ts) HOT;
161  TVector<T>& operator - (const TSVector<T>& ts) HOT;
162  const TVector<T>& operator + (const TVector<T>& a) /*const*/ HOT;
163  const TVector<T>& operator - (const TVector<T>& a) /*const*/ HOT;
164 
165  // fortunately, we don't need friendship here
166 #ifdef TV_OP_FRIENDS
167  friend NOINST const TVector<T>& operator + FGD (const T&, const TVector<T>&);
168  friend NOINST const TVector<T>& operator - FGD (const T&, const TVector<T>&);
169  friend NOINST TSVector<T> operator * FGD (const T&, const TVector<T>&);
170  friend NOINST TSVector<T> operator / FGD (const T&, const TVector<T>&);
171  friend NOINST TVector<T> operator + FGD (const T&, const Vector<T>&);
172  friend NOINST TVector<T> operator - FGD (const T&, const Vector<T>&);
173  friend NOINST TSVector<T> operator * FGD (const T&, const Vector<T>&);
174  //friend NOINST TVector<T> operator / FGD (const T&, const Vector<T>&);
175 
176  // More friendships ...
177  friend NOINST TVector<T> operator + FGD (const T&, const TSVector<T>&);
178  friend NOINST TVector<T> operator - FGD (const T&, const TSVector<T>&);
179 #endif
180 
181  // operator () destroys TVector !!!
182  T operator () (const unsigned long i)
183  { T val = this->vec[i]; /*if (!keep)*/ this->destroy (); return val; }
184  //T& operator () (const unsigned long i) { STD__ cerr << "TVectors may not be written to by the () op\n"; abort(); }
185  T operator [] (const unsigned long i) { return this->operator() (i); }
186  //T& operator [] (const unsigned long i) { return this->operator() (i); }
187 
188  const T& getcref (const unsigned long i) const { return this->vec[i]; }
189  typename tbci_traits<T>::const_refval_type
190  get (const unsigned long i) const { return this->vec[i]; }
191 
192  T& setval (const unsigned long i) const { return this->vec[i]; }
193  T& setval (const T& val, const unsigned long i) const { return this->vec[i] = val; }
194  T& set (const T& val, const unsigned long i) const { return this->vec[i] = val; }
195 
196  TVector<T> slice (const unsigned long, const unsigned long);
197  TVector<T>& incr (const unsigned long, const T = (T)1);
198 
199  TVector<T>& operator + (const T&);
200  TVector<T>& operator - (const T&);
201  TSVector<T> operator * (const T&);
202  TSVector<T> operator / (const T&);
203  TSVector<T> operator - ();
204 
205  T operator * (const Vector<T>& v)
206  { Vector<T> th (*this); return th * v; }
207  T operator * (TVector<T>& tv)
208  { Vector<T> th (*this); Vector<T> v (tv); return th * v; }
209 
210 
211  friend TVector<T> conj FGD (const Vector<T>&);
212  friend TVector<T> real FGD (const Vector<T>&);
213  friend TVector<T> imag FGD (const Vector<T>&);
214 
215  friend TVector<T>& conj FGD (TVector<T>&);
216  friend TVector<T>& real FGD (TVector<T>&);
217  friend TVector<T>& imag FGD (TVector<T>&);
218 
219  static const char* vec_info() { return "TVector"; }
220 
221  friend STD__ ostream& operator << FGD (STD__ ostream&, const TVector<T>&);
222 
223 #ifndef HAVE_PROMOTION_BUG
224  // Promotion (only explicit)
225 //# ifndef HAVE_GCC295_TMPLFRNDCLS_BUG
226 // template <typename U> friend class BVector;
227 //# endif
228  template <typename U> explicit TVector (const BVector<U>& bv) : BVector<T> (bv) {}
229 #endif
230  /* Speed mat-vec-mul */
231  friend NOINST void FRIEND_TBCI2__ do_mat_vec_mult FGD (const unsigned start, const unsigned end, \
232  TVector<T> * res, const Matrix<T> * mat, const Vector<T> * vec);
233  friend NOINST void FRIEND_TBCI2__ do_mat_tsv_mult FGD (const unsigned start, const unsigned end, \
234  TVector<T> * res, const Matrix<T> * mat, const TSVector<T> * rsv);
235  friend NOINST void FRIEND_TBCI2__ do_mat_vec_transmult FGD (const unsigned start, const unsigned end, \
236  TVector<T> * res, const Matrix<T> * mat, const Vector<T> * vec);
237 } ;
238 
239 /* ************************ definitions ********************************** */
240 
241 template <typename T>
242 inline TVector<T>& TVector<T>::operator = (const TVector<T>& tv)
243 {
244  BCHK (this->dim != tv.dim, VecErr, diff size in assignment, tv.dim, *this);
245  if (LIKELY(this->vec != tv.vec)) {
246  this->destroy (); this->vec = tv.vec; this->dim = tv.dim;
247  }
248  return *this;
249 }
250 
251 
252 template <typename T>
253 inline STD__ ostream& operator << (STD__ ostream& os, const TVector<T>& tv)
254 { Vector<T> v(tv); return os << (BVector<T>)v; }
255 
256 template <typename T>
257 inline TVector<T>& TVector<T>::operator = (const TSVector<T>& tsv)
258 {
259  BCHK (this->dim != tsv.dim, VecErr, diff size in assignment, tsv.dim, *this);
260  if (LIKELY(tsv.mut && tsv.fac == (T)1)) {
261  this->destroy (); this->vec = tsv.vec;
262  }
263  else
264  tsv.eval (this->vec);
265  return *this;
266 }
267 
268 /*
269 template <typename T>
270 inline TVector<T>::TVector (const TSVector<T>& ts)
271 {
272  ts.eval (); vec = ts.vec; dim = ts.dim; ts.mut = false;
273 }
274  */
275 
276 #if !defined(SMP) || defined(NOSMP_VECVEC)
277 /* nothing, actually */
278 #else /* SMP */
279 
280 INST(template <typename T> class Vector friend void job_vec_vec_add (struct thr_ctrl *);)
282 template <typename T>
283 void job_vec_vec_add (struct thr_ctrl *tc)
284 {
285  do_vec_vec_add (tc->t_size, (T*)(tc->t_par[0]),
286  (const T*)(tc->t_par[1]), (const T*)(tc->t_par[2]));
287 }
288 INST(template <typename T> class Vector friend void job_vec_vec_sub (struct thr_ctrl *);)
290 template <typename T>
291 void job_vec_vec_sub (struct thr_ctrl *tc)
292 {
293  do_vec_vec_sub (tc->t_size, (T*)(tc->t_par[0]),
294  (const T*)(tc->t_par[1]), (const T*)(tc->t_par[2]));
295 }
296 INST(template <typename T> class Vector friend void job_vec_add_vec (struct thr_ctrl *);)
298 template <typename T>
299 void job_vec_add_vec (struct thr_ctrl *tc)
300 {
301  do_vec_add_vec (tc->t_size, (T*)(tc->t_par[0]),
302  (const T*)(tc->t_par[1]));
303 }
304 INST(template <typename T> class Vector friend void job_vec_sub_vec (struct thr_ctrl *);)
306 template <typename T>
307 void job_vec_sub_vec (struct thr_ctrl *tc)
308 {
309  do_vec_sub_vec (tc->t_size, (T*)(tc->t_par[0]),
310  (const T*)(tc->t_par[1]));
311 }
312 INST(template <typename T> class Vector friend void job_vec_sub_vec_inv (struct thr_ctrl *);)
314 template <typename T>
315 void job_vec_sub_vec_inv (struct thr_ctrl *tc)
316 {
317  do_vec_sub_vec_inv (tc->t_size, (T*)(tc->t_par[0]),
318  (const T*)(tc->t_par[1]));
319 }
320 INST(template <typename T> class Vector friend void job_vec_val_add (struct thr_ctrl *);)
322 template <typename T>
323 void job_vec_val_add (struct thr_ctrl *tc)
324 {
325  do_vec_val_add (tc->t_size, (T*)(tc->t_par[0]),
326  (const T*)(tc->t_par[1]), *(const T*)(tc->t_par[2]));
327 }
328 INST(template <typename T> class Vector friend void job_vec_val_sub (struct thr_ctrl *);)
330 template <typename T>
331 void job_vec_val_sub (struct thr_ctrl *tc)
332 {
333  do_vec_val_sub (tc->t_size, (T*)(tc->t_par[0]),
334  (const T*)(tc->t_par[1]), *(const T*)(tc->t_par[2]));
335 }
336 INST(template <typename T> class Vector friend void job_vec_val_mul (struct thr_ctrl *);)
338 template <typename T>
339 void job_vec_val_mul (struct thr_ctrl *tc)
340 {
341  do_vec_val_mul (tc->t_size, (T*)(tc->t_par[0]),
342  (const T*)(tc->t_par[1]), *(const T*)(tc->t_par[2]));
343 }
344 INST(template <typename T> class Vector friend void job_val_vec_mul (struct thr_ctrl *);)
346 template <typename T>
347 void job_val_vec_mul (struct thr_ctrl *tc)
348 {
349  do_val_vec_mul (tc->t_size, (T*)(tc->t_par[0]),
350  (const T*)(tc->t_par[1]), *(const T*)(tc->t_par[2]));
351 }
352 INST(template <typename T> class Vector friend void job_val_vec_add (struct thr_ctrl *);)
354 template <typename T>
355 void job_val_vec_add (struct thr_ctrl *tc)
356 {
357  do_val_vec_add (tc->t_size, (T*)(tc->t_par[0]),
358  (const T*)(tc->t_par[1]), *(const T*)(tc->t_par[2]));
359 }
360 INST(template <typename T> class Vector friend void job_val_vec_sub (struct thr_ctrl *);)
362 template <typename T>
363 void job_val_vec_sub (struct thr_ctrl *tc)
364 {
365  do_val_vec_sub (tc->t_size, (T*)(tc->t_par[0]),
366  (const T*)(tc->t_par[1]), *(const T*)(tc->t_par[2]));
367 }
368 INST(template <typename T> class Vector friend void job_val_vec_div (struct thr_ctrl *);)
370 template <typename T>
371 void job_val_vec_div (struct thr_ctrl *tc)
372 {
373  do_val_vec_div (tc->t_size, (T*)(tc->t_par[0]),
374  (const T*)(tc->t_par[1]), *(const T*)(tc->t_par[2]));
375 }
376 INST(template <typename T> class Vector friend void job_vec_add_val (struct thr_ctrl *);)
378 template <typename T>
379 void job_vec_add_val (struct thr_ctrl *tc)
380 {
381  do_vec_add_val (tc->t_size, (T*)(tc->t_par[0]),
382  *(const T*)(tc->t_par[1]));
383 }
384 INST(template <typename T> class Vector friend void job_val_add_vec (struct thr_ctrl *);)
386 template <typename T>
387 void job_val_add_vec (struct thr_ctrl *tc)
388 {
389  do_val_add_vec(tc->t_size, (T*)(tc->t_par[0]),
390  *(const T*)(tc->t_par[1]));
391 }
392 INST(template <typename T> class Vector friend void job_vec_sub_val (struct thr_ctrl *);)
394 template <typename T>
395 void job_vec_sub_val (struct thr_ctrl *tc)
396 {
397  do_vec_sub_val (tc->t_size, (T*)(tc->t_par[0]),
398  *(const T*)(tc->t_par[1]));
399 }
400 INST(template <typename T> class Vector friend void job_val_sub_vec (struct thr_ctrl *);)
402 template <typename T>
403 void job_val_sub_vec (struct thr_ctrl *tc)
404 {
405  do_val_sub_vec(tc->t_size, (T*)(tc->t_par[0]),
406  *(const T*)(tc->t_par[1]));
407 }
408 INST(template <typename T> class Vector friend void job_vec_mul_val (struct thr_ctrl *);)
410 template <typename T>
411 void job_vec_mul_val (struct thr_ctrl *tc)
412 {
413  do_vec_mul_val (tc->t_size, (T*)(tc->t_par[0]),
414  *(const T*)(tc->t_par[1]));
415 }
416 INST(template <typename T> class Vector friend void job_vec_div_val (struct thr_ctrl *);)
418 template <typename T>
419 void job_vec_div_val (struct thr_ctrl *tc)
420 {
421  do_vec_div_val (tc->t_size, (T*)(tc->t_par[0]),
422  *(const T*)(tc->t_par[1]));
423 }
424 INST(template <typename T> class Vector friend void job_val_div_vec (struct thr_ctrl *);)
426 template <typename T>
427 void job_val_div_vec (struct thr_ctrl *tc)
428 {
429  do_val_div_vec(tc->t_size, (T*)(tc->t_par[0]),
430  *(const T*)(tc->t_par[1]));
431 }
432 
433 /* TSV variants */
434 INST(template <typename T> class Vector friend void job_vec_svc_add (struct thr_ctrl *);)
436 template <typename T>
437 void job_vec_svc_add (struct thr_ctrl *tc)
438 {
439  do_vec_svc_add (tc->t_size, (T*)(tc->t_par[0]),
440  (const T*)(tc->t_par[1]), (const T*)(tc->t_par[2]),
441  *(const T*)(tc->t_par[3]));
442 }
443 INST(template <typename T> class Vector friend void job_svc_vec_add (struct thr_ctrl *);)
445 template <typename T>
446 void job_svc_vec_add (struct thr_ctrl *tc)
447 {
448  do_svc_vec_add (tc->t_size, (T*)(tc->t_par[0]),
449  (const T*)(tc->t_par[1]), (const T*)(tc->t_par[2]),
450  *(const T*)(tc->t_par[3]));
451 }
452 INST(template <typename T> class Vector friend void job_svc_svc_add (struct thr_ctrl *);)
454 template <typename T>
455 void job_svc_svc_add (struct thr_ctrl *tc)
456 {
457  do_svc_svc_add (tc->t_size, (T*)(tc->t_par[0]),
458  (const T*)(tc->t_par[1]), (const T*)(tc->t_par[2]),
459  *(const T*)(tc->t_par[3]), *(const T*)(tc->t_par[4]));
460 }
462 INST(template <typename T> class Vector friend void job_vec_svc_sub (struct thr_ctrl *);)
463 template <typename T>
464 void job_vec_svc_sub (struct thr_ctrl *tc)
465 {
466  do_vec_svc_sub (tc->t_size, (T*)(tc->t_par[0]),
467  (const T*)(tc->t_par[1]), (const T*)(tc->t_par[2]),
468  *(const T*)(tc->t_par[3]));
469 }
470 INST(template <typename T> class Vector friend void job_svc_vec_sub (struct thr_ctrl *);)
472 template <typename T>
473 void job_svc_vec_sub (struct thr_ctrl *tc)
474 {
475  do_svc_vec_sub (tc->t_size, (T*)(tc->t_par[0]),
476  (const T*)(tc->t_par[1]), (const T*)(tc->t_par[2]),
477  *(const T*)(tc->t_par[3]));
478 }
479 INST(template <typename T> class Vector friend void job_svc_svc_sub (struct thr_ctrl *);)
481 template <typename T>
482 void job_svc_svc_sub (struct thr_ctrl *tc)
483 {
484  do_svc_svc_sub (tc->t_size, (T*)(tc->t_par[0]),
485  (const T*)(tc->t_par[1]), (const T*)(tc->t_par[2]),
486  *(const T*)(tc->t_par[3]), *(const T*)(tc->t_par[4]));
487 }
488 INST(template <typename T> class Vector friend void job_vec_add_svc (struct thr_ctrl *);)
490 template <typename T>
491 void job_vec_add_svc (struct thr_ctrl *tc)
492 {
493  do_vec_add_svc (tc->t_size, (T*)(tc->t_par[0]),
494  (const T*)(tc->t_par[1]),
495  *(const T*)(tc->t_par[2]));
496 }
497 INST(template <typename T> class Vector friend void job_vec_sub_svc (struct thr_ctrl *);)
499 template <typename T>
500 void job_vec_sub_svc (struct thr_ctrl *tc)
501 {
502  do_vec_sub_svc (tc->t_size, (T*)(tc->t_par[0]),
503  (const T*)(tc->t_par[1]),
504  *(const T*)(tc->t_par[2]));
505 }
506 INST(template <typename T> class Vector friend void job_vec_sub_svc_inv (struct thr_ctrl *);)
508 template <typename T>
509 void job_vec_sub_svc_inv (struct thr_ctrl *tc)
510 {
511  do_vec_sub_svc_inv (tc->t_size, (T*)(tc->t_par[0]),
512  (const T*)(tc->t_par[1]),
513  *(const T*)(tc->t_par[2]));
514 }
515 
516 INST(template <typename T> class Vector friend void job_svc_val_add (struct thr_ctrl *);)
518 template <typename T>
519 void job_svc_val_add (struct thr_ctrl *tc)
520 {
521  do_svc_val_add (tc->t_size, (T*)(tc->t_par[0]),
522  (const T*)(tc->t_par[1]),
523  *(const T*)(tc->t_par[2]), *(const T*)(tc->t_par[3]));
524 }
525 INST(template <typename T> class Vector friend void job_svc_val_sub (struct thr_ctrl *);)
527 template <typename T>
528 void job_svc_val_sub (struct thr_ctrl *tc)
529 {
530  do_svc_val_sub (tc->t_size, (T*)(tc->t_par[0]),
531  (const T*)(tc->t_par[1]),
532  *(const T*)(tc->t_par[2]), *(const T*)(tc->t_par[3]));
533 }
534 INST(template <typename T> class Vector friend void job_val_svc_add (struct thr_ctrl *);)
536 template <typename T>
537 void job_val_svc_add (struct thr_ctrl *tc)
538 {
539  do_val_svc_add (tc->t_size, (T*)(tc->t_par[0]),
540  (const T*)(tc->t_par[1]),
541  *(const T*)(tc->t_par[2]), *(const T*)(tc->t_par[3]));
542 }
543 INST(template <typename T> class Vector friend void job_val_svc_sub (struct thr_ctrl *);)
545 template <typename T>
546 void job_val_svc_sub (struct thr_ctrl *tc)
547 {
548  do_val_svc_sub (tc->t_size, (T*)(tc->t_par[0]),
549  (const T*)(tc->t_par[1]),
550  *(const T*)(tc->t_par[2]), *(const T*)(tc->t_par[3]));
551 }
552 INST(template <typename T> class Vector friend void job_val_svc_div (struct thr_ctrl *);)
554 template <typename T>
555 void job_val_svc_div (struct thr_ctrl *tc)
556 {
557  do_val_svc_div (tc->t_size, (T*)(tc->t_par[0]),
558  (const T*)(tc->t_par[1]),
559  *(const T*)(tc->t_par[2]), *(const T*)(tc->t_par[3]));
560 }
561 
562 
563 /* TODO: conj, unary - */
564 
565 
566 #endif
567 
568 /* Naming rules:
569  * 1. first arg is always pointer
570  * 2. we've max 3 pointers and two consts
571  * 3. the last two args are designated by letters, C = const, V = vector
572  * 4. letters that are clear because of rules 2 or 3 are omitted
573  * 5. pointer(vector) args are always before consts
574  */
575 #if !defined(SMP) || defined(NOSMP_VECVEC)
576 #define STD_SMP_TEMPLATE2V(oper, dm, a1, a2) \
577  do_##oper <T> (dm, a1, a2)
578 #define STD_SMP_TEMPLATE2C(oper, dm, a1, a2) \
579  do_##oper <T> (dm, a1, a2)
580 #define STD_SMP_TEMPLATE3VV(oper, dm, a1, a2, a3) \
581  do_##oper <T> (dm, a1, a2, a3)
582 #define STD_SMP_TEMPLATE3VC(oper, dm, a1, a2, a3) \
583  do_##oper <T> (dm, a1, a2, a3)
584 #define STD_SMP_TEMPLATE3CC(oper, dm, a1, a2, a3) \
585  do_##oper <T> (dm, a1, a2, a3)
586 #define STD_SMP_TEMPLATE4V(oper, dm, a1, a2, a3, a4) \
587  do_##oper <T> (dm, a1, a2, a3, a4)
588 #define STD_SMP_TEMPLATE4C(oper, dm, a1, a2, a3, a4) \
589  do_##oper <T> (dm, a1, a2, a3, a4)
590 #define STD_SMP_TEMPLATE5(oper, dm, a1, a2, a3, a4, a5) \
591  do_##oper <T> (dm, a1, a2, a3, a4, a5)
592 
593 #else /* SMP */
594 
595 #define _SMP_TMPL2V(oper,st,en,a1,a2,a3,a4,a5) do_##oper <T> (en-st, a1+st, a2+st)
596 #define _JOB_TMPL2V(oper,st,en,a1,a2,a3,a4,a5) job_##oper <T>, en-st, a1+st, a2+st
597 #define _SMP_TMPL2C(oper,st,en,a1,a2,a3,a4,a5) do_##oper <T> (en-st, a1+st, a2)
598 #define _JOB_TMPL2C(oper,st,en,a1,a2,a3,a4,a5) job_##oper <T>, en-st, a1+st, &a2
599 #define _SMP_TMPL3VV(oper,st,en,a1,a2,a3,a4,a5) do_##oper <T> (en-st, a1+st, a2+st, a3+st)
600 #define _JOB_TMPL3VV(oper,st,en,a1,a2,a3,a4,a5) job_##oper <T>, en-st, a1+st, a2+st, a3+st
601 #define _SMP_TMPL3VC(oper,st,en,a1,a2,a3,a4,a5) do_##oper <T> (en-st, a1+st, a2+st, a3)
602 #define _JOB_TMPL3VC(oper,st,en,a1,a2,a3,a4,a5) job_##oper <T>, en-st, a1+st, a2+st, &a3
603 #define _SMP_TMPL3CC(oper,st,en,a1,a2,a3,a4,a5) do_##oper <T> (en-st, a1+st, a2, a3)
604 #define _JOB_TMPL3CC(oper,st,en,a1,a2,a3,a4,a5) job_##oper <T>, en-st, a1+st, &a2, &a3
605 #define _SMP_TMPL4V(oper,st,en,a1,a2,a3,a4,a5) do_##oper <T> (en-st, a1+st, a2+st, a3+st, a4)
606 #define _JOB_TMPL4V(oper,st,en,a1,a2,a3,a4,a5) job_##oper <T>, en-st, a1+st, a2+st, a3+st, &a4
607 #define _SMP_TMPL4C(oper,st,en,a1,a2,a3,a4,a5) do_##oper <T> (en-st, a1+st, a2+st, a3, a4)
608 #define _JOB_TMPL4C(oper,st,en,a1,a2,a3,a4,a5) job_##oper <T>, en-st, a1+st, a2+st, &a3, &a4
609 #define _SMP_TMPL5(oper,st,en,a1,a2,a3,a4,a5) do_##oper <T> (en-st, a1+st, a2+st, a3+st, a4, a5)
610 #define _JOB_TMPL5(oper,st,en,a1,a2,a3,a4,a5) job_##oper <T>, en-st, a1+st, a2+st, a3+st, &a4, &a5
611 
617 #define STD_SMP_TEMPLATE(oper, dm, a1, a2, a3, a4, a5, SMP_TMPL, JOB_TMPL, p2, p3) \
618  /* use some heuristic to decide for the num of threads */ \
619  const unsigned n_thr = threads_avail (dm / SMP_VECSLICE2); \
620  update_n_thr(n_thr); \
621  if (LIKELY(n_thr < 2)) { \
622  SMP_TMPL(oper, 0, dm, a1, a2, a3, a4, a5); \
623  } else { \
624  PREFETCH_W_MANY(a1, 3); \
625  p2 \
626  p3 \
627  const unsigned long first = slice_offset(2, n_thr, dm, a1); \
628  unsigned long st, en = first; \
629  unsigned _t; \
630  /* smp_barrier(); */ \
631  /* Start threads */ \
632  for (_t = 0; _t < n_thr-1; ++_t) { \
633  st = en; en = slice_offset(_t+2, n_thr, dm, a1); \
634  thread_start (_t, JOB_TMPL(oper, st, en, a1, a2, a3, a4, a5), (void*)0); \
635  } \
636  /* The first slice is handled by the main thread */ \
637  SMP_TMPL(oper, 0UL, first, a1, a2, a3, a4, a5); \
638  /* sched_yield (); */ \
639  /* Wait for the end */ \
640  for (_t = 0; _t < n_thr-1; ++_t) \
641  thread_wait (_t); \
642  }
643 
644 #define STD_SMP_TEMPLATE2V(oper, dm, a1, a2) \
645  STD_SMP_TEMPLATE(oper, dm, a1, a2, NULL, NULL, NULL, _SMP_TMPL2V, _JOB_TMPL2V, PREFETCH_R_MANY(a2,2);,)
646 #define STD_SMP_TEMPLATE2C(oper, dm, a1, a2) \
647  STD_SMP_TEMPLATE(oper, dm, a1, a2, NULL, NULL, NULL, _SMP_TMPL2C, _JOB_TMPL2C,,)
648 #define STD_SMP_TEMPLATE3VV(oper, dm, a1, a2, a3) \
649  STD_SMP_TEMPLATE(oper, dm, a1, a2, a3, NULL, NULL, _SMP_TMPL3VV, _JOB_TMPL3VV, PREFETCH_R_MANY(a2,2);, PREFETCH_R_MANY(a3,2);)
650 #define STD_SMP_TEMPLATE3VC(oper, dm, a1, a2, a3) \
651  STD_SMP_TEMPLATE(oper, dm, a1, a2, a3, NULL, NULL, _SMP_TMPL3VC, _JOB_TMPL3VC, PREFETCH_R_MANY(a2,2);,)
652 #define STD_SMP_TEMPLATE3CC(oper, dm, a1, a2, a3) \
653  STD_SMP_TEMPLATE(oper, dm, a1, a2, a3, NULL, NULL, _SMP_TMPL3CC, _JOB_TMPL3CC,,)
654 #define STD_SMP_TEMPLATE4V(oper, dm, a1, a2, a3, a4) \
655  STD_SMP_TEMPLATE(oper, dm, a1, a2, a3, a4, NULL, _SMP_TMPL4V, _JOB_TMPL4V, PREFETCH_R_MANY(a2,2);,)
656 #define STD_SMP_TEMPLATE4C(oper, dm, a1, a2, a3, a4) \
657  STD_SMP_TEMPLATE(oper, dm, a1, a2, a3, a4, NULL, _SMP_TMPL4C, _JOB_TMPL4C, PREFETCH_R_MANY(a2,2);,)
658 #define STD_SMP_TEMPLATE5(oper, dm, a1, a2, a3, a4, a5) \
659  STD_SMP_TEMPLATE(oper, dm, a1, a2, a3, a4, a5, _SMP_TMPL5, _JOB_TMPL5, PREFETCH_R_MANY(a2,2);, PREFETCH_R_MANY(a3,2);)
660 
661 
662 #endif /* SMP */
663 
665 template <typename T>
666 inline TVector<T>& TVector<T>::operator += (const T& a)
667 {
668  STD_SMP_TEMPLATE2C(vec_add_val, this->dim, this->vec, a);
669  return *this;
670 }
672 template <typename T>
673 inline TVector<T>& TVector<T>::operator -= (const T& a)
674 {
675  STD_SMP_TEMPLATE2C(vec_sub_val, this->dim, this->vec, a);
676  return *this;
677 }
679 template <typename T>
680 inline TVector<T>& TVector<T>::operator *= (const T& a)
681 {
682  STD_SMP_TEMPLATE2C(vec_mul_val, this->dim, this->vec, a);
683  return *this;
684 }
685 
687 template <typename T>
688 inline TVector<T>& TVector<T>::operator /= (const T& a)
689 {
690  BCHK (a == (T)0, VecErr, Divide by zero (TVector), 0, *this);
691  return this->operator *= ((T)1/a);
692 }
693 
694 
696 template <typename T>
697 inline TVector<T>& TVector<T>::operator += (const Vector<T>& a)
698 {
699  BCHK(this->dim!=a.dim,VecErr,"TV += V operator with diff. sizes", a.dim, *this);
700  STD_SMP_TEMPLATE2V(vec_add_vec, this->dim, this->vec, a.vec);
701  return *this;
702 }
704 template <typename T>
705 inline TVector<T>& TVector<T>::operator -= (const Vector<T>& a)
706 {
707  BCHK(this->dim!=a.dim,VecErr,"TV -= V operator with diff. sizes", a.dim, *this);
708  STD_SMP_TEMPLATE2V(vec_sub_vec, this->dim, this->vec, a.vec);
709  return *this;
710 }
711 
713 template <typename T>
714 inline TVector<T>& TVector<T>::operator += (const TSVector<T>& a)
715 {
716  BCHK(this->dim!=a.dim,VecErr,"TV += TSV operator with diff. sizes", a.dim, *this);
717  STD_SMP_TEMPLATE3VC(vec_add_svc, this->dim, this->vec, a.vec, a.fac);
718  if (LIKELY(this->vec != a.vec))
719  a.destroy ();
720  return *this;
721 }
723 template <typename T>
724 inline TVector<T>& TVector<T>::operator -= (const TSVector<T>& a)
725 {
726  BCHK(this->dim!=a.dim,VecErr,"TV -= TSV operator with diff. sizes", a.dim, *this);
727  STD_SMP_TEMPLATE3VC(vec_sub_svc, this->dim, this->vec, a.vec, a.fac);
728  if (LIKELY(this->vec != a.vec))
729  a.destroy ();
730  return *this;
731 }
732 
733 
734 // Binary operations (for TV id. to unary)
735 
736 template <typename T>
737 inline TVector<T>& TVector<T>::operator + (const Vector<T>& a)
738 {
739  //BCHK(dim != a.dim, VecErr, Different sizes in operator TV + V, a.dim, *this);
740  return this->operator += (a);
741 }
742 template <typename T>
743 inline TVector<T>& TVector<T>::operator - (const Vector<T>& a)
744 {
745  //BCHK(dim != a.dim, VecErr, Different sizes in operator TV - V, a.dim, *this);
746  return this->operator -= (a);
747 }
748 
749 
752 template <typename T>
753 const TVector<T>& TVector<T>::operator + (const TVector<T>& tv) /*const*/
754 {
755  BCHK (tv.dim != this->dim, VecErr, Diff. dim in TV + TV, tv.dim, (this->destroy(), tv));
756  //for (REGISTER unsigned long i = 0; i < dim; ++i) tv.vec[i] = vec[i] + tv.vec[i];
757  STD_SMP_TEMPLATE2V(vec_add_vec, this->dim, tv.vec, this->vec);
758  this->destroy (); return tv;
759 }
760 
763 template <typename T>
764 const TVector<T>& TVector<T>::operator - (const TVector<T>& tv) /*const*/
765 {
766  BCHK (tv.dim != this->dim, VecErr, Diff. dim in TV - TV, tv.dim, (this->destroy(), tv));
767  //for (REGISTER unsigned long i = 0; i < dim; ++i) tv.vec[i] = vec[i] - tv.vec[i];
768  STD_SMP_TEMPLATE2V(vec_sub_vec_inv, this->dim, tv.vec, this->vec);
769  this->destroy (); return tv;
770 }
771 
772 
775 INST(template<typename T> class TVector friend const TVector<T>& operator + (const T&, const TVector<T>&);)
776 template <typename T>
777 inline const TVector<T>& operator + (const T& a, const TVector<T>& b)
778 {
779  //for (REGISTER unsigned long i=0; i<b.size(); ++i) b.set(a + b.get(i), i);
780  const unsigned long dim = b.size();
781  T* vec = &b.setval(0);
782  STD_SMP_TEMPLATE2C(val_add_vec, dim, vec, a);
783  return b;
784 }
787 INST(template<typename T> class TVector friend const TVector<T>& operator - (const T&, const TVector<T>&);)
788 template <typename T>
789 inline const TVector<T>& operator - (const T& a, const TVector<T>& b)
790 {
791  //for (REGISTER unsigned long i=0; i<b.size(); ++i) b.set(a - b.get(i), i);
792  const unsigned long dim = b.size();
793  T* vec = &b.setval(0);
794  STD_SMP_TEMPLATE2C(val_sub_vec, dim, vec, a);
795  return b;
796 }
797 
798 
799 INST(template <typename T> class TVector friend TSVector<T> operator * (const T&, const TVector<T>&);)
800 template <typename T>
801 inline TSVector<T> operator * (const T& a, const TVector<T>& b)
802 {
803  return TSVector<T> (b, a);
804 }
805 
806 // a/b = a*b/b^2
807 // We have to take care T not to be an int type!
808 INSTCTL(#ifndef VEC_INT_DIV_SUPPORT)
809 INST(template <typename T> class TVector friend TSVector<T> operator / (const T&, TVector<T>);)
810 INSTCTL(#endif)
811 template <typename T>
812 inline TSVector<T> operator / (const T& a, TVector<T> b)
813 {
814  T dv = a / b.fabssqr();
815  BCHK(dv==(T)0, VecErr, Divide by norm 0 Vector, 0, 0);
816  return TSVector<T> (b, dv);
817 }
818 
819 
820 template <typename T>
821 inline TVector<T>& TVector<T>::operator + (const T& a)
822 { return this->operator += (a); }
823 template <typename T>
824 inline TVector<T>& TVector<T>::operator - (const T& a)
825 { return this->operator -= (a); }
826 
827 
828 template <typename T>
829 inline TSVector<T> TVector<T>::operator * (const T& a)
830 {
831  return TSVector<T> (*this, a);
832 }
833 
834 // We have to take care T not to be an int type!
835 template <typename T>
836 inline TSVector<T> TVector<T>::operator / (const T& a)
837 {
838  return TSVector<T> (*this, (T)1/a);
839 }
840 
841 #if 0
842 template <typename T>
843 inline TVector<T>& TVector<T>::operator - ()
844 {
845  for (REGISTER unsigned long i=0; i<dim; ++i)
846  this->vec[i] = -this->vec[i];
847  return *this;
848 }
849 #else
850 template <typename T>
851 inline TSVector<T> TVector<T>::operator - ()
852 {
853  return TSVector<T> (*this, (T)-1);
854 }
855 #endif
856 
857 //NAMESPACE_END
858 //NAMESPACE_CPLX
859 
860 //INST(template <typename T> class TVector friend TVector<T>& conj (TVector<T>&);)
861 template <typename T>
862 inline TVector<T>& conj (TVector<T>& tv)
863 {
864  for (REGISTER unsigned long i=0; i<tv.size(); ++i)
865  tv.set (CPLX__ conj (tv.get(i)), i);
866  return tv;
867 }
868 
869 //INST(template <typename T> class TVector friend TVector<T>& real (TVector<T>&);)
870 template <typename T>
871 inline TVector<T>& real (TVector<T>& tv)
872 {
873  for (REGISTER unsigned long i=0; i<tv.size(); ++i)
874  tv.set(CPLX__ real (tv.get(i)), i);
875  //tv.vec[i] = (tv.vec[i] + CPLX__ conj (tv.vec[i])) * 0.5;
876  return tv;
877 }
878 
879 //INST(template <typename T> class TVector friend TVector<T>& imag (TVector<T>&);)
880 template <typename T>
881 inline TVector<T>& imag (TVector<T>& tv)
882 {
883  for (REGISTER unsigned long i=0; i<tv.size(); ++i)
884  tv.set(CPLX__ imag (tv.get(i)), i);
885  //tv.vec[i] = (tv.vec[i] - CPLX__ conj (tv.vec[i])) * 0.5;
886  return tv;
887 }
888 
889 //INST(template <typename T> class TVector friend TVector<T> conj (const Vector<T>&);)
890 template <typename T>
891 inline TVector<T> conj (const Vector<T>& v)
892 {
893  TVector<T> tv(v.size());
894  for (REGISTER unsigned long i=0; i<v.size(); ++i)
895  tv.set (CPLX__ conj (v(i)), i);
896  return tv;
897 }
898 
899 //INST(template <typename T> class TVector friend TVector<T> real (const Vector<T>&);)
900 template <typename T>
901 inline TVector<T> real (const Vector<T>& v)
902 {
903  TVector<T> tv (v.size());
904  for (REGISTER unsigned long i=0; i<v.size(); ++i)
905  tv.set (CPLX__ real (v(i)), i);
906  //tv.vec[i] = (v.vec[i] + CPLX__ conj (v.vec[i])) * 0.5;
907  return tv;
908 }
909 
910 //INST(template <typename T> class TVector friend TVector<T> imag (const Vector<T>&);)
911 template <typename T>
912 inline TVector<T> imag (const Vector<T>& v)
913 {
914  TVector<T> tv (v.size());
915  for (REGISTER unsigned long i=0; i<v.size(); ++i)
916  tv.set(CPLX__ imag (v(i)), i);
917  //tv.vec[i] = (v.vec[i] - CPLX__ conj (v.vec[i])) * 0.5;
918  return tv;
919 }
920 
921 //NAMESPACE_CPLX_END
922 //NAMESPACE_TBCI
923 
926 template <typename T>
927 inline TVector<T>& TVector<T>::operator + (const TSVector<T>& ts)
928 {
929  BCHK(this->dim!=ts.dim,VecErr,Wrong dim TV + TSV,ts.dim,*this);
930  //for (REGISTER unsigned long t = 0; t < dim; t++)
931  // vec[t] op##= ts.fac * ts.vec[t];
932  STD_SMP_TEMPLATE3VC(vec_add_svc, this->dim, this->vec, ts.vec, ts.fac);
933  /* included by AA */
934  if (LIKELY(this->vec != ts.vec))
935  ts.destroy();
936  return *this;
937 }
938 
941 template <typename T>
942 inline TVector<T>& TVector<T>::operator - (const TSVector<T>& ts)
943 {
944  BCHK(this->dim!=ts.dim,VecErr,Wrong dim TV - TSV,ts.dim,*this);
945  //for (REGISTER unsigned long t = 0; t < dim; t++)
946  // vec[t] op##= ts.fac * ts.vec[t];
947  STD_SMP_TEMPLATE3VC(vec_sub_svc, this->dim, this->vec, ts.vec, ts.fac);
948  /* included by AA */
949  if (LIKELY(this->vec != ts.vec))
950  ts.destroy();
951  return *this;
952 }
953 
954 
955 template <typename T>
956 inline TVector<T> TVector<T>::slice (const unsigned long i0, const unsigned long i1)
957 {
958  BCHK (i1 < i0, VecErr, Slicing needs ordered arguments, i0, *this);
959  BCHK (i0 >= this->dim, VecErr, Slice: arg1 out of range, i0, *this);
960  BCHK (i1 > this->dim, VecErr, Slice: arg2 out of range, i1, *this);
961  this->dim = i1 - i0;
962  // BEWARE: This relies on the C library freeing the right chunk of
963  // memory, even if it moved within the allocated region!!!
964  if (UNLIKELY(this->dim == 0))
965  TBCIDELETE (T, this->vec, this->dim);
966  else
967  this->vec += i0;
968  return *this;
969 }
970 
971 
972 template <typename T>
973 inline TVector<T>& TVector<T>::incr (const unsigned long wh, const T i)
974 {
975  BCHK (wh >= this->dim, VecErr, Tried to change not-exitsting idx, wh, *this);
976  this->vec[wh] += i;
977  return *this;
978 }
979 
980 //fwd-decl
981 //template <typename T> inline double fabssqr(TVector<T> tv);
982 template <typename T> inline double fabssqr(const Vector<T>& v);
983 
984 NAMESPACE_END
985 
986 NAMESPACE_CSTD
987 
988 INST(template <typename T> class TBCI__ Vector friend double fabs (const TBCI__ Vector<T>&);)
989 template <typename T>
990 inline double fabs (const TBCI__ Vector<T>& v)
991 {
992  return sqrt (TBCI__ fabssqr (v));
993 }
994 
995 INST(template <typename T> class TBCI__ Vector friend T abs (const TBCI__ Vector<T>&);)
996 template <typename T>
997 inline T abs (const TBCI__ Vector<T>& v)
998 {
999  return v.abs();
1000 }
1001 
1002 INST(template <typename T> class TBCI__ TVector friend double fabs (TBCI__ TVector<T>);)
1003 template <typename T>
1004 inline double fabs (TBCI__ TVector<T> tv)
1005 {
1006  TBCI__ Vector<T> v(tv); return fabs (v);
1007 }
1008 
1009 INST(template <typename T> class TBCI__ Vector friend T abs (TBCI__ TVector<T>);)
1010 template <typename T>
1011 inline T abs (TBCI__ TVector<T> tv)
1012 {
1013  TBCI__ Vector<T> v(tv); return abs (v);
1014 }
1015 
1016 NAMESPACE_CSTD_END
1017 
1018 NAMESPACE_TBCI
1019 
1025 template <typename T>
1026 class TSVector : public Vector_Sig<T> //: public TVector<T>
1027 {
1028  protected:
1029  mutable T* vec; mutable unsigned long dim; //bool keep;
1030  mutable T fac ALIGN(MIN_ALIGN2);
1031 
1032  void clone (const bool = false, const T* = 0) const;
1033  public:
1034  void detach (const T* = 0) const;
1035 
1036  public:
1037  mutable bool mut; // mutability (ie private data)
1038  void destroy () const;
1039 
1040  public:
1041  friend class Vector<T>;
1042  friend class TVector<T>;
1043  friend class CRMatrix<T>;
1044  friend NOINST void FRIEND_TBCI2__ do_mat_tsv_mult FGD (const unsigned start, const unsigned end, \
1045  TVector<T> * res, const Matrix<T> * mat, const TSVector<T> * rsv);
1046 
1047  typedef T value_type;
1048  typedef T element_type;
1049  typedef T aligned_value_type TALIGN(MIN_ALIGN2);
1050 
1051  const TSVector<T>& eval (const T* vv = 0) const;
1052 
1053  TSVector () : vec((T*)0), dim(0), fac((T)1), mut(true) {}
1054  ~TSVector () { /* keep = true; */ }
1055  explicit TSVector (const unsigned long d) : vec(new T[d]), dim(d), fac((T)1), mut(true) {}
1056 
1057  TSVector (const TSVector<T>& ts)
1058  : vec(ts.vec), dim(ts.dim), fac(ts.fac), mut(ts.mut) {}
1059 
1060  TSVector (const TVector<T>& tv, const T& f = (T)1)
1061  : vec(tv.vec), dim(tv.dim), fac(f), mut(true) {}
1062  TSVector (const Vector<T>& v, const T& f = (T)1)
1063  : vec(v.vec), dim(v.dim), fac(f), mut(false) {}
1064 
1065 
1066  T operator () (const unsigned long i) const
1067  { T val = vec[i] * fac; destroy (); return val; }
1068  //T& operator () (const unsigned long i) { STD__ cerr << "TSVectors may not be written to by () op.\n"; abort(); }
1069  T operator [] (const unsigned long i) const { return this->operator() (i); }
1070  //T& operator [] (const unsigned long i) { return this->operator() (i); }
1071 
1072  T get (const unsigned long i) const HOT { return vec[i] * fac; }
1073  const T& getfac () const { return fac; }
1074  T& getfacref () const { return fac; }
1075  const T& getcref (const unsigned long i) const HOT { return vec[i]; }
1076 
1077  // alias equ ops
1078  TSVector<T>& operator = (const TSVector<T>& ts)
1079  { destroy (); vec = ts.vec; dim = ts.dim;
1080  fac = ts.fac; mut = ts.mut; return *this; }
1081  TSVector<T>& operator = (const TVector<T>& tv)
1082  { destroy (); vec = tv.vec; dim = tv.dim;
1083  fac = (T)1; mut = true; return *this; }
1084  /*
1085  TSVector<T>& operator = (const Vector<T>& v)
1086  { destroy (); vec = v.vec; dim = v.dim;
1087  fac = (T)1; mut = false; return *this; }
1088  */
1089 
1090  // first the easy to implement operators
1091  TSVector<T>& operator *= (const T& f) { fac *= f; return *this; }
1092  TSVector<T>& operator /= (const T& f) { fac /= f; return *this; }
1093  const TSVector<T>& operator * (const T& f) const { fac *= f; return *this; }
1094  const TSVector<T>& operator / (const T& f) const { fac /= f; return *this; }
1095 
1096  const TSVector<T>& operator - () const { fac = -fac; return *this; }
1097 
1098  // now for the more complicated ones
1099  TVector<T> operator + (const Vector<T>&) const HOT;
1100  TVector<T> operator + (TVector<T>) const HOT;
1101  TVector<T> operator + (const TSVector<T>&) const HOT;
1102  TVector<T> operator + (const T&) const;
1103  TVector<T> operator - (const Vector<T>&) const HOT;
1104  TVector<T> operator - (TVector<T>) const HOT;
1105  TVector<T> operator - (const TSVector<T>&) const HOT;
1106  TVector<T> operator - (const T&) const;
1107 
1108 #ifdef TSV_OP_FRIENDS
1109  friend const TSVector<T>& operator * FGD (const T& f, const TSVector<T>& ts);
1110  friend TVector<T> operator + FGD (const T&, const TSVector<T>&);
1111  friend TVector<T> operator - FGD (const T&, const TSVector<T>&);
1112  // for those "friends"
1113 #else
1114  TVector<T> add_t_tsv (const T&) const;
1115  TVector<T> sub_t_tsv (const T&) const;
1116 #endif
1117 
1118  TVector<T> incr (const unsigned long, const T = (T)1) const;
1119 
1120  inline unsigned long size () const { return dim; }
1121 
1122  T min () { Vector<T> th(*this); return th.min (); }
1123  T max () { Vector<T> th(*this); return th.max (); }
1124  T abs () const /* NOTE: The const IS wrong */
1125  { return ((TVector<T>)*this).abs () * GLBL__ CSTD__ abs(fac); }
1126  double fabs () const
1127  { return ((TVector<T>)*this).fabs () * GLBL__ MATH__ fabs(fac); }
1128  double fabssqr () const
1129  { return ((TVector<T>)*this).fabssqr () * GLBL2__ TBCI__ fabssqr(fac); }
1130 
1131  T sum () { T val = ((TVector<T>)*this).sum () * fac; destroy (); return val; }
1132 
1133  /* Why are those declared but not defined ? */
1134 #ifdef TSVEC_NEED_DOT
1135  T operator * (const TSVector<T>&);
1136  T operator * (const TVector<T>&);
1137  T operator * (const Vector<T>&);
1138 #endif
1139 
1140 #if 0 /* not implemented! */
1141  friend TVector<T> conj FGD (const TSVector<T>&);
1142  friend TVector<T> real FGD (const TSVector<T>&);
1143  friend TVector<T> imag FGD (const TSVector<T>&);
1144 #endif
1145 
1146 // They fortunately also don't need to be friends
1147 #ifndef HAVE_GCC295_FRIEND_BUG
1148  friend NOINST double FRIEND_TBCI__ fabssqr FGDT (const TSVector<T>&);
1149 #endif
1150 #if 0
1151  friend NOINST double MATH__ fabs FGD (const TBCI__ TSVector<T>&);
1152  friend NOINST T CSTD__ abs FGD (const TBCI__ TSVector<T>&);
1153 #endif
1154 
1155  bool operator == (const Vector<T>&) const;
1156  bool operator != (const Vector<T>& v) const { return !(*this == v); }
1157  bool operator == (const TSVector<T>&) const;
1158  bool operator != (const TSVector<T>& ts) const { return !(*this == ts); }
1159  bool operator == (const TVector<T>& tv) const { Vector<T> v(tv); return (*this == v); }
1160  bool operator != (const TVector<T>& tv) const { return !(*this == tv); }
1161 
1162  // Needed for emul et al.
1163  T* const & vecptr () const { return vec; }
1164 
1165  static const char* vec_info() { return "TSVector"; }
1166 
1167  friend STD__ ostream& operator << FGD (STD__ ostream&, const TSVector<T>&);
1168 } ;
1169 
1170 
1171 template <typename T>
1172 inline void TSVector<T>::destroy () const
1173 {
1174  if (LIKELY(!mut || !dim)) return;
1175  TBCIDELETE(T, vec, dim);
1176 }
1177 
1178 template <typename T>
1179 inline void TSVector<T>::detach (const T* vv) const
1180 {
1181  if (LIKELY(!dim || (mut && !vv)))
1182  return;
1183  mut = true;
1184  if (vv)
1185  vec = const_cast<T*> (vv);
1186  else {
1187  vec = NEW(T,dim);
1188  if (UNLIKELY(!vec)) {
1189  dim = 0; mut = false;
1190  }
1191  }
1192 }
1193 
1194 INST(template <typename T> class TSVector<T> friend const TSVector<T>& operator * (const T&, const TSVector<T>&);)
1195 template <typename T>
1196 const TSVector<T>& operator * (const T& f, const TSVector<T>& ts)
1197 { ts.getfacref() = f * ts.getfac(); return ts; }
1198 
1199 template <typename T>
1200 inline void TSVector<T>::clone (const bool evl, const T* vv) const
1201 {
1202  if (LIKELY(!dim || (mut && !vv)))
1203  return;
1204  const T* oldv = vec; bool oldm = mut; detach (vv);
1205  if (LIKELY(evl && fac != (T)1)) {
1206 #if 0 //def __i386__
1207  for (REGISTER unsigned long t = 0; t < dim; t++)
1208  vec[t] = oldv[t] * fac;
1209 #else
1210  STD_SMP_TEMPLATE3VC(vec_val_mul, dim, vec, oldv, fac);
1211 #endif
1212  fac = (T)1;
1213  }
1214  else
1215  TBCICOPY(vec, oldv, T, dim);
1216  if (UNLIKELY(vv && oldm))
1217  TBCIDELETE(T, oldv, dim);
1218 }
1219 
1220 template <typename T>
1221 const TSVector<T>& TSVector<T>::eval (const T* vv) const
1222 {
1223  if (!vv && mut) {
1224  if (fac != (T)1) {
1225 #if 0 //def __i386__
1226  for (REGISTER unsigned long t = 0; t < dim; t++)
1227  vec[t] *= fac;
1228 #else
1229  STD_SMP_TEMPLATE2C(vec_mul_val, dim, vec, fac);
1230 #endif
1231  fac = (T)1;
1232  }
1233  // else do nothing
1234  }
1235  else
1236  clone (true, vv);
1237  return *this;
1238 }
1239 
1240 template <typename T>
1241 bool TSVector<T>::operator == (const Vector<T>& v) const
1242 {
1243  if (LIKELY(dim != v.dim)) {
1244  destroy (); return false;
1245  }
1246  if (LIKELY(dim == 0)) {
1247  /*destroy();*/ return true;
1248  }
1249  if (LIKELY(fac == (T)1)) {
1250  if (LIKELY(vec == v.vec))
1251  return true;
1252  if (TBCICOMP(vec, v.vec, T, dim)) {
1253  destroy (); return false;
1254  }
1255  } else {
1256  if (LIKELY(vec == v.vec))
1257  return false;
1258  for (REGISTER unsigned long i = 0; i < dim; ++i)
1259  if (UNLIKELY(vec[i]*fac != v.vec[i])) {
1260  destroy (); return false;
1261  }
1262  }
1263  destroy ();
1264  return true;
1265 }
1266 
1267 template <typename T>
1268 bool TSVector<T>::operator == (const TSVector<T>& ts) const
1269 {
1270  if (LIKELY(dim != ts.dim)) {
1271  ts.destroy (); destroy (); return false;
1272  }
1273  if (LIKELY(dim == 0)) {
1274  /*destroy(); ts.destroy();*/ return true;
1275  }
1276  if (LIKELY(fac == ts.fac)) {
1277  if (LIKELY(vec == ts.vec)) {
1278  destroy (); return true;
1279  }
1280  if (TBCICOMP(vec, ts.vec, T, dim)) {
1281  ts.destroy (); destroy (); return false;
1282  }
1283  } else {
1284  if (LIKELY(vec == ts.vec)) {
1285  destroy (); return false;
1286  }
1287  for (REGISTER unsigned long i = 0; i < dim; ++i)
1288  if (UNLIKELY(vec[i]*fac != ts.vec[i]*ts.fac)) {
1289  ts.destroy (); destroy (); return false;
1290  }
1291  }
1292  ts.destroy (); destroy (); return true;
1293 }
1294 
1297 template <typename T>
1298 inline TVector<T> TSVector<T>::operator + (TVector<T> tv) const
1299 {
1300  BCHK(dim!=tv.dim,VecErr,Wrong dims in TSV + TV,tv.dim,tv);
1301  //for (REGISTER unsigned long t = 0; t < dim; t++)
1302  // tv.set (fac*vec[t] + tv.vec[t], t);
1303  STD_SMP_TEMPLATE3VC(vec_add_svc, dim, tv.vec, vec, fac);
1304  destroy (); return tv;
1305 }
1308 template <typename T>
1309 inline TVector<T> TSVector<T>::operator - (TVector<T> tv) const
1310 {
1311  BCHK(dim!=tv.dim,VecErr,Wrong dims in TSV - TV,tv.dim,tv);
1312  //for (REGISTER unsigned long t = 0; t < dim; t++)
1313  // tv.set (fac*vec[t] - tv.vec[t], t);
1314  STD_SMP_TEMPLATE3VC(vec_sub_svc_inv, dim, tv.vec, vec, fac);
1315  destroy (); return tv;
1316 }
1317 
1319 template <typename T>
1320 inline TVector<T> TSVector<T>::operator + (const Vector<T>& v) const
1321 {
1322  TVector<T> tv;
1323  if (LIKELY(mut)) {
1324  tv.vec = vec; tv.dim = dim;
1325  } else
1326  tv.resize (dim);
1327  BCHK(dim!=v.dim,VecErr,Wrong dims in TSV + V,v.dim,tv);
1328  //for (REGISTER unsigned long i=0; i<dim; ++i)
1329  // { tv.vec[i] = (vec[i]*fac + v.vec[i]); }
1330  STD_SMP_TEMPLATE4V(svc_vec_add, dim, tv.vec, vec, v.vec, fac);
1331  return (tv);
1332 }
1334 template <typename T>
1335 inline TVector<T> TSVector<T>::operator - (const Vector<T>& v) const
1336 {
1337  TVector<T> tv;
1338  if (LIKELY(mut)) {
1339  tv.vec = vec; tv.dim = dim;
1340  } else
1341  tv.resize (dim);
1342  BCHK(dim!=v.dim,VecErr,Wrong dims in TSV - V,v.dim,tv);
1343  //for (REGISTER unsigned long i=0; i<dim; ++i)
1344  // { tv.vec[i] = (vec[i]*fac - v.vec[i]); }
1345  STD_SMP_TEMPLATE4V(svc_vec_sub, dim, tv.vec, vec, v.vec, fac);
1346  return (tv);
1347 }
1348 
1350 template <typename T>
1351 inline TVector<T> TSVector<T>::operator + (const TSVector<T>& tsv) const
1352 {
1353  TVector<T> tv;
1354  if (LIKELY(mut)) {
1355  tv.vec = vec; tv.dim = dim;
1356  } else if (LIKELY(tsv.mut)) {
1357  tv.vec = tsv.vec; tv.dim = tsv.dim;
1358  } else
1359  tv.resize (dim);
1360  BCHK(dim!=tsv.dim,VecErr,Wrong dims in TSV + TSV,tsv.dim,tv);
1361  //for (REGISTER unsigned long i=0; i<dim; ++i)
1362  // { tv.vec[i] = vec[i]*fac + tsv.vec[i]*tsv.fac; }
1363  STD_SMP_TEMPLATE5(svc_svc_add, dim, tv.vec, vec, tsv.vec, fac, tsv.fac);
1364  if (UNLIKELY(mut && tsv.mut))
1365  tsv.destroy ();
1366  return tv;
1367 }
1369 template <typename T>
1370 inline TVector<T> TSVector<T>::operator - (const TSVector<T>& tsv) const
1371 {
1372  TVector<T> tv;
1373  if (LIKELY(mut)) {
1374  tv.vec = vec; tv.dim = dim;
1375  } else if (LIKELY(tsv.mut)) {
1376  tv.vec = tsv.vec; tv.dim = tsv.dim;
1377  } else
1378  tv.resize (dim);
1379  BCHK(dim!=tsv.dim,VecErr,Wrong dims in TSV - TSV,tsv.dim,tv);
1380  //for (REGISTER unsigned long i=0; i<dim; ++i)
1381  // { tv.vec[i] = vec[i]*fac - tsv.vec[i]*tsv.fac; }
1382  STD_SMP_TEMPLATE5(svc_svc_sub, dim, tv.vec, vec, tsv.vec, fac, tsv.fac);
1383  if (UNLIKELY(mut && tsv.mut))
1384  tsv.destroy ();
1385  return tv;
1386 }
1387 
1389 template <typename T>
1390 inline TVector<T> TSVector<T>::operator + (const T& v) const
1391 {
1392  TVector<T> tv;
1393  if (LIKELY(mut)) {
1394  tv.vec = vec; tv.dim = dim;
1395  } else
1396  tv.resize (dim);
1397  //for (REGISTER unsigned long i=0; i<dim; ++i)
1398  // { tv.vec[i] = vec[i]*fac op v; }
1399  STD_SMP_TEMPLATE4C(svc_val_add, dim, tv.vec, vec, fac, v);
1400  return tv;
1401 }
1403 template <typename T>
1404 inline TVector<T> TSVector<T>::operator - (const T& v) const
1405 {
1406  TVector<T> tv;
1407  if (LIKELY(mut)) {
1408  tv.vec = vec; tv.dim = dim;
1409  } else
1410  tv.resize (dim);
1411  //for (REGISTER unsigned long i=0; i<dim; ++i)
1412  // { tv.vec[i] = vec[i]*fac op v; }
1413  STD_SMP_TEMPLATE4C(svc_val_sub, dim, tv.vec, vec, fac, v);
1414  return tv;
1415 }
1416 
1419 template <typename T>
1420 inline TVector<T> TSVector<T>::add_t_tsv (const T& v) const
1421 {
1422  TVector<T> tv;
1423  if (LIKELY(mut)) {
1424  tv.vec = vec; tv.dim = dim;
1425  }
1426  else
1427  tv.resize (size());
1428  //for (REGISTER unsigned long i=0; i<dim; ++i)
1429  // { tv.vec[i] = (v op vec[i]*fac); }
1430  STD_SMP_TEMPLATE4C(val_svc_add, dim, tv.vec, vec, v, fac);
1431  return tv;
1432 }
1433 INST(template <typename T> class TSVector<T> friend TVector<T> operator + (const T&, const TSVector<T>&);)
1434 template <typename T>
1435 inline TVector<T> operator + (const T& v, const TSVector<T>& tsv)
1436 { return tsv.add_t_tsv (v); }
1437 
1439 template <typename T>
1440 inline TVector<T> TSVector<T>::sub_t_tsv (const T& v) const
1441 {
1442  TVector<T> tv;
1443  if (LIKELY(mut)) {
1444  tv.vec = vec; tv.dim = dim;
1445  }
1446  else
1447  tv.resize (size());
1448  //for (REGISTER unsigned long i=0; i<dim; ++i)
1449  // { tv.vec[i] = (v op vec[i]*fac); }
1450  STD_SMP_TEMPLATE4C(val_svc_sub, dim, tv.vec, vec, v, fac);
1451  return tv;
1452 }
1453 INST(template <typename T> class TSVector<T> friend TVector<T> operator - (const T&, const TSVector<T>&);)
1454 template <typename T>
1455 inline TVector<T> operator - (const T& v, const TSVector<T>& tsv)
1456 { return tsv.sub_t_tsv (v); }
1457 
1458 template <typename T>
1459 STD__ ostream& operator << (STD__ ostream& os, const TSVector<T>& ts)
1460 {
1461  for (unsigned long i = 0; i < ts.dim; ++i)
1462  os << ts.vec[i]*ts.fac << " ";
1463  ts.destroy (); return os;
1464 }
1465 
1466 
1467 template <typename T>
1468 inline TVector<T> TSVector<T>::incr (const unsigned long wh, const T i) const
1469 {
1470  TVector<T> tv (*this);
1471  BCHK (wh >= dim, VecErr, Tried to change not-exitsting idx, wh, tv);
1472  tv.vec[wh] += i; /* destroy(); */
1473  return tv;
1474 }
1475 
1476 
1477 INST(template <typename T> class TSVector friend double fabssqr (const TSVector<T>&);)
1478 template <typename T>
1479 double fabssqr (const TSVector<T>& ts)
1480 { return ts.fabssqr (); }
1481 
1482 NAMESPACE_END
1483 
1484 NAMESPACE_CSTD
1485 
1486 INST(template <typename T> class TBCI__ TSVector friend double fabs (const TBCI__ TSVector<T>&);)
1487 template <typename T>
1488 double inline fabs (const TBCI__ TSVector<T>& ts)
1489 { return ts.fabs(); }
1490 
1491 INST(template <typename T> class TBCI__ TSVector friend T abs (const TBCI__ TSVector<T>&);)
1492 template <typename T>
1493 inline T abs (const TBCI__ TSVector<T>& ts)
1494 { return ts.abs(); }
1495 
1496 NAMESPACE_CSTD_END
1497 
1498 NAMESPACE_TBCI
1499 
1500 //------------------------------------------------------------
1501 
1525 template <typename T>
1526 class Vector : public TVector<T> //, public Vector_Sig<T>
1527 {
1528  public:
1529  //typedef T TALIGN(sizeof(T)) vec_t;
1530  typedef T value_type;
1531  typedef T element_type;
1532  typedef T aligned_value_type TALIGN(MIN_ALIGN2);
1533 
1534  friend class TSVector<T>;
1535  // constructors are NEVER inherited,
1536  // so we inherit them explicitly
1537  // (or do the default constructors do exactly the same ?)
1538 
1539  explicit Vector (const unsigned long d = 0) : TVector<T> (d) {}
1540  Vector (const T& val, const unsigned long d) : TVector<T> (val, d) {}
1541  Vector (const BVector<T>& bv) : TVector<T> (bv) {}
1542  Vector (const Vector<T>& v) : TVector<T> (v) {}
1543  //alias
1544  Vector (const TSVector<T>& ts)
1545  { ts.eval (this->vec); this->vec = ts.vec; this->dim = ts.dim; }
1546  Vector (const TVector<T>& tv) : TVector<T> (tv) {} //remove const ?
1547  Vector (const Mat_Brack<T>&);
1548 #ifndef NO_POD
1549  Vector (const vararg va, ... );
1550 #endif
1551  //disallow automatic conversion to detect errors in the code
1552  // KG, 99/05/19:
1553  // Danger of recursion: int -> cplx<int> -> cplx<cplx<int> > ... compiler death
1554  // Strange enough, egcs can handle it, if no namespaces are used (!)
1555  #if defined(__GNUC__) && defined VEC_CPLX_CONSTR
1556  explicit Vector (const Vector<cplx<T> >& cv);
1557  #endif
1558 
1559  ~Vector() { this->destroy (); }
1560 
1561  /* Read-only access */
1562  typename tbci_traits<T>::const_refval_type
1563  operator () (const unsigned long i) const
1564  { return BVector<T>::operator () (i); }
1565  typename tbci_traits<T>::const_refval_type
1566  operator [] (const unsigned long i) const
1567  { return this->operator() (i); }
1568  // Read-write access
1569  T& operator () (const unsigned long i)
1570  { return BVector<T>::operator () (i); }
1571  T& operator [] (const unsigned long i)
1572  { return this->operator() (i); }
1573 
1574  //T& operator () (const unsigned long i) const
1575  //{ return BVector<T>::operator () (i); }
1576 
1577  Vector<T>& fill (const T& v) { BVector<T>::fill (v); return *this; }
1578  Vector<T>& clear () { return this->fill((T)0); }
1579 
1580  // Ahhh! Slicing!
1581  TVector<T> slice (unsigned long, unsigned long) const;
1582 
1583  TVector<T> incr (const unsigned long, const T = (T)1) const;
1584 
1585  // assignment operators
1586 
1587  Vector<T>& operator = (const T& a)
1588  { this->TVector<T>::operator = (a); return *this; }
1590  Vector<T>& operator = (const Vector<T>& v) HOT
1591  { this->TVector<T>::operator = (v); return *this; }
1592  Vector<T>& operator = (const TSVector<T>& ts) HOT
1593  { TVector<T>::operator = (ts); return *this; }
1594  Vector<T>& operator = (const TVector<T>& tv) HOT
1595  { this->TVector<T>::operator = (tv); return *this; }
1596 #ifdef TBCI_NEW_BRACKET
1597  Vector<T>& operator = (const T* ptr)
1599  { TBCICOPY(this->vec,ptr,T,this->dim); return *this; }
1600 #endif
1601 
1602  // The ones from TVector
1603  bool operator == (const Vector<T>& v) const
1604  { return BVector<T>::operator == ((BVector<T>)v); }
1605  bool operator != (const Vector<T>& v) const { return !(*this == v); }
1606  bool operator == (const TSVector<T>&) const;
1607  bool operator != (const TSVector<T>& ts) const { return !(*this == ts); }
1608  bool operator == (const TVector<T>& tv) const { return (tv == *this); }
1609  bool operator != (const TVector<T>& tv) const { return !(*this == tv); }
1610 
1611  bool operator <= (const Vector<T>& v) const { return BVector<T>::operator <= (v); }
1612  bool operator >= (const Vector<T>& v) const { return BVector<T>::operator >= (v); }
1613  bool operator < (const Vector<T>& v) const { return !((*this) >= v); }
1614  bool operator > (const Vector<T>& v) const { return !((*this) <= v); }
1615 
1616  TVector<T> operator + (const Vector<T>&) const HOT;
1617  TVector<T> operator - (const Vector<T>&) const HOT;
1618  TVector<T> operator + (/*const*/ TVector<T>/*&*/ tv) const HOT;
1619  TVector<T> operator - (/*const*/ TVector<T>/*&*/ tv) const HOT;
1620  TVector<T> operator + (const TSVector<T>&) const HOT;
1621  TVector<T> operator - (const TSVector<T>&) const HOT;
1622 
1623  TVector<T> operator + (const T&) const;
1624  TVector<T> operator - (const T&) const;
1625  TSVector<T> operator * (const T&) const;
1626  TSVector<T> operator / (const T&) const;
1627 
1628  TVector<T> operator - () const;
1629 
1630  // Additional
1631 
1632  friend T dot FGD (const Vector<T>&, const Vector<T>&) HOT;
1633 
1634  T operator * (const Vector<T>&) const;
1635  T operator * (TVector<T>& tv) const
1636  { Vector<T> v (tv); return this->operator * (v); }
1637  T operator / (const Vector<T>&) const; //NB: a/b := a*b/b^2
1638 
1639  T min () const;
1640  T max () const;
1641  T min (unsigned long& pos) const;
1642  T max (unsigned long& pos) const;
1643  T sum () const;
1644  double fabssqr () const HOT;
1645  double fabs () const { return MATH__ sqrt (this->fabssqr()); }
1646  T abs () const { return (T) MATH__ sqrt (this->fabssqr()); }
1647  /* { T res(0.0); do_vec_sumsqr<T> (dim, vec, res); return res; } */
1648 
1649  bool contains (const T& v) const { return BVector<T>::contains (v); }
1650 
1651  //friend TVector<T> operator + FGD (const T&, const Vector<T>&);
1652  //friend TVector<T> operator - FGD (const T&, const Vector<T>&);
1653  //friend TSVector<T> operator * FGD (const T&, const Vector<T>&);
1654 
1655  static const char* vec_info() { return "Vector"; }
1656  // useful?
1657  //Vector<T> apply(T op(T));
1658 
1659 #if 0
1660  friend STD__ ostream& operator << (STD__ ostream& os, const Vector<T>& v);
1661 #endif
1662  /* An inheritance problem in GCC */
1663 #ifdef OLD_GCC
1664  friend STD__ istream& operator >> (STD__ istream& is, Vector<T>& v)
1665  { return is >> (BVector<T>&)v; }
1666 #endif
1667 
1668 #ifndef HAVE_PROMOTION_BUG
1669  // Promotion (only explicit)
1670 //# ifndef HAVE_GCC295_TMPLFRNDCLS_BUG
1671 // template <typename U> friend class BVector;
1672 //# endif
1673  template <typename U> explicit Vector (const BVector<U>& bv) : TVector<T> (bv) {}
1674 #endif
1675 
1676 } ;
1677 
1678 #if 1
1679 INST(template <typename T> friend inline STD__ ostream& operator << (STD__ ostream& os, const Vector<T>& v);)
1680 template <typename T>
1681 inline STD__ ostream& operator << (STD__ ostream& os, const Vector<T>& v)
1682 { return os << (BVector<T>&)v; }
1683 #endif
1684 
1685 #ifndef NO_POD
1686 template <typename T>
1687 inline Vector<T>::Vector (const vararg va, ...)
1688 {
1689  this->vec = (T*)0; this->dim = (unsigned)va; this->keep = false;
1690  if (UNLIKELY(this->dim)) {
1691  this->vec = NEW (T, this->dim);
1692  if (UNLIKELY(!this->vec))
1693  this->dim = 0;
1694  }
1695  va_list vl;
1696  va_start (vl, va);
1697 #if !defined(__clang__) || !defined(CPLX)
1698  for (unsigned long i=0; i < this->dim; ++i)
1699  this->vec[i] = va_arg (vl, T);
1700 #else
1701  throw VecErr("vararg not supported for cplx in clang");
1702 #warning no vararg support with cplx numbers and clang
1703 #endif
1704  va_end (vl);
1705 }
1706 #endif
1707 
1708 #if defined(__GNUC__) && defined VEC_CPLX_CONSTR
1709 template <typename T>
1710 inline Vector<T>::Vector (const Vector<cplx<T> >& vc) : TVector<T> (vc.size())
1711 {
1712  for (unsigned long i = 0; i < this->dim; ++i)
1713  this->vec[i] = CPLX__ real (vc(i));
1714 }
1715 #endif
1716 
1717 template <typename T>
1718 bool Vector<T>::operator == (const TSVector<T>& tsv) const
1719 {
1720  if (LIKELY(this->dim != tsv.dim)) {
1721  tsv.destroy (); return false;
1722  }
1723  if (LIKELY(this->dim == 0)) {
1724  /*tsv.destroy();*/ return true;
1725  }
1726  if (LIKELY(tsv.fac == (T)1)) {
1727  if (LIKELY(this->vec == tsv.vec))
1728  return true;
1729  if (TBCICOMP(this->vec, tsv.vec, T, this->dim)) {
1730  tsv.destroy (); return false;
1731  }
1732  } else {
1733  if (LIKELY(this->vec == tsv.vec))
1734  return false;
1735  for (REGISTER unsigned long i = 0; i < this->dim; ++i)
1736  if (UNLIKELY(tsv.vec[i]*tsv.fac != this->vec[i])) {
1737  tsv.destroy (); return false;
1738  }
1739  }
1740  tsv.destroy (); return true;
1741 }
1742 
1743 
1744 #if 0
1745 template <typename T>
1746 inline Vector<T>& Vector<T>::operator *= (const T& f)
1747 {
1748  for (unsigned long i = 0; i < this->dim; ++i)
1749  this->vec[i] *= f;
1750  return *this;
1751 }
1752 
1753 template <typename T>
1754 inline Vector<T>& Vector<T>::operator /= (const T& f)
1755 {
1756  for (unsigned long i = 0; i < this->dim; ++i)
1757  this->vec[i] /= f;
1758  return *this;
1759 }
1760 #endif
1761 
1762 #if !defined(SMP) || defined(NOSMP_VECSCALAR)
1763 
1764 // Scalar product with conj
1765 template <typename T>
1766 inline T dot (const Vector<T>& a, const Vector<T>& b)
1767 {
1768  BCHK(a.dim != b.dim, VecErr, Scalar dot product with diff. size Vectors, b.size(), T(0));
1769  T s = T(0.0);
1770  if (do_exactsum())
1771  do_vec_dot_exact <T> (a.dim, a.vec, b.vec, s);
1772  else
1773  do_vec_dot_quick <T> (a.dim, a.vec, b.vec, s);
1774  //STD__ cout << "dot: " << s << STD__ endl;
1775 # ifdef SMPVEC_DEBUG
1776  STD__ cerr << "dot: " << STD__ setprecision (16) << s << STD__ endl;
1777 # endif
1778  return s;
1779 }
1780 
1781 // Scalar product withOUT conj
1782 template <typename T>
1783 inline T Vector<T>::operator * (const Vector<T>& a) const
1784 {
1785  BCHK(this->dim != a.dim, VecErr, Scalar product with diff. size Vectors, a.size(), 0);
1786  T s = T(0.0);
1787  if (do_exactsum())
1788  do_vec_mult_exact <T> (this->dim, this->vec, a.vec, s);
1789  else
1790  do_vec_mult_quick <T> (this->dim, this->vec, a.vec, s);
1791 # ifdef SMPVEC_DEBUG
1792  STD__ cerr << "op *: " << STD__ setprecision (16) << s << STD__ endl;
1793 # endif
1794  return s;
1795 }
1796 
1797 
1798 #else /* SMP */
1799 
1800 INST(template <typename T> class TVector friend void job_vec_dot (struct thr_ctrl*);)
1801 template <typename T>
1802 void job_vec_dot (struct thr_ctrl *tc)
1803 {
1804  if (sizeof(T) <= 16) /* tc->t_par[2] == 0) */
1805  tc->t_par[2] = &tc->t_res_ptr;
1806  *(T*)(tc->t_par[2]) = T(0.0);
1807  if (do_exactsum())
1808  do_vec_dot_exact (tc->t_size, (const T*)(tc->t_par[0]),
1809  (const T*)(tc->t_par[1]), *(T*)tc->t_par[2]);
1810  else
1811  do_vec_dot_quick (tc->t_size, (const T*)(tc->t_par[0]),
1812  (const T*)(tc->t_par[1]), *(T*)tc->t_par[2]);
1813 }
1814 
1815 // Scalar product with conj
1816 template <typename T>
1817 inline T dot (const Vector<T>& a, const Vector<T>& b)
1818 {
1819  BCHK(a.size() != b.size(), VecErr, Scalar dot product with diff. size Vectors, b.size(), 0);
1820  T res = T(0); unsigned long sz = a.size();
1821  /* use some heuristic to decide for the num of threads */
1822  const unsigned n_thr = threads_avail (sz / SMP_VECSLICE2);
1823  update_n_thr(n_thr);
1824  if (LIKELY(n_thr < 2)) {
1825  if (do_exactsum())
1826  do_vec_dot_exact<T> (sz, a.vec, b.vec, res);
1827  else
1828  do_vec_dot_quick<T> (sz, a.vec, b.vec, res);
1829  //STD__ cout << "dot: " << res << STD__ endl;
1830  }
1831  else
1832  {
1833  PREFETCH_R(a.vec, 2); PREFETCH_R(b.vec, 2);
1834  const unsigned long first = slice_offset(1, n_thr, sz, a.vec);
1835  unsigned long st, en = first;
1836  T comp = 0.0;
1837  T* results = (T*)0;
1838  if (sizeof(T) > 16)
1839  results = new T[n_thr];
1840  //_tbci_fill<T>(n_thr, results, T(0.0));
1841  /* Start threads */
1842  for (unsigned t = 0; t < n_thr-1; ++t) {
1843  st = en; en = slice_offset(t+2, n_thr, sz, a.vec);
1844  thread_start (t, (thr_job_t)job_vec_dot<T>, en - st,
1845  a.vec+st, b.vec+st, results+t, (void*)0);
1846  }
1847  /* The first slice is handled by the main thread */
1848  if (do_exactsum())
1849  do_vec_dot_exact<T> (first, a.vec, b.vec, res);
1850  else
1851  do_vec_dot_quick<T> (first, a.vec, b.vec, res);
1852  //sched_yield ();
1853  /* Wait for the end */
1854  //STD__ cout << "dot: " << res;
1855  for (unsigned t = 0; t < n_thr-1; ++t) {
1856  job_output out;
1857  thread_wait(t, &out);
1858  T y;
1859  if (sizeof(T) > 16)
1860  y = results[t];
1861  else {
1862  const T* thrres = (const T*)(&out.t_res_ptr);
1863  y = *thrres;
1864  }
1865  T tmp = res + y;
1866  if (do_exactsum())
1867  comp += (tmp - res) - y;
1868  res = tmp;
1869 #ifdef SMPVEC_DEBUG
1870  STD__ cerr << "+"
1871  << (sizeof(T) > THREAD_MAX_RES_LN?
1872  results[t]: *((T*)(threads[t].t_res_dummy)));
1873 #endif
1874  }
1875  //STD__ cout << ":" << res << STD__ endl;
1876  if (results)
1877  delete[] results;
1878  if (do_exactsum())
1879  res -= comp;
1880  }
1881 #ifdef SMPVEC_DEBUG
1882  STD__ cerr << "dot: " << STD__ setprecision (16) << res << STD__ endl;
1883 #endif
1884  return res;
1885 }
1886 
1887 INST(template <typename T> class TVector friend void job_vec_mult (struct thr_ctrl*);)
1888 template <typename T>
1889 void job_vec_mult (struct thr_ctrl *tc)
1890 {
1891  *(T*)(tc->t_par[2]) = T(0.0);
1892  if (do_exactsum())
1893  do_vec_mult_exact (tc->t_size, (const T*)(tc->t_par[0]),
1894  (const T*)(tc->t_par[1]), *(T*)tc->t_par[2]);
1895  else
1896  do_vec_mult_quick (tc->t_size, (const T*)(tc->t_par[0]),
1897  (const T*)(tc->t_par[1]), *(T*)tc->t_par[2]);
1898 }
1899 
1900 // Scalar product withOUT conj
1901 template <typename T>
1902 inline T Vector<T>::operator * (const Vector<T>& a) const
1903 {
1904  BCHK(this->dim != a.size(), VecErr, Scalar product with diff. size Vectors, a.size(), 0);
1905  T res = T(0); unsigned long sz = this->dim;
1906  /* use some heuristic to decide for the num of threads */
1907  const unsigned n_thr = threads_avail (sz / SMP_VECSLICE2);
1908  update_n_thr(n_thr);
1909  if (LIKELY(n_thr < 2)) {
1910  if (do_exactsum())
1911  do_vec_mult_exact<T> (sz, this->vec, a.vec, res);
1912  else
1913  do_vec_mult_quick<T> (sz, this->vec, a.vec, res);
1914  //STD__ cout << "mult: " << res << STD__ endl;
1915  }
1916  else
1917  {
1918  PREFETCH_R(this->vec, 2); PREFETCH_R(a.vec, 2);
1919  const unsigned long first = sz/n_thr - (sz/n_thr)%4;
1920  unsigned long st, en = first;
1921  T comp(0.0);
1922  T* results = new T[n_thr];
1923  //_tbci_fill<T>(n_thr, results, T(0.0));
1924  /* Start threads */
1925  for (unsigned t = 0; t < n_thr-1; ++t) {
1926  st = en; en = (t+2)*sz / n_thr;
1927  if (en < n_thr-2)
1928  en -= en%4;
1929  thread_start (t, (thr_job_t)job_vec_mult<T>, en - st,
1930  this->vec+st, a.vec+st, results+t, (void*)0);
1931  }
1932  /* The last slice is handled by the main thread */
1933  if (do_exactsum())
1934  do_vec_dot_exact<T> (first, this->vec, a.vec, res);
1935  else
1936  do_vec_dot_quick<T> (first, this->vec, a.vec, res);
1937  //sched_yield ();
1938  /* Wait for the end */
1939  //STD__ cout << "dot: " << res;
1940  for (unsigned t = 0; t < n_thr-1; ++t) {
1941  thread_wait (t);
1942  T tmp = res + results[t];
1943  if (do_exactsum())
1944  comp += (tmp - res) - results[t];
1945  res = tmp;
1946 #ifdef SMPVEC_DEBUG
1947  STD__ cerr << "+"
1948  << (sizeof(T) > THREAD_MAX_RES_LN?
1949  results[t]: *((T*)(threads[t].t_res_dummy)));
1950 #endif
1951  }
1952  //STD__ cout << ":" << res << STD__ endl;
1953  delete[] results;
1954  if (do_exactsum())
1955  res -= comp;
1956  }
1957 #ifdef SMPVEC_DEBUG
1958  STD__ cerr << "mlt: " << STD__ setprecision (16) << res << STD__ endl;
1959 #endif
1960  return res;
1961 }
1962 
1963 #endif /* SMP */
1964 
1965 
1966 
1967 HOTDECL(template <typename T>
1968 bool par_comp (const Vector<T>& v1, const Vector<T>& v2))
1969 {
1970  const unsigned long dim = v1.size();
1971  if (LIKELY(dim != v2.size()))
1972  return false;
1973  if (LIKELY(!dim))
1974  return true;
1975  if (LIKELY(v1(0) != v2(0)))
1976  return false;
1977  return !(_par_comp<T>(dim, &v1.getcref(0), &v2.getcref(0)));
1978 }
1979 
1980 HOTDECL(template <typename T>
1981 void par_copy (Vector<T>& v1, const Vector<T>& v2))
1982 {
1983  const unsigned long dim = v1.size();
1984  BCHK(dim != v2.size(), VecErr, Copy with wrong size, dim, v2.size());
1985  if (LIKELY(!dim))
1986  return;
1987  _par_copy<T>(dim, &v1.setval(0), &v2.getcref(0));
1988 }
1989 
1990 INST(template <typename T> class TVector friend bool par_comp (const Vector<T>&, const Vector<T>&);)
1991 INST(template <typename T> class TVector friend bool par_comp (const Vector<T>&, TVector<T>);)
1992 INST(template <typename T> class TVector friend bool par_comp (TVector<T>, const Vector<T>&);)
1993 INST(template <typename T> class TVector friend bool par_comp (TVector<T>, TVector<T>);)
1994 
1995 template <typename T>
1996 inline bool par_comp (const Vector<T>& v1, TVector<T> v2)
1997 {
1998  return par_comp (v1, (const Vector<T>)v2);
1999  // Note: The cast aliases the TV and the mem is freed ;-)
2000 }
2001 
2002 template <typename T>
2003 inline bool par_comp (TVector<T> v1, const Vector<T>& v2)
2004 {
2005  return par_comp ((const Vector<T>)v1, v2);
2006 }
2007 
2008 template <typename T>
2009 inline bool par_comp (TVector<T> v1, TVector<T> v2)
2010 {
2011  return par_comp ((const Vector<T>)v1, (const Vector<T>)v2);
2012 }
2013 
2014 template <typename T>
2015 inline void par_fill(Vector<T>& vec, vec_fill_fn<T> fn, void* par)
2016 {
2017  const unsigned long dim = vec.size();
2018  if (!dim)
2019  return;
2020  _par_fill_fn(dim, &vec.setval(0), fn, par);
2021 }
2022 
2023 template <typename T>
2024 T Vector<T>::operator / (const Vector<T>& a) const
2025 {
2026  double div ALIGN(MIN_ALIGN) = a.fabssqr();
2027  BCHK(div == 0, VecErr, Division by Vector with norm 0, 0, div);
2028  return (T)((*this * a) / div);
2029 }
2030 
2031 template <typename T>
2032 T Vector<T>::min () const
2033 {
2034  if (LIKELY(this->dim == 0))
2035  return 0;
2036  T min ALIGN(MIN_ALIGN) = this->vec[0];
2037  for (REGISTER unsigned long t = 1; t < this->dim; t++)
2038  if (UNLIKELY(this->vec[t] < min))
2039  min = this->vec[t];
2040  return min;
2041 }
2042 
2043 template <typename T>
2044 T Vector<T>::min (unsigned long& pos) const
2045 {
2046  if (LIKELY(this->dim == 0))
2047  return 0;
2048  BCHK (pos >= this->dim, VecErr, pos outside range, pos, 0);
2049  T min ALIGN(MIN_ALIGN) = this->vec[pos];
2050  for (REGISTER unsigned long t = pos+1; t < this->dim; t++)
2051  if (UNLIKELY(this->vec[t] < min)) {
2052  min = this->vec[t]; pos = t;
2053  }
2054  return min;
2055 }
2056 
2057 template <typename T>
2058 T Vector<T>::max () const
2059 {
2060  if (LIKELY(this->dim == 0))
2061  return 0;
2062  T max ALIGN(MIN_ALIGN) = this->vec[0];
2063  for (REGISTER unsigned long t = 1; t < this->dim; t++)
2064  if (UNLIKELY(this->vec[t] > max))
2065  max = this->vec[t];
2066  return max;
2067 }
2068 
2069 template <typename T>
2070 T Vector<T>::max (unsigned long& pos) const
2071 {
2072  if (LIKELY(this->dim == 0))
2073  return 0;
2074  BCHK (pos >= this->dim, VecErr, pos outside range, pos, 0);
2075  T max ALIGN(MIN_ALIGN) = this->vec[pos];
2076  for (REGISTER unsigned long t = pos+1; t < this->dim; t++)
2077  if (UNLIKELY(this->vec[t] > max)) {
2078  max = this->vec[t]; pos = t;
2079  }
2080  return max;
2081 }
2082 
2083 #ifdef SMP
2084 INST(template <typename T> class Vector friend void job_vec_sum(struct thr_ctrl*);)
2085 template <typename T>
2086 void job_vec_sum(struct thr_ctrl *tc)
2087 {
2088  *(T*)tc->t_par[1] = 0.0;
2089  if (do_exactsum())
2090  do_vec_sum_exact<T>(tc->t_size, (const T*)tc->t_par[0], *(T*)tc->t_par[1]);
2091  else
2092  do_vec_sum_quick<T>(tc->t_size, (const T*)tc->t_par[0], *(T*)tc->t_par[1]);
2093 }
2094 
2095 template <typename T>
2096 T Vector<T>::sum() const
2097 {
2098  T res = 0.0; const unsigned long sz = this->dim;
2099  /* use some heuristic to decide for the num of threads */
2100  const unsigned n_thr = threads_avail (sz / SMP_VECSLICE2);
2101  update_n_thr(n_thr);
2102  if (LIKELY(n_thr < 2)) {
2103  if (do_exactsum())
2104  do_vec_sum_exact<T> (sz, this->vec, res);
2105  else
2106  do_vec_sum_quick<T> (sz, this->vec, res);
2107  }
2108  else {
2109  PREFETCH_R(this->vec, 2);
2110  const unsigned long first = slice_offset(1, n_thr, sz, this->vec);
2111  unsigned long st, en = first;
2112  T comp = 0.0;
2113  T* results = new T[n_thr];
2114  /* Start threads */
2115  for (unsigned t = 0; t < n_thr-1; ++t) {
2116  st = en; en = slice_offset(t+2, n_thr, sz, this->vec);
2117  thread_start (t, (thr_job_t)job_vec_sum<T>, en-st,
2118  this->vec+st, results+t, (void*)0);
2119  }
2120  /* The first slice is handled by the main thread */
2121  if (do_exactsum())
2122  do_vec_sum_exact<T> (first, this->vec, res);
2123  else
2124  do_vec_sum_quick<T> (first, this->vec, res);
2125  //sched_yield ();
2126  /* Wait for the end */
2127  for (unsigned t = 0; t < n_thr-1; ++t) {
2128  thread_wait (t);
2129  T tmp = res + results[t];
2130  if (do_exactsum())
2131  comp += (tmp - res) - results[t];
2132  res = tmp;
2133  }
2134  delete[] results;
2135  if (do_exactsum())
2136  res -= comp;
2137  }
2138  return res;
2139 }
2140 
2141 #else
2142 template <typename T>
2143 T Vector<T>::sum () const
2144 {
2145  T sum ALIGN(MIN_ALIGN) = T(0);
2146  if (do_exactsum())
2147  do_vec_sum_exact<T> (this->dim, this->vec, sum);
2148  else
2149  do_vec_sum_quick<T> (this->dim, this->vec, sum);
2150  return sum;
2151 }
2152 #endif
2153 
2154 #if 0
2155 {
2156  if (!this->dim) return 0.0;
2157  double s = CPLX__ real (vec[0] * CPLX__ conj(vec[0]));
2158  for (unsigned long i = 1; i < this->dim; ++i)
2159  s += CPLX__ real ((this->vec[i] * CPLX__ conj(this->vec[i])));
2160  return (sqrt(s));
2161 // STD__ cout << s << STD__ endl;
2162 // STD__ cout << (*this) * (*this) << STD__ endl;
2163 // return MATH__ fabs( MATH__ sqrt (CPLX__ real(((*this) * (*this)))));
2164 }
2165 #endif
2166 
2168 template <typename T>
2169 inline TVector<T> Vector<T>::operator + (const Vector<T>& v) const
2170 {
2171  TVector<T> tv (this->dim);
2172  BCHK(this->dim!=v.dim, VecErr, Unequal size in V + V, v.dim, tv);
2173  //for (REGISTER unsigned long i=0; i<dim; ++i) tv.vec[i] = vec[i] + v.vec[i];
2174  STD_SMP_TEMPLATE3VV(vec_vec_add, this->dim, tv.vec, this->vec, v.vec);
2175  return tv;
2176 }
2178 template <typename T>
2179 inline TVector<T> Vector<T>::operator - (const Vector<T>& v) const
2180 {
2181  TVector<T> tv (this->dim);
2182  BCHK(this->dim!=v.dim, VecErr, Unequal size in V - V, v.dim, tv);
2183  //for (REGISTER unsigned long i=0; i<dim; ++i) tv.vec[i] = vec[i] - v.vec[i];
2184  STD_SMP_TEMPLATE3VV(vec_vec_sub, this->dim, tv.vec, this->vec, v.vec);
2185  return tv;
2186 }
2187 
2190 template <typename T>
2191 inline TVector<T> Vector<T>::operator + (/*const*/ TVector<T>/*&*/ tv) const
2192 {
2193  BCHK(this->dim!=tv.dim, VecErr, Unequal size in V + TV, tv.dim, tv);
2194  //for (REGISTER unsigned long t=0; t<dim; t++) tv.vec[t] = vec[t] + tv.vec[t];
2195  STD_SMP_TEMPLATE2V(vec_add_vec, this->dim, tv.vec, this->vec);
2196  return tv;
2197 }
2200 template <typename T>
2201 inline TVector<T> Vector<T>::operator - (/*const*/ TVector<T>/*&*/ tv) const
2202 {
2203  BCHK(this->dim!=tv.dim, VecErr, Unequal size in V - TV, tv.dim, tv);
2204  //for (REGISTER unsigned long t=0; t<dim; t++) tv.vec[t] = vec[t] - tv.vec[t];
2205  STD_SMP_TEMPLATE2V(vec_sub_vec_inv, this->dim, tv.vec, this->vec);
2206  return tv;
2207 }
2208 
2209 // No need to destroy ts, as it's either an alias to a Vector (=> non-mutable)
2210 // or an alias to the output TVector
2212 template <typename T>
2213 inline TVector<T> Vector<T>::operator + (const TSVector<T>& ts) const
2214 {
2215  TVector<T> tv;
2216  BCHK(this->dim!=ts.dim, VecErr, Wrong dim in V + TSV, ts.dim, tv);
2217  if (LIKELY(ts.mut)) {
2218  tv.vec = ts.vec; tv.dim = ts.dim;
2219  } else
2220  tv.resize (this->dim);
2221  //for (REGISTER unsigned long t = 0; t < dim; t++)
2222  // tv.vec[t] = vec[t] + ts.fac * ts.vec[t];
2223  STD_SMP_TEMPLATE4V(vec_svc_add, this->dim, tv.vec, this->vec, ts.vec, ts.fac);
2224  return tv;
2225 }
2227 template <typename T>
2228 inline TVector<T> Vector<T>::operator - (const TSVector<T>& ts) const
2229 {
2230  TVector<T> tv;
2231  BCHK(this->dim!=ts.dim, VecErr, Wrong dim in V - TSV, ts.dim, tv);
2232  if (LIKELY(ts.mut)) {
2233  tv.vec = ts.vec; tv.dim = ts.dim;
2234  } else
2235  tv.resize (this->dim);
2236  //for (REGISTER unsigned long t = 0; t < dim; t++)
2237  // tv.vec[t] = vec[t] - ts.fac * ts.vec[t];
2238  STD_SMP_TEMPLATE4V(vec_svc_sub, this->dim, tv.vec, this->vec, ts.vec, ts.fac);
2239  return tv;
2240 }
2241 
2243 template <typename T>
2244 inline TVector<T> Vector<T>::operator + (const T& a) const
2245 {
2246  TVector<T> tv (this->dim);
2247  //for (REGISTER unsigned long i=0; i<dim; ++i) tv.vec[i] = vec[i] + a;
2248  STD_SMP_TEMPLATE3VC(vec_val_add, this->dim, tv.vec, this->vec, a);
2249  return tv;
2250 }
2252 template <typename T>
2253 inline TVector<T> Vector<T>::operator - (const T& a) const
2254 {
2255  TVector<T> tv (this->dim);
2256  //for (REGISTER unsigned long i=0; i<dim; ++i) tv.vec[i] = vec[i] - a;
2257  STD_SMP_TEMPLATE3VC(vec_val_sub, this->dim, tv.vec, this->vec, a);
2258  return tv;
2259 }
2260 
2261 
2262 template <typename T>
2263 inline TSVector<T> Vector<T>::operator * (const T& a) const
2264 {
2265  return TSVector<T> (*this, a);
2266 }
2267 
2268 template <typename T>
2269 inline TSVector<T> Vector<T>::operator / (const T& a) const
2270 {
2271  BCHK (a==(T)0, VecErr, Division by zero in TSVec / T, 0, TSVector<T> (*this, 1));
2272  return TSVector<T> (*this, (T)1 / a);
2273 }
2274 
2275 
2276 INST(template <typename T> class Vector friend TVector<T> operator + (const T&, const Vector<T>&);)
2278 template <typename T>
2279 inline TVector<T> operator + (const T& a, const Vector<T>& v)
2280 {
2281  const unsigned long dim = v.size();
2282  TVector<T> tv (dim);
2283  const T* const vvec(&v.getcref(0));
2284  T* const tvvec(&tv.setval(0));
2285  //for (REGISTER unsigned long i=0; i<v.size(); ++i)
2286  // tv.set (a + v.get(i), i);
2287  STD_SMP_TEMPLATE3VC(val_vec_add, dim, tvvec, vvec, a);
2288  return tv;
2289 }
2291 INST(template <typename T> class Vector friend TVector<T> operator - (const T&, const Vector<T>&);)
2292 template <typename T>
2293 inline TVector<T> operator - (const T& a, const Vector<T>& v)
2294 {
2295  const unsigned long dim = v.size();
2296  TVector<T> tv (dim);
2297  const T* const vvec(&v.getcref(0));
2298  T* const tvvec(&tv.setval(0));
2299  //for (REGISTER unsigned long i=0; i<v.size(); ++i)
2300  // tv.set (a - v.get(i), i);
2301  STD_SMP_TEMPLATE3VC(val_vec_sub, dim, tvvec, vvec, a);
2302  return tv;
2303 }
2304 
2305 
2306 INST(template <typename T> class Vector friend TSVector<T> operator * (const T&, const Vector<T>&);)
2307 template <typename T>
2308 inline TSVector<T> operator * (const T& a, const Vector<T>& b)
2309 {
2310  return TSVector<T> (b, a);
2311 }
2312 
2313 // a/b = a*b/b^2
2314 INSTCTL(#ifndef VEC_INT_DIV_SUPPORT)
2315 INST(template <typename T> class Vector friend TSVector<T> operator / (const T&, const Vector<T>&);)
2316 INSTCTL(#endif)
2317 template <typename T>
2318 TSVector<T> operator / (const T& a, const Vector<T>& b)
2319 {
2320  T dv = a / fabssqr (b);
2321  BCHK(dv==(T)0, VecErr, Divide by norm 0 Vector, 0, 0);
2322  return TSVector<T> (b, dv);
2323 }
2324 
2325 template <typename T>
2326 inline TVector<T> Vector<T>::operator - () const
2327 {
2328  TVector<T> tv (this->dim);
2329  do_vec_neg_vec<T> (this->dim, tv.vec, this->vec);
2330  return tv;
2331 }
2332 
2333 
2334 template <typename T>
2335 inline TVector<T> Vector<T>::slice (unsigned long i0, unsigned long i1) const
2336 {
2337  BCHK (i1 < i0, VecErr, Slicing needs ordered arguments, i0, *this);
2338  BCHK (i0 >= this->dim, VecErr, Slice: arg1 out of range, i0, *this);
2339  BCHK (i1 > this->dim, VecErr, Slice: arg2 out of range, i1, *this);
2340  TVector<T> tv (i1 - i0);
2341  if (UNLIKELY(tv.dim))
2342  TBCICOPY(tv.vec, this->vec+i0, T, tv.dim);
2343  return tv;
2344 }
2345 
2346 
2347 template <typename T>
2348 inline TVector<T> Vector<T>::incr (const unsigned long wh, const T i) const
2349 {
2350  TVector<T> tv (*this);
2351  BCHK (wh >= this->dim, VecErr, Tried to change not-exitsting idx, wh, tv);
2352  tv.vec[wh] += i;
2353  return tv;
2354 }
2355 
2356 
2357 /* Cost for these ops with Vecs of dim d */
2358 #define COST_VECFABSSQR(d) (d*(COST_UNIT_LOAD+COST_MULT+COST_ADD+COST_LOOP))
2359 #define COST_VECSCALAR(d) (d*(2*COST_UNIT_LOAD+COST_MULT+COST_ADD+COST_LOOP))
2360 
2361 //INST(template <typename T> class Vector friend void do_vec_fabssqr (const unsigned long, const T*, double &);)
2362 
2363 INST(template <typename T> class Vector friend double fabssqr (const Vector<T>&);)
2364 template <typename T>
2365 inline double fabssqr (const Vector<T>& v)
2366 { return v.fabssqr(); ;}
2367 
2368 
2369 #if !defined(SMP) || defined(NOSMP_VECFABS)
2370 
2371 template <typename T>
2372 inline double Vector<T>::fabssqr () const
2373 {
2374  double res = 0.0;
2375  if (do_exactsum())
2376  do_vec_fabssqr_exact (this->dim, this->vec, res);
2377  else
2378  do_vec_fabssqr_quick (this->dim, this->vec, res);
2379 #ifdef SMPVEC_DEBUG
2380  STD__ cerr << "fabssqr: " << STD__ setprecision (16) << res << STD__ endl;
2381 #endif
2382  return res;
2383 }
2384 
2385 #else /* SMP */
2386 
2387 INST(template <typename T> class Vector friend void job_vec_fabssqr (struct thr_ctrl*);)
2388 template <typename T>
2389 void job_vec_fabssqr (struct thr_ctrl *tc)
2390 {
2391  tc->t_res_d = 0.0;
2392  if (do_exactsum())
2393  do_vec_fabssqr_exact (tc->t_size, (const T*)(tc->t_par[0]), tc->t_res_d);
2394  else
2395  do_vec_fabssqr_quick (tc->t_size, (const T*)(tc->t_par[0]), tc->t_res_d);
2396 }
2397 
2398 template <typename T>
2399 inline double Vector<T>::fabssqr () const
2400 {
2401  double res = 0.0; unsigned long sz = this->dim;
2402  /* use some heuristic to decide for the num of threads */
2403  const unsigned n_thr = threads_avail (sz / SMP_VECSLICE2);
2404  update_n_thr(n_thr);
2405  if (LIKELY(n_thr < 2)) {
2406  if (do_exactsum())
2407  do_vec_fabssqr_exact<T> (sz, this->vec, res);
2408  else
2409  do_vec_fabssqr_quick<T> (sz, this->vec, res);
2410  }
2411  else {
2412  PREFETCH_R(this->vec, 2);
2413  const unsigned long first = slice_offset(1, n_thr, sz, this->vec);
2414  unsigned long st, en = first;
2415  double comp(0.0);
2416  /* Start threads */
2417  for (unsigned t = 0; t < n_thr-1; ++t) {
2418  st = en; en = slice_offset(t+2, n_thr, sz, this->vec);
2419  thread_start (t, (thr_job_t)job_vec_fabssqr<T>, en-st,
2420  this->vec+st, (void*)0);
2421  }
2422  /* The first slice is handled by the main thread */
2423  if (do_exactsum())
2424  do_vec_fabssqr_exact<T> (first, this->vec, res);
2425  else
2426  do_vec_fabssqr_quick<T> (first, this->vec, res);
2427  //sched_yield ();
2428  /* Wait for the end */
2429  for (unsigned t = 0; t < n_thr-1; ++t) {
2430  double y = thread_wait_result (t);
2431  double tmp = res + y;
2432  if (do_exactsum())
2433  comp += (tmp - res) - y;
2434  res = tmp;
2435  }
2436  if (do_exactsum())
2437  res -= comp;
2438  }
2439 #ifdef SMPVEC_DEBUG
2440  STD__ cerr << "fabssqr: " << STD__ setprecision (16) << res << STD__ endl;
2441 #endif
2442  return res;
2443 }
2444 
2445 #endif /* SMP */
2446 
2447 INST(template <typename T> class TVector friend double fabssqr (TVector<T>);)
2448 template <typename T>
2449 inline double fabssqr (TVector<T> tv)
2450 { return tv.fabssqr (); }
2451 
2452 
2453 /*
2454 template <typename T>
2455 Vector<T>::operator Vector<double> () const
2456 {
2457  Vector<double> v(dim); for (unsigned long i = 0; i < dim; ++i)
2458  v.vec[i] = CPLX__ real(vec[i]);
2459  return v;
2460 }
2461 */
2462 
2463 // Here are the emul, cemul, ediv and cediv variations:
2464 #ifdef INC_VEC_EXTRA
2465 # include "vector_extra.h"
2466 #endif
2467 
2468 NAMESPACE_END
2469 
2470 // ---------------------------------------------------------
2471 // Specializations for int types
2472 //----------------------------------------------------------
2473 
2474 #ifdef VEC_INT_DIV_SUPPORT
2475 
2476 NAMESPACE_TBCI
2477 
2478 //INST(template <> class TSVector friend TSVector<int> operator / (const int&, const TVector<int>&);)
2479 template <>
2480 inline TSVector<int> operator / (const int& a, TVector<int> b)
2481 {
2482  double dv = (double)a / b.fabssqr ();
2483  BCHK(!dv, VecErr, Divide by norm 0 Vector, 0, b);
2484  for (unsigned long i = 0; i < b.size (); ++i)
2485  b.set((int)(b.get(i) * dv), i);
2486  return TSVector<int> (b, (int)1);
2487 }
2488 
2489 template <>
2490 inline TVector<int>& TVector<int>::operator /= (const int& a)
2491 {
2492  BCHK(a==0, VecErr, Divide TVector<int> by zero, 0, *this);
2493  for (REGISTER unsigned long i = 0; i < this->dim; ++i)
2494  this->vec[i] /= a;
2495  return *this;
2496 }
2497 
2498 template <> /* Was commented out. Why ??? */
2499 inline TSVector<int> operator / (const int& a, const Vector<int>& b)
2500 {
2501  double dv = (double)a / b.fabssqr ();
2502  TVector<int> r (b.size ());
2503  BCHK(!dv, VecErr, Divide by norm 0 Vector, 0, r);
2504  for (unsigned long i = 0; i < b.size (); ++i)
2505  r.set ((int)(dv * b.get(i)), i);
2506  return TSVector<int>(r, 1);
2507 }
2508 
2509 template <>
2510 inline TSVector<int> Vector<int>::operator / (const int& a) const
2511 {
2512  TVector<int> r (this->dim);
2513  BCHK(a==0, VecErr, Divide Vector<int> by 0, 0, r);
2514  for (unsigned long i = 0; i < this->dim; ++i)
2515  r.vec[i] = this->vec[i] / a;
2516  return TSVector<int> (r, 1);
2517 }
2518 
2519 // ---------------------------------------------------------
2520 // Specializations for unsigned types
2521 //----------------------------------------------------------
2522 
2523 template <>
2524 inline TVector<unsigned>& TVector<unsigned>::operator /= (const unsigned& a)
2525 {
2526  BCHK(a==0, VecErr, Divide TVector<unsigned> by zero, 0, *this);
2527  for (REGISTER unsigned long i = 0; i < this->dim; ++i)
2528  this->vec[i] /= a;
2529  return *this;
2530 }
2531 
2532 
2533 #ifndef HAVE_WIN_32
2534 //does not compile under VC 6.0 because it has been already used
2535 template <>
2536 inline TSVector<unsigned> operator / (const unsigned& a, TVector<unsigned> b)
2537 {
2538  double dv = (double)a / b.fabssqr ();
2539  BCHK(!dv, VecErr, Divide by norm 0 Vector, 0, b);
2540  for (unsigned long i = 0; i < b.size (); ++i)
2541  b.set((unsigned)(b.get(i) * dv), i);
2542  return TSVector<unsigned> (b, (unsigned)1);
2543 }
2544 #endif /* HAVE_WIN_32 */
2545 
2546 template <>
2547 inline TSVector<unsigned> operator / (const unsigned& a, const Vector<unsigned>& b)
2548 {
2549  double dv = (double)a / b.fabssqr ();
2550  BCHK(!dv, VecErr, Divide by norm 0 Vector, 0, b);
2551  TVector<unsigned> r (b.size ());
2552  for (unsigned long i = 0; i < b.size (); ++i)
2553  r.set ((unsigned)(dv * b.get(i)), i);
2554  return TSVector<unsigned>(r, 1);
2555 }
2556 
2557 template <>
2558 inline TSVector<unsigned> Vector<unsigned>::operator / (const unsigned& a) const
2559 {
2560  TVector<unsigned> r (this->dim);
2561  BCHK(a==0, VecErr, Divide Vector<unsigned> by 0, 0, r);
2562  for (unsigned long i = 0; i < this->dim; ++i)
2563  r.vec[i] = this->vec[i] / a;
2564  return TSVector<unsigned> (r, 1);
2565 }
2566 
2567 NAMESPACE_END
2568 
2569 #endif /* VEC_INT_DIV_SUPPORT */
2570 
2571 #endif /* TBCI_VECTOR_H */
TVector::operator-
TSVector< T > operator-()
Definition: vector.h:851
TBCICOPY
#define TBCICOPY(n, o, t, s)
Definition: basics.h:894
BVector::end
iterator end()
Definition: bvector.h:155
TVector::TVector
TVector(const T &val, const unsigned long d)
Definition: vector.h:85
thr_ctrl::t_par
void * t_par[6]
Definition: smp.h:173
BVector::operator<=
bool operator<=(const BVector< T > &bv) const
Definition: bvector.h:509
sum
T sum(const FS_Vector< dims, T > &fv)
Definition: fs_vector.h:599
BVector::operator()
T & operator()(const unsigned long) HOT
Definition: bvector.h:469
Vector::fill
Vector< T > & fill(const T &v)
Definition: vector.h:1577
TSVector::incr
TVector< T > incr(const unsigned long, const T=(T) 1) const
Definition: vector.h:1468
TVector::setval
T & setval(const T &val, const unsigned long i) const
Definition: vector.h:193
TVector::fabs
double fabs()
Definition: vector.h:151
Vector_Sig
abstract base class (signature) for Vectors with arithmetics
Definition: vector_sig.h:24
false
#define false
Definition: bool.h:23
Vector::value_type
T value_type
Definition: vector.h:1530
job_val_div_vec
void job_val_div_vec(struct thr_ctrl *tc)
vec = val/self;
Definition: vector.h:427
do_mat_tsv_mult
void do_mat_tsv_mult(const unsigned start, const unsigned end, TVector< T > *res, const Matrix< T > *mat, const TSVector< T > *vec)
MIN_ALIGN
#define MIN_ALIGN
Definition: basics.h:421
STD_SMP_TEMPLATE3VV
#define STD_SMP_TEMPLATE3VV(oper, dm, a1, a2, a3)
Definition: vector.h:648
ALIGN
#define ALIGN(x)
Definition: basics.h:444
fac
tm fac
Definition: f_matrix.h:1052
Vector::abs
T abs() const
Definition: vector.h:1646
Vector::contains
bool contains(const T &v) const
Definition: vector.h:1649
do_mat_vec_mult
void do_mat_vec_mult(const unsigned start, const unsigned end, TVector< T > *res, const Matrix< T > *mat, const Vector< T > *vec)
job_val_vec_mul
void job_val_vec_mul(struct thr_ctrl *tc)
vec = val * vec;
Definition: vector.h:347
TSVector::FGD
friend NOINST void FRIEND_TBCI2__ do_mat_tsv_mult FGD(const unsigned start, const unsigned end, TVector< T > *res, const Matrix< T > *mat, const TSVector< T > *rsv)
TVector::setval
T & setval(const unsigned long i) const
Definition: vector.h:192
TSVector::TSVector
TSVector()
Definition: vector.h:1053
TVector::min
T min()
Definition: vector.h:146
job_val_vec_div
void job_val_vec_div(struct thr_ctrl *tc)
vec = val / vec;
Definition: vector.h:371
job_svc_vec_sub
void job_svc_vec_sub(struct thr_ctrl *tc)
vec = s*vec - vec;
Definition: vector.h:473
TVector::sum
T sum()
Definition: vector.h:148
BVector
provides basic Vector functionality but arithmetic operators (+=, - , *, /...).
Definition: bvector.h:67
TVector::operator[]
T operator[](const unsigned long i)
Definition: vector.h:185
job_vec_add_svc
void job_vec_add_svc(struct thr_ctrl *tc)
vec += s*vec;
Definition: vector.h:491
TSVector::fabs
double fabs() const
Definition: vector.h:1126
BVector::destroy
void destroy()
Definition: bvector.h:268
Vector::vec_info
static const char * vec_info()
Definition: vector.h:1655
TSVector::TSVector
TSVector(const TSVector< T > &ts)
Definition: vector.h:1057
TSVector::element_type
T element_type
Definition: vector.h:1048
Vector::operator+
TVector< T > operator+(const Vector< T > &) const HOT
TV = V + V.
Definition: vector.h:2169
TSVector::clone
void clone(const bool=false, const T *=0) const
Definition: vector.h:1200
TVector::TVector
TVector(const BVector< U > &bv)
Definition: vector.h:228
VEC_INT_DIV_SUPPORT
#define VEC_INT_DIV_SUPPORT
Definition: veclib.cc:12
BVector::fill
BVector< T > & fill(const T &) HOT
Definition: bvector.h:378
BdMatrix
The class BdMatrix is an implementation to store and do operations on sparse Matrices with a band str...
Definition: band_matrix.h:103
operator>>
STD__ istream & operator>>(STD__ istream &istr, BdMatrix< T > &mat)
Definition: band_matrix.h:2739
BdMatrix::operator+=
BdMatrix< T > & operator+=(const BdMatrix< T > &)
Definition: band_matrix.h:1314
BVector::operator==
bool operator==(const BVector< T > &) const HOT
KG, 2001-06-29: Strange: If we don&#39;t inline this, we seems to get better performance in our solver be...
Definition: bvector.h:495
end
const unsigned end
Definition: matrix_kernels.h:25
Vector::Vector
Vector(const TSVector< T > &ts)
Definition: vector.h:1544
Vector::max
T max() const
Definition: vector.h:2058
cplx
Our own complex class.
Definition: cplx.h:48
REGISTER
#define REGISTER
Definition: basics.h:108
std::fabs
double fabs(const int a)
Definition: basics.h:1215
Vector::operator!=
bool operator!=(const Vector< T > &v) const
Definition: vector.h:1605
TVector::operator+
TVector< T > & operator+(const Vector< T > &) HOT
Definition: vector.h:737
FRIEND_TBCI__
#define FRIEND_TBCI__
Definition: basics.h:334
NAMESPACE_TBCI
#define NAMESPACE_TBCI
Definition: basics.h:317
Mat_Brack
Definition: matrix.h:75
TSVector::FGDT
friend NOINST double FRIEND_TBCI__ fabssqr FGDT(const TSVector< T > &)
TVector::operator+=
TVector< T > & operator+=(const T &)
TV += a.
Definition: vector.h:666
thread_start
void thread_start(const int thr_no, thr_job_t job, const unsigned long sz,...)
Definition: smp.cc:988
job_val_vec_add
void job_val_vec_add(struct thr_ctrl *tc)
vec = val + vec;
Definition: vector.h:355
BdMatrix::dim
unsigned int dim
Definition: band_matrix.h:109
job_vec_div_val
void job_vec_div_val(struct thr_ctrl *tc)
vec /= val;
Definition: vector.h:419
Vector::element_type
T element_type
Definition: vector.h:1531
job_vec_mul_val
void job_vec_mul_val(struct thr_ctrl *tc)
vec *= val;
Definition: vector.h:411
TSVector::vecptr
T *const & vecptr() const
Definition: vector.h:1163
thr_ctrl::t_res_ptr
void * t_res_ptr
Definition: smp.h:180
TVector::abs
T abs()
Definition: vector.h:152
TVector::get
tbci_traits< T >::const_refval_type get(const unsigned long i) const
Definition: vector.h:190
MIN_ALIGN2
#define MIN_ALIGN2
Definition: basics.h:424
Vector::operator-
TVector< T > operator-() const
Definition: vector.h:2326
Vector::Vector
Vector(const BVector< T > &bv)
Definition: vector.h:1541
TSVector::operator()
T operator()(const unsigned long i) const
Definition: vector.h:1066
BVector::dim
unsigned long dim
Definition: bvector.h:74
TVector::vec_info
static const char * vec_info()
Definition: vector.h:219
TSVector::add_t_tsv
TVector< T > add_t_tsv(const T &) const
Helper member fn to prevent friendship TV = T + TSV.
Definition: vector.h:1420
job_svc_val_add
void job_svc_val_add(struct thr_ctrl *tc)
vec = s*vec + val;
Definition: vector.h:519
Vector::Vector
Vector(const T &val, const unsigned long d)
Definition: vector.h:1540
STD_SMP_TEMPLATE5
#define STD_SMP_TEMPLATE5(oper, dm, a1, a2, a3, a4, a5)
Definition: vector.h:658
BCHK
#define BCHK(cond, exc, txt, ind, rtval)
Definition: basics.h:575
par_fill
void par_fill(Matrix< T > &mat, mat_fill_fn< T > fn, void *par)
Definition: matrix.h:2417
TSVector::operator*
const TSVector< T > & operator*(const T &f) const
Definition: vector.h:1093
TVector::element_type
T element_type
Definition: vector.h:81
TVector::value_type
T value_type
Definition: vector.h:80
job_svc_vec_add
void job_svc_vec_add(struct thr_ctrl *tc)
vec = s*vec + vec;
Definition: vector.h:446
BVector::keep
bool keep
Definition: bvector.h:75
Vector::Vector
Vector(const unsigned long d=0)
Definition: vector.h:1539
NAMESPACE_CSTD_END
#define NAMESPACE_CSTD_END
Definition: basics.h:325
TSVector::getfac
const T & getfac() const
Definition: vector.h:1073
max
T max(const FS_Vector< dims, T > &fv)
Definition: fs_vector.h:594
TSVector::operator-
const TSVector< T > & operator-() const
Definition: vector.h:1096
Vector::fabs
double fabs() const
Definition: vector.h:1645
UNLIKELY
#define UNLIKELY(expr)
Definition: basics.h:101
thr_ctrl
Definition: smp.h:168
TSVector::operator/
const TSVector< T > & operator/(const T &f) const
Definition: vector.h:1094
TVector::max
T max()
Definition: vector.h:147
thr_job_t
void(* thr_job_t)(struct thr_ctrl *)
Before the double inclusion guard on purpose!
Definition: smp.h:126
TSVector::value_type
T value_type
Definition: vector.h:1047
job_vec_vec_add
void job_vec_vec_add(struct thr_ctrl *tc)
vec = vec + vec;
Definition: vector.h:283
Vector::slice
TVector< T > slice(unsigned long, unsigned long) const
Definition: vector.h:2335
FRIEND_TBCI2__
#define FRIEND_TBCI2__
Definition: basics.h:335
job_vec_sub_vec
void job_vec_sub_vec(struct thr_ctrl *tc)
vec -= vec;
Definition: vector.h:307
job_val_svc_add
void job_val_svc_add(struct thr_ctrl *tc)
vec = val + s*vec;
Definition: vector.h:537
TVector::TVector
TVector(const BVector< T > &bv)
Definition: vector.h:86
job_svc_svc_sub
void job_svc_svc_sub(struct thr_ctrl *tc)
vec = s*vec - s*vec;
Definition: vector.h:482
TSVector::getcref
const T & getcref(const unsigned long i) const HOT
Definition: vector.h:1075
BVector::operator=
BVector< T > & operator=(const T &a)
Definition: bvector.h:168
Vector::operator>=
bool operator>=(const Vector< T > &v) const
Definition: vector.h:1612
TVector::set
T & set(const T &val, const unsigned long i) const
Definition: vector.h:194
TVector::operator==
bool operator==(const TVector< T > &tv) const
Definition: vector.h:114
F_TSMatrix
Temporary object for scaled matrices.
Definition: cscmatrix.h:50
STD_SMP_TEMPLATE4C
#define STD_SMP_TEMPLATE4C(oper, dm, a1, a2, a3, a4)
Definition: vector.h:656
Vector::operator>
bool operator>(const Vector< T > &v) const
Definition: vector.h:1614
TVector::operator()
T operator()(const unsigned long i)
Definition: vector.h:182
BVector::operator>=
bool operator>=(const BVector< T > &bv) const
Definition: bvector.h:520
job_vec_svc_add
void job_vec_svc_add(struct thr_ctrl *tc)
vec = vec + s*vec;
Definition: vector.h:437
TSVector::destroy
void destroy() const
Definition: vector.h:1172
Vector::clear
Vector< T > & clear()
Definition: vector.h:1578
job_vec_fabssqr
void job_vec_fabssqr(struct thr_ctrl *tc)
Definition: vector.h:2389
NEW
#define NEW(t, s)
Definition: malloc_cache.h:633
std::sqrt
double sqrt(const int a)
Definition: basics.h:1216
Symm_BdMatrix
Matrix class with optimized Matrix-Vector multiplication for symmetrical Matrices.
Definition: symm_bdmatrix.h:25
job_val_add_vec
void job_val_add_vec(struct thr_ctrl *tc)
vec += val;
Definition: vector.h:387
TSVector::TALIGN
T aligned_value_type TALIGN(MIN_ALIGN2)
Definition: vector.h:1049
BdMatrix::size
unsigned int size() const
size: incompatibility to Matrix (!)
Definition: band_matrix.h:243
TSVector::operator*=
TSVector< T > & operator*=(const T &f)
Definition: vector.h:1091
TSVector::operator=
TSVector< T > & operator=(const TSVector< T > &ts)
Definition: vector.h:1078
threads
struct thr_struct * threads
Definition: smp.cc:106
TVector::operator/
TSVector< T > operator/(const T &)
Definition: vector.h:836
CSTD__
#define CSTD__
Definition: basics.h:340
Vector::TALIGN
T aligned_value_type TALIGN(MIN_ALIGN2)
Definition: vector.h:1532
TVector::operator=
TVector< T > & operator=(const T &a)
Definition: vector.h:106
thr_ctrl::t_size
unsigned long t_size
Definition: smp.h:171
PREFETCH_R
#define PREFETCH_R(addr, loc)
In case gcc does not yet support __builtin_prefetch(), we have handcoded assembly with gcc for a few ...
Definition: basics.h:748
Vector::fabssqr
double fabssqr() const HOT
Definition: vector.h:2399
TSVector::min
T min()
Definition: vector.h:1122
job_vec_sub_vec_inv
void job_vec_sub_vec_inv(struct thr_ctrl *tc)
vec -= vec; vec = -vec;
Definition: vector.h:315
job_vec_add_val
void job_vec_add_val(struct thr_ctrl *tc)
vec += val;
Definition: vector.h:379
TVector::TVector
TVector(const unsigned long d=0)
Definition: vector.h:84
BVector::contains
bool contains(const T &, unsigned long *=0) const
Definition: bvector.h:580
job_vec_svc_sub
void job_vec_svc_sub(struct thr_ctrl *tc)
vec = vec - s*vec;
Definition: vector.h:464
CRMatrix
C++ class for sparse matrices using compressed row storage.
Definition: crmatrix.h:63
ts
F_TSMatrix< T > ts
Definition: f_matrix.h:1052
TVector::incr
TVector< T > & incr(const unsigned long, const T=(T) 1)
Definition: vector.h:973
THREAD_MAX_RES_LN
#define THREAD_MAX_RES_LN
Definition: smp.h:130
TBCI__
#define TBCI__
Definition: basics.h:332
TSVector::operator!=
bool operator!=(const Vector< T > &v) const
Definition: vector.h:1156
Vector::operator=
Vector< T > & operator=(const T &a)
Definition: vector.h:1587
TVector::slice
TVector< T > slice(const unsigned long, const unsigned long)
Definition: vector.h:956
job_svc_val_sub
void job_svc_val_sub(struct thr_ctrl *tc)
vec = s*vec - val;
Definition: vector.h:528
BdMatrix::destroy
void destroy()
destroy object explicitly
Definition: band_matrix.h:156
TSVector::max
T max()
Definition: vector.h:1123
TSVector::sum
T sum()
Definition: vector.h:1131
job_vec_sum
void job_vec_sum(struct thr_ctrl *tc)
Definition: vector.h:2086
job_vec_val_mul
void job_vec_val_mul(struct thr_ctrl *tc)
vec = vec * val;
Definition: vector.h:339
job_vec_dot
void job_vec_dot(struct thr_ctrl *tc)
Definition: vector.h:1802
b
F_TMatrix< T > b
Definition: f_matrix.h:736
TVector::operator*
TSVector< T > operator*(const T &)
Definition: vector.h:829
do_mat_vec_transmult
void do_mat_vec_transmult(const unsigned start, const unsigned end, TVector< T > *res, const Matrix< T > *mat, const Vector< T > *vec)
F_Matrix
Definition: bvector.h:51
Vector::operator/
TSVector< T > operator/(const T &) const
Definition: vector.h:2269
CRMatrix::destroy
void destroy()
Definition: crmatrix.h:669
TSVector::vec
T * vec
Definition: vector.h:1029
Vector::sum
T sum() const
Definition: vector.h:2096
GLBL__
#define GLBL__
Definition: basics.h:342
job_val_sub_vec
void job_val_sub_vec(struct thr_ctrl *tc)
vec -= val; vec = -vec;
Definition: vector.h:403
if
if(value==0) return 1
TVector::FGD
friend TVector< T > conj FGD(const Vector< T > &)
INSTCTL
#define INSTCTL(x)
Definition: basics.h:245
Vector::Vector
Vector(const TVector< T > &tv)
Definition: vector.h:1546
std::conj
int conj(const int arg)
conj for elementary types
Definition: basics.h:1055
NAMESPACE_CSTD
#define NAMESPACE_CSTD
Definition: basics.h:319
par_comp
bool par_comp(const Vector< T > &v1, TVector< T > v2)
Definition: vector.h:1996
TSVector::operator/=
TSVector< T > & operator/=(const T &f)
Definition: vector.h:1092
TVector::TVector
TVector(const TVector< T > &tv) HOT
Definition: vector.h:89
TSVector::detach
void detach(const T *=0) const
Definition: vector.h:1179
Vector::Vector
Vector(const Vector< T > &v)
Definition: vector.h:1542
TVector::getcref
const T & getcref(const unsigned long i) const
Definition: vector.h:188
job_vec_vec_sub
void job_vec_vec_sub(struct thr_ctrl *tc)
vec = vec - vec;
Definition: vector.h:291
TSVector::TSVector
TSVector(const unsigned long d)
Definition: vector.h:1055
TSVector::ALIGN
T fac ALIGN(MIN_ALIGN2)
TVector::size
unsigned long size() const
Definition: vector.h:104
STD_SMP_TEMPLATE2C
#define STD_SMP_TEMPLATE2C(oper, dm, a1, a2)
Definition: vector.h:646
TSVector::sub_t_tsv
TVector< T > sub_t_tsv(const T &) const
TV = T - TSV.
Definition: vector.h:1440
TSVector::mut
bool mut
Definition: vector.h:1037
TVector::TALIGN
T aligned_value_type TALIGN(MIN_ALIGN2)
Definition: vector.h:82
Vector::FGD
friend T dot FGD(const Vector< T > &, const Vector< T > &) HOT
Matrix
Definition: bvector.h:49
dims
#define dims
Definition: fsveclib.cc:45
TVector::operator*=
TVector< T > & operator*=(const T &)
TV *= a.
Definition: vector.h:680
detach
tm detach()
CPLX__
#define CPLX__
Definition: basics.h:341
INST
#define INST(x)
Definition: basics.h:238
job_val_svc_div
void job_val_svc_div(struct thr_ctrl *tc)
vec = val / s*vec;
Definition: vector.h:555
TBCI::dot
T dot(const T &a1, const T &a2)
Definition: basics.h:1183
TSVector::getfacref
T & getfacref() const
Definition: vector.h:1074
TVector::operator!=
bool operator!=(const TVector< T > &tv) const
Definition: vector.h:118
F_TMatrix
Temporary Base Class (non referable!) (acc.
Definition: bvector.h:50
i
int i
Definition: LM_fit.h:71
real
#define real
Definition: own_lapack_routines.cpp:20
Vector::incr
TVector< T > incr(const unsigned long, const T=(T) 1) const
Definition: vector.h:2348
BVector::resize
BVector< T > & resize(const BVector< T > &)
Actually it&#39;s a resize and copy (some people would expect the assignment op to do this) ...
Definition: bvector.h:361
CSCMatrix
exception class: Use MatErr from matrix.h
Definition: cscmatrix.h:49
STD_SMP_TEMPLATE2V
#define STD_SMP_TEMPLATE2V(oper, dm, a1, a2)
Definition: vector.h:644
TSVector::eval
const TSVector< T > & eval(const T *vv=0) const
Definition: vector.h:1221
job_svc_svc_add
void job_svc_svc_add(struct thr_ctrl *tc)
vec = s*vec + s*vec;
Definition: vector.h:455
TSVector::vec_info
static const char * vec_info()
Definition: vector.h:1165
STD__
#define STD__
Definition: basics.h:338
TBCIDELETE
#define TBCIDELETE(t, v, sz)
Definition: malloc_cache.h:634
TSVector::operator[]
T operator[](const unsigned long i) const
Definition: vector.h:1069
job_output
Definition: smp.h:141
threads_avail
#define threads_avail(x)
Definition: smp.h:322
TSVector::operator+
TVector< T > operator+(const Vector< T > &) const HOT
TV = TSV + V.
Definition: vector.h:1320
TVector
Temporary Base Class Idiom: Class TVector is used for temporary variables.
Definition: bvector.h:52
TSVector::size
unsigned long size() const
Definition: vector.h:1120
job_vec_sub_val
void job_vec_sub_val(struct thr_ctrl *tc)
vec -= val;
Definition: vector.h:395
job_val_svc_sub
void job_val_svc_sub(struct thr_ctrl *tc)
vec = val + s*vec;
Definition: vector.h:546
TSVector::dim
unsigned long dim
Definition: vector.h:1029
thread_wait
void thread_wait(const int thr_no, struct job_output *out)
Definition: smp.cc:997
thread_wait_result
double thread_wait_result(const int thr_no)
Definition: smp.cc:1017
abs
#define abs(x)
Definition: f2c.h:178
STD_SMP_TEMPLATE3VC
#define STD_SMP_TEMPLATE3VC(oper, dm, a1, a2, a3)
Definition: vector.h:650
CRMatrix::size
unsigned int size() const
Definition: crmatrix.h:193
TSVector::operator==
bool operator==(const Vector< T > &) const
Definition: vector.h:1241
Vector::min
T min() const
Definition: vector.h:2032
NAMESPACE_END
#define NAMESPACE_END
Definition: basics.h:323
BdMatrix::operator*=
BdMatrix< T > & operator*=(const T &)
Definition: band_matrix.h:1334
job_output::t_res_ptr
void * t_res_ptr
Definition: smp.h:149
Vector::Vector
Vector(const BVector< U > &bv)
Definition: vector.h:1673
job_val_vec_sub
void job_val_vec_sub(struct thr_ctrl *tc)
vec = val - vec;
Definition: vector.h:363
real
float real
Definition: f2c.h:31
Vector
Definition: bvector.h:54
SMP_VECSLICE2
#define SMP_VECSLICE2
Definition: bvector.h:626
HOTDECL
#define HOTDECL(x)
Definition: basics.h:497
y
const Vector< T > Vector< T > Vector< T > Vector< T > & y
Definition: LM_fit.h:172
res
const Vector< T > Vector< T > Vector< T > Vector< T > Vector< T > Vector< T > Vector< T > Vector< T > long int res
Definition: LM_fit.h:199
TVector::TVector
TVector(const Vector< T > &v)
Definition: vector.h:87
BVector::vec
T * vec
Definition: bvector.h:73
TBCI::fabssqr
double fabssqr(const double a)
Definition: basics.h:1157
T
#define T
Definition: bdmatlib.cc:20
TBCICOMP
#define TBCICOMP(n, o, t, s)
Definition: basics.h:979
HOT
#define HOT
Definition: basics.h:495
norm
double norm(const TBCI__ cplx< T > &c)
Definition: cplx.h:695
job_vec_add_vec
void job_vec_add_vec(struct thr_ctrl *tc)
vec += vec;
Definition: vector.h:299
true
#define true
Definition: bool.h:24
job_vec_val_sub
void job_vec_val_sub(struct thr_ctrl *tc)
vec = vec - val;
Definition: vector.h:331
operator/
cplx< T > operator/(const T a, const cplx< T > &b)
Definition: cplx.h:342
Vector::operator[]
tbci_traits< T >::const_refval_type operator[](const unsigned long i) const
Definition: vector.h:1566
job_vec_val_add
void job_vec_val_add(struct thr_ctrl *tc)
vec = vec + val;
Definition: vector.h:323
do_exactsum
unsigned int do_exactsum()
Definition: tbci_param.h:150
a
const unsigned TMatrix< T > const Matrix< T > * a
Definition: matrix_kernels.h:25
job_vec_sub_svc_inv
void job_vec_sub_svc_inv(struct thr_ctrl *tc)
vec -= s*vec;
Definition: vector.h:509
job_vec_sub_svc
void job_vec_sub_svc(struct thr_ctrl *tc)
vec -= s*vec;
Definition: vector.h:500
Vector::~Vector
~Vector()
Definition: vector.h:1559
TSVector::abs
T abs() const
Definition: vector.h:1124
TVector::fabssqr
double fabssqr()
Definition: vector.h:150
TVector::operator/=
TVector< T > & operator/=(const T &)
TV /= a.
Definition: vector.h:688
v
const Vector< T > Vector< T > Vector< T > Vector< T > Vector< T > Vector< T > Vector< T > Vector< T > long int int char v
&lt; find minimun of func on grid with resolution res
Definition: LM_fit.h:205
Vector::operator==
bool operator==(const Vector< T > &v) const
Definition: vector.h:1603
min
#define min(a, b)
Definition: f2c.h:180
TVector::operator-=
TVector< T > & operator-=(const T &)
TV -= a.
Definition: vector.h:673
TSVector::~TSVector
~TSVector()
Definition: vector.h:1054
MATH__
#define MATH__
Definition: basics.h:339
Vector::operator()
tbci_traits< T >::const_refval_type operator()(const unsigned long i) const
Definition: vector.h:1563
TBCI::vararg
enum _vararg vararg
Definition: basics.h:1276
BdMatrix::operator-=
BdMatrix< T > & operator-=(const BdMatrix< T > &)
Definition: band_matrix.h:1315
TSVector::fabssqr
double fabssqr() const
Definition: vector.h:1128
LIKELY
#define LIKELY(expr)
branch prediction note that we sometimes on purpose mark the unlikely possibility likely and vice ver...
Definition: basics.h:100
vector_extra.h
contains rarely used functions on TBCI::Vector s
Vector::operator*
TSVector< T > operator*(const T &) const
Definition: vector.h:2263
TVector::TVector
TVector(const TSVector< T > &ts)
Definition: vector.h:91
thr_ctrl::t_res_d
double t_res_d
Definition: smp.h:178
TVector::~TVector
~TVector()
Definition: vector.h:95
GLBL2__
#define GLBL2__
Definition: basics.h:343
std::imag
int imag(const int d)
Definition: basics.h:1068
VecErr
exception class
Definition: bvector.h:31
job_vec_mult
void job_vec_mult(struct thr_ctrl *tc)
Definition: vector.h:1889
STD_SMP_TEMPLATE4V
#define STD_SMP_TEMPLATE4V(oper, dm, a1, a2, a3, a4)
Definition: vector.h:654
Generated by   doxygen 1.8.5