add double precision matrix implementation

6 files changed, 423 insertions, 40 deletions

diff --git a/matrix/generic/cukernel.cu b/matrix/generic/cukernel.cu
new file mode 100644
index 0000000..a37ccf4
--- /dev/null
+++ b/matrix/generic/cukernel.cu
@@ -0,0 +1,184 @@
+#ifdef NERV_GENERIC_CUKERNEL
+#include <assert.h>
+#include <stdio.h>
+#include "matrix.h"
+#include "cuda.h"
+#define CUDA_THREADS_N 16
+#define CUDA_THREADS_NN (16 * 16)
+#define CEIL_DIV(a, b) (((a) + (b) - 1) / (b))
+__global__ void cudak_(sigmoid)(const MATRIX_ELEM *a, MATRIX_ELEM *b,
+                        int nrow, int ncol, int stride) {
+    int j = blockIdx.x * blockDim.x + threadIdx.x;
+    int i = blockIdx.y * blockDim.y + threadIdx.y;
+    long idx;
+    if (i >= nrow || j >= ncol) return;
+    idx = j + i * stride;
+    b[idx] = 1.0 / (1.0 + exp(-a[idx]));
+}
+
+__global__ void cudak_(softmax_final)(const MATRIX_ELEM *a, MATRIX_ELEM *b,
+                        const MATRIX_ELEM *max, const MATRIX_ELEM *deno,
+                        int nrow, int ncol, int stride, int mstride) {
+    int j = blockIdx.x * blockDim.x + threadIdx.x;
+    int i = blockIdx.y * blockDim.y + threadIdx.y;
+    long idx;
+    if (i >= nrow || j >= ncol) return;
+    idx = j + i * stride;
+    b[idx] = exp(a[idx] - max[0 + i * mstride]) / deno[0 + i * mstride];
+}
+
+__global__ void cudak_(block_reduce_sum)(const MATRIX_ELEM *input,
+                                MATRIX_ELEM *output,
+                                const int istride, const int ostride,
+                                const int n) {
+    extern __shared__ MATRIX_ELEM cudak_(arr)[];
+    int j = blockIdx.x * blockDim.x + threadIdx.x;
+    cudak_(arr)[threadIdx.x] = j < n ? input[j + istride * blockIdx.y] : 0;
+    __syncthreads();
+    for (int offset = blockDim.x >> 1;  offset; offset >>= 1)
+    {
+        if (threadIdx.x < offset)
+            cudak_(arr)[threadIdx.x] += cudak_(arr)[threadIdx.x + offset];
+        __syncthreads();
+    }
+    if (threadIdx.x == 0)
+        output[blockIdx.x + ostride * blockIdx.y] = cudak_(arr)[0];
+}
+
+__global__ void cudak_(block_reduce_softmax_sum)(const MATRIX_ELEM *input,
+                                        MATRIX_ELEM *output,
+                                        const MATRIX_ELEM *max,
+                                        const int istride, const int ostride,
+                                        const int mstride, const int n) {
+    extern __shared__ MATRIX_ELEM cudak_(arr)[];
+    int j = blockIdx.x * blockDim.x + threadIdx.x;
+    cudak_(arr)[threadIdx.x] = j < n ? exp(input[j + istride * blockIdx.y] - \
+                                    max[0 + mstride * blockIdx.y]) : 0;
+    __syncthreads();
+    for (int offset = blockDim.x >> 1;  offset; offset >>= 1)
+    {
+        if (threadIdx.x < offset)
+            cudak_(arr)[threadIdx.x] += cudak_(arr)[threadIdx.x + offset];
+        __syncthreads();
+    }
+    if (threadIdx.x == 0)
+        output[blockIdx.x + ostride * blockIdx.y] = cudak_(arr)[0];
+}
+
+__global__ void cudak_(block_reduce_max)(const MATRIX_ELEM *input,
+                                MATRIX_ELEM *output,
+                                const int istride, const int ostride,
+                                const int n) {
+    extern __shared__ MATRIX_ELEM cudak_(arr)[];
+    int j = blockIdx.x * blockDim.x + threadIdx.x;
+    cudak_(arr)[threadIdx.x] = j < n ? input[j + istride * blockIdx.y] : 0;
+    __syncthreads();
+    for (int offset = blockDim.x >> 1;  offset; offset >>= 1)
+    {
+        if (threadIdx.x < offset)
+        {
+            MATRIX_ELEM l = cudak_(arr)[threadIdx.x],
+                  r = cudak_(arr)[threadIdx.x + offset];
+            if (r > l) cudak_(arr)[threadIdx.x] = r;
+        }
+        __syncthreads();
+    }
+    if (threadIdx.x == 0)
+        output[blockIdx.x + ostride * blockIdx.y] = cudak_(arr)[0];
+}
+
+extern "C" {
+#include "../cukernel.h"
+    void cudak_(cuda_sigmoid)(const Matrix *a, Matrix *b) {
+        dim3 threadsPerBlock(CUDA_THREADS_N,
+                CUDA_THREADS_N);
+        dim3 numBlocks(CEIL_DIV(b->ncol, threadsPerBlock.x),
+                CEIL_DIV(b->nrow, threadsPerBlock.y));
+        cudak_(sigmoid)<<<numBlocks, threadsPerBlock>>> \
+            (MATRIX_ELEM_PTR(a), MATRIX_ELEM_PTR(b), b->nrow, b->ncol,
+            b->stride / sizeof(MATRIX_ELEM));
+    }
+
+    void cudak_(cuda_colsum)(const Matrix *a, Matrix *b) {
+        dim3 block(CUDA_THREADS_NN, 1);
+        int ncol = a->ncol;
+        int blocks_per_row = CEIL_DIV(ncol, block.x);
+        dim3 grid(blocks_per_row, a->nrow);
+        MATRIX_ELEM *res;
+        size_t stride;
+        cudaMallocPitch(&res, &stride, blocks_per_row * sizeof(MATRIX_ELEM), a->nrow);
+        cudak_(block_reduce_sum)<<<grid, block, block.x * sizeof(MATRIX_ELEM)>>> \
+            (MATRIX_ELEM_PTR(a), res,
+             a->stride / sizeof(MATRIX_ELEM), stride / sizeof(MATRIX_ELEM),
+             ncol);
+        ncol = blocks_per_row;
+        assert((unsigned long)ncol <= block.x);
+        grid.x = 1;
+        cudak_(block_reduce_sum)<<<grid, block, block.x * sizeof(MATRIX_ELEM)>>> \
+            (res, MATRIX_ELEM_PTR(b),
+             stride / sizeof(MATRIX_ELEM), b->stride / sizeof(MATRIX_ELEM),
+             ncol);
+        cudaFree(res);
+    }
+
+    void cudak_(cuda_softmax_final)(const Matrix *a, const Matrix *max,
+                            const Matrix *deno, Matrix *b) {
+        dim3 threadsPerBlock(CUDA_THREADS_N,
+                CUDA_THREADS_N);
+        dim3 numBlocks(CEIL_DIV(b->ncol, threadsPerBlock.x),
+                CEIL_DIV(b->nrow, threadsPerBlock.y));
+        cudak_(softmax_final)<<<numBlocks, threadsPerBlock>>> \
+                (MATRIX_ELEM_PTR(a), MATRIX_ELEM_PTR(b),
+                MATRIX_ELEM_PTR(max), MATRIX_ELEM_PTR(deno),
+                b->nrow, b->ncol,
+                b->stride / sizeof(MATRIX_ELEM),
+                max->stride / sizeof(MATRIX_ELEM));
+    }
+
+    void cudak_(cuda_softmax_denominator)(const Matrix *a, const Matrix *max, Matrix *b) {
+        dim3 block(CUDA_THREADS_NN, 1);
+        int ncol = a->ncol;
+        int blocks_per_row = CEIL_DIV(ncol, block.x);
+        dim3 grid(blocks_per_row, a->nrow);
+        MATRIX_ELEM *res;
+        size_t stride;
+        assert(max->ncol == 1);
+        cudaMallocPitch(&res, &stride, blocks_per_row * sizeof(MATRIX_ELEM), a->nrow);
+        cudak_(block_reduce_softmax_sum)<<<grid, block, block.x * sizeof(MATRIX_ELEM)>>> \
+            (MATRIX_ELEM_PTR(a), res, MATRIX_ELEM_PTR(max),
+             a->stride / sizeof(MATRIX_ELEM), stride / sizeof(MATRIX_ELEM),
+             max->stride / sizeof(MATRIX_ELEM),
+             ncol);
+        ncol = blocks_per_row;
+        assert((unsigned long)ncol <= block.x);
+        grid.x = 1;
+        cudak_(block_reduce_sum)<<<grid, block, block.x * sizeof(MATRIX_ELEM)>>> \
+            (res, MATRIX_ELEM_PTR(b),
+             stride / sizeof(MATRIX_ELEM), b->stride / sizeof(MATRIX_ELEM),
+             ncol);
+        cudaFree(res);
+    }
+
+    void cudak_(cuda_colmax)(const Matrix *a, Matrix *b) {
+        dim3 block(CUDA_THREADS_NN, 1);
+        int ncol = a->ncol;
+        int blocks_per_row = CEIL_DIV(ncol, block.x);
+        dim3 grid(blocks_per_row, a->nrow);
+        MATRIX_ELEM *res;
+        size_t stride;
+        cudaMallocPitch(&res, &stride, blocks_per_row * sizeof(MATRIX_ELEM), a->nrow);
+        cudak_(block_reduce_max)<<<grid, block, block.x * sizeof(MATRIX_ELEM)>>> \
+            (MATRIX_ELEM_PTR(a), res,
+             a->stride / sizeof(MATRIX_ELEM), stride / sizeof(MATRIX_ELEM),
+             ncol);
+        ncol = blocks_per_row;
+        assert((unsigned long)ncol <= block.x);
+        grid.x = 1;
+        cudak_(block_reduce_max)<<<grid, block, block.x * sizeof(MATRIX_ELEM)>>> \
+            (res, MATRIX_ELEM_PTR(b),
+             stride / sizeof(MATRIX_ELEM), b->stride / sizeof(MATRIX_ELEM),
+             ncol);
+        cudaFree(res);
+    }
+}
+#endif
diff --git a/matrix/generic/cumatrix.c b/matrix/generic/cumatrix.c
new file mode 100644
index 0000000..f0ef99d
--- /dev/null
+++ b/matrix/generic/cumatrix.c
@@ -0,0 +1,143 @@
+#ifdef NERV_GENERIC_CUMATRIX
+#include "matrix.h"
+#include "elem_type.h"
+
+#define MATRIX_DATA_FREE(ptr) cuda_matrix_(free)(ptr)
+#define MATRIX_DATA_ALLOC(dptr, stride, width, height) \
+                            cuda_matrix_(alloc)(dptr, stride, width, height)
+#define MATRIX_DATA_WRITE(data, idx, val) cuda_matrix_(write)(data, idx, val)
+#define MATRIX_DATA_READ(data, idx) cuda_matrix_(read)(data, idx)
+#define MATRIX_INIT(L) cuda_matrix_(init)(L)
+#define NERV_GENERIC_MATRIX
+#define NERV_GENERIC_CUKERNEL
+#include "../../common.h"
+#include "../cukernel.h"
+#include "cuda.h"
+#include "cuda_runtime.h"
+#include "driver_types.h"
+#include "cublas_v2.h"
+
+static cublasHandle_t cublas_handle;
+
+Matrix *nerv_matrix_(new_)(long nrow, long ncol);
+static int nerv_matrix_(add)(lua_State *L) {
+    Matrix *a = luaT_checkudata(L, 1, nerv_matrix_(tname));
+    Matrix *b = luaT_checkudata(L, 2, nerv_matrix_(tname));
+    Matrix *c;
+    long nrow, ncol;
+    if (!(a->nrow == b->nrow && a->ncol == b->ncol))
+        nerv_error(L, "Matrices should be of the same dimension");
+    nrow = a->nrow;
+    ncol = a->ncol;
+    c = nerv_matrix_(new_)(nrow, ncol);
+    MATRIX_ELEM alpha = 1.0f, beta = 1.0f;
+    NERV_CUBLAS_(geam)(cublas_handle, CUBLAS_OP_N, CUBLAS_OP_N,
+                ncol, nrow,
+                &alpha,
+                MATRIX_ELEM_PTR(a), a->stride / sizeof(MATRIX_ELEM),
+                &beta,
+                MATRIX_ELEM_PTR(b), b->stride / sizeof(MATRIX_ELEM),
+                MATRIX_ELEM_PTR(c), c->stride / sizeof(MATRIX_ELEM));
+    luaT_pushudata(L, c, nerv_matrix_(tname));
+    return 1;
+}
+
+static int nerv_matrix_(mul)(lua_State *L) {
+    Matrix *a = luaT_checkudata(L, 1, nerv_matrix_(tname));
+    Matrix *b = luaT_checkudata(L, 2, nerv_matrix_(tname));
+    Matrix *c;
+    if (a->ncol != b->nrow)
+        nerv_error(L, "Wrong dimension of multipliers");
+    c = nerv_matrix_(new_)(a->nrow, b->ncol);
+    MATRIX_ELEM alpha = 1.0f, beta = 0.0f;
+    NERV_CUBLAS_(gemm)(cublas_handle, CUBLAS_OP_N, CUBLAS_OP_N,
+                b->ncol, a->nrow, b->nrow,
+                &alpha,
+                MATRIX_ELEM_PTR(b), b->stride / sizeof(MATRIX_ELEM),
+                MATRIX_ELEM_PTR(a), a->stride / sizeof(MATRIX_ELEM),
+                &beta,
+                MATRIX_ELEM_PTR(c), c->stride / sizeof(MATRIX_ELEM));
+    luaT_pushudata(L, c, nerv_matrix_(tname));
+    return 1;
+}
+
+static int nerv_matrix_(sigmoid)(lua_State *L) {
+    Matrix *a = luaT_checkudata(L, 1, nerv_matrix_(tname));
+    Matrix *b = nerv_matrix_(new_)(a->nrow, a->ncol);
+    cudak_(cuda_sigmoid)(a, b);
+    luaT_pushudata(L, b, nerv_matrix_(tname));
+    return 1;
+}
+
+static int nerv_matrix_(softmax)(lua_State *L) {
+    Matrix *a = luaT_checkudata(L, 1, nerv_matrix_(tname));
+    Matrix *max = nerv_matrix_(new_)(a->nrow, 1);
+    Matrix *dno = nerv_matrix_(new_)(a->nrow, 1);
+    Matrix *b = nerv_matrix_(new_)(a->nrow, a->ncol);
+    cudak_(cuda_colmax)(a, max);
+    cudak_(cuda_softmax_denominator)(a, max, dno);
+    cudak_(cuda_softmax_final)(a, max, dno, b);
+    luaT_pushudata(L, b, nerv_matrix_(tname));
+    return 1;
+}
+
+static int nerv_matrix_(colsum)(lua_State *L) {
+    Matrix *a = luaT_checkudata(L, 1, nerv_matrix_(tname));
+    Matrix *b = nerv_matrix_(new_)(a->nrow, 1);
+    cudak_(cuda_colsum)(a, b);
+    luaT_pushudata(L, b, nerv_matrix_(tname));
+    return 1;
+}
+
+static int nerv_matrix_(colmax)(lua_State *L) {
+    Matrix *a = luaT_checkudata(L, 1, nerv_matrix_(tname));
+    Matrix *b = nerv_matrix_(new_)(a->nrow, 1);
+    cudak_(cuda_colmax)(a, b);
+    luaT_pushudata(L, b, nerv_matrix_(tname));
+    return 1;
+}
+
+static const luaL_Reg nerv_matrix_(extra_methods)[] = {
+    {"__add__", nerv_matrix_(add)},
+    {"__mul__", nerv_matrix_(mul)},
+    {"sigmoid", nerv_matrix_(sigmoid)},
+    {"softmax", nerv_matrix_(softmax)},
+    {"colsum", nerv_matrix_(colsum)},
+    {"colmax", nerv_matrix_(colmax)},
+    {NULL, NULL}
+};
+
+static void cuda_matrix_(init)(lua_State *L) {
+    luaN_append_methods(L, nerv_matrix_(extra_methods));
+    cublasCreate(&cublas_handle);
+}
+
+static void cuda_matrix_(free)(MATRIX_ELEM *ptr) {
+    cudaFree(ptr);
+}
+
+static void cuda_matrix_(alloc)(MATRIX_ELEM **dptr, size_t *stride,
+                                long width, long height) {
+    cudaMallocPitch((void **)dptr, stride, width, height);
+}
+
+static MATRIX_ELEM cuda_matrix_(read)(MATRIX_ELEM *data, int idx) {
+    MATRIX_ELEM res;
+    cudaMemcpy(&res, data + idx, sizeof(MATRIX_ELEM), cudaMemcpyDeviceToHost);
+    return res;
+}
+
+static void cuda_matrix_(write)(MATRIX_ELEM *data, int idx, MATRIX_ELEM val) {
+    cudaMemcpy(data + idx, &val, sizeof(MATRIX_ELEM), cudaMemcpyHostToDevice);
+}
+
+int nerv_matrix_(get_elem)(lua_State *L) {
+    return nerv_error_method_not_implemented(L);
+}
+
+int nerv_matrix_(set_elem)(lua_State *L) {
+    return nerv_error_method_not_implemented(L);
+}
+
+#include "matrix.c"
+#endif
diff --git a/matrix/generic/elem_type.h b/matrix/generic/elem_type.h
new file mode 100644
index 0000000..8f80306
--- /dev/null
+++ b/matrix/generic/elem_type.h
@@ -0,0 +1,11 @@
+#ifdef MATRIX_USE_FLOAT
+
+#define MATRIX_ELEM float
+#define MATRIX_ELEM_PTR(self) ((self)->data.f)
+
+#elif defined(MATRIX_USE_DOUBLE)
+
+#define MATRIX_ELEM double
+#define MATRIX_ELEM_PTR(self) ((self)->data.d)
+
+#endif
diff --git a/matrix/generic/matrix.c b/matrix/generic/matrix.c
index 9ced397..f0f81a9 100644
--- a/matrix/generic/matrix.c
+++ b/matrix/generic/matrix.c
@@ -3,59 +3,61 @@
 #include "matrix.h"
 
 extern const char *nerv_matrix_tname;
-extern const char *nerv_float_matrix_(tname);
+extern const char *nerv_matrix_(tname);
 
-void nerv_float_matrix_(data_free)(Matrix *self) {
+void nerv_matrix_(data_free)(Matrix *self) {
     if (--(*self->data_ref) == 0)
-        MATRIX_DATA_FREE(self->data.f);
+        MATRIX_DATA_FREE(MATRIX_ELEM_PTR(self));
 }
 
-void nerv_float_matrix_(data_retain)(Matrix *self) {
+void nerv_matrix_(data_retain)(Matrix *self) {
     (*self->data_ref)++;
 }
 
-Matrix *nerv_float_matrix_(new_)(long nrow, long ncol) {
+Matrix *nerv_matrix_(new_)(long nrow, long ncol) {
     Matrix *self = (Matrix *)malloc(sizeof(Matrix));
     self->nrow = nrow;
     self->ncol = ncol;
     self->nmax = self->nrow * self->ncol;
-    MATRIX_DATA_ALLOC(&self->data.f, &self->stride, sizeof(float) * self->ncol, self->nrow);
+    MATRIX_DATA_ALLOC(&MATRIX_ELEM_PTR(self), &self->stride,
+                        sizeof(MATRIX_ELEM) * self->ncol, self->nrow);
     self->data_ref = (long *)malloc(sizeof(long));
     *self->data_ref = 0;
-    nerv_float_matrix_(data_retain)(self);
+    nerv_matrix_(data_retain)(self);
     return self;
 }
 
-int nerv_float_matrix_(new)(lua_State *L) {
-    luaT_pushudata(L, nerv_float_matrix_(new_)(luaL_checkinteger(L, 1),
+int nerv_matrix_(new)(lua_State *L) {
+    luaT_pushudata(L, nerv_matrix_(new_)(luaL_checkinteger(L, 1),
                                                 luaL_checkinteger(L, 2)),
-                    nerv_float_matrix_(tname));
+                    nerv_matrix_(tname));
     return 1;
 }
 
-int nerv_float_matrix_(destroy)(lua_State *L) {
-    Matrix *self = luaT_checkudata(L, 1, nerv_float_matrix_(tname));
-    nerv_float_matrix_(data_free)(self);
+int nerv_matrix_(destroy)(lua_State *L) {
+    Matrix *self = luaT_checkudata(L, 1, nerv_matrix_(tname));
+    nerv_matrix_(data_free)(self);
     return 0;
 }
 
-int nerv_float_matrix_(get_elem)(lua_State *L); 
-int nerv_float_matrix_(set_elem)(lua_State *L);
+int nerv_matrix_(get_elem)(lua_State *L); 
+int nerv_matrix_(set_elem)(lua_State *L);
 
-static Matrix *nerv_float_matrix_(getrow)(Matrix *self, int row) {
+static Matrix *nerv_matrix_(getrow)(Matrix *self, int row) {
     Matrix *prow = (Matrix *)malloc(sizeof(Matrix));
     prow->ncol = self->ncol;
     prow->nrow = 1;
     prow->stride = self->stride;
     prow->nmax = prow->ncol;
-    prow->data.f = (float *)((char *)self->data.f + row * self->stride);
+    MATRIX_ELEM_PTR(prow) = \
+        (MATRIX_ELEM *)((char *)MATRIX_ELEM_PTR(self) + row * self->stride);
     prow->data_ref = self->data_ref;
-    nerv_float_matrix_(data_retain)(self);
+    nerv_matrix_(data_retain)(self);
     return prow;
 }
 
-static int nerv_float_matrix_(newindex)(lua_State *L) {
-    Matrix *self = luaT_checkudata(L, 1, nerv_float_matrix_(tname));
+static int nerv_matrix_(newindex)(lua_State *L) {
+    Matrix *self = luaT_checkudata(L, 1, nerv_matrix_(tname));
     if (lua_isnumber(L, 2))
     {
         int idx = luaL_checkinteger(L, 2);
@@ -63,7 +65,8 @@ static int nerv_float_matrix_(newindex)(lua_State *L) {
         {
             if (idx < 0 || idx >= self->ncol)
                 nerv_error(L, "index must be within range [0, %d)", self->ncol);
-            MATRIX_DATA_WRITE(self->data.f, idx, luaL_checknumber(L, 3));
+            MATRIX_DATA_WRITE(MATRIX_ELEM_PTR(self), idx,
+                                luaL_checknumber(L, 3));
         }
         else
             nerv_error(L, "cannot assign a scalar to row vector");
@@ -78,8 +81,8 @@ static int nerv_float_matrix_(newindex)(lua_State *L) {
 }
 
 
-static int nerv_float_matrix_(index)(lua_State *L) {
-    Matrix *self = luaT_checkudata(L, 1, nerv_float_matrix_(tname));
+static int nerv_matrix_(index)(lua_State *L) {
+    Matrix *self = luaT_checkudata(L, 1, nerv_matrix_(tname));
     if (lua_isnumber(L, 2))
     {
         int idx = luaL_checkinteger(L, 2);
@@ -87,13 +90,13 @@ static int nerv_float_matrix_(index)(lua_State *L) {
         {
             if (idx < 0 || idx >= self->ncol)
                 nerv_error(L, "index must be within range [0, %d)", self->ncol);
-            lua_pushnumber(L, MATRIX_DATA_READ(self->data.f, idx));
+            lua_pushnumber(L, MATRIX_DATA_READ(MATRIX_ELEM_PTR(self), idx));
         }
         else
         {
             if (idx < 0 || idx >= self->nrow)
                 nerv_error(L, "index must be within range [0, %d)", self->nrow);
-            luaT_pushudata(L, nerv_float_matrix_(getrow)(self, idx), nerv_float_matrix_(tname));
+            luaT_pushudata(L, nerv_matrix_(getrow)(self, idx), nerv_matrix_(tname));
         }
         lua_pushboolean(L, 1);
         return 2;
@@ -105,33 +108,33 @@ static int nerv_float_matrix_(index)(lua_State *L) {
     }
 }
 
-static int nerv_float_matrix_(ncol)(lua_State *L) {
-    Matrix *self = luaT_checkudata(L, 1, nerv_float_matrix_(tname));
+static int nerv_matrix_(ncol)(lua_State *L) {
+    Matrix *self = luaT_checkudata(L, 1, nerv_matrix_(tname));
     lua_pushinteger(L, self->ncol);
     return 1;
 }
 
-static int nerv_float_matrix_(nrow)(lua_State *L) {
-    Matrix *self = luaT_checkudata(L, 1, nerv_float_matrix_(tname));
+static int nerv_matrix_(nrow)(lua_State *L) {
+    Matrix *self = luaT_checkudata(L, 1, nerv_matrix_(tname));
     lua_pushinteger(L, self->nrow);
     return 1;
 }
 
 
-static const luaL_Reg nerv_float_matrix_(methods)[] = {
-    {"get_elem", nerv_float_matrix_(get_elem)},
-    {"set_elem", nerv_float_matrix_(set_elem)},
-    {"ncol", nerv_float_matrix_(ncol)},
-    {"nrow", nerv_float_matrix_(nrow)},
-    {"__index__", nerv_float_matrix_(index)},
-    {"__newindex__", nerv_float_matrix_(newindex)},
+static const luaL_Reg nerv_matrix_(methods)[] = {
+    {"get_elem", nerv_matrix_(get_elem)},
+    {"set_elem", nerv_matrix_(set_elem)},
+    {"ncol", nerv_matrix_(ncol)},
+    {"nrow", nerv_matrix_(nrow)},
+    {"__index__", nerv_matrix_(index)},
+    {"__newindex__", nerv_matrix_(newindex)},
     {NULL, NULL}
 };
 
-void nerv_float_matrix_(init)(lua_State *L) {
-    luaT_newmetatable(L, nerv_float_matrix_(tname), nerv_matrix_tname,
-                        nerv_float_matrix_(new), nerv_float_matrix_(destroy), NULL);
-    luaL_register(L, NULL, nerv_float_matrix_(methods));
+void nerv_matrix_(init)(lua_State *L) {
+    luaT_newmetatable(L, nerv_matrix_(tname), nerv_matrix_tname,
+                        nerv_matrix_(new), nerv_matrix_(destroy), NULL);
+    luaL_register(L, NULL, nerv_matrix_(methods));
 #ifdef MATRIX_INIT
     MATRIX_INIT(L);
 #endif
diff --git a/matrix/generic/matrix.h b/matrix/generic/matrix.h
index 264859b..276ca5c 100644
--- a/matrix/generic/matrix.h
+++ b/matrix/generic/matrix.h
@@ -1,6 +1,7 @@
 #ifndef NERV_GENERIC_MATRIX_H
 #define NERV_GENERIC_MATRIX_H
 
+#include <stddef.h>
 typedef struct Matrix {
     size_t stride;              /* size of a row */
     long ncol, nrow, nmax;    /* dimension of the matrix */
diff --git a/matrix/generic/mmatrix.c b/matrix/generic/mmatrix.c
new file mode 100644
index 0000000..ac71c3d
--- /dev/null
+++ b/matrix/generic/mmatrix.c
@@ -0,0 +1,41 @@
+#ifdef NERV_GENERIC_MMATRIX
+#include "matrix.h"
+#include "elem_type.h"
+#define MATRIX_DATA_FREE(ptr) free(ptr)
+#define MATRIX_DATA_ALLOC(dptr, stride, width, height) \
+                            host_matrix_(alloc)(dptr, stride, width, height)
+#define MATRIX_DATA_STRIDE(ncol) (sizeof(float) * (ncol))
+#define MATRIX_DATA_WRITE(data, idx, val) (data[idx] = val)
+#define MATRIX_DATA_READ(data, idx) (data[idx])
+#define NERV_GENERIC_MATRIX
+#include "../../common.h"
+
+const char *nerv_matrix_(tname) = "nerv.FloatMMatrix";
+
+static void host_matrix_(alloc)(float **dptr, size_t *stride,
+                                    long width, long height) {
+    *dptr = (float *)malloc(width * height);
+    *stride = width;
+}
+
+int nerv_matrix_(get_elem)(lua_State *L) {
+    Matrix *self = luaT_checkudata(L, 1, nerv_matrix_(tname));
+    int idx = luaL_checkinteger(L, 2);
+    if (idx < 0 || idx >= self->nmax)
+        nerv_error(L, "index must be within range [0, %d)", self->nmax);
+    lua_pushnumber(L, self->data.f[idx]);
+    return 1;
+}
+
+int nerv_matrix_(set_elem)(lua_State *L) {
+    Matrix *self = luaT_checkudata(L, 1, nerv_matrix_(tname));
+    int idx = luaL_checkinteger(L, 2);
+    float v = luaL_checknumber(L, 3);
+    if (idx < 0 || idx >= self->nmax)
+        nerv_error(L, "index must be within range [0, %d)", self->nmax);
+    self->data.f[idx] = v;
+    return 0;
+}
+
+#include "matrix.c"
+#endif