#include <stdio.h>
#include <math.h>

void
initGeMatrix(int m, int n, double *A, int incRowA, int incColA)
{
    int i, j;

    for (i=0; i<m; ++i) {
        for (j=0; j<n; ++j) {
            A[i*incRowA+j*incColA] = i*n + j + 1;
        }
    }
}

void
printGeMatrix(int m, int n, const double *A, int incRowA, int incColA)
{
    int i, j;

    for (i=0; i<m; ++i) {
        for (j=0; j<n; ++j) {
            printf("%10.4lf ", A[i*incRowA+j*incColA]);
        }
        printf("\n");
    }
    printf("\n\n");
}

double
ddot(int n, const double *x, int incX, const double *y, int incY)
{
    int     i;
    double  alpha = 0;

    for (i=0; i<n; ++i) {
        alpha += x[i*incX]*y[i*incY];

    }
    return alpha;
}

void
daxpy(int n, double alpha, const double *x, int incX, double *y, int incY)
{
    int i;

    if (alpha==0) {
        return;
    }
    for (i=0; i<n; ++i) {
        y[i*incY] += alpha*x[i*incX];
    }
}

void
dscal(int n, double alpha, double *x, int incX)
{
    int i;

    if (alpha==1) {
        return;
    }
    for (i=0; i<n; ++i) {
        x[i*incX] *= alpha;
    }
}

void
dcopy(int n, const double *x, int incX, double *y, int incY)
{
    int i;

    for (i=0; i<n; ++i) {
        y[i*incY] = x[i*incX];
    }
}

void
dswap(int n, double *x, int incX, double *y, int incY)
{
    int i;

    for (i=0; i<n; ++i) {
        double tmp;

        tmp = x[i*incX];
        x[i*incX] = y[i*incY];
        y[i*incY] = tmp;
    }
}

double
damax(int n, const double *x, int incX)
{
    int     i;
    double  result = 0;

    for (i=0; i<n; ++i) {
        if (fabs(x[i*incX])>result) {
            result = fabs(x[i*incX]);
        }
    }
    return result;
}

double
dgenrm1(int m, int n, const double *A, int incRowA, int incColA)
{
    int     i, j;
    double  result = 0;

    for (j=0; j<n; ++j) {
        double sum = 0;
        for (i=0; i<m; ++i) {
            sum += fabs(A[i*incRowA+j*incColA]);
        }
        if (sum>result) {
            result = sum;
        }
    }
    return result;
}

#ifndef DIM_M
#define DIM_M 5
#endif

#ifndef DIM_N
#define DIM_N 5
#endif

#ifndef ROWMAJOR
#define ROWMAJOR 0
#endif

#if (ROWMAJOR==1)
#   define INCROW_A  DIM_N
#   define INCCOL_A  1
#else
#   define INCROW_A  1
#   define INCCOL_A  DIM_M
#endif

double A[DIM_M*DIM_N];

int
main()
{
    initGeMatrix(DIM_M, DIM_N, A, INCROW_A, INCCOL_A);

    printf("A =\n");
    printGeMatrix(DIM_M, DIM_N, A, INCROW_A, INCCOL_A);


    printf("A(:,1) * A(:,2) = %lf\n",
           ddot(DIM_M, &A[0*INCCOL_A], INCROW_A, &A[1*INCCOL_A], INCROW_A));


    daxpy(DIM_M, 2.0, &A[1*INCCOL_A], INCROW_A, &A[2*INCCOL_A], INCROW_A);
    printf("A(:,3) = A(:3) + 2*A(:,2)\n");
    printGeMatrix(DIM_M, DIM_N, A, INCROW_A, INCCOL_A);


    dscal(DIM_M, 0.5, &A[3*INCCOL_A], INCROW_A);
    printf("A(:,4) = 0.5*A(:4)\n");
    printGeMatrix(DIM_M, DIM_N, A, INCROW_A, INCCOL_A);

    printf("inf-nrm(A(:,1)) = %lf\n", damax(DIM_M, &A[0*INCCOL_A], INCROW_A));
    printf("inf-nrm(A(2,:)) = %lf\n", damax(DIM_N, &A[1*INCROW_A], INCCOL_A));

    printf("nrm1(A(1:4,1:5)) = %lf\n",
           dgenrm1(4, 5, A, INCROW_A, INCCOL_A));
    printf("nrm1(A(2:4,3:5)) = %lf\n",
           dgenrm1(3, 3, &A[1*INCROW_A+2*INCCOL_A], INCROW_A, INCCOL_A));

    dswap(DIM_M, &A[1*INCCOL_A], INCROW_A, &A[2*INCCOL_A], INCROW_A);
    printf("swap(A(:,2), A(:,3))\n");
    printGeMatrix(DIM_M, DIM_N, A, INCROW_A, INCCOL_A);

    dcopy(DIM_M, &A[2*INCCOL_A], INCROW_A, &A[3*INCCOL_A], INCROW_A);
    printf("A(:,4) = A(:,3)\n");
    printGeMatrix(DIM_M, DIM_N, A, INCROW_A, INCCOL_A);


    return 0;
}