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

//-- setup and print matrices --------------------------------------------------

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

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

    for (i=0; i<m; ++i) {
        for (j=0; j<n; ++j) {
            if (unitDiag && (i==j)) {
                printf("%10.4lf ", 1.0);
            } else if ((lower && (i>=j)) || (!lower && (i<=j))) {
                printf("%10.4lf ", A[i*incRowA+j*incColA]);
            } else {
                printf("%10.4lf ", 0.0);
            }
        }
        printf("\n");
    }
    printf("\n\n");
}


//-- some BLAS Level 1 procedures and functions --------------------------------

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

//-- TRMV implementations ------------------------------------------------------

void
dtrlmv_col(int n, int unitDiag,
           const double *A, int incRowA, int incColA,
           double *x, int incX)
{
    // your code here
}

void
dtrlmv_row(int n, int unitDiag,
           const double *A, int incRowA, int incColA,
           double *x, int incX)
{
    // your code here
}

//------------------------------------------------------------------------------

#ifndef DIM_N
#define DIM_N 7
#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_N
#endif

#ifndef INC_X
#define INC_X 1
#endif

double A[DIM_N*DIM_N];
double x_[DIM_N];
double x[INC_X*DIM_N];

int
main()
{
    initGeMatrix(DIM_N, DIM_N, A, INCROW_A, INCCOL_A);
    initGeMatrix(DIM_N, 1, x_, 1, 0);

    //
    // Lower triangular matrix with unit diagonal
    //
    printf("L_1(A) =\n");
    printTrMatrix(DIM_N, DIM_N, 1, 1, A, INCROW_A, INCCOL_A);

    printf("x =\n");
    printGeMatrix(DIM_N, 1, x_, 1, 0);

    //
    //  Aufruf von trlmv_col
    //
    printf("trlmv_col: A * x =\n\n");

    dcopy(DIM_N, x_, 1, x, INC_X);
    dtrlmv_col(DIM_N, 1,
               A, INCROW_A, INCCOL_A,
               x, INC_X);

    printGeMatrix(DIM_N, 1, x, INC_X, 0);

    //
    //  Aufruf von trlmv_row
    //
    printf("trlmv_row: A * x =\n\n");

    dcopy(DIM_N, x_, 1, x, INC_X);
    dtrlmv_row(DIM_N, 1,
               A, INCROW_A, INCCOL_A,
               x, INC_X);

    printGeMatrix(DIM_N, 1, x, INC_X, 0);

    //
    // Lower triangular matrix with non-unit diagonal
    //
    printf("L(A) =\n");
    printTrMatrix(DIM_N, DIM_N, 0, 1, A, INCROW_A, INCCOL_A);

    printf("x =\n");
    printGeMatrix(DIM_N, 1, x_, 1, 0);

    //
    //  Aufruf von trlmv_col
    //
    printf("trlmv_col: A * x =\n\n");

    dcopy(DIM_N, x_, 1, x, INC_X);
    dtrlmv_col(DIM_N, 0,
               A, INCROW_A, INCCOL_A,
               x, INC_X);

    printGeMatrix(DIM_N, 1, x, INC_X, 0);

    //
    //  Aufruf von trlmv_row
    //
    printf("trlmv_row: A * x =\n\n");

    dcopy(DIM_N, x_, 1, x, INC_X);
    dtrlmv_row(DIM_N, 0,
               A, INCROW_A, INCCOL_A,
               x, INC_X);

    printGeMatrix(DIM_N, 1, x, INC_X, 0);
    return 0;
}