Some BLAS Level 1 Functions

For the sake of simplicity we will only consider matrix/vector elements of type double in this class.
All vectors must be allocated before the function gets called.
An operator description like \(x \rightarrow y\) is read as vector \(y\) gets overwritten with vector \(y\).
Furthermore:
- lower case letters indicate vectors, upper case letters matrices.
- it is assumed that different letters refer to non-overlapping vectors/matrices, i.e. they have no elements in common. Otherwise the result might be undefined.

All BLAS function names get an additional prefix that depends on the matrix/vector element type:

Prefix	Meaning	Element Type in C
s	single precision	float
d	double precision	double
c	complex single precision	float _Complex
z	complex double precision	double _Complex

So functions scopy, dcopy, ccopy and zcopy provide the same functionality for vectors where elements are of type float, double, float _Complex and double _Complex respectively.

Short Description of some BLAS Level 1 Functions

`copy`: Copy Vectors

Operation	Signature
\(x \rightarrow y\)	void dcopy(size_t n, const double x, ptrdiff_t incX, double y, ptrdiff_t incY);

`swap`: Interchanging Vectors

Operation	Signature
\(x \leftrightarrow y\)	void dswap(size_t n, double x, ptrdiff_t incX, double y, ptrdiff_t incY);

`axpy`: Adding Vectors (Alpha x Plus y)

Operation	Signature
\(\alpha x+y\rightarrow y\)	void daxpy(size_t n, double alpha, const double x, ptrdiff_t incX, double y, ptrdiff_t incY);

In the special case \(\alpha=0\) the function returns immediately (NOP case).

`scal`: Scaling Vectors

Operation	Signature
\(\alpha x\rightarrow x\)	void dscal(size_t n, double alpha, double *x, ptrdiff_t incX);

In the special case \(\alpha=1\) the function returns immediately (NOP case).
In the special case \(\alpha=0\) must not multiply elements of \(x\) with zero. Instead \(x\) gets overwritten with zeros. Hence, in this case we allow that \(x\) contains NaN values.

`dot`: dot product

Operation	Signature
\(x^T y \rightarrow alpha\)	double ddot(size_t n, const double x, ptrdiff_t incX, const double y, ptrdiff_t incY);

Exercise

Complete the following code by

implementing the BLAS Level 1 functions and
doing some (very) simple test calls in function main

session04/blas1_ex.c

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

void
initMatrix(size_t m, size_t n,
     double *A,
     ptrdiff_t incRow, ptrdiff_t incCol)
{
    for (size_t i=0; i<m; ++i) {
        for (size_t j=0; j<n; ++j) {
            A[i*incRow+j*incCol] = i*n+j+1;
        }
    }
}

void
printMatrix(size_t m, size_t n,
            const double *A,
            ptrdiff_t incRow, ptrdiff_t incCol)
{
    for (size_t i=0; i<m; ++i) {
        for (size_t j=0; j<n; ++j) {
            printf("%10.3lf", A[i*incRow+j*incCol]);
        }
        printf("\n");
    }
    printf("\n");
}

//------------------------------------------------------------------------------

void
dcopy(size_t n,
      const double *x, ptrdiff_t incX,
      double *y, ptrdiff_t incY)
{
    // *** Your code here ***
}

void
dswap(size_t n,
      double *x, ptrdiff_t incX,
      double *y, ptrdiff_t incY)
{
    // *** Your code here ***
}

void
daxpy(size_t n, double alpha,
      const double *x, ptrdiff_t incX,
      double *y, ptrdiff_t incY)
{
    // *** Your code here ***
}

void
dscal(size_t n, double alpha,
      double *x, ptrdiff_t incX)
{
    // *** Your code here ***
}

double
ddot(size_t n,
     const double *x, ptrdiff_t incX,
     const double *y, ptrdiff_t incY)
{
    // *** Your code here ***
}

int
main()
{
    // Allocate matrix A with dimension m x n
    size_t m = 5;
    size_t n = 8;

    size_t incRow = n;
    size_t incCol = 1;

    double *A = (double *) malloc(m*n*sizeof(double));

    // Initialize matrix A
    initMatrix(m, n, A, incRow, incCol);

    // Print initial matrix A
    printf("A=\n");
    printMatrix(m, n, A, incRow, incCol);

    // Some simple tests for the BLAS Level 1 functions

    // 1) dcopy
    printf("A(:,2) -> A(:, 3)\n");
    // *** Your code here ***
    printMatrix(m, n, A, incRow, incCol);

    // 2) dswap
    printf("A(3,:) <-> A(1,:)\n");
    // *** Your code here ***
    printMatrix(m, n, A, incRow, incCol);

    // 3) daxpy
    printf("2.5*A(0,:) + A(2,:) -> A(2,:)\n");
    // *** Your code here ***
    printMatrix(m, n, A, incRow, incCol);

    // 4) dscal
    printf("0.5*A(:,5) -> A(:,5)\n");
    // *** Your code here ***
    printMatrix(m, n, A, incRow, incCol);

    // Deallocate matrix A
    free(A);
}

Content

Some BLAS Level 1 Functions

Short Description of some BLAS Level 1 Functions

copy: Copy Vectors

swap: Interchanging Vectors

axpy: Adding Vectors (Alpha x Plus y)

scal: Scaling Vectors