# Sample solutions

## 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 #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

#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\$