[Numpy-discussion] cross

[Revaz Yves]

Dear List,

I'm computing the cross product of positions and velocities of n points 
in a 3d space.
Using the numpy function "cross", this can be written as :

n=1000
pos = random.random([n,3])
vel = random.random([n,3])
cross(pos,vel)

I compare the computation time needed with a C-api I wrote (dedicated to 
this operation).
It appears that my api is in average 20 times faster than the cross 
function (for n between 100 and 1000000),
making the latter useless for my purpose :-( .

Is it normal ? or I'm I using the "cross" function the wrong way ?

yves




PS :Here after you can see some lines the of the C-api.



          if (!PyArg_ParseTuple(args, "OO", &pos , &vel))       
              return NULL;
       
          /* create a NumPy object similar to the input */
          int   ld[2];
          ld[0]=pos->dimensions[0];
          ld[1]=pos->dimensions[1];
          lxyz = (PyArrayObject *) 
PyArray_FromDims(pos->nd,ld,pos->descr->type_num); 
      
               
            /* loops over all elements */
            for (i = 0; i < pos->dimensions[0]; i++) {

              x  = (float *) (pos->data + i*(pos->strides[0])            
   );
              y  = (float *) (pos->data + i*(pos->strides[0]) + 
1*pos->strides[1]);
              z  = (float *) (pos->data + i*(pos->strides[0]) + 
2*pos->strides[1]);
       
              vx = (float *) (vel->data + 
i*(vel->strides[0])                    );
              vy = (float *) (vel->data + i*(vel->strides[0]) + 
1*vel->strides[1]);
              vz = (float *) (vel->data + i*(vel->strides[0]) + 
2*vel->strides[1]);   
               
              lx = (*y * *vz - *z * *vy);
              ly = (*z * *vx - *x * *vz);
              lz = (*x * *vy - *y * *vx);
       
              *(float *)(lxyz->data + i*(lxyz->strides[0]) + 
0*lxyz->strides[1])  = lx;
              *(float *)(lxyz->data + i*(lxyz->strides[0]) + 
1*lxyz->strides[1])  = ly;
              *(float *)(lxyz->data + i*(lxyz->strides[0]) + 
2*lxyz->strides[1])  = lz;
            }