#include <ulmaux.h>

#include <stdio.h>          // for printf()
#include <stdlib.h>         // for malloc(), free(), rand(), srand()
#include <math.h>           // for nan(), fabs()
#include <stdbool.h>        // for typedef bool
#include <sys/times.h>      // needed for walltime()
#include <unistd.h>         // needed for walltime()

//-- Functions for benchmarking and testing ------------------------------------

double
walltime()
{
   struct tms    ts;
   static double ClockTick=0.0;

   if (ClockTick==0.0) {
        ClockTick = 1.0 / ((double) sysconf(_SC_CLK_TCK));
   }
   return ((double) times(&ts)) * ClockTick;
}

void
randDGeMatrix(size_t m, size_t n, bool withNan,
              double *A,
              ptrdiff_t incRowA, ptrdiff_t incColA)
{
    for (size_t i=0; i<m; ++i) {
        for (size_t j=0; j<n; ++j) {
            A[i*incRowA + j*incColA] = withNan
                                     ? nan("")
                                     : 2.*(rand()-RAND_MAX/2)/RAND_MAX;
        }
    }
}

void
printfDGeMatrix(const char * fmt, 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(fmt, A[i*incRowA + j*incColA]);
        }
        printf("\n");
    }
    printf("\n");
}

void
printDGeMatrix(size_t m, size_t n,
               const double *A,
               ptrdiff_t incRowA, ptrdiff_t incColA)
{
    printfDGeMatrix("%9.2lf ", m, n, A, incRowA, incColA);
}

void
printIGeMatrix(size_t m, size_t n,
               const size_t *A,
               ptrdiff_t incRowA, ptrdiff_t incColA)
{
    for (size_t i=0; i<m; ++i) {
        for (size_t j=0; j<n; ++j) {
            printf("%6zu ", A[i*incRowA + j*incColA]);
        }
        printf("\n");
    }
    printf("\n");
}

extern inline size_t
ptrdiff_abs(ptrdiff_t x);

void
dfill_nan(size_t bufsize, double *buf)
{
    for (size_t i=0; i<bufsize; ++i) {
        buf[i] = nan("");
    }
}

void
ifill_rand(size_t bufsize, size_t *buf)
{
    for (size_t i=0; i<bufsize; ++i) {
        buf[i] = rand();
    }
}

//-- BLAS extensions -----------------------------------------------------------

void
dgecopy(size_t m, size_t n,
        const double *A, ptrdiff_t incRowA, ptrdiff_t incColA,
        double *B, ptrdiff_t incRowB, ptrdiff_t incColB)
{
    if (incRowB>incColB && incRowA>incColA) {
        dgecopy(n, m, A, incColA, incRowA, B, incColB, incRowB);
        return;
    }

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

void
dgeaxpy(size_t m, size_t n, double alpha,
        const double *X, ptrdiff_t incRowX, ptrdiff_t incColX,
        double *Y, ptrdiff_t incRowY, ptrdiff_t incColY)
{
    if (incRowY>incColY && incRowX>incColX) {
        dgeaxpy(n, m, alpha, X, incColX, incRowX, Y, incColY, incRowY);
        return;
    }

    for (size_t j=0; j<n; ++j) {
        for (size_t i=0; i<m; ++i) {
            Y[i*incRowY + j*incColY] += alpha * X[i*incRowX + j*incColX];
        }
    }
}


double
dgenrm_inf(size_t m, size_t n,
           const double *A,
           ptrdiff_t incRowA, ptrdiff_t incColA)
{
    double result = 0;

    for (size_t i=0; i<m; ++i) {
        double sum = 0;

        for (size_t j=0; j<n; ++j) {
            sum += fabs(A[i*incRowA+j*incColA]);
        }

        if (sum>result) {
            result = sum;
        }
    }

    return result;
}