Sample solutions

Content

Solution for Exercise 3.1

#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
double
dnrm1(size_t n, const double *x, ptrdiff_t incX)
{
    double alpha = 0;
    for (size_t i=0; i<n; ++i) {
        alpha += x[i*incX];
    }
    return alpha;
}


// 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:
    // -------------
    printf("nrm1(A(1,:)) = %lf\n", dnrm1(n, &A[1*incRowA], incColA));
    printf("nrm1(A(:,2)) = %lf\n", dnrm1(m, &A[2*incColA], incRowA));

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

    // 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_sol1 blas1_example_sol1.c
theon$ ./blas1_example_sol1
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 

nrm1(A(1,:)) = 100.000000
nrm1(A(:,2)) = 189.000000
nrm1(diag(A)) = 196.000000
theon$

Solution for Exercise 3.2

session03/blas1_example_sol2.c

#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
double
dnrm1(size_t n, const double *x, ptrdiff_t incX)
{
    double alpha = 0;
    for (size_t i=0; i<n; ++i) {
        alpha += x[i*incX];
    }
    return alpha;
}


// function dswap
void
dswap(size_t n, double *x, ptrdiff_t incX, double *y, ptrdiff_t incY)
{
    for (size_t i=0; i<n; ++i) {
        double t = x[i*incX];
        x[i*incX] = y[i*incY];
        y[i*incY] = t;
    }
}



// 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:
    // -------------
    printf("nrm1(A(1,:)) = %lf\n", dnrm1(n, &A[1*incRowA], incColA));
    printf("nrm1(A(:,2)) = %lf\n", dnrm1(m, &A[2*incColA], incRowA));

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

    // Exercise 3.2:
    // -------------
    dswap(n, &A[1*incRowA], incColA, &A[2*incRowA], incColA);
    printf("After dswap(A(1,:), A(2,:))\n");
    printMatrix(m, n, A, incRowA, incColA);

    dswap(m, &A[1*incColA], incRowA, &A[2*incColA], incRowA);
    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_sol2 blas1_example_sol2.c
theon$ ./blas1_example_sol2
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 

nrm1(A(1,:)) = 100.000000
nrm1(A(:,2)) = 189.000000
nrm1(diag(A)) = 196.000000
After dswap(A(1,:), A(2,:))
   1.000    2.000    3.000    4.000    5.000    6.000    7.000    8.000 
  17.000   18.000   19.000   20.000   21.000   22.000   23.000   24.000 
   9.000   10.000   11.000   12.000   13.000   14.000   15.000   16.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 

After dswap(A(:,1), A(:,2))
   1.000    3.000    2.000    4.000    5.000    6.000    7.000    8.000 
  17.000   19.000   18.000   20.000   21.000   22.000   23.000   24.000 
   9.000   11.000   10.000   12.000   13.000   14.000   15.000   16.000 
  25.000   27.000   26.000   28.000   29.000   30.000   31.000   32.000 
  33.000   35.000   34.000   36.000   37.000   38.000   39.000   40.000 
  41.000   43.000   42.000   44.000   45.000   46.000   47.000   48.000 
  49.000   51.000   50.000   52.000   53.000   54.000   55.000   56.000 
theon$

Solution for Exercise 3.3

session03/blas1_example_sol3.c

#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
double
dnrm1(size_t n, const double *x, ptrdiff_t incX)
{
    double alpha = 0;
    for (size_t i=0; i<n; ++i) {
        alpha += x[i*incX];
    }
    return alpha;
}


// function dswap
void
dswap(size_t n, double *x, ptrdiff_t incX, double *y, ptrdiff_t incY)
{
    for (size_t i=0; i<n; ++i) {
        double t = x[i*incX];
        x[i*incX] = y[i*incY];
        y[i*incY] = t;
    }
}


// function daxpy
void
daxpy(size_t n, double alpha, const double *x, ptrdiff_t incX,
      double *y, ptrdiff_t incY)
{
    for (size_t i=0; i<n; ++i) {
        y[i*incY] += alpha*x[i*incX];
    }
}


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:
    // -------------
    printf("nrm1(A(1,:)) = %lf\n", dnrm1(n, &A[1*incRowA], incColA));
    printf("nrm1(A(:,2)) = %lf\n", dnrm1(m, &A[2*incColA], incRowA));

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

    // Exercise 3.2:
    // -------------
    dswap(n, &A[1*incRowA], incColA, &A[2*incRowA], incColA);
    printf("After dswap(A(1,:), A(2,:))\n");
    printMatrix(m, n, A, incRowA, incColA);

    dswap(m, &A[1*incColA], incRowA, &A[2*incColA], incRowA);
    printf("After dswap(A(:,1), A(:,2))\n");
    printMatrix(m, n, A, incRowA, incColA);

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

    daxpy(n, -2.5, &A[1*incRowA], incColA, &A[3*incRowA], incColA);
    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_sol3 blas1_example_sol3.c
theon$ ./blas1_example_sol3
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 

nrm1(A(1,:)) = 100.000000
nrm1(A(:,2)) = 189.000000
nrm1(diag(A)) = 196.000000
After dswap(A(1,:), A(2,:))
   1.000    2.000    3.000    4.000    5.000    6.000    7.000    8.000 
  17.000   18.000   19.000   20.000   21.000   22.000   23.000   24.000 
   9.000   10.000   11.000   12.000   13.000   14.000   15.000   16.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 

After dswap(A(:,1), A(:,2))
   1.000    3.000    2.000    4.000    5.000    6.000    7.000    8.000 
  17.000   19.000   18.000   20.000   21.000   22.000   23.000   24.000 
   9.000   11.000   10.000   12.000   13.000   14.000   15.000   16.000 
  25.000   27.000   26.000   28.000   29.000   30.000   31.000   32.000 
  33.000   35.000   34.000   36.000   37.000   38.000   39.000   40.000 
  41.000   43.000   42.000   44.000   45.000   46.000   47.000   48.000 
  49.000   51.000   50.000   52.000   53.000   54.000   55.000   56.000 

After daxpy(3.0, A(:,1), A(:,2))
   1.000    3.000   11.000    4.000    5.000    6.000    7.000    8.000 
  17.000   19.000   75.000   20.000   21.000   22.000   23.000   24.000 
   9.000   11.000   43.000   12.000   13.000   14.000   15.000   16.000 
  25.000   27.000  107.000   28.000   29.000   30.000   31.000   32.000 
  33.000   35.000  139.000   36.000   37.000   38.000   39.000   40.000 
  41.000   43.000  171.000   44.000   45.000   46.000   47.000   48.000 
  49.000   51.000  203.000   52.000   53.000   54.000   55.000   56.000 

After daxpy(-2.5, A(1,:), A(3,:))
   1.000    3.000   11.000    4.000    5.000    6.000    7.000    8.000 
  17.000   19.000   75.000   20.000   21.000   22.000   23.000   24.000 
   9.000   11.000   43.000   12.000   13.000   14.000   15.000   16.000 
 -17.500  -20.500  -80.500  -22.000  -23.500  -25.000  -26.500  -28.000 
  33.000   35.000  139.000   36.000   37.000   38.000   39.000   40.000 
  41.000   43.000  171.000   44.000   45.000   46.000   47.000   48.000 
  49.000   51.000  203.000   52.000   53.000   54.000   55.000   56.000 
theon$