The Kitchin Research Group

(vector-shape [[1 2 3]
               [3 4 5]])

(vector-ncols [[1 2 3]
               [3 4 5]])

#include "emacs-module.h"
#include "emacs-module-helpers.h"
#include <mkl.h>

int plugin_is_GPL_compatible;

emacs_value Fmkl_dgemm (emacs_env *env, ptrdiff_t nargs, emacs_value args[], void *data)
{
  double *A, *B, *C;
  int m, n, k, i, j;
  double alpha = 1.0;
  double beta = 0.0;
  
  // These will be 2d vectors
  emacs_value M0 = args[0]; // array 1 - A (m x k)
  emacs_value M1 = args[1]; // array 2 - B (k x n)

  // I need to get the number of rows and columns of each one.
  m = env->vec_size(env, M0);
  k  = 0;
  // We assume that we have a 2d array.
  emacs_value el1 = env->vec_get (env, M0, 0);
  k = env->vec_size(env, el1);
  
  // Now we know A has dimensions (m, k)
 
  emacs_value el2 = env->vec_get (env, M1, 0);
  n = env->vec_size(env, el2);
  
  // Now we know M1 had dimensions (k, n)
  
  // Now we have to build up arrays.
  // We are looking at a * b = c
  A = (double *)mkl_malloc( m*k*sizeof( double ), 64 );
  B = (double *)mkl_malloc( k*n*sizeof( double ), 64 );
  C = (double *)mkl_malloc( m*n*sizeof( double ), 64 );
  if (A == NULL || B == NULL || C == NULL) {
    printf( "\n ERROR: Can't allocate memory for matrices. Aborting... \n\n");
    mkl_free(A);
    mkl_free(B);
    mkl_free(C);
    return 1;
  }

  //populate A
  emacs_value row, ij;
  for (int i = 0; i < m; i++)
    {
      row = env->vec_get(env, M0, i);
      for (int j = 0; j < k; j++)
        {
          // get M0[i, j]
          ij = env->vec_get(env, row, j);
          A[k * i + j] = extract_double(env, ij);
        }
    }

  // populate B
  for (int i = 0; i < k; i++)
    {
      row = env->vec_get(env, M1, i);
      for (int j = 0; j < n; j++)
        {         
          // get M0[i, j]
          ij = env->vec_get(env, row, j);
          B[n * i + j] = extract_double(env, ij);
        }
    }

  // initialize C.  The solution will have dimensions of (rows1, cols2).
  for (int i = 0; i < m; i++)
    {
      for (int j = 0; j < n; j++)
        {
          C[n * i + j] = 0.0;
        }
    }

  // the multiplication is done here.
  cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, 
                m, n, k, alpha, A, k, B, n, beta, C, n);

  // now we build up the vector to return
  emacs_value vector = env->intern(env, "vector");
  emacs_value *array = malloc(sizeof(emacs_value) * m);
  emacs_value *row1;
  emacs_value vec;
  for (int i = 0; i < m; i++)
    {
      row1 = malloc(sizeof(emacs_value) * n);
      for (int j = 0; j < n; j++)
        {
          row1[j] = env->make_float(env, C[j + i * n]);
        }
      vec = env->funcall(env, vector, n, row1);
      array[i] = vec;
      free(row1);
    }

  emacs_value result = env->funcall(env, vector, m, array);
  free(array);
  return result;
}


int emacs_module_init(struct emacs_runtime *ert)
{
  emacs_env *env = ert->get_environment(ert);
  
  DEFUN("mkl-dgemm", Fmkl_dgemm, 2, 2,
        "(mkl-dgemm A B)\n"\
        "Multiply the matrices A and B. A and B must both be 2d vectors.\n" \
        "Returns the product as a vector.",
        NULL);
  provide(env, "mkl");
  
  return 0;
}

sh /opt/intel/mkl/bin/mklvars.sh intel64
gcc -Wall -m64 -I${MKLROOT}/include -fPIC -c mkl.c
gcc -shared -L${MKLROOT}/lib -Wl,-rpath,${MKLROOT}/lib -lmkl_rt -lpthread -lm -ldl -L. -lemacs-module-helpers -o mkl.so mkl.o

(add-to-list 'load-path (expand-file-name "."))
(require 'mkl)
(require 'cl)
(require 'seq)

(defun vector-shape (vec)
  "Return a vector of the shape of VEC."
  (let ((shape (vector (length vec))))
    (if (vectorp (aref vec 0))
        (vconcat shape (vector-shape (aref vec 0)))
      shape)))

(defun vector-ndims (vec)
  "Returns the number of dimensions in VEC."
  (length (vector-shape vec)))


(defun vector-numel (vec)
  "Returns the number of elements in VEC."
  (if (> (length vec) 0)
      (seq-reduce '* (vector-shape vec) 1)
    0))


(defun vector-nrows (vec)
 "Return the number of rows in VEC."
 (aref (vector-shape vec) 0))


(defun vector-ncols (vec)
 "Return the number of columns in VEC."
 (aref (vector-shape vec) 1))


(defun matrix-multiply (A B)
  "Return A * B in the matrix multiply sense."
  (cond
   ;; 1d * 1d i.e. a dot-product
   ((and (= 1 (vector-ndims A))
         (= 1 (vector-ndims B))
         (= (length A) (length B)))
    ;; this is easy to compute so we don't use dgemm.
    (seq-reduce '+ (mapcar* (lambda (a b) (* a b)) A B) 0.0))

   ;; 2d * 1d (m1, n1) * (n2, 1)
   ((and (= 2 (vector-ndims A))
         (= 1 (vector-ndims B))
         ;; ncols-A = len-B
         (= (vector-ncols A) (length B)))
    ;; transform B into a 2d column vector
    (let* ((B2d (apply 'vector (mapcar 'vector B)))
           (result  (mkl-dgemm A B2d)))
      ;; Now call (dgemm A B2d) -> (m2, 1) column vector
      ;; and convert it back to a 1d result
      (cl-map 'vector (lambda (v) (aref v 0)) result)))

   ;; 1d * 2d (1, n1) * (m2, n2) len-A = nrows-B
   ((and (= 1 (vector-ndims A))
         (= 2 (vector-ndims B))
         (= (length A) (vector-nrows B)))
    ;; transform B into a 2d row vector
    (let* ((A2d (vector A))
           (result  (mkl-dgemm A2d B)))
      ;; should be a 2d row vector
      (aref result 0)))

   ;; 2d * 2d (m1, n1) * (m2, n2) rows-A = ncols-B
   ((and (= 2 (vector-ndims A))
         (= 2 (vector-ndims B))
         (= (vector-ncols A)
            (vector-nrows B)))
    ;; call (dgemm A B) and return result
    (mkl-dgemm A B))
   (t
    ;; Error checking, getting here means none of the cases above were caught.
    ;; something is probably wrong.
    (cond
     ((or (> (vector-ndims A) 2)
          (> (vector-ndims B) 2))
      (error "One of your arrays has more than 2 dimensions. Only 1 or 2d arrays are supported"))
     ((and (= 1 (vector-ndims A))
           (= 1 (vector-ndims B))
           (not (= (length A) (length B))))
      (error "A and B must be the same length.
len(A) = %d
len(B) = %d" (length A) (length B)))
     ((and
       (= (vector-ndims A) 2)
       (= (vector-ndims B) 2)
       (not (= (vector-nrows A) (vector-ncols B))))
      (error "Your array shapes are not correct.
The number of rows in array A must equal the number of columns in array B.
There are %d rows in A and %d columns in B" (vector-nrows A) (vector-ncols B)))
     ((and
       (= (vector-ndims A) 2)
       (= (vector-ndims B) 1)
       (not (= (vector-nrows A) (length B))))
      (error "Your array shapes are not correct.
The number of rows in array A must equal the number of columns in array B.
There are %d rows in A and %d columns in B" (vector-nrows A) (length B)))
     (t
      (error "Unknown error"))))))