Vector multiplication with BLAS' catlas_saxpby not working properly

I am trying to have two arbitrary length vectors (typical length will be 2048) and multiply element by element. So Z[n] = X[n] * Y[n] for all n.

The code I have setup to test is fairly basic:

float inputX[4] = { 2, 4, 8, 16 };
float inputY[4] = { 2, 4, 8, 16 };

catlas_saxpby(4, 1, inputX, 1, 1, inputY, 1);

The result goes into inputY, and the result is

4.000000, 8.000000, 16.000000, 32.000000

Which if they were multiplying it should be 4, 16, 64, 256. But it looks like it is adding.

So this is not doing what I expect, and the documentation doesn't give me enough information to figure what it is doing.

Any ideas?

Apple's documentation for BLAS says th开发者_高级运维is: 

Computes the product of two vectors, scaling each one separately (single-precision).

void catlas_saxpby (
   const int N,
   const float alpha,
   const float *X,
   const int incX,
   const float beta,
   float *Y,
   const int incY
);
Parameters
N
Number of elements in the vector.
alpha
Scaling factor for X.
X
Input vector X.
incX
Stride within X. For example, if incX is 7, every 7th element is used.
beta
Scaling factor for Y.
Y
Input vector Y.
incY
Stride within Y. For example, if incY is 7, every 7th element is used.
Discussion
On return, the contents of vector Y are replaced with the result.

As Adam Rosenfield said, the documentation is incorrect. Please file a bug.

That aside, there are some corrections to the rest of his answer. First, saxpby computes alpha * X + beta * Y. Second, and more useful for you: there is no function in BLAS that does what you want, but there is exactly such a function in vDSP, which is also part of the Accelerate.framework: vDSP_vmul.

The Apple documentation is mistaken. The saxpby function computes the expression alpha*X + beta*Y for scalars alpha and beta and vectors X and Y.

I don't think there's a function available to compute the element-wise product of two vectors, since that's not a common operation in linear algebra. You could take the diagonal of the outer product, but that's a severe waste of effort since it computes the entire outer product (N² multiplications instead of N).