1 /***************************************************************************
2 datatypes.cpp - GDL datatypes
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 #define INCLUDE_DATATYPESREF_CPP 1
33 #include "datatypesref.cpp"
34
35 #if defined(USE_PYTHON) || defined(PYTHON_MODULE)
36
37 # define INCLUDE_TOPYTHON_CPP 1
38 # include "topython.cpp"
39
40 # define INCLUDE_GDLPYTHON_CPP 1
41 # include "gdlpython.cpp"
42
43 # ifdef PYTHON_MODULE
44 # define INCLUDE_PYTHONGDL_CPP 1
45 # include "pythongdl.cpp"
46 # endif
47 #endif
48
49 #ifdef TESTTG
50
51 #include "test_template_grouping.cpp"
52 template<class Sp>
TestTemplateGrouping()53 void Data_<Sp>::TestTemplateGrouping()
54 {
55 // Ty ty = Test1();
56 bool b = Test2();
57 }
58
59 #endif
60
61 template<class Sp>
62 FreeListT Data_<Sp>::freeList;
63
64 #ifdef GDLARRAY_CACHE
65
66 #ifdef GDLARRAY_DEBUG
67
TraceCache(SizeT & cacheSize,SizeT sz,bool cacheIsNull,SizeT smallArraySize)68 inline void TraceCache( SizeT& cacheSize, SizeT sz, bool cacheIsNull, SizeT smallArraySize)
69 {
70 // if( cacheSize > smallArraySize && cacheSize == sz && !cacheIsNull)
71 // std::cout << "+++ CACHE HIT\tID: (" << &cacheSize << ") sz: " << cacheSize << std::endl;
72 // else
73 if( sz > smallArraySize)
74 std::cout << "+ New\t\tID: (" << &cacheSize << ") sz: " << sz << " cache size: " << cacheSize <<std::endl;
75 }
76
77 #else
78
79 #define TraceCache( a, b, c, d)
80
81 #endif
82
83 template<class Sp>
84 SizeT GDLArray<Sp>::cacheSize = 0;
85
86 template<>
87 SizeT GDLArray<char>::cacheSize = 0;
88
89 template<class Sp>
90 typename GDLArray<Sp>::Ty* GDLArray<Sp>::cache = NULL;
91
92 template<class Sp>
Cached(SizeT newSize)93 typename GDLArray<Sp>::Ty* GDLArray<Sp>::Cached( SizeT newSize)
94 {
95 assert( newSize > smallArraySize);
96 if( cache != NULL && cacheSize == newSize)
97 {
98 #ifdef GDLARRAY_DEBUG
99 std::cout << "*** CACHE HIT\tID: (" << &cacheSize << ") sz: " << cacheSize << " ***" << std::endl;
100 #endif
101 Ty* hit = cache;
102 cache = NULL;
103 return hit;
104 }
105 return new Ty[ newSize];
106 }
107
GDLArray(const GDLArray<Sp> & cp)108 template<class Sp>GDLArray<Sp>::GDLArray( const GDLArray<Sp>& cp) : sz( cp.size())
109 {
110 TraceCache( cacheSize, sz, cache==NULL, smallArraySize);
111
112 try {
113 buf = (sz > smallArraySize) ? Cached( sz) /* new Ty[ cp.size()]*/ : scalar;
114 } catch (std::bad_alloc&) { ThrowGDLException("Array requires more memory than available"); }
115
116 std::memcpy(buf,cp.buf,sz*sizeof(Ty));
117 }
118
GDLArray(SizeT s,bool b)119 template<class Sp>GDLArray<Sp>:: GDLArray( SizeT s, bool b) : sz( s)
120 {
121 TraceCache( cacheSize, sz, cache==NULL, smallArraySize);
122
123 try {
124 buf = (sz > smallArraySize) ? Cached( sz) /* new Ty[ sz]*/ : scalar;
125 } catch (std::bad_alloc&) { ThrowGDLException("Array requires more memory than available"); }
126 }
127
GDLArray(Ty val,SizeT s)128 template<class Sp>GDLArray<Sp>:: GDLArray( Ty val, SizeT s) : sz( s)
129 {
130 TraceCache( cacheSize, sz, cache==NULL, smallArraySize);
131
132 try {
133 buf = (sz > smallArraySize) ? Cached( sz) /* new Ty[ sz]*/ : scalar;
134 } catch (std::bad_alloc&) { ThrowGDLException("Array requires more memory than available"); }
135 for( SizeT i=0; i<sz; ++i)
136 buf[ i] = val;
137 }
138
GDLArray(const Ty * arr,SizeT s)139 template<class Sp>GDLArray<Sp>:: GDLArray( const Ty* arr, SizeT s) : sz( s)
140 {
141 TraceCache( cacheSize, sz, cache==NULL, smallArraySize);
142
143 try {
144 buf = (sz > smallArraySize) ? Cached( sz) /* new Ty[ sz]*/ : scalar;
145 } catch (std::bad_alloc&) { ThrowGDLException("Array requires more memory than available"); }
146 std::memcpy(buf,arr,sz*sizeof(Ty));
147 }
148
~GDLArray()149 template<class Sp>GDLArray<Sp>::~GDLArray() throw()
150 {
151 #ifdef GDLARRAY_DEBUG
152 if( buf == cache)
153 std::cout << "~~~ recycled cache\tID: (" << &cacheSize << ") sz: " << sz << "\tcacheSize: " << cacheSize << std::endl;
154 #endif
155
156 assert( buf != cache || sz == cacheSize);
157
158 if ( buf != NULL && buf != scalar &&
159 buf != cache // note: assumes cacheSize never changes for a given cache
160 )
161 {
162 assert( sz > smallArraySize);
163 if ( sz <= maxCache )
164 {
165
166 #ifdef GDLARRAY_DEBUG
167 std::cout << "--- free cache\tID: (" << &cacheSize << ") sz: " << cacheSize << "\tnew: " << sz << std::endl;
168 #endif
169 delete cache;
170 cache = buf;
171 cacheSize = sz;
172 }
173 else
174 {
175 delete[] buf;
176 }
177 }
178 }
179 // as strings may occupy arbitrary memory (regardless of the array size), better not cache them
180 // note: Structs are ok, since GDL cleans up the strings they may contain and uses GDLArray as raw memory
181 template<>
~GDLArray()182 GDLArray<DString>::~GDLArray() throw()
183 {
184 if ( buf != scalar )
185 {
186 delete[] buf;
187 }
188 }
189
190
191 template class GDLArray<char>;
192
193 #endif
194
195 template<typename T>
gdlValid(const T & value)196 inline bool gdlValid( const T &value )
197 {
198 T max_value = std::numeric_limits<T>::max();
199 T min_value = - max_value;
200 return ( ( min_value <= value && value <= max_value ) && (value == value));
201 }
gdlValid(const DComplex & value)202 inline bool gdlValid( const DComplex &value )
203 {
204 DFloat max_value = std::numeric_limits<DFloat>::max();
205 DFloat min_value = - max_value;
206 return ( ( min_value <= value.real() && value.real() <= max_value ) && (value.real() == value.real()))&&
207 ( ( min_value <= value.imag() && value.imag() <= max_value ) && (value.imag() == value.imag()));
208 }
gdlValid(const DComplexDbl & value)209 inline bool gdlValid( const DComplexDbl &value )
210 {
211 DDouble max_value = std::numeric_limits<DDouble>::max();
212 DDouble min_value = - max_value;
213 return ( ( min_value <= value.real() && value.real() <= max_value ) && (value.real() == value.real()))&&
214 ( ( min_value <= value.imag() && value.imag() <= max_value ) && (value.imag() == value.imag()));
215 }
216
217
218
operator new(size_t bytes)219 template<class Sp> void* Data_<Sp>::operator new( size_t bytes)
220 {
221 assert( bytes == sizeof( Data_));
222
223 if( freeList.size() > 0)
224 {
225 return freeList.pop_back();
226 // void* res = freeList.back();
227 // freeList.pop_back();
228 // return res;
229 }
230
231 const size_t newSize = multiAlloc - 1;
232
233 static long callCount = 0;
234 ++callCount;
235
236 // reserve space for all instances
237 // note that reserve must do an allocation
238 // this hack divides the number of those allocation
239 // (for the cost of initially larger allocation - but only for pointers)
240 const long allocDivider = 4;
241 freeList.reserve( ((callCount/allocDivider+1)*allocDivider-1)*multiAlloc);
242
243 // resize to what is needed now
244 // freeList.resize( newSize);
245
246 #ifdef USE_EIGEN
247 // we need this allocation here as well (as in typedefs.hpp), because GDLArray needs to be aligned
248 const int alignmentInBytes = 16; // set to multiple of 16 >= sizeof( char*)
249 const size_t realSizeOfType = sizeof( Data_);
250 const SizeT exceed = realSizeOfType % alignmentInBytes;
251 const size_t sizeOfType = realSizeOfType + (alignmentInBytes - exceed);
252 char* res = static_cast< char*>( Eigen::internal::aligned_malloc( sizeOfType * multiAlloc)); // one more than newSize
253 #else
254 const size_t sizeOfType = sizeof( Data_);
255 char* res = static_cast< char*>( malloc( sizeOfType * multiAlloc)); // one more than newSize
256 #endif
257
258 res = freeList.Init( newSize, res, sizeOfType);
259 // freeList[0] = NULL;
260 // for( size_t i=1; i<=newSize; ++i)
261 // {
262 // freeList[ i] = res;
263 // res += sizeOfType;
264 // }
265
266 // the one more
267 return res;
268 }
269
operator delete(void * ptr)270 template<class Sp> void Data_<Sp>::operator delete( void *ptr)
271 {
272 freeList.push_back( ptr);
273 }
274
275
276
277 // destructor
~Data_()278 template<class Sp> Data_<Sp>::~Data_() {}
~Data_()279 template<> Data_<SpDPtr>::~Data_()
280 {
281 if( this->dd.GetBuffer() != NULL)
282 GDLInterpreter::DecRef( this);
283 }
~Data_()284 template<> Data_<SpDObj>::~Data_()
285 {
286 if( this->dd.GetBuffer() != NULL)
287 GDLInterpreter::DecRefObj( this);
288 }
289
290 // default
Data_()291 template<class Sp> Data_<Sp>::Data_(): Sp(), dd() {}
292
293 // scalar
Data_(const Ty & d_)294 template<class Sp> Data_<Sp>::Data_(const Ty& d_): Sp(), dd(d_)
295 {}
296 // template<> Data_<SpDPtr>::Data_(const Ty& d_): SpDPtr(), dd(d_)
297 // {GDLInterpreter::IncRef(d_);}
298 // template<> Data_<SpDObj>::Data_(const Ty& d_): SpDObj(), dd(d_)
299 // {GDLInterpreter::IncRefObj(d_);}
300
301 // new array, zero fields
Data_(const dimension & dim_)302 template<class Sp> Data_<Sp>::Data_(const dimension& dim_):
303 Sp( dim_), dd( Sp::zero, this->dim.NDimElements())
304 {
305 this->dim.Purge();
306 }
307
308 // new one-dim array from Ty*
Data_(const Ty * p,const SizeT nEl)309 template<class Sp> Data_<Sp>::Data_( const Ty* p, const SizeT nEl):
310 Sp( dimension( nEl)), dd( p, nEl)
311 {}
Data_(const Ty * p,const SizeT nEl)312 template<> Data_<SpDPtr>::Data_( const Ty* p, const SizeT nEl):
313 SpDPtr( dimension( nEl)), dd( p, nEl)
314 {GDLInterpreter::IncRef(this);}
Data_(const Ty * p,const SizeT nEl)315 template<> Data_<SpDObj>::Data_( const Ty* p, const SizeT nEl):
316 SpDObj( dimension( nEl)), dd( p, nEl)
317 {GDLInterpreter::IncRefObj(this);}
318
319 // c-i
320 // template<class Sp> Data_<Sp>::Data_(const Data_& d_):
321 // Sp(d_.dim), dd(d_.dd) {}
322
Data_(const dimension & dim_,BaseGDL::InitType iT,DDouble off,DDouble inc)323 template<class Sp> Data_<Sp>::Data_(const dimension& dim_, BaseGDL::InitType iT, DDouble off, DDouble inc):
324 Sp( dim_), dd( (iT == BaseGDL::NOALLOC) ? 0 : this->dim.NDimElements(), false)
325 {
326 this->dim.Purge();
327 if (iT == BaseGDL::NOZERO) return; //very frequent
328 else if (iT == BaseGDL::ZERO) { //rather frequent
329 SizeT sz = dd.size();
330 #pragma omp parallel if (sz >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= sz))
331 {
332 #pragma omp for
333 for (SizeT i = 0; i < sz; ++i) {
334 (*this)[i] = 0;
335 }
336 }
337 }
338 else if (iT == BaseGDL::INDGEN) { //less frequent
339 SizeT sz = dd.size();
340 if (off==0 && inc==1) {
341 #pragma omp parallel if (sz >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= sz))
342 {
343 #pragma omp for
344 for (SizeT i = 0; i < sz; ++i) {
345 (*this)[i]=i;
346 }
347 }
348 } else {
349 #pragma omp parallel if (sz >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= sz))
350 {
351 #pragma omp for
352 for (SizeT i = 0; i < sz; ++i) {
353 (*this)[i]=off+i*inc;
354 }
355 }
356 }
357 }
358 }
359 //IDL uses floats increments and offset for floats and complex (normal) .
Data_(const dimension & dim_,BaseGDL::InitType iT,DDouble off,DDouble inc)360 template<> Data_<SpDFloat>::Data_(const dimension& dim_, BaseGDL::InitType iT, DDouble off, DDouble inc):
361 SpDFloat( dim_), dd( (iT == BaseGDL::NOALLOC) ? 0 : this->dim.NDimElements(), false)
362 {
363 this->dim.Purge();
364 if (iT == BaseGDL::NOZERO) return; //very frequent
365 else if (iT == BaseGDL::ZERO) { //rather frequent
366 SizeT sz = dd.size();
367 #pragma omp parallel if (sz >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= sz))
368 {
369 #pragma omp for
370 for (SizeT i = 0; i < sz; ++i) {
371 (*this)[i] = 0;
372 }
373 }
374 }
375 else if (iT == BaseGDL::INDGEN) { //less frequent
376 SizeT sz = dd.size();
377 if (off==0 && inc==1) {
378 #pragma omp parallel if (sz >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= sz))
379 {
380 #pragma omp for
381 for (SizeT i = 0; i < sz; ++i) {
382 (*this)[i]=i;
383 }
384 }
385 } else {
386 DFloat f_off=off;
387 DFloat f_inc=inc;
388 #pragma omp parallel if (sz >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= sz))
389 {
390 #pragma omp for
391 for (SizeT i = 0; i < sz; ++i) {
392 (*this)[i]=f_off+i*f_inc;
393 }
394 }
395 }
396 }
397 }
398
Data_(const dimension & dim_,BaseGDL::InitType iT,DDouble off,DDouble inc)399 template<> Data_<SpDComplex>::Data_(const dimension& dim_, BaseGDL::InitType iT, DDouble off, DDouble inc):
400 SpDComplex( dim_), dd( (iT == BaseGDL::NOALLOC) ? 0 : this->dim.NDimElements(), false)
401 {
402 this->dim.Purge();
403 if (iT == BaseGDL::NOZERO) return; //very frequent
404 else if (iT == BaseGDL::ZERO) { //rather frequent
405 SizeT sz = dd.size();
406 #pragma omp parallel if (sz >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= sz))
407 {
408 #pragma omp for
409 for (SizeT i = 0; i < sz; ++i) {
410 (*this)[i] = 0;
411 }
412 }
413 }
414 else if (iT == BaseGDL::INDGEN) { //less frequent
415 SizeT sz = dd.size();
416 if (off==0 && inc==1) {
417 #pragma omp parallel if (sz >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= sz))
418 {
419 #pragma omp for
420 for (SizeT i = 0; i < sz; ++i) {
421 (*this)[i]=i;
422 }
423 }
424 } else {
425 DFloat f_off=off;
426 DFloat f_inc=inc;
427 #pragma omp parallel if (sz >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= sz))
428 {
429 #pragma omp for
430 for (SizeT i = 0; i < sz; ++i) {
431 (*this)[i]=f_off+i*f_inc;
432 }
433 }
434 }
435 }
436 }
437
Data_(const dimension & dim_,BaseGDL::InitType iT,DDouble,DDouble)438 template<> Data_<SpDString>::Data_(const dimension& dim_, BaseGDL::InitType iT, DDouble, DDouble):
439 SpDString(dim_), dd( (iT == BaseGDL::NOALLOC) ? 0 : this->dim.NDimElements(), false)
440 {
441 dim.Purge();
442
443 if( iT == BaseGDL::INDGEN)
444 throw GDLException("DStringGDL(dim,InitType=INDGEN) called.");
445 }
Data_(const dimension & dim_,BaseGDL::InitType iT,DDouble,DDouble)446 template<> Data_<SpDPtr>::Data_(const dimension& dim_, BaseGDL::InitType iT, DDouble, DDouble):
447 SpDPtr(dim_), dd( (iT == BaseGDL::NOALLOC) ? 0 : this->dim.NDimElements(), false)
448 {
449 dim.Purge();
450
451 if( iT == BaseGDL::INDGEN)
452 throw GDLException("DPtrGDL(dim,InitType=INDGEN) called.");
453
454 if( iT != BaseGDL::NOALLOC && iT != BaseGDL::NOZERO)
455 {
456 SizeT sz = dd.size();
457 #pragma omp parallel if (sz >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= sz))
458 {
459 #pragma omp for
460 for( SizeT i=0; i<sz; ++i)
461 {
462 (*this)[i]=0;
463 }
464 // val += 1; // no increment operator for floats
465 }
466 }
467 }
Data_(const dimension & dim_,BaseGDL::InitType iT,DDouble,DDouble)468 template<> Data_<SpDObj>::Data_(const dimension& dim_, BaseGDL::InitType iT, DDouble, DDouble):
469 SpDObj(dim_), dd( (iT == BaseGDL::NOALLOC) ? 0 : this->dim.NDimElements(), false)
470 {
471 dim.Purge();
472
473 if( iT == BaseGDL::INDGEN)
474 throw GDLException("DObjGDL(dim,InitType=INDGEN) called.");
475
476 if( iT != BaseGDL::NOALLOC && iT != BaseGDL::NOZERO)
477 {
478 SizeT sz = dd.size();
479 #pragma omp parallel if (sz >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= sz))
480 {
481 #pragma omp for
482 for( SizeT i=0; i<sz; i++)
483 {
484 (*this)[i]=0;
485 }
486 // val += 1; // no increment operator for floats
487 }
488 }
489 }
490
491 // c-i
492 //template<class Sp>
493 //Data_<Sp>::Data_(const Data_& d_): Sp(d_.dim), dd(d_.dd) { }
494
495 //faster, C++ initializer is probably too shy.
496 template<class Sp>
Data_(const Data_ & d_)497 Data_<Sp>::Data_(const Data_& d_) : Sp(d_.dim), dd(this->dim.NDimElements(), false) {
498 this->dim.Purge(); //useful?
499 SizeT sz = dd.size();
500 #pragma omp parallel if (CpuTPOOL_NTHREADS > 1 && sz >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= sz))
501 {
502 #pragma omp for
503 for (SizeT i = 0; i < sz; i++) {
504 dd[i] = d_.dd[i];
505 }
506 }
507 }
508
509 template<>
Data_(const Data_ & d_)510 Data_<SpDPtr>::Data_(const Data_& d_): SpDPtr(d_.dim), dd(d_.dd)
511 {
512 GDLInterpreter::IncRef( this);
513 }
514 template<>
Data_(const Data_ & d_)515 Data_<SpDObj>::Data_(const Data_& d_): SpDObj(d_.dim), dd(d_.dd)
516 {
517 GDLInterpreter::IncRefObj( this);
518 }
519
520
521 template<class Sp>
Dup() const522 Data_<Sp>* Data_<Sp>::Dup() const { return new Data_(*this);}
523
524 // template<>
525 // Data_<SpDPtr>* Data_<SpDPtr>::Dup() const
526 // {
527 // Data_<SpDPtr>* p =new Data_(*this);
528 // GDLInterpreter::IncRef( p);
529 // return p;
530 // }
531 // template<>
532 // Data_<SpDObj>* Data_<SpDObj>::Dup() const
533 // {
534 // Data_<SpDObj>* p =new Data_(*this);
535 // GDLInterpreter::IncRefObj( p);
536 // return p;
537 // }
538
539
540
541 template<class Sp>
Log()542 BaseGDL* Data_<Sp>::Log()
543 {
544 DFloatGDL* res = static_cast<DFloatGDL*>
545 (this->Convert2( GDL_FLOAT, BaseGDL::COPY));
546 res->LogThis();
547 return res;
548 }
549 template<>
Log()550 BaseGDL* Data_<SpDFloat>::Log()
551 {
552 Data_* n = this->New( this->dim, BaseGDL::NOZERO);
553 SizeT nEl = n->N_Elements();
554 if( nEl == 1)
555 {
556 (*n)[ 0] = log( (*this)[ 0]);
557 return n;
558 }
559 #pragma omp parallel for if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl))
560 for( SizeT i=0; i<nEl; ++i) (*n)[ i] = log( (*this)[ i]);
561 return n;
562 }
563 template<>
Log()564 BaseGDL* Data_<SpDDouble>::Log()
565 {
566 Data_* n = this->New( this->dim, BaseGDL::NOZERO);
567 SizeT nEl = n->N_Elements();
568 if( nEl == 1)
569 {
570 (*n)[ 0] = log( (*this)[ 0]);
571 return n;
572 }
573 #pragma omp parallel for if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl))
574 for( SizeT i=0; i<nEl; ++i) (*n)[ i] = log( (*this)[ i]);
575 return n;
576 }
577 template<>
Log()578 BaseGDL* Data_<SpDComplex>::Log()
579 {
580 Data_* n = this->New( this->dim, BaseGDL::NOZERO);
581 SizeT nEl = n->N_Elements();
582 if( nEl == 1)
583 {
584 (*n)[ 0] = log( (*this)[ 0]);
585 return n;
586 }
587 #pragma omp parallel for if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl))
588 for( SizeT i=0; i<nEl; ++i) (*n)[ i] = log( (*this)[ i]);
589 return n;
590 }
591 template<>
Log()592 BaseGDL* Data_<SpDComplexDbl>::Log()
593 {
594 Data_* n = this->New( this->dim, BaseGDL::NOZERO);
595 SizeT nEl = n->N_Elements();
596 if( nEl == 1)
597 {
598 (*n)[ 0] = log( (*this)[ 0]);
599 return n;
600 }
601 #pragma omp parallel for if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl))
602 for( SizeT i=0; i<nEl; ++i) (*n)[ i] = log( (*this)[ i]);
603 return n;
604 }
605
606 // this is actually not a "log" of "this",
607 // but the behaviour is fine with the usage in the library function
608 // the real LogThis is done in the specializations for floating and
609 // complex types
610 template<class Sp>
LogThis()611 BaseGDL* Data_<Sp>::LogThis()
612 {
613 DFloatGDL* res = static_cast<DFloatGDL*>
614 (this->Convert2( GDL_FLOAT, BaseGDL::COPY));
615 res->LogThis(); // calls correct LogThis for float
616 return res;
617 }
618 template<>
LogThis()619 BaseGDL* Data_<SpDFloat>::LogThis()
620 {
621 SizeT nEl = N_Elements();
622 if( nEl == 1)
623 {
624 (*this)[ 0] = log( (*this)[ 0]);
625 return this;
626 }
627 #pragma omp parallel for if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl))
628 for( SizeT i=0; i<nEl; ++i) (*this)[ i] = log( (*this)[ i]);
629 return this;
630 }
631 template<>
LogThis()632 BaseGDL* Data_<SpDDouble>::LogThis()
633 {
634 SizeT nEl = N_Elements();
635 if( nEl == 1)
636 {
637 (*this)[ 0] = log( (*this)[ 0]);
638 return this;
639 }
640 #pragma omp parallel for if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl))
641 for( SizeT i=0; i<nEl; ++i) (*this)[ i] = log( (*this)[ i]);
642 return this;
643 }
644 template<>
LogThis()645 BaseGDL* Data_<SpDComplex>::LogThis()
646 {
647 SizeT nEl = N_Elements();
648 if( nEl == 1)
649 {
650 (*this)[ 0] = log( (*this)[ 0]);
651 return this;
652 }
653 #pragma omp parallel for if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl))
654 for( SizeT i=0; i<nEl; ++i) (*this)[ i] = log( (*this)[ i]);
655 return this;
656 }
657 template<>
LogThis()658 BaseGDL* Data_<SpDComplexDbl>::LogThis()
659 {
660 SizeT nEl = N_Elements();
661 if( nEl == 1)
662 {
663 (*this)[ 0] = log( (*this)[ 0]);
664 return this;
665 }
666 #pragma omp parallel for if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl))
667 for( SizeT i=0; i<nEl; ++i) (*this)[ i] = log( (*this)[ i]);
668 return this;
669 }
670
671 template<class Sp>
Log10()672 BaseGDL* Data_<Sp>::Log10()
673 {
674 DFloatGDL* res = static_cast<DFloatGDL*>
675 (this->Convert2( GDL_FLOAT, BaseGDL::COPY));
676 res->Log10This();
677 return res;
678 }
679 template<>
Log10()680 BaseGDL* Data_<SpDFloat>::Log10()
681 {
682 Data_* n = this->New( this->dim, BaseGDL::NOZERO);
683 SizeT nEl = n->N_Elements();
684 if( nEl == 1)
685 {
686 (*n)[ 0] = log10( (*this)[ 0]);
687 return n;
688 }
689 #pragma omp parallel for if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl))
690 for( SizeT i=0; i<nEl; ++i) (*n)[ i] = log10( (*this)[ i]);
691 return n;
692 }
693 template<>
Log10()694 BaseGDL* Data_<SpDDouble>::Log10()
695 {
696 Data_* n = this->New( this->dim, BaseGDL::NOZERO);
697 SizeT nEl = n->N_Elements();
698 if( nEl == 1)
699 {
700 (*n)[ 0] = log10( (*this)[ 0]);
701 return n;
702 }
703 #pragma omp parallel for if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl))
704 for( SizeT i=0; i<nEl; ++i) (*n)[ i] = log10( (*this)[ i]);
705 return n;
706 }
707 template<>
Log10()708 BaseGDL* Data_<SpDComplex>::Log10()
709 {
710 Data_* n = this->New( this->dim, BaseGDL::NOZERO);
711 SizeT nEl = n->N_Elements();
712 if( nEl == 1)
713 {
714 (*n)[ 0] = log10( (*this)[ 0]);
715 return n;
716 }
717 #pragma omp parallel for if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl))
718 for( SizeT i=0; i<nEl; ++i) (*n)[ i] = log10( (*this)[ i]);
719 return n;
720 }
721 template<>
Log10()722 BaseGDL* Data_<SpDComplexDbl>::Log10()
723 {
724 Data_* n = this->New( this->dim, BaseGDL::NOZERO);
725 SizeT nEl = n->N_Elements();
726 if( nEl == 1)
727 {
728 (*n)[ 0] = log10( (*this)[ 0]);
729 return n;
730 }
731 #pragma omp parallel for if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl))
732 for( SizeT i=0; i<nEl; ++i) (*n)[ i] = log10( (*this)[ i]);
733 return n;
734 }
735
736 // see comment at void Data_<Sp>::LogThis()
737 template<class Sp>
Log10This()738 BaseGDL* Data_<Sp>::Log10This()
739 {
740 DFloatGDL* res = static_cast<DFloatGDL*>
741 (this->Convert2( GDL_FLOAT, BaseGDL::COPY));
742 res->Log10This(); // calls correct Log10This for float
743 return res;
744 }
745 template<>
Log10This()746 BaseGDL* Data_<SpDFloat>::Log10This()
747 {
748 #if 1 || (__GNUC__ == 3) && (__GNUC_MINOR__ == 2) //&& (__GNUC_PATCHLEVEL__ == 2)
749
750 SizeT nEl = N_Elements();
751 if( nEl == 1)
752 {
753 (*this)[ 0] = log10( (*this)[ 0]);
754 return this;
755 }
756 TRACEOMP( __FILE__, __LINE__)
757 #pragma omp parallel for if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl))
758 for( SizeT i=0; i<nEl; ++i)
759 (*this)[ i] = log10( (*this)[ i]);
760 #else
761 dd = log10(dd);
762 #endif
763 return this;
764 }
765 template<>
Log10This()766 BaseGDL* Data_<SpDDouble>::Log10This()
767 {
768 #if 1 || (__GNUC__ == 3) && (__GNUC_MINOR__ == 2) //&& (__GNUC_PATCHLEVEL__ == 2)
769
770 SizeT nEl = N_Elements();
771 if( nEl == 1)
772 {
773 (*this)[ 0] = log10( (*this)[ 0]);
774 return this;
775 }
776 TRACEOMP( __FILE__, __LINE__)
777 #pragma omp parallel for if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl))
778 for( SizeT i=0; i<nEl; ++i)
779 (*this)[ i] = log10( (*this)[ i]);
780 #else
781 dd = log10(dd);
782 #endif
783 return this;
784 }
785 template<>
Log10This()786 BaseGDL* Data_<SpDComplex>::Log10This()
787 {
788 #if 1 || (__GNUC__ == 3) && (__GNUC_MINOR__ == 2) //&& (__GNUC_PATCHLEVEL__ == 2)
789
790 SizeT nEl = N_Elements();
791 if( nEl == 1)
792 {
793 (*this)[ 0] = log10( (*this)[ 0]);
794 return this;
795 }
796 TRACEOMP( __FILE__, __LINE__)
797 #pragma omp parallel for if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl))
798 for( SizeT i=0; i<nEl; ++i)
799 (*this)[ i] = log10( (*this)[ i]);
800 #else
801 dd = log10(dd);
802 #endif
803 return this;
804 }
805 template<>
Log10This()806 BaseGDL* Data_<SpDComplexDbl>::Log10This()
807 {
808 #if 1 || (__GNUC__ == 3) && (__GNUC_MINOR__ == 2) //&& (__GNUC_PATCHLEVEL__ == 2)
809
810 SizeT nEl = N_Elements();
811 if( nEl == 1)
812 {
813 (*this)[ 0] = log10( (*this)[ 0]);
814 return this;
815 }
816 TRACEOMP( __FILE__, __LINE__)
817 #pragma omp parallel for if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl))
818 for( SizeT i=0; i<nEl; ++i)
819 (*this)[ i] = log10( (*this)[ i]);
820 #else
821 dd = log10(dd);
822 #endif
823 return this;
824 }
825
826
827
828 // template<class Sp>
829 // BaseGDL* Data_<Sp>::Abs() const
830 // {
831 // return new Data_( this->dim, dd.abs());
832 // }
833
834 template<class Sp>
Greater(SizeT i1,SizeT i2) const835 inline bool Data_<Sp>::Greater(SizeT i1, SizeT i2) const
836 { return ((*this)[i1] > (*this)[i2]);}
837
838 template<>
Greater(SizeT i1,SizeT i2) const839 inline bool Data_<SpDComplex>::Greater(SizeT i1, SizeT i2) const
840 { return (abs((*this)[i1]) > abs((*this)[i2]));}
841 template<>
Greater(SizeT i1,SizeT i2) const842 inline bool Data_<SpDComplexDbl>::Greater(SizeT i1, SizeT i2) const
843 { return (abs((*this)[i1]) > abs((*this)[i2]));}
844
845
846 template<class Sp>
Equal(SizeT i1,SizeT i2) const847 inline bool Data_<Sp>::Equal(SizeT i1, SizeT i2) const
848 { return ((*this)[i1] == (*this)[i2]);}
849
850
851
852 // calculates the shift to be applied to the destination index
CShiftNormalize(DLong s,SizeT this_dim)853 inline SizeT CShiftNormalize( DLong s, SizeT this_dim)
854 {
855 if ( s >= 0 )
856 return s % this_dim;
857 // s < 0
858 // long dstIx = -(-s % this_dim);
859 long dstIx = -s % this_dim;
860 dstIx = -dstIx;
861 if( dstIx == 0) // if this_dim == 1
862 return 0;
863 assert( dstIx + this_dim > 0);
864 return dstIx + this_dim;
865 }
866
867 template<class Sp>
CShift(DLong d) const868 BaseGDL* Data_<Sp>::CShift( DLong d) const
869 {
870 SizeT nEl = dd.size();
871 SizeT shift = CShiftNormalize( d, nEl);
872
873 if( shift == 0)
874 return this->Dup();
875
876 Data_* sh = new Data_( this->dim, BaseGDL::NOZERO);
877
878 SizeT firstChunk = nEl - shift;
879
880 memcpy( &sh->dd[ shift], &dd[0], firstChunk * sizeof(Ty));
881 memcpy( &sh->dd[ 0], &dd[firstChunk], shift * sizeof(Ty));
882
883 return sh;
884 }
885
886 template<>
CShift(DLong d) const887 BaseGDL* Data_<SpDString>::CShift( DLong d) const
888 {
889 SizeT nEl = dd.size();
890 SizeT shift = CShiftNormalize( d, nEl);
891
892 if( shift == 0)
893 return this->Dup();
894
895 Data_* sh = new Data_( this->dim, BaseGDL::NOZERO);
896
897 SizeT firstChunk = nEl - shift;
898
899 SizeT i=0;
900 for( ; i<firstChunk; ++i)
901 sh->dd[shift++] = dd[ i];
902
903 shift = 0;
904
905 for( ; i<nEl; ++i)
906 sh->dd[shift++] = dd[ i];
907
908 return sh;
909 }
910 template<>
CShift(DLong d) const911 BaseGDL* Data_<SpDPtr>::CShift( DLong d) const
912 {
913 SizeT nEl = dd.size();
914 SizeT shift = CShiftNormalize( d, nEl);
915
916 if( shift == 0)
917 return this->Dup(); // does IncRef
918
919 Data_* sh = new Data_( this->dim, BaseGDL::NOZERO);
920
921 SizeT firstChunk = nEl - shift;
922
923 SizeT i=0;
924 for( ; i<firstChunk; ++i)
925 sh->dd[shift++] = dd[ i];
926
927 shift = 0;
928
929 for( ; i<nEl; ++i)
930 sh->dd[shift++] = dd[ i];
931
932 GDLInterpreter::IncRef( sh);
933 return sh;
934 }
935 template<>
CShift(DLong d) const936 BaseGDL* Data_<SpDObj>::CShift( DLong d) const
937 {
938 SizeT nEl = dd.size();
939 SizeT shift = CShiftNormalize( d, nEl);
940
941 if( shift == 0)
942 return this->Dup(); // does IncRefObj
943
944 Data_* sh = new Data_( this->dim, BaseGDL::NOZERO);
945
946 SizeT firstChunk = nEl - shift;
947
948 SizeT i=0;
949 for( ; i<firstChunk; ++i)
950 sh->dd[shift++] = dd[ i];
951
952 shift = 0;
953
954 for( ; i<nEl; ++i)
955 sh->dd[shift++] = dd[ i];
956
957 GDLInterpreter::IncRefObj( sh);
958 return sh;
959 }
960
961 template<typename Ty>
CShift1(Ty * dst,SizeT & dstLonIx,const Ty * src,SizeT & srcLonIx,SizeT stride_1,SizeT chunk0,SizeT chunk1)962 inline void CShift1( Ty* dst, SizeT& dstLonIx, const Ty* src, SizeT& srcLonIx,
963 SizeT stride_1, SizeT chunk0, SizeT chunk1)
964 {
965 memcpy( &dst[ dstLonIx], &src[ srcLonIx], chunk0 * sizeof(Ty));
966 dstLonIx += chunk0;
967 srcLonIx += chunk0;
968
969 dstLonIx -= stride_1;
970
971 memcpy( &dst[ dstLonIx], &src[ srcLonIx], chunk1 * sizeof(Ty));
972 dstLonIx += chunk1 ;
973 srcLonIx += chunk1;
974
975 dstLonIx += stride_1;
976 }
977
978 #undef TEST_GOOD_OL_VERSION
979
980 template<class Sp>
CShift(DLong s[MAXRANK]) const981 BaseGDL* Data_<Sp>::CShift( DLong s[ MAXRANK]) const {
982 Data_* sh = new Data_(this->dim, BaseGDL::NOZERO);
983
984 SizeT nDim = this->Rank();
985 SizeT nEl = N_Elements();
986
987 SizeT stride[ MAXRANK + 1];
988 this->dim.Stride(stride, nDim);
989
990 long srcIx[ MAXRANK + 1];
991 long dstIx[ MAXRANK + 1];
992 SizeT this_dim[ MAXRANK];
993
994 const Ty* ddP = &(*this)[0];
995 Ty* shP = &(*sh)[0];
996
997 if (nDim == 2) {
998 this_dim[ 0] = this->dim[ 0];
999 this_dim[ 1] = this->dim[ 1];
1000 dstIx[ 0] = CShiftNormalize(s[ 0], this_dim[ 0]);
1001 dstIx[ 1] = CShiftNormalize(s[ 1], this_dim[ 1]);
1002 SizeT dstLonIx = dstIx[ 0] + dstIx[ 1] * stride[ 1];
1003 SizeT freeDstIx_0 = this_dim[ 0] - dstIx[ 0];
1004 SizeT freeDstIx_1 = this_dim[ 1] - dstIx[ 1];
1005
1006 if (Sp::t != GDL_STRING) // strings are not POD all others are
1007 {
1008 SizeT srcLonIx = 0;
1009 SizeT t = 0;
1010
1011 for (; t < freeDstIx_1; ++t) {
1012 CShift1<Ty>(shP, dstLonIx, ddP, srcLonIx, stride[ 1], freeDstIx_0, dstIx[0]);
1013 }
1014 dstLonIx -= stride[ 2];
1015 for (; t < this_dim[ 1]; ++t) {
1016 CShift1<Ty>(shP, dstLonIx, ddP, srcLonIx, stride[ 1], freeDstIx_0, dstIx[0]);
1017 }
1018 } else // Sp::t == GDL_STRING
1019 {
1020 SizeT a = 0;
1021 SizeT t = 0;
1022 for (; t < freeDstIx_1; ++t) {
1023 SizeT s = 0;
1024 for (; s < freeDstIx_0; ++s) {
1025 shP[ dstLonIx++] = ddP[ a++];
1026 }
1027 dstLonIx -= stride[ 1];
1028 for (; s < this_dim[ 0]; ++s) {
1029 shP[ dstLonIx++] = ddP[ a++];
1030 }
1031 dstLonIx += stride[ 1];
1032 }
1033 dstLonIx -= stride[ 2];
1034 for (; t < this_dim[ 1]; ++t) {
1035 SizeT s = 0;
1036 for (; s < freeDstIx_0; ++s) {
1037 shP[ dstLonIx++] = ddP[ a++];
1038 }
1039 dstLonIx -= stride[ 1];
1040 for (; s < this_dim[ 0]; ++s) {
1041 shP[ dstLonIx++] = ddP[ a++];
1042 }
1043 dstLonIx += stride[ 1];
1044 }
1045 // dstLonIx += stride[ 2];
1046 assert(a == nEl);
1047 } // if( Sp::t != GDL_STRING) else
1048
1049 return sh;
1050 }
1051
1052 assert(nDim > 2);
1053
1054 #ifndef TEST_GOOD_OL_VERSION
1055
1056 for (SizeT aSp = 0; aSp < nDim; ++aSp) {
1057 this_dim[ aSp] = this->dim[ aSp];
1058 srcIx[ aSp] = 0;
1059 dstIx[ aSp] = CShiftNormalize(s[ aSp], this_dim[ aSp]);
1060 // dim_stride[ aSp] = this_dim[ aSp] * stride[ aSp];
1061 }
1062 // srcIx[ nDim] = dstIx[ nDim] = 0;
1063 SizeT dstLonIx = dstIx[ 0];
1064 for (SizeT rSp = 1; rSp < nDim; ++rSp)
1065 dstLonIx += dstIx[ rSp] * stride[ rSp];
1066
1067 if (Sp::t != GDL_STRING) {
1068 if (nDim == 3) {
1069 SizeT freeDstIx_0 = this_dim[ 0] - dstIx[ 0];
1070 SizeT freeDstIx_1 = this_dim[ 1] - dstIx[ 1];
1071 SizeT freeDstIx_2 = this_dim[ 2] - dstIx[ 2];
1072
1073 SizeT srcLonIx = 0;
1074 SizeT d2 = 0;
1075 for (; d2 < freeDstIx_2; ++d2) {
1076 SizeT d1 = 0;
1077 for (; d1 < freeDstIx_1; ++d1) {
1078 CShift1<Ty>(shP, dstLonIx, ddP, srcLonIx, stride[ 1], freeDstIx_0, dstIx[0]);
1079 }
1080 dstLonIx -= stride[ 2];
1081 for (; d1 < this_dim[ 1]; ++d1) {
1082 CShift1<Ty>(shP, dstLonIx, ddP, srcLonIx, stride[ 1], freeDstIx_0, dstIx[0]);
1083 }
1084 dstLonIx += stride[ 2];
1085 }
1086 dstLonIx -= stride[ 3];
1087 for (; d2 < this_dim[ 2]; ++d2) {
1088 SizeT d1 = 0;
1089 for (; d1 < freeDstIx_1; ++d1) {
1090 CShift1<Ty>(shP, dstLonIx, ddP, srcLonIx, stride[ 1], freeDstIx_0, dstIx[0]);
1091 }
1092 dstLonIx -= stride[ 2];
1093 for (; d1 < this_dim[ 1]; ++d1) {
1094 CShift1<Ty>(shP, dstLonIx, ddP, srcLonIx, stride[ 1], freeDstIx_0, dstIx[0]);
1095 }
1096 dstLonIx += stride[ 2];
1097 }
1098 assert(srcLonIx == nEl);
1099 return sh;
1100 } // nDim == 3
1101 if (nDim == 4) {
1102 SizeT freeDstIx_0 = this_dim[ 0] - dstIx[ 0];
1103 SizeT freeDstIx_1 = this_dim[ 1] - dstIx[ 1];
1104 SizeT freeDstIx_2 = this_dim[ 2] - dstIx[ 2];
1105 SizeT freeDstIx_3 = this_dim[ 3] - dstIx[ 3];
1106
1107 SizeT srcLonIx = 0;
1108
1109 SizeT d3 = 0;
1110 for (; d3 < freeDstIx_3; ++d3) {
1111 SizeT d2 = 0;
1112 for (; d2 < freeDstIx_2; ++d2) {
1113 SizeT d1 = 0;
1114 for (; d1 < freeDstIx_1; ++d1)
1115 CShift1<Ty>(shP, dstLonIx, ddP, srcLonIx, stride[ 1], freeDstIx_0, dstIx[0]);
1116 dstLonIx -= stride[ 2];
1117 for (; d1 < this_dim[ 1]; ++d1)
1118 CShift1<Ty>(shP, dstLonIx, ddP, srcLonIx, stride[ 1], freeDstIx_0, dstIx[0]);
1119 dstLonIx += stride[ 2];
1120 }
1121 dstLonIx -= stride[ 3];
1122 for (; d2 < this_dim[ 2]; ++d2) {
1123 SizeT d1 = 0;
1124 for (; d1 < freeDstIx_1; ++d1)
1125 CShift1<Ty>(shP, dstLonIx, ddP, srcLonIx, stride[ 1], freeDstIx_0, dstIx[0]);
1126 dstLonIx -= stride[ 2];
1127 for (; d1 < this_dim[ 1]; ++d1)
1128 CShift1<Ty>(shP, dstLonIx, ddP, srcLonIx, stride[ 1], freeDstIx_0, dstIx[0]);
1129 dstLonIx += stride[ 2];
1130 }
1131 dstLonIx += stride[ 3];
1132 }
1133 dstLonIx -= stride[ 4];
1134 for (; d3 < this_dim[ 3]; ++d3) {
1135 SizeT d2 = 0;
1136 for (; d2 < freeDstIx_2; ++d2) {
1137 SizeT d1 = 0;
1138 for (; d1 < freeDstIx_1; ++d1)
1139 CShift1<Ty>(shP, dstLonIx, ddP, srcLonIx, stride[ 1], freeDstIx_0, dstIx[0]);
1140 dstLonIx -= stride[ 2];
1141 for (; d1 < this_dim[ 1]; ++d1)
1142 CShift1<Ty>(shP, dstLonIx, ddP, srcLonIx, stride[ 1], freeDstIx_0, dstIx[0]);
1143 dstLonIx += stride[ 2];
1144 }
1145 dstLonIx -= stride[ 3];
1146 for (; d2 < this_dim[ 2]; ++d2) {
1147 SizeT d1 = 0;
1148 for (; d1 < freeDstIx_1; ++d1)
1149 CShift1<Ty>(shP, dstLonIx, ddP, srcLonIx, stride[ 1], freeDstIx_0, dstIx[0]);
1150 dstLonIx -= stride[ 2];
1151 for (; d1 < this_dim[ 1]; ++d1)
1152 CShift1<Ty>(shP, dstLonIx, ddP, srcLonIx, stride[ 1], freeDstIx_0, dstIx[0]);
1153 dstLonIx += stride[ 2];
1154 }
1155 dstLonIx += stride[ 3];
1156 }
1157 assert(srcLonIx == nEl);
1158 return sh;
1159 } // nDim == 4
1160 } // if( Sp::t != GDL_STRING)
1161
1162 #else
1163
1164 // need to be done earlier within the TEST_GOOD_OL_VERSION section
1165 for (SizeT aSp = 0; aSp < nDim; ++aSp) {
1166 this_dim[ aSp] = this->dim[ aSp];
1167 srcIx[ aSp] = 0;
1168 dstIx[ aSp] = CShiftNormalize(s[ aSp], this_dim[ aSp]);
1169 }
1170 SizeT dstLonIx = dstIx[ 0];
1171 for (SizeT rSp = 1; rSp < nDim; ++rSp)
1172 dstLonIx += dstIx[ rSp] * stride[ rSp];
1173
1174 #endif // TEST_GOOD_OL_VERSION
1175
1176 // good 'ol version RELOADED
1177 SizeT* dim_stride = &stride[1];
1178 SizeT freeDstIx_0 = this_dim[ 0] - dstIx[ 0]; // how many elements till array border is reached (dim 0)
1179 // for( SizeT a=0; a<nEl; ++srcIx[0],++dstIx[0])
1180 for (SizeT a = 0; a < nEl; ++srcIx[1], ++dstIx[1]) {
1181 // for( SizeT aSp=0; aSp<nDim;)
1182 for (SizeT aSp = 1; aSp < nDim;) {
1183 if (dstIx[ aSp] >= this_dim[ aSp]) {
1184 // dstIx[ aSp] -= dim[ aSp];
1185 dstIx[ aSp] = 0;
1186 dstLonIx -= dim_stride[ aSp];
1187 }
1188 if (srcIx[ aSp] < this_dim[ aSp]) break;
1189
1190 srcIx[ aSp] = 0;
1191 if (++aSp >= nDim) break; // ??
1192
1193
1194 ++srcIx[ aSp];
1195 ++dstIx[ aSp];
1196 dstLonIx += stride[ aSp];
1197 }
1198
1199 // code from new version (to avoid the worst :-)
1200 // copy one line
1201 SizeT s = 0;
1202 for (; s < freeDstIx_0; ++s) {
1203 shP[ dstLonIx++] = ddP[ a++];
1204 }
1205 dstLonIx -= stride[ 1];
1206 for (; s < this_dim[ 0]; ++s) {
1207 shP[ dstLonIx++] = ddP[ a++];
1208 }
1209 dstLonIx += stride[ 1];
1210
1211 // copy one element
1212 //shP[ dstLonIx++] = ddP[ a++];
1213 }
1214
1215 return sh;
1216 }
1217
1218
1219
1220 // assumes *perm is already checked according to uniqness and length
1221 // dim[i]_out = dim[perm[i]]_in
1222 // helper function for Transpose()
InitPermDefault()1223 DUInt* InitPermDefault()
1224 {
1225 static DUInt res[ MAXRANK];
1226 for( SizeT i=MAXRANK-1, ii=0; ii<MAXRANK; ++ii, --i)
1227 res[ ii] = i;
1228 return res;
1229 }
1230 template<class Sp>
Transpose(DUInt * perm)1231 BaseGDL* Data_<Sp>::Transpose( DUInt* perm) {
1232 SizeT rank = this->Rank();
1233
1234 if (rank == 1) // special case: vector
1235 {
1236 if (perm != NULL) // must be [0]
1237 {
1238 return Dup();
1239 } else {
1240 Data_* res = Dup();
1241 res->dim >> 1;
1242 return res;
1243 }
1244 }
1245
1246 // 2 - MAXRANK
1247 static DUInt* permDefault = InitPermDefault();
1248 if (perm == NULL) {
1249
1250 // following 2D code is now slower than multi-dim multicore solution below.
1251 // if (rank == 2) {
1252 // SizeT srcDim0 = this->dim[0];
1253 // SizeT srcDim1 = this->dim[1];
1254 // Data_* res = new Data_(dimension(srcDim1, srcDim0), BaseGDL::NOZERO);
1255 //
1256 // SizeT srcIx = 0;
1257 // for (SizeT srcIx1 = 0; srcIx1 < srcDim1; ++srcIx1) // src dim 1
1258 // {
1259 // SizeT resIx = srcIx1;
1260 // SizeT srcLim = srcIx + srcDim0; // src dim 0
1261 // for (; srcIx < srcLim; ++srcIx) {
1262 // (*res)[ resIx] = (*this)[ srcIx];
1263 // resIx += srcDim1;
1264 // }
1265 // }
1266 // return res;
1267 // }
1268
1269 // perm == NULL, rank != 2
1270 perm = &permDefault[ MAXRANK - rank];
1271 }
1272
1273 //new version, parallell the job on a multithreaded machine. Gain is 3/4 number of threads.
1274 SizeT resDim[ MAXRANK]; // permutated!
1275 for (SizeT d = 0; d < rank; ++d) {
1276 resDim[ d] = this->dim[ perm[ d]];
1277 }
1278
1279 Data_* res = new Data_(dimension(resDim, rank), BaseGDL::NOZERO);
1280
1281 // src stride
1282 SizeT srcStride[ MAXRANK+1];
1283 this->dim.Stride(srcStride, rank);
1284
1285 SizeT nElem = dd.size();
1286 long chunksize = nElem;
1287 long nchunk = 1;
1288 if (nElem > CpuTPOOL_MIN_ELTS) { //no use start parallel threading for small numbers.
1289 chunksize = nElem / ((CpuTPOOL_NTHREADS > 32) ? 32 : CpuTPOOL_NTHREADS);
1290 nchunk = nElem / chunksize;
1291 if (chunksize * nchunk < nElem) ++nchunk;
1292 } else {
1293 nchunk = 1;
1294 chunksize = nElem;
1295 }
1296 //compute start parameter for each multiWalk chunks:
1297 // pool of accelerators
1298 SizeT srcDimPool[nchunk][MAXRANK];
1299 for (SizeT i = 0; i < rank; ++i) for (int iloop = 0; iloop < nchunk; ++iloop) srcDimPool[iloop][i] = 0;
1300 //template accelerator
1301 SizeT templateDim[MAXRANK];
1302 for (SizeT i = 0; i < rank; ++i) templateDim[i] = 0;
1303
1304 //compute iloop's accelerator with fast direct method
1305 for (long iloop = 0; iloop < nchunk; ++iloop) {
1306 SizeT e = iloop*chunksize;
1307 SizeT sizeleft = e;
1308 for (long i = 0; i < rank; ++i) {
1309 DUInt pi = perm[i]; //note the transpose effect.
1310 sizeleft /= resDim[i];
1311 templateDim[pi] = e - sizeleft * resDim[i];
1312 e = sizeleft;
1313 }
1314 //memorize current state accelerator for chunk iloop:
1315 for (long j = 0; j < rank; ++j) srcDimPool[iloop][j] = templateDim[j];
1316 }
1317
1318 #pragma omp parallel num_threads(nchunk)
1319 {
1320 #pragma omp for
1321 for (long iloop = 0; iloop < nchunk; ++iloop) {
1322 // populate src multi dim
1323 SizeT srcDim[MAXRANK];
1324 for (SizeT i = 0; i < rank; ++i) srcDim[i] = srcDimPool[iloop][i];
1325 //inner loop
1326 for (SizeT e = iloop * chunksize; (e < (iloop + 1) * chunksize && e < nElem); ++e) {
1327 // src multi dim to one dim src offset index
1328 SizeT ix = 0;
1329 for (SizeT i = 0; i < rank; ++i) ix += srcDim[i] * srcStride[i];
1330 (*res)[ e] = (*this)[ ix];
1331 // update src multi dim for next dest offset index: here is the transpose effect.
1332 for (SizeT i = 0; i < rank; ++i) {
1333 DUInt pi = perm[i];
1334 srcDim[pi]++;
1335 if (srcDim[pi] < resDim[i]) break;
1336 srcDim[pi] = 0;
1337 }
1338 }
1339 }
1340 }
1341 return res;
1342 }
1343
1344 // used by reverse
1345
1346 template<class Sp>
Reverse(DLong dim)1347 void Data_<Sp>::Reverse(DLong dim) {
1348 // SA: based on total_over_dim_template()
1349 // static Data_* tmp = new Data_(dimension(1), BaseGDL::NOZERO);
1350 //Guard<Data_> tmp_guard(tmp);
1351 SizeT nEl = N_Elements();
1352 SizeT revStride = this->dim.Stride(dim);
1353 SizeT outerStride = this->dim.Stride(dim + 1);
1354 SizeT revLimit = this->dim[dim] * revStride;
1355 #pragma omp parallel //if ((nEl/outerStride)*revStride >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= (nEl/outerStride)*revStride))
1356 {
1357 #pragma omp for
1358 for (SizeT o = 0; o < nEl; o += outerStride) {
1359 for (SizeT i = 0; i < revStride; ++i) {
1360 SizeT oi = o + i;
1361 SizeT last_plus_oi = revLimit + oi - revStride + oi;
1362 SizeT half = ((revLimit / revStride) / 2) * revStride + oi;
1363 for (SizeT s = oi; s < half; s += revStride) {
1364 SizeT opp = last_plus_oi - s;
1365 Ty tmp = (*this)[s];
1366 (*this)[s] = (*this)[opp];
1367 (*this)[opp] = tmp;
1368 }
1369 }
1370 }
1371 }
1372 }
1373
1374 template<class Sp>
DupReverse(DLong dim)1375 BaseGDL* Data_<Sp>::DupReverse(DLong dim) {
1376 // SA: based on total_over_dim_template()
1377 Data_* res = new Data_(this->dim, BaseGDL::NOZERO);
1378 Guard<Data_> res_guard(res);
1379 SizeT nEl = N_Elements();
1380 SizeT revStride = this->dim.Stride(dim);
1381 SizeT outerStride = this->dim.Stride(dim + 1);
1382 SizeT revLimit = this->dim[dim] * revStride;
1383 #pragma omp parallel //if ((nEl/outerStride)*revStride >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= (nEl/outerStride)*revStride))
1384 {
1385 #pragma omp for
1386 for (SizeT o = 0; o < nEl; o += outerStride) {
1387 for (SizeT i = 0; i < revStride; ++i) {
1388 SizeT oi = o + i;
1389 SizeT last_plus_oi = revLimit + oi - revStride + oi;
1390 SizeT half = ((revLimit / revStride) / 2) * revStride + oi;
1391 for (SizeT s = oi; s < half + 1; s += revStride) {
1392 SizeT opp = last_plus_oi - s;
1393 // cout << s <<" "<< opp << " " << (*this)[s] << " " << (*this)[opp] << endl;
1394 (*res)[s] = (*this)[opp];
1395 (*res)[opp] = (*this)[s];
1396 }
1397 }
1398 }
1399 }
1400 return res_guard.release();
1401 }
1402 template<>
DupReverse(DLong dim)1403 BaseGDL* Data_<SpDPtr>::DupReverse( DLong dim) {
1404 // SA: based on total_over_dim_template()
1405 Data_* res = new Data_(this->dim, BaseGDL::NOZERO);
1406 Guard<Data_> res_guard(res);
1407 SizeT nEl = N_Elements();
1408 SizeT revStride = this->dim.Stride(dim);
1409 SizeT outerStride = this->dim.Stride(dim + 1);
1410 SizeT revLimit = this->dim[dim] * revStride;
1411 #pragma omp parallel //if ((nEl/outerStride)*revStride >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= (nEl/outerStride)*revStride))
1412 {
1413 #pragma omp for
1414 for (SizeT o = 0; o < nEl; o += outerStride) {
1415 for (SizeT i = 0; i < revStride; ++i) {
1416 SizeT oi = o + i;
1417 SizeT last_plus_oi = revLimit + oi - revStride + oi;
1418 SizeT half = ((revLimit / revStride) / 2) * revStride + oi;
1419 for (SizeT s = oi; s < half + 1; s += revStride) {
1420 SizeT opp = last_plus_oi - s;
1421 (*res)[s] = (*this)[opp];
1422 (*res)[opp] = (*this)[s];
1423 }
1424 }
1425 }
1426 }
1427 GDLInterpreter::IncRef(res);
1428 return res_guard.release();
1429 }
1430 template<>
DupReverse(DLong dim)1431 BaseGDL* Data_<SpDObj>::DupReverse(DLong dim) {
1432 // SA: based on total_over_dim_template()
1433 Data_* res = new Data_(this->dim, BaseGDL::NOZERO);
1434 Guard<Data_> res_guard(res);
1435 SizeT nEl = N_Elements();
1436 SizeT revStride = this->dim.Stride(dim);
1437 SizeT outerStride = this->dim.Stride(dim + 1);
1438 SizeT revLimit = this->dim[dim] * revStride;
1439 #pragma omp parallel //if ((nEl/outerStride)*revStride >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= (nEl/outerStride)*revStride))
1440 {
1441 #pragma omp for
1442 for (SizeT o = 0; o < nEl; o += outerStride) {
1443 for (SizeT i = 0; i < revStride; ++i) {
1444 SizeT oi = o + i;
1445 SizeT last_plus_oi = revLimit + oi - revStride + oi;
1446 SizeT half = ((revLimit / revStride) / 2) * revStride + oi;
1447 for (SizeT s = oi; s < half + 1; s += revStride) {
1448 SizeT opp = last_plus_oi - s;
1449 (*res)[s] = (*this)[opp];
1450 (*res)[opp] = (*this)[s];
1451 }
1452 }
1453 }
1454 }
1455 GDLInterpreter::IncRefObj(res);
1456 return res_guard.release();
1457 }
1458
1459 // rank must be 1 or 2 (already checked)
1460 template<class Sp>
Rotate(DLong dir)1461 BaseGDL* Data_<Sp>::Rotate( DLong dir) {
1462 dir = (dir % 8 + 8) % 8; // bring into 0..7 range
1463
1464 if (dir == 0) return Dup();
1465 if (dir == 2) {
1466 Data_* res = new Data_(this->dim, BaseGDL::NOZERO);
1467 SizeT nEl = N_Elements();
1468 //no pragma: better optimized directly by the compiler!
1469 {
1470 for (SizeT i = 0; i < nEl; ++i)
1471 (*res)[i] = (*this)[ nEl - 1 - i];
1472 }
1473 return res;
1474 }
1475
1476 if (this->Rank() == 1) {
1477 if (dir == 7) return Dup();
1478
1479 if (dir == 1 || dir == 4) {
1480 return new Data_(dimension(1, N_Elements()), dd);
1481 }
1482 if (dir == 5) // || dir == 2
1483 {
1484 Data_* res = new Data_(this->dim, BaseGDL::NOZERO);
1485 SizeT nEl = N_Elements();
1486 //no pragma: better optimized directly by the compiler!
1487 {
1488 for (SizeT i = 0; i < nEl; ++i)
1489 (*res)[ i] = (*this)[ nEl - 1 - i];
1490 }
1491 return res;
1492 }
1493 // 3 || 6
1494 Data_* res = new Data_(dimension(1, N_Elements()), BaseGDL::NOZERO);
1495 SizeT nEl = N_Elements();
1496 //no pragma: better optimized directly by the compiler!
1497 {
1498 for (SizeT i = 0; i < nEl; ++i)
1499 (*res)[ i] = (*this)[ nEl - 1 - i];
1500 }
1501 return res;
1502 }
1503
1504 // rank == 2, dir == 0 and dir == 2 already handled
1505 bool keepDim = (dir == 5) || (dir == 7);
1506
1507 Data_* res;
1508 if (keepDim) {
1509 res = new Data_(this->dim, BaseGDL::NOZERO);
1510 } else {
1511 res = new Data_(dimension(this->dim[1], this->dim[0]), BaseGDL::NOZERO);
1512 }
1513
1514 // bool flipX = dir == 3 || dir == 5 || dir == 6;
1515 // bool flipY = dir == 1 || dir == 6 || dir == 7;
1516
1517 SizeT xEl = this->dim[0];
1518 SizeT yEl = this->dim[1];
1519 SizeT i = 0;
1520 // enable fast optimzation by removing ifs outside loops.
1521 if (dir == 1) { // flipY
1522 for (SizeT y = 0; y < yEl; ++y) for (SizeT x = 0; x < xEl; ++x) (*res)[ x * yEl + yEl - 1 - y] = (*this)[ i++];
1523 } else if (dir == 3) { // flipX
1524 for (SizeT y = 0; y < yEl; ++y) for (SizeT x = 0; x < xEl; ++x) (*res)[ (xEl - 1 - x) * yEl + y] = (*this)[ i++];
1525 } else if (dir == 4) {
1526 for (SizeT y = 0; y < yEl; ++y) for (SizeT x = 0; x < xEl; ++x) (*res)[x * yEl + y] = (*this)[ i++];
1527 } else if (dir == 5) { // flipX, keepDim
1528 for (SizeT y = 0; y < yEl; ++y) for (SizeT x = 0; x < xEl; ++x) (*res)[y * xEl + xEl - 1 - x] = (*this)[ i++];
1529 } else if (dir == 6) { // flipX, flipY
1530 for (SizeT y = 0; y < yEl; ++y) for (SizeT x = 0; x < xEl; ++x) (*res)[(xEl - 1 - x )* yEl + yEl - 1 - y] = (*this)[ i++];
1531 } else if (dir == 7) { // flipY, keepDim
1532 for (SizeT y = 0; y < yEl; ++y) for (SizeT x = 0; x < xEl; ++x) (*res)[(yEl - 1 - y) * xEl + x ] = (*this)[ i++];
1533 }
1534 //was:
1535 // for (SizeT y = 0; y < yEl; ++y) {
1536 // SizeT yR = flipY ? yEl - 1 - y : y;
1537 // for (SizeT x = 0; x < xEl; ++x) {
1538 // SizeT xR = flipX ? xEl - 1 - x : x;
1539 //
1540 // SizeT ix = keepDim ? yR * xEl + xR : xR * yEl + yR;
1541 //
1542 // (*res)[ix] = (*this)[ i++];
1543 // }
1544 // }
1545 return res;
1546 }
1547
1548 template<class Sp>
Sum() const1549 typename Data_<Sp>::Ty Data_<Sp>::Sum() const
1550 {
1551 Ty s= dd[ 0];
1552 SizeT nEl = dd.size();
1553 TRACEOMP( __FILE__, __LINE__)
1554 #pragma omp parallel if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl)) shared( s)
1555 {
1556 #pragma omp for reduction(+:s)
1557 for( SizeT i=1; i<nEl; ++i)
1558 {
1559 s += dd[ i];
1560 }
1561 }
1562 return s;
1563 }
1564
1565 template<>
Sum() const1566 Data_<SpDString>::Ty Data_<SpDString>::Sum() const
1567 {
1568 Ty s= dd[ 0];
1569 SizeT nEl = dd.size();
1570 for( SizeT i=1; i<nEl; ++i)
1571 {
1572 s += dd[ i];
1573 }
1574 return s;
1575 }
1576 template<>
Sum() const1577 Data_<SpDComplexDbl>::Ty Data_<SpDComplexDbl>::Sum() const
1578 {
1579 DDouble sr= dd[ 0].real();
1580 DDouble si= dd[ 0].imag();
1581 SizeT nEl = dd.size();
1582 TRACEOMP( __FILE__, __LINE__)
1583 #pragma omp parallel if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl)) shared( sr,si)
1584 {
1585 #pragma omp for reduction(+:si,sr)
1586 for( SizeT i=1; i<nEl; ++i)
1587 {
1588 sr += dd[i].real();
1589 si += dd[i].imag();
1590 }
1591 }
1592 return std::complex<double>(sr,si);
1593 }
1594 template<>
Sum() const1595 Data_<SpDComplex>::Ty Data_<SpDComplex>::Sum() const
1596 {
1597 DFloat sr= dd[ 0].real();
1598 DFloat si= dd[ 0].imag();
1599 SizeT nEl = dd.size();
1600 TRACEOMP( __FILE__, __LINE__)
1601 #pragma omp parallel if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl)) shared( sr,si)
1602 {
1603 #pragma omp for reduction(+:si,sr)
1604 for( SizeT i=1; i<nEl; ++i)
1605 {
1606 sr += dd[i].real();
1607 si += dd[i].imag();
1608 }
1609 }
1610 return std::complex<float>(sr,si);
1611 }
1612
1613 // template<class Sp>
1614 // typename Data_<Sp>::DataT& Data_<Sp>:: Resize( SizeT n)
1615 // {
1616 // if( n > dd.size())
1617 // throw GDLException("Internal error: Data_::Resize(...) tried to grow.");
1618 // if( dd.size() != n)
1619 // {
1620 // DataT rsArr( n);
1621 // std::copy( &(*this)[0], &(*this)[n], &rsArr[0]);
1622 // dd.resize( n); // discards data
1623 // std::copy( &rsArr[0], &rsArr[n], &(*this)[0]);
1624 // }
1625 // return dd;
1626 // }
1627
1628 // template<class Sp>
1629 // typename Data_<Sp>::Ty& Data_<Sp>::operator[] (const SizeT d1)
1630 // { return (*this)[d1];}
1631
1632 // only used from DStructGDL
1633 template<class Sp>
operator =(const BaseGDL & r)1634 Data_<Sp>& Data_<Sp>::operator=(const BaseGDL& r)
1635 {
1636 assert( r.Type() == this->Type());
1637 const Data_<Sp>& right = static_cast<const Data_<Sp>&>( r);
1638 assert( &right != this);
1639 // if( &right == this) return *this; // self assignment
1640 this->dim = right.dim;
1641 dd = right.dd;
1642 return *this;
1643 }
1644 // only used from DStructGDL
1645 template<>
operator =(const BaseGDL & r)1646 Data_<SpDPtr>& Data_<SpDPtr>::operator=(const BaseGDL& r)
1647 {
1648 assert( r.Type() == this->Type());
1649 const Data_& right = static_cast<const Data_&>( r);
1650 assert( &right != this);
1651 // if( &right == this) return *this; // self assignment
1652 this->dim = right.dim;
1653 GDLInterpreter::DecRef( this);
1654 dd = right.dd;
1655 GDLInterpreter::IncRef( this);
1656 return *this;
1657 }
1658 // only used from DStructGDL
1659 template<>
operator =(const BaseGDL & r)1660 Data_<SpDObj>& Data_<SpDObj>::operator=(const BaseGDL& r)
1661 {
1662 assert( r.Type() == this->Type());
1663 const Data_& right = static_cast<const Data_&>( r);
1664 assert( &right != this);
1665 // if( &right == this) return *this; // self assignment
1666 this->dim = right.dim;
1667 GDLInterpreter::DecRefObj( this);
1668 dd = right.dd;
1669 GDLInterpreter::IncRefObj( this);
1670 return *this;
1671 }
1672
1673
1674 template<class Sp>
InitFrom(const BaseGDL & r)1675 void Data_<Sp>::InitFrom(const BaseGDL& r)
1676 {
1677 assert( r.Type() == this->Type());
1678 const Data_<Sp>& right = static_cast<const Data_<Sp>&>( r);
1679 assert( &right != this);
1680 // if( &right == this) return *this; // self assignment
1681 this->dim = right.dim;
1682 dd.InitFrom( right.dd);
1683 // return *this;
1684 }
1685 // only used from DStructGDL::DStructGDL(const DStructGDL& d_)
1686 template<>
InitFrom(const BaseGDL & r)1687 void Data_<SpDPtr>::InitFrom(const BaseGDL& r)
1688 {
1689 assert( r.Type() == this->Type());
1690 const Data_& right = static_cast<const Data_&>( r);
1691 assert( &right != this);
1692 // if( &right == this) return *this; // self assignment
1693 this->dim = right.dim;
1694 // GDLInterpreter::DecRef( this);
1695 dd.InitFrom( right.dd);
1696 GDLInterpreter::IncRef( this);
1697 // return *this;
1698 }
1699 // only used from DStructGDL::DStructGDL(const DStructGDL& d_)
1700 template<>
InitFrom(const BaseGDL & r)1701 void Data_<SpDObj>::InitFrom(const BaseGDL& r)
1702 {
1703 assert( r.Type() == this->Type());
1704 const Data_& right = static_cast<const Data_&>( r);
1705 assert( &right != this);
1706 // if( &right == this) return *this; // self assignment
1707 this->dim = right.dim;
1708 // GDLInterpreter::DecRefObj( this);
1709 dd.InitFrom( right.dd);
1710 GDLInterpreter::IncRefObj( this);
1711 // return *this;
1712 }
1713
1714 template< class Sp>
EqType(const BaseGDL * r) const1715 bool Data_<Sp>::EqType( const BaseGDL* r) const
1716 { return (this->Type() == r->Type());}
1717
1718 template< class Sp>
DataAddr()1719 void* Data_<Sp>::DataAddr()// SizeT elem)
1720 { return &(*this)[0];}//elem];}
1721
1722 // template<>
1723 // void* Data_<SpDString>::DataAddr()// SizeT elem)
1724 // {
1725 // throw GDLException("STRING not allowed in this context.");
1726 // return NULL;
1727 // }
1728 // template<>
1729 // void* Data_<SpDPtr>::DataAddr()// SizeT elem)
1730 // {
1731 // throw GDLException("PTR not allowed in this context.");
1732 // return NULL;
1733 // }
1734 // template<>
1735 // void* Data_<SpDObj>::DataAddr()// SizeT elem)
1736 // {
1737 // throw GDLException("Object expression not allowed in this context.");
1738 // return NULL;
1739 // }
1740
1741 template< class Sp>
N_Elements() const1742 SizeT Data_<Sp>::N_Elements() const
1743 { return dd.size();}
1744
1745 template<>
N_Elements() const1746 SizeT Data_<SpDObj>::N_Elements() const
1747 {
1748 if( !this->StrictScalar())
1749 return dd.size();
1750 #if 1 // GJ change this and see how many problems it solves, creates. 2016/5/5
1751 else return 1;
1752 #elif 0
1753 DObj s = dd[0]; // is StrictScalar()
1754 if( s == 0) // no overloads for null object
1755 return 1;
1756
1757 DStructGDL* oStructGDL= GDLInterpreter::GetObjHeapNoThrow( s);
1758 if( oStructGDL == NULL) // if object not valid -> default behaviour
1759 return 1;
1760
1761 DStructDesc* desc = oStructGDL->Desc();
1762
1763 if( desc->IsParent("LIST"))
1764 {
1765 // no static here, might vary in derived object
1766 unsigned nListTag = desc->TagIndex( "NLIST");
1767 SizeT listSize = (*static_cast<DLongGDL*>(oStructGDL->GetTag( nListTag, 0)))[0];
1768 return listSize;
1769 }
1770 if( desc->IsParent("HASH"))
1771 {
1772 // no static here, might vary in derived object
1773 unsigned nListTag = desc->TagIndex( "TABLE_COUNT");
1774 SizeT listSize = (*static_cast<DLongGDL*>(oStructGDL->GetTag( nListTag, 0)))[0];
1775 return listSize;
1776 }
1777
1778 return 1;
1779 #endif
1780 }
1781
1782
1783 template< class Sp>
Size() const1784 SizeT Data_<Sp>::Size() const
1785 { return dd.size();}
1786 template< class Sp>
Sizeof() const1787 SizeT Data_<Sp>::Sizeof() const
1788 { return sizeof(Ty);}
1789
1790 template< class Sp>
Clear()1791 void Data_<Sp>::Clear()
1792 {
1793 SizeT nEl = dd.size();
1794 /*#pragma omp parallel if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl))
1795 {
1796 #pragma omp for*/
1797 for( SizeT i = 0; i<nEl; ++i) (*this)[ i] = Sp::zero;
1798 }//}
1799
1800 // first time initialization (construction)
1801 template< class Sp>
Construct()1802 void Data_<Sp>::Construct()
1803 {
1804 // note that this is not possible in cases where an operation
1805 // (here: 'new' which is ok) isn't defined for any POD
1806 // (although this code never executes and should be optimized away anyway)
1807 const bool isPOD = Sp::IS_POD;
1808 // do nothing for POD
1809 if( !isPOD)
1810 {
1811 SizeT nEl = dd.size();
1812 // for( SizeT i = 0; i<nEl; ++i) new (&(*this)[ i]) Ty;
1813 /*#pragma omp parallel if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl))
1814 {
1815 #pragma omp for*/
1816 for( SizeT i = 0; i<nEl; ++i) new (&(dd[ i])) Ty;
1817 }
1818 }
1819 template<>
Construct()1820 void Data_<SpDPtr>::Construct()
1821 {
1822 SizeT nEl = dd.size();
1823 // for( SizeT i = 0; i<nEl; ++i) new (&(*this)[ i]) Ty;
1824 /*#pragma omp parallel if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl))
1825 {
1826 #pragma omp for*/
1827 for( SizeT i = 0; i<nEl; ++i) dd[ i] = 0;
1828 }//}
1829 template<>
Construct()1830 void Data_<SpDObj>::Construct()
1831 {
1832 SizeT nEl = dd.size();
1833 // for( SizeT i = 0; i<nEl; ++i) new (&(*this)[ i]) Ty;
1834 /*#pragma omp parallel if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl))
1835 {
1836 #pragma omp for*/
1837 for( SizeT i = 0; i<nEl; ++i) dd[ i] = 0;
1838 }//}
1839 // // non POD - use placement new
1840 // template<>
1841 // void Data_< SpDString>::Construct()
1842 // {
1843 // SizeT nEl = dd.size();
1844 // // for( SizeT i = 0; i<nEl; ++i) new (&(*this)[ i]) Ty;
1845 // /*#pragma omp parallel if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl))
1846 // {
1847 // #pragma omp for*/
1848 // for( SizeT i = 0; i<nEl; ++i) new (&(dd[ i])) Ty;
1849 // }//}
1850 // template<>
1851 // void Data_< SpDComplex>::Construct()
1852 // {
1853 // SizeT nEl = dd.size();
1854 // /*#pragma omp parallel if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl))
1855 // {
1856 // #pragma omp for*/
1857 // for( SizeT i = 0; i<nEl; ++i) new (&(*this)[ i]) Ty;
1858 // }//}
1859 // template<>
1860 // void Data_< SpDComplexDbl>::Construct()
1861 // {
1862 // SizeT nEl = dd.size();
1863 // /*#pragma omp parallel if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl))
1864 // {
1865 // #pragma omp for*/
1866 // for( SizeT i = 0; i<nEl; ++i) new (&(*this)[ i]) Ty;
1867 // }//}
1868
1869 // construction and initalization to zero
1870 template< class Sp>
ConstructTo0()1871 void Data_<Sp>::ConstructTo0()
1872 {
1873 if( Sp::IS_POD)
1874 {
1875 SizeT nEl = dd.size();
1876 /*#pragma omp parallel if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl))
1877 {
1878 #pragma omp for*/
1879 for( SizeT i = 0; i<nEl; ++i) (*this)[ i] = Sp::zero;
1880 }
1881 else
1882 {
1883 SizeT nEl = dd.size();
1884 /*#pragma omp parallel if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl))
1885 {
1886 #pragma omp for*/
1887 for( SizeT i = 0; i<nEl; ++i) new (&(*this)[ i]) Ty( Sp::zero);
1888 }
1889 }//}
1890 // // non POD - use placement new
1891 // template<>
1892 // void Data_< SpDString>::ConstructTo0()
1893 // {
1894 // SizeT nEl = dd.size();
1895 // /*#pragma omp parallel if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl))
1896 // {
1897 // #pragma omp for*/
1898 // for( SizeT i = 0; i<nEl; ++i) new (&(*this)[ i]) Ty( zero);
1899 // }//}
1900 // template<>
1901 // void Data_< SpDComplex>::ConstructTo0()
1902 // {
1903 // SizeT nEl = dd.size();
1904 // /*#pragma omp parallel if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl))
1905 // {
1906 // #pragma omp for*/
1907 // for( SizeT i = 0; i<nEl; ++i) new (&(*this)[ i]) Ty( zero);
1908 // }//}
1909 // template<>
1910 // void Data_< SpDComplexDbl>::ConstructTo0()
1911 // {
1912 // SizeT nEl = dd.size();
1913 // /*#pragma omp parallel if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl))
1914 // {
1915 // #pragma omp for*/
1916 // for( SizeT i = 0; i<nEl; ++i) new (&(*this)[ i]) Ty( zero);
1917 // }//}
1918
1919 template< class Sp>
Destruct()1920 void Data_<Sp>::Destruct()
1921 {
1922 // no destruction for POD
1923 if( !Sp::IS_POD)
1924 {
1925 SizeT nEl = dd.size();
1926 /*#pragma omp parallel if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl))
1927 {
1928 #pragma omp for*/
1929 for( SizeT i = 0; i<nEl; ++i)
1930 (*this)[ i].~Ty();
1931 }
1932 }
1933 template<>
Destruct()1934 void Data_< SpDPtr>::Destruct()
1935 {
1936 GDLInterpreter::DecRef( this);
1937 }
1938 template<>
Destruct()1939 void Data_< SpDObj>::Destruct()
1940 {
1941 GDLInterpreter::DecRefObj( this);
1942 }
1943 // template<>
1944 // void Data_< SpDString>::Destruct()
1945 // {
1946 // SizeT nEl = dd.size();
1947 // /*#pragma omp parallel if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl))
1948 // {
1949 // #pragma omp for*/
1950 // for( SizeT i = 0; i<nEl; ++i)
1951 // (*this)[ i].~DString();
1952 // }//}
1953 // template<>
1954 // void Data_< SpDComplex>::Destruct()
1955 // {
1956 // SizeT nEl = dd.size();
1957 // /*#pragma omp parallel if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl))
1958 // {
1959 // #pragma omp for*/
1960 // for( SizeT i = 0; i<nEl; ++i)
1961 // (*this)[ i].~DComplex();
1962 // }//}
1963 // template<>
1964 // void Data_< SpDComplexDbl>::Destruct()
1965 // {
1966 // SizeT nEl = dd.size();
1967 // /*#pragma omp parallel if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl))
1968 // {
1969 // #pragma omp for*/
1970 // for( SizeT i = 0; i<nEl; ++i)
1971 // (*this)[ i].~DComplexDbl();
1972 // }//}
1973
1974 template< class Sp>
SetBuffer(const void * b)1975 BaseGDL* Data_<Sp>::SetBuffer( const void* b)
1976 {
1977 dd.SetBuffer( static_cast< Ty*>(const_cast<void*>( b)));
1978 return this;
1979 }
1980 template< class Sp>
SetBufferSize(SizeT s)1981 void Data_<Sp>::SetBufferSize( SizeT s)
1982 {
1983 dd.SetBufferSize( s);
1984 }
1985
1986 // template< class Sp>
1987 // Data_<Sp>* Data_<Sp>::Dup()
1988 // { return new Data_(*this);}
1989
1990
1991 template< class Sp>
New(const dimension & dim_,BaseGDL::InitType noZero) const1992 Data_<Sp>* Data_<Sp>::New( const dimension& dim_, BaseGDL::InitType noZero) const
1993 {
1994 if( noZero == BaseGDL::NOZERO) return new Data_(dim_, BaseGDL::NOZERO);
1995 if( noZero == BaseGDL::INIT)
1996 {
1997 Data_* res = new Data_(dim_, BaseGDL::NOZERO);
1998 SizeT nEl = res->dd.size();
1999 /*#pragma omp parallel if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl))
2000 {
2001 #pragma omp for*/
2002 for( SizeT i=0; i<nEl; ++i) (*res)[ i] = (*this)[ 0]; // set all to scalar
2003 //}
2004 return res;
2005 }
2006 return new Data_(dim_); // zero data
2007 }
2008
2009
2010 template< class Sp>
NewResult() const2011 Data_<Sp>* Data_<Sp>::NewResult() const
2012 {
2013 return new Data_(this->dim, BaseGDL::NOZERO);
2014 }
2015
2016 // template< class Sp>
2017 // bool Data_<Sp>::Scalar() const
2018 // {
2019 // return (dd.size() == 1);
2020 // }
2021
2022 // template< class Sp>
2023 // bool Data_<Sp>::Scalar(Ty& s) const
2024 // {
2025 // if( dd.size() != 1) return false;
2026 // s=(*this)[0];
2027 // return true;
2028 // }
2029
2030
2031 template<class Sp>
NBytes() const2032 SizeT Data_<Sp>::NBytes() const
2033 { return (dd.size() * sizeof(Ty));}
2034
2035 // string, ptr, obj cannot calculate their bytes
2036 // only used by assoc function
NBytes() const2037 template<> SizeT Data_<SpDString>::NBytes() const
2038 {
2039 SizeT nEl = dd.size();
2040 SizeT nB = 0;
2041 #pragma omp parallel for if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl)) reduction(+:nB)
2042 for( SizeT i=0; i<nEl; ++i) nB += (*this)[i].size();
2043 return nB;
2044 }
2045 // template<> SizeT Data_<SpDObj>::NBytes() const
2046 // {
2047 // throw GDLException("object references");
2048 // }
2049 // template<> SizeT Data_<SpDPtr>::NBytes() const
2050 // {
2051 // throw GDLException("pointers");
2052 // }
2053
2054 template<class Sp>
ToTransfer() const2055 SizeT Data_<Sp>::ToTransfer() const
2056 {
2057 return dd.size();
2058 }
2059 // complex has 2 elements to transfer
ToTransfer() const2060 template<> SizeT Data_<SpDComplex>::ToTransfer() const
2061 {
2062 return N_Elements() * 2;
2063 }
ToTransfer() const2064 template<> SizeT Data_<SpDComplexDbl>::ToTransfer() const
2065 {
2066 return N_Elements() * 2;
2067 }
2068
2069 // // note that min and max are not defined in BaseGDL
2070 // template<class Sp>
2071 // typename Data_<Sp>::Ty Data_<Sp>::min() const
2072 // { return dd.min();}
2073 // template<class Sp>
2074 // typename Data_<Sp>::Ty Data_<Sp>::max() const
2075 // { return dd.max();}
2076 // template<>
2077 // Data_<SpDComplex>::Ty Data_<SpDComplex>::min() const
2078 // {
2079 // throw GDLException("COMPLEX expression not allowed in this context.");
2080 // }
2081 // template<>
2082 // Data_<SpDComplex>::Ty Data_<SpDComplex>::max() const
2083 // {
2084 // throw GDLException("COMPLEX expression not allowed in this context.");
2085 // }
2086 // template<>
2087 // Data_<SpDComplexDbl>::Ty Data_<SpDComplexDbl>::min() const
2088 // {
2089 // throw GDLException("COMPLEXDBL expression not allowed in this context.");
2090 // }
2091 // template<>
2092 // Data_<SpDComplexDbl>::Ty Data_<SpDComplexDbl>::max() const
2093 // {
2094 // throw GDLException("COMPLEXDBL expression not allowed in this context.");
2095 // }
2096
2097
2098 // for HASH objects
2099 template<class Sp>
HashValue() const2100 DDouble Data_<Sp>::HashValue() const
2101 {
2102 return static_cast<DDouble>((*this)[0]);
2103 }
2104 template<>
HashValue() const2105 DDouble Data_<SpDComplex>::HashValue() const
2106 {
2107 return real((*this)[0]);
2108 }
2109 template<>
HashValue() const2110 DDouble Data_<SpDComplexDbl>::HashValue() const
2111 {
2112 return real((*this)[0]);
2113 }
2114 template<>
HashValue() const2115 DDouble Data_<SpDString>::HashValue() const
2116 {
2117 throw GDLException("STRING expression not allowed as index. Please report.");
2118 return 0; // get rid of warning
2119 }
2120 template<>
HashValue() const2121 DDouble Data_<SpDPtr>::HashValue() const
2122 {
2123 throw GDLException("PTR expression not allowed as index. Please report.");
2124 return 0; // get rid of warning
2125 }
2126
2127 template<>
HashValue() const2128 DDouble Data_<SpDObj>::HashValue() const
2129 {
2130 throw GDLException("Object expression not allowed as index. Please report.");
2131 return 0; // get rid of warning
2132 }
2133
2134
2135 // -1 -> p2 is greater
2136 // 0 -> equal
2137 // 1 -> this is greater
2138
2139 // note: this is for internal use only (for HASH objects)
2140 // this should not be called on non-numeric types (also for p2)
2141 template<class Sp>
HashCompare(BaseGDL * p2) const2142 int Data_<Sp>::HashCompare( BaseGDL* p2) const
2143 {
2144 assert( dd.size() == 1);
2145 assert( p2->N_Elements() == 1);
2146 if( p2->Type() == GDL_STRING)
2147 return 1; // strings 1st (smallest)
2148
2149 assert( NumericType(p2->Type()));
2150
2151 if( this->IS_INTEGER)
2152 {
2153 if( IntType( p2->Type())) // make full use of data type
2154 {
2155 RangeT thisValue = this->LoopIndex();
2156 RangeT p2Value = p2->LoopIndex();
2157 if( thisValue == p2Value)
2158 return 0;
2159 if( thisValue < p2Value)
2160 return -1;
2161 return 1;
2162 }
2163 }
2164 DDouble thisValue = this->HashValue();
2165 DDouble p2Value = p2->HashValue();
2166 if( thisValue == p2Value)
2167 return 0;
2168 if( thisValue < p2Value)
2169 return -1;
2170 return 1;
2171 }
2172
2173 template<>
HashCompare(BaseGDL * p2) const2174 int Data_<SpDString>::HashCompare( BaseGDL* p2) const
2175 {
2176 assert( dd.size() == 1);
2177 assert( p2->N_Elements() == 1);
2178 if( p2->Type() != this->Type())
2179 return -1; // strings 1st (smallest)
2180
2181 Data_* p2String = static_cast<Data_*>(p2);
2182 if( dd[0].length() == (*p2String)[0].length())
2183 {
2184 if( dd[0] == (*p2String)[0])
2185 return 0;
2186 if( dd[0] < (*p2String)[0])
2187 return -1;
2188 return 1;
2189 }
2190 else if( dd[0].length() < (*p2String)[0].length())
2191 {
2192 return -1;
2193 }
2194 return 1;
2195 }
2196
2197
2198 // Scalar2Index
2199 // used by the interpreter
2200 // -2 < 0 array
2201 // -1 < 0 scalar
2202 // 0 empty or array
2203 // 1 scalar
2204 // 2 one-element array
2205 template<class Sp>
Scalar2Index(SizeT & st) const2206 int Data_<Sp>::Scalar2Index( SizeT& st) const
2207 {
2208 if( dd.size() != 1) return 0;
2209
2210 // the next statement gives a warning for unsigned integer types:
2211 // "comparison is always false due to limited range of data type"
2212 // This is because the same template is used here for signed and
2213 // unsigned data. A template specialization for the unsigned integer
2214 // types would result in three identical specializations, which is bad
2215 // for maintainability. And as any modern C++ compiler will optimize
2216 // away the superflous (for unsigned data) statement anyway, it is
2217 // better to keep the code this way here.
2218 if( (*this)[0] < 0) {
2219 if( this->dim.Rank() != 0)
2220 return -2;
2221 else
2222 return -1;
2223 }
2224
2225 st= static_cast<SizeT>((*this)[0]);
2226 if( this->dim.Rank() != 0) return 2;
2227 return 1;
2228 }
2229
2230 template<class Sp>
Scalar2RangeT(RangeT & st) const2231 int Data_<Sp>::Scalar2RangeT( RangeT& st) const
2232 {
2233 if( dd.size() != 1) return 0;
2234
2235 st= static_cast<RangeT>((*this)[0]);
2236 if( this->dim.Rank() != 0) return 2;
2237 return 1;
2238 }
2239
2240 template<>
Scalar2Index(SizeT & st) const2241 int Data_<SpDComplex>::Scalar2Index( SizeT& st) const
2242 {
2243 if( dd.size() != 1) return 0;
2244 float r=real((*this)[0]);
2245 if( r < 0.0) return -1;
2246 st= static_cast<SizeT>(r);
2247 if( this->dim.Rank() != 0) return 2;
2248 return 1;
2249 }
2250
2251 template<>
Scalar2RangeT(RangeT & st) const2252 int Data_<SpDComplex>::Scalar2RangeT( RangeT& st) const
2253 {
2254 if( dd.size() != 1) return 0;
2255 float r=real((*this)[0]);
2256 st= static_cast<RangeT>(r);
2257 if( this->dim.Rank() != 0) return 2;
2258 return 1;
2259 }
2260
2261 template<>
Scalar2Index(SizeT & st) const2262 int Data_<SpDComplexDbl>::Scalar2Index( SizeT& st) const
2263 {
2264 if( dd.size() != 1) return 0;
2265 double r=real((*this)[0]);
2266 if( r < 0.0) return -1;
2267 st= static_cast<SizeT>(r);
2268 if( this->dim.Rank() != 0) return 2;
2269 return 1;
2270 }
2271
2272 template<>
Scalar2RangeT(RangeT & st) const2273 int Data_<SpDComplexDbl>::Scalar2RangeT( RangeT& st) const
2274 {
2275 if( dd.size() != 1) return 0;
2276 double r=real((*this)[0]);
2277 st= static_cast<RangeT>(r);
2278 if( this->dim.Rank() != 0) return 2;
2279 return 1;
2280 }
2281
2282
2283 template<>
Scalar2Index(SizeT & st) const2284 int Data_<SpDString>::Scalar2Index( SizeT& st) const
2285 {
2286 if( dd.size() != 1) return 0;
2287
2288 SizeT tSize=(*this)[0].size();
2289
2290 if( tSize == 0)
2291 {
2292 st=0;
2293 }
2294 else
2295 {
2296 long int number = Str2L( (*this)[0].c_str());
2297 if( number < 0) return -1;
2298 st=number;
2299 }
2300 if( dim.Rank() != 0) return 2;
2301 return 1;
2302 }
2303
2304 template<>
Scalar2RangeT(RangeT & st) const2305 int Data_<SpDString>::Scalar2RangeT( RangeT& st) const
2306 {
2307 if( dd.size() != 1) return 0;
2308
2309 SizeT tSize=(*this)[0].size();
2310
2311 if( tSize == 0)
2312 {
2313 st=0;
2314 }
2315 else
2316 {
2317 long int number = Str2L( (*this)[0].c_str());
2318 st=number;
2319 }
2320 if( dim.Rank() != 0) return 2;
2321 return 1;
2322 }
2323
2324
2325 template<>
Scalar2Index(SizeT & st) const2326 int Data_<SpDPtr>::Scalar2Index( SizeT& st) const
2327 {
2328 throw GDLException("PTR expression not allowed in this context.");
2329 return 0; // get rid of warning
2330 }
2331 template<>
Scalar2RangeT(RangeT & st) const2332 int Data_<SpDPtr>::Scalar2RangeT( RangeT& st) const
2333 {
2334 throw GDLException("PTR expression not allowed in this context.");
2335 return 0; // get rid of warning
2336 }
2337
2338 template<>
Scalar2Index(SizeT & st) const2339 int Data_<SpDObj>::Scalar2Index( SizeT& st) const
2340 {
2341 throw GDLException("Object expression not allowed in this context.");
2342 return 0; // get rid of warning
2343 }
2344 template<>
Scalar2RangeT(RangeT & st) const2345 int Data_<SpDObj>::Scalar2RangeT( RangeT& st) const
2346 {
2347 throw GDLException("Object expression not allowed in this context.");
2348 return 0; // get rid of warning
2349 }
2350
2351
2352
2353
2354
2355
2356 // for FOR loop *indices*
2357 template<class Sp>
LoopIndex() const2358 RangeT Data_<Sp>::LoopIndex() const
2359 {
2360 // if( dd.size() != 1) return 0;
2361
2362 // the next statement gives a warning for unsigned integer types:
2363 // "comparison is always false due to limited range of data type"
2364 // This is because the same template is used here for signed and
2365 // unsigned data. A template specialization for the unsigned integer
2366 // types would result in three identical specializations, which is bad
2367 // for maintainability. And as any modern C++ compiler will optimize
2368 // away the superflous (for unsigned data) statement anyway, it is
2369 // better to keep the code this way here.
2370 // if( (*this)[0] < 0)
2371 // throw GDLException( "Index variable <0.");
2372
2373 return static_cast<RangeT>((*this)[0]);
2374 }
2375 template<>
LoopIndex() const2376 RangeT Data_<SpDFloat>::LoopIndex() const
2377 {
2378 // if( (*this)[0] < 0.0f)
2379 // //if( (*this)[0] <= 1.0f)
2380 // throw GDLException( "Index variable <0.");
2381 //else
2382 // return 0;
2383
2384 return static_cast<RangeT>((*this)[0]);
2385 }
2386 template<>
LoopIndex() const2387 RangeT Data_<SpDDouble>::LoopIndex() const
2388 {
2389 // if( (*this)[0] < 0.0)
2390 // //if( (*this)[0] <= 1.0)
2391 // throw GDLException( "Index variable <0.");
2392 //else
2393 // return 0;
2394
2395 return static_cast<RangeT>((*this)[0]);
2396 }
2397 template<>
LoopIndex() const2398 RangeT Data_<SpDComplex>::LoopIndex() const
2399 {
2400 throw GDLException( "Complex expression not allowed as index.");
2401 return 0;
2402 }
2403 template<>
LoopIndex() const2404 RangeT Data_<SpDComplexDbl>::LoopIndex() const
2405 {
2406 throw GDLException( "Complex expression not allowed as index.");
2407 return 0;
2408 }
2409 template<>
LoopIndex() const2410 RangeT Data_<SpDString>::LoopIndex() const
2411 {
2412 if( (*this)[0] == "")
2413 return 0;
2414
2415 const char* cStart=(*this)[0].c_str();
2416 char* cEnd;
2417 RangeT ix=strtol(cStart,&cEnd,10);
2418 if( cEnd == cStart)
2419 {
2420 Warning( "Type conversion error: "
2421 "Unable to convert given STRING: '"+
2422 (*this)[0]+"' to index.");
2423 return 0;
2424 }
2425 return ix;
2426 }
2427
2428 template<>
LoopIndex() const2429 RangeT Data_<SpDPtr>::LoopIndex() const
2430 {
2431 throw GDLException("PTR expression not allowed as index.");
2432 return 0; // get rid of warning
2433 }
2434
2435 template<>
LoopIndex() const2436 RangeT Data_<SpDObj>::LoopIndex() const
2437 {
2438 throw GDLException("Object expression not allowed as index.");
2439 return 0; // get rid of warning
2440 }
2441
2442
2443 // True
2444 template<class Sp>
True()2445 bool Data_<Sp>::True()
2446 {
2447 Ty s;
2448 if( !Scalar( s))
2449 throw GDLException("Expression must be a scalar or 1 element array in this context.",true,false);
2450 return (s % 2);
2451 }
2452
2453 template<>
True()2454 bool Data_<SpDFloat>::True()
2455 {
2456 Ty s;
2457 if( !Scalar( s))
2458 throw GDLException("Expression must be a scalar or 1 element array in this context.",true,false);
2459 return (s != 0.0f);
2460 }
2461
2462 template<>
True()2463 bool Data_<SpDDouble>::True()
2464 {
2465 Ty s;
2466 if( !Scalar( s))
2467 throw GDLException("Expression must be a scalar or 1 element array in this context.",true,false);
2468 return (s != 0.0);
2469 }
2470
2471 template<>
True()2472 bool Data_<SpDString>::True()
2473 {
2474 Ty s;
2475 if( !Scalar( s))
2476 throw GDLException("Expression must be a scalar or 1 element array in this context.",true,false);
2477 return (s != "");
2478 }
2479
2480 template<>
True()2481 bool Data_<SpDComplex>::True()
2482 {
2483 Ty s;
2484 if( !Scalar( s))
2485 throw GDLException("Expression must be a scalar or 1 element array in this context.",true,false);
2486 return (real(s) != 0.0 || imag(s) != 0.0);
2487 }
2488 template<>
True()2489 bool Data_<SpDComplexDbl>::True()
2490 {
2491 Ty s;
2492 if( !Scalar( s))
2493 throw GDLException("Expression must be a scalar or 1 element array in this context.",true,false);
2494 return (real(s) != 0.0 || imag(s) != 0.0);
2495 }
2496
2497 template<>
True()2498 bool Data_<SpDPtr>::True()
2499 {
2500 Ty s;
2501 if( !Scalar( s))
2502 throw GDLException("Expression must be a scalar or 1 element array in this context.",true,false);
2503 return (s != 0);
2504 }
2505
2506 template<>
True()2507 bool Data_<SpDObj>::True()
2508 {
2509 Ty s;
2510 if( !Scalar( s))
2511 throw GDLException("Expression must be a scalar or 1 element array in this context.",true,false);
2512 if( s == 0)
2513 return false; // on overloads for null object -> default is 'false'
2514
2515 // DStructGDL* oStructGDL= GDLInterpreter::GetObjHeapNoThrow( s);
2516 // if( oStructGDL == NULL) // object not valid -> default behaviour
2517 // return true; // true is ok here: Default behaviour is to just checks for null object
2518 //
2519 // DStructDesc* desc = oStructGDL->Desc();
2520 //
2521 // DFun* isTrueOverload = static_cast<DFun*>(desc->GetOperator( OOIsTrue));
2522
2523 DSubUD* isTrueOverload = static_cast<DSubUD*>(GDLInterpreter::GetObjHeapOperator( s, OOIsTrue));
2524 if( isTrueOverload == NULL)
2525 return true; // not overloaded, false case for default already returned (s. a.)
2526
2527 ProgNodeP callingNode = interpreter->GetRetTree();
2528
2529 // BaseGDL* self = this->Dup();
2530 // Guard<BaseGDL> selfGuard( self);
2531 // EnvUDT* newEnv= new EnvUDT( callingNode, isTrueOverload, &self);
2532 // no parameters
2533
2534 EnvUDT* newEnv;
2535 DObjGDL* self;
2536 Guard<BaseGDL> selfGuard;
2537 // Dup() here is not optimal
2538 // avoid at least for internal overload routines (which do/must not change SELF or r)
2539 bool internalDSubUD = isTrueOverload->GetTree()->IsWrappedNode();
2540 if( internalDSubUD)
2541 {
2542 self = this;
2543 newEnv= new EnvUDT( callingNode, isTrueOverload, &self);
2544 }
2545 else
2546 {
2547 self = this->Dup();
2548 selfGuard.Init( self);
2549 newEnv= new EnvUDT( callingNode, isTrueOverload, &self);
2550 }
2551
2552 StackGuard<EnvStackT> guard(interpreter->CallStack());
2553
2554 interpreter->CallStack().push_back( newEnv);
2555
2556 // make the call
2557 BaseGDL* res=interpreter->call_fun(static_cast<DSubUD*>(newEnv->GetPro())->GetTree());
2558
2559 if( !internalDSubUD && self != selfGuard.Get())
2560 {
2561 // always put out warning first, in case of a later crash
2562 Warning( "WARNING: " + isTrueOverload->ObjectName() +
2563 ": Assignment to SELF detected (GDL session still ok).");
2564 // assignment to SELF -> self was deleted and points to new variable
2565 // which it owns
2566 selfGuard.Release();
2567 if( (BaseGDL*)self != NullGDL::GetSingleInstance())
2568 selfGuard.Reset(self);
2569 }
2570 if( NullGDL::IsNULLorNullGDL( res))
2571 {
2572 throw GDLException( isTrueOverload->ObjectName() + " returned an undefined value.",true,false);
2573 }
2574
2575 Guard<BaseGDL> resGuard( res);
2576
2577 // prevent recursion
2578 if( res->Type() == GDL_OBJ)
2579 {
2580 std::ostringstream os;
2581 res->ToStream(os);
2582 throw GDLException( isTrueOverload->ObjectName() + ": Object reference expression not allowed in this context: " +
2583 os.str(),true,false);
2584 }
2585
2586 return res->LogTrue();
2587 }
2588
2589 // False
2590 template<class Sp>
False()2591 bool Data_<Sp>::False()
2592 {
2593 return !True();
2594 }
2595
2596 // Sgn
2597 template<class Sp>
Sgn()2598 int Data_<Sp>::Sgn() // -1,0,1
2599 {
2600 Ty s;
2601 if( !Scalar( s))
2602 throw GDLException("Expression must be a scalar or 1 element array in this context.",true,false);
2603 if( s > 0) return 1;
2604 if( s == 0) return 0;
2605 return -1;
2606 }
2607
2608 template<>
Sgn()2609 int Data_<SpDString>::Sgn() // -1,0,1
2610 {
2611 throw GDLException("String expression not allowed in this context.");
2612 return 0;
2613 }
2614 template<>
Sgn()2615 int Data_<SpDComplex>::Sgn() // -1,0,1
2616 {
2617 throw GDLException("Complex expression not allowed in this context.");
2618 return 0;
2619 }
2620 template<>
Sgn()2621 int Data_<SpDComplexDbl>::Sgn() // -1,0,1
2622 {
2623 throw GDLException("Complex expression not allowed in this context.");
2624 return 0;
2625 }
2626
2627 template<>
Sgn()2628 int Data_<SpDPtr>::Sgn() // -1,0,1
2629 {
2630 throw GDLException("Ptr expression not allowed in this context.");
2631 return 0;
2632 }
2633
2634 template<>
Sgn()2635 int Data_<SpDObj>::Sgn() // -1,0,1
2636 {
2637 throw GDLException("Object expression not allowed in this context.");
2638 return 0;
2639 }
2640
2641 // Equal (deletes r) only used in ForCheck(...)
2642 template<class Sp>
Equal(BaseGDL * r) const2643 bool Data_<Sp>::Equal( BaseGDL* r) const
2644 {
2645 assert( r->StrictScalar());
2646 // if( !r->Scalar())
2647 // {
2648 // GDLDelete(r);
2649 // throw GDLException("Expression must be a scalar in this context.");
2650 // }
2651 assert( r->Type() == this->t);
2652 Data_* rr=static_cast<Data_*>(r);
2653 // Data_* rr=static_cast<Data_*>(r->Convert2( this->t));
2654 bool ret= ((*this)[0] == (*rr)[0]);
2655 GDLDelete(rr);
2656 return ret;
2657 }
2658 template<>
Equal(BaseGDL * r) const2659 bool Data_<SpDFloat>::Equal( BaseGDL* r) const
2660 {
2661 assert( r->StrictScalar());
2662 // if( !r->Scalar())
2663 // {
2664 // GDLDelete(r);
2665 // throw GDLException("Expression must be a scalar in this context.");
2666 // }
2667 assert( r->Type() == this->t);
2668 Data_* rr=static_cast<Data_*>(r);
2669 // Data_* rr=static_cast<Data_*>(r->Convert2( this->t));
2670 bool ret= fabs((*this)[0] - (*rr)[0]) < 1.0f;
2671 GDLDelete(rr);
2672 return ret;
2673 }
2674 template<>
Equal(BaseGDL * r) const2675 bool Data_<SpDDouble>::Equal( BaseGDL* r) const
2676 {
2677 assert( r->StrictScalar());
2678 // if( !r->Scalar())
2679 // {
2680 // GDLDelete(r);
2681 // throw GDLException("Expression must be a scalar in this context.");
2682 // }
2683 assert( r->Type() == this->t);
2684 Data_* rr=static_cast<Data_*>(r);
2685 // Data_* rr=static_cast<Data_*>(r->Convert2( this->t));
2686 bool ret= fabs((*this)[0] - (*rr)[0]) < 1.0;
2687 GDLDelete(rr);
2688 return ret;
2689 }
2690
2691 // Equal (does not delete r)
2692 template<class Sp>
EqualNoDelete(const BaseGDL * r) const2693 bool Data_<Sp>::EqualNoDelete( const BaseGDL* r) const
2694 {
2695 if( !r->Scalar())
2696 {
2697 throw GDLException("Expression must be a scalar in this context.");
2698 }
2699 bool ret;
2700 if( r->Type() != this->t)
2701 {
2702 Data_* rr=static_cast<Data_*>(const_cast<BaseGDL*>(r)->Convert2( this->t, BaseGDL::COPY));
2703 ret= ((*this)[0] == (*rr)[0]);
2704 GDLDelete(rr);
2705 }
2706 else
2707 {
2708 const Data_* rr=static_cast<const Data_*>(r);
2709 ret= ((*this)[0] == (*rr)[0]);
2710 }
2711 return ret;
2712 }
2713
Equal(BaseGDL * r) const2714 bool DStructGDL::Equal( BaseGDL* r) const
2715 {
2716 GDLDelete(r);
2717 throw GDLException("Struct expression not allowed in this context.");
2718 return false;
2719 }
2720
2721 // For array_equal r must be of same type
2722 template<class Sp>
ArrayEqual(BaseGDL * rIn)2723 bool Data_<Sp>::ArrayEqual( BaseGDL* rIn) {
2724 Data_<Sp>* r = static_cast<Data_<Sp>*> (rIn);
2725 SizeT nEl = N_Elements();
2726 SizeT rEl = r->N_Elements();
2727 if (rEl == 1) {
2728 for (SizeT i = 0; i < nEl; ++i)
2729 if ((*this)[i] != (*r)[0]) return false;
2730 return true;
2731 }
2732 if (nEl == 1) {
2733 for (SizeT i = 0; i < rEl; ++i)
2734 if ((*this)[0] != (*r)[i]) return false;
2735 return true;
2736 }
2737 if (nEl != rEl) return false;
2738 for (SizeT i = 0; i < nEl; ++i)
2739 if ((*this)[i] != (*r)[i]) return false;
2740 return true;
2741 }
2742
ArrayEqual(BaseGDL * r)2743 bool DStructGDL::ArrayEqual( BaseGDL* r)
2744 {
2745 throw GDLException("Struct expression not allowed in this context.");
2746 return false;
2747 }
2748
2749 // For array_equal r must be of same type
2750 template<class Sp>
ArrayNeverEqual(BaseGDL * rIn)2751 bool Data_<Sp>::ArrayNeverEqual( BaseGDL* rIn)
2752 {
2753 Data_<Sp>* r = static_cast< Data_<Sp>*>( rIn);
2754 SizeT nEl = N_Elements();
2755 SizeT rEl = r->N_Elements();
2756 if( rEl == 1)
2757 {
2758 for( SizeT i=0; i<nEl; ++i)
2759 if( (*this)[i] == (*r)[0]) return false;
2760 return true;
2761 }
2762 if( nEl == 1)
2763 {
2764 for( SizeT i=0; i<rEl; ++i)
2765 if( (*this)[0] == (*r)[i]) return false;
2766 return true;
2767 }
2768 if( nEl != rEl) return true;
2769 for( SizeT i=0; i<nEl; ++i)
2770 if( (*this)[i] == (*r)[i]) return false;
2771 return true;
2772 }
2773
ArrayNeverEqual(BaseGDL * r)2774 bool DStructGDL::ArrayNeverEqual( BaseGDL* r)
2775 {
2776 throw GDLException("Struct expression not allowed in this context.");
2777 return false;
2778 }
2779
2780
2781 //Not used
2782 //template<class Sp>
2783 //bool Data_<Sp>::OutOfRangeOfInt() const
2784 //{
2785 // assert( this->StrictScalar());
2786 // return (*this)[0] > std::numeric_limits< DInt>::max() || (*this)[0] < std::numeric_limits< DInt>::min();
2787 //}
2788 //
2789 //template<>
2790 //bool Data_<SpDString>::OutOfRangeOfInt() const
2791 //{
2792 // return false;
2793 //}
2794 //template<>
2795 //bool Data_<SpDByte>::OutOfRangeOfInt() const
2796 //{
2797 // return false;
2798 //}
2799 //template<>
2800 //bool Data_<SpDComplex>::OutOfRangeOfInt() const
2801 //{
2802 // return false;
2803 //}
2804 //template<>
2805 //bool Data_<SpDComplexDbl>::OutOfRangeOfInt() const
2806 //{
2807 // return false;
2808 //}
2809
2810 // for statement compliance (int types , float types scalar only)
2811 // (convert strings to floats here (not for first argument)
2812
2813 template<class Sp>
ForCheck(BaseGDL ** lEnd,BaseGDL ** lStep)2814 bool Data_<Sp>::ForCheck(BaseGDL** lEnd, BaseGDL** lStep) {
2815 // all scalars?
2816 if (!StrictScalar())
2817 throw GDLException("Loop INIT must be a scalar in this context.");
2818
2819 if (!(*lEnd)->StrictScalar())
2820 throw GDLException("Loop LIMIT must be a scalar in this context.");
2821
2822 if (lStep != NULL && !(*lStep)->StrictScalar())
2823 throw GDLException("Loop INCREMENT must be a scalar in this context.");
2824
2825 // only proper types?
2826 if (this->t == GDL_UNDEF)
2827 throw GDLException("Expression is undefined.");
2828 if (this->t == GDL_COMPLEX || this->t == GDL_COMPLEXDBL)
2829 throw GDLException("Complex expression not allowed in this context.");
2830 if (this->t == GDL_PTR)
2831 throw GDLException("Pointer expression not allowed in this context.");
2832 if (this->t == GDL_OBJ)
2833 throw GDLException("Object expression not allowed in this context.");
2834 if (this->t == GDL_STRING)
2835 throw GDLException("String expression not allowed in this context.");
2836
2837 DType lType = (*lEnd)->Type();
2838 if (lType == GDL_COMPLEX || lType == GDL_COMPLEXDBL)
2839 throw GDLException("Complex expression not allowed in this context.");
2840
2841 //ForCheck() is donce only once so the penalty is minimal for these complicated tests.
2842 //BYTE is quite special and (see #816) there are dubious cases. The following code
2843 //avoids all traps at the expense of being complicated and promote to INT in a case where IDL does not:
2844 //"for i=255b,255,1 do help,i"
2845 if (this->t == GDL_BYTE) {
2846 //early exit with simple check is needed for Bytes beacuse the arithmetic is not the same as for other types.
2847 // simple check: end vs. start
2848 (*lEnd) = (*lEnd)->Convert2(GDL_BYTE); //convert to byte!
2849 DByte* endval = static_cast<DByte*> ((*lEnd)->DataAddr());
2850 DByte* startval = static_cast<DByte*> (this->DataAddr());
2851
2852 if (lStep != NULL) { //overflow case NOT handled in forAddCondUp
2853 (*lStep) = (*lStep)->Convert2(GDL_LONG); //check a few things with a long argument
2854 //now what if end+step > 255 ?
2855 DLong* step = static_cast<DLong*> ((*lStep)->DataAddr());
2856 if (std::signbit(step[0])) {
2857 // 1) check lEnd is not already smaller than lStart
2858 if (endval[0] > startval[0]) return false;
2859 // 2) We must convert to INTs, but this may be a problem.
2860 (*lEnd) = (*lEnd)->Convert2(GDL_INT);
2861 } else { //step >= 0
2862 if (endval[0] < startval[0]) return false;
2863 int final = endval[0] + step[0];
2864 if (final > 255) { //whatever the start value, the end value must be 255 when the final step is reached.
2865 //We must convert to INTs, but this may be a problem.
2866 (*lEnd) = (*lEnd)->Convert2(GDL_INT);
2867 }
2868 }
2869 //Must have lStep always as LEnd
2870 (*lStep) = (*lStep)->Convert2((*lEnd)->Type());
2871 } else if (endval[0] < startval[0]) return false; //lStep==NULL case
2872 return true;
2873 }
2874
2875 // Check for promotion as we have to be robust to lEnd==std::numeric_limits< DInt>::max() and lEnd+Step > numeric limit. (bug #816)
2876 // (*lend)+/-(*lstep) must not overflow an INT numeric limit, else the loop will fail.If (*lStep) is not defined,
2877 // forAddCondUp is used and will save the day as the test is made in the right order. this is not the general case, where the
2878 //addition (or substraction) is made in the prognode Run() loop. In that case we *NEED* to promote the loop variable correctly.
2879
2880 (*lEnd) = (*lEnd)->Convert2(GDL_LONG64); // upgrade to safe encoding
2881 DLong64* endval = static_cast<DLong64*> ((*lEnd)->DataAddr());
2882
2883 DLong64 testVal=endval[0];
2884
2885 if (lStep != NULL) (*lStep) = (*lStep)->Convert2(GDL_LONG64); //idem
2886 DLong64* step = NULL;
2887 if (lStep != NULL) {
2888 step = static_cast<DLong64*> ((*lStep)->DataAddr());
2889 testVal+=step[0];
2890 }
2891
2892 if (this->t == GDL_INT) {
2893 if (testVal < std::numeric_limits< DInt>::max() && testVal > std::numeric_limits< DInt>::min()) *lEnd = (*lEnd)->Convert2(GDL_INT);
2894 else if (testVal < std::numeric_limits< DLong>::max() && testVal > std::numeric_limits< DLong>::min()) *lEnd = (*lEnd)->Convert2(GDL_LONG);
2895 //Must have lStep always as LEnd
2896 if (lStep != NULL) (*lStep) = (*lStep)->Convert2((*lEnd)->Type());
2897 return true; // finished for GDL_INT
2898 }
2899 if (this->t == GDL_LONG) {
2900 if (testVal < std::numeric_limits< DLong>::max() && testVal > std::numeric_limits< DLong>::min()) *lEnd = (*lEnd)->Convert2(GDL_LONG);
2901 //Must have lStep always as LEnd
2902 if (lStep != NULL) (*lStep) = (*lStep)->Convert2((*lEnd)->Type());
2903 return true; // finished for GDL_LONG
2904 }
2905 // other cases: no promotion
2906 *lEnd = (*lEnd)->Convert2(this->t);
2907 if (lStep != NULL) *lStep = (*lStep)->Convert2(this->t);
2908 return true;
2909 }
2910
ForCheck(BaseGDL ** lEnd,BaseGDL ** lStep)2911 bool DStructGDL::ForCheck( BaseGDL** lEnd, BaseGDL** lStep)
2912 {
2913 throw GDLException("Struct expression not allowed in this context.");
2914 return false;
2915 }
2916
2917 // ForCheck must have been called before
2918 template<class Sp>
ForAddCondUp(BaseGDL * endLoopVar)2919 bool Data_<Sp>::ForAddCondUp( BaseGDL* endLoopVar)
2920 {
2921 if( endLoopVar->Type() != this->t) throw GDLException("Type of FOR index variable changed.");
2922 Data_* lEnd=static_cast<Data_*>(endLoopVar);
2923 bool what=true;
2924 if ((*this)[0] == (*lEnd)[0]) //This way, loop will stop on good end value and loop index will be incremented...
2925 //but only afterwards as to avoid bug #816
2926 {
2927 what=false;
2928 } else what=((*this)[0] < (*lEnd)[0]);
2929 (*this)[0] += 1;
2930 return what;
2931 }
2932 // ForCheck must have been called before
2933 template<class Sp>
ForCondUp(BaseGDL * lEndIn)2934 bool Data_<Sp>::ForCondUp( BaseGDL* lEndIn)
2935 {
2936 if( lEndIn->Type() != this->t) throw GDLException("Type of FOR index variable changed.");
2937 Data_* lEnd=static_cast<Data_*>(lEndIn);
2938 return (*this)[0] <= (*lEnd)[0];
2939 }
2940 template<class Sp>
ForCondDown(BaseGDL * lEndIn)2941 bool Data_<Sp>::ForCondDown( BaseGDL* lEndIn)
2942 {
2943 if( lEndIn->Type() != this->t) throw GDLException("Type of FOR index variable changed.");
2944 Data_* lEnd=static_cast<Data_*>(lEndIn);
2945 return (*this)[0] >= (*lEnd)[0];
2946 }
2947
2948 // error if the type of the loop variable changed
ForCondUp(BaseGDL *)2949 bool DStructGDL::ForCondUp( BaseGDL*)
2950 {
2951 throw GDLException("Type of FOR index variable changed to STRUCT.");
2952 return false;
2953 }
ForCondDown(BaseGDL *)2954 bool DStructGDL::ForCondDown( BaseGDL*)
2955 {
2956 throw GDLException("Type of FOR index variable changed to STRUCT.");
2957 return false;
2958 }
2959 template<>
ForCondUp(BaseGDL *)2960 bool Data_<SpDComplex>::ForCondUp( BaseGDL*)
2961 {
2962 throw GDLException("Type of FOR index variable changed to COMPLEX.");
2963 return false;
2964 }
2965
2966 template<>
ForAddCondUp(BaseGDL * loopInfo)2967 bool Data_<SpDComplex>::ForAddCondUp( BaseGDL* loopInfo)
2968 // bool Data_<SpDComplex>::ForAddCondUp( ForLoopInfoT& loopInfo)
2969 {
2970 throw GDLException("Type of FOR index variable changed to COMPLEX.");
2971 return false;
2972 }
2973 template<>
ForCondDown(BaseGDL *)2974 bool Data_<SpDComplex>::ForCondDown( BaseGDL*)
2975 {
2976 throw GDLException("Type of FOR index variable changed to COMPLEX.");
2977 return false;
2978 }
2979
2980 template<>
ForAddCondUp(BaseGDL * loopInfo)2981 bool Data_<SpDComplexDbl>::ForAddCondUp( BaseGDL* loopInfo)
2982 // bool Data_<SpDComplexDbl>::ForAddCondUp( ForLoopInfoT& loopInfo)
2983 {
2984 throw GDLException("Type of FOR index variable changed to DCOMPLEX.");
2985 return false;
2986 }
2987 template<>
ForCondUp(BaseGDL *)2988 bool Data_<SpDComplexDbl>::ForCondUp( BaseGDL*)
2989 {
2990 throw GDLException("Type of FOR index variable changed to DCOMPLEX.");
2991 return false;
2992 }
2993 template<>
ForCondDown(BaseGDL *)2994 bool Data_<SpDComplexDbl>::ForCondDown( BaseGDL*)
2995 {
2996 throw GDLException("Type of FOR index variable changed to DCOMPLEX.");
2997 return false;
2998 }
2999
3000 // ForCheck must have been called before
3001 // general version
3002 template<class Sp>
ForAdd(BaseGDL * addIn)3003 void Data_<Sp>::ForAdd( BaseGDL* addIn)
3004 {
3005 if( addIn == NULL)
3006 {
3007 (*this)[0] += 1;
3008 return;
3009 }
3010 Data_* add=static_cast<Data_*>(addIn);
3011 (*this)[0] += (*add)[0];
3012 }
3013 // cannnot be called, just to make the compiler shut-up
ForAdd(BaseGDL * addIn)3014 void DStructGDL::ForAdd( BaseGDL* addIn) {}
3015
3016
3017 //NOT USED (GD)
3018 //// normal (+1) version
3019 //template<class Sp>
3020 //void Data_<Sp>::ForAdd()
3021 //{
3022 // (*this)[0] += 1;
3023 //}
3024 // cannnot be called, just to make the compiler shut-up
3025 //void DStructGDL::ForAdd() {}
3026
3027 template<class Sp>
AssignAtIx(RangeT ixR,BaseGDL * srcIn)3028 void Data_<Sp>::AssignAtIx( RangeT ixR, BaseGDL* srcIn)
3029 {
3030 if( ixR < 0)
3031 {
3032 SizeT nEl = this->N_Elements();
3033
3034 if( -ixR > nEl)
3035 throw GDLException("Subscript out of range: " + i2s(ixR));
3036
3037 SizeT ix = nEl + ixR;
3038
3039 if( srcIn->Type() != this->Type())
3040 {
3041 Data_* rConv = static_cast<Data_*>(srcIn->Convert2( this->Type(), BaseGDL::COPY_BYTE_AS_INT));
3042 // Data_* rConv = static_cast<Data_*>(srcIn->Convert2( this->Type(), BaseGDL::COPY));
3043 Guard<Data_> conv_guard( rConv);
3044 (*this)[ix] = (*rConv)[0];
3045 }
3046 else
3047 (*this)[ix] = (*static_cast<Data_*>(srcIn))[0];
3048
3049 return;
3050 } // ixR >= 0
3051 if( srcIn->Type() != this->Type())
3052 {
3053 Data_* rConv = static_cast<Data_*>(srcIn->Convert2( this->Type(), BaseGDL::COPY_BYTE_AS_INT));
3054 // Data_* rConv = static_cast<Data_*>(srcIn->Convert2( this->Type(), BaseGDL::COPY));
3055 Guard<Data_> conv_guard( rConv);
3056 (*this)[ixR] = (*rConv)[0];
3057 }
3058 else
3059 (*this)[ixR] = (*static_cast<Data_*>(srcIn))[0];
3060 }
3061
3062 // assigns srcIn to this at ixList, if ixList is NULL does linear copy
3063 // assumes: ixList has this already set as variable
3064 // used by DotAccessDescT::DoAssign
3065 // GDLInterpreter::l_array_expr
3066 template<class Sp>
AssignAt(BaseGDL * srcIn,ArrayIndexListT * ixList,SizeT offset)3067 void Data_<Sp>::AssignAt( BaseGDL* srcIn, ArrayIndexListT* ixList, SizeT offset)
3068 {
3069 // breakpoint(); // gdbg can not handle breakpoints in template functions
3070 Data_* src = static_cast<Data_*>(srcIn);
3071
3072 SizeT srcElem= src->N_Elements();
3073 // bool isScalar= (srcElem == 1);
3074 bool isScalar= (srcElem == 1) && (src->Rank() == 0);
3075 if( isScalar)
3076 { // src is scalar
3077 Ty scalar=(*src)[0];
3078
3079 if( ixList == NULL)
3080 {
3081 SizeT nCp=Data_::N_Elements();
3082
3083 /*#pragma omp parallel if (nCp >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nCp))
3084 {
3085 #pragma omp for*/
3086 for( int c=0; c<nCp; ++c)
3087 (*this)[ c]=scalar;
3088 // }
3089 }
3090 else
3091 {
3092 SizeT nCp=ixList->N_Elements();
3093
3094 AllIxBaseT* allIx = ixList->BuildIx();
3095 /*#pragma omp parallel if (nCp >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nCp))
3096 {
3097 #pragma omp for*/
3098 (*this)[ allIx->InitSeqAccess()]=scalar;
3099 for( SizeT c=1; c<nCp; ++c)
3100 (*this)[ allIx->SeqAccess()]=scalar;
3101 //} // (*this)[ ixList->GetIx( c)]=scalar;
3102 }
3103 }
3104 else
3105 {
3106 if( ixList == NULL)
3107 {
3108 SizeT nCp=Data_::N_Elements();
3109
3110 // if (non-indexed) src is smaller -> just copy its number of elements
3111 if( nCp > (srcElem-offset)) {
3112 if( offset == 0)
3113 nCp=srcElem;
3114 else
3115 throw GDLException("Source expression contains not enough elements.");
3116 }
3117 /*#pragma omp parallel if (nCp >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nCp))
3118 {
3119 #pragma omp for*/
3120 for( int c=0; c<nCp; ++c)
3121 (*this)[ c]=(*src)[c+offset];
3122 }//}
3123 else
3124 {
3125 // crucial part
3126 SizeT nCp=ixList->N_Elements();
3127
3128 if( nCp == 1)
3129 {
3130 SizeT destStart = ixList->LongIx();
3131 // len = 1;
3132 SizeT rStride = srcIn->Stride(this->Rank());
3133 (*this)[ destStart] = (*src)[ offset/rStride];
3134
3135 // InsAt( src, ixList, offset);
3136 }
3137 else
3138 {
3139 if( offset == 0)
3140 {
3141 if( srcElem < nCp)
3142 throw GDLException("Array subscript must have same size as"
3143 " source expression.");
3144
3145 AllIxBaseT* allIx = ixList->BuildIx();
3146 /*#pragma omp parallel if (nCp >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nCp))
3147 {
3148 #pragma omp for*/
3149 (*this)[ allIx->InitSeqAccess()]=(*src)[0];
3150 for( SizeT c=1; c<nCp; ++c)
3151 (*this)[ allIx->SeqAccess()]=(*src)[c];
3152 // } // (*this)[ ixList->GetIx( c)]=(*src)[c+offset];
3153 }
3154 else
3155 {
3156 if( (srcElem-offset) < nCp)
3157 throw GDLException("Array subscript must have same size as"
3158 " source expression.");
3159
3160 AllIxBaseT* allIx = ixList->BuildIx();
3161 /*#pragma omp parallel if (nCp >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nCp))
3162 {
3163 #pragma omp for*/
3164 (*this)[ allIx->InitSeqAccess()]=(*src)[offset];
3165 for( SizeT c=1; c<nCp; ++c)
3166 (*this)[ allIx->SeqAccess()]=(*src)[c+offset];
3167 // } // (*this)[ ixList->GetIx( c)]=(*src)[c+offset];
3168 }
3169 }
3170 }
3171 }
3172 }
3173 template<class Sp>
AssignAt(BaseGDL * srcIn,ArrayIndexListT * ixList)3174 void Data_<Sp>::AssignAt( BaseGDL* srcIn, ArrayIndexListT* ixList)
3175 {
3176 assert( ixList != NULL);
3177
3178 // breakpoint(); // gdbg can not handle breakpoints in template functions
3179 Data_* src = static_cast<Data_*>(srcIn);
3180
3181 SizeT srcElem= src->N_Elements();
3182 bool isScalar= (srcElem == 1);
3183 if( isScalar)
3184 { // src is scalar
3185 SizeT nCp=ixList->N_Elements();
3186
3187 if( nCp == 1)
3188 {
3189 (*this)[ ixList->LongIx()] = (*src)[0];
3190 }
3191 else
3192 {
3193 Ty scalar=(*src)[0];
3194 AllIxBaseT* allIx = ixList->BuildIx();
3195 (*this)[ allIx->InitSeqAccess()]=scalar;
3196 /*#pragma omp parallel if (nCp >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nCp))
3197 {
3198 #pragma omp for*/
3199 for( int c=1; c<nCp; ++c)
3200 (*this)[ allIx->SeqAccess()]=scalar;
3201 // (*this)[ (*allIx)[ c]]=scalar;
3202
3203
3204 // } // (*this)[ ixList->GetIx( c)]=scalar;
3205 }
3206 }
3207 else
3208 {
3209 // crucial part
3210 SizeT nCp=ixList->N_Elements();
3211
3212 if( nCp == 1)
3213 {
3214 InsAt( src, ixList);
3215 }
3216 else
3217 {
3218 if( srcElem < nCp)
3219 throw GDLException("Array subscript must have same size as"
3220 " source expression.");
3221
3222 AllIxBaseT* allIx = ixList->BuildIx();
3223 (*this)[ allIx->InitSeqAccess()]=(*src)[0];
3224 /*#pragma omp parallel if (nCp >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nCp))
3225 {
3226 #pragma omp for*/
3227 for( int c=1; c<nCp; ++c)
3228 (*this)[ allIx->SeqAccess()]=(*src)[c];
3229 // (*this)[ (*allIx)[ c]]=(*src)[c];
3230 // } // (*this)[ ixList->GetIx( c)]=(*src)[c+offset];
3231 }
3232 }
3233 }
3234 template<class Sp>
AssignAt(BaseGDL * srcIn)3235 void Data_<Sp>::AssignAt( BaseGDL* srcIn)
3236 {
3237 // breakpoint(); // gdbg can not handle breakpoints in template functions
3238 Data_* src = static_cast<Data_*>(srcIn);
3239
3240 SizeT srcElem= src->N_Elements();
3241 bool isScalar= (srcElem == 1);
3242 if( isScalar)
3243 { // src is scalar
3244 Ty scalar=(*src)[0];
3245
3246 /* dd = scalar;*/
3247 SizeT nCp=Data_::N_Elements();
3248
3249
3250 /*#pragma omp parallel if (nCp >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nCp))
3251 {
3252 #pragma omp for*/
3253 for( int c=0; c<nCp; ++c)
3254 (*this)[ c]=scalar;
3255 // }
3256 // SizeT nCp=Data_::N_Elements();
3257
3258 // for( SizeT c=0; c<nCp; ++c)
3259 // (*this)[ c]=scalar;
3260 }
3261 else
3262 {
3263 SizeT nCp=Data_::N_Elements();
3264
3265 // if (non-indexed) src is smaller -> just copy its number of elements
3266 if( nCp > srcElem) nCp=srcElem;
3267
3268 /*#pragma omp parallel if (nCp >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nCp))
3269 {
3270 #pragma omp for*/
3271 for( int c=0; c<nCp; ++c)
3272 (*this)[ c]=(*src)[c];
3273 // }
3274 }
3275 }
3276
3277 // increment/decrement operators
3278 // integers
3279 template<class Sp>
DecAt(ArrayIndexListT * ixList)3280 void Data_<Sp>::DecAt( ArrayIndexListT* ixList)
3281 {
3282 if( ixList == NULL)
3283 {
3284 dd -= 1;
3285 // SizeT nCp=Data_::N_Elements();
3286
3287 // for( SizeT c=0; c<nCp; ++c)
3288 // (*this)[ c]--;
3289 }
3290 else
3291 {
3292 SizeT nCp=ixList->N_Elements();
3293
3294 AllIxBaseT* allIx = ixList->BuildIx();
3295 /*#pragma omp parallel if (nCp >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nCp))
3296 {
3297 #pragma omp for*/
3298 (*this)[ allIx->InitSeqAccess()]--;
3299 for( SizeT c=1; c<nCp; ++c)
3300 (*this)[ allIx->SeqAccess()]--;
3301 }// }
3302 }
3303 template<class Sp>
IncAt(ArrayIndexListT * ixList)3304 void Data_<Sp>::IncAt( ArrayIndexListT* ixList)
3305 {
3306 if( ixList == NULL)
3307 {
3308 dd += 1;
3309 // SizeT nCp=Data_::N_Elements();
3310
3311 // for( SizeT c=0; c<nCp; ++c)
3312 // (*this)[ c]++;
3313 }
3314 else
3315 {
3316 SizeT nCp=ixList->N_Elements();
3317
3318 AllIxBaseT* allIx = ixList->BuildIx();
3319 /*#pragma omp parallel if (nCp >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nCp))
3320 {
3321 #pragma omp for*/
3322 (*this)[ allIx->InitSeqAccess()]++;
3323 for( SizeT c=1; c<nCp; ++c)
3324 (*this)[ allIx->SeqAccess()]++;
3325 }// }
3326 }
3327 // float, double
3328 template<>
DecAt(ArrayIndexListT * ixList)3329 void Data_<SpDFloat>::DecAt( ArrayIndexListT* ixList)
3330 {
3331 if( ixList == NULL)
3332 {
3333 dd -= 1.0f;
3334
3335 // SizeT nCp=Data_::N_Elements();
3336
3337 // for( SizeT c=0; c<nCp; ++c)
3338 // (*this)[ c] -= 1.0;
3339 }
3340 else
3341 {
3342 SizeT nCp=ixList->N_Elements();
3343
3344 AllIxBaseT* allIx = ixList->BuildIx();
3345 /*#pragma omp parallel if (nCp >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nCp))
3346 {
3347 #pragma omp for*/
3348 (*this)[ allIx->InitSeqAccess()] -= 1.0f;
3349 for( SizeT c=1; c<nCp; ++c)
3350 (*this)[ allIx->SeqAccess()] -=1.0f;
3351 }// }
3352 }
3353 template<>
IncAt(ArrayIndexListT * ixList)3354 void Data_<SpDFloat>::IncAt( ArrayIndexListT* ixList)
3355 {
3356 if( ixList == NULL)
3357 {
3358 dd += 1.0f;
3359
3360 // SizeT nCp=Data_::N_Elements();
3361
3362 // for( SizeT c=0; c<nCp; ++c)
3363 // (*this)[ c] += 1.0;
3364 }
3365 else
3366 {
3367 SizeT nCp=ixList->N_Elements();
3368
3369 AllIxBaseT* allIx = ixList->BuildIx();
3370 /*#pragma omp parallel if (nCp >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nCp))
3371 {
3372 #pragma omp for*/
3373 (*this)[ allIx->InitSeqAccess()] += 1.0f;
3374 for( SizeT c=1; c<nCp; ++c)
3375 (*this)[ allIx->SeqAccess()] +=1.0f;
3376 }// }
3377 }
3378 template<>
DecAt(ArrayIndexListT * ixList)3379 void Data_<SpDDouble>::DecAt( ArrayIndexListT* ixList)
3380 {
3381 if( ixList == NULL)
3382 {
3383 dd -= 1.0;
3384
3385 // SizeT nCp=Data_::N_Elements();
3386
3387 // for( SizeT c=0; c<nCp; ++c)
3388 // (*this)[ c] -= 1.0;
3389 }
3390 else
3391 {
3392 SizeT nCp=ixList->N_Elements();
3393
3394 AllIxBaseT* allIx = ixList->BuildIx();
3395 /*#pragma omp parallel if (nCp >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nCp))
3396 {
3397 #pragma omp for*/
3398 (*this)[ allIx->InitSeqAccess()] -= 1.0;
3399 for( SizeT c=1; c<nCp; ++c)
3400 (*this)[ allIx->SeqAccess()] -=1.0;
3401 }// }
3402 }
3403 template<>
IncAt(ArrayIndexListT * ixList)3404 void Data_<SpDDouble>::IncAt( ArrayIndexListT* ixList)
3405 {
3406 if( ixList == NULL)
3407 {
3408 dd += 1.0;
3409
3410 // SizeT nCp=Data_::N_Elements();
3411
3412 // for( SizeT c=0; c<nCp; ++c)
3413 // (*this)[ c] += 1.0;
3414 }
3415 else
3416 {
3417 SizeT nCp=ixList->N_Elements();
3418
3419 AllIxBaseT* allIx = ixList->BuildIx();
3420 /*#pragma omp parallel if (nCp >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nCp))
3421 {
3422 #pragma omp for*/
3423 (*this)[ allIx->InitSeqAccess()] += 1.0f;
3424 for( SizeT c=1; c<nCp; ++c)
3425 (*this)[ allIx->SeqAccess()] +=1.0f;
3426 }// }
3427 }
3428 // complex
3429 template<>
DecAt(ArrayIndexListT * ixList)3430 void Data_<SpDComplex>::DecAt( ArrayIndexListT* ixList)
3431 {
3432 if( ixList == NULL)
3433 {
3434 // dd -= 1.0f;
3435
3436 SizeT nCp=Data_::N_Elements();
3437
3438 /*#pragma omp parallel if (nCp >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nCp))
3439 {
3440 #pragma omp for*/
3441 for( int c=0; c<nCp; ++c)
3442 (*this)[ c] -= 1.0;
3443 }// }
3444 else
3445 {
3446 SizeT nCp=ixList->N_Elements();
3447
3448 AllIxBaseT* allIx = ixList->BuildIx();
3449 /*#pragma omp parallel if (nCp >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nCp))
3450 {
3451 #pragma omp for*/
3452 (*this)[ allIx->InitSeqAccess()] -= 1.0;
3453 for( SizeT c=1; c<nCp; ++c)
3454 (*this)[ allIx->SeqAccess()] -=1.0;
3455 // for( SizeT c=0; c<nCp; ++c)
3456 // (*this)[ (*allIx)[ c]] -= 1.0;
3457 }// }
3458 }
3459 template<>
IncAt(ArrayIndexListT * ixList)3460 void Data_<SpDComplex>::IncAt( ArrayIndexListT* ixList)
3461 {
3462 if( ixList == NULL)
3463 {
3464 // dd += 1.0f;
3465
3466 SizeT nCp=Data_::N_Elements();
3467
3468 /*#pragma omp parallel if (nCp >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nCp))
3469 {
3470 #pragma omp for*/
3471 for( int c=0; c<nCp; ++c)
3472 (*this)[ c] += 1.0;
3473 }// }
3474 else
3475 {
3476 SizeT nCp=ixList->N_Elements();
3477
3478 AllIxBaseT* allIx = ixList->BuildIx();
3479 /*#pragma omp parallel if (nCp >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nCp))
3480 {
3481 #pragma omp for*/
3482 (*this)[ allIx->InitSeqAccess()] += 1.0;
3483 for( SizeT c=1; c<nCp; ++c)
3484 (*this)[ allIx->SeqAccess()] +=1.0;
3485 // for( SizeT c=0; c<nCp; ++c)
3486 // (*this)[ (*allIx)[ c]] += 1.0;
3487 }// }
3488 }
3489 template<>
DecAt(ArrayIndexListT * ixList)3490 void Data_<SpDComplexDbl>::DecAt( ArrayIndexListT* ixList)
3491 {
3492 if( ixList == NULL)
3493 {
3494 // dd -= 1.0;
3495
3496 SizeT nCp=Data_::N_Elements();
3497
3498 /*#pragma omp parallel if (nCp >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nCp))
3499 {
3500 #pragma omp for*/
3501 for( int c=0; c<nCp; ++c)
3502 (*this)[ c] -= 1.0;
3503 }// }
3504 else
3505 {
3506 SizeT nCp=ixList->N_Elements();
3507
3508 AllIxBaseT* allIx = ixList->BuildIx();
3509 /*#pragma omp parallel if (nCp >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nCp))
3510 {
3511 #pragma omp for*/
3512 (*this)[ allIx->InitSeqAccess()] -= 1.0;
3513 for( SizeT c=1; c<nCp; ++c)
3514 (*this)[ allIx->SeqAccess()] -=1.0;
3515 // for( SizeT c=0; c<nCp; ++c)
3516 // (*this)[ (*allIx)[ c]] -= 1.0;
3517 } // }
3518 }
3519 template<>
IncAt(ArrayIndexListT * ixList)3520 void Data_<SpDComplexDbl>::IncAt( ArrayIndexListT* ixList)
3521 {
3522 if( ixList == NULL)
3523 {
3524 // dd += 1.0;
3525
3526 SizeT nCp=Data_::N_Elements();
3527
3528 /*#pragma omp parallel if (nCp >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nCp))
3529 {
3530 #pragma omp for*/
3531 for( int c=0; c<nCp; ++c)
3532 (*this)[ c] += 1.0;
3533 }// }
3534 else
3535 {
3536 SizeT nCp=ixList->N_Elements();
3537
3538 AllIxBaseT* allIx = ixList->BuildIx();
3539 /*#pragma omp parallel if (nCp >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nCp))
3540 {
3541 #pragma omp for*/
3542 (*this)[ allIx->InitSeqAccess()] += 1.0;
3543 for( SizeT c=1; c<nCp; ++c)
3544 (*this)[ allIx->SeqAccess()] +=1.0;
3545 // for( SizeT c=0; c<nCp; ++c)
3546 // (*this)[ (*allIx)[ c]] += 1.0;
3547 }// }
3548 }
3549 // forbidden types
3550 template<>
DecAt(ArrayIndexListT * ixList)3551 void Data_<SpDString>::DecAt( ArrayIndexListT* ixList)
3552 {
3553 throw GDLException("String expression not allowed in this context.");
3554 }
3555 template<>
DecAt(ArrayIndexListT * ixList)3556 void Data_<SpDPtr>::DecAt( ArrayIndexListT* ixList)
3557 {
3558 throw GDLException("Pointer expression not allowed in this context.");
3559 }
3560 template<>
DecAt(ArrayIndexListT * ixList)3561 void Data_<SpDObj>::DecAt( ArrayIndexListT* ixList)
3562 {
3563 throw GDLException("Object expression not allowed in this context.");
3564 }
3565 template<>
IncAt(ArrayIndexListT * ixList)3566 void Data_<SpDString>::IncAt( ArrayIndexListT* ixList)
3567 {
3568 throw GDLException("String expression not allowed in this context.");
3569 }
3570 template<>
IncAt(ArrayIndexListT * ixList)3571 void Data_<SpDPtr>::IncAt( ArrayIndexListT* ixList)
3572 {
3573 throw GDLException("Pointer expression not allowed in this context.");
3574 }
3575 template<>
IncAt(ArrayIndexListT * ixList)3576 void Data_<SpDObj>::IncAt( ArrayIndexListT* ixList)
3577 {
3578 throw GDLException("Object expression not allowed in this context.");
3579 }
3580
3581
3582 // used by AccessDescT for resolving, no checking is done
3583 // inserts srcIn[ ixList] at offset
3584 // used by DotAccessDescT::DoResolve
3585 template<class Sp>
InsertAt(SizeT offset,BaseGDL * srcIn,ArrayIndexListT * ixList)3586 void Data_<Sp>::InsertAt( SizeT offset, BaseGDL* srcIn,
3587 ArrayIndexListT* ixList)
3588 {
3589 Data_* src=static_cast<Data_* >(srcIn);
3590 if( ixList == NULL)
3591 {
3592 SizeT nCp=src->N_Elements();
3593
3594 /*#pragma omp parallel if (nCp >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nCp))
3595 {
3596 #pragma omp for*/
3597 for( int c=0; c<nCp; ++c)
3598 (*this)[ c+offset]=(*src)[c];
3599 }// }
3600 else
3601 {
3602 SizeT nCp=ixList->N_Elements();
3603
3604 AllIxBaseT* allIx = ixList->BuildIx();
3605 /*#pragma omp parallel if (nCp >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nCp))
3606 {
3607 #pragma omp for*/
3608 (*this)[ offset]=(*src)[ allIx->InitSeqAccess()];
3609 for( SizeT c=1; c<nCp; ++c)
3610 (*this)[ c+offset]=(*src)[ allIx->SeqAccess()];
3611 //} // (*this)[ c+offset]=(*src)[ ixList->GetIx( c)];
3612 }
3613 }
3614
3615
3616 // used for array concatenation
3617 template<class Sp>
CatArray(ExprListT & exprList,const SizeT catRankIx,const SizeT rank)3618 Data_<Sp>* Data_<Sp>::CatArray( ExprListT& exprList,
3619 const SizeT catRankIx,
3620 const SizeT rank)
3621 {
3622 // breakpoint();
3623 SizeT rankIx = RankIx( rank);
3624 SizeT maxIx = (catRankIx > rankIx)? catRankIx : rankIx;
3625
3626 dimension catArrDim(this->dim); // list contains at least one element
3627
3628 catArrDim.MakeRank( maxIx+1);
3629 catArrDim.SetOneDim(catRankIx,0); // clear rank which is added up
3630
3631 SizeT dimSum=0;
3632 ExprListIterT i=exprList.begin();
3633 for(; i != exprList.end(); ++i)
3634 {
3635 // conversion done already here to throw if type is Assoc_<>
3636 (*i)=(*i)->Convert2( this->t);
3637
3638 for( SizeT dIx=0; dIx<=maxIx; dIx++)
3639 {
3640 if( dIx != catRankIx)
3641 {
3642 if( catArrDim[dIx] == (*i)->Dim(dIx)) continue;
3643 if( (catArrDim[dIx] > 1) || ((*i)->Dim(dIx) > 1))
3644 throw GDLException("Unable to concatenate variables "
3645 "because the dimensions do not agree");
3646 }
3647 else
3648 {
3649 SizeT add=(*i)->Dim(dIx);
3650 dimSum+=(add)?add:1;
3651 }
3652 }
3653 }
3654
3655 catArrDim.SetOneDim(catRankIx,dimSum);
3656
3657 // special case: identity (c=[a]) equivalent to c=a, e.g. Dup(). Speedup
3658 if (exprList.size()==1 && (this->dim==catArrDim)) return this->Dup();
3659
3660 // the concatenated array
3661 Data_<Sp>* catArr=New(catArrDim, BaseGDL::NOZERO);
3662
3663 SizeT at=0;
3664
3665 // accelerated in CatInsert. Note that "at" is computed serially in CatInsert but can be computed externally (if one wanted to precompute it for parallel calls to CatInsert:
3666 // SizeT add=static_cast<Data_<Sp>*>( (*i))->dim[catRankIx]; // update 'at'
3667 // for( i=exprList.begin(); i != exprList.end(); ++i) at += (add > 1)? add : 1;
3668
3669 for( i=exprList.begin(); i != exprList.end(); ++i)
3670 {
3671 catArr->CatInsert(static_cast<Data_<Sp>*>( (*i)),
3672 catRankIx,at); // advances 'at'
3673 }
3674
3675 return catArr;
3676 }
3677
3678 // returns (*this)[ ixList]
3679 template<class Sp>
Index(ArrayIndexListT * ixList)3680 Data_<Sp>* Data_<Sp>::Index( ArrayIndexListT* ixList)
3681 {
3682 // ixList->SetVariable( this);
3683
3684 Data_* res=Data_::New( ixList->GetDim(), BaseGDL::NOZERO);
3685
3686 SizeT nCp=ixList->N_Elements();
3687
3688 // cout << "nCP = " << nCp << endl;
3689 // cout << "dim = " << this->dim << endl;
3690
3691 // DataT& res_dd = res->dd;
3692 AllIxBaseT* allIx = ixList->BuildIx();
3693
3694 if( nCp == 1)
3695 {
3696 (*res)[0]=(*this)[ (*allIx)[ 0]];
3697 return res;
3698 }
3699 // else
3700 // {
3701 /*#pragma omp parallel if (nCp >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nCp))
3702 {
3703 #pragma omp for*/
3704 (*res)[0]=(*this)[ allIx->InitSeqAccess()];
3705 for( SizeT c=1; c<nCp; ++c)
3706 (*res)[c]=(*this)[ allIx->SeqAccess()];
3707 //} // res_(*this)[c]=(*this)[ (*allIx)[ c]];
3708 // (*res)[c]=(*this)[ ixList->GetIx(c)];
3709 // }
3710 return res;
3711 }
3712
3713 // inserts srcIn at index ixDim
3714 // respects the exact structure of srcIn
3715 template<class Sp>
InsAt(Data_ * srcIn,ArrayIndexListT * ixList,SizeT offset)3716 void Data_<Sp>::InsAt( Data_* srcIn, ArrayIndexListT* ixList, SizeT offset)
3717 {
3718 // max. number of dimensions to copy
3719 SizeT nDim = ixList->NDim();
3720
3721 if( nDim == 1)
3722 // {
3723 // SizeT destStart = ixList->LongIx();
3724
3725 // SizeT len = srcIn->Dim( 0); // length of segment to copy
3726 // // check if in bounds of a
3727 // if( (destStart+len) > this->N_Elements()) //dim[0])
3728 // throw GDLException("Out of range subscript encountered (1).");
3729
3730 // DataT& srcIn_dd = srcIn->dd;
3731 // SizeT srcIx = 0; // this one simply runs from 0 to N_Elements(srcIn)
3732
3733 // SizeT destEnd = destStart + len;
3734 // for( SizeT destIx = destStart; destIx < destEnd; ++destIx)
3735 // (*this)[ destIx] = (*srcIn)[ srcIx++];
3736
3737 // return;
3738 // }
3739 {
3740 SizeT destStart = ixList->LongIx();
3741
3742 //SizeT len;
3743 if( this->N_Elements() == 1)
3744 {
3745 // len = 1;
3746 SizeT rStride = srcIn->Stride(this->Rank());
3747 (*this)[ destStart] = (*srcIn)[ offset/rStride];
3748 }
3749 else
3750 {
3751 SizeT len = srcIn->Dim( 0); // length of segment to copy
3752 // TODO: IDL reports here (and probably in the insert-dimension case below as well)
3753 // the name of a variable, e.g.:
3754 // IDL> a=[0,0,0] & a[2]=[2,2,2]
3755 // % Out of range subscript encountered: A.
3756 if( (destStart+len) > this->N_Elements()) //dim[0])
3757 throw GDLException("Out of range subscript encountered (length of insert exceeds array boundaries).");
3758
3759 // DataT& srcIn_dd = srcIn->dd;
3760 SizeT srcIx = 0; // this one simply runs from 0 to N_Elements(srcIn)
3761
3762 SizeT destEnd = destStart + len;
3763 for( SizeT destIx = destStart; destIx < destEnd; ++destIx)
3764 (*this)[ destIx] = (*srcIn)[ srcIx++];
3765 }
3766
3767 return;
3768 }
3769
3770 SizeT destStart; // 1-dim starting index
3771 // ATTENTION: dimension is used as an index here
3772 dimension ixDim = ixList->GetDimIx0( destStart);
3773 nDim--;
3774
3775 dimension srcDim=srcIn->Dim();
3776 srcDim.Purge(); // newly inserted: ticket #675
3777 SizeT len=srcDim[0]; // length of one segment to copy (one line of srcIn)
3778
3779 // SizeT nDim =RankIx(ixDim.Rank());
3780 SizeT srcNDim=RankIx(srcDim.Rank()); // number of source dimensions
3781 if( srcNDim < nDim) nDim=srcNDim;
3782
3783 // check limits (up to Rank to consider)
3784 for( SizeT dIx=0; dIx <= nDim; ++dIx)
3785 // check if in bounds of a
3786 if( (ixDim[dIx]+srcDim[dIx]) > this->dim[dIx])
3787 throw GDLException("Out of range subscript encountered (dimension of insert exceeds array boundaries for dimension " + i2s(dIx +1) + ").");
3788
3789 SizeT nCp=srcIn->Stride(nDim+1)/len; // number of OVERALL copy actions
3790
3791 // as lines are copied, we need the stride from 2nd dim on
3792 SizeT retStride[MAXRANK+1];
3793 for( SizeT a=0; a <= nDim; ++a) retStride[a]=srcDim.Stride(a+1)/len;
3794
3795 // a magic number, to reset destStart for this dimension
3796 SizeT resetStep[MAXRANK+1];
3797 for( SizeT a=1; a <= nDim; ++a)
3798 resetStep[a]=(retStride[a]-1)/retStride[a-1]*this->dim.Stride(a);
3799
3800 // SizeT destStart=this->dim.LongIndex(ixDim); // starting pos
3801
3802 // DataT& srcIn_dd = srcIn->dd;
3803
3804 SizeT srcIx=0; // this one simply runs from 0 to N_Elements(srcIn)
3805 for( SizeT c=1; c<=nCp; ++c) // linearized verison of nested loops
3806 {
3807 // copy one segment
3808 SizeT destEnd=destStart+len;
3809 for( SizeT destIx=destStart; destIx<destEnd; ++destIx)
3810 (*this)[destIx] = (*srcIn)[ srcIx++];
3811
3812 // update destStart for all dimensions
3813 if( c < nCp)
3814 for( SizeT a=1; a<=nDim; ++a)
3815 {
3816 if( c % retStride[a])
3817 {
3818 // advance to next
3819 destStart += this->dim.Stride(a);
3820 break;
3821 }
3822 else
3823 {
3824 // reset
3825 destStart -= resetStep[a];
3826 }
3827 }
3828 }
3829 }
3830
3831
3832 // used for concatenation, called from CatArray
3833 // assumes that everything is checked (see CatInfo)
3834 template<class Sp>
CatInsert(const Data_ * srcArr,const SizeT atDim,SizeT & at)3835 void Data_<Sp>::CatInsert (const Data_* srcArr, const SizeT atDim, SizeT& at)
3836 {
3837 // length of one segment to copy
3838 SizeT len = srcArr->dim.Stride (atDim + 1); // src array
3839
3840 SizeT nEl = srcArr->N_Elements ();
3841 // number of copy actions
3842 SizeT nCp = nEl / len;
3843
3844 // initial offset
3845 SizeT destStart = this->dim.Stride (atDim) * at; // dest array
3846 SizeT destEnd = destStart + len;
3847
3848 // number of elements to skip
3849 SizeT gap = this->dim.Stride (atDim + 1); // dest array
3850
3851 #ifdef _OPENMP
3852 //GD: speed up by using indexing that permit parallel and collapse.
3853 #pragma omp parallel for collapse(2) if (len*nCp >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= len*nCp))
3854 for (OMPInt c = 0; c < nCp; ++c)
3855 {
3856 // // copy one segment
3857 // SizeT destStartLoop = destStart + c * gap;
3858 // SizeT destEndLoop = destStartLoop + len;
3859 // SizeT srcIxLoop = c * len;
3860 // std::cerr<<"ncp="<<nCp<<",c="<<c<<": start: "<<destStart+c*gap<<" end: "<<destStart + c * gap + len<<" srcIxLoop="<<c * len<<std::endl;
3861 for (SizeT destIx = 0; destIx < len; destIx++) (*this)[destIx + destStart + c * gap] = (*srcArr)[ destIx + c * len];
3862 }
3863 #else // #ifdef _OPENMP
3864 SizeT srcIx = 0;
3865 for (SizeT c = 0; c < nCp; ++c)
3866 {
3867 // copy one segment
3868 for (SizeT destIx = destStart; destIx < destEnd; destIx++)
3869 (*this)[destIx] = (*srcArr)[ srcIx++];
3870
3871 // set new destination pointer
3872 destStart += gap;
3873 destEnd += gap;
3874 }
3875 #endif
3876
3877 SizeT add = srcArr->dim[atDim]; // update 'at'
3878 at += (add > 1) ? add : 1;
3879 }
3880
3881 // Logical True
3882 // integers
3883 template<class Sp>
LogTrue()3884 bool Data_<Sp>::LogTrue()
3885 {
3886 Ty s;
3887 if( !Scalar( s))
3888 throw GDLException("Expression must be a scalar or 1 element array in this context.",true,false);
3889 return (s != 0);
3890 }
3891 template<>
LogTrue()3892 bool Data_<SpDFloat>::LogTrue()
3893 {
3894 Ty s;
3895 if( !Scalar( s))
3896 throw GDLException("Expression must be a scalar or 1 element array in this context.",true,false);
3897 return (s != 0.0f);
3898 }
3899 template<>
LogTrue()3900 bool Data_<SpDDouble>::LogTrue()
3901 {
3902 Ty s;
3903 if( !Scalar( s))
3904 throw GDLException("Expression must be a scalar or 1 element array in this context.",true,false);
3905 return (s != 0.0);
3906 }
3907 template<>
LogTrue()3908 bool Data_<SpDString>::LogTrue()
3909 {
3910 Ty s;
3911 if( !Scalar( s))
3912 throw GDLException("Expression must be a scalar or 1 element array in this context.",true,false);
3913 return (s != "");
3914 }
3915 template<>
LogTrue()3916 bool Data_<SpDComplex>::LogTrue()
3917 {
3918 Ty s;
3919 if( !Scalar( s))
3920 throw GDLException("Expression must be a scalar or 1 element array in this context.",true,false);
3921 return (real(s) != 0.0 || imag(s) != 0.0);
3922 }
3923 template<>
LogTrue()3924 bool Data_<SpDComplexDbl>::LogTrue()
3925 {
3926 Ty s;
3927 if( !Scalar( s))
3928 throw GDLException("Expression must be a scalar or 1 element array in this context.",true,false);
3929 return (real(s) != 0.0 || imag(s) != 0.0);
3930 }
3931 template<>
LogTrue()3932 bool Data_<SpDPtr>::LogTrue()
3933 {
3934 Ty s;
3935 if( !Scalar( s))
3936 throw GDLException("Expression must be a scalar or 1 element array in this context.",true,false);
3937 return (s != 0);
3938 }
3939 template<>
LogTrue()3940 bool Data_<SpDObj>::LogTrue()
3941 {
3942 // ::_overloadIsTrue is handled in True()
3943
3944 return this->True();
3945 }
3946 // structs are not allowed
3947
3948 // indexed version
3949 // integers
3950 template<class Sp>
LogTrue(SizeT i)3951 bool Data_<Sp>::LogTrue(SizeT i)
3952 {
3953 return ((*this)[i] != 0);
3954 }
3955 template<>
LogTrue(SizeT i)3956 bool Data_<SpDFloat>::LogTrue(SizeT i)
3957 {
3958 return ((*this)[i] != 0.0f);
3959 }
3960 template<>
LogTrue(SizeT i)3961 bool Data_<SpDDouble>::LogTrue(SizeT i)
3962 {
3963 return ((*this)[i] != 0.0);
3964 }
3965 template<>
LogTrue(SizeT i)3966 bool Data_<SpDString>::LogTrue(SizeT i)
3967 {
3968 return ((*this)[i] != "");
3969 }
3970 template<>
LogTrue(SizeT i)3971 bool Data_<SpDComplex>::LogTrue(SizeT i)
3972 {
3973 return ((*this)[i].real() != 0.0 || (*this)[i].imag() != 0.0);
3974 }
3975 template<>
LogTrue(SizeT i)3976 bool Data_<SpDComplexDbl>::LogTrue(SizeT i)
3977 {
3978 return ((*this)[i].real() != 0.0 || (*this)[i].imag() != 0.0);
3979 }
3980 template<>
LogTrue(SizeT i)3981 bool Data_<SpDPtr>::LogTrue(SizeT i)
3982 {
3983 return (*this)[i] != 0;
3984 }
3985 template<>
LogTrue(SizeT i)3986 bool Data_<SpDObj>::LogTrue(SizeT i)
3987 {
3988 return (*this)[i] != 0;
3989 }
3990
3991 template<>
Rebin(const dimension & newDim,bool sample)3992 BaseGDL* Data_<SpDString>::Rebin( const dimension& newDim, bool sample)
3993 {
3994 throw GDLException("String expression not allowed in this context.");
3995 }
3996 template<>
Rebin(const dimension & newDim,bool sample)3997 BaseGDL* Data_<SpDComplex>::Rebin( const dimension& newDim, bool sample)
3998 {
3999 throw GDLException("Complex expression not allowed in this context.");
4000 }
4001 template<>
Rebin(const dimension & newDim,bool sample)4002 BaseGDL* Data_<SpDComplexDbl>::Rebin( const dimension& newDim, bool sample)
4003 {
4004 throw GDLException("Double complex expression not allowed in this context.");
4005 }
4006
4007
4008 // rebin over dimIx, new value: newDim
4009 // newDim != srcDim[ dimIx] -> compress or expand
4010 template< typename T>
Rebin1(T * src,const dimension & srcDim,SizeT dimIx,SizeT newDim,bool sample)4011 T* Rebin1( T* src,
4012 const dimension& srcDim,
4013 SizeT dimIx, SizeT newDim, bool sample)
4014 {
4015 SizeT nEl = src->N_Elements();
4016
4017 if( newDim == 0) newDim = 1;
4018
4019 // get dest dim and number of summations
4020 dimension destDim = srcDim;
4021
4022 destDim.MakeRank( dimIx + 1);
4023
4024 SizeT srcDimIx = destDim[ dimIx];
4025
4026 destDim.SetOneDim( dimIx, newDim);
4027
4028 SizeT resStride = destDim.Stride( dimIx);
4029
4030 // dimStride is also the number of linear src indexing
4031 SizeT dimStride = srcDim.Stride( dimIx);
4032 SizeT outerStride = srcDim.Stride( dimIx + 1);
4033
4034 SizeT rebinLimit = srcDimIx * dimStride;
4035
4036 if( newDim < srcDimIx) // compress
4037 {
4038
4039 SizeT ratio = srcDimIx / newDim;
4040
4041 if( sample)
4042 {
4043 T* res = new T( destDim, BaseGDL::NOZERO);
4044
4045 for( SizeT o=0; o < nEl; o += outerStride) // outer dim
4046 for( SizeT i=0; i < dimStride; ++i) // src element offset (lower dim)
4047 {
4048 SizeT oi = o+i;
4049 SizeT oiLimit = oi + rebinLimit;
4050
4051 for( SizeT s=oi; s<oiLimit; s += dimStride*ratio) // run over dim
4052 {
4053 (*res)[ (s / dimStride / ratio) * dimStride + i] = (*src)[ s];
4054 }
4055 }
4056
4057 return res;
4058 }
4059 else
4060 {
4061 T* res = new T( destDim); // zero fields
4062
4063 for( SizeT o=0; o < nEl; o += outerStride) // outer dim
4064 for( SizeT i=0; i < dimStride; ++i) // src element offset (lower dim)
4065 {
4066 SizeT oi = o+i;
4067 SizeT oiLimit = oi + rebinLimit;
4068
4069 for( SizeT s=oi; s<oiLimit; s += dimStride) // run over dim
4070 {
4071 // the way s indexes:
4072 // assume over b (compress index)
4073 // src[ a, b, c]
4074 // [ 0, 0, 0] [ 0, 1, 0] [ 0, 2, 0] ...
4075 // [ 1, 0, 0] [ 1, 1, 0] [ 1, 2, 0] ...
4076 // [ 2, 0, 0] [ 2, 1, 0] [ 2, 2, 0] ...
4077
4078 (*res)[ (s / dimStride / ratio) * dimStride + i] += (*src)[ s];
4079 }
4080 }
4081
4082 SizeT resEl = res->N_Elements();
4083 /*#pragma omp parallel if (resEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= resEl))
4084 {
4085 #pragma omp for*/
4086 for( int r=0; r < resEl; ++r)
4087 (*res)[ r] /= ratio;
4088 // }
4089 return res;
4090 }
4091 }
4092 else // expand
4093 {
4094 T* res = new T( destDim, BaseGDL::NOZERO);
4095
4096 if( sample)
4097 {
4098 SizeT ratio = newDim / srcDimIx;
4099
4100 for( SizeT o=0; o < nEl; o += outerStride) // outer dim
4101 for( SizeT i=0; i < dimStride; ++i) // src element offset (lower dim)
4102 {
4103 SizeT oi = o+i;
4104 SizeT oiLimit = oi + rebinLimit;
4105
4106 for( SizeT s=oi; s<oiLimit; s += dimStride) // run over dim
4107 {
4108 typename T::Ty src_s = (*src)[ s];
4109
4110 SizeT s_dimStride_ratio_dimStride_i =
4111 (s / dimStride) * ratio * dimStride + i;
4112
4113 for( SizeT r=0; r<ratio; ++r)
4114 {
4115 (*res)[ s_dimStride_ratio_dimStride_i + r * dimStride] =
4116 src_s;
4117 }
4118 }
4119 }
4120 }
4121 else
4122 {
4123 DLong64 ratio = newDim / srcDimIx; // make sure 32 bit integers are working also
4124
4125 for( SizeT o=0; o < nEl; o += outerStride) // outer dim
4126 for( SizeT i=0; i < dimStride; ++i) // src element offset (lower dim)
4127 {
4128 SizeT oi = o+i;
4129 SizeT oiLimit = oi + rebinLimit;
4130
4131 for( SizeT s=oi; s<oiLimit; s += dimStride) // run over dim
4132 {
4133 typename T::Ty first = (*src)[ s];
4134 typename T::Ty next =
4135 (*src)[ (s+dimStride)<oiLimit?s+dimStride:s];
4136
4137 SizeT s_dimStride_ratio_dimStride_i =
4138 (s / dimStride) * ratio * dimStride + i;
4139 for( DLong64 r=0; r<ratio; ++r)
4140 {
4141 (*res)[ s_dimStride_ratio_dimStride_i + r * dimStride] =
4142 (first * (ratio - r) + next * r) / ratio; // 64 bit temporary
4143 }
4144 }
4145 }
4146 }
4147
4148 return res;
4149 }
4150 }
4151
4152 // for integer
4153 template< typename T, typename TNext>
4154 T* Rebin1Int( T* src,
4155 const dimension& srcDim,
4156 SizeT dimIx, SizeT newDim, bool sample)
4157 {
4158 SizeT nEl = src->N_Elements();
4159
4160 if( newDim == 0) newDim = 1;
4161
4162 // get dest dim and number of summations
4163 dimension destDim = srcDim;
4164
4165 destDim.MakeRank( dimIx + 1);
4166
4167 SizeT srcDimIx = destDim[ dimIx];
4168
4169 destDim.SetOneDim( dimIx, newDim);
4170
4171 SizeT resStride = destDim.Stride( dimIx);
4172
4173 // dimStride is also the number of linear src indexing
4174 SizeT dimStride = srcDim.Stride( dimIx);
4175 SizeT outerStride = srcDim.Stride( dimIx + 1);
4176
4177 SizeT rebinLimit = srcDimIx * dimStride;
4178
4179 if( newDim < srcDimIx) // compress
4180 {
4181
4182 SizeT ratio = srcDimIx / newDim;
4183
4184 if( sample)
4185 {
4186 T* res = new T( destDim, BaseGDL::NOZERO);
4187
4188 SizeT ratio = srcDimIx / newDim;
4189
4190 for( SizeT o=0; o < nEl; o += outerStride) // outer dim
4191 for( SizeT i=0; i < dimStride; ++i) // src element offset (lower dim)
4192 {
4193 SizeT oi = o+i;
4194 SizeT oiLimit = oi + rebinLimit;
4195
4196 for( SizeT s=oi; s<oiLimit; s += dimStride*ratio) // run over dim
4197 {
4198 (*res)[ (s / dimStride / ratio) * dimStride + i] = (*src)[ s];
4199 }
4200 }
4201
4202 return res;
4203 }
4204 else
4205 {
4206 T* res = new T( destDim); // zero fields
4207
4208 for( SizeT o=0; o < nEl; o += outerStride) // outer dim
4209 for( SizeT i=0; i < dimStride; ++i) // src element offset (lower dim)
4210 {
4211 SizeT oi = o+i;
4212 SizeT oiLimit = oi + rebinLimit;
4213
4214 TNext tmp = 0;
4215 for( SizeT s=oi; s<oiLimit; s += dimStride) // run over dim
4216 {
4217 tmp += (*src)[ s];
4218
4219 if( (s / dimStride) % ratio == (ratio - 1))
4220 {
4221 (*res)[ (s / dimStride / ratio) * dimStride + i] = tmp / ratio;
4222 tmp = 0;
4223 }
4224 }
4225 }
4226
4227 return res;
4228 }
4229 }
4230 else // expand
4231 {
4232 T* res = new T( destDim, BaseGDL::NOZERO);
4233
4234 if( sample)
4235 {
4236 SizeT ratio = newDim / srcDimIx;
4237
4238 for( SizeT o=0; o < nEl; o += outerStride) // outer dim
4239 for( SizeT i=0; i < dimStride; ++i) // src element offset (lower dim)
4240 {
4241 SizeT oi = o+i;
4242 SizeT oiLimit = oi + rebinLimit;
4243
4244 for( SizeT s=oi; s<oiLimit; s += dimStride) // run over dim
4245 {
4246 typename T::Ty src_s = (*src)[ s];
4247
4248 SizeT s_dimStride_ratio_dimStride_i =
4249 (s / dimStride) * ratio * dimStride + i;
4250
4251 for( SizeT r=0; r<ratio; ++r)
4252 {
4253 (*res)[ s_dimStride_ratio_dimStride_i + r * dimStride] =
4254 src_s;
4255 }
4256 }
4257 }
4258 }
4259 else
4260 {
4261 DLong64 ratio = newDim / srcDimIx; // make sure 32 bit integers are working also
4262
4263 for( SizeT o=0; o < nEl; o += outerStride) // outer dim
4264 for( SizeT i=0; i < dimStride; ++i) // src element offset (lower dim)
4265 {
4266 SizeT oi = o+i;
4267 SizeT oiLimit = oi + rebinLimit;
4268
4269 for( SizeT s=oi; s<oiLimit; s += dimStride) // run over dim
4270 {
4271 typename T::Ty first = (*src)[ s];
4272 typename T::Ty next =
4273 (*src)[ (s+dimStride)<oiLimit?s+dimStride:s];
4274
4275 SizeT s_dimStride_ratio_dimStride_i =
4276 (s / dimStride) * ratio * dimStride + i;
4277 for( DLong64 r=0; r<ratio; ++r)
4278 {
4279 (*res)[ s_dimStride_ratio_dimStride_i + r * dimStride] =
4280 (first * (ratio - r) + next * r) / ratio; // 64 bit temporary
4281 }
4282 }
4283 }
4284 }
4285
4286 return res;
4287 }
4288 }
4289
4290 // for float, double
4291 template<class Sp>
Rebin(const dimension & newDim,bool sample)4292 BaseGDL* Data_<Sp>::Rebin( const dimension& newDim, bool sample)
4293 {
4294 SizeT resRank = newDim.Rank();
4295 SizeT srcRank = this->Rank();
4296
4297 SizeT nDim;
4298 if( resRank < srcRank)
4299 nDim = srcRank;
4300 else
4301 nDim = resRank;
4302
4303 dimension actDim( this->dim);
4304 Data_* actIn = this;
4305
4306 // 1st compress
4307 for( SizeT d=0; d<nDim; ++d)
4308 if( newDim[d] < this->dim[d])
4309 { // compress
4310
4311 Data_* act = Rebin1( actIn, actDim, d, newDim[d], sample);
4312 actDim = act->Dim();
4313
4314 if( actIn != this) GDLDelete(actIn);
4315 actIn = act;
4316 }
4317
4318 // 2nd expand
4319 for( SizeT d=0; d<nDim; ++d)
4320 if( newDim[ d] > this->dim[d])
4321 { // expand
4322
4323 Data_* act = Rebin1( actIn, actDim, d, newDim[d], sample);
4324 actDim = act->Dim();
4325
4326 if( actIn != this) GDLDelete(actIn);
4327 actIn = act;
4328 }
4329
4330 if( actIn == this) return actIn->Dup();
4331 return actIn;
4332 }
4333
4334 // integer types
4335 template<>
Rebin(const dimension & newDim,bool sample)4336 BaseGDL* Data_<SpDByte>::Rebin( const dimension& newDim, bool sample)
4337 {
4338 SizeT resRank = newDim.Rank();
4339 SizeT srcRank = Rank();
4340
4341 SizeT nDim;
4342 if( resRank < srcRank)
4343 nDim = srcRank;
4344 else
4345 nDim = resRank;
4346
4347 dimension actDim( dim);
4348 Data_* actIn = this;
4349
4350 // 1st compress
4351 for( SizeT d=0; d<nDim; ++d)
4352 if( newDim[d] < dim[d])
4353 { // compress
4354
4355 Data_* act = Rebin1Int<DByteGDL, DULong64>( actIn, actDim, d, newDim[d], sample);
4356 actDim = act->Dim();
4357
4358 if( actIn != this) GDLDelete(actIn);
4359 actIn = act;
4360 }
4361
4362 // 2nd expand
4363 for( SizeT d=0; d<nDim; ++d)
4364 if( newDim[ d] > dim[d])
4365 { // expand
4366
4367 Data_* act = Rebin1Int<DByteGDL, DULong64>( actIn, actDim, d, newDim[d], sample);
4368 actDim = act->Dim();
4369
4370 if( actIn != this) GDLDelete(actIn);
4371 actIn = act;
4372 }
4373
4374 if( actIn == this) return actIn->Dup();
4375 return actIn;
4376 }
4377 template<>
Rebin(const dimension & newDim,bool sample)4378 BaseGDL* Data_<SpDInt>::Rebin( const dimension& newDim, bool sample)
4379 {
4380 SizeT resRank = newDim.Rank();
4381 SizeT srcRank = Rank();
4382
4383 SizeT nDim;
4384 if( resRank < srcRank)
4385 nDim = srcRank;
4386 else
4387 nDim = resRank;
4388
4389 dimension actDim( dim);
4390 Data_* actIn = this;
4391
4392 // 1st compress
4393 for( SizeT d=0; d<nDim; ++d)
4394 if( newDim[d] < dim[d])
4395 { // compress
4396
4397 Data_* act = Rebin1Int<DIntGDL, DLong64>( actIn, actDim, d, newDim[d], sample);
4398 actDim = act->Dim();
4399
4400 if( actIn != this) GDLDelete(actIn);
4401 actIn = act;
4402 }
4403
4404 // 2nd expand
4405 for( SizeT d=0; d<nDim; ++d)
4406 if( newDim[ d] > dim[d])
4407 { // expand
4408
4409 Data_* act = Rebin1Int<DIntGDL, DLong64>( actIn, actDim, d, newDim[d], sample);
4410 actDim = act->Dim();
4411
4412 if( actIn != this) GDLDelete(actIn);
4413 actIn = act;
4414 }
4415
4416 if( actIn == this) return actIn->Dup();
4417 return actIn;
4418 }
4419 template<>
Rebin(const dimension & newDim,bool sample)4420 BaseGDL* Data_<SpDUInt>::Rebin( const dimension& newDim, bool sample)
4421 {
4422 SizeT resRank = newDim.Rank();
4423 SizeT srcRank = Rank();
4424
4425 SizeT nDim;
4426 if( resRank < srcRank)
4427 nDim = srcRank;
4428 else
4429 nDim = resRank;
4430
4431 dimension actDim( dim);
4432 Data_* actIn = this;
4433
4434 // 1st compress
4435 for( SizeT d=0; d<nDim; ++d)
4436 if( newDim[d] < dim[d])
4437 { // compress
4438
4439 Data_* act = Rebin1Int<DUIntGDL, DULong64>( actIn, actDim, d, newDim[d], sample);
4440 actDim = act->Dim();
4441
4442 if( actIn != this) GDLDelete(actIn);
4443 actIn = act;
4444 }
4445
4446 // 2nd expand
4447 for( SizeT d=0; d<nDim; ++d)
4448 if( newDim[ d] > dim[d])
4449 { // expand
4450
4451 Data_* act = Rebin1Int<DUIntGDL, DULong64>( actIn, actDim, d, newDim[d], sample);
4452 actDim = act->Dim();
4453
4454 if( actIn != this) GDLDelete(actIn);
4455 actIn = act;
4456 }
4457
4458 if( actIn == this) return actIn->Dup();
4459 return actIn;
4460 }
4461 template<>
Rebin(const dimension & newDim,bool sample)4462 BaseGDL* Data_<SpDLong>::Rebin( const dimension& newDim, bool sample)
4463 {
4464 SizeT resRank = newDim.Rank();
4465 SizeT srcRank = Rank();
4466
4467 SizeT nDim;
4468 if( resRank < srcRank)
4469 nDim = srcRank;
4470 else
4471 nDim = resRank;
4472
4473 dimension actDim( dim);
4474 Data_* actIn = this;
4475
4476 // 1st compress
4477 for( SizeT d=0; d<nDim; ++d)
4478 if( newDim[d] < dim[d])
4479 { // compress
4480
4481 Data_* act = Rebin1Int<DLongGDL, DLong64>( actIn, actDim, d, newDim[d], sample);
4482 actDim = act->Dim();
4483
4484 if( actIn != this) GDLDelete(actIn);
4485 actIn = act;
4486 }
4487
4488 // 2nd expand
4489 for( SizeT d=0; d<nDim; ++d)
4490 if( newDim[ d] > dim[d])
4491 { // expand
4492
4493 Data_* act = Rebin1Int<DLongGDL, DLong64>( actIn, actDim, d, newDim[d], sample);
4494 actDim = act->Dim();
4495
4496 if( actIn != this) GDLDelete(actIn);
4497 actIn = act;
4498 }
4499
4500 if( actIn == this) return actIn->Dup();
4501 return actIn;
4502 }
4503 template<>
Rebin(const dimension & newDim,bool sample)4504 BaseGDL* Data_<SpDULong>::Rebin( const dimension& newDim, bool sample)
4505 {
4506 SizeT resRank = newDim.Rank();
4507 SizeT srcRank = Rank();
4508
4509 SizeT nDim;
4510 if( resRank < srcRank)
4511 nDim = srcRank;
4512 else
4513 nDim = resRank;
4514
4515 dimension actDim( dim);
4516 Data_* actIn = this;
4517
4518 // 1st compress
4519 for( SizeT d=0; d<nDim; ++d)
4520 if( newDim[d] < dim[d])
4521 { // compress
4522
4523 Data_* act = Rebin1Int<DULongGDL, DULong64>( actIn, actDim, d, newDim[d], sample);
4524 actDim = act->Dim();
4525
4526 if( actIn != this) GDLDelete(actIn);
4527 actIn = act;
4528 }
4529
4530 // 2nd expand
4531 for( SizeT d=0; d<nDim; ++d)
4532 if( newDim[ d] > dim[d])
4533 { // expand
4534
4535 Data_* act = Rebin1Int<DULongGDL, DULong64>( actIn, actDim, d, newDim[d], sample);
4536 actDim = act->Dim();
4537
4538 if( actIn != this) GDLDelete(actIn);
4539 actIn = act;
4540 }
4541
4542 if( actIn == this) return actIn->Dup();
4543 return actIn;
4544 }
4545
4546 // plain copy of nEl from src
4547 // no checking
4548 template<class Sp>
Assign(BaseGDL * src,SizeT nEl)4549 void Data_<Sp>::Assign( BaseGDL* src, SizeT nEl)
4550 {
4551 Data_* srcT; // = dynamic_cast<Data_*>( src);
4552
4553 Guard< Data_> srcTGuard;
4554 if( src->Type() != Data_::t)
4555 {
4556 srcT = static_cast<Data_*>( src->Convert2( Data_::t, BaseGDL::COPY));
4557 srcTGuard.Init( srcT);
4558 }
4559 else
4560 {
4561 srcT = static_cast<Data_*>( src);
4562 }
4563 /*#pragma omp parallel if (nEl >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nEl))
4564 {
4565 #pragma omp for*/
4566 for(long k=0; k < nEl; ++k)
4567 {
4568 (*this)[ k] = (*srcT)[ k];
4569 }// }
4570 }
4571
4572
4573
4574
4575
4576 // return a new type of itself (only for one dimensional case)
4577 template<class Sp>
NewIx(SizeT ix)4578 BaseGDL* Data_<Sp>::NewIx( SizeT ix)
4579 {
4580 return new Data_( (*this)[ ix]);
4581 }
4582 template<class Sp>
NewIx(AllIxBaseT * ix,const dimension * dIn)4583 Data_<Sp>* Data_<Sp>::NewIx( AllIxBaseT* ix, const dimension* dIn)
4584 {
4585 SizeT nCp = ix->size();
4586 Data_* res=Data_::New( *dIn, BaseGDL::NOZERO);
4587 /*#pragma omp parallel if (nCp >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nCp))
4588 {
4589 #pragma omp for*/
4590 for( int c=0; c<nCp; ++c)
4591 (*res)[c]=(*this)[ (*ix)[ c]];
4592 // }
4593 return res;
4594 }
4595 template<class Sp>
NewIxFrom(SizeT s)4596 Data_<Sp>* Data_<Sp>::NewIxFrom( SizeT s)
4597 {
4598 SizeT nCp = dd.size() - s;
4599 Data_* res=Data_::New( dimension( nCp), BaseGDL::NOZERO);
4600 /*#pragma omp parallel if (nCp >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nCp))
4601 {
4602 #pragma omp for*/
4603 for( int c=0; c<nCp; ++c)
4604 (*res)[c]=(*this)[ s+c];
4605 // }
4606 return res;
4607 }
4608 template<class Sp>
NewIxFrom(SizeT s,SizeT e)4609 Data_<Sp>* Data_<Sp>::NewIxFrom( SizeT s, SizeT e)
4610 {
4611 SizeT nCp = e - s + 1;
4612 Data_* res=Data_::New( dimension( nCp), BaseGDL::NOZERO);
4613 /*#pragma omp parallel if (nCp >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nCp))
4614 {
4615 #pragma omp for*/
4616 for( int c=0; c<nCp; ++c)
4617 (*res)[c]=(*this)[ s+c];
4618 // }
4619 return res;
4620 }
4621 template<class Sp>
NewIxFromStride(SizeT s,SizeT stride)4622 Data_<Sp>* Data_<Sp>::NewIxFromStride( SizeT s, SizeT stride)
4623 {
4624 SizeT nCp = (dd.size() - s + stride - 1)/stride;
4625 Data_* res=Data_::New( dimension( nCp), BaseGDL::NOZERO);
4626 for( SizeT c=0; c<nCp; ++c, s += stride)
4627 (*res)[c]=(*this)[ s];
4628 return res;
4629 }
4630 template<class Sp>
NewIxFromStride(SizeT s,SizeT e,SizeT stride)4631 Data_<Sp>* Data_<Sp>::NewIxFromStride( SizeT s, SizeT e, SizeT stride)
4632 {
4633 SizeT nCp = (e - s + stride)/stride;
4634 Data_* res=Data_::New( dimension( nCp), BaseGDL::NOZERO);
4635 for( SizeT c=0; c<nCp; ++c, s += stride)
4636 (*res)[c]=(*this)[ s];
4637 return res;
4638 }
4639
4640
4641
4642 template<class Sp>
NewIx(BaseGDL * ix,bool strict)4643 Data_<Sp>* Data_<Sp>::NewIx( BaseGDL* ix, bool strict)
4644 {
4645 assert( ix->Type() != GDL_UNDEF);
4646
4647 // no type checking needed here: GetAsIndex() will fail with grace
4648 // int typeCheck = DTypeOrder[ dType];
4649 // if( typeCheck >= 100)
4650 // throw GDLException("Type "+ix->TypeStr()+" not allowed as subscript.");
4651
4652 SizeT nElem = ix->N_Elements();
4653
4654 Data_* res = New( ix->Dim(), BaseGDL::NOZERO);
4655 Guard<Data_> guard( res);
4656
4657 SizeT upper = dd.size() - 1;
4658 Ty upperVal = (*this)[ upper];
4659 if( strict)
4660 {
4661 for(SizeT i = 0 ; i < nElem; ++i)
4662 {
4663 SizeT actIx = ix->GetAsIndexStrict( i);
4664 if( actIx > upper)
4665 throw GDLException("Array used to subscript array "
4666 "contains out of range (>) subscript (at index: "+i2s(i)+").");
4667 (*res)[i]= (*this)[ actIx];
4668 }
4669 }
4670 else // not strict
4671 {
4672 for(SizeT i = 0 ; i < nElem; ++i)
4673 {
4674 SizeT actIx = ix->GetAsIndex( i);
4675 if( actIx >= upper)
4676 (*res)[i] = upperVal;
4677 else
4678 (*res)[i]= (*this)[ actIx];
4679 }
4680 }
4681 return guard.release();
4682 }
4683
4684
4685 // unsigned types
GetAsIndex(SizeT i) const4686 template<class Sp> SizeT Data_<Sp>::GetAsIndex( SizeT i) const
4687 {
4688 return (*this)[ i];
4689 }
GetAsIndexStrict(SizeT i) const4690 template<class Sp> SizeT Data_<Sp>::GetAsIndexStrict( SizeT i) const
4691 {
4692 return (*this)[ i]; // good for unsigned types
4693 }
4694
4695 template<>
GetAsIndex(SizeT i) const4696 SizeT Data_<SpDInt>::GetAsIndex( SizeT i) const
4697 {
4698 if( (*this)[i] < 0)
4699 return 0;
4700 return (*this)[i];
4701 }
4702 template<>
GetAsIndexStrict(SizeT i) const4703 SizeT Data_<SpDInt>::GetAsIndexStrict( SizeT i) const
4704 {
4705 if( (*this)[i] < 0)
4706 throw GDLException(-1,NULL,"Array used to subscript array "
4707 "contains out of range (<0) subscript (at index: " + i2s(i) + ").",true,false);
4708 return (*this)[i];
4709 }
4710 template<>
GetAsIndex(SizeT i) const4711 SizeT Data_<SpDLong>::GetAsIndex( SizeT i) const
4712 {
4713 if( (*this)[i] < 0)
4714 return 0;
4715 return (*this)[i];
4716 }
4717 template<>
GetAsIndexStrict(SizeT i) const4718 SizeT Data_<SpDLong>::GetAsIndexStrict( SizeT i) const
4719 {
4720 if( (*this)[i] < 0)
4721 throw GDLException(-1,NULL,"Array used to subscript array "
4722 "contains out of range (<0) subscript (at index: " + i2s(i) + ").",true,false);
4723 return (*this)[i];
4724 }
4725 template<>
GetAsIndex(SizeT i) const4726 SizeT Data_<SpDLong64>::GetAsIndex( SizeT i) const
4727 {
4728 if( (*this)[i] < 0)
4729 return 0;
4730 return (*this)[i];
4731 }
4732 template<>
GetAsIndexStrict(SizeT i) const4733 SizeT Data_<SpDLong64>::GetAsIndexStrict( SizeT i) const
4734 {
4735 if( (*this)[i] < 0)
4736 throw GDLException(-1,NULL,"Array used to subscript array "
4737 "contains out of range (<0) subscript (at index: " + i2s(i) + ").",true,false);
4738 return (*this)[i];
4739 }
4740 template<>
GetAsIndex(SizeT i) const4741 SizeT Data_<SpDFloat>::GetAsIndex( SizeT i) const
4742 {
4743 if( (*this)[i] <= 0.0)
4744 return 0;
4745 return Real2Int<SizeT,float>((*this)[i]);
4746 }
4747 template<>
GetAsIndexStrict(SizeT i) const4748 SizeT Data_<SpDFloat>::GetAsIndexStrict( SizeT i) const
4749 {
4750 if( (*this)[i] <= -1.0)
4751 throw GDLException(-1,NULL,"Array used to subscript array "
4752 "contains out of range (<0) subscript (at index: " + i2s(i) + ").",true,false);
4753 if( (*this)[i] <= 0.0)
4754 return 0;
4755 return Real2Int<SizeT,float>((*this)[i]);
4756 }
4757 template<>
GetAsIndex(SizeT i) const4758 SizeT Data_<SpDDouble>::GetAsIndex( SizeT i) const
4759 {
4760 if( (*this)[i] <= 0.0)
4761 return 0;
4762 return Real2Int<SizeT,double>((*this)[i]);
4763 }
4764 template<>
GetAsIndexStrict(SizeT i) const4765 SizeT Data_<SpDDouble>::GetAsIndexStrict( SizeT i) const
4766 {
4767 if( (*this)[i] <= -1.0)
4768 throw GDLException(-1,NULL,"Array used to subscript array "
4769 "contains out of range (<0) subscript (at index: " + i2s(i) + ").",true,false);
4770 if( (*this)[i] <= 0.0)
4771 return 0;
4772 return Real2Int<SizeT,double>((*this)[i]);
4773 }
4774 template<>
GetAsIndex(SizeT i) const4775 SizeT Data_<SpDString>::GetAsIndex( SizeT i) const
4776 {
4777 const char* cStart=(*this)[i].c_str();
4778 char* cEnd;
4779 long l=strtol(cStart,&cEnd,10);
4780 if( cEnd == cStart)
4781 {
4782 Warning("Type conversion error: Unable to convert given STRING to LONG (at index: " + i2s(i) + ")");
4783 return 0;
4784 }
4785 if( l < 0)
4786 return 0;
4787 return l;
4788 }
4789 template<>
GetAsIndexStrict(SizeT i) const4790 SizeT Data_<SpDString>::GetAsIndexStrict( SizeT i) const
4791 {
4792 const char* cStart=(*this)[i].c_str();
4793 char* cEnd;
4794 long l=strtol(cStart,&cEnd,10);
4795 if( cEnd == cStart)
4796 {
4797 Warning("Type conversion error: Unable to convert given STRING to LONG (at index: " + i2s(i) + ")");
4798 return 0;
4799 }
4800 if( l < 0)
4801 throw GDLException(-1,NULL,"Array used to subscript array "
4802 "contains out of range (<0) subscript.",true,false);
4803 return l;
4804 }
4805
4806 template<>
GetAsIndex(SizeT i) const4807 SizeT Data_<SpDComplex>::GetAsIndex( SizeT i) const
4808 {
4809 if( real((*this)[i]) <= 0.0)
4810 return 0;
4811 return Real2Int<SizeT,float>(real((*this)[i]));
4812 }
4813 template<>
GetAsIndexStrict(SizeT i) const4814 SizeT Data_<SpDComplex>::GetAsIndexStrict( SizeT i) const
4815 {
4816 if( real((*this)[i]) <= -1.0)
4817 throw GDLException(-1,NULL,"Array used to subscript array "
4818 "contains out of range (<0) subscript (at index: " + i2s(i) + ").",true,false);
4819 if( real((*this)[i]) <= 0.0)
4820 return 0;
4821 return Real2Int<SizeT,float>(real((*this)[i]));
4822 }
4823 template<>
GetAsIndex(SizeT i) const4824 SizeT Data_<SpDComplexDbl>::GetAsIndex( SizeT i) const
4825 {
4826 if( real((*this)[i]) <= 0.0)
4827 return 0;
4828 return Real2Int<SizeT,double>(real((*this)[i]));
4829 }
4830 template<>
GetAsIndexStrict(SizeT i) const4831 SizeT Data_<SpDComplexDbl>::GetAsIndexStrict( SizeT i) const
4832 {
4833 if( real((*this)[i]) <= -1.0)
4834 throw GDLException(-1,NULL,"Array used to subscript array "
4835 "contains out of range (<0) subscript (at index: " + i2s(i) + ").",true,false);
4836 if( real((*this)[i]) <= 0.0)
4837 return 0;
4838 return Real2Int<SizeT,double>(real((*this)[i]));
4839 }
4840
4841 #include "instantiate_templates.hpp"
4842