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

    if (!A) {
        abort();
    }

    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;
}