1 /***************************************************************************
2 basic_op_new.cpp - GDL operators
3 -------------------
4 begin : July 22 2002
5 copyright : (C) 2002 by Marc Schellens
6 email : m_schellens@users.sf.net
7 ***************************************************************************/
8
9 /***************************************************************************
10 * *
11 * This program is free software; you can redistribute it and/or modify *
12 * it under the terms of the GNU General Public License as published by *
13 * the Free Software Foundation; either version 2 of the License, or *
14 * (at your option) any later version. *
15 * *
16 ***************************************************************************/
17 #ifdef HAVE_CONFIG_H
18 #include <config.h>
19 #endif
20
21 #ifdef _OPENMP
22 #include <omp.h>
23 #endif
24
25 //#include "datatypes.hpp" // for friend declaration
26 #include "nullgdl.hpp"
27 #include "dinterpreter.hpp"
28
29 // needed with gcc-3.3.2
30 #include <cassert>
31
32 #include "typetraits.hpp"
33
34 #include "sigfpehandler.hpp"
35
36 using namespace std;
37
38 // binary operators
39
40 // in basic_op.cpp:
41 // 1. operators that always return a new result (basic_op.cpp)
42 // 2. operators which operate on 'this' (basic_op.cpp)
43
44 // here:
45 // 3. same operators as under 2. that always return a new result
46
47 // AndOp &&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
48 // Ands right and itself into a new DataT_
49 // right must always have more or same number of elements
50 // for integers
51 template<class Sp>
AndOpNew(BaseGDL * r)52 Data_<Sp>* Data_<Sp>::AndOpNew( BaseGDL* r)
53 {
54 Data_* right=static_cast<Data_*>(r);
55
56 ULong nEl=N_Elements();
57 assert( nEl);
58 assert(right->N_Elements());
59
60 Data_* res = NewResult();
61 if( nEl == 1)
62 {
63 (*res)[0] = (*this)[0] & (*right)[0];
64 return res;
65 }
66 TRACEOMP( __FILE__, __LINE__)
67 #pragma omp parallel if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl))
68 {
69 #pragma omp for
70 for( OMPInt i=0; i < nEl; ++i)
71 (*res)[i] = (*this)[i] & (*right)[i]; // & Ty(1);
72 }
73 return res;
74 }
75 // different for floats
76 template<class Sp>
AndOpInvNew(BaseGDL * right)77 Data_<Sp>* Data_<Sp>::AndOpInvNew( BaseGDL* right)
78 {
79 return AndOpNew( right);
80 }
81 // for floats
82 template<>
AndOpNew(BaseGDL * r)83 Data_<SpDFloat>* Data_<SpDFloat>::AndOpNew( BaseGDL* r)
84 {
85 Data_* right=static_cast<Data_*>(r);
86
87 // ULong rEl=right->N_Elements();
88 ULong nEl=N_Elements();
89 Data_* res = NewResult();
90
91 // assert( rEl);
92 assert( nEl);
93 if( nEl == 1)
94 {
95 if ( (*right)[0] == zero ) (*res)[0] = zero; else (*res)[0] = (*this)[0];
96 return res;
97 }
98 TRACEOMP( __FILE__, __LINE__)
99 #pragma omp parallel if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl))
100 {
101 #pragma omp for
102 for ( OMPInt i=0; i < nEl; ++i )
103 if ( (*right)[i] == zero ) (*res)[i] = zero; else (*res)[i] = (*this)[i];
104 }
105 return res;
106 }
107 template<>
AndOpInvNew(BaseGDL * r)108 Data_<SpDFloat>* Data_<SpDFloat>::AndOpInvNew( BaseGDL* r)
109 {
110 Data_* right=static_cast<Data_*>(r);
111
112 ULong nEl=N_Elements();
113 assert( nEl);
114 assert( right->N_Elements());
115
116 Data_* res = NewResult();
117 if( nEl == 1)
118 {
119 if( (*this)[0] != zero) (*res)[0] = (*right)[0]; else (*res)[0] = zero;
120 return res;
121 }
122 TRACEOMP( __FILE__, __LINE__)
123 #pragma omp parallel if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl))
124 {
125 #pragma omp for
126 for( OMPInt i=0; i < nEl; ++i)
127 if( (*this)[i] != zero) (*res)[i] = (*right)[i]; else (*res)[i] = zero;
128 }
129 return res;
130 }
131 // for doubles
132 template<>
AndOpNew(BaseGDL * r)133 Data_<SpDDouble>* Data_<SpDDouble>::AndOpNew( BaseGDL* r)
134 {
135 Data_* right=static_cast<Data_*>(r);
136
137 ULong nEl=N_Elements();
138 assert( nEl);
139 assert( right->N_Elements());
140
141 Data_* res = NewResult();
142 if( nEl == 1)
143 {
144 if ( (*right)[0] == zero ) (*res)[0] = zero; else (*res)[0] = (*this)[0];
145 return res;
146 }
147 TRACEOMP( __FILE__, __LINE__)
148 #pragma omp parallel if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl))
149 {
150 #pragma omp for
151 for( OMPInt i=0; i < nEl; ++i)
152 if ( (*right)[i] == zero ) (*res)[i] = zero; else (*res)[i] = (*this)[i];
153 }
154 return res;
155 }
156 template<>
AndOpInvNew(BaseGDL * r)157 Data_<SpDDouble>* Data_<SpDDouble>::AndOpInvNew( BaseGDL* r)
158 {
159 Data_* right=static_cast<Data_*>(r);
160
161 ULong nEl=N_Elements();
162 assert( nEl);
163 assert( right->N_Elements());
164
165 Data_* res = NewResult();
166 if( nEl == 1)
167 {
168 if( (*this)[0] != zero) (*res)[0] = (*right)[0]; else (*res)[0] = zero;
169 return res;
170 }
171 TRACEOMP( __FILE__, __LINE__)
172 #pragma omp parallel if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl))
173 {
174 #pragma omp for
175 for( OMPInt i=0; i < nEl; ++i)
176 if( (*this)[i] != zero) (*res)[i] = (*right)[i]; else (*res)[i] = zero;
177 }
178 return res;
179 }
180 // invalid types
181 template<>
AndOpNew(BaseGDL * r)182 Data_<SpDString>* Data_<SpDString>::AndOpNew( BaseGDL* r)
183 {
184 throw GDLException("Cannot apply operation to datatype STRING.",true,false);
185 return NULL;
186 }
187 template<>
AndOpNew(BaseGDL * r)188 Data_<SpDComplex>* Data_<SpDComplex>::AndOpNew( BaseGDL* r)
189 {
190 throw GDLException("Cannot apply operation to datatype "+str+".",true,false);
191 return NULL;
192 }
193 template<>
AndOpNew(BaseGDL * r)194 Data_<SpDComplexDbl>* Data_<SpDComplexDbl>::AndOpNew( BaseGDL* r)
195 {
196 throw GDLException("Cannot apply operation to datatype "+str+".",true,false);
197 return NULL;
198 }
199 // template<>
200 // Data_<SpDString>* Data_<SpDString>::AndOpInvNew( BaseGDL* r)
201 // {
202 // throw GDLException("Cannot apply operation to datatype STRING.",true,false);
203 // return res;
204 // }
205 template<>
AndOpNew(BaseGDL * r)206 Data_<SpDPtr>* Data_<SpDPtr>::AndOpNew( BaseGDL* r)
207 {
208 throw GDLException("Cannot apply operation to datatype PTR.",true,false);
209 return NULL;
210 }
211 // template<>
212 // Data_<SpDPtr>* Data_<SpDPtr>::AndOpInvNew( BaseGDL* r)
213 // {
214 // throw GDLException("Cannot apply operation to datatype PTR.",true,false);
215 // return res;
216 // }
217 template<>
AndOpNew(BaseGDL * r)218 Data_<SpDObj>* Data_<SpDObj>::AndOpNew( BaseGDL* r)
219 {
220 throw GDLException("Cannot apply operation to datatype OBJECT.",true,false);
221 return NULL;
222 }
223 // template<>
224 // Data_<SpDPtr>* Data_<SpDPtr>::AndOpInvNew( BaseGDL* r)
225 // {
226 // throw GDLException("Cannot apply operation to datatype PTR.",true,false);
227 // return res;
228 // }
229
230 // scalar versions
231 template<class Sp>
AndOpSNew(BaseGDL * r)232 Data_<Sp>* Data_<Sp>::AndOpSNew( BaseGDL* r)
233 {
234 Data_* right=static_cast<Data_*>(r);
235
236 ULong nEl=N_Elements();
237 assert( nEl);
238
239 Ty s = (*right)[0];
240
241 Data_* res = NewResult();
242 if( nEl == 1)
243 {
244 (*res)[0] = (*this)[0] & s;
245 return res;
246 }
247 TRACEOMP( __FILE__, __LINE__)
248 #pragma omp parallel if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl)) shared(s)
249 {
250 #pragma omp for
251 for( OMPInt i=0; i < nEl; ++i)
252 (*res)[i] = (*this)[i] & s;
253 }
254 return res;
255 }
256 // different for floats
257 template<class Sp>
AndOpInvSNew(BaseGDL * right)258 Data_<Sp>* Data_<Sp>::AndOpInvSNew( BaseGDL* right)
259 {
260 return AndOpSNew( right);
261 }
262 // for floats
263 template<>
AndOpSNew(BaseGDL * r)264 Data_<SpDFloat>* Data_<SpDFloat>::AndOpSNew( BaseGDL* r)
265 {
266 Data_* right=static_cast<Data_*>(r);
267 if( (*right)[0] == zero)
268 {
269 return New( this->dim, BaseGDL::ZERO);
270 }
271 return this->Dup();
272 }
273 template<>
AndOpInvSNew(BaseGDL * r)274 Data_<SpDFloat>* Data_<SpDFloat>::AndOpInvSNew( BaseGDL* r)
275 {
276 Data_* right=static_cast<Data_*>(r);
277
278 ULong nEl=N_Elements();
279 assert( nEl);
280 Ty s = (*right)[0];
281 if( s == zero)
282 {
283 return New( this->dim, BaseGDL::ZERO);
284 }
285 else
286 {
287 Data_* res = NewResult();
288 if( nEl == 1)
289 {
290 if( (*this)[0] != zero) (*res)[0] = s; else (*res)[0] = zero;
291 return res;
292 }
293 TRACEOMP( __FILE__, __LINE__)
294 #pragma omp parallel if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl))
295 {
296 #pragma omp for
297 for( OMPInt i=0; i < nEl; ++i)
298 if( (*this)[i] != zero) (*res)[i] = s; else (*res)[i] = zero;
299 }
300 return res;
301 }
302 }
303 // for doubles
304 template<>
AndOpSNew(BaseGDL * r)305 Data_<SpDDouble>* Data_<SpDDouble>::AndOpSNew( BaseGDL* r)
306 {
307 Data_* right=static_cast<Data_*>(r);
308 if( (*right)[0] == zero)
309 {
310 return New( this->dim, BaseGDL::ZERO);
311 }
312 return this->Dup();
313 }
314 template<>
AndOpInvSNew(BaseGDL * r)315 Data_<SpDDouble>* Data_<SpDDouble>::AndOpInvSNew( BaseGDL* r)
316 {
317 Data_* right=static_cast<Data_*>(r);
318
319 ULong nEl=N_Elements();
320 assert( nEl);
321 Ty s = (*right)[0];
322 if( s == zero)
323 {
324 return New( this->dim, BaseGDL::ZERO);
325 }
326 else
327 {
328 Data_* res = NewResult();
329 if( nEl == 1)
330 {
331 if( (*this)[0] != zero) (*res)[0] = s; else (*res)[0] = zero;
332 return res;
333 }
334 TRACEOMP( __FILE__, __LINE__)
335 #pragma omp parallel if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl))
336 {
337 #pragma omp for
338 for( OMPInt i=0; i < nEl; ++i)
339 if( (*this)[i] != zero) (*res)[i] = s; else (*res)[i] = zero;
340 }
341 return res;
342 }
343 }
344 // invalid types
345 template<>
AndOpSNew(BaseGDL * r)346 Data_<SpDString>* Data_<SpDString>::AndOpSNew( BaseGDL* r)
347 {
348 throw GDLException("Cannot apply operation to datatype STRING.",true,false);
349 return NULL;
350 }
351 template<>
AndOpSNew(BaseGDL * r)352 Data_<SpDComplex>* Data_<SpDComplex>::AndOpSNew( BaseGDL* r)
353 {
354 throw GDLException("Cannot apply operation to datatype "+str+".",true,false);
355 return NULL;
356 }
357 template<>
AndOpSNew(BaseGDL * r)358 Data_<SpDComplexDbl>* Data_<SpDComplexDbl>::AndOpSNew( BaseGDL* r)
359 {
360 throw GDLException("Cannot apply operation to datatype "+str+".",true,false);
361 return NULL;
362 }
363 // template<>
364 // Data_<SpDString>* Data_<SpDString>::AndOpInvSNew( BaseGDL* r)
365 // {
366 // throw GDLException("Cannot apply operation to datatype STRING.",true,false);
367 // return res;
368 // }
369 template<>
AndOpSNew(BaseGDL * r)370 Data_<SpDPtr>* Data_<SpDPtr>::AndOpSNew( BaseGDL* r)
371 {
372 throw GDLException("Cannot apply operation to datatype PTR.",true,false);
373 return NULL;
374 }
375 template<>
AndOpSNew(BaseGDL * r)376 Data_<SpDObj>* Data_<SpDObj>::AndOpSNew( BaseGDL* r)
377 {
378 throw GDLException("Cannot apply operation to datatype OBJECT.",true,false);
379 return NULL;
380 }
381
382
383
384 // OrOp ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
385 // Ors right to itself returns new result
386 // right must always have more or same number of elements
387 // for integers
388 template<class Sp>
OrOpNew(BaseGDL * r)389 Data_<Sp>* Data_<Sp>::OrOpNew( BaseGDL* r)
390 {
391 Data_* right=static_cast<Data_*>(r);
392
393 // ULong rEl=right->N_Elements();
394 ULong nEl=N_Elements();
395 Data_* res = NewResult();
396 // assert( rEl);
397 assert( nEl);
398 //if( !rEl || !nEl) throw GDLException("Variable is undefined.");
399 if( nEl == 1)
400 {
401 (*res)[0] = (*this)[0] | (*right)[0]; // | Ty(1);
402 return res;
403 }
404 TRACEOMP( __FILE__, __LINE__)
405 #pragma omp parallel if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl))
406 {
407 #pragma omp for
408 for( OMPInt i=0; i < nEl; ++i)
409 (*res)[i] = (*this)[i] | (*right)[i]; // | Ty(1);
410 }
411 //C delete right;
412 return res;
413 }
414 // different for floats
415 template<class Sp>
OrOpInvNew(BaseGDL * right)416 Data_<Sp>* Data_<Sp>::OrOpInvNew( BaseGDL* right)
417 {
418 return OrOpNew( right);
419 }
420 // for floats
421 template<>
OrOpNew(BaseGDL * r)422 Data_<SpDFloat>* Data_<SpDFloat>::OrOpNew( BaseGDL* r)
423 {
424 Data_* right=static_cast<Data_*>(r);
425
426 // ULong rEl=right->N_Elements();
427 ULong nEl=N_Elements();
428 Data_* res = NewResult();
429 // assert( rEl);
430 assert( nEl);
431 // if( !rEl || !nEl) throw GDLException("Variable is undefined.");
432 if( nEl == 1)
433 {
434 if( (*this)[0] == zero) (*res)[0] = (*right)[0]; else (*res)[0] = (*this)[0];
435 return res;
436 }
437 TRACEOMP( __FILE__, __LINE__)
438 #pragma omp parallel if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl))
439 {
440 #pragma omp for
441 for( OMPInt i=0; i < nEl; ++i)
442 if( (*this)[i] == zero) (*res)[i] = (*right)[i]; else (*res)[i] = (*this)[i];
443 } //C delete right;
444 return res;
445 }
446 template<>
OrOpInvNew(BaseGDL * r)447 Data_<SpDFloat>* Data_<SpDFloat>::OrOpInvNew( BaseGDL* r)
448 {
449 Data_* right=static_cast<Data_*>(r);
450
451 // ULong rEl=right->N_Elements();
452 ULong nEl=N_Elements(); Data_* res = NewResult();
453 // assert( rEl);
454 assert( nEl);
455 // if( !rEl || !nEl) throw GDLException("Variable is undefined.");
456 if( nEl == 1)
457 {
458 if( (*right)[0] != zero) (*res)[0] = (*right)[0]; else (*res)[0] = (*this)[0];
459 return res;
460 }
461 TRACEOMP( __FILE__, __LINE__)
462 #pragma omp parallel if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl))
463 {
464 #pragma omp for
465 for( OMPInt i=0; i < nEl; ++i)
466 if( (*right)[i] != zero) (*res)[i] = (*right)[i]; else (*res)[i] = (*this)[i];
467 } //C delete right;
468 return res;
469 }
470 // for doubles
471 template<>
OrOpNew(BaseGDL * r)472 Data_<SpDDouble>* Data_<SpDDouble>::OrOpNew( BaseGDL* r)
473 {
474 Data_* right=static_cast<Data_*>(r);
475
476 // ULong rEl=right->N_Elements();
477 ULong nEl=N_Elements();
478 Data_* res = NewResult();
479 // assert( rEl);
480 assert( nEl);
481 // if( !rEl || !nEl) throw GDLException("Variable is undefined.");
482 if( nEl == 1)
483 {
484 if( (*this)[0] == zero) (*res)[0] = (*right)[0]; else (*res)[0] = (*this)[0];
485 return res;
486 }
487 TRACEOMP( __FILE__, __LINE__)
488 #pragma omp parallel if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl))
489 {
490 #pragma omp for
491 for( OMPInt i=0; i < nEl; ++i)
492 if( (*this)[i] == zero) (*res)[i] = (*right)[i]; else (*res)[i] = (*this)[i];
493 } //C delete right;
494 return res;
495 }
496 template<>
OrOpInvNew(BaseGDL * r)497 Data_<SpDDouble>* Data_<SpDDouble>::OrOpInvNew( BaseGDL* r)
498 {
499 Data_* right=static_cast<Data_*>(r);
500
501 // ULong rEl=right->N_Elements();
502 ULong nEl=N_Elements(); Data_* res = NewResult();
503 // assert( rEl);
504 assert( nEl);
505 // if( !rEl || !nEl) throw GDLException("Variable is undefined.");
506 if( nEl == 1)
507 {
508 if( (*right)[0] != zero) (*res)[0] = (*right)[0]; else (*res)[0] = (*this)[0];
509 return res;
510 }
511 TRACEOMP( __FILE__, __LINE__)
512 #pragma omp parallel if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl))
513 {
514 #pragma omp for
515 for( OMPInt i=0; i < nEl; ++i)
516 if( (*right)[i] != zero) (*res)[i] = (*right)[i]; else (*res)[i] = (*this)[i];
517 } //C delete right;
518 return res;
519 }
520 // invalid types
521 template<>
OrOpNew(BaseGDL * r)522 Data_<SpDString>* Data_<SpDString>::OrOpNew( BaseGDL* r)
523 {
524 throw GDLException("Cannot apply operation to datatype STRING.",true,false);
525 return NULL;
526 }
527 template<>
OrOpNew(BaseGDL * r)528 Data_<SpDComplex>* Data_<SpDComplex>::OrOpNew( BaseGDL* r)
529 {
530 throw GDLException("Cannot apply operation to datatype "+str+".",true,false);
531 return NULL;
532 }
533 template<>
OrOpNew(BaseGDL * r)534 Data_<SpDComplexDbl>* Data_<SpDComplexDbl>::OrOpNew( BaseGDL* r)
535 {
536 throw GDLException("Cannot apply operation to datatype "+str+".",true,false);
537 return NULL;
538 }
539 template<>
OrOpNew(BaseGDL * r)540 Data_<SpDPtr>* Data_<SpDPtr>::OrOpNew( BaseGDL* r)
541 {
542 throw GDLException("Cannot apply operation to datatype PTR.",true,false);
543 return NULL;
544 }
545 template<>
OrOpNew(BaseGDL * r)546 Data_<SpDObj>* Data_<SpDObj>::OrOpNew( BaseGDL* r)
547 {
548 throw GDLException("Cannot apply operation to datatype OBJECT.",true,false);
549 return NULL;
550 }
551
552 // OrOpS
553 // for integers
554 template<class Sp>
OrOpSNew(BaseGDL * r)555 Data_<Sp>* Data_<Sp>::OrOpSNew( BaseGDL* r)
556 {
557 Data_* right=static_cast<Data_*>(r);
558
559 ULong nEl=N_Elements();
560 Data_* res = NewResult();
561 assert( nEl);
562 Ty s = (*right)[0];
563 // right->Scalar(s);
564 //s &= Ty(1);
565 // dd |= s;
566 if( nEl == 1)
567 {
568 (*res)[0] = (*this)[0] | s;
569 return res;
570 }
571 TRACEOMP( __FILE__, __LINE__)
572 #pragma omp parallel if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl))
573 {
574 #pragma omp for
575 for( OMPInt i=0; i < nEl; ++i)
576 (*res)[i] = (*this)[i] | s;
577 } //C delete right;
578 return res;
579 }
580 // different for floats
581 template<class Sp>
OrOpInvSNew(BaseGDL * right)582 Data_<Sp>* Data_<Sp>::OrOpInvSNew( BaseGDL* right)
583 {
584 return OrOpSNew( right);
585 }
586 // for floats
587 template<>
OrOpSNew(BaseGDL * r)588 Data_<SpDFloat>* Data_<SpDFloat>::OrOpSNew( BaseGDL* r)
589 {
590 Data_* right=static_cast<Data_*>(r);
591
592 ULong nEl=N_Elements();
593 assert( nEl);
594
595 Ty s = (*right)[0];
596 if( s != zero)
597 {
598 Data_* res = NewResult();
599 if( nEl == 1)
600 {
601 if( (*this)[0] == zero) (*res)[0] = s; else (*res)[0] = (*this)[0];
602 return res;
603 }
604 TRACEOMP( __FILE__, __LINE__)
605 #pragma omp parallel if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl))
606 {
607 #pragma omp for
608 for( OMPInt i=0; i < nEl; ++i)
609 if( (*this)[i] == zero) (*res)[i] = s; else (*res)[i] = (*this)[i];
610 }
611 return res;
612 }
613 else // s == zero
614 {
615 return this->Dup();
616 }
617 }
618 template<>
OrOpInvSNew(BaseGDL * r)619 Data_<SpDFloat>* Data_<SpDFloat>::OrOpInvSNew( BaseGDL* r)
620 {
621 Data_* right=static_cast<Data_*>(r);
622
623 ULong nEl=N_Elements();
624 Data_* res = NewResult();
625 assert( nEl);
626 Ty s = (*right)[0];
627 if( s != zero)
628 {
629 for( SizeT i=0; i < nEl; ++i)
630 (*res)[i] = s;
631 return res;
632 }
633 else
634 {
635 if( nEl == 1)
636 {
637 if( (*this)[0] != zero) (*res)[0] = s; else (*res)[0] = zero;
638 return res;
639 }
640 TRACEOMP( __FILE__, __LINE__)
641 #pragma omp parallel if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl))
642 {
643 #pragma omp for
644 for( OMPInt i=0; i < nEl; ++i)
645 if( (*this)[i] != zero) (*res)[i] = s; else (*res)[i] = zero;
646 }
647 return res;
648 }
649 }
650 // for doubles
651 template<>
OrOpSNew(BaseGDL * r)652 Data_<SpDDouble>* Data_<SpDDouble>::OrOpSNew( BaseGDL* r)
653 {
654 Data_* right=static_cast<Data_*>(r);
655
656 ULong nEl=N_Elements(); Data_* res = NewResult();
657 assert( nEl);
658 Ty s = (*right)[0];
659 // right->Scalar(s);
660 if( s != zero)
661 {
662 if( nEl == 1)
663 {
664 if( (*this)[0] == zero) (*res)[0] = s; else (*res)[0] = (*this)[0];
665 return res;
666 }
667 TRACEOMP( __FILE__, __LINE__)
668 #pragma omp parallel if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl))
669 {
670 #pragma omp for
671 for( OMPInt i=0; i < nEl; ++i)
672 if( (*this)[i] == zero) (*res)[i] = s; else (*res)[i] = (*this)[i];
673 }
674 return res;
675 }
676 // s == zero
677 return this->Dup();
678 }
679 template<>
OrOpInvSNew(BaseGDL * r)680 Data_<SpDDouble>* Data_<SpDDouble>::OrOpInvSNew( BaseGDL* r)
681 {
682 Data_* right=static_cast<Data_*>(r);
683
684 ULong nEl=N_Elements();
685 Data_* res = NewResult();
686 assert( nEl);
687 Ty s = (*right)[0];
688 if( s != zero)
689 {
690 for( SizeT i=0; i < nEl; ++i)
691 (*res)[i] = s;
692 return res;
693 }
694 else
695 {
696 if( nEl == 1)
697 {
698 if( (*this)[0] != zero) (*res)[0] = s; else (*res)[0] = zero;
699 return res;
700 }
701 TRACEOMP( __FILE__, __LINE__)
702 #pragma omp parallel if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl))
703 {
704 #pragma omp for
705 for( OMPInt i=0; i < nEl; ++i)
706 if( (*this)[i] != zero) (*res)[i] = s; else (*res)[i] = zero;
707 }
708 return res;
709 }
710 }
711 // invalid types
712 template<>
OrOpSNew(BaseGDL * r)713 Data_<SpDString>* Data_<SpDString>::OrOpSNew( BaseGDL* r)
714 {
715 throw GDLException("Cannot apply operation to datatype STRING.",true,false);
716 return NULL;
717 }
718 template<>
OrOpSNew(BaseGDL * r)719 Data_<SpDComplex>* Data_<SpDComplex>::OrOpSNew( BaseGDL* r)
720 {
721 throw GDLException("Cannot apply operation to datatype "+str+".",true,false);
722 return NULL;
723 }
724 template<>
OrOpSNew(BaseGDL * r)725 Data_<SpDComplexDbl>* Data_<SpDComplexDbl>::OrOpSNew( BaseGDL* r)
726 {
727 throw GDLException("Cannot apply operation to datatype "+str+".",true,false);
728 return NULL;
729 }
730 template<>
OrOpSNew(BaseGDL * r)731 Data_<SpDPtr>* Data_<SpDPtr>::OrOpSNew( BaseGDL* r)
732 {
733 throw GDLException("Cannot apply operation to datatype PTR.",true,false);
734 return NULL;
735 }
736 template<>
OrOpSNew(BaseGDL * r)737 Data_<SpDObj>* Data_<SpDObj>::OrOpSNew( BaseGDL* r)
738 {
739 throw GDLException("Cannot apply operation to datatype OBJECT.",true,false);
740 return NULL;
741 }
742
743
744
745 // XorOp ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
746 // Xors right to itself, //C deletes right
747 // right must always have more or same number of elements
748 // for integers
749 template<class Sp>
XorOpNew(BaseGDL * r)750 Data_<Sp>* Data_<Sp>::XorOpNew( BaseGDL* r)
751 {
752 Data_* right=static_cast<Data_*>(r);
753
754 ULong nEl=N_Elements();
755 assert( nEl);
756
757 if( nEl == 1)
758 {
759 Data_* res = NewResult();
760 (*res)[0] = (*this)[0] ^ (*right)[0];
761 return res;
762 }
763
764 Ty s;
765 if( right->StrictScalar(s))
766 {
767 if( s == Sp::zero)
768 return this->Dup();
769
770 Data_* res = NewResult();
771 TRACEOMP( __FILE__, __LINE__)
772 #pragma omp parallel if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl))
773 {
774 #pragma omp for
775 for( OMPInt i=0; i < nEl; ++i)
776 (*res)[i] = (*this)[i] ^ s;
777 }
778 return res;
779 }
780 else
781 {
782 Data_* res = NewResult();
783 TRACEOMP( __FILE__, __LINE__)
784 #pragma omp parallel if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl))
785 {
786 #pragma omp for
787 for( OMPInt i=0; i < nEl; ++i)
788 (*res)[i] = (*this)[i] ^ (*right)[i];
789 }
790 return res;
791 }
792 }
793 // invalid types
794 template<>
XorOpNew(BaseGDL * r)795 Data_<SpDFloat>* Data_<SpDFloat>::XorOpNew( BaseGDL* r)
796 {
797 throw GDLException("Cannot apply operation to datatype FLOAT.",true,false);
798 return NULL;
799 }
800 template<>
XorOpNew(BaseGDL * r)801 Data_<SpDDouble>* Data_<SpDDouble>::XorOpNew( BaseGDL* r)
802 {
803 throw GDLException("Cannot apply operation to datatype DOUBLE.",true,false);
804 return NULL;
805 }
806 template<>
XorOpNew(BaseGDL * r)807 Data_<SpDString>* Data_<SpDString>::XorOpNew( BaseGDL* r)
808 {
809 throw GDLException("Cannot apply operation to datatype STRING.",true,false);
810 return NULL;
811 }
812 template<>
XorOpNew(BaseGDL * r)813 Data_<SpDComplex>* Data_<SpDComplex>::XorOpNew( BaseGDL* r)
814 {
815 throw GDLException("Cannot apply operation to datatype "+str+".",true,false);
816 return NULL;
817 }
818 template<>
XorOpNew(BaseGDL * r)819 Data_<SpDComplexDbl>* Data_<SpDComplexDbl>::XorOpNew( BaseGDL* r)
820 {
821 throw GDLException("Cannot apply operation to datatype "+str+".",true,false);
822 return NULL;
823 }
824 template<>
XorOpNew(BaseGDL * r)825 Data_<SpDPtr>* Data_<SpDPtr>::XorOpNew( BaseGDL* r)
826 {
827 throw GDLException("Cannot apply operation to datatype PTR.",true,false);
828 return NULL;
829 }
830 template<>
XorOpNew(BaseGDL * r)831 Data_<SpDObj>* Data_<SpDObj>::XorOpNew( BaseGDL* r)
832 {
833 throw GDLException("Cannot apply operation to datatype OBJECT.",true,false);
834 return NULL;
835 }
836 template<class Sp>
XorOpSNew(BaseGDL * r)837 Data_<Sp>* Data_<Sp>::XorOpSNew( BaseGDL* r)
838 {
839 Data_* right=static_cast<Data_*>(r);
840
841 ULong nEl=N_Elements();
842 assert( nEl);
843 if( nEl == 1)
844 {
845 Data_* res = NewResult();
846 (*res)[0] = (*this)[0] ^ (*right)[0];
847 return res;
848 }
849 Ty s = (*right)[0];
850 if( s == this->zero)
851 return this->Dup();
852
853 Data_* res = NewResult();
854 TRACEOMP( __FILE__, __LINE__)
855 #pragma omp parallel if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl))
856 {
857 #pragma omp for
858 for( OMPInt i=0; i < nEl; ++i)
859 (*res)[i] = (*this)[i] ^ s;
860 }
861 return res;
862 }
863 // different for floats
864 // for floats
865 template<>
XorOpSNew(BaseGDL * r)866 Data_<SpDFloat>* Data_<SpDFloat>::XorOpSNew( BaseGDL* r)
867 {
868 throw GDLException("Cannot apply operation to datatype FLOAT.",true,false);
869 return NULL;
870 }
871 // for doubles
872 template<>
XorOpSNew(BaseGDL * r)873 Data_<SpDDouble>* Data_<SpDDouble>::XorOpSNew( BaseGDL* r)
874 {
875 throw GDLException("Cannot apply operation to datatype DOUBLE.",true,false);
876 return NULL;
877 }
878 // invalid types
879 template<>
XorOpSNew(BaseGDL * r)880 Data_<SpDString>* Data_<SpDString>::XorOpSNew( BaseGDL* r)
881 {
882 throw GDLException("Cannot apply operation to datatype STRING.",true,false);
883 return NULL;
884 }
885 template<>
XorOpSNew(BaseGDL * r)886 Data_<SpDComplex>* Data_<SpDComplex>::XorOpSNew( BaseGDL* r)
887 {
888 throw GDLException("Cannot apply operation to datatype "+str+".",true,false);
889 return NULL;
890 }
891 template<>
XorOpSNew(BaseGDL * r)892 Data_<SpDComplexDbl>* Data_<SpDComplexDbl>::XorOpSNew( BaseGDL* r)
893 {
894 throw GDLException("Cannot apply operation to datatype "+str+".",true,false);
895 return NULL;
896 }
897 template<>
XorOpSNew(BaseGDL * r)898 Data_<SpDPtr>* Data_<SpDPtr>::XorOpSNew( BaseGDL* r)
899 {
900 throw GDLException("Cannot apply operation to datatype PTR.",true,false);
901 return NULL;
902 }
903 template<>
XorOpSNew(BaseGDL * r)904 Data_<SpDObj>* Data_<SpDObj>::XorOpSNew( BaseGDL* r)
905 {
906 throw GDLException("Cannot apply operation to datatype OBJECT.",true,false);
907 return NULL;
908 }
909
910 // Add ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
911 // Adds right to itself returns new result
912 // right must always have more or same number of elements
913 template<class Sp>
AddNew(BaseGDL * r)914 BaseGDL* Data_<Sp>::AddNew( BaseGDL* r)
915 {
916 Data_* right=static_cast<Data_*>(r);
917
918 // ULong rEl=right->N_Elements();
919 ULong nEl=N_Elements();
920 assert( nEl);
921 assert( right->N_Elements());
922
923 Data_* res = NewResult();
924 if( nEl == 1)
925 {
926 (*res)[0] = (*this)[0] + (*right)[0];
927 return res;
928 }
929
930 #ifdef USE_EIGEN
931
932 Eigen::Map<Eigen::Array<Ty,Eigen::Dynamic,1> ,Eigen::Aligned> mThis(&(*this)[0], nEl);
933 Eigen::Map<Eigen::Array<Ty,Eigen::Dynamic,1> ,Eigen::Aligned> mRight(&(*right)[0], nEl);
934 Eigen::Map<Eigen::Array<Ty,Eigen::Dynamic,1> ,Eigen::Aligned> mRes(&(*res)[0], nEl);
935 mRes = mThis + mRight;
936 return res;
937 #else
938
939 TRACEOMP( __FILE__, __LINE__)
940 #pragma omp parallel if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl))
941 {
942 #pragma omp for
943 for( OMPInt i=0; i < nEl; ++i)
944 (*res)[i] = (*this)[i] + (*right)[i];
945 } //C delete right;
946 return res;
947 #endif
948 }
949
950 template<class Sp>
AddInvNew(BaseGDL * r)951 BaseGDL* Data_<Sp>::AddInvNew( BaseGDL* r)
952 {
953 return AddNew( r);
954 }
955
956 template<>
AddInvNew(BaseGDL * r)957 BaseGDL* Data_<SpDString>::AddInvNew( BaseGDL* r)
958 {
959 Data_* right=static_cast<Data_*>(r);
960
961 ULong nEl=N_Elements();
962 assert( nEl);
963
964 Data_* res = NewResult();
965 if( nEl == 1)
966 {
967 (*res)[0] = (*right)[0] + (*this)[0] ;
968 return res;
969 }
970
971 TRACEOMP( __FILE__, __LINE__)
972 #pragma omp parallel if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl))
973 {
974 #pragma omp for
975 for( OMPInt i=0; i < nEl; ++i)
976 (*res)[i] = (*right)[i] + (*this)[i];
977 } //C delete right;
978 return res;
979 }
980
981 // invalid types
982 template<>
AddNew(BaseGDL * r)983 BaseGDL* Data_<SpDPtr>::AddNew( BaseGDL* r)
984 {
985 throw GDLException("Cannot apply operation to datatype PTR.",true,false);
986 return NULL;
987 }
988 template<>
AddNew(BaseGDL * r)989 BaseGDL* Data_<SpDObj>::AddNew( BaseGDL* r)
990 {
991 return Add( r);
992 }
993 template<>
AddInvNew(BaseGDL * r)994 BaseGDL* Data_<SpDObj>::AddInvNew( BaseGDL* r)
995 {
996 return AddInv( r);
997 }
998
999 // scalar versions
1000 template<class Sp>
AddSNew(BaseGDL * r)1001 BaseGDL* Data_<Sp>::AddSNew( BaseGDL* r)
1002 {
1003 Data_* right=static_cast<Data_*>(r);
1004
1005 ULong nEl=N_Elements();
1006 assert( nEl);
1007
1008 Data_* res = NewResult();
1009 if( nEl == 1)
1010 {
1011 (*res)[0] = (*this)[0] + (*right)[0];
1012 return res;
1013 }
1014 Ty s = (*right)[0];
1015 #ifdef USE_EIGEN
1016
1017 Eigen::Map<Eigen::Array<Ty,Eigen::Dynamic,1> ,Eigen::Aligned> mThis(&(*this)[0], nEl);
1018 Eigen::Map<Eigen::Array<Ty,Eigen::Dynamic,1> ,Eigen::Aligned> mRes(&(*res)[0], nEl);
1019 mRes = mThis + s;
1020 return res;
1021 #else
1022
1023 TRACEOMP( __FILE__, __LINE__)
1024 #pragma omp parallel if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl))
1025 {
1026 #pragma omp for
1027 for( OMPInt i=0; i < nEl; ++i)
1028 (*res)[i] = (*this)[i] + s;
1029 } //C delete right;
1030 return res;
1031 #endif
1032
1033 }
1034 template<class Sp>
AddInvSNew(BaseGDL * r)1035 BaseGDL* Data_<Sp>::AddInvSNew( BaseGDL* r)
1036 {
1037 return AddSNew( r);
1038 }
1039 template<>
AddInvSNew(BaseGDL * r)1040 BaseGDL* Data_<SpDString>::AddInvSNew( BaseGDL* r)
1041 {
1042 Data_* right=static_cast<Data_*>(r);
1043
1044 ULong nEl=N_Elements(); Data_* res = NewResult();
1045 assert( nEl);
1046 if( nEl == 1)
1047 {
1048 (*res)[0] = (*right)[0] + (*this)[0];
1049 return res;
1050 }
1051 Ty s = (*right)[0];
1052 TRACEOMP( __FILE__, __LINE__)
1053 #pragma omp parallel if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl))
1054 {
1055 #pragma omp for
1056 for( OMPInt i=0; i < nEl; ++i)
1057 (*res)[i] = s + (*this)[i];
1058 } //C delete right;
1059 return res;
1060 }
1061
1062 // invalid types
1063 template<>
AddSNew(BaseGDL * r)1064 BaseGDL* Data_<SpDPtr>::AddSNew( BaseGDL* r)
1065 {
1066 throw GDLException("Cannot apply operation to datatype PTR.",true,false);
1067 return NULL;
1068 }
1069 template<>
AddSNew(BaseGDL * r)1070 BaseGDL* Data_<SpDObj>::AddSNew( BaseGDL* r)
1071 {
1072 return Add( r);
1073 }
1074 template<>
AddInvSNew(BaseGDL * r)1075 BaseGDL* Data_<SpDObj>::AddInvSNew( BaseGDL* r)
1076 {
1077 return AddInv( r);
1078 }
1079
1080
1081
1082 // Sub ----------------------------------------------------------------------
1083 // substraction: res=left-right
1084 template<class Sp>
SubNew(BaseGDL * r)1085 BaseGDL* Data_<Sp>::SubNew( BaseGDL* r)
1086 {
1087 Data_* right=static_cast<Data_*>(r);
1088
1089 ULong rEl=right->N_Elements();
1090 ULong nEl=N_Elements();
1091 assert( rEl);
1092 assert( nEl);
1093
1094 Data_* res = NewResult();
1095
1096 if( nEl == 1)// && rEl == 1)
1097 {
1098 (*res)[0] = (*this)[0] - (*right)[0];
1099 return res;
1100 }
1101
1102 Ty s;
1103 if( right->StrictScalar(s))
1104 {
1105 #ifdef USE_EIGEN
1106
1107 Eigen::Map<Eigen::Array<Ty,Eigen::Dynamic,1> ,Eigen::Aligned> mThis(&(*this)[0], nEl);
1108 Eigen::Map<Eigen::Array<Ty,Eigen::Dynamic,1> ,Eigen::Aligned> mRes(&(*res)[0], nEl);
1109 mRes = mThis - s;
1110 return res;
1111 #else
1112 TRACEOMP( __FILE__, __LINE__)
1113 #pragma omp parallel if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl))
1114 {
1115 #pragma omp for
1116 for( OMPInt i=0; i < nEl; ++i)
1117 (*res)[i] = (*this)[i] - s;
1118 }
1119 return res;
1120 #endif
1121
1122 }
1123 else
1124 {
1125 #ifdef USE_EIGEN
1126
1127 Eigen::Map<Eigen::Array<Ty,Eigen::Dynamic,1> ,Eigen::Aligned> mThis(&(*this)[0], nEl);
1128 Eigen::Map<Eigen::Array<Ty,Eigen::Dynamic,1> ,Eigen::Aligned> mRight(&(*right)[0], nEl);
1129 Eigen::Map<Eigen::Array<Ty,Eigen::Dynamic,1> ,Eigen::Aligned> mRes(&(*res)[0], nEl);
1130 mRes = mThis - mRight;
1131 return res;
1132 #else
1133 TRACEOMP( __FILE__, __LINE__)
1134 #pragma omp parallel if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl))
1135 {
1136 #pragma omp for
1137 for( OMPInt i=0; i < nEl; ++i)
1138 (*res)[i] = (*this)[i] - (*right)[i];
1139 }
1140 return res;
1141 #endif
1142 }
1143 }
1144 // inverse substraction: left=right-left
1145 template<class Sp>
SubInvNew(BaseGDL * r)1146 BaseGDL* Data_<Sp>::SubInvNew( BaseGDL* r)
1147 {
1148 Data_* right=static_cast<Data_*>(r);
1149
1150 ULong rEl=right->N_Elements();
1151 ULong nEl=N_Elements();
1152 assert( rEl);
1153 assert( nEl);
1154 Data_* res = NewResult();
1155 if( nEl == 1)
1156 {
1157 (*res)[0] = (*right)[0] - (*this)[0];
1158 return res;
1159 }
1160 #ifdef USE_EIGEN
1161
1162 Eigen::Map<Eigen::Array<Ty,Eigen::Dynamic,1> ,Eigen::Aligned> mThis(&(*this)[0], nEl);
1163 Eigen::Map<Eigen::Array<Ty,Eigen::Dynamic,1> ,Eigen::Aligned> mRight(&(*right)[0], nEl);
1164 Eigen::Map<Eigen::Array<Ty,Eigen::Dynamic,1> ,Eigen::Aligned> mRes(&(*res)[0], nEl);
1165 mRes = mRight - mThis;
1166 return res;
1167 #else
1168 TRACEOMP( __FILE__, __LINE__)
1169 #pragma omp parallel if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl))
1170 {
1171 #pragma omp for
1172 for( OMPInt i=0; i < nEl; ++i)
1173 (*res)[i] = (*right)[i] - (*this)[i];
1174 }
1175 return res;
1176 #endif
1177 }
1178 // invalid types
1179 template<>
SubNew(BaseGDL * r)1180 BaseGDL* Data_<SpDString>::SubNew( BaseGDL* r)
1181 {
1182 throw GDLException("Cannot apply operation to datatype STRING.",true,false);
1183 return NULL;
1184 }
1185 template<>
SubInvNew(BaseGDL * r)1186 BaseGDL* Data_<SpDString>::SubInvNew( BaseGDL* r)
1187 {
1188 throw GDLException("Cannot apply operation to datatype STRING.",true,false);
1189 return NULL;
1190 }
1191 template<>
SubNew(BaseGDL * r)1192 BaseGDL* Data_<SpDPtr>::SubNew( BaseGDL* r)
1193 {
1194 throw GDLException("Cannot apply operation to datatype PTR.",true,false);
1195 return NULL;
1196 }
1197 template<>
SubInvNew(BaseGDL * r)1198 BaseGDL* Data_<SpDPtr>::SubInvNew( BaseGDL* r)
1199 {
1200 throw GDLException("Cannot apply operation to datatype PTR.",true,false);
1201 return NULL;
1202 }
1203 template<>
SubNew(BaseGDL * r)1204 BaseGDL* Data_<SpDObj>::SubNew( BaseGDL* r)
1205 {
1206 return Sub( r);
1207 }
1208 template<>
SubInvNew(BaseGDL * r)1209 BaseGDL* Data_<SpDObj>::SubInvNew( BaseGDL* r)
1210 {
1211 return SubInv( r);
1212 }
1213
1214 // scalar versions
1215 template<class Sp>
SubSNew(BaseGDL * r)1216 BaseGDL* Data_<Sp>::SubSNew( BaseGDL* r)
1217 {
1218 Data_* right=static_cast<Data_*>(r);
1219
1220 ULong nEl=N_Elements();
1221 assert( nEl);
1222
1223 Data_* res = NewResult();
1224 if( nEl == 1)
1225 {
1226 (*res)[0] = (*this)[0] - (*right)[0];
1227 return res;
1228 }
1229
1230 Ty s = (*right)[0];
1231 #ifdef USE_EIGEN
1232
1233 Eigen::Map<Eigen::Array<Ty,Eigen::Dynamic,1> ,Eigen::Aligned> mThis(&(*this)[0], nEl);
1234 Eigen::Map<Eigen::Array<Ty,Eigen::Dynamic,1> ,Eigen::Aligned> mRes(&(*res)[0], nEl);
1235 mRes = mThis - s;
1236 return res;
1237 #else
1238
1239 TRACEOMP( __FILE__, __LINE__)
1240 #pragma omp parallel if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl))
1241 {
1242 #pragma omp for
1243 for( OMPInt i=0; i < nEl; ++i)
1244 (*res)[i] = (*this)[i] - s;
1245 }
1246 return res;
1247 #endif
1248
1249 }
1250 // inverse substraction: left=right-left
1251 template<class Sp>
SubInvSNew(BaseGDL * r)1252 BaseGDL* Data_<Sp>::SubInvSNew( BaseGDL* r)
1253 {
1254 Data_* right=static_cast<Data_*>(r);
1255
1256 ULong nEl=N_Elements();
1257 assert( nEl);
1258
1259 Data_* res = NewResult();
1260 if( nEl == 1)
1261 {
1262 (*res)[0] = (*right)[0] - (*this)[0];
1263 return res;
1264 }
1265
1266 Ty s = (*right)[0];
1267 // right->Scalar(s);
1268 // dd = s - dd;
1269 #ifdef USE_EIGEN
1270
1271 Eigen::Map<Eigen::Array<Ty,Eigen::Dynamic,1> ,Eigen::Aligned> mThis(&(*this)[0], nEl);
1272 Eigen::Map<Eigen::Array<Ty,Eigen::Dynamic,1> ,Eigen::Aligned> mRes(&(*res)[0], nEl);
1273 mRes = s - mThis;
1274 return res;
1275 #else
1276 TRACEOMP( __FILE__, __LINE__)
1277 #pragma omp parallel if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl))
1278 {
1279 #pragma omp for
1280 for( OMPInt i=0; i < nEl; ++i)
1281 (*res)[i] = s - (*this)[i];
1282 }
1283 return res;
1284 #endif
1285
1286 }
1287 // invalid types
1288 template<>
SubSNew(BaseGDL * r)1289 BaseGDL* Data_<SpDString>::SubSNew( BaseGDL* r)
1290 {
1291 throw GDLException("Cannot apply operation to datatype STRING.",true,false);
1292 return NULL;
1293 }
1294 template<>
SubInvSNew(BaseGDL * r)1295 BaseGDL* Data_<SpDString>::SubInvSNew( BaseGDL* r)
1296 {
1297 throw GDLException("Cannot apply operation to datatype STRING.",true,false);
1298 return NULL;
1299 }
1300 template<>
SubSNew(BaseGDL * r)1301 BaseGDL* Data_<SpDPtr>::SubSNew( BaseGDL* r)
1302 {
1303 throw GDLException("Cannot apply operation to datatype PTR.",true,false);
1304 return NULL;
1305 }
1306 template<>
SubInvSNew(BaseGDL * r)1307 BaseGDL* Data_<SpDPtr>::SubInvSNew( BaseGDL* r)
1308 {
1309 throw GDLException("Cannot apply operation to datatype PTR.",true,false);
1310 return NULL;
1311 }
1312 template<>
SubSNew(BaseGDL * r)1313 BaseGDL* Data_<SpDObj>::SubSNew( BaseGDL* r)
1314 {
1315 return Sub( r);
1316 }
1317 template<>
SubInvSNew(BaseGDL * r)1318 BaseGDL* Data_<SpDObj>::SubInvSNew( BaseGDL* r)
1319 {
1320 return SubInv( r);
1321 }
1322
1323 // LtMark <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
1324 // LtMarks right to itself, //C deletes right
1325 // right must always have more or same number of elements
1326 template<class Sp>
LtMarkNew(BaseGDL * r)1327 Data_<Sp>* Data_<Sp>::LtMarkNew( BaseGDL* r)
1328 {
1329 Data_* right=static_cast<Data_*>(r);
1330
1331 // ULong rEl=right->N_Elements();
1332 ULong nEl=N_Elements(); Data_* res = NewResult();
1333 // assert( rEl);
1334 assert( nEl);
1335 // if( !rEl || !nEl) throw GDLException("Variable is undefined.");
1336 if( nEl == 1)
1337 {
1338 if( (*this)[0] > (*right)[0]) (*res)[0] = (*right)[0]; else (*res)[0] = (*this)[0];
1339 return res;
1340 }
1341 TRACEOMP( __FILE__, __LINE__)
1342 #pragma omp parallel if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl))
1343 {
1344 #pragma omp for
1345 for( OMPInt i=0; i < nEl; ++i)
1346 if( (*this)[i] > (*right)[i]) (*res)[i] = (*right)[i]; else (*res)[i] = (*this)[i];
1347 } //C delete right;
1348 return res;
1349 }
1350 // invalid types
1351 template<>
LtMarkNew(BaseGDL * r)1352 Data_<SpDString>* Data_<SpDString>::LtMarkNew( BaseGDL* r)
1353 {
1354 throw GDLException("Cannot apply operation to datatype STRING.",true,false);
1355 return NULL;
1356 }
1357 template<>
LtMarkNew(BaseGDL * r)1358 Data_<SpDComplex>* Data_<SpDComplex>::LtMarkNew( BaseGDL* r)
1359 {
1360 throw GDLException("Cannot apply operation to datatype "+str+".",true,false);
1361 return NULL;
1362 }
1363 template<>
LtMarkNew(BaseGDL * r)1364 Data_<SpDComplexDbl>* Data_<SpDComplexDbl>::LtMarkNew( BaseGDL* r)
1365 {
1366 throw GDLException("Cannot apply operation to datatype "+str+".",true,false);
1367 return NULL;
1368 }
1369 template<>
LtMarkNew(BaseGDL * r)1370 Data_<SpDPtr>* Data_<SpDPtr>::LtMarkNew( BaseGDL* r)
1371 {
1372 throw GDLException("Cannot apply operation to datatype PTR.",true,false);
1373 return NULL;
1374 }
1375 template<>
LtMarkNew(BaseGDL * r)1376 Data_<SpDObj>* Data_<SpDObj>::LtMarkNew( BaseGDL* r)
1377 {
1378 throw GDLException("Cannot apply operation to datatype OBJECT.",true,false);
1379 return NULL;
1380 }
1381
1382 // scalar versions
1383 template<class Sp>
LtMarkSNew(BaseGDL * r)1384 Data_<Sp>* Data_<Sp>::LtMarkSNew( BaseGDL* r)
1385 {
1386 Data_* right=static_cast<Data_*>(r);
1387
1388 ULong nEl=N_Elements();
1389 assert( nEl);
1390 Data_* res = NewResult();
1391 if( nEl == 1)
1392 {
1393 if( (*this)[0] > (*right)[0]) (*res)[0] = (*right)[0]; else (*res)[0] = (*this)[0];
1394 return res;
1395 }
1396 Ty s = (*right)[0];
1397 // right->Scalar(s);
1398 TRACEOMP( __FILE__, __LINE__)
1399 #pragma omp parallel if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl))
1400 {
1401 #pragma omp for
1402 for( OMPInt i=0; i < nEl; ++i)
1403 if( (*this)[i] > s) (*res)[i] = s; else (*res)[i] = (*this)[i];
1404 } //C delete right;
1405 return res;
1406 }
1407 // invalid types
1408 template<>
LtMarkSNew(BaseGDL * r)1409 Data_<SpDString>* Data_<SpDString>::LtMarkSNew( BaseGDL* r)
1410 {
1411 throw GDLException("Cannot apply operation to datatype STRING.",true,false);
1412 return NULL;
1413 }
1414 template<>
LtMarkSNew(BaseGDL * r)1415 Data_<SpDComplex>* Data_<SpDComplex>::LtMarkSNew( BaseGDL* r)
1416 {
1417 throw GDLException("Cannot apply operation to datatype "+str+".",true,false);
1418 return NULL;
1419 }
1420 template<>
LtMarkSNew(BaseGDL * r)1421 Data_<SpDComplexDbl>* Data_<SpDComplexDbl>::LtMarkSNew( BaseGDL* r)
1422 {
1423 throw GDLException("Cannot apply operation to datatype "+str+".",true,false);
1424 return NULL;
1425 }
1426 template<>
LtMarkSNew(BaseGDL * r)1427 Data_<SpDPtr>* Data_<SpDPtr>::LtMarkSNew( BaseGDL* r)
1428 {
1429 throw GDLException("Cannot apply operation to datatype PTR.",true,false);
1430 return NULL;
1431 }
1432 template<>
LtMarkSNew(BaseGDL * r)1433 Data_<SpDObj>* Data_<SpDObj>::LtMarkSNew( BaseGDL* r)
1434 {
1435 throw GDLException("Cannot apply operation to datatype OBJECT.",true,false);
1436 return NULL;
1437 }
1438
1439 // GtMark >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
1440 // GtMarks right to itself returns new result
1441 // right must always have more or same number of elements
1442 template<class Sp>
GtMarkNew(BaseGDL * r)1443 Data_<Sp>* Data_<Sp>::GtMarkNew( BaseGDL* r)
1444 {
1445 Data_* right=static_cast<Data_*>(r);
1446
1447 // ULong rEl=right->N_Elements();
1448 ULong nEl=N_Elements(); Data_* res = NewResult();
1449 // assert( rEl);
1450 assert( nEl);
1451 // if( !rEl || !nEl) throw GDLException("Variable is undefined.");
1452 if( nEl == 1)
1453 {
1454 if( (*this)[0] < (*right)[0]) (*res)[0] = (*right)[0]; else (*res)[0] = (*this)[0];
1455 return res;
1456 }
1457 TRACEOMP( __FILE__, __LINE__)
1458 #pragma omp parallel if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl))
1459 {
1460 #pragma omp for
1461 for( OMPInt i=0; i < nEl; ++i)
1462 if( (*this)[i] < (*right)[i]) (*res)[i] = (*right)[i]; else (*res)[i] = (*this)[i];
1463 } //C delete right;
1464 return res;
1465 }
1466 // invalid types
1467 template<>
GtMarkNew(BaseGDL * r)1468 Data_<SpDString>* Data_<SpDString>::GtMarkNew( BaseGDL* r)
1469 {
1470 throw GDLException("Cannot apply operation to datatype STRING.",true,false);
1471 return NULL;
1472 }
1473 template<>
GtMarkNew(BaseGDL * r)1474 Data_<SpDComplex>* Data_<SpDComplex>::GtMarkNew( BaseGDL* r)
1475 {
1476 throw GDLException("Cannot apply operation to datatype "+str+".",true,false);
1477 return NULL;
1478 }
1479 template<>
GtMarkNew(BaseGDL * r)1480 Data_<SpDComplexDbl>* Data_<SpDComplexDbl>::GtMarkNew( BaseGDL* r)
1481 {
1482 throw GDLException("Cannot apply operation to datatype "+str+".",true,false);
1483 return NULL;
1484 }
1485 template<>
GtMarkNew(BaseGDL * r)1486 Data_<SpDPtr>* Data_<SpDPtr>::GtMarkNew( BaseGDL* r)
1487 {
1488 throw GDLException("Cannot apply operation to datatype PTR.",true,false);
1489 return NULL;
1490 }
1491 template<>
GtMarkNew(BaseGDL * r)1492 Data_<SpDObj>* Data_<SpDObj>::GtMarkNew( BaseGDL* r)
1493 {
1494 throw GDLException("Cannot apply operation to datatype OBJECT.",true,false);
1495 return NULL;
1496 }
1497
1498 // scalar versions
1499 template<class Sp>
GtMarkSNew(BaseGDL * r)1500 Data_<Sp>* Data_<Sp>::GtMarkSNew( BaseGDL* r)
1501 {
1502 Data_* right=static_cast<Data_*>(r);
1503
1504 ULong nEl=N_Elements(); Data_* res = NewResult();
1505 assert( nEl);
1506 if( nEl == 1)
1507 {
1508 if( (*this)[0] < (*right)[0]) (*res)[0] = (*right)[0]; else (*res)[0] = (*this)[0];
1509 return res;
1510 }
1511
1512 Ty s = (*right)[0];
1513 // right->Scalar(s);
1514 TRACEOMP( __FILE__, __LINE__)
1515 #pragma omp parallel if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl))
1516 {
1517 #pragma omp for
1518 for( OMPInt i=0; i < nEl; ++i)
1519 if( (*this)[i] < s) (*res)[i] = s; else (*res)[i] = (*this)[i];
1520 } ;
1521 return res;
1522 }
1523 // invalid types
1524 template<>
GtMarkSNew(BaseGDL * r)1525 Data_<SpDString>* Data_<SpDString>::GtMarkSNew( BaseGDL* r)
1526 {
1527 throw GDLException("Cannot apply operation to datatype STRING.",true,false);
1528 return NULL;
1529 }
1530 template<>
GtMarkSNew(BaseGDL * r)1531 Data_<SpDComplex>* Data_<SpDComplex>::GtMarkSNew( BaseGDL* r)
1532 {
1533 throw GDLException("Cannot apply operation to datatype "+str+".",true,false);
1534 return NULL;
1535 }
1536 template<>
GtMarkSNew(BaseGDL * r)1537 Data_<SpDComplexDbl>* Data_<SpDComplexDbl>::GtMarkSNew( BaseGDL* r)
1538 {
1539 throw GDLException("Cannot apply operation to datatype "+str+".",true,false);
1540 return NULL;
1541 }
1542 template<>
GtMarkSNew(BaseGDL * r)1543 Data_<SpDPtr>* Data_<SpDPtr>::GtMarkSNew( BaseGDL* r)
1544 {
1545 throw GDLException("Cannot apply operation to datatype PTR.",true,false);
1546 return NULL;
1547 }
1548 template<>
GtMarkSNew(BaseGDL * r)1549 Data_<SpDObj>* Data_<SpDObj>::GtMarkSNew( BaseGDL* r)
1550 {
1551 throw GDLException("Cannot apply operation to datatype OBJECT.",true,false);
1552 return NULL;
1553 }
1554
1555 // Mult *********************************************************************
1556 // Mults right to itself, //C deletes right
1557 // right must always have more or same number of elements
1558 template<class Sp>
MultNew(BaseGDL * r)1559 Data_<Sp>* Data_<Sp>::MultNew( BaseGDL* r)
1560 {
1561 Data_* right=static_cast<Data_*>(r);
1562
1563 Data_* res=NewResult();
1564
1565 // ULong rEl=right->N_Elements();
1566 ULong nEl=N_Elements();
1567 // assert( rEl);
1568 assert( nEl);
1569 // if( !rEl || !nEl) throw GDLException("Variable is undefined.");
1570 if( nEl == 1)
1571 {
1572 (*res)[0] = (*this)[0] * (*right)[0];
1573 return res;
1574 }
1575 #ifdef USE_EIGEN
1576
1577 Eigen::Map<Eigen::Array<Ty,Eigen::Dynamic,1> ,Eigen::Aligned> mThis(&(*this)[0], nEl);
1578 Eigen::Map<Eigen::Array<Ty,Eigen::Dynamic,1> ,Eigen::Aligned> mRight(&(*right)[0], nEl);
1579 Eigen::Map<Eigen::Array<Ty,Eigen::Dynamic,1> ,Eigen::Aligned> mRes(&(*res)[0], nEl);
1580 mRes = mThis * mRight;
1581 return res;
1582 #else
1583 TRACEOMP( __FILE__, __LINE__)
1584 #pragma omp parallel if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl))
1585 {
1586 #pragma omp for
1587 for( OMPInt i=0; i < nEl; ++i)
1588 (*res)[i] = (*this)[i] * (*right)[i];
1589 } //C delete right;
1590 return res;
1591 #endif
1592
1593 }
1594 // invalid types
1595 template<>
MultNew(BaseGDL * r)1596 Data_<SpDString>* Data_<SpDString>::MultNew( BaseGDL* r)
1597 {
1598 throw GDLException("Cannot apply operation to datatype STRING.",true,false);
1599 return NULL;
1600 }
1601 template<>
MultNew(BaseGDL * r)1602 Data_<SpDPtr>* Data_<SpDPtr>::MultNew( BaseGDL* r)
1603 {
1604 throw GDLException("Cannot apply operation to datatype PTR.",true,false);
1605 return NULL;
1606 }
1607 template<>
MultNew(BaseGDL * r)1608 Data_<SpDObj>* Data_<SpDObj>::MultNew( BaseGDL* r)
1609 {
1610 throw GDLException("Cannot apply operation to datatype OBJECT.",true,false);
1611 return NULL;
1612 }
1613
1614 // scalar versions
1615 template<class Sp>
MultSNew(BaseGDL * r)1616 Data_<Sp>* Data_<Sp>::MultSNew( BaseGDL* r )
1617 {
1618 Data_* right=static_cast<Data_*> ( r );
1619
1620 ULong nEl=N_Elements();
1621 assert ( nEl );
1622
1623 Data_* res = NewResult();
1624 if ( nEl == 1 )
1625 {
1626 ( *res )[0] = ( *this )[0] * ( *right )[0];
1627 return res;
1628 }
1629 Ty s = ( *right ) [0];
1630 #ifdef USE_EIGEN
1631
1632 Eigen::Map<Eigen::Array<Ty,Eigen::Dynamic,1> ,Eigen::Aligned> mThis(&(*this)[0], nEl);
1633 Eigen::Map<Eigen::Array<Ty,Eigen::Dynamic,1> ,Eigen::Aligned> mRes(&(*res)[0], nEl);
1634 mRes = mThis * s;
1635 return res;
1636 #else
1637 TRACEOMP ( __FILE__, __LINE__ )
1638 #pragma omp parallel if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl))
1639 {
1640 #pragma omp for
1641 for ( OMPInt i=0; i < nEl; ++i )
1642 (*res ) [i] = (*this )[i] * s;
1643 }
1644 return res;
1645 #endif
1646
1647 }
1648 // invalid types
1649 template<>
MultSNew(BaseGDL * r)1650 Data_<SpDString>* Data_<SpDString>::MultSNew( BaseGDL* r)
1651 {
1652 throw GDLException("Cannot apply operation to datatype STRING.",true,false);
1653 return NULL;
1654 }
1655 template<>
MultSNew(BaseGDL * r)1656 Data_<SpDPtr>* Data_<SpDPtr>::MultSNew( BaseGDL* r)
1657 {
1658 throw GDLException("Cannot apply operation to datatype PTR.",true,false);
1659 return NULL;
1660 }
1661 template<>
MultSNew(BaseGDL * r)1662 Data_<SpDObj>* Data_<SpDObj>::MultSNew( BaseGDL* r)
1663 {
1664 throw GDLException("Cannot apply operation to datatype OBJECT.",true,false);
1665 return NULL;
1666 }
1667
1668
1669
1670 // Div //////////////////////////////////////////////////////////////////////
1671 // division: left=left/right
1672 template<class Sp>
DivNew(BaseGDL * r)1673 Data_<Sp>* Data_<Sp>::DivNew( BaseGDL* r)
1674 {
1675 Data_* right=static_cast<Data_*>(r);
1676
1677 // ULong rEl=right->N_Elements();
1678 ULong nEl=N_Elements();
1679 assert( nEl);
1680
1681 Data_* res = NewResult();
1682
1683 SizeT i = 0;
1684
1685 if( sigsetjmp( sigFPEJmpBuf, 1) == 0)
1686 {
1687 // TODO: Check if we can use OpenMP here (is longjmp allowed?)
1688 // if yes: need to run the full loop after the longjmp
1689 for( ; i < nEl; ++i)
1690 (*res)[i] = (*this)[i] / (*right)[i];
1691 return res;
1692 }
1693 else
1694 {
1695 TRACEOMP( __FILE__, __LINE__)
1696 #pragma omp parallel if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl))
1697 {
1698 #pragma omp for
1699 for( OMPInt ix=i; ix < nEl; ++ix)
1700 if( (*right)[ix] != this->zero)
1701 (*res)[ix] = (*this)[ix] / (*right)[ix];
1702 else
1703 (*res)[ix] = (*this)[ix];
1704 }
1705 return res;
1706 }
1707 }
1708 // inverse division: left=right/left
1709 template<class Sp>
DivInvNew(BaseGDL * r)1710 Data_<Sp>* Data_<Sp>::DivInvNew( BaseGDL* r)
1711 {
1712 Data_* right=static_cast<Data_*>(r);
1713
1714 // ULong rEl=right->N_Elements();
1715 ULong nEl=N_Elements(); Data_* res = NewResult();
1716 // assert( rEl);
1717 assert( nEl);
1718
1719 SizeT i = 0;
1720
1721 // if( !rEl || !nEl) throw GDLException("Variable is undefined.");
1722 if( sigsetjmp( sigFPEJmpBuf, 1) == 0)
1723 {
1724 for( /*SizeT i=0*/; i < nEl; ++i)
1725 (*res)[i] = (*right)[i] / (*this)[i];
1726 return res;
1727 }
1728 else
1729 {
1730 TRACEOMP( __FILE__, __LINE__)
1731 #pragma omp parallel if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl))
1732 {
1733 #pragma omp for
1734 for( OMPInt ix=i; ix < nEl; ++ix)
1735 if( (*this)[ix] != this->zero)
1736 (*res)[ix] = (*right)[ix] / (*this)[ix];
1737 else
1738 (*res)[ix] = (*right)[ix];
1739 } //C delete right;
1740 return res;
1741 }
1742 }
1743 // invalid types
1744 template<>
DivNew(BaseGDL * r)1745 Data_<SpDString>* Data_<SpDString>::DivNew( BaseGDL* r)
1746 {
1747 throw GDLException("Cannot apply operation to datatype STRING.",true,false);
1748 return NULL;
1749 }
1750 template<>
DivInvNew(BaseGDL * r)1751 Data_<SpDString>* Data_<SpDString>::DivInvNew( BaseGDL* r)
1752 {
1753 throw GDLException("Cannot apply operation to datatype STRING.",true,false);
1754 return NULL;
1755 }
1756 template<>
DivNew(BaseGDL * r)1757 Data_<SpDPtr>* Data_<SpDPtr>::DivNew( BaseGDL* r)
1758 {
1759 throw GDLException("Cannot apply operation to datatype PTR.",true,false);
1760 return NULL;
1761 }
1762 template<>
DivInvNew(BaseGDL * r)1763 Data_<SpDPtr>* Data_<SpDPtr>::DivInvNew( BaseGDL* r)
1764 {
1765 throw GDLException("Cannot apply operation to datatype PTR.",true,false);
1766 return NULL;
1767 }
1768 template<>
DivNew(BaseGDL * r)1769 Data_<SpDObj>* Data_<SpDObj>::DivNew( BaseGDL* r)
1770 {
1771 throw GDLException("Cannot apply operation to datatype OBJECT.",true,false);
1772 return NULL;
1773 }
1774 template<>
DivInvNew(BaseGDL * r)1775 Data_<SpDObj>* Data_<SpDObj>::DivInvNew( BaseGDL* r)
1776 {
1777 throw GDLException("Cannot apply operation to datatype OBJECT.",true,false);
1778 return NULL;
1779 }
1780
1781 // scalar versions
1782 template<class Sp>
DivSNew(BaseGDL * r)1783 Data_<Sp>* Data_<Sp>::DivSNew( BaseGDL* r)
1784 {
1785 Data_* right=static_cast<Data_*>(r);
1786
1787 ULong nEl=N_Elements();
1788 assert( nEl);
1789 Ty s = (*right)[0];
1790 SizeT i=0;
1791 Data_* res = NewResult();
1792 if( s != this->zero)
1793 {
1794 for( SizeT i=0; i < nEl; ++i)
1795 (*res)[i] = (*this)[i] / s;
1796 return res;
1797 }
1798
1799 if( sigsetjmp( sigFPEJmpBuf, 1) == 0)
1800 {
1801 for( SizeT i=0; i < nEl; ++i)
1802 (*res)[i] = (*this)[i] / s;
1803 }
1804 else
1805 {
1806 for( SizeT i=0; i < nEl; ++i)
1807 (*res)[i] = (*this)[i];
1808 // }
1809 }
1810 return res;
1811 }
1812
1813 // inverse division: left=right/left
1814 template<class Sp>
DivInvSNew(BaseGDL * r)1815 Data_<Sp>* Data_<Sp>::DivInvSNew( BaseGDL* r)
1816 {
1817 Data_* right=static_cast<Data_*>(r);
1818
1819 ULong nEl=N_Elements(); Data_* res = NewResult();
1820 assert( nEl);
1821 // if( !rEl || !nEl) throw GDLException("Variable is undefined.");
1822
1823 if( nEl == 1 && (*this)[0] != this->zero)
1824 {
1825 (*res)[0] = (*right)[0] / (*this)[0];
1826 return res;
1827 }
1828
1829 Ty s = (*right)[0];
1830 SizeT i=0;
1831 if( sigsetjmp( sigFPEJmpBuf, 1) == 0)
1832 {
1833 for( ; i < nEl; ++i)
1834 (*res)[i] = s / (*this)[i];
1835 return res;
1836 }
1837 else
1838 {
1839 TRACEOMP( __FILE__, __LINE__)
1840 #pragma omp parallel if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl))
1841 {
1842 #pragma omp for
1843 for( OMPInt ix=i; ix < nEl; ++ix)
1844 if( (*this)[ix] != this->zero)
1845 (*res)[ix] = s / (*this)[ix];
1846 else
1847 (*res)[ix] = s;
1848 }
1849 return res;
1850 }
1851 }
1852 // invalid types
1853 template<>
DivSNew(BaseGDL * r)1854 Data_<SpDString>* Data_<SpDString>::DivSNew( BaseGDL* r)
1855 {
1856 throw GDLException("Cannot apply operation to datatype STRING.",true,false);
1857 return NULL;
1858 }
1859 template<>
DivInvSNew(BaseGDL * r)1860 Data_<SpDString>* Data_<SpDString>::DivInvSNew( BaseGDL* r)
1861 {
1862 throw GDLException("Cannot apply operation to datatype STRING.",true,false);
1863 return NULL;
1864 }
1865 template<>
DivSNew(BaseGDL * r)1866 Data_<SpDPtr>* Data_<SpDPtr>::DivSNew( BaseGDL* r)
1867 {
1868 throw GDLException("Cannot apply operation to datatype PTR.",true,false);
1869 return NULL;
1870 }
1871 template<>
DivInvSNew(BaseGDL * r)1872 Data_<SpDPtr>* Data_<SpDPtr>::DivInvSNew( BaseGDL* r)
1873 {
1874 throw GDLException("Cannot apply operation to datatype PTR.",true,false);
1875 return NULL;
1876 }
1877 template<>
DivSNew(BaseGDL * r)1878 Data_<SpDObj>* Data_<SpDObj>::DivSNew( BaseGDL* r)
1879 {
1880 throw GDLException("Cannot apply operation to datatype OBJECT.",true,false);
1881 return NULL;
1882 }
1883 template<>
DivInvSNew(BaseGDL * r)1884 Data_<SpDObj>* Data_<SpDObj>::DivInvSNew( BaseGDL* r)
1885 {
1886 throw GDLException("Cannot apply operation to datatype OBJECT.",true,false);
1887 return NULL;
1888 }
1889
1890
1891
1892 // Mod %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1893 // modulo division: left=left % right
1894
1895 // float modulo division: left=left % right
Modulo(const DFloat & l,const DFloat & r)1896 inline DFloat Modulo( const DFloat& l, const DFloat& r)
1897 {
1898 // float t=abs(l/r);
1899 // if( l < 0.0) return t=(floor(t)-t)*abs(r);
1900 // return (t-floor(t))*abs(r);
1901 return fmod(l,r);
1902 }
1903 // in basic_op.cpp
1904 // double modulo division: left=left % right
DModulo(const DDouble & l,const DDouble & r)1905 inline DDouble DModulo( const DDouble& l, const DDouble& r)
1906 {
1907 // DDouble t=abs(l/r);
1908 // if( l < 0.0) return t=(floor(t)-t)*abs(r);
1909 // return (t-floor(t))*abs(r);
1910 return fmod(l,r);
1911 }
1912
1913 template<class Sp>
ModNew(BaseGDL * r)1914 Data_<Sp>* Data_<Sp>::ModNew( BaseGDL* r)
1915 {
1916 Data_* right=static_cast<Data_*>(r);
1917
1918 // ULong rEl=right->N_Elements();
1919 ULong nEl=N_Elements();
1920 Data_* res = NewResult();
1921 // assert( rEl);
1922 assert( nEl);
1923
1924 SizeT i=0;
1925 if( sigsetjmp( sigFPEJmpBuf, 1) == 0)
1926 {
1927 for( ; i < nEl; ++i)
1928 (*res)[i] = (*this)[i] % (*right)[i];
1929 return res;
1930 }
1931 else
1932 {
1933 TRACEOMP( __FILE__, __LINE__)
1934 #pragma omp parallel if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl))
1935 {
1936 #pragma omp for
1937 for( OMPInt ix=i; ix < nEl; ++ix)
1938 if( (*right)[ix] != this->zero)
1939 (*res)[ix] = (*this)[ix] % (*right)[ix];
1940 else
1941 (*res)[ix] = this->zero;
1942 }
1943 return res;
1944 }
1945 }
1946 // inverse modulo division: left=right % left
1947 template<class Sp>
ModInvNew(BaseGDL * r)1948 Data_<Sp>* Data_<Sp>::ModInvNew( BaseGDL* r)
1949 {
1950 Data_* right=static_cast<Data_*>(r);
1951
1952 ULong nEl=N_Elements();
1953 Data_* res = NewResult();
1954 assert( nEl);
1955
1956 SizeT i=0;
1957
1958 if( sigsetjmp( sigFPEJmpBuf, 1) == 0)
1959 {
1960 for( ; i < nEl; ++i)
1961 (*res)[i] = (*right)[i] % (*this)[i];
1962 return res;
1963 }
1964 else
1965 {
1966 TRACEOMP( __FILE__, __LINE__)
1967 #pragma omp parallel if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl))
1968 {
1969 #pragma omp for
1970 for( OMPInt ix=i; ix < nEl; ++ix)
1971 if( (*this)[ix] != this->zero)
1972 (*res)[ix] = (*right)[ix] % (*this)[ix];
1973 else
1974 (*res)[ix] = this->zero;
1975 }
1976 return res;
1977 }
1978 }
1979 // in basic_op.cpp
1980 // float modulo division: left=left % right
1981
1982 template<>
ModNew(BaseGDL * r)1983 Data_<SpDFloat>* Data_<SpDFloat>::ModNew( BaseGDL* r)
1984 {
1985 Data_* right=static_cast<Data_*>(r);
1986
1987 ULong nEl=N_Elements();
1988 Data_* res = NewResult();
1989
1990 assert( nEl);
1991 if( nEl == 1)
1992 {
1993 (*res)[0] = Modulo((*this)[0],(*right)[0]);
1994 return res;
1995 }
1996
1997 TRACEOMP( __FILE__, __LINE__)
1998 #pragma omp parallel if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl))
1999 {
2000 #pragma omp for
2001 for( OMPInt i=0; i < nEl; ++i)
2002 (*res)[i] = Modulo((*this)[i],(*right)[i]);
2003 }
2004 return res;
2005 }
2006 // float inverse modulo division: left=right % left
2007 template<>
ModInvNew(BaseGDL * r)2008 Data_<SpDFloat>* Data_<SpDFloat>::ModInvNew( BaseGDL* r)
2009 {
2010 Data_* right=static_cast<Data_*>(r);
2011
2012 // ULong rEl=right->N_Elements();
2013 ULong nEl=N_Elements();
2014 Data_* res = NewResult();
2015
2016 assert( nEl);
2017 if( nEl == 1)
2018 {
2019 (*res)[0] = Modulo((*right)[0],(*this)[0]);
2020 return res;
2021 }
2022
2023 TRACEOMP( __FILE__, __LINE__)
2024 #pragma omp parallel if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl))
2025 {
2026 #pragma omp for
2027 for( OMPInt i=0; i < nEl; ++i)
2028 (*res)[i] = Modulo((*right)[i],(*this)[i]);
2029 }
2030 return res;
2031 }
2032 // double modulo division: left=left % right
2033
2034 template<>
ModNew(BaseGDL * r)2035 Data_<SpDDouble>* Data_<SpDDouble>::ModNew( BaseGDL* r)
2036 {
2037 Data_* right=static_cast<Data_*>(r);
2038
2039 // ULong rEl=right->N_Elements();
2040 ULong nEl=N_Elements();
2041 assert( nEl);
2042
2043 Data_* res = NewResult();
2044 if( nEl == 1)
2045 {
2046 (*res)[0] = DModulo((*this)[0],(*right)[0]);
2047 return res;
2048 }
2049
2050 TRACEOMP( __FILE__, __LINE__)
2051 #pragma omp parallel if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl))
2052 {
2053 #pragma omp for
2054 for( OMPInt i=0; i < nEl; ++i)
2055 (*res)[i] = DModulo((*this)[i],(*right)[i]);
2056 }
2057 return res;
2058 }
2059 // double inverse modulo division: left=right % left
2060 template<>
ModInvNew(BaseGDL * r)2061 Data_<SpDDouble>* Data_<SpDDouble>::ModInvNew( BaseGDL* r)
2062 {
2063 Data_* right=static_cast<Data_*>(r);
2064
2065 // ULong rEl=right->N_Elements();
2066 ULong nEl=N_Elements();
2067 Data_* res = NewResult();
2068
2069 assert( nEl);
2070 if( nEl == 1)
2071 {
2072 (*res)[0] = DModulo((*right)[0],(*this)[0]);
2073 return res;
2074 }
2075
2076 TRACEOMP( __FILE__, __LINE__)
2077 #pragma omp parallel if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl))
2078 {
2079 #pragma omp for
2080 for( OMPInt i=0; i < nEl; ++i)
2081 (*res)[i] = DModulo((*right)[i],(*this)[i]);
2082 }
2083 return res;
2084 }
2085 // invalid types
2086 template<>
ModNew(BaseGDL * r)2087 Data_<SpDString>* Data_<SpDString>::ModNew( BaseGDL* r)
2088 {
2089 throw GDLException("Cannot apply operation to datatype STRING.",true,false);
2090 return NULL;
2091 }
2092 template<>
ModInvNew(BaseGDL * r)2093 Data_<SpDString>* Data_<SpDString>::ModInvNew( BaseGDL* r)
2094 {
2095 throw GDLException("Cannot apply operation to datatype STRING.",true,false);
2096 return NULL;
2097 }
2098 template<>
ModNew(BaseGDL * r)2099 Data_<SpDComplex>* Data_<SpDComplex>::ModNew( BaseGDL* r)
2100 {
2101 throw GDLException("Cannot apply operation to datatype "+str+".",true,false);
2102 return NULL;
2103 }
2104 template<>
ModNew(BaseGDL * r)2105 Data_<SpDComplexDbl>* Data_<SpDComplexDbl>::ModNew( BaseGDL* r)
2106 {
2107 throw GDLException("Cannot apply operation to datatype "+str+".",true,false);
2108 return NULL;
2109 }
2110 template<>
ModInvNew(BaseGDL * r)2111 Data_<SpDComplex>* Data_<SpDComplex>::ModInvNew( BaseGDL* r)
2112 {
2113 throw GDLException("Cannot apply operation to datatype "+str+".",true,false);
2114 return NULL;
2115 }
2116 template<>
ModInvNew(BaseGDL * r)2117 Data_<SpDComplexDbl>* Data_<SpDComplexDbl>::ModInvNew( BaseGDL* r)
2118 {
2119 throw GDLException("Cannot apply operation to datatype "+str+".",true,false);
2120 return NULL;
2121 }
2122 template<>
ModNew(BaseGDL * r)2123 Data_<SpDPtr>* Data_<SpDPtr>::ModNew( BaseGDL* r)
2124 {
2125 throw GDLException("Cannot apply operation to datatype PTR.",true,false);
2126 return NULL;
2127 }
2128 template<>
ModInvNew(BaseGDL * r)2129 Data_<SpDPtr>* Data_<SpDPtr>::ModInvNew( BaseGDL* r)
2130 {
2131 throw GDLException("Cannot apply operation to datatype PTR.",true,false);
2132 return NULL;
2133 }
2134 template<>
ModNew(BaseGDL * r)2135 Data_<SpDObj>* Data_<SpDObj>::ModNew( BaseGDL* r)
2136 {
2137 throw GDLException("Cannot apply operation to datatype OBJECT.",true,false);
2138 return NULL;
2139 }
2140 template<>
ModInvNew(BaseGDL * r)2141 Data_<SpDObj>* Data_<SpDObj>::ModInvNew( BaseGDL* r)
2142 {
2143 throw GDLException("Cannot apply operation to datatype OBJECT.",true,false);
2144 return NULL;
2145 }
2146
2147 // scalar versions
2148 template<class Sp>
ModSNew(BaseGDL * r)2149 Data_<Sp>* Data_<Sp>::ModSNew( BaseGDL* r)
2150 {
2151 Data_* right=static_cast<Data_*>(r);
2152
2153 ULong nEl=N_Elements();
2154 assert( nEl);
2155
2156 Ty s = (*right)[0];
2157
2158 Data_* res = NewResult();
2159 if( s != this->zero)
2160 {
2161 for( SizeT i=0; i < nEl; ++i)
2162 (*res)[i] = (*this)[i] % s;
2163 return res;
2164 }
2165
2166
2167 if( sigsetjmp( sigFPEJmpBuf, 1) == 0)
2168 {
2169 for( SizeT i=0; i < nEl; ++i)
2170 (*res)[i] = (*this)[i] % s;
2171 return res;
2172 }
2173 else
2174 {
2175 assert( s == this->zero);
2176 for( SizeT i=0; i < nEl; ++i)
2177 (*res)[i] = this->zero;
2178 return res;
2179 }
2180 }
2181 // inverse modulo division: left=right % left
2182 template<class Sp>
ModInvSNew(BaseGDL * r)2183 Data_<Sp>* Data_<Sp>::ModInvSNew( BaseGDL* r)
2184 {
2185 Data_* right=static_cast<Data_*>(r);
2186
2187 ULong nEl=N_Elements();
2188 assert( nEl);
2189
2190 Data_* res = NewResult();
2191 if( nEl == 1 && (*this)[0] != this->zero)
2192 {
2193 (*res)[0] = (*right)[0] % (*this)[0];
2194 return res;
2195 }
2196
2197 Ty s = (*right)[0];
2198 SizeT i=0;
2199
2200 if( sigsetjmp( sigFPEJmpBuf, 1) == 0)
2201 {
2202 for( /*SizeT i=0*/; i < nEl; ++i)
2203 {
2204 (*res)[i] = s % (*this)[i];
2205 }
2206 return res;
2207 }
2208 else
2209 {
2210 TRACEOMP( __FILE__, __LINE__)
2211 #pragma omp parallel if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl))
2212 {
2213 #pragma omp for
2214 for( OMPInt ix=i; ix < nEl; ++ix)
2215 if( (*this)[ix] != this->zero)
2216 (*res)[ix] = s % (*this)[ix];
2217 else
2218 (*res)[ix] = this->zero;
2219 }
2220 return res;
2221 }
2222 }
2223 template<>
ModSNew(BaseGDL * r)2224 Data_<SpDFloat>* Data_<SpDFloat>::ModSNew( BaseGDL* r)
2225 {
2226 Data_* right=static_cast<Data_*>(r);
2227
2228 ULong nEl=N_Elements();
2229 assert( nEl);
2230
2231 Data_* res = NewResult();
2232 if( nEl == 1)
2233 {
2234 (*res)[0] = Modulo((*this)[0],(*right)[0]);
2235 return res;
2236 }
2237
2238 Ty s = (*right)[0];
2239 TRACEOMP( __FILE__, __LINE__)
2240 #pragma omp parallel if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl))
2241 {
2242 #pragma omp for
2243 for( OMPInt i=0; i < nEl; ++i)
2244 (*res)[i] = Modulo((*this)[i],s);
2245 }
2246 return res;
2247 }
2248 // float inverse modulo division: left=right % left
2249 template<>
ModInvSNew(BaseGDL * r)2250 Data_<SpDFloat>* Data_<SpDFloat>::ModInvSNew( BaseGDL* r)
2251 {
2252 Data_* right=static_cast<Data_*>(r);
2253
2254 ULong nEl=N_Elements();
2255 assert( nEl);
2256
2257 Data_* res = NewResult();
2258 if( nEl == 1)
2259 {
2260 (*res)[0] = Modulo((*right)[0],(*this)[0]);
2261 return res;
2262 }
2263
2264 Ty s = (*right)[0];
2265 TRACEOMP( __FILE__, __LINE__)
2266 #pragma omp parallel if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl))
2267 {
2268 #pragma omp for
2269 for( OMPInt i=0; i < nEl; ++i)
2270 (*res)[i] = Modulo(s,(*this)[i]);
2271 }
2272 return res;
2273 }
2274 template<>
ModSNew(BaseGDL * r)2275 Data_<SpDDouble>* Data_<SpDDouble>::ModSNew( BaseGDL* r)
2276 {
2277 Data_* right=static_cast<Data_*>(r);
2278
2279 ULong nEl=N_Elements();
2280 assert( nEl);
2281
2282 Data_* res = NewResult();
2283 if( nEl == 1)
2284 {
2285 (*res)[0] = DModulo((*this)[0],(*right)[0]);
2286 return res;
2287 }
2288
2289 Ty s = (*right)[0];
2290 TRACEOMP( __FILE__, __LINE__)
2291 #pragma omp parallel if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl))
2292 {
2293 #pragma omp for
2294 for( OMPInt i=0; i < nEl; ++i)
2295 (*res)[i] = DModulo((*this)[i],s);
2296 }
2297 return res;
2298 }
2299 // double inverse modulo division: left=right % left
2300 template<>
ModInvSNew(BaseGDL * r)2301 Data_<SpDDouble>* Data_<SpDDouble>::ModInvSNew( BaseGDL* r)
2302 {
2303 Data_* right=static_cast<Data_*>(r);
2304
2305 ULong nEl=N_Elements();
2306 assert( nEl);
2307
2308 Data_* res = NewResult();
2309 if( nEl == 1)
2310 {
2311 (*res)[0] = DModulo((*right)[0],(*this)[0]);
2312 return res;
2313 }
2314
2315 Ty s = (*right)[0];
2316 TRACEOMP( __FILE__, __LINE__)
2317 #pragma omp parallel if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl))
2318 {
2319 #pragma omp for
2320 for( OMPInt i=0; i < nEl; ++i)
2321 (*res)[i] = DModulo(s,(*this)[i]);
2322 }
2323 return res;
2324 }
2325 // invalid types
2326 template<>
ModSNew(BaseGDL * r)2327 Data_<SpDString>* Data_<SpDString>::ModSNew( BaseGDL* r)
2328
2329 {
2330 throw GDLException("Cannot apply operation to datatype STRING.",true,false);
2331 return NULL;
2332 }
2333 template<>
ModInvSNew(BaseGDL * r)2334 Data_<SpDString>* Data_<SpDString>::ModInvSNew( BaseGDL* r)
2335 {
2336 throw GDLException("Cannot apply operation to datatype STRING.",true,false);
2337 return NULL;
2338 }
2339 template<>
ModSNew(BaseGDL * r)2340 Data_<SpDComplex>* Data_<SpDComplex>::ModSNew( BaseGDL* r)
2341 {
2342 throw GDLException("Cannot apply operation to datatype "+str+".",true,false);
2343 return NULL;
2344 }
2345 template<>
ModSNew(BaseGDL * r)2346 Data_<SpDComplexDbl>* Data_<SpDComplexDbl>::ModSNew( BaseGDL* r)
2347 {
2348 throw GDLException("Cannot apply operation to datatype "+str+".",true,false);
2349 return NULL;
2350 }
2351 template<>
ModInvSNew(BaseGDL * r)2352 Data_<SpDComplex>* Data_<SpDComplex>::ModInvSNew( BaseGDL* r)
2353 {
2354 throw GDLException("Cannot apply operation to datatype "+str+".",true,false);
2355 return NULL;
2356 }
2357 template<>
ModInvSNew(BaseGDL * r)2358 Data_<SpDComplexDbl>* Data_<SpDComplexDbl>::ModInvSNew( BaseGDL* r)
2359 {
2360 throw GDLException("Cannot apply operation to datatype "+str+".",true,false);
2361 return NULL;
2362 }
2363 template<>
ModSNew(BaseGDL * r)2364 Data_<SpDPtr>* Data_<SpDPtr>::ModSNew( BaseGDL* r)
2365 {
2366 throw GDLException("Cannot apply operation to datatype PTR.",true,false);
2367 return NULL;
2368 }
2369 template<>
ModInvSNew(BaseGDL * r)2370 Data_<SpDPtr>* Data_<SpDPtr>::ModInvSNew( BaseGDL* r)
2371 {
2372 throw GDLException("Cannot apply operation to datatype PTR.",true,false);
2373 return NULL;
2374 }
2375 template<>
ModSNew(BaseGDL * r)2376 Data_<SpDObj>* Data_<SpDObj>::ModSNew( BaseGDL* r)
2377 {
2378 throw GDLException("Cannot apply operation to datatype OBJECT.",true,false);
2379 return NULL;
2380 }
2381 template<>
ModInvSNew(BaseGDL * r)2382 Data_<SpDObj>* Data_<SpDObj>::ModInvSNew( BaseGDL* r)
2383 {
2384 throw GDLException("Cannot apply operation to datatype OBJECT.",true,false);
2385 return NULL;
2386 }
2387
2388
2389
2390 // Pow pow pow pow pow pow pow pow pow pow pow pow pow pow pow pow pow pow
2391 // C++ defines pow only for floats and doubles
2392 // in basic_op.cpp:
2393 // template <typename T> T pow( const T r, const T l)
2394 // {
2395 // typedef T TT;
2396 //
2397 // if( l == 0) return 1;
2398 // if( l < 0) return 0;
2399 //
2400 // const int nBits = sizeof(TT) * 8;
2401 //
2402 // T arr = r;
2403 // T res = 1;
2404 // TT mask = 1;
2405 // for( SizeT i=0; i<nBits; ++i)
2406 // {
2407 // if( l & mask) res *= arr;
2408 // mask <<= 1;
2409 // if( l < mask) return res;
2410 // arr *= arr;
2411 // }
2412 //
2413 // return res;
2414 // }
2415
2416 // power of value: left=left ^ right
2417 // integral types
2418 template<class Sp>
PowNew(BaseGDL * r)2419 Data_<Sp>* Data_<Sp>::PowNew( BaseGDL* r)
2420 {
2421 Data_* right=static_cast<Data_*>(r);
2422
2423 ULong nEl=N_Elements();
2424 Data_* res = NewResult();
2425 assert( nEl);
2426 if( nEl == 1)
2427 {
2428 (*res)[0] = pow( (*this)[0], (*right)[0]);
2429 return res;
2430 }
2431 TRACEOMP( __FILE__, __LINE__)
2432 #pragma omp parallel if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl))
2433 {
2434 #pragma omp for
2435 for( OMPInt i=0; i < nEl; ++i)
2436 (*res)[i] = pow( (*this)[i], (*right)[i]); // valarray
2437 } //C delete right;
2438 return res;
2439 }
2440 // inverse power of value: left=right ^ left
2441 template<class Sp>
PowInvNew(BaseGDL * r)2442 Data_<Sp>* Data_<Sp>::PowInvNew( BaseGDL* r)
2443 {
2444 Data_* right=static_cast<Data_*>(r);
2445
2446 ULong nEl=N_Elements();
2447 assert( nEl);
2448 Data_* res = NewResult();
2449 if( nEl == 1)
2450 {
2451 (*res)[0] = pow( (*right)[0], (*this)[0]);
2452 return res;
2453 }
2454
2455 TRACEOMP( __FILE__, __LINE__)
2456 #pragma omp parallel if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl))
2457 {
2458 #pragma omp for
2459 for( OMPInt i=0; i < nEl; ++i)
2460 (*res)[i] = pow( (*right)[i], (*this)[i]);
2461 }
2462 return res;
2463 }
2464
2465 // PowInt and PowIntNew can only be called for GDL_FLOAT and GDL_DOUBLE
2466 template<class Sp>
PowIntNew(BaseGDL * r)2467 Data_<Sp>* Data_<Sp>::PowIntNew( BaseGDL* r)
2468 {
2469 assert( 0);
2470 throw GDLException("Internal error: Data_::PowIntNew called.",true,false);
2471 return NULL;
2472 }
2473 // floats power of value with GDL_LONG: left=left ^ right
2474 template<>
PowIntNew(BaseGDL * r)2475 Data_<SpDFloat>* Data_<SpDFloat>::PowIntNew( BaseGDL* r)
2476 {
2477 DLongGDL* right=static_cast<DLongGDL*>(r);
2478
2479 ULong rEl=right->N_Elements();
2480 ULong nEl=N_Elements();
2481 assert( rEl);
2482 assert( nEl);
2483 if( r->StrictScalar())
2484 {
2485 Data_* res = new Data_( Dim(), BaseGDL::NOZERO);
2486 DLong r0 = (*right)[0];
2487 TRACEOMP( __FILE__, __LINE__)
2488 #pragma omp parallel if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl))
2489 {
2490 #pragma omp for
2491 for( OMPInt i=0; i < nEl; ++i)
2492 (*res)[i] = pow( (*this)[i], r0);
2493 } return res;
2494 }
2495 if( StrictScalar())
2496 {
2497 Data_* res = new Data_( right->Dim(), BaseGDL::NOZERO);
2498 Ty s0 = (*this)[ 0];
2499 TRACEOMP( __FILE__, __LINE__)
2500 #pragma omp parallel if (rEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= rEl))
2501 {
2502 #pragma omp for
2503 for( OMPInt i=0; i < rEl; ++i)
2504 (*res)[ i] = pow( s0, (*right)[ i]);
2505 } return res;
2506 }
2507 if( nEl <= rEl)
2508 {
2509 Data_* res = new Data_( Dim(), BaseGDL::NOZERO);
2510 TRACEOMP( __FILE__, __LINE__)
2511 #pragma omp parallel if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl))
2512 {
2513 #pragma omp for
2514 for( OMPInt i=0; i < nEl; ++i)
2515 (*res)[i] = pow( (*this)[i], (*right)[i]);
2516 } return res;
2517 }
2518 else
2519 {
2520 Data_* res = new Data_( right->Dim(), BaseGDL::NOZERO);
2521 TRACEOMP( __FILE__, __LINE__)
2522 #pragma omp parallel if (rEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= rEl))
2523 {
2524 #pragma omp for
2525 for( OMPInt i=0; i < rEl; ++i)
2526 (*res)[i] = pow( (*this)[i], (*right)[i]);
2527 } return res;
2528 }
2529 }
2530 template<>
PowIntNew(BaseGDL * r)2531 Data_<SpDDouble>* Data_<SpDDouble>::PowIntNew( BaseGDL* r)
2532 {
2533 DLongGDL* right=static_cast<DLongGDL*>(r);
2534
2535 ULong rEl=right->N_Elements();
2536 ULong nEl=N_Elements();
2537 assert( rEl);
2538 assert( nEl);
2539 if( r->StrictScalar())
2540 {
2541 Data_* res = new Data_( Dim(), BaseGDL::NOZERO);
2542 DLong r0 = (*right)[0];
2543 TRACEOMP( __FILE__, __LINE__)
2544 #pragma omp parallel if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl))
2545 {
2546 #pragma omp for
2547 for( OMPInt i=0; i < nEl; ++i)
2548 (*res)[i] = pow( (*this)[i], r0);
2549 } return res;
2550 }
2551 if( StrictScalar())
2552 {
2553 Data_* res = new Data_( right->Dim(), BaseGDL::NOZERO);
2554 Ty s0 = (*this)[ 0];
2555 TRACEOMP( __FILE__, __LINE__)
2556 #pragma omp parallel if (rEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= rEl))
2557 {
2558 #pragma omp for
2559 for( OMPInt i=0; i < rEl; ++i)
2560 (*res)[ i] = pow( s0, (*right)[ i]);
2561 } return res;
2562 }
2563 if( nEl <= rEl)
2564 {
2565 Data_* res = new Data_( Dim(), BaseGDL::NOZERO);
2566 TRACEOMP( __FILE__, __LINE__)
2567 #pragma omp parallel if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl))
2568 {
2569 #pragma omp for
2570 for( OMPInt i=0; i < nEl; ++i)
2571 (*res)[i] = pow( (*this)[i], (*right)[i]);
2572 } return res;
2573 }
2574 else
2575 {
2576 Data_* res = new Data_( right->Dim(), BaseGDL::NOZERO);
2577 TRACEOMP( __FILE__, __LINE__)
2578 #pragma omp parallel if (rEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= rEl))
2579 {
2580 #pragma omp for
2581 for( OMPInt i=0; i < rEl; ++i)
2582 (*res)[i] = pow( (*this)[i], (*right)[i]);
2583 } return res;
2584 }
2585 }
2586
2587 // floats power of value: left=left ^ right
2588 template<>
PowNew(BaseGDL * r)2589 Data_<SpDFloat>* Data_<SpDFloat>::PowNew( BaseGDL* r)
2590 {
2591 Data_* right=static_cast<Data_*>(r);
2592
2593 ULong nEl=N_Elements();
2594 assert( nEl);
2595 Data_* res = NewResult();
2596 if( nEl == 1)
2597 {
2598 (*res)[0] = pow( (*this)[0], (*right)[0]);
2599 return res;
2600 }
2601 TRACEOMP( __FILE__, __LINE__)
2602 #pragma omp parallel if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl))
2603 {
2604 #pragma omp for
2605 for( OMPInt i=0; i < nEl; ++i)
2606 (*res)[i] = pow( (*this)[i], (*right)[i]);
2607 }
2608 return res;
2609 }
2610 // floats inverse power of value: left=right ^ left
2611 template<>
PowInvNew(BaseGDL * r)2612 Data_<SpDFloat>* Data_<SpDFloat>::PowInvNew( BaseGDL* r)
2613 {
2614 Data_* right=static_cast<Data_*>(r);
2615
2616 ULong nEl=N_Elements();
2617 assert( nEl);
2618 Data_* res = NewResult();
2619 if( nEl == 1)
2620 {
2621 (*res)[0] = pow( (*right)[0], (*this)[0]);
2622 return res;
2623 }
2624 TRACEOMP( __FILE__, __LINE__)
2625 #pragma omp parallel if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl))
2626 {
2627 #pragma omp for
2628 for( OMPInt i=0; i < nEl; ++i)
2629 (*res)[i] = pow( (*right)[i], (*this)[i]);
2630 } //C delete right;
2631 return res;
2632 }
2633 // doubles power of value: left=left ^ right
2634 template<>
PowNew(BaseGDL * r)2635 Data_<SpDDouble>* Data_<SpDDouble>::PowNew( BaseGDL* r)
2636 {
2637 Data_* right=static_cast<Data_*>(r);
2638
2639 ULong nEl=N_Elements();
2640 assert( nEl);
2641 Data_* res = NewResult();
2642 if( nEl == 1)
2643 {
2644 (*res)[0] = pow( (*this)[0], (*right)[0]);
2645 return res;
2646 }
2647 TRACEOMP( __FILE__, __LINE__)
2648 #pragma omp parallel if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl))
2649 {
2650 #pragma omp for
2651 for( OMPInt i=0; i < nEl; ++i)
2652 (*res)[i] = pow( (*this)[i], (*right)[i]);
2653 }
2654 return res;
2655 }
2656 // doubles inverse power of value: left=right ^ left
2657 template<>
PowInvNew(BaseGDL * r)2658 Data_<SpDDouble>* Data_<SpDDouble>::PowInvNew( BaseGDL* r)
2659 {
2660 Data_* right=static_cast<Data_*>(r);
2661
2662 ULong nEl=N_Elements();
2663 Data_* res = NewResult();
2664 assert( nEl);
2665 if( nEl == 1)
2666 {
2667 (*res)[0] = pow( (*right)[0], (*this)[0]);
2668 return res;
2669 }
2670
2671 TRACEOMP( __FILE__, __LINE__)
2672 #pragma omp parallel if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl))
2673 {
2674 #pragma omp for
2675 for( OMPInt i=0; i < nEl; ++i)
2676 (*res)[i] = pow( (*right)[i], (*this)[i]);
2677 } //C delete right;
2678 return res;
2679 }
2680 // complex power of value: left=left ^ right
2681 // complex is special here
2682 template<>
PowNew(BaseGDL * r)2683 Data_<SpDComplex>* Data_<SpDComplex>::PowNew( BaseGDL* r)
2684 {
2685 SizeT nEl = N_Elements();
2686
2687 assert( nEl > 0);
2688 assert( r->N_Elements() > 0);
2689
2690 if( r->Type() == GDL_FLOAT)
2691 {
2692 Data_<SpDFloat>* right=static_cast<Data_<SpDFloat>* >(r);
2693
2694 DFloat s;
2695 // note: changes here have to be reflected in POWNCNode::Eval() (dnode.cpp)
2696 // (concerning when a new variable is created vs. using this)
2697 // (must also be consistent with ComplexDbl)
2698 if( right->StrictScalar(s))
2699 {
2700 DComplexGDL* res = new DComplexGDL( this->Dim(),
2701 BaseGDL::NOZERO);
2702 TRACEOMP( __FILE__, __LINE__)
2703 #pragma omp parallel if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl))
2704 {
2705 #pragma omp for
2706 for( OMPInt i=0; i<nEl; ++i)
2707 (*res)[ i] = pow( (*this)[ i], s);
2708 } //C delete right;
2709 return res;
2710 }
2711 else
2712 {
2713 SizeT rEl = right->N_Elements();
2714 if( nEl < rEl)
2715 {
2716 DComplex s;
2717 if( StrictScalar(s))
2718 {
2719 DComplexGDL* res = new DComplexGDL( right->Dim(),
2720 BaseGDL::NOZERO);
2721 TRACEOMP( __FILE__, __LINE__)
2722 #pragma omp parallel if (rEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= rEl))
2723 {
2724 #pragma omp for
2725 for( OMPInt i=0; i<rEl; ++i)
2726 (*res)[ i] = pow( s, (*right)[ i]);
2727 } //C delete right;
2728 return res;
2729 }
2730
2731 DComplexGDL* res = new DComplexGDL( this->Dim(),
2732 BaseGDL::NOZERO);
2733 TRACEOMP( __FILE__, __LINE__)
2734 #pragma omp parallel if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl))
2735 {
2736 #pragma omp for
2737 for( OMPInt i=0; i<nEl; ++i)
2738 (*res)[ i] = pow( (*this)[ i], (*right)[ i]);
2739 } //C delete right;
2740 return res;
2741 }
2742 else
2743 {
2744 DComplexGDL* res = new DComplexGDL( right->Dim(),
2745 BaseGDL::NOZERO);
2746 TRACEOMP( __FILE__, __LINE__)
2747 #pragma omp parallel if (rEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= rEl))
2748 {
2749 #pragma omp for
2750 for( OMPInt i=0; i<rEl; ++i)
2751 (*res)[ i] = pow( (*this)[ i], (*right)[ i]);
2752 } //C delete right;
2753 //C delete this;
2754 return res;
2755 }
2756 }
2757 }
2758 if( r->Type() == GDL_LONG)
2759 {
2760 Data_<SpDLong>* right=static_cast<Data_<SpDLong>* >(r);
2761
2762 DLong s;
2763 // note: changes here have to be reflected in POWNCNode::Eval() (dnode.cpp)
2764 // (concerning when a new variable is created vs. using this)
2765 // (must also be consistent with ComplexDbl)
2766 if( right->StrictScalar(s))
2767 {
2768 DComplexGDL* res = new DComplexGDL( this->Dim(),
2769 BaseGDL::NOZERO);
2770 TRACEOMP( __FILE__, __LINE__)
2771 #pragma omp parallel if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl))
2772 {
2773 #pragma omp for
2774 for( OMPInt i=0; i<nEl; ++i)
2775 (*res)[ i] = pow( (*this)[ i], s);
2776 } //C delete right;
2777 return res;
2778 }
2779 else
2780 {
2781 SizeT rEl = right->N_Elements();
2782 if( nEl < rEl)
2783 {
2784 DComplex s;
2785 if( StrictScalar(s))
2786 {
2787 DComplexGDL* res = new DComplexGDL( right->Dim(),
2788 BaseGDL::NOZERO);
2789 TRACEOMP( __FILE__, __LINE__)
2790 #pragma omp parallel if (rEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= rEl))
2791 {
2792 #pragma omp for
2793 for( OMPInt i=0; i<rEl; ++i)
2794 (*res)[ i] = pow( s, (*right)[ i]);
2795 } //C delete right;
2796 return res;
2797 }
2798
2799 DComplexGDL* res = new DComplexGDL( this->Dim(),
2800 BaseGDL::NOZERO);
2801 TRACEOMP( __FILE__, __LINE__)
2802 #pragma omp parallel if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl))
2803 {
2804 #pragma omp for
2805 for( OMPInt i=0; i<nEl; ++i)
2806 (*res)[ i] = pow( (*this)[ i], (*right)[ i]);
2807 } //C delete right;
2808 return res;
2809 }
2810 else
2811 {
2812 DComplexGDL* res = new DComplexGDL( right->Dim(),
2813 BaseGDL::NOZERO);
2814 TRACEOMP( __FILE__, __LINE__)
2815 #pragma omp parallel if (rEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= rEl))
2816 {
2817 #pragma omp for
2818 for( OMPInt i=0; i<rEl; ++i)
2819 (*res)[ i] = pow( (*this)[ i], (*right)[ i]);
2820 } //C delete right;
2821 //C delete this;
2822 return res;
2823 }
2824 }
2825 }
2826
2827 Data_* right=static_cast<Data_*>(r);
2828
2829 // ULong rEl=right->N_Elements();
2830 // ULong nEl=N_Elements(); Data_* res = NewResult();
2831 // if( !rEl || !nEl) throw GDLException("Variable is undefined.");
2832 Ty s;
2833 if( right->StrictScalar(s))
2834 {
2835 DComplexGDL* res = new DComplexGDL( this->Dim(),
2836 BaseGDL::NOZERO);
2837 //#if (__GNUC__ == 3) && (__GNUC_MINOR__ <= 2)
2838 TRACEOMP( __FILE__, __LINE__)
2839 #pragma omp parallel if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl))
2840 {
2841 #pragma omp for
2842 for( OMPInt i=0; i<nEl; ++i)
2843 (*res)[ i] = pow( (*this)[ i], s);
2844 }
2845 return res;
2846 }
2847 else
2848 {
2849 DComplexGDL* res = new DComplexGDL( this->Dim(),
2850 BaseGDL::NOZERO);
2851 TRACEOMP( __FILE__, __LINE__)
2852 #pragma omp parallel if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl))
2853 {
2854 #pragma omp for
2855 for( OMPInt i=0; i < nEl; ++i)
2856 (*res)[i] = pow( (*this)[i], (*right)[i]);
2857 }
2858 return res;
2859 }
2860 }
2861 // complex inverse power of value: left=right ^ left
2862 template<>
PowInvNew(BaseGDL * r)2863 Data_<SpDComplex>* Data_<SpDComplex>::PowInvNew( BaseGDL* r)
2864 {
2865 Data_* right=static_cast<Data_*>(r);
2866
2867 ULong nEl=N_Elements();
2868 assert( nEl);
2869
2870 Data_* res = NewResult();
2871 TRACEOMP( __FILE__, __LINE__)
2872 #pragma omp parallel if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl))
2873 {
2874 #pragma omp for
2875 for( OMPInt i=0; i < nEl; ++i)
2876 (*res)[i] = pow( (*right)[i], (*this)[i]);
2877 }
2878 return res;
2879 }
2880 // double complex power of value: left=left ^ right
2881 template<>
PowNew(BaseGDL * r)2882 Data_<SpDComplexDbl>* Data_<SpDComplexDbl>::PowNew( BaseGDL* r)
2883 {
2884 SizeT nEl = N_Elements();
2885
2886 assert( nEl > 0);
2887 assert( r->N_Elements() > 0);
2888
2889 if( r->Type() == GDL_DOUBLE)
2890 {
2891 Data_<SpDDouble>* right=static_cast<Data_<SpDDouble>* >(r);
2892
2893 DDouble s;
2894 // note: changes here have to be reflected in POWNCNode::Eval() (dnode.cpp)
2895 // (concerning when a new variable is created vs. using this)
2896 // (must also be consistent with ComplexDbl)
2897 if( right->StrictScalar(s))
2898 {
2899 DComplexDblGDL* res = new DComplexDblGDL( this->Dim(),
2900 BaseGDL::NOZERO);
2901 TRACEOMP( __FILE__, __LINE__)
2902 #pragma omp parallel if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl))
2903 {
2904 #pragma omp for
2905 for( OMPInt i=0; i<nEl; ++i)
2906 (*res)[ i] = pow( (*this)[ i], s);
2907 }
2908 return res;
2909 }
2910 else
2911 {
2912 SizeT rEl = right->N_Elements();
2913 if( nEl < rEl)
2914 {
2915 DComplexDbl s;
2916 if( StrictScalar(s))
2917 {
2918 DComplexDblGDL* res = new DComplexDblGDL( right->Dim(),
2919 BaseGDL::NOZERO);
2920 TRACEOMP( __FILE__, __LINE__)
2921 #pragma omp parallel if (rEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= rEl))
2922 {
2923 #pragma omp for
2924 for( OMPInt i=0; i<rEl; ++i)
2925 (*res)[ i] = pow( s, (*right)[ i]);
2926 }
2927 return res;
2928 }
2929
2930 DComplexDblGDL* res = new DComplexDblGDL( this->Dim(),
2931 BaseGDL::NOZERO);
2932 TRACEOMP( __FILE__, __LINE__)
2933 #pragma omp parallel if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl))
2934 {
2935 #pragma omp for
2936 for( OMPInt i=0; i<nEl; ++i)
2937 (*res)[ i] = pow( (*this)[ i], (*right)[ i]);
2938 }
2939 return res;
2940 }
2941 else
2942 {
2943 DComplexDblGDL* res = new DComplexDblGDL( right->Dim(),
2944 BaseGDL::NOZERO);
2945 TRACEOMP( __FILE__, __LINE__)
2946 #pragma omp parallel if (rEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= rEl))
2947 {
2948 #pragma omp for
2949 for( OMPInt i=0; i<rEl; ++i)
2950 (*res)[ i] = pow( (*this)[ i], (*right)[ i]);
2951 }
2952 return res;
2953 }
2954 }
2955 }
2956 if( r->Type() == GDL_LONG)
2957 {
2958 Data_<SpDLong>* right=static_cast<Data_<SpDLong>* >(r);
2959
2960 DLong s;
2961 // note: changes here have to be reflected in POWNCNode::Eval() (dnode.cpp)
2962 // (concerning when a new variable is created vs. using this)
2963 // (must also be consistent with ComplexDbl)
2964 if( right->StrictScalar(s))
2965 {
2966 DComplexDblGDL* res = new DComplexDblGDL( this->Dim(),
2967 BaseGDL::NOZERO);
2968 TRACEOMP( __FILE__, __LINE__)
2969 #pragma omp parallel if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl))
2970 {
2971 #pragma omp for
2972 for( OMPInt i=0; i<nEl; ++i)
2973 (*res)[ i] = pow( (*this)[ i], s);
2974 }
2975 return res;
2976 }
2977 else
2978 {
2979 SizeT rEl = right->N_Elements();
2980 if( nEl < rEl)
2981 {
2982 DComplexDbl s;
2983 if( StrictScalar(s))
2984 {
2985 DComplexDblGDL* res = new DComplexDblGDL( right->Dim(),
2986 BaseGDL::NOZERO);
2987 TRACEOMP( __FILE__, __LINE__)
2988 #pragma omp parallel if (rEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= rEl))
2989 {
2990 #pragma omp for
2991 for( OMPInt i=0; i<rEl; ++i)
2992 (*res)[ i] = pow( s, (*right)[ i]);
2993 }
2994 return res;
2995 }
2996
2997 DComplexDblGDL* res = new DComplexDblGDL( this->Dim(),
2998 BaseGDL::NOZERO);
2999 TRACEOMP( __FILE__, __LINE__)
3000 #pragma omp parallel if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl))
3001 {
3002 #pragma omp for
3003 for( OMPInt i=0; i<nEl; ++i)
3004 (*res)[ i] = pow( (*this)[ i], (*right)[ i]);
3005 } //C delete right;
3006 return res;
3007 }
3008 else
3009 {
3010 DComplexDblGDL* res = new DComplexDblGDL( right->Dim(),
3011 BaseGDL::NOZERO);
3012 TRACEOMP( __FILE__, __LINE__)
3013 #pragma omp parallel if (rEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= rEl))
3014 {
3015 #pragma omp for
3016 for( OMPInt i=0; i<rEl; ++i)
3017 (*res)[ i] = pow( (*this)[ i], (*right)[ i]);
3018 }
3019 return res;
3020 }
3021 }
3022 }
3023
3024 Data_* right=static_cast<Data_*>(r);
3025
3026 Ty s;
3027 if( right->StrictScalar(s))
3028 {
3029 DComplexDblGDL* res = new DComplexDblGDL( this->Dim(),
3030 BaseGDL::NOZERO);
3031 TRACEOMP( __FILE__, __LINE__)
3032 #pragma omp parallel if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl))
3033 {
3034 #pragma omp for
3035 for( OMPInt i=0; i<nEl; ++i)
3036 (*res)[ i] = pow( (*this)[ i], s);
3037 }
3038 return res;
3039 }
3040 else
3041 {
3042 DComplexDblGDL* res = new DComplexDblGDL( this->Dim(),
3043 BaseGDL::NOZERO);
3044 TRACEOMP( __FILE__, __LINE__)
3045 #pragma omp parallel if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl))
3046 {
3047 #pragma omp for
3048 for( OMPInt i=0; i < nEl; ++i)
3049 (*res)[i] = pow( (*this)[i], (*right)[i]);
3050 }
3051 return res;
3052 }
3053 }
3054 // double complex inverse power of value: left=right ^ left
3055 template<>
PowInvNew(BaseGDL * r)3056 Data_<SpDComplexDbl>* Data_<SpDComplexDbl>::PowInvNew( BaseGDL* r)
3057 {
3058 Data_* right=static_cast<Data_*>(r);
3059
3060 ULong rEl=right->N_Elements();
3061 ULong nEl=N_Elements();
3062 assert( rEl);
3063 assert( nEl);
3064 Data_* res = NewResult();
3065 TRACEOMP( __FILE__, __LINE__)
3066 #pragma omp parallel if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl))
3067 {
3068 #pragma omp for
3069 for( OMPInt i=0; i < nEl; ++i)
3070 (*res)[i] = pow( (*right)[i], (*this)[i]);
3071 }
3072 return res;
3073 }
3074 // invalid types
3075 template<>
PowNew(BaseGDL * r)3076 Data_<SpDString>* Data_<SpDString>::PowNew( BaseGDL* r)
3077 {
3078 throw GDLException("Cannot apply operation to datatype STRING.",true,false);
3079 return NULL;
3080 }
3081 template<>
PowInvNew(BaseGDL * r)3082 Data_<SpDString>* Data_<SpDString>::PowInvNew( BaseGDL* r)
3083 {
3084 throw GDLException("Cannot apply operation to datatype STRING.",true,false);
3085 return NULL;
3086 }
3087 template<>
PowNew(BaseGDL * r)3088 Data_<SpDPtr>* Data_<SpDPtr>::PowNew( BaseGDL* r)
3089 {
3090 throw GDLException("Cannot apply operation to datatype PTR.",true,false);
3091 return NULL;
3092 }
3093 template<>
PowInvNew(BaseGDL * r)3094 Data_<SpDPtr>* Data_<SpDPtr>::PowInvNew( BaseGDL* r)
3095 {
3096 throw GDLException("Cannot apply operation to datatype PTR.",true,false);
3097 return NULL;
3098 }
3099 template<>
PowNew(BaseGDL * r)3100 Data_<SpDObj>* Data_<SpDObj>::PowNew( BaseGDL* r)
3101 {
3102 throw GDLException("Cannot apply operation to datatype OBJECT.",true,false);
3103 return NULL;
3104 }
3105 template<>
PowInvNew(BaseGDL * r)3106 Data_<SpDObj>* Data_<SpDObj>::PowInvNew( BaseGDL* r)
3107 {
3108 throw GDLException("Cannot apply operation to datatype OBJECT.",true,false);
3109 return NULL;
3110 }
3111
3112 // scalar versions
3113 template<class Sp>
PowSNew(BaseGDL * r)3114 Data_<Sp>* Data_<Sp>::PowSNew( BaseGDL* r)
3115 {
3116 Data_* right=static_cast<Data_*>(r);
3117
3118 ULong nEl=N_Elements();
3119 Data_* res = NewResult();
3120 assert( nEl);
3121 Ty s = (*right)[0];
3122 if( nEl == 1)
3123 {
3124 (*res)[0] = pow( (*this)[0], s);
3125 return res;
3126 }
3127 TRACEOMP( __FILE__, __LINE__)
3128 #pragma omp parallel if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl))
3129 {
3130 #pragma omp for
3131 for( OMPInt i=0; i < nEl; ++i)
3132 (*res)[i] = pow( (*this)[i], s);
3133 }
3134 //C delete right;
3135 return res;
3136 }
3137 // inverse power of value: left=right ^ left
3138 template<class Sp>
PowInvSNew(BaseGDL * r)3139 Data_<Sp>* Data_<Sp>::PowInvSNew( BaseGDL* r)
3140 {
3141 Data_* right=static_cast<Data_*>(r);
3142
3143 ULong nEl=N_Elements();
3144 assert( nEl);
3145 Ty s = (*right)[0];
3146 Data_* res = NewResult();
3147 if( nEl == 1)
3148 {
3149 (*res)[0] = pow( s, (*this)[0]);
3150 return res;
3151 }
3152 TRACEOMP( __FILE__, __LINE__)
3153 #pragma omp parallel if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl))
3154 {
3155 #pragma omp for
3156 for( OMPInt i=0; i < nEl; ++i)
3157 (*res)[i] = pow( s, (*this)[i]);
3158 } //C delete right;
3159 return res;
3160 }
3161 // complex power of value: left=left ^ right
3162 // complex is special here
3163 template<>
PowSNew(BaseGDL * r)3164 Data_<SpDComplex>* Data_<SpDComplex>::PowSNew( BaseGDL* r)
3165 {
3166 SizeT nEl = N_Elements();
3167 assert( nEl > 0);
3168 assert( r->N_Elements() > 0);
3169
3170 if( r->Type() == GDL_FLOAT)
3171 {
3172 Data_<SpDFloat>* right=static_cast<Data_<SpDFloat>* >(r);
3173
3174 DFloat s;
3175 // note: changes here have to be reflected in POWNCNode::Eval() (dnode.cpp)
3176 // (concerning when a new variable is created vs. using this)
3177 // (must also be consistent with ComplexDbl)
3178 if( right->StrictScalar(s))
3179 {
3180 Data_* res = NewResult();
3181 TRACEOMP( __FILE__, __LINE__)
3182 #pragma omp parallel if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl))
3183 {
3184 #pragma omp for
3185 for( OMPInt i=0; i<nEl; ++i)
3186 (*res)[i] = pow( (*this)[ i], s);
3187 }
3188 return res;
3189 }
3190 else
3191 {
3192 SizeT rEl = right->N_Elements();
3193 if( nEl < rEl)
3194 {
3195 DComplex s;
3196 if( StrictScalar(s))
3197 {
3198 DComplexGDL* res = new DComplexGDL( right->Dim(),
3199 BaseGDL::NOZERO);
3200 TRACEOMP( __FILE__, __LINE__)
3201 #pragma omp parallel if (rEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= rEl))
3202 {
3203 #pragma omp for
3204 for( OMPInt i=0; i<rEl; ++i)
3205 (*res)[ i] = pow( s, (*right)[ i]);
3206 }
3207 return res;
3208 }
3209
3210 Data_* res = NewResult();
3211 TRACEOMP( __FILE__, __LINE__)
3212 #pragma omp parallel if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl))
3213 {
3214 #pragma omp for
3215 for( OMPInt i=0; i<nEl; ++i)
3216 (*res)[i] = pow( (*this)[ i], (*right)[ i]);
3217 } //C delete right;
3218 return res;
3219 }
3220 else
3221 {
3222 DComplexGDL* res = new DComplexGDL( right->Dim(),
3223 BaseGDL::NOZERO);
3224 TRACEOMP( __FILE__, __LINE__)
3225 #pragma omp parallel if (rEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= rEl))
3226 {
3227 #pragma omp for
3228 for( OMPInt i=0; i<rEl; ++i)
3229 (*res)[ i] = pow( (*this)[ i], (*right)[ i]);
3230 } //C delete right;
3231 //C delete this;
3232 return res;
3233 }
3234 }
3235 }
3236 if( r->Type() == GDL_LONG)
3237 {
3238 Data_<SpDLong>* right=static_cast<Data_<SpDLong>* >(r);
3239
3240 DLong s;
3241 // note: changes here have to be reflected in POWNCNode::Eval() (dnode.cpp)
3242 // (concerning when a new variable is created vs. using this)
3243 // (must also be consistent with ComplexDbl)
3244 if( right->StrictScalar(s))
3245 {
3246 Data_* res = NewResult();
3247 TRACEOMP( __FILE__, __LINE__)
3248 #pragma omp parallel if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl))
3249 {
3250 #pragma omp for
3251 for( OMPInt i=0; i<nEl; ++i)
3252 (*res)[i] = pow( (*this)[ i], s);
3253 } //C delete right;
3254 return res;
3255 }
3256 else
3257 {
3258 SizeT rEl = right->N_Elements();
3259 if( nEl < rEl)
3260 {
3261 DComplex s;
3262 if( StrictScalar(s))
3263 {
3264 DComplexGDL* res = new DComplexGDL( right->Dim(),
3265 BaseGDL::NOZERO);
3266 TRACEOMP( __FILE__, __LINE__)
3267 #pragma omp parallel if (rEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= rEl))
3268 {
3269 #pragma omp for
3270 for( OMPInt i=0; i<rEl; ++i)
3271 (*res)[ i] = pow( s, (*right)[ i]);
3272 } //C delete right;
3273 return res;
3274 }
3275
3276 Data_* res = NewResult();
3277 TRACEOMP( __FILE__, __LINE__)
3278 #pragma omp parallel if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl))
3279 {
3280 #pragma omp for
3281 for( OMPInt i=0; i<nEl; ++i)
3282 (*res)[i] = pow( (*this)[ i], (*right)[ i]);
3283 } //C delete right;
3284 return res;
3285 }
3286 else
3287 {
3288 DComplexGDL* res = new DComplexGDL( right->Dim(),
3289 BaseGDL::NOZERO);
3290 TRACEOMP( __FILE__, __LINE__)
3291 #pragma omp parallel if (rEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= rEl))
3292 {
3293 #pragma omp for
3294 for( OMPInt i=0; i<rEl; ++i)
3295 (*res)[ i] = pow( (*this)[ i], (*right)[ i]);
3296 }
3297 return res;
3298 }
3299 }
3300 }
3301
3302 // r->Type() == GDL_COMPLEX
3303 Data_* right=static_cast<Data_*>(r);
3304
3305 Ty s = (*right)[0];
3306 Data_* res = NewResult();
3307 TRACEOMP( __FILE__, __LINE__)
3308 #pragma omp parallel if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl))
3309 {
3310 #pragma omp for
3311 for( OMPInt i=0; i<nEl; ++i)
3312 (*res)[i] = pow( (*this)[ i], s);
3313 }
3314
3315 return res;
3316 }
3317 // complex inverse power of value: left=right ^ left
3318 template<>
PowInvSNew(BaseGDL * r)3319 Data_<SpDComplex>* Data_<SpDComplex>::PowInvSNew( BaseGDL* r)
3320 {
3321 Data_* right=static_cast<Data_*>(r);
3322 ULong nEl=N_Elements();
3323 assert( nEl);
3324 assert( right->N_Elements());
3325 Ty s = (*right)[0];
3326 Data_* res = NewResult();
3327 if( nEl == 1)
3328 {
3329 (*res)[0] = pow( s, (*this)[0]);
3330 return res;
3331 }
3332 TRACEOMP( __FILE__, __LINE__)
3333 #pragma omp parallel if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl))
3334 {
3335 #pragma omp for
3336 for( OMPInt i=0; i<nEl; ++i)
3337 (*res)[i] = pow( s, (*this)[ i]);
3338 }
3339 return res;
3340 }
3341 // double complex power of value: left=left ^ right
3342 template<>
PowSNew(BaseGDL * r)3343 Data_<SpDComplexDbl>* Data_<SpDComplexDbl>::PowSNew( BaseGDL* r)
3344 {
3345 SizeT nEl = N_Elements();
3346
3347 assert( nEl > 0);
3348
3349 if( r->Type() == GDL_DOUBLE)
3350 {
3351 Data_<SpDDouble>* right=static_cast<Data_<SpDDouble>* >(r);
3352
3353 assert( right->N_Elements() > 0);
3354
3355 DDouble s;
3356
3357 // note: changes here have to be reflected in POWNCNode::Eval() (dnode.cpp)
3358 // (concerning when a new variable is created vs. using this)
3359 if( right->StrictScalar(s))
3360 {
3361 Data_* res = NewResult();
3362 TRACEOMP( __FILE__, __LINE__)
3363 #pragma omp parallel if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl))
3364 {
3365 #pragma omp for
3366 for( OMPInt i=0; i<nEl; ++i)
3367 (*res)[i] = pow( (*this)[ i], s);
3368 }
3369 return res;
3370 }
3371 else
3372 {
3373 SizeT rEl = right->N_Elements();
3374 if( nEl < rEl)
3375 {
3376 DComplexDbl s;
3377 if( StrictScalar(s))
3378 {
3379 DComplexDblGDL* res = new DComplexDblGDL( right->Dim(),
3380 BaseGDL::NOZERO);
3381 TRACEOMP( __FILE__, __LINE__)
3382 #pragma omp parallel if (rEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= rEl))
3383 {
3384 #pragma omp for
3385 for( OMPInt i=0; i<rEl; ++i)
3386 (*res)[ i] = pow( s, (*right)[ i]);
3387 }
3388 return res;
3389 }
3390
3391 Data_* res = NewResult();
3392 TRACEOMP( __FILE__, __LINE__)
3393 #pragma omp parallel if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl))
3394 {
3395 #pragma omp for
3396 for( OMPInt i=0; i<nEl; ++i)
3397 (*res)[i] = pow( (*this)[ i], (*right)[ i]);
3398 } //C delete right;
3399 return res;
3400 }
3401 else
3402 {
3403 DComplexDblGDL* res = new DComplexDblGDL( right->Dim(),
3404 BaseGDL::NOZERO);
3405 TRACEOMP( __FILE__, __LINE__)
3406 #pragma omp parallel if (rEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= rEl))
3407 {
3408 #pragma omp for
3409 for( OMPInt i=0; i<rEl; ++i)
3410 (*res)[ i] = pow( (*this)[ i], (*right)[ i]);
3411 }
3412 return res;
3413 }
3414 }
3415 }
3416 if( r->Type() == GDL_LONG)
3417 {
3418 Data_<SpDLong>* right=static_cast<Data_<SpDLong>* >(r);
3419
3420 assert( right->N_Elements() > 0);
3421
3422 DLong s;
3423
3424 // note: changes here have to be reflected in POWNCNode::Eval() (dnode.cpp)
3425 // (concerning when a new variable is created vs. using this)
3426 if( right->StrictScalar(s))
3427 {
3428 Data_* res = NewResult();
3429 TRACEOMP( __FILE__, __LINE__)
3430 #pragma omp parallel if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl))
3431 {
3432 #pragma omp for
3433 for( OMPInt i=0; i<nEl; ++i)
3434 (*res)[i] = pow( (*this)[ i], s);
3435 }
3436 return res;
3437 }
3438 else
3439 {
3440 SizeT rEl = right->N_Elements();
3441 if( nEl < rEl)
3442 {
3443 DComplexDbl s;
3444 if( StrictScalar(s))
3445 {
3446 DComplexDblGDL* res = new DComplexDblGDL( right->Dim(),
3447 BaseGDL::NOZERO);
3448 TRACEOMP( __FILE__, __LINE__)
3449 #pragma omp parallel if (rEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= rEl))
3450 {
3451 #pragma omp for
3452 for( OMPInt i=0; i<rEl; ++i)
3453 (*res)[ i] = pow( s, (*right)[ i]);
3454 }
3455 return res;
3456 }
3457
3458 Data_* res = NewResult();
3459 TRACEOMP( __FILE__, __LINE__)
3460 #pragma omp parallel if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl))
3461 {
3462 #pragma omp for
3463 for( OMPInt i=0; i<nEl; ++i)
3464 (*res)[i] = pow( (*this)[ i], (*right)[ i]);
3465 } //C delete right;
3466 return res;
3467 }
3468 else
3469 {
3470 DComplexDblGDL* res = new DComplexDblGDL( right->Dim(),
3471 BaseGDL::NOZERO);
3472 TRACEOMP( __FILE__, __LINE__)
3473 #pragma omp parallel if (rEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= rEl))
3474 {
3475 #pragma omp for
3476 for( OMPInt i=0; i<rEl; ++i)
3477 (*res)[ i] = pow( (*this)[ i], (*right)[ i]);
3478 }
3479 return res;
3480 }
3481 }
3482 }
3483
3484 Data_* right=static_cast<Data_*>(r);
3485 const Ty s = (*right)[0];
3486 Data_* res = NewResult();
3487 TRACEOMP( __FILE__, __LINE__)
3488 #pragma omp parallel if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl))
3489 {
3490 #pragma omp for
3491 for( OMPInt i=0; i<nEl; ++i)
3492 (*res)[i] = pow( (*this)[ i], s);
3493 }
3494 return res;
3495 }
3496 // double complex inverse power of value: left=right ^ left
3497 template<>
PowInvSNew(BaseGDL * r)3498 Data_<SpDComplexDbl>* Data_<SpDComplexDbl>::PowInvSNew( BaseGDL* r)
3499 {
3500 Data_* right=static_cast<Data_*>(r);
3501
3502 ULong nEl=N_Elements();
3503 assert( nEl);
3504 Ty s = (*right)[0];
3505 Data_* res = NewResult();
3506 TRACEOMP( __FILE__, __LINE__)
3507 #pragma omp parallel if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl))
3508 {
3509 #pragma omp for
3510 for( OMPInt i=0; i<nEl; ++i)
3511 (*res)[i] = pow( s, (*this)[ i]);
3512 }
3513 return res;
3514 }
3515 // invalid types
3516 template<>
PowSNew(BaseGDL * r)3517 Data_<SpDString>* Data_<SpDString>::PowSNew( BaseGDL* r)
3518 {
3519 throw GDLException("Cannot apply operation to datatype STRING.",true,false);
3520 return NULL;
3521 }
3522 template<>
PowInvSNew(BaseGDL * r)3523 Data_<SpDString>* Data_<SpDString>::PowInvSNew( BaseGDL* r)
3524 {
3525 throw GDLException("Cannot apply operation to datatype STRING.",true,false);
3526 return NULL;
3527 }
3528 template<>
PowSNew(BaseGDL * r)3529 Data_<SpDPtr>* Data_<SpDPtr>::PowSNew( BaseGDL* r)
3530 {
3531 throw GDLException("Cannot apply operation to datatype PTR.",true,false);
3532 return NULL;
3533 }
3534 template<>
PowInvSNew(BaseGDL * r)3535 Data_<SpDPtr>* Data_<SpDPtr>::PowInvSNew( BaseGDL* r)
3536 {
3537 throw GDLException("Cannot apply operation to datatype PTR.",true,false);
3538 return NULL;
3539 }
3540 template<>
PowSNew(BaseGDL * r)3541 Data_<SpDObj>* Data_<SpDObj>::PowSNew( BaseGDL* r)
3542 {
3543 throw GDLException("Cannot apply operation to datatype OBJECT.",true,false);
3544 return NULL;
3545 }
3546 template<>
PowInvSNew(BaseGDL * r)3547 Data_<SpDObj>* Data_<SpDObj>::PowInvSNew( BaseGDL* r)
3548 {
3549 throw GDLException("Cannot apply operation to datatype OBJECT.",true,false);
3550 return NULL;
3551 }
3552
3553 #include "instantiate_templates.hpp"
3554