[Numpy-discussion] cross

Revaz Yves yves.revaz@obspm...
Mon Mar 3 13:20:54 CST 2008

```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;
}

```