[Numpy-discussion] How can CDOUBLE_to_CDOUBLE work correctly?

[Pearu Peterson]

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 */