#
# Create pdlcore.c
# - needed for bad-value handling in whichdatatype
#
use strict;
use Config;
use File::Basename qw(&basename &dirname);
require 'Dev.pm'; PDL::Core::Dev->import;
use vars qw( %PDL_DATATYPES );
# check for bad value support
use vars qw( $bvalflag $usenan );
require "badsupport.p";
require 'Types.pm';
PDL::Types->import(':All');
# This forces PL files to create target in same directory as PL file.
# This is so that make depend always knows where to find PL derivatives.
chdir(dirname($0));
my $file;
($file = basename($0)) =~ s/\.PL$//;
$file =~ s/\.pl$//
if ($Config{'osname'} eq 'VMS' or
$Config{'osname'} eq 'OS2'); # "case-forgiving"
if ( $bvalflag ) {
print "Extracting $file (WITH bad value support)\n";
} else {
print "Extracting $file (NO bad value support)\n";
}
open OUT,">$file" or die "Can't create $file: $!";
chmod 0644, $file;
print OUT <<"!WITH!SUBS!";
/* pdlcore.c - generated automatically by pdlcore.c.PL */
/* - bad value support = $bvalflag */
!WITH!SUBS!
;
print OUT <<'!HEADER!'
#undef FOODEB
#define PDL_CORE /* For certain ifdefs */
#include "pdl.h" /* Data structure declarations */
#include "pdlcore.h" /* Core declarations */
/*** Turn on definitions to print lots of fencepost information in the constructor ***/
//#define DEBUG_SETAV_TYPE 1
//#define DEBUG_KLUDGE_COPY 1
/***************
* Paranoid check is commented out because SvPOK breaks some kinds of perl scalars.
* (blessed scalars like '$#$foo' get set to zero in perl 5.8.x)
*/
/* #define sv_undef(sv) ( (!(sv) || ((sv)==&PL_sv_undef) ) || !(SvNIOK(sv) || SvPOK(sv) || SvROK(sv)) ) */
#define sv_undef(sv) ( (!(sv) || ((sv)==&PL_sv_undef)) || !(SvNIOK(sv) || (SvTYPE(sv)==SVt_PVMG) || SvPOK(sv) || SvROK(sv)))
!HEADER!
;
if($Config{cc} eq 'cl') {
# _finite in CV++ 4.0
print OUT <<'FOO';
#define finite _finite
#include <float.h>
FOO
;
}
my $finite_inc;
foreach my $inc ( qw/ math.h ieeefp.h / )
{
if ( trylink ("finite: $inc", "#include <$inc>", 'finite(3.2);', '' ) ) {
$finite_inc = $inc;
last;
}
}
if ( defined $finite_inc )
{
print OUT "#include <$finite_inc>\n";
}
else
{
print OUT <<'!NO!SUBS!'
#ifndef finite
#ifdef isfinite
#define finite isfinite
#else
#define finite(a) (((a) * 0) == (0))
#endif
#endif
!NO!SUBS!
}
print OUT <<'!NO!SUBS!'
static SV *getref_pdl(pdl *it) {
SV *newref;
if(!it->sv) {
SV *ref;
HV *stash = gv_stashpv("PDL",TRUE);
SV *psv = newSViv(PTR2IV(it));
it->sv = psv;
newref = newRV_noinc(it->sv);
(void)sv_bless(newref,stash);
} else {
newref = newRV_inc(it->sv);
SvAMAGIC_on(newref);
}
return newref;
}
void SetSV_PDL ( SV *sv, pdl *it ) {
SV *newref = getref_pdl(it); /* YUCK!!!! */
sv_setsv(sv,newref);
SvREFCNT_dec(newref);
}
/* Size of data type information */
int pdl_howbig (int datatype) {
switch (datatype) {
!NO!SUBS!
;
# generate the cases for the various types
for my $type (typesrtkeys()) {
my ($sym,$ctype) = map {typefld($type,$_)} qw/sym ctype/;
print OUT << "!WITH!SUBS!";
case $sym:
return sizeof($ctype);
!WITH!SUBS!
}
print OUT <<'!NO!SUBS!';
default:
croak("Unknown datatype code = %d",datatype);
}
}
/* Check minimum datatype required to represent number */
/* Microsoft compilers do some unbelievable things - hence
some #ifdef's inserted by Sisyphus */
#if defined _MSC_VER && _MSC_VER < 1400
#define TESTTYPE(b,a) {a foo = nv; a bar = foo; foo += 0; if(nv == bar) return b;}
#else
#define TESTTYPE(b,a) {a foo = nv; if(nv == foo) return b;}
#endif
#if defined _MSC_VER && _MSC_VER < 1400
#pragma optimize("", off)
#endif
int pdl_whichdatatype (double nv) {
!NO!SUBS!
# generate the cases for the various types
for my $type (typesrtkeys()) {
my ($sym,$ctype) = map {typefld($type,$_)} qw/sym ctype/;
print OUT << "!WITH!SUBS!";
TESTTYPE($sym,$ctype)
!WITH!SUBS!
}
print OUT <<'!NO!SUBS!';
if( !finite(nv) ) { return PDL_D; }
croak("Something's gone wrong: %lf cannot be converted by whichdatatype",
nv);
}
/* Check minimum, at least float, datatype required to represent number */
int pdl_whichdatatype_double (double nv) {
TESTTYPE(PDL_F,PDL_Float)
TESTTYPE(PDL_D,PDL_Double)
if( !finite(nv) ) { return PDL_D; }
croak("Something's gone wrong: %lf cannot be converted by whichdatatype_double",
nv);
}
#if defined _MSC_VER && _MSC_VER < 1400
#pragma optimize("", on)
#endif
/* Make a scratch data existence for a pdl */
void pdl_makescratchhash(pdl *ret,double data, int datatype) {
STRLEN n_a;
HV *hash;
SV *dat; PDL_Indx fake[1];
/* Compress to smallest available type. This may have strange
results sometimes :( */
ret->datatype = datatype;
ret->data = pdl_malloc(pdl_howbig(ret->datatype)); /* Wasteful */
dat = newSVpv(ret->data,pdl_howbig(ret->datatype));
ret->data = SvPV(dat,n_a);
ret->datasv = dat;
#ifdef FOO
/* Refcnt should be 1 already... */
SvREFCNT_inc(ret->datasv); /* XXX MEMLEAK */
#endif
/* This is an important point: it makes this whole piddle mortal
* so destruction will happen at the right time.
* If there are dangling references, pdlapi.c knows not to actually
* destroy the C struct. */
sv_2mortal(getref_pdl(ret));
pdl_setdims(ret, fake, 0); /* However, there are 0 dims in scalar */
ret->nvals = 1;
/* NULLs should be ok because no dimensions. */
pdl_set(ret->data, ret->datatype, NULL, NULL, NULL, 0, 0, data);
}
/*
"Convert" a perl SV into a pdl (alright more like a mapping as
the data block is not actually copied in the most common case
of a single scalar) scalars are automatically converted to PDLs.
*/
pdl* SvPDLV ( SV* sv ) {
pdl* ret;
PDL_Indx fake[1];
SV *sv2;
if ( !SvROK(sv) ) {
/* The scalar is not a ref, so we can use direct conversion. */
SV *dat;
double data;
int datatype;
ret = pdl_new(); /* Scratch pdl */
/* Scratch hash for the pdl :( - slow but safest. */
/* handle undefined values */
if( sv_undef(sv) ) {
sv = get_sv("PDL::undefval",1);
if(SvIV(get_sv("PDL::debug",1))){
fprintf(stderr,"Warning: SvPDLV converted undef to $PDL::undefval (%g).\n",SvNV(sv));
}
}
/* Figure datatype to use */
if ( !SvIOK(sv) && SvNOK(sv) && SvNIOK(sv) ) {/* Perl Double (e.g. 2.0) */
data = SvNV(sv);
!NO!SUBS!
##########
# If badvals and usenan are true, we default NaNs to type double.
if ( $bvalflag and $usenan ) {
print OUT q{
/*
* default NaN/Infs to double
*
*/
if ( finite(data) == 0 ) {
datatype = PDL_D;
} else {
datatype = pdl_whichdatatype_double(data);
}
}
} else {
print OUT q{
datatype = pdl_whichdatatype_double(data);
}
} # if: $bvalflag
print OUT <<'!NO!SUBS!';
} /* end of Perl Double case for data type */
else { /* Perl Int (e.g. 2) */
data = SvNV(sv);
datatype = pdl_whichdatatype(data);
}
pdl_makescratchhash(ret,data,datatype);
return ret;
} /* End of scalar case */
/* If execution reaches here, then sv is NOT a scalar
* (i.e. it is a ref).
*/
if(SvTYPE(SvRV(sv)) == SVt_PVHV) {
HV *hash = (HV*)SvRV(sv);
SV **svp = hv_fetch(hash,"PDL",3,0);
if(svp == NULL) {
croak("Hash given as a pdl - but not {PDL} key!");
}
if(*svp == NULL) {
croak("Hash given as a pdl - but not {PDL} key (*svp)!");
}
/* This is the magic hook which checks to see if {PDL}
is a code ref, and if so executes it. It should
return a standard piddle. This allows
all kinds of funky objects to be derived from PDL,
and allow normal PDL functions to still work so long
as the {PDL} code returns a standard piddle on
demand - KGB */
if (SvROK(*svp) && SvTYPE(SvRV(*svp)) == SVt_PVCV) {
dSP;
int count;
ENTER ;
SAVETMPS ;
PUSHMARK(sp) ;
count = perl_call_sv(*svp, G_SCALAR|G_NOARGS);
SPAGAIN ;
if (count != 1)
croak("Execution of PDL structure failed to return one value\n") ;
sv=newSVsv(POPs);
PUTBACK ;
FREETMPS ;
LEAVE ;
}
else {
sv = *svp;
}
if(SvGMAGICAL(sv)) {
mg_get(sv);
}
if ( !SvROK(sv) ) { /* Got something from a hash but not a ref */
croak("Hash given as pdl - but PDL key is not a ref!");
}
}
if(SvTYPE(SvRV(sv)) == SVt_PVAV) {
/* This is similar to pdl_avref in Core.xs.PL -- we do the same steps here. */
AV *dims, *av;
int i, depth;
int datalevel = -1;
pdl *p;
av = (AV *)SvRV(sv);
dims = (AV *)sv_2mortal((SV *)newAV());
av_store(dims,0,newSViv( (IV) av_len(av)+1 ) );
/* Pull sizes using av_ndcheck */
depth = 1 + av_ndcheck(av,dims,0,&datalevel);
return pdl_from_array(av,dims,PDL_D,NULL);
} /* end of AV code */
if (SvTYPE(SvRV(sv)) != SVt_PVMG)
croak("Error - tried to use an unknown data structure as a PDL");
else if( !( sv_derived_from( sv, "PDL") ) )
croak("Error - tried to use an unknown Perl object type as a PDL");
sv2 = (SV*) SvRV(sv);
/* Return the pdl * pointer */
ret = INT2PTR(pdl *, SvIV(sv2));
/* Final check -- make sure it has the right magic number */
if(ret->magicno != PDL_MAGICNO) {
croak("Fatal error: argument is probably not a piddle, or\
magic no overwritten. You're in trouble, guv: %p %p %lu\n",sv2,ret,ret->magicno);
}
return ret;
}
/* Make a new pdl object as a copy of an old one and return - implement by
callback to perl method "copy" or "new" (for scalar upgrade) */
SV* pdl_copy( pdl* a, char* option ) {
SV* retval;
char meth[20];
dSP ; int count ;
retval = newSVpv("",0); /* Create the new SV */
ENTER ; SAVETMPS ; PUSHMARK(sp) ;
/* Push arguments */
#ifdef FOOBAR
if (sv_isobject((SV*)a->hash)) {
#endif
XPUSHs(sv_2mortal(getref_pdl(a)));
strcpy(meth,"copy");
XPUSHs(sv_2mortal(newSVpv(option, 0))) ;
#ifdef FOOBAR
}
else{
XPUSHs(perl_get_sv("PDL::name",FALSE)); /* Default object */
XPUSHs(sv_2mortal(getref_pdl(a)));
strcpy(meth,"new");
}
#endif
PUTBACK ;
count = perl_call_method(meth, G_SCALAR); /* Call Perl */
SPAGAIN;
if (count !=1)
croak("Error calling perl function\n");
sv_setsv( retval, POPs ); /* Save the perl returned value */
PUTBACK ; FREETMPS ; LEAVE ;
return retval;
}
/* Pack dims array - returns dims[] (pdl_malloced) and ndims */
PDL_Indx* pdl_packdims ( SV* sv, int *ndims ) {
SV* bar;
AV* array;
int i;
PDL_Indx *dims;
if (!(SvROK(sv) && SvTYPE(SvRV(sv))==SVt_PVAV)) /* Test */
return NULL;
array = (AV *) SvRV(sv); /* dereference */
*ndims = (int) av_len(array) + 1; /* Number of dimensions */
dims = (PDL_Indx *) pdl_malloc( (*ndims) * sizeof(*dims) ); /* Array space */
if (dims == NULL)
croak("Out of memory");
for(i=0; i<(*ndims); i++) {
bar = *(av_fetch( array, i, 0 )); /* Fetch */
dims[i] = (PDL_Indx) SvIV(bar);
}
return dims;
}
/* unpack dims array into PDL SV* */
void pdl_unpackdims ( SV* sv, PDL_Indx *dims, int ndims ) {
AV* array;
HV* hash;
int i;
hash = (HV*) SvRV( sv );
array = newAV();
(void)hv_store(hash, "Dims", strlen("Dims"), newRV( (SV*) array), 0 );
if (ndims==0 )
return;
for(i=0; i<ndims; i++)
av_store( array, i, newSViv( (IV)dims[i] ) );
}
PDL_Indx pdl_safe_indterm( PDL_Indx dsz, PDL_Indx at, char *file, int lineno)
{
if (at >= 0 && at < dsz) return at;
pdl_barf("access [%d] out of range [0..%d] (inclusive) at %s line %d",
at, dsz-1, file?file:"?", lineno);
return at; // This can never happen - placed to avoid a compiler warning.
}
/*
pdl_malloc - utility to get temporary memory space. Uses
a mortal *SV for this so it is automatically freed when the current
context is terminated without having to call free(). Naughty but
nice!
*/
void* pdl_malloc ( STRLEN nbytes ) {
STRLEN n_a;
SV* work;
work = sv_2mortal(newSVpv("", 0));
SvGROW( work, nbytes);
return (void *) SvPV(work, n_a);
}
/*********** Stuff for barfing *************/
/*
This routine barfs/warns in a thread-safe manner. If we're in the main thread,
this calls the perl-level barf/warn. If in a worker thread, we save the
message to barf/warn in the main thread later
*/
static void pdl_barf_or_warn(const char* pat, int iswarn, va_list* args)
{
/* If we're in a worker thread, we queue the
* barf/warn for later, and exit the thread ...
*/
if( pdl_pthread_barf_or_warn(pat, iswarn, args) )
return;
/* ... otherwise we fall through and barf by calling
* the perl-level PDL::barf() or PDL::cluck()
*/
{ /* scope block for C89 compatibility */
SV * sv;
dSP;
ENTER;
SAVETMPS;
PUSHMARK(SP);
sv = sv_2mortal(newSV(0));
sv_vsetpvfn(sv, pat, strlen(pat), args, Null(SV**), 0, Null(bool*));
va_end(*args);
XPUSHs(sv);
PUTBACK;
if(iswarn) call_pv("PDL::cluck", G_DISCARD);
else call_pv("PDL::barf", G_DISCARD);
FREETMPS;
LEAVE;
} /* end C89 compatibility scope block */
}
#define GEN_PDL_BARF_OR_WARN_I_STDARG(type, iswarn) \
void pdl_##type(const char* pat, ...) \
{ \
va_list args; \
va_start(args, pat); \
pdl_barf_or_warn(pat, iswarn, &args); \
}
#define GEN_PDL_BARF_OR_WARN_LEGACY(type, iswarn) \
void pdl_##type(pat, va_alist) \
char *pat; \
va_dcl \
{ \
va_list args; \
va_start(args); \
pdl_barf_or_warn(pat, iswarn, &args); \
}
#ifdef I_STDARG
GEN_PDL_BARF_OR_WARN_I_STDARG(barf, 0)
GEN_PDL_BARF_OR_WARN_I_STDARG(warn, 1)
#else
GEN_PDL_BARF_OR_WARN_LEGACY(barf, 0)
GEN_PDL_BARF_OR_WARN_LEGACY(warn, 1)
#endif
/**********************************************************************
*
* CONSTRUCTOR/INGESTION HELPERS
*
* The following routines assist with the permissive constructor,
* which is designed to build a PDL out of basically anything thrown at it.
*
* They are all called by pdl_avref in Core.xs, which in turn is called by the constructors
* in Core.pm.PL.
*
*/
/******************************
* av_ndcheck -
* traverse a Perl array ref recursively, following down any number of
* levels of references, and generate a minimal PDL dim list that can
* encompass them all according to permissive-constructor rules.
*
* Scalars, array refs, and PDLs may be mixed in the incoming AV.
*
* The routine works out the dimensions of a corresponding
* piddle (in the AV dims) in reverse notation (vs PDL conventions).
*
* It does not enforce a rectangular array, the idea being that
* omitted values will be set to zero in the resulting piddle,
* i.e. we can make piddles from 'sparse' array refs.
*
* Empty PDLs are treated like any other dimension -- i.e. their
* 0-length dimensions are thrown into the mix just like nonzero
* dimensions would be.
*
* The possible presence of empty PDLs forces us to pad out dimensions
* to unity explicitly in cases like
* [ Empty[2x0x2], 5 ]
* where simple parsing would yield a dimlist of
* [ 2,0,2,2 ]
* which is still Empty.
*/
PDL_Indx av_ndcheck(AV* av, AV* dims, int level, int *datalevel)
{
PDL_Indx i, len, oldlen;
int newdepth, depth = 0;
int n_scalars = 0;
SV *el, **elp;
pdl *pdl; /* Stores PDL argument */
if(dims==NULL) {
pdl_barf("av_ndcheck - got a null dim array! This is a bug in PDL.");
}
/* Start with a clean slate */
if(level==0) {
av_clear(dims);
}
len = av_len(av); /* Loop over elements of the AV */
for (i=0; i<= len; i++) {
newdepth = 0; /* Each element - find depth */
elp = av_fetch(av,i,0);
el = elp ? *elp : 0; /* Get the ith element */
if (el && SvROK(el)) { /* It is a reference */
if (SvTYPE(SvRV(el)) == SVt_PVAV) { /* It is an array reference */
/* Recurse to find depth inside the array reference */
newdepth = 1 + av_ndcheck((AV *) SvRV(el), dims, level+1, datalevel);
} else if ( (pdl = SvPDLV(el)) ) {
/* It is a PDL - walk down its dimension list, exactly as if it
* were a bunch of nested array refs. We pull the ndims and dims
* fields out to local variables so that nulls can be treated specially.
*/
int j;
short pndims;
PDL_Indx *pdims;
PDL_Indx pnvals;
#ifdef DEBUG_SETAV_TYPE
printf("av_ndcheck - found a PDL....\n");
#endif
pdl_make_physdims(pdl);
pndims = pdl->ndims;
pdims = pdl->dims;
pnvals = pdl->nvals;
#ifdef DEBUG_SETAV_TYPE
{
printf("av_ndcheck: nvals is %d; ndims is %d; dims are (",pdl->nvals, pdl->ndims);
for(j=0;j<pdl->ndims; j++) {
printf("%s%d",j?", ":"",pdl->dims[j]);
}
printf("); level is %d\n",level);
}
#endif
for(j=0;j<pndims;j++) {
int jl = pndims-j+level;
PDL_Indx siz = pdims[j];
if( av_len(dims) >= jl &&
av_fetch(dims,jl,0) != NULL &&
SvIOK(*(av_fetch(dims,jl,0)))) {
/* We have already found something that specifies this dimension -- so */
/* we keep the size if possible, or enlarge if necessary. */
oldlen=(PDL_Indx)SvIV(*(av_fetch(dims,jl,0)));
if(siz > oldlen) {
sv_setiv(*(av_fetch(dims,jl,0)),(IV)(pdims[j]));
}
} else {
/* Breaking new dimensional ground here -- if this is the first element */
/* in the arg list, then we can keep zero elements -- but if it is not */
/* the first element, we have to pad zero dims to unity (because the */
/* prior object had implicit size of 1 in all implicit dimensions) */
av_store(dims, jl, newSViv((IV)(siz?siz:(i?1:0))));
}
}
/* We have specified all the dims in this PDL. Now pad out the implicit */
/* dims of size unity, to wipe out any dims of size zero we have already */
/* marked. */
for(j=pndims+1; j <= av_len(dims); j++) {
SV **svp = av_fetch(dims,j,0);
if(!svp){
av_store(dims, j, newSViv((IV)1));
} else if( (int)SvIV(*svp) == 0 ) {
sv_setiv(*svp, (IV)1);
}
}
newdepth= pndims;
} else {
croak("av_ndcheck: non-array, non-PDL ref in structure\n\t(this is usually a problem with a pdl() call)");
}
} else {
/* got a scalar (not a ref) */
n_scalars++;
}
if (newdepth > depth)
depth = newdepth;
}
len++; // convert from funky av_len return value to real count
#ifdef DEBUG_SETAV_TYPE
{
int i,dim;
if(dims != NULL) {
dim = av_len(dims) + 1;
} else {
dim =0 ;
}
printf("av_ndcheck: depth=%d, length is %d; derived dim list is [",depth, dim);
fflush(stdout);
for( i=0; i<dim; i++ ) {
SV **svp = av_fetch(dims, i, 0);
PDL_Indx k;
if(svp != NULL) {
k = SvIV (*svp);
} else {
k = -1;
}
printf(" %d",k);
}
printf(" ]\n");
}
#endif
if (av_len(dims) >= level && av_fetch(dims, level, 0) != NULL
&& SvIOK(*(av_fetch(dims, level, 0)))) {
oldlen = (PDL_Indx) SvIV(*(av_fetch(dims, level, 0)));
if (len > oldlen)
sv_setiv(*(av_fetch(dims, level, 0)), (IV) len);
}
else
av_store(dims,level,newSViv((IV) len));
/* We found at least one element -- so pad dims to unity at levels earlier than this one */
#ifdef DEBUG_SETAV_TYPE
printf("n_scalars=%d\n",n_scalars);
#endif
if(n_scalars) {
for(i=0;i<level;i++) {
SV **svp = av_fetch(dims, i, 0);
if(!svp) {
av_store(dims, i, newSViv((IV)1));
} else if( (PDL_Indx)SvIV(*svp) == 0) {
sv_setiv(*svp, (IV)1);
}
}
for(i=level+1; i <= av_len(dims); i++) {
SV **svp = av_fetch(dims, i, 0);
if(!svp) {
av_store(dims, i, newSViv((IV)1));
} else if( (PDL_Indx)SvIV(*svp) == 0) {
sv_setiv(*svp, (IV)1);
}
}
}
#ifdef DEBUG_SETAV_TYPE
{
int i,dim;
if(dims != NULL) {
dim = av_len(dims) + 1;
} else {
dim =0 ;
}
printf("av_ndcheck: depth=%d, length is %d; derived dim list is [",depth, dim);
fflush(stdout);
for( i=0; i<dim; i++ ) {
SV **svp = av_fetch(dims, i, 0);
PDL_Indx k;
if(svp != NULL) {
k = SvIV (*svp);
} else {
k = -1;
}
printf(" %d",k);
}
printf(" ]\n");
}
#endif
return depth;
}
pdl* pdl_from_array(AV* av, AV* dims, int type, pdl* p)
{
int ndims, i, level=0;
PDL_Indx *pdims;
double undefval;
ndims = av_len(dims)+1;
pdims = (PDL_Indx *) pdl_malloc( (ndims) * sizeof(*pdims) );
for (i=0; i<ndims; i++) {
pdims[i] = SvIV(*(av_fetch(dims, ndims-1-i, 0))); /* reverse order */
}
#ifdef DEBUG_SETAV_TYPE
{
int ii;
printf("pdl_from_array: dim list is: (");
for(ii=0; ii<ndims; ii++) {
printf("%s%d",ii?", ":"",pdims[ii]);
}
printf(")\n");
}
#endif
if (p == NULL)
p = pdl_new();
pdl_setdims (p, pdims, ndims);
p->datatype = type;
pdl_allocdata (p);
pdl_make_physical(p);
/* this one assigns the data */
{
/******
* Copy the undefval to fill empty spots in the piddle...
*/
SV *sv = get_sv("PDL::undefval",0);
undefval = ((!sv) || (sv==&PL_sv_undef)) ? 0 : (double)SvNV(sv);
}
switch (type) {
!NO!SUBS!
##########
# Perl snippet autogenerates switch statement to distribute
# pdl_setav calls...
#
for my $type(sort keys %PDL_DATATYPES){
my $t2 = $PDL_DATATYPES{$type};
$t2 =~ s/PDL_//;
print OUT <<"!WITH!SUBS!";
case $type:
pdl_setav_$t2(p->data,av,pdims,ndims,level, undefval);
break;
!WITH!SUBS!
}
#
# Back to your regularly scheduled C code emission...
########
print OUT <<'!NO!SUBS!';
default:
croak("pdl_from_array: internal error: got type %d",type);
break;
}
p->state &= ~PDL_NOMYDIMS;
return p;
}
!NO!SUBS!
######################################################################
# these are helper functions for fast assignment from array refs
# mainly used by pdl_avref in Core.xs which implements converting
# array refs to piddles
for my $in ( keys %PDL_DATATYPES ) {
(my $type = $PDL_DATATYPES{$in}) =~ s/^PDL_//;
print OUT <<"!WITH!SUBS!";
/*
* pdl_kludge_copy - copy a PDL into a part of a being-formed PDL.
* It is only used by pdl_setav, to handle the case where a PDL is part
* of the argument list. Ideally this would use the existing threadloop
* code but that seems too hard.
*
* kludge_copy recursively walks down the dim list of both the source and dest
* pdls, copying values in as we go. It differs from PP copy in that it operates
* on only a portion of the output pdl.
*
* (If I were Lazier I would have popped up into the perl level and used threadloops to
* assign to a slice of the output pdl -- but this is probably a little faster.)
*
* -CED 17-Jun-2004
*
* Arguments:
* poff is an integer indicating which element along the current direction is being treated (for padding accounting)
* pdata is a pointer into the destination PDL's data;
* pdims is a pointer to the destination PDL's dim list;
* ndims is the size of the destination PDL's dimlist;
* level is the conjugate dimension along which copying is happening (indexes pdims).
* "conjugate" means that it counts backward through the dimension array.
* stride is the increment in the data array corresponding to this dimension;
*
* pdl is the input PDL.
* plevel is the dim number for the input PDL, which works in the same sense as level.
* It is offset to account for the difference in dimensionality between the input and
* output PDLs. It is allowed to be negative (which is equivalent to the "permissive
* slicing" that treats missing dimensions as present and having size 1), but should
* not match or exceed pdl->ndims.
* pptr is the current offset data pointer into pdl->data.
*
* Kludge-copy works backward through the dim lists, so that padding is simpler: if undefval
* padding is required at any particular dimension level, the padding occupies a contiguous
* block of memory.
*/
PDL_Indx pdl_kludge_copy_$type(PDL_Indx poff,
PDL_$type* pdata,
PDL_Indx* pdims,
PDL_Indx ndims,
int level,
PDL_Indx stride,
pdl* pdl,
int plevel,
void* pptr,
double undefval
) {
PDL_Indx i;
PDL_Indx undef_count = 0;
#ifdef DEBUG_KLUDGE_COPY
printf("entering pdl_kludge_copy: level=%d, ndims=%d, plevel=%d; pdl->ndims=%d\\n",level,ndims,plevel,pdl->ndims);
#endif
if(level > ndims ) {
fprintf(stderr,"pdl_kludge_copy: level=%d; ndims=%d\\n",level,ndims);
croak("Internal error - please submit a bug report at http://sourceforge.net/projects/pdl/:\\n pdl_kludge_copy: Assertion failed; ndims-1-level (%d) < 0!.",ndims-1-level);
}
if(level >= ndims - 1) {
/* In this case we are in as far as we can go in the destination PDL, so direct copying is in order. */
int pdldim = pdl->ndims - 1 - plevel;
PDL_Indx pdlsiz;
int oob = (ndims-1-level < 0);
/* Do bounds checking on the source dimension -- if we wander off the end, we
* are doing permissive-slicing kind of stuff; if we wander off the beginning, we
* are doing dimensional padding. In either case, we just iterate once.
*/
if(pdldim < 0 || pdldim >= pdl->ndims) {
pdldim = (pdldim < 0) ? (0) : (pdl->ndims - 1);
pdlsiz = 1;
} else {
pdlsiz = pdl->dims[pdldim];
}
#ifdef DEBUG_KLUDGE_COPY
fprintf(stderr," pdldim_expr=%d, pdldim=%d, pdlsiz=%d, pdims[0]=%d\\n",pdl->ndims-1-plevel,pdldim, pdlsiz,pdims[0]);
#endif
switch(pdl->datatype) {
!WITH!SUBS!
# perl loop to emit code for all the PDL types
#
foreach my $switch_type (keys %PDL::Types::typehash) {
my $ctype = $PDL::Types::typehash{$switch_type}{ctype};
print OUT <<"!WITH!SUBS!";
case ${switch_type}:
#ifdef DEBUG_KLUDGE_COPY
if(pptr && pdata) {
fprintf(stderr,"*** kludge Assigning to %d (num - %g); pdlsiz=%d\\n",pdata, (double)( *((${ctype} *)pptr)),pdlsiz);
} else {
fprintf(stderr,"*** kludge Skipping assignment (null pointer in source)\\n");
}
#endif
/* copy data (unless the source pointer is null) */
i=0;
if(pptr && pdata && pdlsiz) {
for(; i<pdlsiz; i++)
pdata[i] = (PDL_$type) ((${ctype} *)pptr)[i];
} else {
if(pdata)
pdata[i] = undefval;
}
/* pad out, in the innermost dimension */
#ifdef DEBUG_KLUDGE_COPY
fprintf(stderr,"padding; "); fflush(stderr);
fprintf(stderr," ndims-1-level=%d; ndims=%d, level=%d; plevel=%d; pdl->ndims=%d; poff=%ld; oob=%d\\n", ndims-1-level,ndims, level, plevel, pdl->ndims, poff, oob);
#endif
if( !oob ) {
for(; i< pdims[0]-poff; i++) {
undef_count++;
pdata[i] = undefval;
}
#ifdef DEBUG_KLUDGE_COPY
fprintf(stderr," filled in row: ");
for(i=0;i<pdims[0] - poff; i++)
fprintf(stderr,"%g ",pdata[i]);
fprintf(stderr,"\\n");
#endif
}
break;
!WITH!SUBS!
} # end of foreach in the perl generator code
print OUT <<"!WITH!SUBS!";
default:
croak("Internal error - please submit a bug report at http://sourceforge.net/projects/pdl/:\\n pdl_kludge_copy: unknown type of %d.",pdl->datatype);
break;
}
return undef_count;
}
/* If we're here we're not at the bottom level yet... */
#ifdef DEBUG_KLUDGE_COPY
printf("Entering PKC recursion loop. imax is %d (plevel=%d; pdl->dims[%d-1-%d] = %d)\\n",
( (plevel >= 0 && (pdl->ndims - 1 - plevel >= 0) && (pdl->ndims - 1 - plevel < pdl->ndims)) ? (pdl->dims[ pdl->ndims - 1 - plevel ]) : 1 ),
plevel,
ndims,
plevel,
pdl->dims[pdl->ndims-1-plevel]
);
#endif
for(i=0;
i < ( (plevel >= 0 && (pdl->ndims - 1 - plevel >= 0) && (pdl->ndims - 1 - plevel < pdl->ndims)) ? (pdl->dims[ pdl->ndims-1-plevel ]) : 1 );
i++) {
#ifdef DEBUG_KLUDGE_COPY
{
char buf[10240];
char sb[1024];
int ii;
*buf=0;
for(ii=0;ii<ndims;ii++) {
sprintf(sb,"%d",pdims[ii]);
if(ii) strcat(buf,",");
strcat(buf,sb);
}
printf("pdl_kludge_copy: pushing... level=%d, i=%d, pdata=%d, pdims=(%s), ndims=%d, stride=%d, plevel=%d, pptr=%d\\n",
level, i, pdata, buf, ndims, stride, plevel, pptr);
printf(" (pdims[ndims-2-level] is %d)\\n",pdims[ndims-2-level]);
}
#endif
undef_count += pdl_kludge_copy_$type(0, pdata + stride * i,
pdims,
ndims,
level+1,
stride / ((pdims[ndims-2-level]) ? (pdims[ndims-2-level]) : 1),
pdl,
plevel+1,
((PDL_Byte *) pptr) + pdl->dimincs[pdl->ndims-1-plevel] * i * pdl_howbig(pdl->datatype),
undefval
);
#ifdef DEBUG_KLUDGE_COPY
printf("(level=%d, back from recursion)\\n",level);
#endif
} /* end of kludge_copy recursion loop */
/* pad tree to zero if there are not enough elements... */
#ifdef DEBUG_KLUDGE_COPY
printf(" pdl_kludge_copy - finished recursion loop at level %d... i is %d, pdims[%d]=%d...\\n",
level, i, ndims - 1 - level,pdims[ndims - 1 - level]);
#endif
if(i < pdims[ndims - 1 - level]) {
int cursor, target;
cursor = i * stride;
target = pdims[ndims-1-level]*stride;
undef_count += target - cursor;
#ifdef DEBUG_KLUDGE_COPY
printf("i=%d, pdims[%d - 1 - %d]=%d, stride=%d, cursor=%d, target=%d....\\n",i,ndims,level,pdims[ndims - 1 - level],stride,cursor,target);
{
int ii;
printf(" pdims is:[");
for(ii=0; ii<ndims;ii++) {
printf("%s%d",(ii?", ":""),pdims[ii]);
}
printf("]\\n");
}
#endif
for(;
cursor < target;
cursor++) {
pdata[cursor] = undefval;
}
} /* end of padding IF statement */
return undef_count;
}
/*
* pdl_setav_type loads a new PDL with values from a Perl AV, another PDL, or
* a mix of both. Heterogeneous sizes are handled by padding the new PDLs
* values out to size with the undefval. It is called by pdl_setav only.
*
*
* Recent changes to setav_$type:
* - Look for PDLs and deep copy them with pdl_kludge_copy just as if they were
* array refs.
* - Allow multiple depths in different elements. The max depth should have been
* determined by pdl_av_ndcheck (it comes in as the ndims parameter), so other
* depths are descended to and extra elements are filled in with zeroes.
* In the Best of All Possible Worlds this would be badval compliant and
* the extra elements would get filled in with BAD.
* --CED, 17-Jun-2004
*
* - Check for undef values and set the PDL element to PDL::undefval for those
* elements.
* - Keep track of how many undef values were encountered, for debugging
* (I am wary of action-at-a-distance)
*
* -CED, 28-Jul-2004
*
*
* - pdata is the data pointer from a PDL
* - av is the array ref (or PDL) to use to fill the data with,
* - pdims is the dimlist
* - ndims is the size of the dimlist
* - level is the recursion level, which is also the dimension that we are filling
*/
PDL_Indx pdl_setav_$type(PDL_$type* pdata, AV* av,
PDL_Indx* pdims, PDL_Long ndims, int level, double undefval)
{
PDL_Indx cursz = pdims[ndims-1-level]; /* we are starting from the highest dim
inwards */
PDL_Indx len = av_len(av);
PDL_Indx i,stride=1;
SV *el, **elp;
PDL_Indx undef_count = 0;
fflush(stdout);
for (i=0;i<ndims-1-level;i++) {
stride *= pdims[i];
}
#ifdef DEBUG_SETAV_TYPE
printf("entering pdl_setav_$type: pdata=%d\\n",pdata);
{
int i;
printf("ndims=%d; level=%d; ndims-1-level is %d; pdims factor is (",ndims,level,ndims-1-level);
for(i=0;i<ndims-1-level; i++) {
printf("%s%d",(i?", ":""),pdims[i]);
}
printf("); initial stride is %d\\n",stride);
}
#endif
for (i=0;i<=len;i++,pdata += stride) { /* note len is actually the highest index in the array */
int foo;
#ifdef DEBUG_SETAV_TYPE
{
/* Generate some debugging info */
char buf[10240];
char sb[1024];
int ii;
*buf=0;
for(ii=0;ii<ndims; ii++) {
sprintf(sb,"%d",pdims[ii]);
if(ii) strcat(buf,",");
strcat(buf,sb);
}
printf("pdl_setav_$type: level=%d, i=%d, pdata=%d, pdims=(%s), ndims=%d\\t",
level, i, pdata, buf, ndims);
}
#endif
elp = av_fetch(av,i,0);
el = (elp ? *elp : 0);
foo = el ? SVavref(el) : 0;
if (foo) {
#ifdef DEBUG_SETAV_TYPE
printf("found an array ref -- recursing\\n");
#endif
undef_count += pdl_setav_$type(pdata, (AV *) SvRV(el), pdims, ndims, level+1, undefval);
} else if( el && SvROK(el) ) {
pdl *pdl;
if( (pdl = SvPDLV(el)) ) {
pdl_make_physical(pdl);
#ifdef DEBUG_SETAV_TYPE
printf("source pdl has %d vals...\\n",pdl->nvals);
#endif
{ // convenience block.. recursively copy/pad each PDL.
#ifdef DEBUG_SETAV_TYPE
{
/* Generate some debugging info */
char buf[10240];
char sb[1024];
int ii;
*buf=0;
for(ii=0;ii<ndims; ii++) {
sprintf(sb,"%d",pdims[ii]);
if(ii) strcat(buf,",");
strcat(buf,sb);
}
printf("Calling pdl_kludge_copy - pdata is %d, pdims is (%s), ndims is %d, level is %d, stride is %d, pdl->data is %d\\n",
pdata, buf, ndims, level, stride, pdl->data
);
}
#endif
{
PDL_Indx pd;
int pddex;
pddex = ndims - 2 - level;
pd = (pddex >= 0 && pddex < ndims ? pdims[ pddex ] : 0);
if(!pd)
pd = 1;
undef_count += pdl_kludge_copy_$type(0, pdata,pdims,ndims, level+1, stride / pd , pdl, 0, pdl->data,undefval);
}
}
} else {
croak("Non-array, non-PDL element in list");
}
} else { /* SvROK(el)==0; scalar element */
if( sv_undef(el) ) {
#ifdef DEBUG_SETAV_TYPE
printf("undefined scalar element\\n");
#endif
*pdata = (PDL_$type) undefval; /* undef case */
undef_count++;
} else {
#ifdef DEBUG_SETAV_TYPE
printf("defined scalar element %g\\n",(double)SvNV(el));
#endif
*pdata = (PDL_$type) SvNV(el);
}
/* Pad dim if we are not deep enough */
if(level < ndims-1) {
PDL_$type *cursor = pdata;
PDL_$type *target = pdata + stride;
#ifdef DEBUG_SETAV_TYPE
printf("\\tPadding: level=%d; ndims-1=%d. pdata is %d, stride is %d, target is %d\\n",level,ndims-1,pdata,stride, target);
#endif
for( cursor++; cursor < target; cursor++ ) {
*cursor = (PDL_$type)undefval;
undef_count++;
}
}
}
} /* end of i loop */
/* in case this dim is incomplete set remaining elements to the undefval */
#ifdef DEBUG_SETAV_TYPE
printf("\\tloop is complete. len is %d, cursz-1 is %d, stride is %d\\n",len,cursz-1, stride);
#endif
if(len < cursz-1 ) {
PDL_$type *target = pdata + stride * (cursz - 1 - len);
#ifdef DEBUG_SETAV_TYPE
printf("\\tpadding %d elements with the undefval...\\n",target-pdata);
#endif
for( ;
pdata < target;
pdata++
) {
*pdata = (PDL_$type) undefval;
undef_count++;
}
}
if(level==0 && undef_count) {
char debug_flag;
SV *sv;
sv = get_sv("PDL::debug",0);
debug_flag = (sv_undef(sv)) ? 0 : (char)SvIV(sv);
if(debug_flag) {
fprintf(stderr,"Warning: pdl_setav_$type converted undef to $PDL::undefval (%g) %ld time%s\\n",undefval,undef_count,undef_count==1?"":"s");
}
}
return undef_count;
}
!WITH!SUBS!
} # end type loop