#include <mkl.h>
#include <stdio.h>

// column major index within a 5x5 matrix
#define I(r,c) ((c)*5 + (r))

static void
print_mat(const double* vals, const char * name)
{
    printf("%s=[\n", name);
    for (int r = 0; r < 5; r++) {
        for (int c = 0; c < 5; c++) {
            printf("%s %.2f", c==0 ? "\t[" : ",", vals[I(r,c)]);
        }
        puts(r==4 ? "]]" : "],");
    }
}

int
main(void)
{
    // 5x5 matrix from http://software.intel.com/sites/products/documentation/hpc/mkl/mklman/GUID-9FCEB1C4-670D-4738-81D2-F378013412B0.htm#MKL_APPA_SMSF_3
    double values[13] = {1, -1, -3, -2, 5, 4, 6, 4, -4, 2, 7, 8, -5};
    MKL_INT columns[13] = {0, 1, 3, 0, 1, 2, 3, 4, 0, 2, 3, 1, 4};
    MKL_INT pointerB[5] = {0, 3, 5, 8, 11};
    MKL_INT pointerE[5] = {3, 5, 8, 11, 13};

    // copy the CSR representation into a dense matrix
    double a[25] = { };
    for (int r = 0; r < 5; ++r) {
        const MKL_INT begin = pointerB[r];
        const MKL_INT end = pointerE[r];
        for (int i = begin; i < end; ++i)
            a[I(columns[i],r)] = values[i];
    }

    double b[25] = {
            [I(0,0)] = 1,
            [I(1,0)] = 1,
            [I(2,0)] = 1,
            [I(3,0)] = 1,
            [I(4,0)] = 1,
    };

    double c[25] = { };

    print_mat(a, "A");
    print_mat(b, "B");

    char N = 'N';
    MKL_INT five = 5;
    double alpha = 1.0;
    double beta = 0.0;
    dgemm(&N, &N, &five, &five, &five, &alpha,
        a, &five,
        b, &five, &beta,
        c, &five);

    print_mat(c, "C(dgemm)");

    char matdescra[6] = {
            [0]='G', // general
            [3]='C', // zero based
    };
    mkl_dcsrmm(&N, &five, &five, &five, &alpha,
        matdescra,
        values, columns, pointerB, pointerE,
        b, &five, &beta,
        c, &five);

    print_mat(c, "C(mkl_dcsrmm)");

    return 0;
}