TBCI Numerical high perf. C++ Library 2.8.0
vector.h
Go to the documentation of this file.
1
12
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
39NAMESPACE_TBCI
40
41//template <typename T> class cplx;
42
43template <typename T> class Vector;
44template <typename T> class TVector;
45template <typename T> class TSVector;
46template <typename T> class Matrix;
47template <typename T> class TMatrix;
48template <typename T> class TSMatrix;
49template <typename T> class F_Matrix;
50template <typename T> class F_TMatrix;
51template <typename T> class F_TSMatrix;
52template <typename T> class Mat_Brack;
53template <typename T> class CRMatrix;
54template <typename T> class CSCMatrix;
55template <typename T> class BdMatrix;
56template <typename T> class Symm_BdMatrix;
57
58#ifdef PRAGMA_I
59# pragma interface "vector.h"
60#endif
61
62//------------------------------------------------------------------------
63
70
71template <typename T>
72class 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; }
110
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
241template <typename T>
242inline 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
252template <typename T>
253inline STD__ ostream& operator << (STD__ ostream& os, const TVector<T>& tv)
254{ Vector<T> v(tv); return os << (BVector<T>)v; }
255
256template <typename T>
257inline 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/*
269template <typename T>
270inline 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
280INST(template <typename T> class Vector friend void job_vec_vec_add (struct thr_ctrl *);)
282template <typename T>
283void 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}
288INST(template <typename T> class Vector friend void job_vec_vec_sub (struct thr_ctrl *);)
290template <typename T>
291void 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}
296INST(template <typename T> class Vector friend void job_vec_add_vec (struct thr_ctrl *);)
298template <typename T>
299void 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}
304INST(template <typename T> class Vector friend void job_vec_sub_vec (struct thr_ctrl *);)
306template <typename T>
307void 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}
312INST(template <typename T> class Vector friend void job_vec_sub_vec_inv (struct thr_ctrl *);)
314template <typename T>
315void 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}
320INST(template <typename T> class Vector friend void job_vec_val_add (struct thr_ctrl *);)
322template <typename T>
323void 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}
328INST(template <typename T> class Vector friend void job_vec_val_sub (struct thr_ctrl *);)
330template <typename T>
331void 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}
336INST(template <typename T> class Vector friend void job_vec_val_mul (struct thr_ctrl *);)
338template <typename T>
339void 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}
344INST(template <typename T> class Vector friend void job_val_vec_mul (struct thr_ctrl *);)
346template <typename T>
347void 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}
352INST(template <typename T> class Vector friend void job_val_vec_add (struct thr_ctrl *);)
354template <typename T>
355void 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}
360INST(template <typename T> class Vector friend void job_val_vec_sub (struct thr_ctrl *);)
362template <typename T>
363void 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}
368INST(template <typename T> class Vector friend void job_val_vec_div (struct thr_ctrl *);)
370template <typename T>
371void 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}
376INST(template <typename T> class Vector friend void job_vec_add_val (struct thr_ctrl *);)
378template <typename T>
379void 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}
384INST(template <typename T> class Vector friend void job_val_add_vec (struct thr_ctrl *);)
386template <typename T>
387void 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}
392INST(template <typename T> class Vector friend void job_vec_sub_val (struct thr_ctrl *);)
394template <typename T>
395void 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}
400INST(template <typename T> class Vector friend void job_val_sub_vec (struct thr_ctrl *);)
402template <typename T>
403void 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}
408INST(template <typename T> class Vector friend void job_vec_mul_val (struct thr_ctrl *);)
410template <typename T>
411void 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}
416INST(template <typename T> class Vector friend void job_vec_div_val (struct thr_ctrl *);)
418template <typename T>
419void 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}
424INST(template <typename T> class Vector friend void job_val_div_vec (struct thr_ctrl *);)
426template <typename T>
427void 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 */
434INST(template <typename T> class Vector friend void job_vec_svc_add (struct thr_ctrl *);)
436template <typename T>
437void 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}
443INST(template <typename T> class Vector friend void job_svc_vec_add (struct thr_ctrl *);)
445template <typename T>
446void 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}
452INST(template <typename T> class Vector friend void job_svc_svc_add (struct thr_ctrl *);)
454template <typename T>
455void 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}
461
462INST(template <typename T> class Vector friend void job_vec_svc_sub (struct thr_ctrl *);)
463template <typename T>
464void 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}
470INST(template <typename T> class Vector friend void job_svc_vec_sub (struct thr_ctrl *);)
472template <typename T>
473void 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}
479INST(template <typename T> class Vector friend void job_svc_svc_sub (struct thr_ctrl *);)
481template <typename T>
482void 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}
488INST(template <typename T> class Vector friend void job_vec_add_svc (struct thr_ctrl *);)
490template <typename T>
491void 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}
497INST(template <typename T> class Vector friend void job_vec_sub_svc (struct thr_ctrl *);)
499template <typename T>
500void 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}
506INST(template <typename T> class Vector friend void job_vec_sub_svc_inv (struct thr_ctrl *);)
508template <typename T>
509void 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
516INST(template <typename T> class Vector friend void job_svc_val_add (struct thr_ctrl *);)
518template <typename T>
519void 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}
525INST(template <typename T> class Vector friend void job_svc_val_sub (struct thr_ctrl *);)
527template <typename T>
528void 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}
534INST(template <typename T> class Vector friend void job_val_svc_add (struct thr_ctrl *);)
536template <typename T>
537void 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}
543INST(template <typename T> class Vector friend void job_val_svc_sub (struct thr_ctrl *);)
545template <typename T>
546void 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}
552INST(template <typename T> class Vector friend void job_val_svc_div (struct thr_ctrl *);)
554template <typename T>
555void 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
616
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
665template <typename T>
666inline 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}
671
672template <typename T>
673inline 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}
678
679template <typename T>
680inline 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
687template <typename T>
688inline 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
696template <typename T>
697inline 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}
703
704template <typename T>
705inline 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
713template <typename T>
714inline 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}
722
723template <typename T>
724inline 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
736template <typename T>
737inline 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}
742template <typename T>
743inline 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
752template <typename T>
753const 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
763template <typename T>
764const 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
775INST(template<typename T> class TVector friend const TVector<T>& operator + (const T&, const TVector<T>&);)
776template <typename T>
777inline 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}
785
787INST(template<typename T> class TVector friend const TVector<T>& operator - (const T&, const TVector<T>&);)
788template <typename T>
789inline 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
799INST(template <typename T> class TVector friend TSVector<T> operator * (const T&, const TVector<T>&);)
800template <typename T>
801inline 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!
808INSTCTL(#ifndef VEC_INT_DIV_SUPPORT)
809INST(template <typename T> class TVector friend TSVector<T> operator / (const T&, TVector<T>);)
810INSTCTL(#endif)
811template <typename T>
812inline 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
820template <typename T>
821inline TVector<T>& TVector<T>::operator + (const T& a)
822{ return this->operator += (a); }
823template <typename T>
824inline TVector<T>& TVector<T>::operator - (const T& a)
825{ return this->operator -= (a); }
826
827
828template <typename T>
829inline 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!
835template <typename T>
836inline TSVector<T> TVector<T>::operator / (const T& a)
837{
838 return TSVector<T> (*this, (T)1/a);
839}
840
841#if 0
842template <typename T>
843inline 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
850template <typename T>
851inline 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>&);)
861template <typename T>
862inline 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>&);)
870template <typename T>
871inline 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>&);)
880template <typename T>
881inline 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>&);)
890template <typename T>
891inline 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>&);)
900template <typename T>
901inline 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>&);)
911template <typename T>
912inline 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
926template <typename T>
927inline 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
941template <typename T>
942inline 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
955template <typename T>
956inline 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
972template <typename T>
973inline 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);
982template <typename T> inline double fabssqr(const Vector<T>& v);
983
984NAMESPACE_END
985
986NAMESPACE_CSTD
987
988INST(template <typename T> class TBCI__ Vector friend double fabs (const TBCI__ Vector<T>&);)
989template <typename T>
990inline double fabs (const TBCI__ Vector<T>& v)
991{
992 return sqrt (TBCI__ fabssqr (v));
993}
994
995INST(template <typename T> class TBCI__ Vector friend T abs (const TBCI__ Vector<T>&);)
996template <typename T>
997inline T abs (const TBCI__ Vector<T>& v)
998{
999 return v.abs();
1000}
1001
1002INST(template <typename T> class TBCI__ TVector friend double fabs (TBCI__ TVector<T>);)
1003template <typename T>
1004inline double fabs (TBCI__ TVector<T> tv)
1005{
1006 TBCI__ Vector<T> v(tv); return fabs (v);
1007}
1008
1009INST(template <typename T> class TBCI__ Vector friend T abs (TBCI__ TVector<T>);)
1010template <typename T>
1011inline T abs (TBCI__ TVector<T> tv)
1012{
1013 TBCI__ Vector<T> v(tv); return abs (v);
1014}
1015
1016NAMESPACE_CSTD_END
1017
1018NAMESPACE_TBCI
1019
1025template <typename T>
1026class 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
1171template <typename T>
1172inline void TSVector<T>::destroy () const
1173{
1174 if (LIKELY(!mut || !dim)) return;
1175 TBCIDELETE(T, vec, dim);
1176}
1177
1178template <typename T>
1179inline 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
1194INST(template <typename T> class TSVector<T> friend const TSVector<T>& operator * (const T&, const TSVector<T>&);)
1195template <typename T>
1196const TSVector<T>& operator * (const T& f, const TSVector<T>& ts)
1197{ ts.getfacref() = f * ts.getfac(); return ts; }
1198
1199template <typename T>
1200inline 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
1220template <typename T>
1221const 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
1240template <typename T>
1241bool 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
1267template <typename T>
1268bool 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
1297template <typename T>
1298inline 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}
1306
1308template <typename T>
1309inline 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
1319template <typename T>
1320inline 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}
1333
1334template <typename T>
1335inline 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
1350template <typename T>
1351inline 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}
1368
1369template <typename T>
1370inline 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
1389template <typename T>
1390inline 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}
1402
1403template <typename T>
1404inline 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
1419template <typename T>
1420inline 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}
1433INST(template <typename T> class TSVector<T> friend TVector<T> operator + (const T&, const TSVector<T>&);)
1434template <typename T>
1435inline TVector<T> operator + (const T& v, const TSVector<T>& tsv)
1436{ return tsv.add_t_tsv (v); }
1437
1439template <typename T>
1440inline 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}
1453INST(template <typename T> class TSVector<T> friend TVector<T> operator - (const T&, const TSVector<T>&);)
1454template <typename T>
1455inline TVector<T> operator - (const T& v, const TSVector<T>& tsv)
1456{ return tsv.sub_t_tsv (v); }
1457
1458template <typename T>
1459STD__ 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
1467template <typename T>
1468inline 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
1477INST(template <typename T> class TSVector friend double fabssqr (const TSVector<T>&);)
1478template <typename T>
1479double fabssqr (const TSVector<T>& ts)
1480{ return ts.fabssqr (); }
1481
1482NAMESPACE_END
1483
1484NAMESPACE_CSTD
1485
1486INST(template <typename T> class TBCI__ TSVector friend double fabs (const TBCI__ TSVector<T>&);)
1487template <typename T>
1488double inline fabs (const TBCI__ TSVector<T>& ts)
1489{ return ts.fabs(); }
1490
1491INST(template <typename T> class TBCI__ TSVector friend T abs (const TBCI__ TSVector<T>&);)
1492template <typename T>
1493inline T abs (const TBCI__ TSVector<T>& ts)
1494{ return ts.abs(); }
1495
1496NAMESPACE_CSTD_END
1497
1498NAMESPACE_TBCI
1499
1500//------------------------------------------------------------
1501
1524
1525template <typename T>
1526class 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; }
1589
1590 Vector<T>& operator = (const Vector<T>& v) HOT
1591 { this->TVector<T>::operator = (v); return *this; }
1594 Vector<T>& operator = (const TVector<T>& tv) HOT
1595 { this->TVector<T>::operator = (tv); return *this; }
1596#ifdef TBCI_NEW_BRACKET
1598 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
1679INST(template <typename T> friend inline STD__ ostream& operator << (STD__ ostream& os, const Vector<T>& v);)
1680template <typename T>
1681inline STD__ ostream& operator << (STD__ ostream& os, const Vector<T>& v)
1682{ return os << (BVector<T>&)v; }
1683#endif
1684
1685#ifndef NO_POD
1686template <typename T>
1687inline 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
1709template <typename T>
1710inline 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
1717template <typename T>
1718bool 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
1745template <typename T>
1746inline 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
1753template <typename T>
1754inline 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
1765template <typename T>
1766inline 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
1782template <typename T>
1783inline 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
1800INST(template <typename T> class TVector friend void job_vec_dot (struct thr_ctrl*);)
1801template <typename T>
1802void 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
1816template <typename T>
1817inline 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
1887INST(template <typename T> class TVector friend void job_vec_mult (struct thr_ctrl*);)
1888template <typename T>
1889void 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
1901template <typename T>
1902inline 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
1967HOTDECL(template <typename T>
1968bool 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
1980HOTDECL(template <typename T>
1981void 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
1990INST(template <typename T> class TVector friend bool par_comp (const Vector<T>&, const Vector<T>&);)
1991INST(template <typename T> class TVector friend bool par_comp (const Vector<T>&, TVector<T>);)
1992INST(template <typename T> class TVector friend bool par_comp (TVector<T>, const Vector<T>&);)
1993INST(template <typename T> class TVector friend bool par_comp (TVector<T>, TVector<T>);)
1994
1995template <typename T>
1996inline 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
2002template <typename T>
2003inline bool par_comp (TVector<T> v1, const Vector<T>& v2)
2004{
2005 return par_comp ((const Vector<T>)v1, v2);
2006}
2007
2008template <typename T>
2009inline bool par_comp (TVector<T> v1, TVector<T> v2)
2010{
2011 return par_comp ((const Vector<T>)v1, (const Vector<T>)v2);
2012}
2013
2014template <typename T>
2015inline 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
2023template <typename T>
2024T 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
2031template <typename T>
2032T 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
2043template <typename T>
2044T 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
2057template <typename T>
2058T 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
2069template <typename T>
2070T 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
2084INST(template <typename T> class Vector friend void job_vec_sum(struct thr_ctrl*);)
2085template <typename T>
2086void 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
2095template <typename T>
2096T 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
2142template <typename T>
2143T 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
2168template <typename T>
2169inline 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}
2177
2178template <typename T>
2179inline 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
2190template <typename T>
2191inline 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}
2198
2200template <typename T>
2201inline 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
2212template <typename T>
2213inline 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}
2226
2227template <typename T>
2228inline 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
2243template <typename T>
2244inline 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}
2251
2252template <typename T>
2253inline 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
2262template <typename T>
2263inline TSVector<T> Vector<T>::operator * (const T& a) const
2264{
2265 return TSVector<T> (*this, a);
2266}
2267
2268template <typename T>
2269inline 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
2276INST(template <typename T> class Vector friend TVector<T> operator + (const T&, const Vector<T>&);)
2278template <typename T>
2279inline 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}
2290
2291INST(template <typename T> class Vector friend TVector<T> operator - (const T&, const Vector<T>&);)
2292template <typename T>
2293inline 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
2306INST(template <typename T> class Vector friend TSVector<T> operator * (const T&, const Vector<T>&);)
2307template <typename T>
2308inline 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
2314INSTCTL(#ifndef VEC_INT_DIV_SUPPORT)
2315INST(template <typename T> class Vector friend TSVector<T> operator / (const T&, const Vector<T>&);)
2316INSTCTL(#endif)
2317template <typename T>
2318TSVector<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
2325template <typename T>
2326inline 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
2334template <typename T>
2335inline 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
2347template <typename T>
2348inline 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
2363INST(template <typename T> class Vector friend double fabssqr (const Vector<T>&);)
2364template <typename T>
2365inline double fabssqr (const Vector<T>& v)
2366{ return v.fabssqr(); ;}
2367
2368
2369#if !defined(SMP) || defined(NOSMP_VECFABS)
2370
2371template <typename T>
2372inline 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
2387INST(template <typename T> class Vector friend void job_vec_fabssqr (struct thr_ctrl*);)
2388template <typename T>
2389void 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
2398template <typename T>
2399inline 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
2447INST(template <typename T> class TVector friend double fabssqr (TVector<T>);)
2448template <typename T>
2449inline double fabssqr (TVector<T> tv)
2450{ return tv.fabssqr (); }
2451
2452
2453/*
2454template <typename T>
2455Vector<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
2468NAMESPACE_END
2469
2470// ---------------------------------------------------------
2471// Specializations for int types
2472//----------------------------------------------------------
2473
2474#ifdef VEC_INT_DIV_SUPPORT
2475
2476NAMESPACE_TBCI
2477
2478//INST(template <> class TSVector friend TSVector<int> operator / (const int&, const TVector<int>&);)
2479template <>
2480inline 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
2489template <>
2490inline 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
2498template <> /* Was commented out. Why ??? */
2499inline 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
2509template <>
2510inline 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
2523template <>
2524inline 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
2535template <>
2536inline 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
2546template <>
2547inline 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
2557template <>
2558inline 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
2567NAMESPACE_END
2568
2569#endif /* VEC_INT_DIV_SUPPORT */
2570
2571#endif /* TBCI_VECTOR_H */
i
int i
Definition LM_fit.h:71
y
doublereal y
Definition TOMS_707.C:27
operator>>
STD__ istream & operator>>(STD__ istream &istr, BdMatrix< T > &mat)
Definition band_matrix.h:2739
CSTD__
#define CSTD__
Definition basics.h:340
STD__
#define STD__
Definition basics.h:338
HOT
#define HOT
Definition basics.h:495
BCHK
#define BCHK(cond, exc, txt, ind, rtval)
Definition basics.h:575
LIKELY
#define LIKELY(expr)
branch prediction note that we sometimes on purpose mark the unlikely possibility likely and vice ver...
Definition basics.h:100
FRIEND_TBCI__
#define FRIEND_TBCI__
Definition basics.h:334
NAMESPACE_CSTD
#define NAMESPACE_CSTD
Definition basics.h:319
MIN_ALIGN2
#define MIN_ALIGN2
Definition basics.h:424
FRIEND_TBCI2__
#define FRIEND_TBCI2__
Definition basics.h:335
NAMESPACE_END
#define NAMESPACE_END
Definition basics.h:323
TBCICOMP
#define TBCICOMP(n, o, t, s)
Definition basics.h:981
HOTDECL
#define HOTDECL(x)
Definition basics.h:497
INST
#define INST(x)
Definition basics.h:238
NAMESPACE_CSTD_END
#define NAMESPACE_CSTD_END
Definition basics.h:325
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
NAMESPACE_TBCI
#define NAMESPACE_TBCI
Definition basics.h:317
INSTCTL
#define INSTCTL(x)
Definition basics.h:245
TBCI__
#define TBCI__
Definition basics.h:332
UNLIKELY
#define UNLIKELY(expr)
Definition basics.h:101
REGISTER
#define REGISTER
Definition basics.h:108
MATH__
#define MATH__
Definition basics.h:339
NOINST
#define NOINST
Definition basics.h:244
GLBL__
#define GLBL__
Definition basics.h:342
ALIGN
#define ALIGN(x)
Definition basics.h:444
TBCICOPY
#define TBCICOPY(n, o, t, s)
Definition basics.h:895
GLBL2__
#define GLBL2__
Definition basics.h:343
MIN_ALIGN
#define MIN_ALIGN
Definition basics.h:421
FGD
#define FGD
Definition basics.h:144
T
#define T
Definition bdmatlib.cc:20
true
#define true
Definition bool.h:24
false
#define false
Definition bool.h:23
SMP_VECSLICE2
#define SMP_VECSLICE2
Definition bvector.h:626
BVector
provides basic Vector functionality but arithmetic operators (+=, - , *, /...).
Definition bvector.h:68
BVector::dim
unsigned long dim
Definition bvector.h:74
BVector::keep
bool keep
Definition bvector.h:75
BVector::operator>=
bool operator>=(const BVector< T > &bv) const
Definition bvector.h:520
BVector::end
iterator end()
Definition bvector.h:155
BVector::operator=
BVector< T > & operator=(const T &a)
Definition bvector.h:168
BVector::contains
bool contains(const T &, unsigned long *=0) const
Definition bvector.h:580
BVector::BVector
friend class BVector
Definition bvector.h:200
BVector::fill
BVector< T > & fill(const T &) HOT
Definition bvector.h:378
BVector::destroy
void destroy()
Definition bvector.h:268
BVector::resize
BVector< T > & resize(const BVector< T > &)
Actually it's a resize and copy (some people would expect the assignment op to do this).
Definition bvector.h:361
BVector::vec
T * vec
Definition bvector.h:73
BVector::operator==
bool operator==(const BVector< T > &) const HOT
KG, 2001-06-29: Strange: If we don't inline this, we seems to get better performance in our solver be...
Definition bvector.h:495
BVector::operator()
T & operator()(const unsigned long) HOT
Definition bvector.h:469
BVector::operator<=
bool operator<=(const BVector< T > &bv) const
Definition bvector.h:509
BdMatrix
The class BdMatrix is an implementation to store and do operations on sparse Matrices with a band str...
Definition band_matrix.h:104
BdMatrix::operator-=
BdMatrix< T > & operator-=(const BdMatrix< T > &)
Definition band_matrix.h:1315
BdMatrix::operator+=
BdMatrix< T > & operator+=(const BdMatrix< T > &)
Definition band_matrix.h:1314
BdMatrix::operator*=
BdMatrix< T > & operator*=(const T &)
Definition band_matrix.h:1334
BdMatrix::dim
unsigned int dim
Definition band_matrix.h:109
BdMatrix::destroy
void destroy()
destroy object explicitly
Definition band_matrix.h:156
CRMatrix
C++ class for sparse matrices using compressed row storage.
Definition crmatrix.h:64
CRMatrix::size
unsigned int size() const
Definition crmatrix.h:193
CRMatrix::destroy
void destroy()
Definition crmatrix.h:669
CSCMatrix
exception class: Use MatErr from matrix.h
Definition cscmatrix.h:66
F_Matrix
Definition f_matrix.h:1172
F_TMatrix
Temporary Base Class (non referable!) (acc.
Definition f_matrix.h:71
F_TSMatrix
Temporary object for scaled matrices.
Definition f_matrix.h:784
Mat_Brack
Definition matrix.h:2311
Symm_BdMatrix
Matrix class with optimized Matrix-Vector multiplication for symmetrical Matrices.
Definition symm_bdmatrix.h:64
TSVector::element_type
T element_type
Definition vector.h:1048
TSVector::detach
void detach(const T *=0) const
Definition vector.h:1179
TSVector::vec
T * vec
Definition vector.h:1029
TSVector::destroy
void destroy() const
Definition vector.h:1172
TSVector::operator==
bool operator==(const Vector< T > &) const
Definition vector.h:1241
TSVector::operator/=
TSVector< T > & operator/=(const T &f)
Definition vector.h:1092
TSVector::TSVector
TSVector(const TVector< T > &tv, const T &f=(T) 1)
Definition vector.h:1060
TSVector::operator=
TSVector< T > & operator=(const TSVector< T > &ts)
Definition vector.h:1078
TSVector::abs
T abs() const
Definition vector.h:1124
TSVector::~TSVector
~TSVector()
Definition vector.h:1054
TSVector::max
T max()
Definition vector.h:1123
TSVector::TSVector
TSVector()
Definition vector.h:1053
TSVector::operator/
const TSVector< T > & operator/(const T &f) const
Definition vector.h:1094
TSVector::dim
unsigned long dim
Definition vector.h:1029
TSVector::ALIGN
T fac ALIGN(MIN_ALIGN2)
TSVector::getcref
const T & getcref(const unsigned long i) const HOT
Definition vector.h:1075
TSVector::fabssqr
double fabssqr() const
Definition vector.h:1128
TSVector::min
T min()
Definition vector.h:1122
TSVector::get
T get(const unsigned long i) const HOT
Definition vector.h:1072
TSVector::clone
void clone(const bool=false, const T *=0) const
Definition vector.h:1200
TSVector::eval
const TSVector< T > & eval(const T *vv=0) const
Definition vector.h:1221
TSVector::operator*=
TSVector< T > & operator*=(const T &f)
Definition vector.h:1091
TSVector::vec_info
static const char * vec_info()
Definition vector.h:1165
TSVector::sub_t_tsv
TVector< T > sub_t_tsv(const T &) const
TV = T - TSV.
Definition vector.h:1440
TSVector::operator+
TVector< T > operator+(const Vector< T > &) const HOT
TV = TSV + V.
Definition vector.h:1320
TSVector::getfacref
T & getfacref() const
Definition vector.h:1074
TSVector::incr
TVector< T > incr(const unsigned long, const T=(T) 1) const
Definition vector.h:1468
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)
TSVector::operator()
T operator()(const unsigned long i) const
Definition vector.h:1066
TSVector::FGDT
friend NOINST double FRIEND_TBCI__ fabssqr FGDT(const TSVector< T > &)
TSVector::fabs
double fabs() const
Definition vector.h:1126
TSVector::sum
T sum()
Definition vector.h:1131
TSVector::mut
bool mut
Definition vector.h:1037
TSVector::vecptr
T *const & vecptr() const
Definition vector.h:1163
TSVector::getfac
const T & getfac() const
Definition vector.h:1073
TSVector::TSVector
TSVector(const Vector< T > &v, const T &f=(T) 1)
Definition vector.h:1062
TSVector::size
unsigned long size() const
Definition vector.h:1120
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
TSVector::operator[]
T operator[](const unsigned long i) const
Definition vector.h:1069
TSVector::TALIGN
T aligned_value_type TALIGN(MIN_ALIGN2)
Definition vector.h:1049
TSVector::TSVector
TSVector(const unsigned long d)
Definition vector.h:1055
TSVector::operator!=
bool operator!=(const Vector< T > &v) const
Definition vector.h:1156
TSVector::TSVector
TSVector(const TSVector< T > &ts)
Definition vector.h:1057
TSVector::value_type
T value_type
Definition vector.h:1047
TSVector::operator-
const TSVector< T > & operator-() const
Definition vector.h:1096
TSVector::operator*
const TSVector< T > & operator*(const T &f) const
Definition vector.h:1093
TVector
Temporary Base Class Idiom: Class TVector is used for temporary variables.
Definition vector.h:73
TVector::operator/
TSVector< T > operator/(const T &)
Definition vector.h:836
TVector::FGD
friend TVector< T > conj FGD(const Vector< T > &)
TVector::get
tbci_traits< T >::const_refval_type get(const unsigned long i) const
Definition vector.h:190
TVector::FGD
friend TVector< T > imag FGD(const Vector< T > &)
TVector::TVector
TVector(const T &val, const unsigned long d)
Definition vector.h:85
TVector::getcref
const T & getcref(const unsigned long i) const
Definition vector.h:188
TVector::operator==
bool operator==(const TVector< T > &tv) const
Definition vector.h:114
TVector::operator+=
TVector< T > & operator+=(const T &)
TV += a.
Definition vector.h:666
TVector::min
T min()
Definition vector.h:146
TVector::FGD
friend TVector< T > &real FGD(TVector< T > &)
TVector::operator-=
TVector< T > & operator-=(const T &)
TV -= a.
Definition vector.h:673
TVector::size
unsigned long size() const
Definition vector.h:104
TVector::slice
TVector< T > slice(const unsigned long, const unsigned long)
Definition vector.h:956
TVector::value_type
T value_type
Definition vector.h:80
TVector::TVector
TVector(const BVector< T > &bv)
Definition vector.h:86
TVector::operator*
TSVector< T > operator*(const T &)
Definition vector.h:829
TVector::operator/=
TVector< T > & operator/=(const T &)
TV /= a.
Definition vector.h:688
TVector::set
T & set(const T &val, const unsigned long i) const
Definition vector.h:194
TVector::operator-
TSVector< T > operator-()
Definition vector.h:851
TVector::operator!=
bool operator!=(const TVector< T > &tv) const
Definition vector.h:118
TVector::TALIGN
T aligned_value_type TALIGN(MIN_ALIGN2)
Definition vector.h:82
TVector::TVector
TVector(const TSVector< T > &ts)
Definition vector.h:91
TVector::FGD
friend TVector< T > &imag FGD(TVector< T > &)
TVector::incr
TVector< T > & incr(const unsigned long, const T=(T) 1)
Definition vector.h:973
TVector::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::element_type
T element_type
Definition vector.h:81
TVector::vec_info
static const char * vec_info()
Definition vector.h:219
TVector::max
T max()
Definition vector.h:147
TVector::FGD
friend TVector< T > &conj FGD(TVector< T > &)
TVector::setval
T & setval(const T &val, const unsigned long i) const
Definition vector.h:193
TVector::operator=
TVector< T > & operator=(const T &a)
Definition vector.h:106
TVector::operator[]
T operator[](const unsigned long i)
Definition vector.h:185
TVector::fabs
double fabs()
Definition vector.h:151
TVector::operator*=
TVector< T > & operator*=(const T &)
TV *= a.
Definition vector.h:680
TVector::TVector
TVector(const unsigned long d=0)
Definition vector.h:84
TVector::FGD
friend NOINST void FRIEND_TBCI2__ do_mat_vec_mult FGD(const unsigned start, const unsigned end, TVector< T > *res, const Matrix< T > *mat, const Vector< T > *vec)
TVector::FGD
friend NOINST void FRIEND_TBCI2__ do_mat_vec_transmult FGD(const unsigned start, const unsigned end, TVector< T > *res, const Matrix< T > *mat, const Vector< T > *vec)
TVector::setval
T & setval(const unsigned long i) const
Definition vector.h:192
TVector::TVector
TVector(const TVector< T > &tv) HOT
Definition vector.h:89
TVector::TVector
TVector(const BVector< U > &bv)
Definition vector.h:228
TVector::fabssqr
double fabssqr()
Definition vector.h:150
TVector::~TVector
~TVector()
Definition vector.h:95
TVector::FGD
friend TVector< T > real FGD(const Vector< T > &)
TVector::operator+
TVector< T > & operator+(const Vector< T > &) HOT
Definition vector.h:737
TVector::sum
T sum()
Definition vector.h:148
TVector::TVector
TVector(const Vector< T > &v)
Definition vector.h:87
TVector::abs
T abs()
Definition vector.h:152
TVector::operator()
T operator()(const unsigned long i)
Definition vector.h:182
VecErr
exception class
Definition bvector.h:32
Vector_Sig::Vector_Sig
Vector_Sig()
Definition vector_sig.h:76
Vector::TALIGN
T aligned_value_type TALIGN(MIN_ALIGN2)
Definition vector.h:1532
Vector::operator>
bool operator>(const Vector< T > &v) const
Definition vector.h:1614
Vector::Vector
Vector(const BVector< T > &bv)
Definition vector.h:1541
Vector::operator()
tbci_traits< T >::const_refval_type operator()(const unsigned long i) const
Definition vector.h:1563
Vector::operator>=
bool operator>=(const Vector< T > &v) const
Definition vector.h:1612
Vector::vec_info
static const char * vec_info()
Definition vector.h:1655
Vector::operator/
TSVector< T > operator/(const T &) const
Definition vector.h:2269
Vector::operator!=
bool operator!=(const Vector< T > &v) const
Definition vector.h:1605
Vector::Vector
Vector(const TSVector< T > &ts)
Definition vector.h:1544
Vector::Vector
Vector(const TVector< T > &tv)
Definition vector.h:1546
Vector::sum
T sum() const
Definition vector.h:2096
Vector::operator[]
tbci_traits< T >::const_refval_type operator[](const unsigned long i) const
Definition vector.h:1566
Vector::operator=
Vector< T > & operator=(const T &a)
Definition vector.h:1587
Vector::min
T min() const
Definition vector.h:2032
Vector::operator==
bool operator==(const Vector< T > &v) const
Definition vector.h:1603
Vector::operator*
TSVector< T > operator*(const T &) const
Definition vector.h:2263
Vector::value_type
T value_type
Definition vector.h:1530
Vector::abs
T abs() const
Definition vector.h:1646
Vector::Vector
Vector(const T &val, const unsigned long d)
Definition vector.h:1540
Vector::fabs
double fabs() const
Definition vector.h:1645
Vector::element_type
T element_type
Definition vector.h:1531
Vector::slice
TVector< T > slice(unsigned long, unsigned long) const
Definition vector.h:2335
Vector::Vector
Vector(const Vector< T > &v)
Definition vector.h:1542
Vector::fabssqr
double fabssqr() const HOT
Definition vector.h:2399
Vector::fill
Vector< T > & fill(const T &v)
Definition vector.h:1577
Vector::FGD
friend T dot FGD(const Vector< T > &, const Vector< T > &) HOT
Vector::operator-
TVector< T > operator-() const
Definition vector.h:2326
Vector::max
T max() const
Definition vector.h:2058
Vector::incr
TVector< T > incr(const unsigned long, const T=(T) 1) const
Definition vector.h:2348
Vector::Vector
Vector(const unsigned long d=0)
Definition vector.h:1539
Vector::contains
bool contains(const T &v) const
Definition vector.h:1649
Vector::Vector
Vector(const BVector< U > &bv)
Definition vector.h:1673
Vector::clear
Vector< T > & clear()
Definition vector.h:1578
Vector::operator<
bool operator<(const Vector< T > &v) const
Definition vector.h:1613
Vector::~Vector
~Vector()
Definition vector.h:1559
Vector::operator<=
bool operator<=(const Vector< T > &v) const
Definition vector.h:1611
Vector::operator+
TVector< T > operator+(const Vector< T > &) const HOT
TV = V + V.
Definition vector.h:2169
cplx
Our own complex class.
Definition cplx.h:56
imag
T imag(const TBCI__ cplx< T > &z)
Definition cplx.h:674
fabs
double fabs(const TBCI__ cplx< T > &c)
Definition cplx.h:746
dot
cplx< T > dot(const cplx< T > &a, const cplx< T > &b)
Definition cplx.h:300
fabssqr
double fabssqr(const cplx< T > &c)
Definition cplx.h:390
conj
NAMESPACE_END NAMESPACE_CPLX TBCI__ cplx< T > conj(const TBCI__ cplx< T > &c)
Definition cplx.h:663
norm
double norm(const TBCI__ cplx< T > &c)
Definition cplx.h:695
sqrt
NAMESPACE_END NAMESPACE_CSTD TBCI__ cplx< T > sqrt(const TBCI__ cplx< T > &z)
Definition cplx.h:751
abs
#define abs(x)
Definition f2c.h:178
ts
F_TSMatrix< T > ts
Definition f_matrix.h:1052
fac
tm fac
Definition f_matrix.h:1052
b
F_TMatrix< T > b
Definition f_matrix.h:736
detach
tm detach()
sum
T sum(const FS_Vector< dims, T > &fv)
Definition fs_vector.h:599
dims
#define dims
Definition fsveclib.cc:24
TBCIDELETE
#define TBCIDELETE(t, v, sz)
Definition malloc_cache.h:634
NEW
#define NEW(t, s)
Definition malloc_cache.h:633
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)
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)
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)
a
const unsigned TMatrix< T > const Matrix< T > * a
Definition matrix_kernels.h:26
end
const unsigned end
Definition matrix_kernels.h:25
res
const unsigned TMatrix< T > * res
Definition matrix_kernels.h:26
CPLX__
real
#define real
Definition own_lapack_routines.cpp:20
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
thread_start
void thread_start(const int thr_no, thr_job_t job, const unsigned long sz,...)
Definition smp.cc:988
threads
struct thr_struct * threads
Definition smp.cc:106
thr_job_t
void(* thr_job_t)(struct thr_ctrl *)
Before the double inclusion guard on purpose!
Definition smp.h:126
THREAD_MAX_RES_LN
#define THREAD_MAX_RES_LN
Definition smp.h:130
threads_avail
#define threads_avail(x)
Definition smp.h:322
job_output
Definition smp.h:141
job_output::t_res_ptr
void * t_res_ptr
Definition smp.h:149
thr_ctrl
Definition smp.h:168
thr_ctrl::t_res_d
double t_res_d
Definition smp.h:178
thr_ctrl::t_res_ptr
void * t_res_ptr
Definition smp.h:180
thr_ctrl::t_size
unsigned long t_size
Definition smp.h:171
thr_ctrl::t_par
void * t_par[6]
Definition smp.h:173
do_exactsum
unsigned int do_exactsum()
Definition tbci_param.h:150
VEC_INT_DIV_SUPPORT
#define VEC_INT_DIV_SUPPORT
Definition veclib.cc:12
STD_SMP_TEMPLATE5
#define STD_SMP_TEMPLATE5(oper, dm, a1, a2, a3, a4, a5)
Definition vector.h:658
job_vec_sum
void job_vec_sum(struct thr_ctrl *tc)
Definition vector.h:2086
job_svc_svc_add
void job_svc_svc_add(struct thr_ctrl *tc)
vec = s*vec + s*vec;
Definition vector.h:455
fabssqr
double fabssqr(const Vector< T > &v)
Definition vector.h:2365
STD_SMP_TEMPLATE2V
#define STD_SMP_TEMPLATE2V(oper, dm, a1, a2)
Definition vector.h:644
job_vec_sub_vec
void job_vec_sub_vec(struct thr_ctrl *tc)
vec -= vec;
Definition vector.h:307
job_vec_add_val
void job_vec_add_val(struct thr_ctrl *tc)
vec += val;
Definition vector.h:379
par_fill
void par_fill(Vector< T > &vec, vec_fill_fn< T > fn, void *par)
Definition vector.h:2015
job_vec_div_val
void job_vec_div_val(struct thr_ctrl *tc)
vec /= val;
Definition vector.h:419
fabs
NAMESPACE_END NAMESPACE_CSTD double fabs(const TBCI__ Vector< T > &v)
Definition vector.h:990
job_svc_val_sub
void job_svc_val_sub(struct thr_ctrl *tc)
vec = s*vec - val;
Definition vector.h:528
job_svc_val_add
void job_svc_val_add(struct thr_ctrl *tc)
vec = s*vec + val;
Definition vector.h:519
STD_SMP_TEMPLATE4V
#define STD_SMP_TEMPLATE4V(oper, dm, a1, a2, a3, a4)
Definition vector.h:654
job_vec_mult
void job_vec_mult(struct thr_ctrl *tc)
Definition vector.h:1889
job_val_vec_mul
void job_val_vec_mul(struct thr_ctrl *tc)
vec = val * vec;
Definition vector.h:347
job_val_div_vec
void job_val_div_vec(struct thr_ctrl *tc)
vec = val/self;
Definition vector.h:427
job_vec_add_svc
void job_vec_add_svc(struct thr_ctrl *tc)
vec += s*vec;
Definition vector.h:491
job_vec_sub_vec_inv
void job_vec_sub_vec_inv(struct thr_ctrl *tc)
vec -= vec; vec = -vec;
Definition vector.h:315
par_comp
bool par_comp(const Vector< T > &v1, TVector< T > v2)
Definition vector.h:1996
job_val_vec_add
void job_val_vec_add(struct thr_ctrl *tc)
vec = val + vec;
Definition vector.h:355
job_vec_sub_val
void job_vec_sub_val(struct thr_ctrl *tc)
vec -= val;
Definition vector.h:395
job_svc_vec_add
void job_svc_vec_add(struct thr_ctrl *tc)
vec = s*vec + vec;
Definition vector.h:446
job_val_svc_add
void job_val_svc_add(struct thr_ctrl *tc)
vec = val + s*vec;
Definition vector.h:537
operator<<
STD__ ostream & operator<<(STD__ ostream &os, const TVector< T > &tv)
Definition vector.h:253
job_vec_sub_svc
void job_vec_sub_svc(struct thr_ctrl *tc)
vec -= s*vec;
Definition vector.h:500
job_val_svc_div
void job_val_svc_div(struct thr_ctrl *tc)
vec = val / s*vec;
Definition vector.h:555
job_vec_val_sub
void job_vec_val_sub(struct thr_ctrl *tc)
vec = vec - val;
Definition vector.h:331
job_vec_svc_add
void job_vec_svc_add(struct thr_ctrl *tc)
vec = vec + s*vec;
Definition vector.h:437
job_vec_add_vec
void job_vec_add_vec(struct thr_ctrl *tc)
vec += vec;
Definition vector.h:299
job_vec_sub_svc_inv
void job_vec_sub_svc_inv(struct thr_ctrl *tc)
vec -= s*vec;
Definition vector.h:509
imag
TVector< T > & imag(TVector< T > &tv)
Definition vector.h:881
dot
T dot(const Vector< T > &a, const Vector< T > &b)
Definition vector.h:1817
job_val_svc_sub
void job_val_svc_sub(struct thr_ctrl *tc)
vec = val + s*vec;
Definition vector.h:546
job_vec_vec_sub
void job_vec_vec_sub(struct thr_ctrl *tc)
vec = vec - vec;
Definition vector.h:291
STD_SMP_TEMPLATE3VV
#define STD_SMP_TEMPLATE3VV(oper, dm, a1, a2, a3)
Definition vector.h:648
real
TVector< T > & real(TVector< T > &tv)
Definition vector.h:871
job_val_vec_div
void job_val_vec_div(struct thr_ctrl *tc)
vec = val / vec;
Definition vector.h:371
job_val_sub_vec
void job_val_sub_vec(struct thr_ctrl *tc)
vec -= val; vec = -vec;
Definition vector.h:403
STD_SMP_TEMPLATE3VC
#define STD_SMP_TEMPLATE3VC(oper, dm, a1, a2, a3)
Definition vector.h:650
STD_SMP_TEMPLATE2C
#define STD_SMP_TEMPLATE2C(oper, dm, a1, a2)
Definition vector.h:646
job_vec_svc_sub
void job_vec_svc_sub(struct thr_ctrl *tc)
vec = vec - s*vec;
Definition vector.h:464
job_vec_val_mul
void job_vec_val_mul(struct thr_ctrl *tc)
vec = vec * val;
Definition vector.h:339
conj
TVector< T > & conj(TVector< T > &tv)
Definition vector.h:862
job_val_vec_sub
void job_val_vec_sub(struct thr_ctrl *tc)
vec = val - vec;
Definition vector.h:363
job_svc_vec_sub
void job_svc_vec_sub(struct thr_ctrl *tc)
vec = s*vec - vec;
Definition vector.h:473
job_vec_val_add
void job_vec_val_add(struct thr_ctrl *tc)
vec = vec + val;
Definition vector.h:323
STD_SMP_TEMPLATE4C
#define STD_SMP_TEMPLATE4C(oper, dm, a1, a2, a3, a4)
Definition vector.h:656
job_vec_vec_add
void job_vec_vec_add(struct thr_ctrl *tc)
vec = vec + vec;
Definition vector.h:283
job_vec_dot
void job_vec_dot(struct thr_ctrl *tc)
Definition vector.h:1802
job_val_add_vec
void job_val_add_vec(struct thr_ctrl *tc)
vec += val;
Definition vector.h:387
operator/
TSVector< T > operator/(const T &a, TVector< T > b)
Definition vector.h:812
job_vec_fabssqr
void job_vec_fabssqr(struct thr_ctrl *tc)
Definition vector.h:2389
job_svc_svc_sub
void job_svc_svc_sub(struct thr_ctrl *tc)
vec = s*vec - s*vec;
Definition vector.h:482
job_vec_mul_val
void job_vec_mul_val(struct thr_ctrl *tc)
vec *= val;
Definition vector.h:411
vector_extra.h
contains rarely used functions on TBCI::Vector s
Generated by doxygen 1.16.1