#include <ulmblas/level3.h>
#include <ulmblas/level2.h>
#include <ulmblas/level1.h>

//
// BLAS Level 3
//

//-- GEMM ----------------------------------------------------------------------

void
dgemm_mv(int m, int n, int k,
         double alpha,
         const double *A, int incRowA, int incColA,
         const double *B, int incRowB, int incColB,
         double beta,
         double *C, int incRowC, int incColC)
{
    int j;

    for (j=0; j<n; ++j) {
        dgemv(m, k,
              alpha,
              A, incRowA, incColA,
              &B[j*incColB], incRowB,
              beta,
              &C[j*incColC], incRowC);
    }
}

void
dgemm_jil(int m, int n, int k,
          double alpha,
          const double *A, int incRowA, int incColA,
          const double *B, int incRowB, int incColB,
          double beta,
          double *C, int incRowC, int incColC)
{
    int i, j, l;

    for (j=0; j<n; ++j) {
        for (i=0; i<m; ++i) {
            C[i*incRowC+j*incColC] *= beta;
            for (l=0; l<k; ++l) {
                C[i*incRowC+j*incColC] += alpha*A[i*incRowA+l*incColA]
                                               *B[l*incRowB+j*incColB];
            }
        }
    }
}

void
dgemm_jli(int m, int n, int k,
          double alpha,
          const double *A, int incRowA, int incColA,
          const double *B, int incRowB, int incColB,
          double beta,
          double *C, int incRowC, int incColC)
{
    int i, j, l;

    for (j=0; j<n; ++j) {
        for (l=0; l<k; ++l) {
            for (i=0; i<m; ++i) {
                if (l==0) {
                    C[i*incRowC+j*incColC] *= beta;
                }
                C[i*incRowC+j*incColC] += alpha*A[i*incRowA+l*incColA]
                                               *B[l*incRowB+j*incColB];
            }
        }
    }
}

void
dgemm_vm(int m, int n, int k,
         double alpha,
         const double *A, int incRowA, int incColA,
         const double *B, int incRowB, int incColB,
         double beta,
         double *C, int incRowC, int incColC)
{
    int i;

    for (i=0; i<m; ++i) {
        dgemv(n, k,
              alpha,
              B, incColB, incRowB,
              &A[i*incRowA], incColA,
              beta,
              &C[i*incRowC], incColC);
    }
}

void
dgemm_ilj(int m, int n, int k,
          double alpha,
          const double *A, int incRowA, int incColA,
          const double *B, int incRowB, int incColB,
          double beta,
          double *C, int incRowC, int incColC)
{
    int i, j, l;

    for (i=0; i<m; ++i) {
        for (l=0; l<k; ++l) {
            for (j=0; j<n; ++j) {
                if (l==0) {
                    C[i*incRowC+j*incColC] *= beta;
                }
                C[i*incRowC+j*incColC] += alpha*A[i*incRowA+l*incColA]
                                               *B[l*incRowB+j*incColB];
            }
        }
    }
}

void
dgemm_ijl(int m, int n, int k,
          double alpha,
          const double *A, int incRowA, int incColA,
          const double *B, int incRowB, int incColB,
          double beta,
          double *C, int incRowC, int incColC)
{
    int i, j, l;

    for (i=0; i<m; ++i) {
        for (j=0; j<n; ++j) {
            C[i*incRowC+j*incColC] *= beta;
            for (l=0; l<k; ++l) {
                C[i*incRowC+j*incColC] += alpha*A[i*incRowA+l*incColA]
                                               *B[l*incRowB+j*incColB];
            }
        }
    }
}

void
dgemm_vv(int m, int n, int k,
         double alpha,
         const double *A, int incRowA, int incColA,
         const double *B, int incRowB, int incColB,
         double beta,
         double *C, int incRowC, int incColC)
{
    int l;

    dgescal(m, n, alpha, C, incRowC, incColC);

    for (l=0; l<k; ++l) {
        dger(m, n, alpha,
             &A[l*incColA], incRowA,
             &B[l*incRowB], incColB,
             C, incRowC, incColC);
    }
}

void
dgemm_lij(int m, int n, int k,
          double alpha,
          const double *A, int incRowA, int incColA,
          const double *B, int incRowB, int incColB,
          double beta,
          double *C, int incRowC, int incColC)
{
    int i, j, l;

    for (j=0; j<n; ++j) {
        for (i=0; i<m; ++i) {
            C[i*incRowC+j*incColC] *= beta;
        }
    }
    for (l=0; l<k; ++l) {
        for (i=0; i<m; ++i) {
            for (j=0; j<n; ++j) {
                C[i*incRowC+j*incColC] += alpha*A[i*incRowA+l*incColA]
                                               *B[l*incRowB+j*incColB];
            }
        }
    }
}

void
dgemm_lji(int m, int n, int k,
          double alpha,
          const double *A, int incRowA, int incColA,
          const double *B, int incRowB, int incColB,
          double beta,
          double *C, int incRowC, int incColC)
{
    int i, j, l;

    for (j=0; j<n; ++j) {
        for (i=0; i<m; ++i) {
            C[i*incRowC+j*incColC] *= beta;
        }
    }
    for (l=0; l<k; ++l) {
        for (j=0; j<n; ++j) {
            for (i=0; i<m; ++i) {
                C[i*incRowC+j*incColC] += alpha*A[i*incRowA+l*incColA]
                                               *B[l*incRowB+j*incColB];
            }
        }
    }
}