DAXPY - A BLAS Example

Linear Algebra

A BLAS Example

axpy

In this section we will re-implement axpy as it is in BLAS.

Add $αx$ to $y$ .

$N$ is the number of elements in the input vectors
$D A$ is the constant alpha
$D X$ this is the X array
$I NCX$ spacing between elements X
$D Y$ this is the Y array
$I NC Y$ this is the spacing between elements

Remark: $I NCX$ and $I NC Y$ let you stride through the array. i.e. you could multiply every other element or every fifth element. This will cause a lot of cache misses, which is why this is not advisable to use. $I NCX$ and $I NC Y$ can be negative... this lets you perform the operation backwards.

Here we show the code for the a modified C version of the BLAS routine DAXPY, which is the double version of the axpy routine. In the code below we allocate memory using malloc and then we subsequently free the memory. The memory is sized to be $N \times sizeof(double)$ , this means that there is enough memory to store $N$ doubles.

#include <stdio.h>
#include <stdlib.h>

/* BLAS-style DAXPY: y <- alpha*x + y */
void daxpy(int N, double alpha, const double *x, int incx, double *y, int incy)
{
    if (N <= 0 || alpha == 0.0)
        return;

    int ix = 0;
    int iy = 0;

    if (incx < 0) ix = (1 - N) * incx;
    if (incy < 0) iy = (1 - N) * incy;

    for (int i = 0; i < N; ++i) {
        y[iy] += alpha * x[ix];
        ix += incx;
        iy += incy;
    }
}

int main(void)
{
    int N = 8;
    double alpha = 2.0;

    /* allocate vectors */
    double *x = malloc(N * sizeof(double));
    double *y = malloc(N * sizeof(double));
    if (!x || !y) {
        perror("malloc failed");
        free(x);
        free(y);
        return 1;
    }

    /* initialize vectors */
    for (int i = 0; i < N; ++i) {
        x[i] = (double)(i + 1); // 1,2,3,...
        y[i] = 1.0;
    }

    /* perform daxpy: y = alpha*x + y */
    daxpy(N, alpha, x, 1, y, 1);

    /* print results */
    for (int i = 0; i < N; ++i)
        printf("y[%d] = %f\n", i, y[i]);

    free(x);
    free(y);
    return 0;
}

The mallocs and frees can be quite fun to work with if you are not careful. Add another free to the end of the file and you will begin to see what I mean.