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

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

void
printMatrix(size_t m, size_t n,
            const double *A,
            ptrdiff_t incRow, ptrdiff_t incCol)
{
    for (size_t i=0; i<m; ++i) {
        for (size_t j=0; j<n; ++j) {
            printf("%10.3lf", A[i*incRow+j*incCol]);
        }
        printf("\n");
    }
    printf("\n");
}

//-- BLAS Level 1 --------------------------------------------------------------

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

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

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

void
dscal(size_t n, double alpha,
      double *x, ptrdiff_t incX)
{
    if (alpha==1) {
        return;
    }
    if (alpha==0) {
        for (size_t i=0; i<n; ++i) {
            x[i*incX] = 0;
        }
    } else {
        for (size_t i=0; i<n; ++i) {
            x[i*incX] *= alpha;
        }
    }
}

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

//-- BLAS Level 2 --------------------------------------------------------------

void
dger(size_t m, size_t n, double alpha,
     const double *x, ptrdiff_t incX,
     const double *y, ptrdiff_t incY,
     double *A, ptrdiff_t incRowA, ptrdiff_t incColA)
{
    if (alpha==0 || m==0 || n==0) {
        return;
    }
    if (incRowA<incColA) {
        // col major case
        for (size_t j=0; j<n; ++j) {
            daxpy(m, alpha*y[j*incY], x, incX, &A[j*incColA], incRowA);
        }
    } else {
        // row major case
        for (size_t i=0; i<m; ++i) {
            daxpy(n, alpha*x[i*incX], y, incY, &A[i*incRowA], incColA);
        }
    }
}

int
main()
{
    // Allocate matrix A_ with dimension m_ x n_
    size_t m_ = 5;
    size_t n_ = 8;

    size_t incRow = n_;
    size_t incCol = 1;

    double *A_ = (double *) malloc(m_*n_*sizeof(double));

    // Initialize matrix A_
    initMatrix(m_, n_, A_, incRow, incCol);

    // Print initial matrix A_
    printf("A_=\n");
    printMatrix(m_, n_, A_, incRow, incCol);

    // Matrix A refers to A_(1:m-1, 1:n-1)
    size_t    m       = m_ - 1;
    size_t    n       = n_ - 1;
    double    *A      = &A_[1*incRow+1*incCol];
    ptrdiff_t incRowA = incRow;
    ptrdiff_t incColA = incCol;

    // Vector x refers to A_(0, 1:n-1)
    double    *x      = &A_[0*incRow+1*incCol];
    ptrdiff_t incX    = incCol;

    // Vector y refers to A_(1:m-1, 0)
    double    *y      = &A_[1*incRow+0*incCol];
    ptrdiff_t incY    = incRow;

    // print A, x, y (for saving space: display vectors as row vectors)
    printf("A=\n");
    printMatrix(m, n, A, incRowA, incColA);
    printf("x=\n");
    printMatrix(1, n, x, 1, incX);
    printf("y=\n");
    printMatrix(1, m, y, 1, incY);

    // Some simple tests for the BLAS Level 2 functions

    // 1) dger
    printf("1.5*y * x^T + A -> A\n");
    dger(m, n, 1.5, y, incY, x, incX, A, incRowA, incColA);
    printMatrix(m, n, A, incRow, incCol);

    // 2) dger
    printf("1.5*x * y^T + A^T -> A^T\n");
    dger(n, m, 1.5, x, incX, y, incY, A, incColA, incRowA);
    printMatrix(m, n, A, incRow, incCol);

    // Deallocate matrix A_
    free(A_);
}