Some BLAS Level 1 functions

Content

BLAS Level 1 functions operate on vectors. In practical cases these vectors are often rows or columns of a matrix. Hence, it is important to recall that in memory the vector elements are separated from each other by a constant increment that can be larger than one.

Skeleton for the exercises

You can use the following skeleton to solve the exercises below step by step:

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

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

void
printMatrix(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("%8.3lf ", A[i*incRowA+j*incColA]);
        }
        printf("\n");
    }
    printf("\n");
}


// function dnrm1
// TODO: implement drnm1 here


// function dswap
// TODO: implement dswap here


// function daxpy
// TODO: implement daxpy here


int
main()
{
    size_t m = 7;
    size_t n = 8;

    ptrdiff_t incRowA = 1;
    ptrdiff_t incColA = m;

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

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

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

    // Exercise 3.1:
    // -------------
    /*
    //TODO/FIXME: printf("nrm1(A(1,:)) = %lf\n", dnrm1(???));
    //TODO/FIXME: printf("nrm1(A(:,2)) = %lf\n", dnrm1(???));

    size_t mn = (m<n) ? m : n;  // length of diagonal
    //TODO/FIXME: printf("nrm1(diag(A)) = %lf\n", dnrm1(???));
    */

    // Exercise 3.2:
    // -------------
    /*
    // TODO/FIXME: dswap(???);
    printf("After dswap(A(1,:), A(2,:))\n");
    printMatrix(m, n, A, incRowA, incColA);

    // TODO/FIXME: dswap(???);
    printf("After dswap(A(:,1), A(:,2))\n");
    printMatrix(m, n, A, incRowA, incColA);
    */

    // Exercise 3.3:
    // -------------
    /*
    // TODO/FIXME: daxpy(???);
    printf("After daxpy(3.0, A(:,1), A(:,2))\n");
    printMatrix(m, n, A, incRowA, incColA);

    // TODO/FIXME: daxpy(???);
    printf("After daxpy(-2.5, A(1,:), A(3,:))\n");
    printMatrix(m, n, A, incRowA, incColA);
    */

    free(A);

    return 0;
}

theon$ gcc -Wall -std=c99 -o blas1_example blas1_example.c
theon$ ./blas1_example
A = 
   1.000    2.000    3.000    4.000    5.000    6.000    7.000    8.000 
   9.000   10.000   11.000   12.000   13.000   14.000   15.000   16.000 
  17.000   18.000   19.000   20.000   21.000   22.000   23.000   24.000 
  25.000   26.000   27.000   28.000   29.000   30.000   31.000   32.000 
  33.000   34.000   35.000   36.000   37.000   38.000   39.000   40.000 
  41.000   42.000   43.000   44.000   45.000   46.000   47.000   48.000 
  49.000   50.000   51.000   52.000   53.000   54.000   55.000   56.000 
theon$ 

Exercise: One-Norm of a vector dnrm1

Write a function dnrm1 that computes and returns

\[\alpha = \| x \|_1 = \left\| \left(\begin{array}{c} x_1 \\ \vdots \\ x_n \end{array}\right) \right\|_1 = |x_1| + \dots + |x_n|\]

Adjust the previous test program such that it computes the one-norm

Exercise: Swap vectors dswap

Write a function dswap that swaps (i.e. exchanges) the elements of two vectors.

Extend the test program, such that it

Exercise: Adding a scaled vector daxpy

The so called axpy-operation (alpha x plus y) computes

\[y \leftarrow \alpha x + y\]

Implement a function daxpy for this operation. Test it by