[Numpy-discussion] How can CDOUBLE_to_CDOUBLE work correctly?
Pearu Peterson
pearu at cens.ioc.ee
Sun Mar 3 04:41:15 CST 2002
Hi!
I am trying to copy a 2-d complex array to another 2-d complex array in an
extension module. Both arrays may be noncontiguous. But using
the same routine (namely Travis's copy_ND_array, you can find it at the
end of this messsage) as for real arrays seems not work.
After some playing around and reading docs about strides and stuff,
I found that the reason might be in how Numeric (20.3) defines
the CDOUBLE_to_CDOUBLE function:
static void CDOUBLE_to_CDOUBLE(double *ip, int ipstep,
double *op, int opstep, int n)
{int i; for(i=0;i<2*n;i++,ip+=ipstep,op+=opstep) {*op = (double)*ip;}}
It seems not to take into account that real and imaginary part are
always contiguous in memory, even if an array itself is not. Actually, I
don't understand how this code can work (unless some magic is done in
places where this code is used). I would have expected that the code
for CDOUBLE_to_CDOUBLE to be analoguous to relative functions but for the
real data. For example, DOUBLE_to_DOUBLE is defined as
static void DOUBLE_to_DOUBLE(double *ip, int ipstep,
double *op, int opstep, int n)
{int i; for(i=0;i<n;i++,ip+=ipstep,op+=opstep) {*op = (double)*ip;}}
So, I would have expected CDOUBLE_to_CDOUBLE to be
static void CDOUBLE_to_CDOUBLE(double *ip, int ipstep,
double *op, int opstep, int n)
{ int i;
for(i=0;i<n;i++,ip+=ipstep,op+=opstep) {
*op = (double)*ip; /* copy real part */
*(op+1) = (double)*(ip+1); /* copy imaginary part that always
follows the real part in memory */
}
}
Could someone explain how Numeric can work with the current
CDOUBLE_to_CDOUBLE?
Because I don't understand which one is broken, the Numeric's
CDOUBLE_to_CDOUBLE (and relative) functions, or my code.
The latter may be the case, but to fix it, I need some clarification on
the issue. Can you help me?
Thanks,
Pearu
/************************* copy_ND_array *******************************/
#define INCREMENT(ret_ind, nd, max_ind) \
{ \
int k; \
k = (nd) - 1; \
if (k<0) (ret_ind)[0] = (max_ind)[0]; else \
if (++(ret_ind)[k] >= (max_ind)[k]) { \
while (k >= 0 && ((ret_ind)[k] >= (max_ind)[k]-1)) \
(ret_ind)[k--] = 0; \
if (k >= 0) (ret_ind)[k]++; \
else (ret_ind)[0] = (max_ind)[0]; \
} \
}
#define CALCINDEX(indx, nd_index, strides, ndim) \
{ \
int i; \
indx = 0; \
for (i=0; i < (ndim); i++) \
indx += nd_index[i]*strides[i]; \
}
extern
int copy_ND_array(const PyArrayObject *in, PyArrayObject *out)
{
/* This routine copies an N-D array in to an N-D array out where both
can be discontiguous. An appropriate (raw) cast is made on the data.
*/
/* It works by using an N-1 length vector to hold the N-1 first indices
into the array. This counter is looped through copying (and casting)
the entire last dimension at a time.
*/
int *nd_index, indx1;
int indx2, last_dim;
int instep, outstep;
if (0 == in->nd) {
in->descr->cast[out->descr->type_num]((void *)in->data,1,
(void*)out->data,1,1);
return 0;
}
if (1 == in->nd) {
in->descr->cast[out->descr->type_num]((void *)in->data,1,
(void*)out->data,1,in->dimensions[0]);
return 0;
}
nd_index = (int *)calloc(in->nd-1,sizeof(int));
last_dim = in->nd - 1;
instep = in->strides[last_dim] / in->descr->elsize;
outstep = out->strides[last_dim] / out->descr->elsize;
if (NULL == nd_index ) {
fprintf(stderr,"Could not allocate memory for index array.\n");
return -1;
}
while(nd_index[0] != in->dimensions[0]) {
CALCINDEX(indx1,nd_index,in->strides,in->nd-1);
CALCINDEX(indx2,nd_index,out->strides,out->nd-1);
/* Copy (with an appropriate cast) the last dimension of the array */
(in->descr->cast[out->descr->type_num])((void*)(in->data+indx1),instep,
(void*)(out->data+indx2),outstep,in->dimensions[last_dim]);
INCREMENT(nd_index,in->nd-1,in->dimensions);
}
free(nd_index);
return 0;
}
/* EOF copy_ND_array */
More information about the Numpy-discussion
mailing list