1 #ifndef WIN32
2 #include <unistd.h>
3 #include <sys/mman.h>
4 #include <fcntl.h>
5 #define USE_MMAP
6 #endif
7
8 #include "EXTERN.h" /* std perl include */
9 #include "perl.h" /* std perl include */
10 #include "XSUB.h" /* XSUB include */
11
12 #if defined(CONTEXT)
13 #undef CONTEXT
14 #endif
15
16 #define PDL_CORE /* For certain ifdefs */
17 #include "pdl.h" /* Data structure declarations */
18 #include "pdlcore.h" /* Core declarations */
19
20 #if BADVAL
21 # if !BADVAL_USENAN
22 #include <float.h>
23 # endif
24 #include <limits.h>
25 #endif
26
27 /* Return a integer or numeric scalar as approroate */
28
29 #define setflag(reg,flagval,val) (val?(reg |= flagval):(reg &= ~flagval))
30
31 Core PDL; /* Struct holding pointers to shared C routines */
32
33 #ifdef FOO
pdl__Core_get_Core()34 Core *pdl__Core_get_Core() /* INTERNAL TO CORE! DONT CALL FROM OUTSIDE */
35 {
36 return PDL;
37 }
38 #endif
39
40 int pdl_debugging=0;
41 int pdl_autopthread_targ = 0; /* No auto-pthreading unless set using the set_autopthread_targ */
42 int pdl_autopthread_actual = 0;
43 int pdl_autopthread_size = 1;
44
45 #define CHECKP(p) if ((p) == NULL) croak("Out of memory")
46
pdl_packint(SV * sv,int * ndims)47 static PDL_Indx* pdl_packint( SV* sv, int *ndims ) {
48
49 SV* bar;
50 AV* array;
51 int i;
52 PDL_Indx *dims;
53
54 if (!(SvROK(sv) && SvTYPE(SvRV(sv))==SVt_PVAV)) /* Test */
55 return NULL;
56 array = (AV *) SvRV(sv); /* dereference */
57 *ndims = (int) av_len(array) + 1; /* Number of dimensions */
58 /* Array space */
59 dims = (PDL_Indx *) pdl_malloc( (*ndims) * sizeof(*dims) );
60 CHECKP(dims);
61
62 for(i=0; i<(*ndims); i++) {
63 bar = *(av_fetch( array, i, 0 )); /* Fetch */
64 dims[i] = (PDL_Indx) SvIV(bar);
65 }
66 return dims;
67 }
68
pdl_unpackint(PDL_Indx * dims,int ndims)69 static SV* pdl_unpackint ( PDL_Indx *dims, int ndims ) {
70
71 AV* array;
72 int i;
73
74 array = newAV();
75
76 for(i=0; i<ndims; i++) /* if ndims == 0, nothing stored -> ok */
77 av_store( array, i, newSViv( (IV)dims[i] ) );
78
79 return (SV*) array;
80 }
81
82 /*
83 * Free the data if possible; used by mmapper
84 * Moved from pdlhash.c July 10 2006 DJB
85 */
pdl_freedata(pdl * a)86 static void pdl_freedata (pdl *a) {
87 if(a->datasv) {
88 SvREFCNT_dec(a->datasv);
89 a->datasv=0;
90 a->data=0;
91 } else if(a->data) {
92 die("Trying to free data of untouchable (mmapped?) pdl");
93 }
94 }
95
96 #if BADVAL
97
98 #ifdef FOOFOO_PROPAGATE_BADFLAG
99
100 /*
101 * this seems to cause an infinite loop in between tests 42 & 43 of
102 * t/bad.t - ie
103 *
104 * $a = sequence( byte, 2, 3 );
105 * $b = $a->slice("(1),:");
106 * my $mask = sequence( byte, 2, 3 );
107 * $mask = $mask->setbadif( ($mask % 3) == 2 );
108 * print "a,b == ", $a->badflag, ",", $b->badflag, "\n";
109 * $a->inplace->copybad( $mask ); <-- think this is the call
110 * print "a,b == ", $a->badflag, ",", $b->badflag, "\n";
111 * print "$a $b\n";
112 * ok( $b->badflag, 1 );
113 *
114 */
115
116 /* used by propagate_badflag() */
117
propagate_badflag_children(pdl * it,int newval)118 void propagate_badflag_children( pdl *it, int newval ) {
119 PDL_DECL_CHILDLOOP(it)
120 PDL_START_CHILDLOOP(it)
121 {
122 pdl_trans *trans = PDL_CHILDLOOP_THISCHILD(it);
123 int i;
124
125 for( i = trans->vtable->nparents;
126 i < trans->vtable->npdls;
127 i++ ) {
128
129 pdl *child = trans->pdls[i];
130
131 if ( newval ) child->state |= PDL_BADVAL;
132 else child->state &= ~PDL_BADVAL;
133
134 /* make sure we propagate to grandchildren, etc */
135 propagate_badflag_children( child, newval );
136
137 } /* for: i */
138 }
139 PDL_END_CHILDLOOP(it)
140 } /* propagate_badflag_children */
141
142 /* used by propagate_badflag() */
143
propagate_badflag_parents(pdl * it)144 void propagate_badflag_parents( pdl *it ) {
145 PDL_DECL_CHILDLOOP(it)
146 PDL_START_CHILDLOOP(it)
147 {
148 pdl_trans *trans = PDL_CHILDLOOP_THISCHILD(it);
149 int i;
150
151 for( i = 0;
152 i < trans->vtable->nparents;
153 i++ ) {
154
155 pdl *parent = trans->pdls[i];
156
157 /* only sets allowed here */
158 parent->state |= PDL_BADVAL;
159
160 /* make sure we propagate to grandparents, etc */
161 propagate_badflag_parents( parent );
162
163 } /* for: i */
164 }
165 PDL_END_CHILDLOOP(it)
166 } /* propagate_badflag_parents */
167
168 /*
169 * we want to change the bad flag of the children
170 * (newval = 1 means set flag, 0 means clear it).
171 * If newval == 1, then we also loop through the
172 * parents, setting their bad flag
173 *
174 * thanks to Christian Soeller for this
175 */
176
propagate_badflag(pdl * it,int newval)177 void propagate_badflag( pdl *it, int newval ) {
178 /* only do anything if the flag has changed - do we need this check ? */
179 if ( newval ) {
180 if ( (it->state & PDL_BADVAL) == 0 ) {
181 propagate_badflag_parents( it );
182 propagate_badflag_children( it, newval );
183 }
184 } else {
185 if ( (it->state & PDL_BADVAL) > 0 ) {
186 propagate_badflag_children( it, newval );
187 }
188
189 }
190
191 } /* propagate_badflag */
192
193 #else /* FOOFOO_PROPAGATE_BADFLAG */
194
195 /* newval = 1 means set flag, 0 means clear it */
196 /* thanks to Christian Soeller for this */
197
propagate_badflag(pdl * it,int newval)198 void propagate_badflag( pdl *it, int newval ) {
199 PDL_DECL_CHILDLOOP(it)
200 PDL_START_CHILDLOOP(it)
201 {
202 pdl_trans *trans = PDL_CHILDLOOP_THISCHILD(it);
203 int i;
204
205 for( i = trans->vtable->nparents;
206 i < trans->vtable->npdls; i++ ) {
207
208 pdl *child = trans->pdls[i];
209
210 if ( newval ) child->state |= PDL_BADVAL;
211 else child->state &= ~PDL_BADVAL;
212
213 /* make sure we propagate to grandchildren, etc */
214 propagate_badflag( child, newval );
215
216 } /* for: i */
217 }
218 PDL_END_CHILDLOOP(it)
219 } /* propagate_badflag */
220
221 #endif /* FOOFOO_PROPAGATE_BADFLAG */
222
propagate_badvalue(pdl * it)223 void propagate_badvalue( pdl *it ) {
224 PDL_DECL_CHILDLOOP(it)
225 PDL_START_CHILDLOOP(it)
226 {
227 pdl_trans *trans = PDL_CHILDLOOP_THISCHILD(it);
228 int i;
229
230 for( i = trans->vtable->nparents;
231 i < trans->vtable->npdls; i++ ) {
232
233 pdl *child = trans->pdls[i];
234
235 child->has_badvalue = 1;
236 child->badvalue = it->badvalue;
237
238 /* make sure we propagate to grandchildren, etc */
239 propagate_badvalue( child );
240
241 } /* for: i */
242 }
243 PDL_END_CHILDLOOP(it)
244 } /* propagate_badvalue */
245
246
247 /* this is horrible - the routines from bad should perhaps be here instead ? */
pdl_get_badvalue(int datatype)248 PDL_Anyval pdl_get_badvalue( int datatype ) {
249 PDL_Anyval retval = { -1, 0 };
250 switch ( datatype ) {
251
252 #include "pdldataswitch.c"
253
254 default:
255 croak("Unknown type sent to pdl_get_badvalue\n");
256 }
257 return retval;
258 } /* pdl_get_badvalue() */
259
260
pdl_get_pdl_badvalue(pdl * it)261 PDL_Anyval pdl_get_pdl_badvalue( pdl *it ) {
262 PDL_Anyval retval = { -1, 0 };
263 int datatype;
264
265 #if BADVAL_PER_PDL
266 if (it->has_badvalue) {
267 retval = it->badvalue;
268 } else {
269 datatype = it->datatype;
270 retval = pdl_get_badvalue( datatype );
271 }
272 #else
273 datatype = it->datatype;
274 retval = pdl_get_badvalue( datatype );
275 #endif
276 return retval;
277 } /* pdl_get_pdl_badvalue() */
278
279 #endif
280
281 MODULE = PDL::Core PACKAGE = PDL
282
283
284 # Destroy a PDL - note if a hash do nothing, the $$x{PDL} component
285 # will be destroyed anyway on a separate call
286
287 void
288 DESTROY(sv)
289 SV * sv;
290 PREINIT:
291 pdl *self;
292 CODE:
293 if ( !( (SvROK(sv) && SvTYPE(SvRV(sv)) == SVt_PVHV) ) ) {
294 self = SvPDLV(sv);
295 PDLDEBUG_f(printf("DESTROYING %p\n",(void*)self);)
296 if (self != NULL)
297 pdl_destroy(self);
298 }
299
300 # Return the transformation object or an undef otherwise.
301
302 SV *
303 get_trans(self)
304 pdl *self;
305 CODE:
306 ST(0) = sv_newmortal();
307 if(self->trans) {
308 sv_setref_pv(ST(0), "PDL::Trans", (void*)(self->trans));
309 } else {
310 ST(0) = &PL_sv_undef;
311 }
312
313
314 MODULE = PDL::Core PACKAGE = PDL
315
316 int
317 iscontig(x)
318 pdl* x
319 CODE:
320 RETVAL = 1;
321 pdl_make_physvaffine( x );
322 if PDL_VAFFOK(x) {
323 int i;
324 PDL_Indx inc=1;
325 PDLDEBUG_f(printf("vaff check...\n");)
326 for (i=0;i<x->ndims;i++) {
327 if (PDL_REPRINC(x,i) != inc) {
328 RETVAL = 0;
329 break;
330 }
331 inc *= x->dims[i];
332 }
333 }
334 OUTPUT:
335 RETVAL
336
337 # using "perl" not $^X because that doesn't work on "perl in space"
338 # TODO: switching back to $^X since using "perl" is not a viable fix
339 INCLUDE_COMMAND: $^X -e "require q{./Dev.pm}; PDL::Core::Dev::generate_core_flags()"
340
341 #if 0
342 =begin windows_mmap
343
344 I found this at http://mollyrocket.com/forums/viewtopic.php?p=2529&sid=973b8e0a1e639e3008d7ef05f686c6fa
345 and thougt we might consider using it to make windows mmapping possible.
346
347 -David Mertens
348
349 /*
350 This code was placed in the public domain by the author,
351 Sean Barrett, in November 2007. Do with it as you will.
352 (Seee the page for stb_vorbis or the mollyrocket source
353 page for a longer description of the public domain non-license).
354 */
355
356 #define WIN32_LEAN_AND_MEAN
357 #include <windows.h>
358
359 typedef struct
360 {
361 HANDLE f;
362 HANDLE m;
363 void *p;
364 } SIMPLE_UNMMAP;
365
366 // map 'filename' and return a pointer to it. fill out *length and *un if not-NULL
367 void *simple_mmap(const char *filename, int *length, SIMPLE_UNMMAP *un)
368 {
369 HANDLE f = CreateFile(filename, GENERIC_READ, FILE_SHARE_READ, NULL, OPEN_EXISTING, FILE_ATTRIBUTE_NORMAL, NULL);
370 HANDLE m;
371 void *p;
372 if (!f) return NULL;
373 m = CreateFileMapping(f, NULL, PAGE_READONLY, 0,0, NULL);
374 if (!m) { CloseHandle(f); return NULL; }
375 p = MapViewOfFile(m, FILE_MAP_READ, 0,0,0);
376 if (!p) { CloseHandle(m); CloseHandle(f); return NULL; }
377 if (n) *n = GetFileSize(f, NULL);
378 if (un) {
379 un->f = f;
380 un->m = m;
381 un->p = p;
382 }
383 return p;
384 }
385
386 void simple_unmmap(SIMPLE_UNMMAP *un)
387 {
388 UnmapViewOfFile(un->p);
389 CloseHandle(un->m);
390 CloseHandle(un->f);
391 }
392
393 =end windows_mmap
394
395 =cut
396
397 #endif /* 0 - commented out */
398
399 void
400 set_inplace(self,val)
401 pdl *self;
402 int val;
403 CODE:
404 setflag(self->state,PDL_INPLACE,val);
405
406 IV
407 address(self)
408 pdl *self;
409 CODE:
410 RETVAL = PTR2IV(self);
411 OUTPUT:
412 RETVAL
413
414 pdl *
415 pdl_hard_copy(src)
416 pdl *src;
417
418 pdl *
419 sever(src)
420 pdl *src;
421 CODE:
422 if(src->trans) {
423 pdl_make_physvaffine(src);
424 pdl_destroytransform(src->trans,1);
425 }
426 RETVAL=src;
427 OUTPUT:
428 RETVAL
429
430 int
set_data_by_mmap(it,fname,len,shared,writable,creat,mode,trunc)431 set_data_by_mmap(it,fname,len,shared,writable,creat,mode,trunc)
432 pdl *it
433 char *fname
434 STRLEN len
435 int writable
436 int shared
437 int creat
438 int mode
439 int trunc
440 CODE:
441 #ifdef USE_MMAP
442 int fd;
443 pdl_freedata(it);
444 fd = open(fname,(writable && shared ? O_RDWR : O_RDONLY)|
445 (creat ? O_CREAT : 0),mode);
446 if(fd < 0) {
447 croak("Error opening file");
448 }
449 if(trunc) {
450 int error = ftruncate(fd,0); /* Clear all previous data */
451
452 if(error)
453 {
454 fprintf(stderr,"Failed to set length of '%s' to %d. errno=%d",fname,(int)len,(int)error);
455 croak("set_data_by_mmap: first ftruncate failed");
456 }
457
458 error = ftruncate(fd,len); /* And make it long enough */
459
460 if(error)
461 {
462 fprintf(stderr,"Failed to set length of '%s' to %d. errno=%d",fname,(int)len,(int)error);
463 croak("set_data_by_mmap: second ftruncate failed");
464 }
465 }
466 if(len) {
467 it->data = mmap(0,len,PROT_READ | (writable ?
468 PROT_WRITE : 0),
469 (shared ? MAP_SHARED : MAP_PRIVATE),
470 fd,0);
471 if(!it->data)
472 croak("Error mmapping!");
473 } else {
474 /* Special case: zero-length file */
475 it->data = NULL;
476 }
477 PDLDEBUG_f(printf("PDL::MMap: mapped to %p\n",it->data);)
478 it->state |= PDL_DONTTOUCHDATA | PDL_ALLOCATED;
479 pdl_add_deletedata_magic(it, pdl_delete_mmapped_data, len);
480 close(fd);
481 #else
482 croak("mmap not supported on this architecture");
483 #endif
484 RETVAL = 1;
485 OUTPUT:
486 RETVAL
487
488 int
489 set_state_and_add_deletedata_magic(it,len)
490 pdl *it
491 STRLEN len
492 CODE:
493 it->state |= PDL_DONTTOUCHDATA | PDL_ALLOCATED;
494 pdl_add_deletedata_magic(it, pdl_delete_mmapped_data, len);
495 RETVAL = 1;
496 OUTPUT:
497 RETVAL
498
499 int
500 set_data_by_offset(it,orig,offset)
501 pdl *it
502 pdl *orig
503 STRLEN offset
504 CODE:
505 pdl_freedata(it);
506 it->data = ((char *) orig->data) + offset;
507 it->datasv = orig->sv;
508 (void)SvREFCNT_inc(it->datasv);
509 it->state |= PDL_DONTTOUCHDATA | PDL_ALLOCATED;
510 RETVAL = 1;
511 OUTPUT:
512 RETVAL
513
514 PDL_Indx
515 nelem(x)
516 pdl *x
517 CODE:
518 pdl_make_physdims(x);
519 RETVAL = x->nvals;
520 OUTPUT:
521 RETVAL
522
523 # Convert PDL to new datatype (called by float(), int() etc.)
524
525 # SV *
526 # convert(a,datatype)
527 # pdl* a
528 # int datatype
529 # CODE:
530 # pdl* b;
531 # pdl_make_physical(a);
532 # RETVAL = pdl_copy(a,""); /* Init value to return */
533 # b = SvPDLV(RETVAL); /* Map */
534 # pdl_converttype( &b, datatype, PDL_PERM );
535 # PDLDEBUG_f(printf("converted %d, %d, %d, %d\n",a, b, a->datatype, b->datatype));
536
537 # OUTPUT:
538 # RETVAL
539
540
541 # Call my howbig function
542
543 int
544 howbig_c(datatype)
545 int datatype
546 CODE:
547 RETVAL = pdl_howbig(datatype);
548 OUTPUT:
549 RETVAL
550
551
552 int
553 set_autopthread_targ(i)
554 int i;
555 CODE:
556 RETVAL = i;
557 pdl_autopthread_targ = i;
558 OUTPUT:
559 RETVAL
560
561 int
562 get_autopthread_targ()
563 CODE:
564 RETVAL = pdl_autopthread_targ;
565 OUTPUT:
566 RETVAL
567
568
569 int
570 set_autopthread_size(i)
571 int i;
572 CODE:
573 RETVAL = i;
574 pdl_autopthread_size = i;
575 OUTPUT:
576 RETVAL
577
578 int
579 get_autopthread_size()
580 CODE:
581 RETVAL = pdl_autopthread_size;
582 OUTPUT:
583 RETVAL
584
585 int
586 get_autopthread_actual()
587 CODE:
588 RETVAL = pdl_autopthread_actual;
589 OUTPUT:
590 RETVAL
591
592 MODULE = PDL::Core PACKAGE = PDL::Core
593
594 unsigned int
595 is_scalar_SvPOK(arg)
596 SV* arg;
597 CODE:
598 RETVAL = SvPOK(arg);
599 OUTPUT:
600 RETVAL
601
602
603 int
604 set_debugging(i)
605 int i;
606 CODE:
607 RETVAL = pdl_debugging;
608 pdl_debugging = i;
609 OUTPUT:
610 RETVAL
611
612
613
614 SV *
615 sclr_c(it)
616 pdl* it
617 PREINIT:
618 PDL_Indx nullp = 0;
619 PDL_Indx dummyd = 1;
620 PDL_Indx dummyi = 1;
621 PDL_Anyval result = { -1, 0 };
622 CODE:
623 /* get the first element of a piddle and return as
624 * Perl scalar (autodetect suitable type IV or NV)
625 */
626 pdl_make_physvaffine( it );
627 if (it->nvals < 1)
628 croak("piddle must have at least one element");
629 /* offs = PDL_REPROFFS(it); */
630 /* result = pdl_get_offs(PDL_REPRP(it),offs); */
631 result=pdl_at(PDL_REPRP(it), it->datatype, &nullp, &dummyd,
632 &dummyi, PDL_REPROFFS(it),1);
633 ANYVAL_TO_SV(RETVAL, result);
634
635 OUTPUT:
636 RETVAL
637
638
639 SV *
640 at_c(x,position)
641 pdl* x
642 SV* position
643 PREINIT:
644 PDL_Indx * pos;
645 int npos;
646 int ipos;
647 PDL_Anyval result = { -1, 0 };
648 CODE:
649 pdl_make_physvaffine( x );
650
651 pos = pdl_packdims( position, &npos);
652
653 if (pos == NULL || npos < x->ndims)
654 croak("Invalid position");
655
656 /* allow additional trailing indices
657 * which must be all zero, i.e. a
658 * [3,1,5] piddle is treated as an [3,1,5,1,1,1,....]
659 * infinite dim piddle
660 */
661 for (ipos=x->ndims; ipos<npos; ipos++)
662 if (pos[ipos] != 0)
663 croak("Invalid position");
664
665 result=pdl_at(PDL_REPRP(x), x->datatype, pos, x->dims,
666 (PDL_VAFFOK(x) ? x->vafftrans->incs : x->dimincs), PDL_REPROFFS(x),
667 x->ndims);
668
669 ANYVAL_TO_SV(RETVAL, result);
670
671 OUTPUT:
672 RETVAL
673
674 SV *
675 at_bad_c(x,position)
676 pdl* x
677 SV * position
678 PREINIT:
679 PDL_Indx * pos;
680 int npos;
681 int ipos;
682 int badflag;
683 PDL_Anyval result = { -1, 0 };
684 CODE:
685 pdl_make_physvaffine( x );
686
687 pos = pdl_packdims( position, &npos);
688
689 if (pos == NULL || npos < x->ndims)
690 croak("Invalid position");
691
692 /* allow additional trailing indices
693 * which must be all zero, i.e. a
694 * [3,1,5] piddle is treated as an [3,1,5,1,1,1,....]
695 * infinite dim piddle
696 */
697 for (ipos=x->ndims; ipos<npos; ipos++)
698 if (pos[ipos] != 0)
699 croak("Invalid position");
700
701 result=pdl_at(PDL_REPRP(x), x->datatype, pos, x->dims,
702 (PDL_VAFFOK(x) ? x->vafftrans->incs : x->dimincs), PDL_REPROFFS(x),
703 x->ndims);
704 #if BADVAL
705 badflag = (x->state & PDL_BADVAL) > 0;
706 # if BADVAL_USENAN
707 /* do we have to bother about NaN's? */
708 if ( badflag &&
709 ( ( x->datatype < PDL_F && ANYVAL_EQ_ANYVAL(result, pdl_get_badvalue(x->datatype)) ) ||
710 ( x->datatype == PDL_F && finite(result.value.F) == 0 ) ||
711 ( x->datatype == PDL_D && finite(result.value.D) == 0 )
712 )
713 ) {
714 RETVAL = newSVpvn( "BAD", 3 );
715 } else
716 # else
717 if ( badflag &&
718 ANYVAL_EQ_ANYVAL( result, pdl_get_badvalue( x->datatype ) )
719 ) {
720 RETVAL = newSVpvn( "BAD", 3 );
721 } else
722 # endif
723 #endif
724
725 ANYVAL_TO_SV(RETVAL, result);
726
727 OUTPUT:
728 RETVAL
729
730
731 void
732 list_c(x)
733 pdl *x
734 PREINIT:
735 PDL_Indx *inds;
736 PDL_Indx *incs;
737 PDL_Indx offs;
738 void *data;
739 int ind;
740 int stop = 0;
741 SV *sv;
742 PPCODE:
743 pdl_make_physvaffine( x );
744 inds = pdl_malloc(sizeof(PDL_Indx) * x->ndims); /* GCC -> on stack :( */
745
746 data = PDL_REPRP(x);
747 incs = (PDL_VAFFOK(x) ? x->vafftrans->incs : x->dimincs);
748 offs = PDL_REPROFFS(x);
749 EXTEND(sp,x->nvals);
750 for(ind=0; ind < x->ndims; ind++) inds[ind] = 0;
751 while(!stop) {
752 PDL_Anyval pdl_val = { -1, 0 };
753 pdl_val = pdl_at( data, x->datatype, inds, x->dims, incs, offs, x->ndims);
754 ANYVAL_TO_SV(sv,pdl_val);
755 PUSHs(sv_2mortal(sv));
756 stop = 1;
757 for(ind = 0; ind < x->ndims; ind++)
758 if(++(inds[ind]) >= x->dims[ind])
759 inds[ind] = 0;
760 else
761 {stop = 0; break;}
762 }
763
764 # returns the string 'BAD' if an element is bad
765 #
766
767 SV *
768 listref_c(x)
769 pdl *x
770 PREINIT:
771 PDL_Indx * inds;
772 PDL_Indx * incs;
773 PDL_Indx offs;
774 void *data;
775 int ind;
776 int lind;
777 int stop = 0;
778 AV *av;
779 SV *sv;
780 PDL_Anyval pdl_val = { -1, 0 };
781 PDL_Anyval pdl_badval = { -1, 0 };
782 CODE:
783 #if BADVAL
784 /*
785 # note:
786 # the badvalue is stored in a PDL_Anyval, but that's what pdl_at()
787 # returns
788 */
789
790 int badflag = (x->state & PDL_BADVAL) > 0;
791 # if BADVAL_USENAN
792 /* do we have to bother about NaN's? */
793 if ( badflag && x->datatype < PDL_F ) {
794 pdl_badval = pdl_get_pdl_badvalue( x );
795 }
796 # else
797 if ( badflag ) {
798 pdl_badval = pdl_get_pdl_badvalue( x );
799 }
800 # endif
801 #endif
802
803 pdl_make_physvaffine( x );
804 inds = pdl_malloc(sizeof(PDL_Indx) * x->ndims); /* GCC -> on stack :( */
805 data = PDL_REPRP(x);
806 incs = (PDL_VAFFOK(x) ? x->vafftrans->incs : x->dimincs);
807 offs = PDL_REPROFFS(x);
808 av = newAV();
809 av_extend(av,x->nvals);
810 lind=0;
811 for(ind=0; ind < x->ndims; ind++) inds[ind] = 0;
812 while(!stop) {
813 #if BADVAL
814 pdl_val = pdl_at( data, x->datatype, inds, x->dims, incs, offs, x->ndims );
815 if ( badflag &&
816 # if BADVAL_USENAN
817 ( (x->datatype < PDL_F && ANYVAL_EQ_ANYVAL(pdl_val, pdl_badval)) ||
818 (x->datatype == PDL_F && finite(pdl_val.value.F) == 0) ||
819 (x->datatype == PDL_D && finite(pdl_val.value.D) == 0) )
820 # else
821 ANYVAL_EQ_ANYVAL(pdl_val, pdl_badval)
822 # endif
823 ) {
824 sv = newSVpvn( "BAD", 3 );
825 } else {
826 ANYVAL_TO_SV(sv, pdl_val);
827 }
828 av_store( av, lind, sv );
829 #else
830 pdl_val = pdl_at( data, x->datatype, inds, x->dims, incs, offs, x->ndims );
831 ANYVAL_TO_SV(sv, pdl_val);
832 av_store(av, lind, sv);
833 #endif
834
835 lind++;
836 stop = 1;
837 for(ind = 0; ind < x->ndims; ind++) {
838 if(++(inds[ind]) >= x->dims[ind]) {
839 inds[ind] = 0;
840 } else {
841 stop = 0; break;
842 }
843 }
844 }
845 RETVAL = newRV_noinc((SV *)av);
846 OUTPUT:
847 RETVAL
848
849 void
850 set_c(x,position,value)
851 pdl* x
852 SV* position
853 PDL_Anyval value
854 PREINIT:
855 PDL_Indx * pos;
856 int npos;
857 int ipos;
858 CODE:
859 pdl_make_physvaffine( x );
860
861 pos = pdl_packdims( position, &npos);
862 if (pos == NULL || npos < x->ndims)
863 croak("Invalid position");
864
865 /* allow additional trailing indices
866 * which must be all zero, i.e. a
867 * [3,1,5] piddle is treated as an [3,1,5,1,1,1,....]
868 * infinite dim piddle
869 */
870 for (ipos=x->ndims; ipos<npos; ipos++)
871 if (pos[ipos] != 0)
872 croak("Invalid position");
873
874 pdl_children_changesoon( x , PDL_PARENTDATACHANGED );
875 pdl_set(PDL_REPRP(x), x->datatype, pos, x->dims,
876 (PDL_VAFFOK(x) ? x->vafftrans->incs : x->dimincs), PDL_REPROFFS(x),
877 x->ndims,value);
878 if (PDL_VAFFOK(x))
879 pdl_vaffinechanged(x, PDL_PARENTDATACHANGED);
880 else
881 pdl_changed( x , PDL_PARENTDATACHANGED , 0 );
882
883 BOOT:
884 {
885 #if NVSIZE > 8
886 fprintf(stderr, "Your perl NV has more precision than PDL_Double. There will be loss of floating point precision!\n");
887 #endif
888
889 /* Initialize structure of pointers to core C routines */
890
891 PDL.Version = PDL_CORE_VERSION;
892 PDL.SvPDLV = SvPDLV;
893 PDL.SetSV_PDL = SetSV_PDL;
894 PDL.create = pdl_create;
895 PDL.pdlnew = pdl_external_new;
896 PDL.tmp = pdl_external_tmp;
897 PDL.destroy = pdl_destroy;
898 PDL.null = pdl_null;
899 PDL.copy = pdl_copy;
900 PDL.hard_copy = pdl_hard_copy;
901 PDL.converttype = pdl_converttype;
902 PDL.twod = pdl_twod;
903 PDL.smalloc = pdl_malloc;
904 PDL.howbig = pdl_howbig;
905 PDL.packdims = pdl_packdims;
906 PDL.unpackdims = pdl_unpackdims;
907 PDL.setdims = pdl_setdims;
908 PDL.grow = pdl_grow;
909 PDL.flushcache = NULL;
910 PDL.reallocdims = pdl_reallocdims;
911 PDL.reallocthreadids = pdl_reallocthreadids;
912 PDL.resize_defaultincs = pdl_resize_defaultincs;
913 PDL.get_threadoffsp = pdl_get_threadoffsp;
914 PDL.thread_copy = pdl_thread_copy;
915 PDL.clearthreadstruct = pdl_clearthreadstruct;
916 PDL.initthreadstruct = pdl_initthreadstruct;
917 PDL.startthreadloop = pdl_startthreadloop;
918 PDL.iterthreadloop = pdl_iterthreadloop;
919 PDL.freethreadloop = pdl_freethreadloop;
920 PDL.thread_create_parameter = pdl_thread_create_parameter;
921 PDL.add_deletedata_magic = pdl_add_deletedata_magic;
922
923 PDL.setdims_careful = pdl_setdims_careful;
924 PDL.put_offs = pdl_put_offs;
925 PDL.get_offs = pdl_get_offs;
926 PDL.get = pdl_get;
927 PDL.set_trans_childtrans = pdl_set_trans_childtrans;
928 PDL.set_trans_parenttrans = pdl_set_trans_parenttrans;
929
930 PDL.get_convertedpdl = pdl_get_convertedpdl;
931
932 PDL.make_trans_mutual = pdl_make_trans_mutual;
933 PDL.trans_mallocfreeproc = pdl_trans_mallocfreeproc;
934 PDL.make_physical = pdl_make_physical;
935 PDL.make_physdims = pdl_make_physdims;
936 PDL.make_physvaffine = pdl_make_physvaffine;
937 PDL.pdl_barf = pdl_barf;
938 PDL.pdl_warn = pdl_warn;
939 PDL.allocdata = pdl_allocdata;
940 PDL.safe_indterm = pdl_safe_indterm;
941 PDL.children_changesoon = pdl_children_changesoon;
942 PDL.changed = pdl_changed;
943 PDL.vaffinechanged = pdl_vaffinechanged;
944
945 PDL.NaN_float = union_nan_float.f;
946 PDL.NaN_double = union_nan_double.d;
947 #if BADVAL
948 PDL.propagate_badflag = propagate_badflag;
949 PDL.propagate_badvalue = propagate_badvalue;
950 PDL.get_pdl_badvalue = pdl_get_pdl_badvalue;
951 #include "pdlbadvalinit.c"
952 #endif
953 /*
954 "Publish" pointer to this structure in perl variable for use
955 by other modules
956 */
957 sv_setiv(get_sv("PDL::SHARE",TRUE|GV_ADDMULTI), PTR2IV(&PDL));
958 }
959
960 # make piddle belonging to 'class' and of type 'type'
961 # from avref 'array_ref' which is checked for being
962 # rectangular first
963
964 SV*
965 pdl_avref(array_ref, class, type)
966 SV* array_ref
967 char* class
968 int type
969 PREINIT:
970 AV *dims, *av;
971 int i, depth;
972 int datalevel = -1;
973 SV* psv;
974 pdl* p;
975 CODE:
976 /* make a piddle from a Perl array ref */
977
978 if (!SvROK(array_ref))
979 croak("pdl_avref: not a reference");
980
981
982 if (SvTYPE(SvRV(array_ref)) != SVt_PVAV)
983 croak("pdl_avref: not an array reference");
984
985 // Expand the array ref to a list, and allocate a Perl list to hold the dimlist
986 av = (AV *) SvRV(array_ref);
987 dims = (AV *) sv_2mortal( (SV *) newAV());
988
989 av_store(dims,0,newSViv((IV) av_len(av)+1));
990
991 /* even if we contain nothing depth is one */
992 depth = 1 + av_ndcheck(av,dims,0,&datalevel);
993
994 /* printf("will make type %s\n",class); */
995 /*
996 at this stage start making a piddle and populate it with
997 values from the array (which has already been checked in av_check)
998 */
999 if (strcmp(class,"PDL") == 0) {
1000 p = pdl_from_array(av,dims,type,NULL); /* populate with data */
1001 ST(0) = sv_newmortal();
1002 SetSV_PDL(ST(0),p);
1003 } else {
1004 /* call class->initialize method */
1005 PUSHMARK(SP);
1006 XPUSHs(sv_2mortal(newSVpv(class, 0)));
1007 PUTBACK;
1008 perl_call_method("initialize", G_SCALAR);
1009 SPAGAIN;
1010 psv = POPs;
1011 PUTBACK;
1012 p = SvPDLV(psv); /* and get piddle from returned object */
1013 ST(0) = psv;
1014 pdl_from_array(av,dims,type,p); /* populate ;) */
1015 }
1016
1017 MODULE = PDL::Core PACKAGE = PDL
1018
1019 # pdl_null is created/imported with no PREFIX as pdl_null.
1020 # 'null' is supplied in Core.pm that calls 'initialize' which calls
1021 # the pdl_null here
1022
1023 pdl *
1024 pdl_null(...)
1025
1026
1027 MODULE = PDL::Core PACKAGE = PDL::Core PREFIX = pdl_
1028
1029 int
1030 pdl_pthreads_enabled()
1031
1032 MODULE = PDL::Core PACKAGE = PDL PREFIX = pdl_
1033
1034 int
1035 isnull(self)
1036 pdl *self;
1037 CODE:
1038 RETVAL= !!(self->state & PDL_NOMYDIMS);
1039 OUTPUT:
1040 RETVAL
1041
1042 pdl *
1043 make_physical(self)
1044 pdl *self;
1045 CODE:
1046 pdl_make_physical(self);
1047 RETVAL = self;
1048 OUTPUT:
1049 RETVAL
1050
1051 pdl *
1052 make_physvaffine(self)
1053 pdl *self;
1054 CODE:
1055 pdl_make_physvaffine(self);
1056 RETVAL = self;
1057 OUTPUT:
1058 RETVAL
1059
1060
1061 pdl *
1062 make_physdims(self)
1063 pdl *self;
1064 CODE:
1065 pdl_make_physdims(self);
1066 RETVAL = self;
1067 OUTPUT:
1068 RETVAL
1069
1070 void
1071 pdl_dump(x)
1072 pdl *x;
1073
1074 void
1075 pdl_add_threading_magic(it,nthdim,nthreads)
1076 pdl *it
1077 int nthdim
1078 int nthreads
1079
1080 void
1081 pdl_remove_threading_magic(it)
1082 pdl *it
1083 CODE:
1084 pdl_add_threading_magic(it,-1,-1);
1085
1086 MODULE = PDL::Core PACKAGE = PDL
1087
1088 SV *
initialize(class)1089 initialize(class)
1090 SV *class
1091 PREINIT:
1092 HV *bless_stash;
1093 PPCODE:
1094 if (SvROK(class)) { /* a reference to a class */
1095 bless_stash = SvSTASH(SvRV(class));
1096 } else { /* a class name */
1097 bless_stash = gv_stashsv(class, 0);
1098 }
1099 ST(0) = sv_newmortal();
1100 SetSV_PDL(ST(0),pdl_null()); /* set a null PDL to this SV * */
1101 ST(0) = sv_bless(ST(0), bless_stash); /* bless appropriately */
1102 XSRETURN(1);
1103
1104 SV *
1105 get_dataref(self)
1106 pdl *self
1107 CODE:
1108 if(self->state & PDL_DONTTOUCHDATA) {
1109 croak("Trying to get dataref to magical (mmaped?) pdl");
1110 }
1111 pdl_make_physical(self); /* XXX IS THIS MEMLEAK WITHOUT MORTAL? */
1112 RETVAL = (newRV(self->datasv));
1113 OUTPUT:
1114 RETVAL
1115
1116 int
1117 get_datatype(self)
1118 pdl *self
1119 CODE:
1120 RETVAL = self->datatype;
1121 OUTPUT:
1122 RETVAL
1123
1124 int
upd_data(self)1125 upd_data(self)
1126 pdl *self
1127 PREINIT:
1128 STRLEN n_a;
1129 CODE:
1130 if(self->state & PDL_DONTTOUCHDATA) {
1131 croak("Trying to touch dataref of magical (mmaped?) pdl");
1132 }
1133 self->data = SvPV((SV*)self->datasv,n_a);
1134 XSRETURN(0);
1135
1136 void
1137 set_dataflow_f(self,value)
1138 pdl *self;
1139 int value;
1140 CODE:
1141 if(value)
1142 self->state |= PDL_DATAFLOW_F;
1143 else
1144 self->state &= ~PDL_DATAFLOW_F;
1145
1146 void
1147 set_dataflow_b(self,value)
1148 pdl *self;
1149 int value;
1150 CODE:
1151 if(value)
1152 self->state |= PDL_DATAFLOW_B;
1153 else
1154 self->state &= ~PDL_DATAFLOW_B;
1155
1156 int
1157 getndims(x)
1158 pdl *x
1159 ALIAS:
1160 PDL::ndims = 1
1161 CODE:
1162 pdl_make_physdims(x);
1163 RETVAL = x->ndims;
1164 OUTPUT:
1165 RETVAL
1166
1167 PDL_Indx
1168 getdim(x,y)
1169 pdl *x
1170 int y
1171 ALIAS:
1172 PDL::dim = 1
1173 CODE:
1174 pdl_make_physdims(x);
1175 if (y < 0) y = x->ndims + y;
1176 if (y < 0) croak("negative dim index too large");
1177 if (y < x->ndims)
1178 RETVAL = x->dims[y];
1179 else
1180 RETVAL = 1; /* return size 1 for all other dims */
1181 OUTPUT:
1182 RETVAL
1183
1184 int
1185 getnthreadids(x)
1186 pdl *x
1187 CODE:
1188 pdl_make_physdims(x);
1189 RETVAL = x->nthreadids;
1190 OUTPUT:
1191 RETVAL
1192
1193 int
1194 getthreadid(x,y)
1195 pdl *x
1196 int y
1197 CODE:
1198 RETVAL = x->threadids[y];
1199 OUTPUT:
1200 RETVAL
1201
1202 void
1203 setdims(x,dims_arg)
1204 pdl *x
1205 SV * dims_arg
1206 PREINIT:
1207 PDL_Indx * dims;
1208 int ndims;
1209 int i;
1210 CODE:
1211 {
1212 /* This mask avoids all kinds of subtle dereferencing bugs (CED 11/2015) */
1213 if(x->trans || x->vafftrans || x->children.next ) {
1214 pdl_barf("Can't setdims on a PDL that already has children");
1215 }
1216
1217 /* not sure if this is still necessary with the mask above... (CED 11/2015) */
1218 pdl_children_changesoon(x,PDL_PARENTDIMSCHANGED|PDL_PARENTDATACHANGED);
1219 dims = pdl_packdims(dims_arg,&ndims);
1220 pdl_reallocdims(x,ndims);
1221 for(i=0; i<ndims; i++) x->dims[i] = dims[i];
1222 pdl_resize_defaultincs(x);
1223 x->threadids[0] = ndims;
1224 x->state &= ~PDL_NOMYDIMS;
1225 pdl_changed(x,PDL_PARENTDIMSCHANGED|PDL_PARENTDATACHANGED,0);
1226 }
1227
1228 void
1229 dowhenidle()
1230 CODE:
1231 pdl_run_delayed_magic();
1232 XSRETURN(0);
1233
1234 void
1235 bind(p,c)
1236 pdl *p
1237 SV *c
1238 PROTOTYPE: $&
1239 CODE:
1240 pdl_add_svmagic(p,c);
1241 XSRETURN(0);
1242
1243 void
sethdr(p,h)1244 sethdr(p,h)
1245 pdl *p
1246 SV *h
1247 PREINIT:
1248 HV* hash;
1249 CODE:
1250 if(p->hdrsv == NULL) {
1251 p->hdrsv = &PL_sv_undef; /*(void*) newSViv(0);*/
1252 }
1253
1254 /* Throw an error if we're not either undef or hash */
1255 if ( (h != &PL_sv_undef && h != NULL) &&
1256 ( !SvROK(h) || SvTYPE(SvRV(h)) != SVt_PVHV )
1257 )
1258 croak("Not a HASH reference");
1259
1260 /* Clear the old header */
1261 SvREFCNT_dec(p->hdrsv);
1262
1263 /* Put the new header (or undef) in place */
1264 if(h == &PL_sv_undef || h == NULL)
1265 p->hdrsv = NULL;
1266 else
1267 p->hdrsv = (void*) newRV( (SV*) SvRV(h) );
1268
1269 SV *
1270 hdr(p)
1271 pdl *p
1272 CODE:
1273 pdl_make_physdims(p);
1274
1275 /* Make sure that in the undef case we return not */
1276 /* undef but an empty hash ref. */
1277
1278 if((p->hdrsv==NULL) || (p->hdrsv == &PL_sv_undef)) {
1279 p->hdrsv = (void*) newRV_noinc( (SV*)newHV() );
1280 }
1281
1282 RETVAL = newRV( (SV*) SvRV((SV*)p->hdrsv) );
1283
1284 OUTPUT:
1285 RETVAL
1286
1287 # fhdr(p) is implemented in perl; see Core.pm.PL if you're looking for it
1288 # --CED 9-Feb-2003
1289 #
1290
1291 SV *
1292 gethdr(p)
1293 pdl *p
1294 CODE:
1295 pdl_make_physdims(p);
1296
1297 if((p->hdrsv==NULL) || (p->hdrsv == &PL_sv_undef)) {
1298 RETVAL = &PL_sv_undef;
1299 } else {
1300 RETVAL = newRV( (SV*) SvRV((SV*)p->hdrsv) );
1301 }
1302
1303 OUTPUT:
1304 RETVAL
1305
1306 void
1307 set_datatype(a,datatype)
1308 pdl *a
1309 int datatype
1310 CODE:
1311 pdl_make_physical(a);
1312 if(a->trans)
1313 pdl_destroytransform(a->trans,1);
1314 /* if(! (a->state && PDL_NOMYDIMS)) { */
1315 pdl_converttype( &a, datatype, PDL_PERM );
1316 /* } */
1317
1318 void
1319 threadover_n(...)
1320 PREINIT:
1321 int npdls;
1322 SV *sv;
1323 CODE:
1324 {
1325 npdls = items - 1;
1326 if(npdls <= 0)
1327 croak("Usage: threadover_n(pdl[,pdl...],sub)");
1328 {
1329 int i,sd;
1330 pdl **pdls = malloc(sizeof(pdl *) * npdls);
1331 PDL_Indx *realdims = malloc(sizeof(PDL_Indx) * npdls);
1332 pdl_thread pdl_thr;
1333 SV *code = ST(items-1);
1334 for(i=0; i<npdls; i++) {
1335 pdls[i] = SvPDLV(ST(i));
1336 /* XXXXXXXX Bad */
1337 pdl_make_physical(pdls[i]);
1338 realdims[i] = 0;
1339 }
1340 PDL_THR_CLRMAGIC(&pdl_thr);
1341 pdl_initthreadstruct(0,pdls,realdims,realdims,npdls,NULL,&pdl_thr,NULL, 1);
1342 pdl_startthreadloop(&pdl_thr,NULL,NULL);
1343 sd = pdl_thr.ndims;
1344 do {
1345 dSP;
1346 PUSHMARK(sp);
1347 EXTEND(sp,items);
1348 PUSHs(sv_2mortal(newSViv((sd-1))));
1349 for(i=0; i<npdls; i++) {
1350 PDL_Anyval pdl_val = { -1, 0 };
1351 pdl_val = pdl_get_offs(pdls[i],pdl_thr.offs[i]);
1352 ANYVAL_TO_SV(sv, pdl_val);
1353 PUSHs(sv_2mortal(sv));
1354 }
1355 PUTBACK;
1356 perl_call_sv(code,G_DISCARD);
1357 } while( (sd = pdl_iterthreadloop(&pdl_thr,0)) );
1358 pdl_freethreadloop(&pdl_thr);
1359 free(pdls);
1360 free(realdims);
1361 }
1362 }
1363
1364 void
1365 threadover(...)
1366 PREINIT:
1367 int npdls;
1368 int targs;
1369 int nothers = -1;
1370 CODE:
1371 {
1372 targs = items - 4;
1373 if (items > 0) nothers = SvIV(ST(0));
1374 if(targs <= 0 || nothers < 0 || nothers >= targs)
1375 croak("Usage: threadover(nothers,pdl[,pdl...][,otherpars..],realdims,creating,sub)");
1376 npdls = targs-nothers;
1377 {
1378 int i,nd1,nd2,dtype=0;
1379 PDL_Indx j,nc=npdls;
1380 SV* rdimslist = ST(items-3);
1381 SV* cdimslist = ST(items-2);
1382 SV *code = ST(items-1);
1383 pdl_thread pdl_thr;
1384 pdl **pdls = malloc(sizeof(pdl *) * npdls);
1385 pdl **child = malloc(sizeof(pdl *) * npdls);
1386 SV **csv = malloc(sizeof(SV *) * npdls);
1387 SV **dims = malloc(sizeof(SV *) * npdls);
1388 SV **incs = malloc(sizeof(SV *) * npdls);
1389 SV **others = malloc(sizeof(SV *) * nothers);
1390 PDL_Indx *creating = pdl_packint(cdimslist,&nd2);
1391 PDL_Indx *realdims = pdl_packint(rdimslist,&nd1);
1392 CHECKP(pdls); CHECKP(child); CHECKP(dims);
1393 CHECKP(incs); CHECKP(csv);
1394
1395 if (nd1 != npdls || nd2 < npdls)
1396 croak("threadover: need one realdim and creating flag "
1397 "per pdl!");
1398 for(i=0; i<npdls; i++) {
1399 pdls[i] = SvPDLV(ST(i+1));
1400 if (creating[i])
1401 nc += realdims[i];
1402 else {
1403 pdl_make_physical(pdls[i]); /* is this what we want?XXX */
1404 dtype = PDLMAX(dtype,pdls[i]->datatype);
1405 }
1406 }
1407 for (i=npdls+1; i<=targs; i++)
1408 others[i-npdls-1] = ST(i);
1409 if (nd2 < nc)
1410 croak("Not enough dimension info to create pdls");
1411 #ifdef DEBUG_PTHREAD
1412 for (i=0;i<npdls;i++) { /* just for debugging purposes */
1413 printf("pdl %d Dims: [",i);
1414 for (j=0;j<realdims[i];j++)
1415 printf("%d ",pdls[i]->dims[j]);
1416 printf("] Incs: [");
1417 for (j=0;j<realdims[i];j++)
1418 printf("%d ",PDL_REPRINC(pdls[i],j));
1419 printf("]\n");
1420 }
1421 #endif
1422 PDL_THR_CLRMAGIC(&pdl_thr);
1423 pdl_initthreadstruct(0,pdls,realdims,creating,npdls,
1424 NULL,&pdl_thr,NULL, 1);
1425 for(i=0, nc=npdls; i<npdls; i++) /* create as necessary */
1426 if (creating[i]) {
1427 PDL_Indx *cp = creating+nc;
1428 pdls[i]->datatype = dtype;
1429 pdl_thread_create_parameter(&pdl_thr,i,cp,0);
1430 nc += realdims[i];
1431 pdl_make_physical(pdls[i]);
1432 PDLDEBUG_f(pdl_dump(pdls[i]));
1433 /* And make it nonnull, now that we've created it */
1434 pdls[i]->state &= (~PDL_NOMYDIMS);
1435 }
1436 pdl_startthreadloop(&pdl_thr,NULL,NULL);
1437 for(i=0; i<npdls; i++) { /* will the SV*'s be properly freed? */
1438 dims[i] = newRV(pdl_unpackint(pdls[i]->dims,realdims[i]));
1439 incs[i] = newRV(pdl_unpackint(PDL_VAFFOK(pdls[i]) ?
1440 pdls[i]->vafftrans->incs: pdls[i]->dimincs,realdims[i]));
1441 /* need to make sure we get the vaffine (grand)parent */
1442 if (PDL_VAFFOK(pdls[i]))
1443 pdls[i] = pdls[i]->vafftrans->from;
1444 child[i]=pdl_null();
1445 /* instead of pdls[i] its vaffine parent !!!XXX */
1446 PDL.affine_new(pdls[i],child[i],pdl_thr.offs[i],dims[i],
1447 incs[i]);
1448 pdl_make_physical(child[i]); /* make sure we can get at
1449 the vafftrans */
1450 csv[i] = sv_newmortal();
1451 SetSV_PDL(csv[i], child[i]); /* pdl* into SV* */
1452 }
1453 do { /* the actual threadloop */
1454 pdl_trans_affine *traff;
1455 dSP;
1456 PUSHMARK(sp);
1457 EXTEND(sp,npdls);
1458 for(i=0; i<npdls; i++) {
1459 /* just twiddle the offset - quick and dirty */
1460 /* we must twiddle both !! */
1461 traff = (pdl_trans_affine *) child[i]->trans;
1462 traff->offs = pdl_thr.offs[i];
1463 child[i]->vafftrans->offs = pdl_thr.offs[i];
1464 child[i]->state |= PDL_PARENTDATACHANGED;
1465 PUSHs(csv[i]);
1466 }
1467 for (i=0; i<nothers; i++)
1468 PUSHs(others[i]); /* pass the OtherArgs onto the stack */
1469 PUTBACK;
1470 perl_call_sv(code,G_DISCARD);
1471 } while (pdl_iterthreadloop(&pdl_thr,0));
1472 pdl_freethreadloop(&pdl_thr);
1473 free(pdls); /* should all these be done with pdl_malloc */
1474 free(dims); /* in case the sub barfs ? XXXX */
1475 free(child);
1476 free(csv);
1477 free(incs);
1478 free(others);
1479 }
1480 }
1481