add CuContext/MContextwrapped-handles

27 files changed, 730 insertions, 319 deletions

diff --git a/nerv/Makefile b/nerv/Makefile
index c0db53a..a2155b9 100644
--- a/nerv/Makefile
+++ b/nerv/Makefile
@@ -6,14 +6,15 @@ LIB_PATH := $(LUA_BINDIR)/../lib
 INC_PATH := $(LUA_BINDIR)/../include/nerv
 LUA_DIR = $(INST_LUADIR)/nerv
 OBJ_DIR := $(BUILD_DIR)/objs
-ISUBDIR := io matrix luaT
+ISUBDIR := lib matrix lib/io lib/matrix lib/luaT
 SUBDIR := matrix io layer examples nn tnn lib/io lib/luaT lib/matrix
 
 INC_SUBDIR := $(addprefix $(INC_PATH)/,$(ISUBDIR))
 OBJ_SUBDIR := $(addprefix $(OBJ_DIR)/,$(SUBDIR))
 LUA_SUBDIR := $(addprefix $(LUA_DIR)/,$(SUBDIR))
 
-INCS := common.h matrix/matrix.h io/chunk_file.h luaT/luaT.h
+INCS := lib/common.h lib/matrix/matrix.h lib/matrix/mmatrix.h lib/io/chunk_file.h lib/luaT/luaT.h \
+		matrix/matrix.h
 CORE_OBJS := lib/common.o lib/io/chunk_file.o \
 			lib/matrix/mmatrix.o lib/matrix/cumatrix.o lib/matrix/cukernel.o
 NERV_OBJS := nerv.o \
@@ -82,5 +83,5 @@ clean:
 
 install: $(LIBS) $(LUA_DIR) $(LUA_SUBDIR) $(LUA_LIBS) $(INC_SUBDIR) $(INCS)
 
-$(INC_PATH)/%.h: lib/%.h
+$(INC_PATH)/%.h: %.h
 	cp $< $@
diff --git a/nerv/init.lua b/nerv/init.lua
index 6312df1..e7d668c 100644
--- a/nerv/init.lua
+++ b/nerv/init.lua
@@ -54,7 +54,7 @@ end
 function nerv.warning(fmt, ...)
     nerv.printf(
         string.format("(%s)[nerv] warning: %s\n",
-            os.date("%H:%M:%S %F"), fmt), ...)
+            os.date("%H:%M:%S.%N %F"), fmt), ...)
 end
 
 --- Create a class (Torch-compatible).
diff --git a/nerv/lib/common.c b/nerv/lib/common.c
index d977f8d..879ae9d 100644
--- a/nerv/lib/common.c
+++ b/nerv/lib/common.c
@@ -56,7 +56,7 @@ void luaN_append_methods(lua_State *L, const luaL_Reg *mlist) {
     }
 }
 
-HashMap *hashmap_create(size_t size, HashKey_t hfunc, HashMapCmp_t cmp) {
+HashMap *nerv_hashmap_create(size_t size, HashKey_t hfunc, HashMapCmp_t cmp) {
     HashMap *res = (HashMap *)malloc(sizeof(HashMap));
     res->bucket = calloc(size, sizeof(HashNode));
     res->cmp = cmp;
@@ -65,7 +65,7 @@ HashMap *hashmap_create(size_t size, HashKey_t hfunc, HashMapCmp_t cmp) {
     return res;
 }
 
-void *hashmap_getval(HashMap *h, const char *key) {
+void *nerv_hashmap_getval(HashMap *h, const char *key) {
     size_t idx = h->hfunc(key) % h->size;
     HashNode *ptr;
     for (ptr = h->bucket[idx]; ptr; ptr = ptr->next)
@@ -76,7 +76,7 @@ void *hashmap_getval(HashMap *h, const char *key) {
     return NULL;
 }
 
-void hashmap_setval(HashMap *h, const char *key, void *val) {
+void nerv_hashmap_setval(HashMap *h, const char *key, void *val) {
     size_t idx = h->hfunc(key) % h->size;
     HashNode *ptr = malloc(sizeof(HashNode));
     ptr->next = h->bucket[idx];
@@ -85,7 +85,7 @@ void hashmap_setval(HashMap *h, const char *key, void *val) {
     ptr->val = val;
 }
 
-void hashmap_clear(HashMap *h) {
+void nerv_hashmap_clear(HashMap *h) {
     size_t i;
     for (i = 0; i < h->size; i++)
     {
@@ -100,6 +100,11 @@ void hashmap_clear(HashMap *h) {
     }
 }
 
+void nerv_hashmap_destroy(HashMap *h) {
+    nerv_hashmap_clear(h);
+    free(h);
+}
+
 size_t bkdr_hash(const char *key) {
     unsigned int seed = 131;
     unsigned int res = 0;
diff --git a/nerv/lib/common.h b/nerv/lib/common.h
index 1c588d1..3283ac1 100644
--- a/nerv/lib/common.h
+++ b/nerv/lib/common.h
@@ -7,6 +7,8 @@
 #include <stdio.h>
 #include <stdlib.h>
 
+#define PROFILE_HASHMAP_SIZE 123457
+
 typedef enum ErrCode {
     NERV_NORMAL,
     /* matrix err */
@@ -75,10 +77,11 @@ typedef struct HashMap {
     size_t size;
 } HashMap;
 
-HashMap *hashmap_create(size_t size, HashKey_t hfunc, HashMapCmp_t cmp);
-void *hashmap_getval(HashMap *h, const char *key);
-void hashmap_setval(HashMap *h, const char *key, void *val);
-void hashmap_clear(HashMap *h);
+HashMap *nerv_hashmap_create(size_t size, HashKey_t hfunc, HashMapCmp_t cmp);
+void *nerv_hashmap_getval(HashMap *h, const char *key);
+void nerv_hashmap_setval(HashMap *h, const char *key, void *val);
+void nerv_hashmap_clear(HashMap *h);
+void nerv_hashmap_destroy(HashMap *h);
 
 size_t bkdr_hash(const char *key);
 
diff --git a/nerv/lib/matrix/cuda_helper.h b/nerv/lib/matrix/cuda_helper.h
index 13d5728..5c75e38 100644
--- a/nerv/lib/matrix/cuda_helper.h
+++ b/nerv/lib/matrix/cuda_helper.h
@@ -54,6 +54,28 @@
         cudaDeviceSynchronize(); \
     } while (0)
 
+#define CURAND_SAFE_SYNC_CALL(call, status) \
+    do { \
+        curandStatus_t  err = (call); \
+        if (err != CURAND_STATUS_SUCCESS) \
+        { \
+            NERV_SET_STATUS(status, MAT_CUBLAS_ERR, curandGetErrorString(err)); \
+            return; \
+        } \
+        cudaDeviceSynchronize(); \
+    } while (0)
+
+#define CURAND_SAFE_SYNC_CALL_RET(call, status) \
+    do { \
+        curandStatus_t  err = (call); \
+        if (err != CURAND_STATUS_SUCCESS) \
+        { \
+            NERV_SET_STATUS(status, MAT_CUBLAS_ERR, curandGetErrorString(err)); \
+            return 0; \
+        } \
+        cudaDeviceSynchronize(); \
+    } while (0)
+
 #define CHECK_SAME_DIMENSION(a, b, status) \
     do { \
         if (!(a->nrow == b->nrow && a->ncol == b->ncol)) \
@@ -96,15 +118,46 @@ static const char *cublasGetErrorString(cublasStatus_t err) {
     return "<unknown>";
 }
 
+static const char *curandGetErrorString(curandStatus_t err) {
+    switch (err)
+    {
+        case CURAND_STATUS_VERSION_MISMATCH:
+            return "Header file and linked library version do not match";
+        case CURAND_STATUS_NOT_INITIALIZED:
+            return "Generator not initialized";
+        case CURAND_STATUS_ALLOCATION_FAILED:
+            return "Memory allocation failed";
+        case CURAND_STATUS_TYPE_ERROR:
+            return "Generator is wrong type";
+        case CURAND_STATUS_OUT_OF_RANGE:
+            return "Argument out of range";
+        case CURAND_STATUS_LENGTH_NOT_MULTIPLE:
+            return "Length requested is not a multple of dimension";
+        case CURAND_STATUS_DOUBLE_PRECISION_REQUIRED:
+            return "GPU does not have double precision required by MRG32k3a";
+        case CURAND_STATUS_LAUNCH_FAILURE:
+            return "Kernel launch failure";
+        case CURAND_STATUS_PREEXISTING_FAILURE:
+            return "Preexisting failure on library entry";
+        case CURAND_STATUS_INITIALIZATION_FAILED:
+            return "Initialization of CUDA failed";
+        case CURAND_STATUS_ARCH_MISMATCH:
+            return "Architecture mismatch, GPU does not support requested feature";
+        case CURAND_STATUS_INTERNAL_ERROR:
+            return "Internal library error";
+    }
+    return "<unknown>";
+}
 #define PROFILE_START \
     do { \
-        cudaEventRecord(profile_start, 0);
+        cudaEventRecord(context->profile_start, 0);
 #define PROFILE_STOP \
-        cudaEventRecord(profile_stop, 0); \
-        cudaEventSynchronize(profile_stop); \
+        cudaEventRecord(context->profile_stop, 0); \
+        cudaEventSynchronize(context->profile_stop); \
         float milliseconds = 0; \
-        cudaEventElapsedTime(&milliseconds, profile_start, profile_stop); \
-        accu_profile(__func__, milliseconds / 1000); \
+        cudaEventElapsedTime(&milliseconds, context->profile_start, \
+                                            context->profile_stop); \
+        nerv_cuda_context_accu_profile(context, __func__, milliseconds / 1000); \
     } while (0);
 
 #define PROFILE_END
diff --git a/nerv/lib/matrix/cukernel.h b/nerv/lib/matrix/cukernel.h
index c84200e..d59a070 100644
--- a/nerv/lib/matrix/cukernel.h
+++ b/nerv/lib/matrix/cukernel.h
@@ -3,7 +3,7 @@ void cudak_(cuda_mul_elem)(const Matrix *a, const Matrix *b, Matrix *c);
 void cudak_(cuda_log_elem)(const Matrix *a, Matrix *b);
 void cudak_(cuda_sigmoid)(const Matrix *a, Matrix *b);
 void cudak_(cuda_sigmoid_grad)(const Matrix *output, const Matrix *err, Matrix *nerr);
-void cudak_(cuda_rand_uniform)(const Matrix *a); /* a's curand_gen may be modified */
+void cudak_(cuda_rand_uniform)(const Matrix *a, CuContext *context); /* a's curand_gen may be modified */
 void cudak_(cuda_thres_mask)(const Matrix *a, const Matrix *b, double thres, double low, double high);
 void cudak_(cuda_tanh)(const Matrix *a, Matrix *b);
 void cudak_(cuda_tanh_grad)(const Matrix *output, const Matrix *err, Matrix *nerr);
diff --git a/nerv/lib/matrix/cumatrix.c b/nerv/lib/matrix/cumatrix.c
index d998871..2fbe7d8 100644
--- a/nerv/lib/matrix/cumatrix.c
+++ b/nerv/lib/matrix/cumatrix.c
@@ -1,23 +1,12 @@
 #define NERV_GENERIC_CUMATRIX
+#define MATRIX_CONTEXT CuContext
 #include "cumatrix.h"
 #include "cuda_helper.h"
 #include <string.h>
 #include <time.h>
-#define PROFILE_HASHMAP_SIZE 123457
-static cublasHandle_t cublas_handle;
-static cudaEvent_t profile_start, profile_stop;
-curandGenerator_t curand_gen;
-static HashMap *profile;
 
-void nerv_cumatrix_select_gpu(int dev, Status *status) {
-    fprintf(stderr, "** selecting GPU %d\n", dev);
-    NERV_SET_STATUS(status, NERV_NORMAL, 0);
-    CUDA_SAFE_SYNC_CALL(cudaSetDevice(dev), status);
-    CUDA_SAFE_SYNC_CALL(cublasDestroy(cublas_handle), status);
-    CUDA_SAFE_SYNC_CALL(cublasCreate(&cublas_handle), status);
-}
-
-void nerv_cumatrix_print_profile() {
+void nerv_cuda_context_print_profile(CuContext *context) {
+    HashMap *profile = context->profile;
     size_t i;
     fprintf(stderr, "*** [nerv cumatrix profile] **\n");
     for (i = 0; i < profile->size; i++)
@@ -30,28 +19,72 @@ void nerv_cumatrix_print_profile() {
     }
 }
 
-void nerv_cumatrix_clear_profile() {
-    hashmap_clear(profile);
+void nerv_cuda_context_clear_profile(CuContext *context) {
+    nerv_hashmap_clear(context->profile);
 }
 
-void accu_profile(const char *name, float delta) {
-    float *val = hashmap_getval(profile, name);
+void nerv_cuda_context_accu_profile(CuContext *context,
+                                    const char *name, float delta) {
+    HashMap *profile = context->profile;
+    float *val = nerv_hashmap_getval(profile, name);
     if (!val)
     {
         val = malloc(sizeof(float));
         *val = 0;
-        hashmap_setval(profile, name, val);
+        nerv_hashmap_setval(profile, name, val);
     }
     *val += delta;
 }
 
-void nerv_cumatrix_init() {
-    cublasCreate(&cublas_handle);
-    curandCreateGenerator(&curand_gen, CURAND_RNG_PSEUDO_DEFAULT);
-    curandSetPseudoRandomGeneratorSeed(curand_gen, time(NULL));
-    cudaEventCreate(&profile_start);
-    cudaEventCreate(&profile_stop);
-    profile = hashmap_create(PROFILE_HASHMAP_SIZE, bkdr_hash, strcmp);
+static void new_cuda_handles(CuContext *context, Status *status) {
+    CUBLAS_SAFE_SYNC_CALL(cublasCreate(&(context->cublas_handle)), status);
+    CURAND_SAFE_SYNC_CALL(curandCreateGenerator(&(context->curand_gen),
+                            CURAND_RNG_PSEUDO_DEFAULT), status);
+    CURAND_SAFE_SYNC_CALL(
+            curandSetPseudoRandomGeneratorSeed(context->curand_gen, time(NULL)),
+            status);
+    CUDA_SAFE_SYNC_CALL(cudaEventCreate(&(context->profile_start)), status);
+    CUDA_SAFE_SYNC_CALL(cudaEventCreate(&(context->profile_stop)), status);
+    NERV_SET_STATUS(status, NERV_NORMAL, 0);
+}
+
+static void free_cuda_handles(CuContext *context, Status *status) {
+    CUBLAS_SAFE_SYNC_CALL(cublasDestroy(context->cublas_handle), status);
+    CURAND_SAFE_SYNC_CALL(curandDestroyGenerator(context->curand_gen), status);
+    CUDA_SAFE_SYNC_CALL(cudaEventDestroy(context->profile_start), status);
+    CUDA_SAFE_SYNC_CALL(cudaEventDestroy(context->profile_stop), status);
+    NERV_SET_STATUS(status, NERV_NORMAL, 0);
+}
+
+CuContext *nerv_cuda_context_create(Status *status) {
+    CuContext *context = (CuContext *)malloc(sizeof(CuContext));
+    new_cuda_handles(context, status);
+    if (status->err_code != NERV_NORMAL)
+        return NULL;
+    context->profile = nerv_hashmap_create(PROFILE_HASHMAP_SIZE, bkdr_hash, strcmp);
+    NERV_SET_STATUS(status, NERV_NORMAL, 0);
+    return context;
+}
+
+void nerv_cuda_context_destroy(CuContext *context, Status *status) {
+    free_cuda_handles(context, status);
+    if (status->err_code != NERV_NORMAL)
+        return;
+    nerv_hashmap_destroy(context->profile);
+    free(context);
+    NERV_SET_STATUS(status, NERV_NORMAL, 0);
+}
+
+void nerv_cuda_context_select_gpu(CuContext *context,
+                                    int dev, Status *status) {
+    free_cuda_handles(context, status);
+    if (status->err_code != NERV_NORMAL)
+        return;
+    CUDA_SAFE_SYNC_CALL(cudaSetDevice(dev), status);
+    new_cuda_handles(context, status);
+    if (status->err_code != NERV_NORMAL)
+        return;
+    NERV_SET_STATUS(status, NERV_NORMAL, 0);
 }
 
 #define MATRIX_USE_FLOAT
@@ -59,7 +92,6 @@ void nerv_cumatrix_init() {
 #define nerv_matrix_(NAME) nerv_matrix_cuda_float_##NAME
 #define cudak_(NAME) cudak_float_ ## NAME
 #define NERV_CUBLAS_(NAME) cublasS##NAME
-#define MATRIX_CUMATRIX_HOST_TNAME nerv_matrix_host_float_tname
 #include "generic/cumatrix.c"
 
 #undef NERV_CUBLAS_
@@ -72,12 +104,10 @@ void nerv_cumatrix_init() {
 #undef MATRIX_ELEM_PTR_BASE
 #undef MATRIX_ELEM_FMT
 #undef MATRIX_ELEM_WRITE_FMT
-#undef MATRIX_CUMATRIX_HOST_TNAME
 
 #define MATRIX_USE_DOUBLE
 #define cuda_matrix_(NAME) cuda_matrix_double_##NAME
 #define nerv_matrix_(NAME) nerv_matrix_cuda_double_##NAME
 #define cudak_(NAME) cudak_double_ ## NAME
 #define NERV_CUBLAS_(NAME) cublasD##NAME
-#define MATRIX_CUMATRIX_HOST_TNAME nerv_matrix_host_double_tname
 #include "generic/cumatrix.c"
diff --git a/nerv/lib/matrix/cumatrix.h b/nerv/lib/matrix/cumatrix.h
index b47e14b..280035b 100644
--- a/nerv/lib/matrix/cumatrix.h
+++ b/nerv/lib/matrix/cumatrix.h
@@ -2,8 +2,19 @@
 #define NERV_CUMATRIX_H
 #include "matrix.h"
 #include "../common.h"
-void nerv_cumatrix_print_profile();
-void nerv_cumatrix_clear_profile();
-void nerv_cumatrix_init();
-void nerv_cumatrix_select_gpu(int dev, Status *status);
+#include "cuda_helper.h"
+
+typedef struct CuContext {
+    cublasHandle_t cublas_handle;
+    cudaEvent_t profile_start, profile_stop;
+    curandGenerator_t curand_gen;
+    HashMap *profile;
+} CuContext;
+
+void nerv_cuda_context_print_profile(CuContext *context);
+void nerv_cuda_context_clear_profile(CuContext *context);
+void nerv_cuda_context_accu_profile(CuContext *context, const char *name, float delta);
+void nerv_cuda_context_select_gpu(CuContext *context, int dev, Status *status);
+CuContext *nerv_cuda_context_create(Status *status);
+void nerv_cuda_context_destroy(CuContext *contex, Status *status);
 #endif
diff --git a/nerv/lib/matrix/generic/cukernel.cu b/nerv/lib/matrix/generic/cukernel.cu
index 51e3b6a..0e09cfa 100644
--- a/nerv/lib/matrix/generic/cukernel.cu
+++ b/nerv/lib/matrix/generic/cukernel.cu
@@ -445,13 +445,12 @@ extern "C" {
         cudaStreamSynchronize(0);
     }
 
-    extern curandGenerator_t curand_gen;
-    void cudak_(cuda_rand_uniform)(const Matrix *a) {
+    void cudak_(cuda_rand_uniform)(const Matrix *a, CuContext *context) {
         #ifdef MATRIX_USE_FLOAT
-        curandGenerateUniform(curand_gen, MATRIX_ELEM_PTR(a), a->nrow * a->stride / sizeof(MATRIX_ELEM));
+        curandGenerateUniform(context->curand_gen, MATRIX_ELEM_PTR(a), a->nrow * a->stride / sizeof(MATRIX_ELEM));
         #endif
         #ifdef MATRIX_USE_DOUBLE
-        curandGenerateUniformDouble(curand_gen, MATRIX_ELEM_PTR(a), a->nrow * a->stride / sizeof(MATRIX_ELEM));
+        curandGenerateUniformDouble(context->curand_gen, MATRIX_ELEM_PTR(a), a->nrow * a->stride / sizeof(MATRIX_ELEM));
         #endif
     }
 
diff --git a/nerv/lib/matrix/generic/cumatrix.c b/nerv/lib/matrix/generic/cumatrix.c
index 7b70607..6342d90 100644
--- a/nerv/lib/matrix/generic/cumatrix.c
+++ b/nerv/lib/matrix/generic/cumatrix.c
@@ -1,10 +1,11 @@
 #ifdef NERV_GENERIC_CUMATRIX
 #include "matrix.h"
 #include "elem_type.h"
-#define MATRIX_DATA_FREE(ptr, status) cuda_matrix_(free)(ptr, status)
-#define MATRIX_DATA_ALLOC(dptr, stride, width, height, status) \
-                            cuda_matrix_(alloc)(dptr, stride, width, height, status)
-
+#define MATRIX_DATA_FREE(ptr, context, status) \
+                            cuda_matrix_(free)(ptr, context, status)
+#define MATRIX_DATA_ALLOC(dptr, stride, width, height, context, status) \
+                            cuda_matrix_(alloc)(dptr, stride, width, height, \
+                                                context, status)
 #define NERV_GENERIC_MATRIX
 #define NERV_GENERIC_CUKERNEL
 #include "../../common.h"
@@ -14,12 +15,13 @@
 
 void nerv_matrix_(add)(Matrix *c, const Matrix *a, const Matrix *b,
                             MATRIX_ELEM alpha, MATRIX_ELEM beta,
+                            CuContext *context,
                             Status *status) {
     CHECK_SAME_DIMENSION(a, b, status);
     CHECK_SAME_DIMENSION(a, c, status);
     PROFILE_START
     CUBLAS_SAFE_SYNC_CALL(
-            NERV_CUBLAS_(geam)(cublas_handle, CUBLAS_OP_N, CUBLAS_OP_N,
+            NERV_CUBLAS_(geam)(context->cublas_handle, CUBLAS_OP_N, CUBLAS_OP_N,
                 a->ncol, a->nrow,
                 &alpha,
                 MATRIX_ELEM_PTR(a), a->stride / sizeof(MATRIX_ELEM),
@@ -33,7 +35,8 @@ void nerv_matrix_(add)(Matrix *c, const Matrix *a, const Matrix *b,
 
 void nerv_matrix_(mul)(Matrix *c, const Matrix *a, const Matrix *b,
                             MATRIX_ELEM alpha, MATRIX_ELEM beta,
-                            int ta, int tb, Status *status) {
+                            int ta, int tb,
+                            CuContext *context, Status *status) {
 #define SWAP(a, b) \
     do { int t = (a); (a) = (b); (b) = t; } while (0)
 
@@ -46,7 +49,7 @@ void nerv_matrix_(mul)(Matrix *c, const Matrix *a, const Matrix *b,
     /* Because matrix in Nerv is row-major, here b comes first */
     PROFILE_START
     CUBLAS_SAFE_SYNC_CALL(
-            NERV_CUBLAS_(gemm)(cublas_handle, tb, ta,
+            NERV_CUBLAS_(gemm)(context->cublas_handle, tb, ta,
                 bn, am, bm,
                 &alpha,
                 MATRIX_ELEM_PTR(b), b->stride / sizeof(MATRIX_ELEM),
@@ -58,7 +61,8 @@ void nerv_matrix_(mul)(Matrix *c, const Matrix *a, const Matrix *b,
     NERV_SET_STATUS(status, NERV_NORMAL, 0);
 }
 
-void nerv_matrix_(sigmoid)(Matrix *a, const Matrix *b, Status *status) {
+void nerv_matrix_(sigmoid)(Matrix *a, const Matrix *b,
+                            CuContext *context, Status *status) {
     CHECK_SAME_DIMENSION(a, b, status);
     PROFILE_START
     cudak_(cuda_sigmoid)(b, a);
@@ -67,7 +71,8 @@ void nerv_matrix_(sigmoid)(Matrix *a, const Matrix *b, Status *status) {
 }
 
 void nerv_matrix_(sigmoid_grad)(Matrix *nerr, const Matrix *err,
-                                const Matrix *output, Status *status) {
+                                const Matrix *output,
+                                CuContext *context, Status *status) {
     CHECK_SAME_DIMENSION(nerr, err, status);
     CHECK_SAME_DIMENSION(nerr, output, status);
     PROFILE_START
@@ -76,14 +81,16 @@ void nerv_matrix_(sigmoid_grad)(Matrix *nerr, const Matrix *err,
     NERV_SET_STATUS(status, NERV_NORMAL, 0);
 }
 
-void nerv_matrix_(rand_uniform)(Matrix *a, Status *status) {
+void nerv_matrix_(rand_uniform)(Matrix *a, CuContext *context, Status *status) {
     PROFILE_START
-    cudak_(cuda_rand_uniform)(a);
+    cudak_(cuda_rand_uniform)(a, context);
     PROFILE_STOP
     NERV_SET_STATUS(status, NERV_NORMAL, 0);
 }
 
-void nerv_matrix_(thres_mask)(Matrix *a, Matrix *b, double thres, double low, double high, Status *status) {
+void nerv_matrix_(thres_mask)(Matrix *a, Matrix *b, double thres,
+                                double low, double high,
+                                CuContext *context, Status *status) {
     CHECK_SAME_DIMENSION(a, b, status);
     PROFILE_START
     cudak_(cuda_thres_mask)(a, b, thres, low, high);
@@ -91,7 +98,8 @@ void nerv_matrix_(thres_mask)(Matrix *a, Matrix *b, double thres, double low, do
     NERV_SET_STATUS(status, NERV_NORMAL, 0);
 }
 
-void nerv_matrix_(tanh)(Matrix *a, const Matrix *b, Status *status) {
+void nerv_matrix_(tanh)(Matrix *a, const Matrix *b,
+                        CuContext *context, Status *status) {
     CHECK_SAME_DIMENSION(a, b, status);
     PROFILE_START
     cudak_(cuda_tanh)(b, a);
@@ -99,8 +107,8 @@ void nerv_matrix_(tanh)(Matrix *a, const Matrix *b, Status *status) {
     NERV_SET_STATUS(status, NERV_NORMAL, 0);
 }
 
-void nerv_matrix_(tanh_grad)(Matrix *nerr, const Matrix *err,
-                                const Matrix *output, Status *status) {
+void nerv_matrix_(tanh_grad)(Matrix *nerr, const Matrix *err, const Matrix *output,
+                                CuContext *context, Status *status) {
     CHECK_SAME_DIMENSION(nerr, err, status);
     CHECK_SAME_DIMENSION(nerr, output, status);
     PROFILE_START
@@ -109,24 +117,25 @@ void nerv_matrix_(tanh_grad)(Matrix *nerr, const Matrix *err,
     NERV_SET_STATUS(status, NERV_NORMAL, 0);
 }
 
-Matrix *nerv_matrix_(softmax)(Matrix *b, const Matrix *a, Status *status) {
+Matrix *nerv_matrix_(softmax)(Matrix *b, const Matrix *a,
+                                CuContext *context, Status *status) {
     Matrix *max, *max_idx;
     Matrix *dno;
     CHECK_SAME_DIMENSION_RET(a, b, status);
-    max = nerv_matrix_(create)(a->nrow, 1, status);
+    max = nerv_matrix_(create)(a->nrow, 1, context, status);
     if (status->err_code != NERV_NORMAL)
         return NULL;
-    max_idx = nerv_matrix_(create)(a->nrow, 1, status);
+    max_idx = nerv_matrix_(create)(a->nrow, 1, context, status);
     if (status->err_code != NERV_NORMAL)
     {
-        nerv_matrix_(destroy)(max, status);
+        nerv_matrix_(destroy)(max, context, status);
         return NULL;
     }
-    dno = nerv_matrix_(create)(a->nrow, 1, status);
+    dno = nerv_matrix_(create)(a->nrow, 1, context, status);
     if (status->err_code != NERV_NORMAL)
     {   /* FIXME: destroy may also fail? */
-        nerv_matrix_(destroy)(max, status);
-        nerv_matrix_(destroy)(max_idx, status);
+        nerv_matrix_(destroy)(max, context, status);
+        nerv_matrix_(destroy)(max_idx, context, status);
         return NULL;
     }
     PROFILE_START
@@ -134,14 +143,14 @@ Matrix *nerv_matrix_(softmax)(Matrix *b, const Matrix *a, Status *status) {
     cudak_(cuda_softmax_denominator)(a, max, dno);
     cudak_(cuda_softmax_final)(a, max, dno, b);
     PROFILE_STOP
-    nerv_matrix_(destroy)(max, status);
-    nerv_matrix_(destroy)(dno, status);
+    nerv_matrix_(destroy)(max, context, status);
+    nerv_matrix_(destroy)(dno, context, status);
     NERV_SET_STATUS(status, NERV_NORMAL, 0);
     return max_idx;
 }
 
-Matrix *nerv_matrix_(rowsum)(Matrix *a, Status *status) {
-    Matrix *b = nerv_matrix_(create)(a->nrow, 1, status);
+Matrix *nerv_matrix_(rowsum)(Matrix *a, CuContext *context, Status *status) {
+    Matrix *b = nerv_matrix_(create)(a->nrow, 1, context, status);
     if (status->err_code != NERV_NORMAL)
         return NULL;
     PROFILE_START
@@ -151,8 +160,8 @@ Matrix *nerv_matrix_(rowsum)(Matrix *a, Status *status) {
     return b;
 }
 
-Matrix *nerv_matrix_(colsum)(Matrix *a, Status *status) {
-    Matrix *b = nerv_matrix_(create)(1, a->ncol, status);
+Matrix *nerv_matrix_(colsum)(Matrix *a, CuContext *context, Status *status) {
+    Matrix *b = nerv_matrix_(create)(1, a->ncol, context, status);
     if (status->err_code != NERV_NORMAL)
         return NULL;
     PROFILE_START
@@ -163,8 +172,8 @@ Matrix *nerv_matrix_(colsum)(Matrix *a, Status *status) {
 }
 
 Matrix *nerv_matrix_(colsame)(Matrix *a, const Matrix *ref,
-                                Status *status) {
-    Matrix *b = nerv_matrix_(create)(1, a->ncol, status);
+                                CuContext *context, Status *status) {
+    Matrix *b = nerv_matrix_(create)(1, a->ncol, context, status);
     if (status->err_code != NERV_NORMAL)
         return NULL;
     CHECK_SAME_DIMENSION_RET(a, ref, status);
@@ -175,8 +184,8 @@ Matrix *nerv_matrix_(colsame)(Matrix *a, const Matrix *ref,
     return b;
 }
 
-Matrix *nerv_matrix_(rowmax)(Matrix *a, Status *status) {
-    Matrix *b = nerv_matrix_(create)(a->nrow, 1, status);
+Matrix *nerv_matrix_(rowmax)(Matrix *a, CuContext *context, Status *status) {
+    Matrix *b = nerv_matrix_(create)(a->nrow, 1, context, status);
     if (status->err_code != NERV_NORMAL)
         return NULL;
     PROFILE_START
@@ -187,15 +196,15 @@ Matrix *nerv_matrix_(rowmax)(Matrix *a, Status *status) {
 }
 
 void nerv_matrix_(rowmax_idx)(Matrix *a, Matrix **b, Matrix **idx,
-                                Status *status) {
-    *b = nerv_matrix_(create)(a->nrow, 1, status);
+                                CuContext *context, Status *status) {
+    *b = nerv_matrix_(create)(a->nrow, 1, context, status);
     if (status->err_code != NERV_NORMAL)
         return;
-    *idx = nerv_matrix_(create)(a->nrow, 1, status);
+    *idx = nerv_matrix_(create)(a->nrow, 1, context, status);
     if (status->err_code != NERV_NORMAL)
     {
         /* FIXME: destroy may also fail? */
-        nerv_matrix_(destroy)(*b, status);
+        nerv_matrix_(destroy)(*b, context, status);
         return;
     }
     PROFILE_START
@@ -205,7 +214,7 @@ void nerv_matrix_(rowmax_idx)(Matrix *a, Matrix **b, Matrix **idx,
 }
 
 void nerv_matrix_(add_row)(Matrix *b, const Matrix *a, double beta,
-                            Status *status) {
+                            CuContext *context, Status *status) {
     if (a->ncol != b->ncol)
         NERV_EXIT_STATUS(status, MAT_MISMATCH_DIM, 0);
     if (a->nrow != 1)
@@ -216,23 +225,25 @@ void nerv_matrix_(add_row)(Matrix *b, const Matrix *a, double beta,
     NERV_SET_STATUS(status, NERV_NORMAL, 0);
 }
 
-void nerv_matrix_(fill)(Matrix *self, double val, Status *status) {
+void nerv_matrix_(fill)(Matrix *self, double val,
+                        CuContext *context, Status *status) {
     PROFILE_START
     cudak_(cuda_fill)(self, val);
     PROFILE_STOP
     NERV_SET_STATUS(status, NERV_NORMAL, 0);
 }
 
-void nerv_matrix_(clip)(Matrix *self, double val_1, double val_2, Status *status) {
+void nerv_matrix_(clip)(Matrix *self, double val1, double val2,
+                        CuContext *context, Status *status) {
     PROFILE_START
-    cudak_(cuda_clip)(self, val_1, val_2);
+    cudak_(cuda_clip)(self, val1, val2);
     PROFILE_STOP
     NERV_SET_STATUS(status, NERV_NORMAL, 0);
 }
 
 void nerv_matrix_(copy_fromd)(Matrix *a, const Matrix *b,
                             int a_begin, int b_begin, int b_end,
-                            Status *status) {
+                            CuContext *context, Status *status) {
     if (!(0 <= b_begin && b_begin < b_end && b_end <= b->nrow &&
             a_begin + b_end - b_begin <= a->nrow))
         NERV_EXIT_STATUS(status, MAT_INVALID_COPY_INTERVAL, 0);
@@ -251,7 +262,7 @@ void nerv_matrix_(copy_fromd)(Matrix *a, const Matrix *b,
 
 void nerv_matrix_(copy_fromh)(Matrix *a, const Matrix *b,
                             int a_begin, int b_begin, int b_end,
-                            Status *status) { 
+                            CuContext *context, Status *status) {
     if (!(0 <= b_begin && b_begin < b_end && b_end <= b->nrow &&
             a_begin + b_end - b_begin <= a->nrow))
         NERV_EXIT_STATUS(status, MAT_INVALID_COPY_INTERVAL, 0);
@@ -270,7 +281,7 @@ void nerv_matrix_(copy_fromh)(Matrix *a, const Matrix *b,
 
 void nerv_matrix_(copy_toh)(Matrix *a, const Matrix *b,
                             int a_begin, int a_end, int b_begin,
-                            Status *status) {
+                            CuContext *context, Status *status) {
     if (!(0 <= a_begin && a_begin < a_end && a_end <= a->nrow &&
             b_begin + a_end - a_begin <= b->nrow))
         NERV_EXIT_STATUS(status, MAT_INVALID_COPY_INTERVAL, 0);
@@ -287,15 +298,15 @@ void nerv_matrix_(copy_toh)(Matrix *a, const Matrix *b,
     NERV_SET_STATUS(status, NERV_NORMAL, 0);
 }
 
-Matrix *nerv_matrix_(trans)(Matrix *a, Status *status) {
+Matrix *nerv_matrix_(trans)(Matrix *a, CuContext *context, Status *status) {
     MATRIX_ELEM alpha = 1, beta = 0;
-    Matrix *b = nerv_matrix_(create)(a->ncol, a->nrow, status);
+    Matrix *b = nerv_matrix_(create)(a->ncol, a->nrow, context, status);
     if (status->err_code != NERV_NORMAL)
         return NULL;
     /* FIXME: possible memory leak when lua error is raised */
     PROFILE_START
     CUBLAS_SAFE_SYNC_CALL_RET(
-            NERV_CUBLAS_(geam)(cublas_handle, CUBLAS_OP_T, CUBLAS_OP_T,
+            NERV_CUBLAS_(geam)(context->cublas_handle, CUBLAS_OP_T, CUBLAS_OP_T,
                 a->nrow, a->ncol,
                 &alpha,
                 MATRIX_ELEM_PTR(a), a->stride / sizeof(MATRIX_ELEM),
@@ -309,7 +320,7 @@ Matrix *nerv_matrix_(trans)(Matrix *a, Status *status) {
 }
 
 void nerv_matrix_(mul_elem)(Matrix *c, const Matrix *a, const Matrix *b,
-                            Status *status) {
+                            CuContext *context, Status *status) {
     CHECK_SAME_DIMENSION(a, b, status);
     CHECK_SAME_DIMENSION(a, c, status);
     PROFILE_START
@@ -318,7 +329,8 @@ void nerv_matrix_(mul_elem)(Matrix *c, const Matrix *a, const Matrix *b,
     NERV_SET_STATUS(status, NERV_NORMAL, 0);
 }
 
-void nerv_matrix_(log_elem)(Matrix *b, const Matrix *a, Status *status) {
+void nerv_matrix_(log_elem)(Matrix *b, const Matrix *a,
+                            CuContext *context, Status *status) {
     CHECK_SAME_DIMENSION(a, b, status);
     PROFILE_START
     cudak_(cuda_log_elem)(a, b);
@@ -326,14 +338,15 @@ void nerv_matrix_(log_elem)(Matrix *b, const Matrix *a, Status *status) {
     NERV_SET_STATUS(status, NERV_NORMAL, 0);
 }
 
-Matrix *ner