Some BLAS Level 1 Functions

Content

Some BLAS Level 1 Functions
- Exercise: BLAS Level 1
Possible Solution
Simple Benchmark for Matrix Initialization
Possible Solution

We will implement the following BLAS Level 1 operations:

Function dcopy

Copying vectors denoted mathematically as $x \to y$.
Function daxpy

Scaled vector update denoted mathematically as $\alpha x + y \to y$. Hence, axpy for alpha x plus y.
Function ddot

Computing the dot product of two vector, i.e. $x^T \cdot y$.
Function dscal
- If $\alpha \neq 0$ it computed the scaled vector $x \leftarrow \alpha x$.
- If $\alpha = 0$ it overwrites $x$ with zeros, i.e. $x \leftarrow \mathcal{0}$.
Note that these cases need to be handled separately because in the case $\alpha=0$ vector $x$ is allowed to contain Not-a-Number (NaN) values.

If a matrix has NaN entries and $\alpha \neq 0$ the result of scal is allowed to be undefined.

Exercise: BLAS Level 1

The signature of function dcopy is given:
- Add the implementation
- Fix the function call to match the operation indicated by the preceding printf statement
How function daxpy is supposed to be used is given in main:
- Derive a proper declaration for the function and
- implement the function.
For function ddot we give fewer hints. Specify a signature that indicates:
- The function expects two vectors of same length.
- The vectors are not modified by the function.
- The function returns the result.
Then implement the function and use it in main.
For function dscal we provide a skeleton:
- In case $\alpha=1$ the function makes a quick return.
- Implement the cases for
  - $\alpha \neq 1$ and \alpha \neq 0$
  - $\alpha \neq 1$ and \alpha = 0$
Note that main contains a test where NaNs are inserted through function nan (declared in math.h).

#include <stdlib.h>         // for malloc(), free()
#include <stdio.h>          // for printf()
#include <stddef.h>         // for size_t, ptrdiff_t
#include <math.h>           // for nan()

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

void
printGeMatrix(size_t m, size_t n,
              const double *A,
              ptrdiff_t incRowA, ptrdiff_t incColA)
{
    for (size_t i=0; i<m; ++i) {
        for (size_t j=0; j<n; ++j) {
            printf("%9.2lf ", A[i*incRowA + j*incColA]);
        }
        printf("\n");
    }
    printf("\n");
}

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

// function daxpy
/*
   Your code here
*/

// function ddot
/*
   Your code here
*/

#ifndef COLMAJOR
#define COLMAJOR 1
#endif

int
main()
{
    printf("COLMAJOR = %d\n", COLMAJOR);

    size_t m = 5, n = 10;

    ptrdiff_t incRowA = COLMAJOR ? 1 : n;
    ptrdiff_t incColA = COLMAJOR ? m : 1;

    double *A = malloc(m*n*sizeof(double));
    if (!A) {
        abort();
    }

    initGeMatrix(m, n, A, incRowA, incColA);

    printf("A =\n");
    printGeMatrix(m, n, A, incRowA, incColA);

    printf("ddot(A(:,1), A(:,4)) = %lf\n",
           ddot(...));                                              // <- FIXME

    printf("compute: A(:,1) -= 0.5*A(:,4)\n");
    daxpy(m, -0.5, &A[4*incColA], incRowA, &A[1*incColA], incRowA);

    printf("A =\n");
    printGeMatrix(m, n, A, incRowA, incColA);

    printf("copy: A(1,:) = A(2,:)\n");
    dcopy(m, A, incRowA, A, incColA);                               // <- FIXME

    printf("A =\n");
    printGeMatrix(m, n, A, incRowA, incColA);

    printf("scal: A(1,:) = 1.5*A(1,:)\n");
    dscal(n, 1.5, &A[1*incRowA], incColA);

    printf("A =\n");
    printGeMatrix(m, n, A, incRowA, incColA);

    //
    // We manually insert some NaNs into A
    //
    A[1*incRowA+0*incColA] = nan("");
    A[1*incRowA+3*incColA] = nan("");
    A[1*incRowA+9*incColA] = nan("");

    printf("\n\nAfter inserting some NaNs into A:\n");
    printf("A =\n");
    printGeMatrix(m, n, A, incRowA, incColA);

    printf("scal: A(1,:) = 1.5*A(1,:)       <- will give wrong results\n");
    dscal(n, 1.5, &A[1*incRowA], incColA);

    printf("A =\n");
    printGeMatrix(m, n, A, incRowA, incColA);

    printf("scal: A(1,:) = 0*A(1,:)\n");
    dscal(n, 0.0, &A[1*incRowA], incColA);

    printf("A =\n");
    printGeMatrix(m, n, A, incRowA, incColA);

    free(A);
}