TBCI Numerical high perf. C++ Library  2.8.0
vector_extra.h
Go to the documentation of this file.
1 
6 /* KG: 1998/08/17 */
7 /* $Id: vector_extra.h,v 1.9.2.17 2022/11/03 17:28:12 garloff Exp $ */
8 
9 #ifndef TBCI_VECTOR_EXTRA_H
10 #define TBCI_VECTOR_EXTRA_H
11 
12 #include "tbci/vector.h"
13 #include "tbci/matrix.h"
14 
15 // avoid -fguiding-decls (here: unnecessary)
16 #if !defined(NO_GD) && !defined(AUTO_DECL)
17 # include "tbci/vector_extra_gd.h"
18 #endif
19 
20 NAMESPACE_TBCI
21 
22 // elementwise multiplication
23 INST(template <typename T> class TVector friend TVector<T> emul (const Vector<T>&, const Vector<T>&);)
24 template <typename T>
25 inline TVector<T> emul (const Vector<T>& a, const Vector<T>& b)
26 {
27  const unsigned long dim = b.size();
28  TVector<T> s (dim);
29  BCHK(a.size() != dim, VecErr, product with diff. size Vectors, dim, s);
30  do_vec_vec_mul<T> (dim, &s.setval(0), &a.getcref(0), &b.getcref(0));
31  return s;
32 }
33 
34 INST(template <typename T> class TVector friend TVector<T> emul (TVector<T>, const Vector<T>&);)
35 template <typename T>
36 inline TVector<T> emul (TVector<T> a, const Vector<T>& b)
37 {
38  const unsigned long dim = b.size();
39  BCHK(a.size() != dim, VecErr, product with diff. size Vectors, dim, a);
40  do_vec_mul_vec<T> (dim, &a.setval(0), &b.getcref(0));
41  return a;
42 }
43 
44 INST(template <typename T> class TVector friend TVector<T> emul (const Vector<T>&, TVector<T>);)
45 template <typename T>
46 inline TVector<T> emul (const Vector<T>& a, TVector<T> b)
47 {
48  const unsigned long dim = b.size();
49  BCHK(a.size() != dim, VecErr, product with diff. size Vectors, dim, b);
50  do_vec_mul_vec<T> (dim, &b.setval(0), &a.getcref(0));
51  return b;
52 }
53 
54 INST(template <typename T> class TVector friend TVector<T> emul (TVector<T>, TVector<T>);)
55 template <typename T>
56 inline TVector<T> emul (TVector<T> a, TVector<T> b)
57 {
58  const unsigned long dim = b.size();
59  BCHK(a.size() != dim, VecErr, product with diff. size Vectors, dim, a);
60  do_vec_mul_vec<T> (dim, &a.setval(0), &b.getcref(0));
61  b.destroy (); return a;
62 }
63 
64 INST(template <typename T> class TVector friend TVector<T> cemul (const Vector<T>&, const Vector<T>&);)
65 template <typename T>
66 inline TVector<T> cemul (const Vector<T>& a, const Vector<T>& b)
67 {
68  const unsigned long dim = b.size();
69  TVector<T> s(dim);
70  BCHK(a.size() != dim, VecErr, product with diff. size Vectors, dim, s);
71  do_vec_vec_cmul<T> (dim, &s.setval(0), &a.getcref(0), &b.getcref(0));
72  return s;
73 }
74 
75 INST(template <typename T> class TVector friend TVector<T> cemul (TVector<T>, const Vector<T>&);)
76 template <typename T>
77 inline TVector<T> cemul (TVector<T> a, const Vector<T>& b)
78 {
79  const unsigned long dim = b.size();
80  BCHK(a.size() != dim, VecErr, product with diff. size Vectors, dim, a);
81  do_vec_cmul_vec<T> (dim, &a.setval(0), &b.getcref(0));
82  return a;
83 }
84 
85 INST(template <typename T> class TVector friend TVector<T> cemul (const Vector<T>&, TVector<T>);)
86 template <typename T>
87 inline TVector<T> cemul (const Vector<T>& a, TVector<T> b)
88 {
89  const unsigned long dim = b.size();
90  BCHK(a.size() != dim, VecErr, product with diff. size Vectors, dim, b);
91  do_vec_cmul_vec_inv<T> (dim, &b.setval(0), &a.getcref(0));
92  return b;
93 }
94 
95 INST(template <typename T> class TVector friend TVector<T> cemul (TVector<T>, TVector<T>);)
96 template <typename T>
97 inline TVector<T> cemul (TVector<T> a, TVector<T> b)
98 {
99  const unsigned long dim = b.size();
100  BCHK(a.size() != dim, VecErr, product with diff. size Vectors, dim, a);
101  do_vec_cmul_vec<T> (dim, &a.setval(0), &b.getcref(0));
102  b.destroy (); return a;
103 }
104 
105 
106 // elementwise division
107 
108 INST(template <typename T> class TVector friend TVector<T> ediv (const Vector<T>&, const Vector<T>&);)
109 template <typename T>
110 inline TVector<T> ediv (const Vector<T>& a, const Vector<T>& b)
111 {
112  TVector<T> s (b.size());
113  BCHK(a.size() != b.size(), VecErr, division with diff. size Vectors, b.size(), s);
114  do_vec_vec_div<T> (b.size(), &s.setval(0), &a.getcref(0), &b.getcref(0));
115  return s;
116 }
117 
118 INST(template <typename T> class TVector friend TVector<T> ediv (TVector<T>, const Vector<T>&);)
119 template <typename T>
120 inline TVector<T> ediv (TVector<T> a, const Vector<T>& b)
121 {
122  BCHK(a.size() != b.size(), VecErr, division with diff. size Vectors, b.size(), a);
123  do_vec_div_vec<T> (b.size(), &a.setval(0), &b.getcref(0));
124  return a;
125 }
126 
127 INST(template <typename T> class TVector friend TVector<T> ediv (const Vector<T>&, TVector<T>);)
128 template <typename T>
129 inline TVector<T> ediv (const Vector<T>& a, TVector<T> b)
130 {
131  BCHK(a.size() != b.size(), VecErr, division with diff. size Vectors, b.size(), b);
132  do_vec_div_vec_inv<T> (b.size(), &b.setval(0), &a.getcref(0));
133  return b;
134 }
135 
136 INST(template <typename T> class TVector friend TVector<T> ediv (TVector<T>, TVector<T>);)
137 template <typename T>
138 inline TVector<T> ediv (TVector<T> a, TVector<T> b)
139 {
140  BCHK(a.size() != b.size(), VecErr, division with diff. size Vectors, b.size(), a);
141  do_vec_div_vec<T> (b.size(), &a.setval(0), &b.getcref(0));
142  b.destroy (); return a;
143 }
144 
145 INST(template <typename T> class TVector friend TVector<T> cediv (const Vector<T>&, const Vector<T>&);)
146 template <typename T>
147 inline TVector<T> cediv (const Vector<T>& a, const Vector<T>& b)
148 {
149  TVector<T> s (b.size());
150  BCHK(a.size() != b.size(), VecErr, division with diff. size Vectors, b.size(), s);
151  do_vec_vec_cdiv<T> (b.size(), &s.setval(0), &a.getcref(0), &b.getcref(0));
152  return s;
153 }
154 
155 INST(template <typename T> class TVector friend TVector<T> cediv (TVector<T>, const Vector<T>&);)
156 template <typename T>
157 inline TVector<T> cediv (TVector<T> a, const Vector<T>& b)
158 {
159  BCHK(a.size() != b.size(), VecErr, division with diff. size Vectors, b.size(), a);
160  do_vec_cdiv_vec<T> (b.size(), &a.setval(0), &b.getcref(0));
161  return a;
162 }
163 
164 INST(template <typename T> class TVector friend TVector<T> cediv (const Vector<T>&, TVector<T>);)
165 template <typename T>
166 inline TVector<T> cediv (const Vector<T>& a, TVector<T> b)
167 {
168  BCHK(a.size() != b.size(), VecErr, division with diff. size Vectors, b.size(), b);
169  do_vec_cdiv_vec_inv<T> (a.size(), &b.setval(0), &a.getcref(0));
170  return b;
171 }
172 
173 INST(template <typename T> class TVector friend TVector<T> cediv (TVector<T>, TVector<T>);)
174 template <typename T>
175 inline TVector<T> cediv (TVector<T> a, TVector<T> b)
176 {
177  BCHK(a.size() != b.size(), VecErr, division with diff. size Vectors, b.size(), a);
178  do_vec_cdiv_vec<T> (b.size(), &a.setval(0), &b.getcref(0));
179  b.destroy (); return a;
180 }
181 
182 
183 // Now for the combinations with TSVector<T>
184 // Not yet optimized with the unrolled kernels from vec_kern_unr_pref.h
185 
186 INST(template <typename T> class TSVector friend TVector<T> emul (TVector<T>, const TSVector<T>&);)
187 template <typename T>
188 TVector<T> emul (TVector<T> tv, const TSVector<T>& ts)
189 {
190  BCHK(tv.size() != ts.size(), VecErr, Wrong dim in emul, ts.size(), tv);
191  for (REGISTER unsigned long i = 0; i < tv.size(); i++)
192  tv.setval(i) *= ts.get(i);
193  ts.destroy (); return tv;
194 }
195 
196 INST(template <typename T> class TSVector friend TVector<T> emul (const TSVector<T>&, TVector<T>);)
197 template <typename T>
198 TVector<T> emul (const TSVector<T>& ts, TVector<T> tv)
199 {
200  BCHK(tv.size() != ts.size(), VecErr, Wrong dim in emul, tv.size(), tv);
201  for (REGISTER unsigned long i = 0; i < tv.size(); i++)
202  tv.set (ts.get(i) * tv.get(i), i);
203  ts.destroy (); return tv;
204 }
205 
206 
207 INST(template <typename T> class TSVector friend TVector<T> emul (const Vector<T>&, const TSVector<T>&);)
208 template <typename T>
209 TVector<T> emul (const Vector<T>& v, const TSVector<T>& ts)
210 {
211  TVector<T> tv;
212  if (ts.mut) { tv.setptr(ts.vecptr()); tv.setsize(ts.size()); }
213  else tv.resize (ts.size());
214  BCHK(v.size() != ts.size(), VecErr, Wrong dim in emul, ts.size(), tv);
215  for (REGISTER unsigned long i = 0; i < tv.size(); i++)
216  tv.set (v.get(i) * ts.get(i), i);
217  return tv;
218 }
219 
220 INST(template <typename T> class TSVector friend TVector<T> emul (const TSVector<T>&, const Vector<T>&);)
221 template <typename T>
222 TVector<T> emul (const TSVector<T>& ts, const Vector<T>& v)
223 {
224  TVector<T> tv;
225  if (ts.mut) { tv.setptr(ts.vecptr()); tv.setsize(ts.size()); }
226  else tv.resize (ts.size());
227  BCHK(v.size() != ts.size(), VecErr, Wrong dim in emul, v.size(), tv);
228  for (REGISTER unsigned long i = 0; i < tv.size(); i++)
229  tv.set (ts.get(i) * v.get(i), i);
230  return tv;
231 }
232 
233 INST(template <typename T> class TSVector friend TVector<T> emul (const TSVector<T>&, const TSVector<T>&);)
234 template <typename T>
235 TVector<T> emul (const TSVector<T>& ts1, const TSVector<T>& ts2)
236 {
237  TVector<T> tv;
238  if (ts1.mut) { tv.setptr(ts1.vecptr()); tv.setsize(ts1.size()); }
239  else if (ts2.mut) { tv.setptr(ts2.vecptr()); tv.setsize(ts2.size()); }
240  else tv.resize (ts1.size());
241  BCHK(ts1.size() != ts2.size(), VecErr, Wrong dim in emul, ts2.size(), tv);
242  for (REGISTER unsigned long i = 0; i < tv.size(); i++)
243  tv.set (ts1.get(i), i);
244  if (ts1.mut && ts2.mut) ts2.destroy ();
245  return tv;
246 }
247 
248 
249 INST(template <typename T> class TSVector friend TVector<T> cemul (TVector<T>, const TSVector<T>&);)
250 template <typename T>
251 TVector<T> cemul (TVector<T> tv, const TSVector<T>& ts)
252 {
253  BCHK(tv.size() != ts.size(), VecErr, Wrong dim in emul, ts.size(), tv);
254  for (REGISTER unsigned long i = 0; i < tv.size(); i++)
255  tv.set (CPLX__ conj(tv.get(i)) * ts.get(i), i);
256  ts.destroy (); return tv;
257 }
258 
259 INST(template <typename T> class TSVector friend TVector<T> cemul (const TSVector<T>&, TVector<T>);)
260 template <typename T>
261 TVector<T> cemul (const TSVector<T>& ts, TVector<T> tv)
262 {
263  BCHK(tv.size() != ts.size(), VecErr, Wrong dim in emul, tv.size(), tv);
264  for (REGISTER unsigned long i = 0; i < tv.size(); i++)
265  tv.set (CPLX__ conj (ts.get(i)) * tv.get(i), i);
266  ts.destroy (); return tv;
267 }
268 
269 INST(template <typename T> class TSVector friend TVector<T> cemul (const Vector<T>&, const TSVector<T>&);)
270 template <typename T>
271 TVector<T> cemul (const Vector<T>& v, const TSVector<T>& ts)
272 {
273  TVector<T> tv;
274  if (ts.mut) { tv.setptr(ts.vecptr()); tv.setsize(ts.size()); }
275  else tv.resize (ts.size());
276  BCHK(v.size() != ts.size(), VecErr, Wrong dim in emul, ts.size(), tv);
277  for (REGISTER unsigned long i = 0; i < tv.size(); i++)
278  tv.set (CPLX__ conj(v.get(i)) * ts.get(i), i);
279  return tv;
280 }
281 
282 
283 INST(template <typename T> class TSVector friend TVector<T> cemul (const TSVector<T>&, const Vector<T>&);)
284 template <typename T>
285 TVector<T> cemul (const TSVector<T>& ts, const Vector<T>& v)
286 {
287  TVector<T> tv;
288  if (ts.mut) { tv.setptr(ts.vecptr()); tv.setsize(ts.size()); }
289  else tv.resize (ts.size());
290  BCHK(v.size() != ts.size(), VecErr, Wrong dim in emul, v.size(), tv);
291  for (REGISTER unsigned long i = 0; i < tv.size(); i++)
292  tv.set (CPLX__ conj (ts.get(i)) * v.get(i), i);
293  return tv;
294 }
295 
296 INST(template <typename T> class TSVector friend TVector<T> cemul (const TSVector<T>&, const TSVector<T>&);)
297 template <typename T>
298 TVector<T> cemul (const TSVector<T>& ts1, const TSVector<T>& ts2)
299 {
300  TVector<T> tv;
301  if (ts1.mut) { tv.setptr(ts1.vecptr()); tv.setsize(ts1.size()); }
302  else if (ts2.mut) { tv.setptr(ts2.vecptr()); tv.setsize(ts2.size()); }
303  else tv.resize (ts1.size());
304  BCHK(ts1.size() != ts2.size(), VecErr, Wrong dim in emul, ts2.size(), tv);
305  for (REGISTER unsigned long i = 0; i < tv.size(); i++)
306  tv.set (CPLX__ conj (ts1.get(i)), i);
307  if (ts1.mut && ts2.mut) ts2.destroy ();
308  return tv;
309 }
310 
311 
312 INST(template <typename T> class TSVector friend TVector<T> ediv (TVector<T>, const TSVector<T>&);)
313 template <typename T>
314 TVector<T> ediv (TVector<T> tv, const TSVector<T>& ts)
315 {
316  BCHK(tv.size() != ts.size(), VecErr, Wrong dim in ediv, ts.size(), tv);
317  for (REGISTER unsigned long i = 0; i < tv.size(); i++)
318  tv.setval(i) /= ts.get(i);
319  ts.destroy (); return tv;
320 }
321 
322 INST(template <typename T> class TSVector friend TVector<T> ediv (const Vector<T>&, const TSVector<T>&);)
323 template <typename T>
324 TVector<T> ediv (const Vector<T>& v, const TSVector<T>& ts)
325 {
326  TVector<T> tv;
327  if (ts.mut) { tv.setptr(ts.vecptr()); tv.setsize(ts.size()); }
328  else tv.resize (ts.size());
329  BCHK(v.size() != ts.size(), VecErr, Wrong dim in ediv, ts.size(), tv);
330  for (REGISTER unsigned long i = 0; i < tv.size(); i++)
331  tv.set (v.get(i) / ts.get(i), i);
332  return tv;
333 }
334 
335 INST(template <typename T> class TSVector friend TVector<T> ediv (const TSVector<T>&, const TSVector<T>&);)
336 template <typename T>
337 TVector<T> ediv (const TSVector<T>& ts1, const TSVector<T>& ts2)
338 {
339  TVector<T> tv;
340  if (ts1.mut) { tv.setptr(ts1.vecptr()); tv.setsize(ts1.size()); }
341  else if (ts2.mut) { tv.setptr(ts2.vecptr()); tv.setsize(ts2.size()); }
342  else tv.resize (ts1.size());
343  BCHK(ts1.size() != ts2.size(), VecErr, Wrong dim in ediv, ts2.size(), tv);
344  for (REGISTER unsigned long i = 0; i < tv.size(); i++)
345  tv.set (ts1.get(i), i);
346  if (ts1.mut && ts2.mut) ts2.destroy ();
347  return tv;
348 }
349 
350 INST(template <typename T> class TSVector friend TVector<T> ediv (const TSVector<T>&, const Vector<T>&);)
351 template <typename T>
352 TVector<T> ediv (const TSVector<T>& ts, const Vector<T>& v)
353 {
354  TVector<T> tv;
355  if (ts.mut) { tv.setptr(ts.vecptr()); tv.setsize(ts.size()); }
356  else tv.resize (ts.size());
357  BCHK(v.size() != ts.size(), VecErr, Wrong dim in ediv, v.size(), tv);
358  for (REGISTER unsigned long i = 0; i < tv.size(); i++)
359  tv.set (ts.get(i) / v.get(i), i);
360  return tv;
361 }
362 
363 INST(template <typename T> class TSVector friend TVector<T> ediv (const TSVector<T>&, TVector<T>);)
364 template <typename T>
365 TVector<T> ediv (const TSVector<T>& ts, TVector<T> tv)
366 {
367  BCHK(tv.size() != ts.size(), VecErr, Wrong dim in ediv, tv.size(), tv);
368  for (REGISTER unsigned long i = 0; i < tv.size(); i++)
369  tv.set (ts.get(i) / tv.get(i), i);
370  ts.destroy (); return tv;
371 }
372 
373 
374 INST(template <typename T> class TSVector friend TVector<T> cediv (TVector<T>, const TSVector<T>&);)
375 template <typename T>
376 TVector<T> cediv (TVector<T> tv, const TSVector<T>& ts)
377 {
378  BCHK(tv.size() != ts.size(), VecErr, Wrong dim in ediv, ts.size(), tv);
379  for (REGISTER unsigned long i = 0; i < tv.size(); i++)
380  tv.set (CPLX__ conj(tv.get(i)) / ts.get(i), i);
381  ts.destroy (); return tv;
382 }
383 
384 INST(template <typename T> class TSVector friend TVector<T> cediv (const TSVector<T>&, TVector<T>);)
385 template <typename T>
386 TVector<T> cediv (const TSVector<T>& ts, TVector<T> tv)
387 {
388  BCHK(tv.size() != ts.size(), VecErr, Wrong dim in ediv, tv.size(), tv);
389  for (REGISTER unsigned long i = 0; i < tv.size(); i++)
390  tv.set (CPLX__ conj (ts.get(i)) / tv.get(i), i);
391  ts.destroy (); return tv;
392 }
393 
394 INST(template <typename T> class TSVector friend TVector<T> cediv (const Vector<T>&, const TSVector<T>&);)
395 template <typename T>
396 TVector<T> cediv (const Vector<T>& v, const TSVector<T>& ts)
397 {
398  TVector<T> tv;
399  if (ts.mut) { tv.setptr(ts.vecptr()); tv.setsize(ts.size()); }
400  else tv.resize (ts.size());
401  BCHK(v.size() != ts.size(), VecErr, Wrong dim in ediv, ts.size(), tv);
402  for (REGISTER unsigned long i = 0; i < tv.size(); i++)
403  tv.set (CPLX__ conj(v.get(i)) / ts.get(i), i);
404  return tv;
405 }
406 
407 
408 INST(template <typename T> class TSVector friend TVector<T> cediv (const TSVector<T>&, const Vector<T>&);)
409 template <typename T>
410 TVector<T> cediv (const TSVector<T>& ts, const Vector<T>& v)
411 {
412  TVector<T> tv;
413  if (ts.mut) { tv.setptr(ts.vecptr()); tv.setsize(ts.size()); }
414  else tv.resize (ts.size());
415  BCHK(v.size() != ts.size(), VecErr, Wrong dim in ediv, v.size(), tv);
416  for (REGISTER unsigned long i = 0; i < tv.size(); i++)
417  tv.set (CPLX__ conj (ts.get(i)) / v.get(i), i);
418  return tv;
419 }
420 
421 INST(template <typename T> class TSVector friend TVector<T> cediv (const TSVector<T>&, const TSVector<T>&);)
422 template <typename T>
423 TVector<T> cediv (const TSVector<T>& ts1, const TSVector<T>& ts2)
424 {
425  TVector<T> tv;
426  if (ts1.mut) { tv.setptr(ts1.vecptr()); tv.setsize(ts1.size()); }
427  else if (ts2.mut) { tv.setptr(ts2.vecptr()); tv.setsize(ts2.size()); }
428  else tv.resize (ts1.size());
429  BCHK(ts1.size() != ts2.size(), VecErr, Wrong dim in ediv, ts2.size(), tv);
430  for (REGISTER unsigned long i = 0; i < tv.size(); i++)
431  tv.set (CPLX__ conj (ts1.get(i)), i);
432  if (ts1.mut && ts2.mut) ts2.destroy ();
433  return tv;
434 }
435 
436 
437 
438 INST(template <typename T> class TVector friend TVector<T> ediv (const T&, TVector<T>);)
439 template <typename T>
440 inline TVector<T> ediv (const T& a, TVector<T> b)
441 {
442  STD_SMP_TEMPLATE2C(val_div_vec, b.size(), &b.setval(0), a);
443  return b;
444 }
445 
446 INST(template <typename T> class TVector friend TVector<T> ediv (const T&, const Vector<T>&);)
447 template <typename T>
448 inline TVector<T> ediv (const T& a, const Vector<T>& b)
449 {
450  const unsigned d = b.size();
451  TVector<T> res (d);
452  STD_SMP_TEMPLATE3VC(val_vec_div, d, &res.setval(0), &b.getcref(0), a);
453  return res;
454 }
455 
456 INST(template <typename T> class TSVector friend TVector<T> ediv (const T&, const TSVector<T>&);)
457 template <typename T>
458 inline TVector<T> ediv (const T& a, const TSVector<T>& ts)
459 {
460  TSVector<T> tm (ts); tm.detach ();
461  STD_SMP_TEMPLATE4C(val_svc_div, tm.size(), tm.vecptr(), &ts.getcref(0), a, ts.getfac());
462  tm.getfacref() = (T)1; return TVector<T> (tm);
463 }
464 
465 /* Summing up a vector is seldomly relevant.
466  * TODO: Figure a way to use the algorithms below for other reductions,
467  * such as dot products, fabs, Mat-Vec and Mat-Mat products.
468  * Also create SIMD and parallel versions of it.
469  */
470 
471 /* http://en.wikipedia.org/wiki/Kahan_summation_algorithm */
472 INST(template <typename T> class Vector friend T kahan_sum(const Vector<T>&);)
473 template <typename T> T kahan_sum(const Vector<T>& v)
474 {
475  T sum = T(0.0); T comp = T(0.0);
476  const unsigned long sz = v.size();
477  for (REGISTER unsigned long i = 0; i < sz; ++i) {
478  T y = v(i) - comp;
479  T t = sum + y;
480  comp = (t - sum) - y;
481  sum = t;
482  }
483  return sum /* - comp */;
484 }
485 
486 /* Improved algo */
487 INST(template <typename T> class Vector friend T exact_sum(const Vector<T>&);)
488 template <typename T> T exact_sum(const Vector<T>& v)
489 {
490  T sum = T(0.0); T comp = T(0.0);
491  const unsigned long sz = v.size();
492  for (REGISTER unsigned long i = 0; i < sz; ++i) {
493  T y = v(i);
494  T t = sum + y;
495  comp += (t - sum) - y;
496  sum = t;
497  }
498  return sum - comp;
499 }
500 
501 // elementwise multiplication / division (with and without taking complex conj.)
502 // for matrices
503 
504 # define _VEC getvec()
505 # define _ENDVEC getendvec()
506 # define _DIM size()
507 # define _ROW rows()
508 # define _COL columns()
509 # define _FAC getfac()
510 
511 INST(template <typename T> class TMatrix friend TMatrix<T> emul (const Matrix<T>&, const Matrix<T>&);)
512 template <typename T>
513 TMatrix<T> emul (const Matrix<T>& a, const Matrix<T>& b)
514 {
515  const unsigned int rw = a.rows(), cl = a.columns();
516  BCHK(rw != b.rows(), MatErr, rows differ in emul, rw, TMatrix<T>(0,0));
517  BCHK(cl != b.columns(), MatErr, columns differ in emul, cl, TMatrix<T>(0,0));
518  TMatrix<T> res (rw, cl);
519  do_vec_vec_mul<T> (a._DIM, res._VEC, a._VEC, b._VEC);
520  return res;
521 }
522 
523 INST(template <typename T> class TMatrix friend TMatrix<T> cemul (const Matrix<T>&, const Matrix<T>&);)
524 template <typename T>
525 TMatrix<T> cemul (const Matrix<T>& a, const Matrix<T>& b)
526 {
527  const unsigned int rw = a .rows(), cl = a.columns();
528  BCHK(rw != b.rows(), MatErr, rows differ in cemul, rw, TMatrix<T>(0,0));
529  BCHK(cl != b.columns(), MatErr, columns differ in cemul, cl, TMatrix<T>(0,0));
530  TMatrix<T> res (rw, cl);
531  do_vec_vec_cmul<T> (a._DIM, res._VEC, a._VEC, b._VEC);
532  return res;
533 }
534 
535 INST(template <typename T> class TMatrix friend TMatrix<T> ediv (const Matrix<T>&, const Matrix<T>&);)
536 template <typename T>
537 TMatrix<T> ediv (const Matrix<T>& a, const Matrix<T>& b)
538 {
539  const unsigned int rw = a .rows(), cl = a.columns();
540  BCHK(rw != b.rows(), MatErr, rows differ in ediv, rw, TMatrix<T>(0,0));
541  BCHK(cl != b.columns(), MatErr, columns differ in ediv, cl, TMatrix<T>(0,0));
542  TMatrix<T> res (rw, cl);
543  do_vec_vec_div<T> (a._DIM, res._VEC, a._VEC, b._VEC);
544  return res;
545 }
546 
547 INST(template <typename T> class TMatrix friend TMatrix<T> cediv (const Matrix<T>&, const Matrix<T>&);)
548 template <typename T>
549 TMatrix<T> cediv (const Matrix<T>& a, const Matrix<T>& b)
550 {
551  const unsigned int rw = a .rows(), cl = a.columns();
552  BCHK(rw != b.rows(), MatErr, rows differ in cediv, rw, TMatrix<T>(0,0));
553  BCHK(cl != b.columns(), MatErr, columns differ in cediv, cl, TMatrix<T>(0,0));
554  TMatrix<T> res (rw, cl);
555  do_vec_vec_cdiv<T> (a._DIM, res._VEC, a._VEC, b._VEC);
556  return res;
557 }
558 
559 
560 INST(template <typename T> class TMatrix friend TMatrix<T> ediv (const T&, TMatrix<T>);)
561 template <typename T>
562 inline TMatrix<T> ediv (const T& a, TMatrix<T> b)
563 {
564  STD_SMP_TEMPLATE2C(val_div_vec, b._DIM, b._VEC, a);
565  return b;
566 }
567 
568 INST(template <typename T> class TMatrix friend TMatrix<T> ediv (const T&, const Matrix<T>&);)
569 template <typename T>
570 inline TMatrix<T> ediv (const T& a, const Matrix<T>& b)
571 {
572  const unsigned r = b.rows(), c = b.columns();
573  TMatrix<T> res (r, c);
574  STD_SMP_TEMPLATE3VC(val_vec_div, b._DIM, res._VEC, b._VEC, a);
575  return res;
576 }
577 
578 INST(template <typename T> class TSMatrix friend TMatrix<T> ediv (const T&, const TSMatrix<T>&);)
579 template <typename T>
580 inline TMatrix<T> ediv (const T& a, const TSMatrix<T>& ts)
581 {
582  TSMatrix<T> tm (ts); tm.detach ();
583  STD_SMP_TEMPLATE4C(val_svc_div, tm._DIM, tm._VEC, ts._VEC, a, ts._FAC);
584  tm._FAC = (T)1; return TMatrix<T> (tm);
585 }
586 
587 #undef _DIM
588 #undef _VEC
589 #undef _ENDVEC
590 #undef _ROW
591 #undef _COL
592 #undef _FAC
593 
594 NAMESPACE_END
595 
596 #endif /* TBCI_VECTOR_EXTRA_H */
sum
T sum(const FS_Vector< dims, T > &fv)
Definition: fs_vector.h:599
CPLX__
TSMatrix::detach
void detach(TMatrix< T > *=0)
Definition: matrix.h:1239
TVector::setval
T & setval(const unsigned long i) const
Definition: vector.h:192
BVector::destroy
void destroy()
Definition: bvector.h:268
cemul
TVector< T > cemul(const Vector< T > &a, const Vector< T > &b)
Definition: vector_extra.h:66
BVector::setsize
BVector< T > & setsize(const unsigned long size)
Definition: bvector.h:162
REGISTER
#define REGISTER
Definition: basics.h:108
c
return c
Definition: f_matrix.h:760
NAMESPACE_TBCI
#define NAMESPACE_TBCI
Definition: basics.h:317
TSVector::vecptr
T *const & vecptr() const
Definition: vector.h:1163
TVector::get
tbci_traits< T >::const_refval_type get(const unsigned long i) const
Definition: vector.h:190
res
const unsigned TMatrix< T > * res
Definition: matrix_kernels.h:25
BCHK
#define BCHK(cond, exc, txt, ind, rtval)
Definition: basics.h:575
TSVector::getfac
const T & getfac() const
Definition: vector.h:1073
emul
NAMESPACE_TBCI TVector< T > emul(const Vector< T > &a, const Vector< T > &b)
Definition: vector_extra.h:25
TSVector::getcref
const T & getcref(const unsigned long i) const HOT
Definition: vector.h:1075
TVector::set
T & set(const T &val, const unsigned long i) const
Definition: vector.h:194
STD_SMP_TEMPLATE4C
#define STD_SMP_TEMPLATE4C(oper, dm, a1, a2, a3, a4)
Definition: vector.h:656
TSVector::destroy
void destroy() const
Definition: vector.h:1172
TSVector::get
T get(const unsigned long i) const HOT
Definition: vector.h:1072
cediv
TVector< T > cediv(const Vector< T > &a, const Vector< T > &b)
Definition: vector_extra.h:147
TMatrix::columns
unsigned int columns() const
number of columns
Definition: matrix.h:315
kahan_sum
T kahan_sum(const Vector< T > &v)
Definition: vector_extra.h:473
ts
F_TSMatrix< T > ts
Definition: f_matrix.h:1052
b
F_TMatrix< T > b
Definition: f_matrix.h:736
BVector::setptr
BVector< T > & setptr(T *pointer)
Definition: bvector.h:161
MatErr
exception class
Definition: matrix.h:53
std::conj
int conj(const int arg)
conj for elementary types
Definition: basics.h:1055
TSVector::detach
void detach(const T *=0) const
Definition: vector.h:1179
TVector::getcref
const T & getcref(const unsigned long i) const
Definition: vector.h:188
TVector::size
unsigned long size() const
Definition: vector.h:104
STD_SMP_TEMPLATE2C
#define STD_SMP_TEMPLATE2C(oper, dm, a1, a2)
Definition: vector.h:646
TSVector::mut
bool mut
Definition: vector.h:1037
Matrix
Definition: bvector.h:49
exact_sum
T exact_sum(const Vector< T > &v)
Definition: vector_extra.h:488
INST
#define INST(x)
Definition: basics.h:238
TSVector::getfacref
T & getfacref() const
Definition: vector.h:1074
i
int i
Definition: LM_fit.h:71
BVector::resize
BVector< T > & resize(const BVector< T > &)
Actually it&#39;s a resize and copy (some people would expect the assignment op to do this) ...
Definition: bvector.h:361
y
doublereal y
Definition: TOMS_707.C:27
TVector
Temporary Base Class Idiom: Class TVector is used for temporary variables.
Definition: bvector.h:52
TSVector::size
unsigned long size() const
Definition: vector.h:1120
TSVector::dim
unsigned long dim
Definition: vector.h:1029
STD_SMP_TEMPLATE3VC
#define STD_SMP_TEMPLATE3VC(oper, dm, a1, a2, a3)
Definition: vector.h:650
NAMESPACE_END
#define NAMESPACE_END
Definition: basics.h:323
Vector
Definition: bvector.h:54
T
#define T
Definition: bdmatlib.cc:20
a
const unsigned TMatrix< T > const Matrix< T > * a
Definition: matrix_kernels.h:25
ediv
TVector< T > ediv(const Vector< T > &a, const Vector< T > &b)
Definition: vector_extra.h:110
TMatrix::rows
unsigned int rows() const
number of rows
Definition: matrix.h:317
VecErr
exception class
Definition: bvector.h:31
Generated by   doxygen 1.8.11