Skip to content

Instantly share code, notes, and snippets.

@bayerj

bayerj/gist:1012340

Created June 7, 2011 14:18
Show Gist options

  • Select an option

    Learn more about clone URLs
  • Save bayerj/1012340 to your computer and use it in GitHub Desktop.

Select an option

Learn more about clone URLs
Save bayerj/1012340 to your computer and use it in GitHub Desktop.
Download ZIP
theano run warning
Raw
gistfile1.txt
mod.cu: In function ‘PyObject* CudaNdarray_ptr_int_size(PyObject*, PyObject*)’:
mod.cu:1928: warning: comparison between ‘enum cudaError’ and ‘enum cublasStatus_t’
mod.cu: In function ‘int CudaNdarray_gemm(float, const CudaNdarray*, const CudaNdarray*, float, CudaNdarray*)’:
mod.cu:2753: warning: comparison between ‘enum cublasStatus_t’ and ‘enum cudaError_t’
mod.cu: In function ‘int CudaNdarray_sger(float, CudaNdarray*, CudaNdarray*, CudaNdarray*)’:
mod.cu:2795: warning: comparison between ‘enum cublasStatus_t’ and ‘enum cudaError_t’
i686-apple-darwin10-gcc-4.2.1: python: No such file or directory
i686-apple-darwin10-gcc-4.2.1: warning: '-x c++' after last input file has no effect
i686-apple-darwin10-gcc-4.2.1: no input files
===============================
1 #include <Python.h>
2 #include <structmember.h>
3
4 #include <numpy/arrayobject.h>
5 #include <iostream>
6
7 #include "cuda_ndarray.cuh"
8
9 //If true, when there is a gpu malloc or free error, we print the size of allocated memory on the device.
10 #define COMPUTE_GPU_MEM_USED 0
11
12 //If true, we fill with NAN allocated device memory.
13 #define ALLOC_MEMSET 0
14
15 /////////////////////////
16 // Alloc and Free
17 /////////////////////////
18
19 static int g_gpu_context_active = 0;
20
21
22 PyObject *
23 CudaNdarray_Dimshuffle(PyObject* _unused, PyObject* args);
24
25 /**
26 *
27 * In the test program I'm using, the _outstanding_mallocs decreases with every call.
28 * This suggests there are more free() calls being made than alloc(), but I can't figure out why.
29 *
30 */
31 int _outstanding_mallocs[] = {0,0};
32 #if COMPUTE_GPU_MEM_USED
33 int _allocated_size = 0;
34 const int TABLE_SIZE = 10000;
35 struct table_struct{
36 void* ptr;
37 int size;
38 };
39 table_struct _alloc_size_table[TABLE_SIZE];
40 #endif
41 void * device_malloc(size_t size)
42 {
43 void * rval=NULL;
44 cudaError_t err = cudaMalloc(&rval, size);
45 if (cudaSuccess != err)
46 {
47 #if COMPUTE_GPU_MEM_USED
48 fprintf(stderr, "Error allocating %li bytes of device memory (%s). %d already allocated\n", (long)size, cudaGetErrorString(err),_allocated_size);
49 #else
50 fprintf(stderr, "Error allocating %li bytes of device memory (%s).\n", (long)size, cudaGetErrorString(err));
51 #endif
52 PyErr_Format(PyExc_MemoryError, "Error allocating %li bytes of device memory (%s).", (long)size, cudaGetErrorString(err));
53 return NULL;
54 }
55 _outstanding_mallocs[0] += (rval != NULL);
56 #if COMPUTE_GPU_MEM_USED
57 for(int i=0;i<TABLE_SIZE;i++){
58 if(NULL==_alloc_size_table[i].ptr){
59 _alloc_size_table[i].ptr=rval;
60 _alloc_size_table[i].size=size;
61 break;
62 }
63 }
64 _allocated_size += size;
65 #endif
66 //fprintf(stderr, "allocated %li bytes of device memory (%s). %d already allocated, ptr: %p\n", (long)size, cudaGetErrorString(err),_allocated_size,rval);
67
68 if(ALLOC_MEMSET){
69 //We init them to nan to make sure we catch more debug case.
70 cudaMemset(rval, 0xFF, size);
71 //printf("MEMSET\n");
72 }
73 return rval;
74 }
75 int device_free(void *ptr)
76 {
77 // if there is no gpu context, the call to cudaFree will fail; skip it entirely
78 if(!g_gpu_context_active) {
79 return 0;
80 }
81 cudaError_t err = cudaFree(ptr);
82 if (cudaSuccess != err)
83 {
84 #if COMPUTE_GPU_MEM_USED
85 fprintf(stderr, "Error freeing device pointer %p (%s).%d byte already allocated\n", ptr, cudaGetErrorString(err), _allocated_size);
86 #else
87 fprintf(stderr, "Error freeing device pointer %p (%s).\n", ptr, cudaGetErrorString(err));
88 #endif
89 PyErr_Format(PyExc_MemoryError, "error freeing device pointer %p (%s)", ptr, cudaGetErrorString(err));
90 return -1;
91 }
92 _outstanding_mallocs[0] -= (ptr != NULL);
93 #if COMPUTE_GPU_MEM_USED
94 int i=0;
95 size_t total_freed = 0;
96 for(;i<TABLE_SIZE;i++)
97 if(_alloc_size_table[i].ptr==ptr){
98 _allocated_size -= _alloc_size_table[i].size;
99 total_freed += _alloc_size_table[i].size;
100 _alloc_size_table[i].ptr=0;
101 _alloc_size_table[i].size=0;
102
103 break;
104 }
105 if(i==TABLE_SIZE)
106 printf("Unallocated unknow size!\n");
107 //fprintf(stderr, "freed %li bytes of device memory (%s). %d already allocated, ptr=%p\n", (long)total_freed, cudaGetErrorString(err),_allocated_size,ptr);
108 #endif
109 return 0;
110 }
111 static PyObject *
112 outstanding_mallocs(PyObject* self, PyObject * args)
113 {
114 return PyInt_FromLong(_outstanding_mallocs[0]);
115 }
116
117 /////////////////////////
118 // Static helper methods
119 /////////////////////////
120
121 static void
122 CudaNdarray_null_init(CudaNdarray*self)
123 {
124 self->base = NULL;
125 self->nd = -1;
126 self->host_structure = NULL;
127 self->data_allocated = 0;
128 self->dev_structure_fresh = 1;
129 self->dev_structure = NULL;
130 self->devdata = NULL;
131 }
132
133 static int
134 CudaNdarray_uninit(CudaNdarray*self)
135 {
136 int rval = 0;
137 if (self->data_allocated) {
138 assert(self->devdata);
139 if (device_free(self->devdata))
140 {
141 fprintf(stderr,
142 "!!!! error freeing device memory %p (self=%p)\n",
143 self->devdata, self);
144 rval = -1;
145 }
146 self->devdata = NULL;
147 self->data_allocated = 0;
148 }
149 if (self->dev_structure)
150 {
151 if (device_free(self->dev_structure))
152 {
153 fprintf(stderr,
154 "!!!! error freeing dev_structure memory %p (self=%p)\n",
155 self->dev_structure, self);
156 rval = -1;
157 }
158 self->dev_structure = NULL;
159 }
160 if (self->host_structure)
161 {
162 free(self->host_structure);
163 self->host_structure = NULL;
164 }
165 self->nd = -1;
166 Py_XDECREF(self->base);
167 self->base = NULL;
168 return rval;
169 }
170
171
172 //make the rightmost coords change fastest
173 //TODO: why does a downward for-loop not work????
174 //TODO: use the log2_dims and driver code to remove / and %
175 //TODO: skip the last division (when d == 0)
176 #define decl_k_elemwise_unary_rowmajor(name, F) \
177 __global__ void name (unsigned int numEls, \
178 unsigned int nd, \
179 const int * dim, \
180 const float * a_data, const int * a_str, \
181 float * z_data, const int * z_str) \
182 { \
183 const unsigned int idx = blockIdx.x * blockDim.x + threadIdx.x; \
184 const unsigned int numThreads = blockDim.x * gridDim.x; \
185 \
186 for (unsigned int i = idx; i < numEls; i += numThreads) \
187 { \
188 unsigned int ii = i; \
189 const float * a_i = a_data; \
190 float * z_i = z_data; \
191 for (unsigned int _d = 0; _d < nd; ++_d) \
192 { \
193 unsigned int d = nd - _d-1; \
194 /* i_d used to be unsigned, but their is a bug in nvcc 3.0. making it signed fix the bug.*/\
195 int i_d = ii % dim[d]; /* i_d is our position in the d'th dimension */ \
196 ii = ii / dim[d]; \
197 a_i += i_d * a_str[d]; /* increment our a and z pointers by i_d elements */ \
198 z_i += i_d * z_str[d]; \
199 } \
200 z_i[0] = F(a_i[0]); \
201 } \
202 }
203
204 template<typename T> __device__ T unary_copy(T a) { return a; }
205 decl_k_elemwise_unary_rowmajor(k_elemwise_unary_rowmajor_copy, unary_copy<float>)
206
207 template<typename T> __device__ T unary_exp(T a) { return exp(a); }
208 decl_k_elemwise_unary_rowmajor(k_elemwise_unary_rowmajor_exp, unary_exp<float>)
209
210 /////////////////////////////
211 // Satisfying reqs to be Type
212 /////////////////////////////
213
214 //DON'T use directly(if their is other CudaNdarray that point to it, it will cause problem)! use Py_DECREF() instead
215 static void
216 CudaNdarray_dealloc(CudaNdarray* self)
217 {
218 if (0) std::cerr << "CudaNdarray dealloc " << self << " " << self->devdata << '\n';
219 if(self->ob_refcnt>1)
220 printf("WARNING:CudaNdarray_dealloc called when their is still active reference to it.\n");
221 CudaNdarray_uninit(self);
222 self->ob_type->tp_free((PyObject*)self);
223 --_outstanding_mallocs[1];
224 if (0)
225 {
226 fprintf(stderr, "device_malloc_counts: (device) %i (obj) %i\n",
227 _outstanding_mallocs[0],
228 _outstanding_mallocs[1]);
229 }
230 }
231
232 static PyObject *
233 CudaNdarray_new(PyTypeObject *type, PyObject *args, PyObject *kwds)
234 {
235 CudaNdarray *self;
236
237 self = (CudaNdarray *)type->tp_alloc(type, 0);
238 if (self != NULL)
239 {
240 CudaNdarray_null_init(self);
241 ++_outstanding_mallocs[1];
242 }
243 return (PyObject *)self;
244 }
245 static int
246 CudaNdarray_init(CudaNdarray *self, PyObject *args, PyObject *kwds)
247 {
248 PyObject *arr=NULL;
249
250 if (! PyArg_ParseTuple(args, "O", &arr))
251 return -1;
252 if (! PyArray_Check(arr))
253 {
254 PyErr_SetString(PyExc_TypeError, "PyArray arg required");
255 return -1;
256 }
257 int rval = CudaNdarray_CopyFromArray(self, (PyArrayObject*)arr);
258 return rval;
259 }
260 static PyMemberDef CudaNdarray_members[] =
261 {
262 /*
263 {"first", T_OBJECT_EX, offsetof(CudaNdarray, first), 0,
264 "first name"},
265 {"last", T_OBJECT_EX, offsetof(CudaNdarray, last), 0,
266 "last name"},
267 {"number", T_INT, offsetof(CudaNdarray, number), 0,
268 "noddy number"},
269 */
270 {NULL} /* Sentinel */
271 };
272
273 PyObject * CudaNdarray_CreateArrayObj(CudaNdarray * self)
274 {
275 int verbose = 0;
276 if(self->nd>=0 && CudaNdarray_SIZE(self)==0){
277 npy_intp * npydims = (npy_intp*)malloc(self->nd * sizeof(npy_intp));
278 assert (npydims);
279 for (int i = 0; i < self->nd; ++i) npydims[i] = (npy_intp)(CudaNdarray_HOST_DIMS(self)[i]);
280 PyObject * rval = PyArray_SimpleNew(self->nd, npydims, REAL_TYPENUM);
281 free(npydims);
282 if (!rval){
283 return NULL;
284 }
285 assert (PyArray_ITEMSIZE(rval) == sizeof(real));
286 return rval;
287 }
288 if ((self->nd < 0) || (self->devdata == 0))
289 {
290 PyErr_SetString(PyExc_ValueError, "can't copy from un-initialized CudaNdarray");
291 return NULL;
292 }
293 CudaNdarray * contiguous_self = NULL;
294 if (CudaNdarray_is_c_contiguous(self))
295 {
296 contiguous_self = self;
297 Py_INCREF(contiguous_self);
298 if (verbose) std::cerr << "CreateArrayObj already contiguous" << contiguous_self << '\n';
299 }
300 else
301 {
302 contiguous_self = (CudaNdarray*)CudaNdarray_Copy(self);
303 if (verbose) std::cerr << "CreateArrayObj created contiguous" << contiguous_self << '\n';
304 }
305 if (!contiguous_self)
306 {
307 return NULL;
308 }
309
310 npy_intp * npydims = (npy_intp*)malloc(self->nd * sizeof(npy_intp));
311 assert (npydims);
312 for (int i = 0; i < self->nd; ++i)
313 npydims[i] = (npy_intp)(CudaNdarray_HOST_DIMS(self)[i]);
314 PyObject * rval = PyArray_SimpleNew(self->nd, npydims, REAL_TYPENUM);
315 free(npydims);
316 if (!rval)
317 {
318 Py_DECREF(contiguous_self);
319 return NULL;
320 }
321
322 assert (PyArray_ITEMSIZE(rval) == sizeof(real));
323
324 cublasGetVector(PyArray_SIZE(rval), sizeof(real),
325 contiguous_self->devdata, 1,
326 PyArray_DATA(rval), 1);
327 CNDA_THREAD_SYNC;
328
329 if (CUBLAS_STATUS_SUCCESS != cublasGetError())
330 {
331 PyErr_SetString(PyExc_RuntimeError, "error copying data to host");
332 Py_DECREF(rval);
333 rval = NULL;
334 }
335
336 Py_DECREF(contiguous_self);
337 return rval;
338 }
339
340 // TODO-- we have two functions here, ZEROS and Zeros.
341 // ZEROS is meant to be called just from C code (you don't need to pass it PyObject * s)
342 // but this naming is very weird, makes it look like a macro
343 // we should figure out the correct convention and change to that
344 PyObject* CudaNdarray_ZEROS(int n, int * dims)
345 {
346
347 int total_elements = 1;
348 for(int i=0;i<n;i++)
349 total_elements*=dims[i];
350
351 // total_elements now contains the size of the array, in reals
352 int total_size = total_elements * sizeof(real);
353
354 CudaNdarray* rval = (CudaNdarray*)CudaNdarray_New();
355 if (!rval)
356 {
357 PyErr_SetString(PyExc_RuntimeError, "CudaNdarray_ZEROS: call to New failed");
358 return NULL;
359 }
360
361 if (CudaNdarray_alloc_contiguous(rval, n, dims))
362 {
363 PyErr_SetString(PyExc_RuntimeError, "CudaNdarray_ZEROS: allocation failed.");
364 Py_DECREF(rval);
365 return NULL;
366 }
367
368 // Fill with zeros
369 //fprintf(stdout, "Sizeof: %d\n", total_size);
370 if (cudaSuccess != cudaMemset(rval->devdata, 0, total_size))
371 {
372 PyErr_Format(PyExc_MemoryError, "CudaNdarray_ZEROS: Error memsetting %d bytes of device memory.", total_size);
373 Py_DECREF(rval);
374 return NULL;
375 }
376
377 if (cnda_copy_structure_to_device(rval))
378 {
379 PyErr_SetString(PyExc_RuntimeError, "CudaNdarray_ZEROS: syncing structure to device failed");
380 Py_DECREF(rval);
381 return NULL;
382 }
383 return (PyObject*) rval;
384 }
385
386 // declared as a static method (hence "dummy" is not used)
387 // Based on _Copy and _dimshuffle
388 PyObject* CudaNdarray_Zeros(PyObject* dummy, PyObject* shape)
389 {
390 if(!PySequence_Check(shape))
391 {
392 PyErr_SetString(PyExc_TypeError, "shape argument must be a sequence");
393 return NULL;
394 }
395
396 int shplen = PySequence_Length(shape);
397
398 if (shplen == 0)
399 {
400 PyErr_SetString(PyExc_ValueError,
401 "CudaNdarray_Zeros: empty shape not allowed");
402 return NULL;
403 }
404
405 int* newdims = (int *)malloc(sizeof(int) * shplen);
406
407 if (!newdims)
408 {
409 PyErr_SetString(PyExc_MemoryError,
410 "CudaNdarray_Zeros: Failed to allocate temporary space");
411 return NULL;
412 }
413
414 // start from the end to compute strides
415 for (int i = shplen-1; i >= 0; --i)
416 {
417 PyObject* shp_el_obj = PySequence_GetItem(shape, i);
418 if(shp_el_obj == NULL)
419 {
420 // shouldn't happen since we checked length before...
421 PyErr_SetString(PyExc_RuntimeError, "CudaNdarray_Zeros: Index out of bound in sequence");
422 free(newdims);
423 return NULL;
424 }
425
426 int shp_el = PyInt_AsLong(shp_el_obj);
427 Py_DECREF(shp_el_obj);
428
429 if (shp_el <= 0)
430 {
431 PyErr_SetString(PyExc_ValueError, "CudaNdarray_Zeros: shape must not contain 0 (or negative value) for size of a dimension");
432 free(newdims);
433 return NULL;
434 }
435
436 newdims[i] = shp_el;
437 }
438
439 PyObject* rval = CudaNdarray_ZEROS(shplen,newdims);
440
441 free(newdims);
442
443 return (PyObject*)rval;
444 }
445
446
447
448
449
450 PyObject * CudaNdarray_Copy(CudaNdarray * self)
451 {
452 PyObject * rval = CudaNdarray_New();
453 if ((!rval) || (-1 == self->nd))
454 {
455 return rval;
456 }
457 if (CudaNdarray_alloc_contiguous((CudaNdarray*)rval, self->nd, CudaNdarray_HOST_DIMS(self)))
458 {
459 Py_DECREF(rval);
460 return NULL;
461 }
462 if (CudaNdarray_CopyFromCudaNdarray((CudaNdarray*)rval, self))
463 {
464 Py_DECREF(rval);
465 return NULL;
466 }
467 return rval;
468 }
469 PyObject * CudaNdarray_DeepCopy(CudaNdarray * self, PyObject * memo)
470 {
471 assert(PyDict_Check(memo));
472 PyObject * selfkey = PyInt_FromLong((long)self);
473 assert(selfkey);
474 if (PyDict_Contains(memo, selfkey))
475 {
476 PyObject * rval = PyDict_GetItem(memo, selfkey);
477 Py_DECREF(selfkey);
478 Py_XINCREF(rval);
479 return rval;
480 }
481 else
482 {
483 PyObject * rval = CudaNdarray_Copy(self);
484 if (0) std::cerr << "DeepCopy created " << rval << " devdata " << ((CudaNdarray*)rval)->devdata << "\n";
485 if (NULL == rval)
486 {
487 Py_DECREF(selfkey);
488 return NULL;
489 }
490 if (PyDict_SetItem(memo, selfkey, rval))
491 {
492 Py_DECREF(rval);
493 Py_DECREF(selfkey);
494 return NULL;
495 }
496 Py_DECREF(selfkey);
497 return rval;
498 }
499 }
500 PyObject * CudaNdarray_ReduceSum(CudaNdarray * self, PyObject * py_reduce_mask)
501 {
502 if (!PySequence_Check(py_reduce_mask))
503 {
504 PyErr_SetString(PyExc_TypeError, "reduce_mask must be sequence of ints");
505 return NULL;
506 }
507 int len = PySequence_Length(py_reduce_mask);
508 if (len != self->nd)
509 {
510 PyErr_SetString(PyExc_TypeError, "length of reduce_mask must match self->nd");
511 return NULL;
512 }
513 CudaNdarray * self_sum = (CudaNdarray*)CudaNdarray_New();
514 if (!self_sum)
515 {
516 return NULL;
517 }
518 //TODO: allocate a fixed size dimshuffle_pattern_cache on the stack,
519 // and use it if it is big enough.
520 int * dimshuffle_pattern = (int*)malloc(len * 2 * sizeof(int));
521 int * sum_dims = dimshuffle_pattern + len;
522 int n_remaining_dims = 0;
523 if (!dimshuffle_pattern)
524 {
525 Py_DECREF(self_sum);
526 PyErr_SetString(PyExc_MemoryError, "failed to alloc internal storage");
527 return NULL;
528 }
529 for (int i = 0; i < len; ++i)
530 {
531 PyObject *o_i = PySequence_GetItem(py_reduce_mask, i);
532 int o_i_int = PyInt_AsLong(o_i);
533 Py_XDECREF(o_i);
534 if (PyErr_Occurred())
535 {
536 Py_DECREF(self_sum);
537 free(dimshuffle_pattern);
538 return NULL;
539 }
540 if (o_i_int) // this is a dimension over which we are reducing
541 {
542 sum_dims[i] = 1;
543 }
544 else
545 {
546 sum_dims[i] = CudaNdarray_HOST_DIMS(self)[i];
547 dimshuffle_pattern[n_remaining_dims++] = i;
548 }
549 }
550 if (0 || CudaNdarray_alloc_contiguous(self_sum, len, sum_dims)
551 || CudaNdarray_reduce_sum(self_sum, self)
552 || CudaNdarray_dimshuffle(self_sum, n_remaining_dims, dimshuffle_pattern))
553 {
554 Py_DECREF(self_sum);
555 free(dimshuffle_pattern);
556 return NULL;
557 }
558 free(dimshuffle_pattern);
559 return (PyObject*)self_sum;
560 }
561
562 __global__ void k_copy_reshape_rowmajor(unsigned int numEls,
563 unsigned int a_nd, const float * a_data, const int * a_dim, const int * a_str,
564 unsigned int z_nd, float * z_data, const int * z_dim, const int * z_str)
565 {
566 const unsigned int idx = blockIdx.x * blockDim.x + threadIdx.x;
567 const unsigned int numThreads = blockDim.x * gridDim.x;
568
569 for (unsigned int i = idx; i < numEls; i += numThreads)
570 {
571 const float * a_i = a_data;
572 unsigned int a_ii = i;
573 for (unsigned int _d = 0; _d < a_nd; ++_d) //make the rightmost coords change fastest
574 {
575 unsigned int d = a_nd - _d-1;
576 unsigned int a_i_d = a_ii % a_dim[d];
577 a_ii = a_ii / a_dim[d];
578 a_i += a_i_d * a_str[d];
579 }
580 unsigned int z_ii = i;
581 float * z_i = z_data;
582 for (unsigned int _d = 0; _d < z_nd; ++_d) //make the rightmost coords change fastest
583 {
584 unsigned int d = z_nd - _d-1;
585 //i tried to make the for loop count down, but it didn't work!?
586 unsigned int z_i_d = z_ii % z_dim[d];
587 z_i += z_i_d * z_str[d];
588 z_ii = z_ii / z_dim[d];
589 }
590 z_i[0] = a_i[0]; //copy one lousy float!
591 }
592 }
593 PyObject * CudaNdarray_Reshape(CudaNdarray * self, PyObject * shape)
594 {
595 // check shape tuple
596 unsigned int rval_nd;
597 unsigned int * rval_dims;
598 unsigned int rval_size = 1;
599
600 if (PyTuple_Check(shape)){
601 // copy shape to integer array
602 rval_nd = PyTuple_Size(shape);
603 }else if (PyInt_Check(shape)){
604 rval_nd = 1;
605 }else{
606 PyErr_SetString(PyExc_TypeError, "shape must be tuple of integers or an integer");
607 return NULL;
608 }
609 rval_dims = (unsigned int*)malloc(rval_nd * sizeof(int));
610
611 if(PyTuple_Check(shape)){
612 for (int i = 0; i < rval_nd; ++i)
613 {
614 rval_dims[i] = PyInt_AsLong(PyTuple_GetItem(shape, i)); //GetItem returns borrowed reference
615 if (PyErr_Occurred()) //error in AsLong
616 {
617 free(rval_dims);
618 return NULL;
619 }
620 if(rval_dims[i]<=0){
621 PyErr_Format(PyExc_ValueError, "Reshape has invalid dimension %i (must be >0)",rval_dims[i]);
622 free(rval_dims);
623 return NULL;
624 }
625 rval_size = rval_size * rval_dims[i];
626 }
627 }else{
628 rval_size = PyInt_AsLong(shape);
629 rval_dims[0] = rval_size;
630 }
631 // calculate new size, assert same as old size
632 if (rval_size != CudaNdarray_SIZE(self))
633 {
634 PyErr_Format(PyExc_ValueError, "size must remain unchanged, changed from %i to %i", CudaNdarray_SIZE(self), rval_size);
635 free(rval_dims);
636 return NULL;
637 }
638 if (rval_size==0)
639 {
640 PyObject * rval = CudaNdarray_NewDims(rval_nd, rval_dims);
641 free(rval_dims);
642 return rval;
643 }
644
645 if(CudaNdarray_is_c_contiguous(self))
646 {
647 //return a view, not a copy
648 CudaNdarray * rval = (CudaNdarray * )CudaNdarray_New(rval_nd);
649
650 if (!rval || 0 != rval->data_allocated
651 ||CudaNdarray_set_device_data(rval, CudaNdarray_DEV_DATA(self), self))
652 {
653 Py_XDECREF(rval);
654 free(rval_dims);
655 return NULL;
656 }
657 //set dim and stride
658 int size = 1;
659 for (int i = rval_nd-1; i >= 0; --i)
660 {
661 CudaNdarray_set_stride(rval, i, (rval_dims[i] == 1) ? 0 : size);
662 CudaNdarray_set_dim(rval, i, rval_dims[i]);
663 size = size * rval_dims[i];
664 }
665 free(rval_dims);
666 return (PyObject*)rval;
667 }
668
669 // allocate new space (TODO: test to see if we can re-use old one)
670 CudaNdarray * rval = (CudaNdarray * )CudaNdarray_New();
671 if (!rval || CudaNdarray_alloc_contiguous(rval, rval_nd, rval_dims)){
672 Py_XDECREF(rval);
673 free(rval_dims);
674 return NULL;
675 }
676
677 // call worker routine
678 unsigned int threads_per_block = std::min(rval_size, (unsigned int)NUM_VECTOR_OP_THREADS_PER_BLOCK);
679 unsigned int n_blocks = std::min(ceil_intdiv(rval_size,threads_per_block), (unsigned int)NUM_VECTOR_OP_BLOCKS);
680 k_copy_reshape_rowmajor<<<n_blocks,threads_per_block>>>(
681 rval_size,
682 self->nd,
683 CudaNdarray_DEV_DATA(self), CudaNdarray_DEV_DIMS(self), CudaNdarray_DEV_STRIDES(self),
684 rval->nd,
685 CudaNdarray_DEV_DATA(rval), CudaNdarray_DEV_DIMS(rval), CudaNdarray_DEV_STRIDES(rval));
686
687 CNDA_THREAD_SYNC;
688 cudaError_t err = cudaGetLastError();
689 if( cudaSuccess != err)
690 {
691 Py_DECREF(rval);
692 PyErr_Format(PyExc_RuntimeError, "Cuda error: %s: %s.\n", "k_copy_reshape_rowmajor", cudaGetErrorString(err));
693 free(rval_dims);
694 return NULL;
695 }
696 free(rval_dims);
697 return (PyObject*)rval;
698 }
699 PyObject * CudaNdarray_View(CudaNdarray * self)
700 {
701 CudaNdarray * rval = (CudaNdarray*)CudaNdarray_New(self->nd);
702 if (!rval || CudaNdarray_set_device_data(rval, CudaNdarray_DEV_DATA(self), self))
703 {
704 Py_XDECREF(rval);
705 rval = NULL;
706 }
707 else
708 {
709 for (int i = 0; i < self->nd; ++i)
710 {
711 CudaNdarray_set_dim(rval, i, CudaNdarray_HOST_DIMS(self)[i]);
712 CudaNdarray_set_stride(rval, i, CudaNdarray_HOST_STRIDES(self)[i]);
713 }
714 }
715 return (PyObject*)rval;
716 }
717 PyObject * CudaNdarray_SetStride(CudaNdarray * self, PyObject *args)
718 {
719 int pos, stride;
720 if (! PyArg_ParseTuple(args, "ii", &pos, &stride))
721 return NULL;
722 if ((pos < 0) || (pos >= self->nd))
723 {
724 PyErr_Format(PyExc_ValueError, "position argument out of legal range [0, %i)", self->nd);
725 return NULL;
726 }
727 CudaNdarray_set_stride(self, pos, stride);
728 if (cnda_copy_structure_to_device(self))
729 {
730 return NULL;
731 }
732 Py_INCREF(Py_None);
733 return Py_None;
734 }
735 PyObject * CudaNdarray_SetShapeI(CudaNdarray * self, PyObject *args)
736 {
737 int pos, dim;
738 if (! PyArg_ParseTuple(args, "ii", &pos, &dim))
739 return NULL;
740 if ((pos < 0) || (pos >= self->nd))
741 {
742 PyErr_Format(PyExc_ValueError, "position argument out of legal range [0, %i)", self->nd);
743 return NULL;
744 }
745 CudaNdarray_set_dim(self, pos, dim);
746 if (cnda_copy_structure_to_device(self))
747 {
748 return NULL;
749 }
750 Py_INCREF(Py_None);
751 return Py_None;
752 }
753
754 static PyObject *
755 CudaNdarray_exp(CudaNdarray* self)
756 {
757 CudaNdarray * rval = (CudaNdarray *)CudaNdarray_New();
758 if ((NULL == rval) || CudaNdarray_alloc_contiguous(rval, self->nd, CudaNdarray_HOST_DIMS(self)))
759 {
760 Py_XDECREF(rval);
761 return NULL;
762 }
763 unsigned int size = 1;
764 for (int i = 0; i < self->nd; i++)
765 {
766 size *= (unsigned int) CudaNdarray_HOST_DIMS(self)[i];
767 }
768 unsigned int threads_per_block = std::min(size, (unsigned int)NUM_VECTOR_OP_THREADS_PER_BLOCK);
769 unsigned int n_blocks = std::min(ceil_intdiv(size,threads_per_block), (unsigned int)NUM_VECTOR_OP_BLOCKS);
770 k_elemwise_unary_rowmajor_exp<<<n_blocks,threads_per_block>>>(size, self->nd, CudaNdarray_DEV_DIMS(self),
771 CudaNdarray_DEV_DATA(self), CudaNdarray_DEV_STRIDES(self),
772 CudaNdarray_DEV_DATA(rval), CudaNdarray_DEV_STRIDES(rval));
773
774 //TODO: don't do this right away, do it when we need the result
775 CNDA_THREAD_SYNC;
776 cudaError_t err = cudaGetLastError();
777 if( cudaSuccess != err)
778 {
779 Py_DECREF(rval);
780 PyErr_Format(PyExc_RuntimeError, "Cuda error: %s: %s.\n", "kExp", cudaGetErrorString(err));
781 return NULL;
782 }
783
784 return (PyObject*)rval;
785 }
786 static PyMethodDef CudaNdarray_methods[] =
787 {
788 {"__array__",
789 (PyCFunction)CudaNdarray_CreateArrayObj, METH_NOARGS,
790 "Copy from the device to a numpy ndarray"},
791 {"__copy__",
792 (PyCFunction)CudaNdarray_View, METH_NOARGS,
793 "Create a shallow copy of this object. used by module copy"},
794 {"__deepcopy__",
795 (PyCFunction)CudaNdarray_DeepCopy, METH_O,
796 "Create a copy of this object"},
797 {"zeros",
798 (PyCFunction)CudaNdarray_Zeros, METH_STATIC,
799 "Create a new CudaNdarray with specified shape, filled with zeros."},
800 {"copy",
801 (PyCFunction)CudaNdarray_Copy, METH_NOARGS,
802 "Create a copy of this object"},
803 {"reduce_sum",
804 (PyCFunction)CudaNdarray_ReduceSum, METH_O,
805 "Reduce over the given dimensions by summation"},
806 {"exp",
807 (PyCFunction)CudaNdarray_exp, METH_NOARGS,
808 "Return the exponential of all elements"},
809 {"reshape",
810 (PyCFunction)CudaNdarray_Reshape, METH_O,
811 "Return a reshaped view (or copy) of this ndarray\n\
812 The required argument is a tuple of integers specifying the shape of the new ndarray."},
813 {"view",
814 (PyCFunction)CudaNdarray_View, METH_NOARGS,
815 "Return an alias of this ndarray"},
816 {"_set_stride",
817 (PyCFunction)CudaNdarray_SetStride, METH_VARARGS,
818 "For integer arguments (i, s), set the 'i'th stride to 's'"},
819 {"_set_shape_i",
820 (PyCFunction)CudaNdarray_SetShapeI, METH_VARARGS,
821 "For integer arguments (i, s), set the 'i'th shape to 's'"},
822 {NULL, NULL, NULL, NULL} /* Sentinel */
823 };
824
825
826 ////////////////////
827 // Number protocol
828 ////////////////////
829
830 __global__ void kAdd_contiguous(float* a, float* b, float* dest, unsigned int numEls) {
831 const unsigned int idx = blockIdx.x * blockDim.x + threadIdx.x;
832 const unsigned int numThreads = blockDim.x * gridDim.x;
833
834 for (unsigned int i = idx; i < numEls; i += numThreads) {
835 dest[i] = a[i] + b[i];
836 }
837 }
838
839 // Will be called by __add__ in Python
840 static PyObject *
841 CudaNdarray_add(PyObject* py_self, PyObject * py_other)
842 {
843 if (! CudaNdarray_Check(py_self)) {
844 PyErr_SetString(PyExc_TypeError, "need a CudaNdarray on left");
845 return NULL;
846 }
847 if (! CudaNdarray_Check(py_other)) {
848 PyErr_SetString(PyExc_TypeError, "need a CudaNdarray on right");
849 return NULL;
850 }
851 CudaNdarray * self = (CudaNdarray *)py_self;
852 CudaNdarray * other = (CudaNdarray *)py_other;
853 if(!CudaNdarray_is_c_contiguous(self) || !CudaNdarray_is_c_contiguous(other)){
854 PyErr_SetString(PyExc_TypeError, "We have implementet only the c_contiguous version for now.");
855 return NULL;
856 }
857
858 //standard elemwise size checks
859 if (self->nd != other->nd)
860 {
861 PyErr_SetString(PyExc_TypeError, "CudaNdarray_add: need same number of dims");
862 return NULL;
863 }
864 //standard elemwise dim checks
865 unsigned int size = 1;
866 for (int i = 0; i< self->nd; ++i)
867 {
868 if (CudaNdarray_HOST_DIMS(self)[i] != CudaNdarray_HOST_DIMS(other)[i])
869 {
870 PyErr_SetString(PyExc_TypeError, "need same dimensions");
871 return NULL;
872 }
873 size *= (unsigned int) CudaNdarray_HOST_DIMS(self)[i];
874 }
875 CudaNdarray * rval = (CudaNdarray *)CudaNdarray_New();
876 if (!rval || CudaNdarray_alloc_contiguous(rval, self->nd, CudaNdarray_HOST_DIMS(self)))
877 {
878 Py_XDECREF(rval);
879 return NULL;
880 }
881
882 if(CudaNdarray_SIZE((CudaNdarray *)py_self)==0 && CudaNdarray_SIZE((CudaNdarray *)py_other)==0){
883 return (PyObject *) rval;
884 }
885
886 int threads_per_block = std::min(size, (unsigned int)NUM_VECTOR_OP_THREADS_PER_BLOCK);
887 int n_blocks = std::min(ceil_intdiv(size,(unsigned int)threads_per_block), (unsigned int)NUM_VECTOR_OP_BLOCKS);
888 kAdd_contiguous<<<n_blocks,threads_per_block>>>(
889 self->devdata, other->devdata, rval->devdata, size);
890 CNDA_THREAD_SYNC;
891 cudaError_t err = cudaGetLastError();
892 if( cudaSuccess != err)
893 {
894 PyErr_Format(PyExc_RuntimeError, "Cuda error: %s: %s.\n", "kAdd", cudaGetErrorString(err));
895 Py_DECREF(rval);
896 return NULL;
897 }
898 return (PyObject *) rval;
899 }
900
901 /*
902 #define decl_k_elemwise_binary_inplace_rowmajor_3(name, F) \
903 __global__ void name(const int d0, const int d1, const int d2,\
904 float* a, const int sA0, const int sA1, const int sA2,\
905 const float* b, const int sB0, const int sB1, const int sB2){\
906 for (int i0 = blockIdx.x; i0 < d0; i0 += gridDim.x){\
907 for (int i1 = blockIdx.y; i1 < d1; i1 += gridDim.y){\
908 for (int i2 = threadIdx.x; i2 < d2; i2 += blockDim.x){\
909 F(a[i0*sA0 + i1*sA1 + i2*sA2], b[i0*sB0 + i1*sB1 + i2*sB2]); \
910 }\
911 }\
912 }\
913 }
914
915 #define decl_k_elemwise_binary_inplace_rowmajor_4(name, F) \
916 __global__ void name(const int d0, const int d1, const int d2, const int d3,\
917 float* a, const int sA0, const int sA1,\
918 const int sA2, const int sA3,\
919 const float* b, const int sB0, const int sB1,\
920 const int sB2, const int sB3){\
921 for (int i0 = blockIdx.x; i0 < d0; i0 += gridDim.x){\
922 for (int i1 = blockIdx.y; i1 < d1; i1 += gridDim.y){\
923 for (int i2 = threadIdx.x; i2 < d2; i2 += blockDim.x){\
924 for (int i3 = threadIdx.y; i3 < d3; i3 += blockDim.y){\
925 F(a[i0*sA0 + i1*sA1 + i2*sA2 + i3*sA3], b[i0*sB0 + i1*sB1 + i2*sB2 + i3*sB3]); \
926 }\
927 }\
928 }\
929 }\
930 }
931
932 template<typename T> __device__ T binary_iadd(T a, T b) { a = a+b; }
933 template<typename T> __device__ T binary_idiv(T a, T b) { a = a/b; }
934
935 decl_k_elemwise_binary_inplace_rowmajor_3(k_iAdd_3, binary_iadd<float>)
936 decl_k_elemwise_binary_inplace_rowmajor_4(k_iAdd_4, binary_iadd<float>)
937 decl_k_elemwise_binary_inplace_rowmajor_3(k_iDiv_3, binary_idiv<float>)
938 decl_k_elemwise_binary_inplace_rowmajor_4(k_iDiv_4, binary_idiv<float>)
939 */
940 __global__ void k_iAdd_3(const int d0, const int d1, const int d2,\
941 float* a, const int sA0, const int sA1, const int sA2,\
942 const float* b, const int sB0, const int sB1, const int sB2){\
943 for (int i0 = blockIdx.x; i0 < d0; i0 += gridDim.x){\
944 for (int i1 = blockIdx.y; i1 < d1; i1 += gridDim.y){\
945 for (int i2 = threadIdx.x; i2 < d2; i2 += blockDim.x){\
946 a[i0*sA0 + i1*sA1 + i2*sA2]+= b[i0*sB0 + i1*sB1 + i2*sB2]; \
947 }\
948 }\
949 }\
950 }
951
952 __global__ void k_iAdd_4(const int d0, const int d1, const int d2, const int d3,\
953 float* a, const int sA0, const int sA1,\
954 const int sA2, const int sA3,\
955 const float* b, const int sB0, const int sB1,\
956 const int sB2, const int sB3){\
957 for (int i0 = blockIdx.x; i0 < d0; i0 += gridDim.x){\
958 for (int i1 = blockIdx.y; i1 < d1; i1 += gridDim.y){\
959 for (int i2 = threadIdx.x; i2 < d2; i2 += blockDim.x){\
960 for (int i3 = threadIdx.y; i3 < d3; i3 += blockDim.y){\
961 a[i0*sA0 + i1*sA1 + i2*sA2 + i3*sA3] += b[i0*sB0 + i1*sB1 + i2*sB2 + i3*sB3]; \
962 }\
963 }\
964 }\
965 }\
966 }
967
968 __global__ void k_iDiv_3(const int d0, const int d1, const int d2,\
969 float* a, const int sA0, const int sA1, const int sA2,\
970 const float* b, const int sB0, const int sB1, const int sB2){\
971 for (int i0 = blockIdx.x; i0 < d0; i0 += gridDim.x){\
972 for (int i1 = blockIdx.y; i1 < d1; i1 += gridDim.y){\
973 for (int i2 = threadIdx.x; i2 < d2; i2 += blockDim.x){\
974 a[i0*sA0 + i1*sA1 + i2*sA2]/= b[i0*sB0 + i1*sB1 + i2*sB2]; \
975 }\
976 }\
977 }\
978 }
979
980 __global__ void k_iDiv_4(const int d0, const int d1, const int d2, const int d3,\
981 float* a, const int sA0, const int sA1,\
982 const int sA2, const int sA3,\
983 const float* b, const int sB0, const int sB1,\
984 const int sB2, const int sB3){\
985 for (int i0 = blockIdx.x; i0 < d0; i0 += gridDim.x){\
986 for (int i1 = blockIdx.y; i1 < d1; i1 += gridDim.y){\
987 for (int i2 = threadIdx.x; i2 < d2; i2 += blockDim.x){\
988 for (int i3 = threadIdx.y; i3 < d3; i3 += blockDim.y){\
989 a[i0*sA0 + i1*sA1 + i2*sA2 + i3*sA3] /= b[i0*sB0 + i1*sB1 + i2*sB2 + i3*sB3]; \
990 }\
991 }\
992 }\
993 }\
994 }
995
996 static PyObject *
997 CudaNdarray_inplace_add_div(PyObject* py_self, PyObject * py_other, int fct_nb)
998 {
999 int verbose = 0;
1000 if (! CudaNdarray_Check(py_self)) {
1001 PyErr_SetString(PyExc_TypeError, "CudaNdarray_inplace_add_div need a CudaNdarray on left");
1002 return NULL;
1003 }
1004 if (! CudaNdarray_Check(py_other)) {
1005 PyErr_SetString(PyExc_TypeError, "CudaNdarray_inplace_add_div need a CudaNdarray on right");
1006 return NULL;
1007 }
1008 if (fct_nb<0 || fct_nb>1){
1009 PyErr_SetString(PyExc_TypeError, "CudaNdarray_inplace_add_div fct_nb param supported are only 0 and 1.");
1010 return NULL;
1011 }
1012
1013 CudaNdarray * self = (CudaNdarray *)py_self;
1014 CudaNdarray * other = (CudaNdarray *)py_other;
1015
1016 if (verbose) fprintf(stderr, "INPLACE ADD/DIV for self->nd=%d other->nd=%d\n",
1017 self->nd, other->nd);
1018
1019 //standard elemwise size checks
1020 if (self->nd != other->nd)
1021 {
1022 PyErr_Format(PyExc_TypeError, "CudaNdarray_inplace_add_div: need same number of dims. Got %d and %d", self->nd, other->nd);
1023 return NULL;
1024 }
1025 //standard elemwise dim checks
1026 unsigned int size = 1;
1027 for (int i = 0; i< self->nd; ++i)
1028 {
1029 if ((CudaNdarray_HOST_DIMS(self)[i] != CudaNdarray_HOST_DIMS(other)[i])
1030 && (CudaNdarray_HOST_DIMS(other)[i] != 1))
1031 {
1032 PyErr_SetString(PyExc_TypeError, "need same dimensions (or broadcastable dimension)");
1033 return NULL;
1034 }
1035 size *= (unsigned int) CudaNdarray_HOST_DIMS(self)[i];
1036 }
1037
1038 if(CudaNdarray_SIZE((CudaNdarray *)py_self)==0 && CudaNdarray_SIZE((CudaNdarray *)py_other)==0){
1039 Py_INCREF(py_self);
1040 return py_self;
1041 }
1042 void (*k_iop_3)(const int, const int, const int,
1043 float*, const int, const int, const int,
1044 const float*, const int, const int, const int);
1045 void (*k_iop_4)(const int, const int, const int, const int,
1046 float*, const int, const int,
1047 const int, const int,
1048 const float*, const int, const int,
1049 const int, const int);
1050 if(fct_nb == 0){
1051 k_iop_3 = k_iAdd_3;
1052 k_iop_4 = k_iAdd_4;
1053 }else if(fct_nb == 1){
1054 k_iop_3 = k_iDiv_3;
1055 k_iop_4 = k_iDiv_4;
1056 }
1057
1058 switch(self->nd)
1059 {
1060 case 0:
1061 {
1062 dim3 n_blocks(1, 1, 1);
1063 dim3 n_threads(1);
1064 k_iop_3<<<n_blocks, n_threads>>>(1,
1065 1, //CudaNdarray_HOST_DIMS(self)[0],
1066 1, //CudaNdarray_HOST_DIMS(self)[0],
1067 CudaNdarray_DEV_DATA(self),
1068 1,
1069 1, //CudaNdarray_HOST_STRIDES(self)[0],
1070 CudaNdarray_HOST_STRIDES(self)[0],
1071 CudaNdarray_DEV_DATA(other),
1072 1,
1073 1, //CudaNdarray_HOST_STRIDES(other)[0],
1074 CudaNdarray_HOST_STRIDES(other)[0]);
1075 CNDA_THREAD_SYNC;
1076 cudaError_t err = cudaGetLastError();
1077 if( cudaSuccess != err)
1078 {
1079 PyErr_Format(PyExc_RuntimeError, "Cuda error: %s: %s.\n", "k_iop_3", cudaGetErrorString(err));
1080 return NULL;
1081 }
1082 Py_INCREF(py_self);
1083 return py_self;
1084 }
1085 case 1:
1086 {
1087 dim3 n_blocks(1, 1, 1);
1088 dim3 n_threads(
1089 std::min(CudaNdarray_HOST_DIMS(self)[0], NUM_VECTOR_OP_THREADS_PER_BLOCK)
1090 );
1091 k_iop_3<<<n_blocks, n_threads>>>(1,
1092 1, //CudaNdarray_HOST_DIMS(self)[0],
1093 CudaNdarray_HOST_DIMS(self)[0],
1094 CudaNdarray_DEV_DATA(self),
1095 1,
1096 1, //CudaNdarray_HOST_STRIDES(self)[0],
1097 CudaNdarray_HOST_STRIDES(self)[0],
1098 CudaNdarray_DEV_DATA(other),
1099 1,
1100 1, //CudaNdarray_HOST_STRIDES(other)[0],
1101 CudaNdarray_HOST_STRIDES(other)[0]);
1102 CNDA_THREAD_SYNC;
1103 cudaError_t err = cudaGetLastError();
1104 if( cudaSuccess != err)
1105 {
1106 PyErr_Format(PyExc_RuntimeError, "Cuda error: %s: %s.\n", "k_iop_3", cudaGetErrorString(err));
1107 return NULL;
1108 }
1109 Py_INCREF(py_self);
1110 return py_self;
1111 }
1112 case 2:
1113 {
1114 dim3 n_blocks(1,
1115 std::min(CudaNdarray_HOST_DIMS(self)[0], NUM_VECTOR_OP_BLOCKS)
1116 );
1117 dim3 n_threads(
1118 std::min(CudaNdarray_HOST_DIMS(self)[1], NUM_VECTOR_OP_THREADS_PER_BLOCK)
1119 );
1120 k_iop_3<<<n_blocks, n_threads>>>(1,
1121 CudaNdarray_HOST_DIMS(self)[0],
1122 CudaNdarray_HOST_DIMS(self)[1],
1123 CudaNdarray_DEV_DATA(self),
1124 1,
1125 CudaNdarray_HOST_STRIDES(self)[0],
1126 CudaNdarray_HOST_STRIDES(self)[1],
1127 CudaNdarray_DEV_DATA(other),
1128 1,
1129 CudaNdarray_HOST_STRIDES(other)[0],
1130 CudaNdarray_HOST_STRIDES(other)[1]);
1131 CNDA_THREAD_SYNC;
1132 cudaError_t err = cudaGetLastError();
1133 if( cudaSuccess != err)
1134 {
1135 PyErr_Format(PyExc_RuntimeError, "Cuda error: %s: %s.\n", "k_iop_3", cudaGetErrorString(err));
1136 return NULL;
1137 }
1138 Py_INCREF(py_self);
1139 return py_self;
1140 }
1141 case 3:
1142 {
1143 dim3 n_blocks(
1144 std::min(CudaNdarray_HOST_DIMS(self)[0], NUM_VECTOR_OP_BLOCKS),
1145 CudaNdarray_HOST_DIMS(self)[1]
1146 );
1147 while (n_blocks.x * n_blocks.y > NUM_VECTOR_OP_BLOCKS) n_blocks.y /= 2;
1148 dim3 n_threads(
1149 std::min(CudaNdarray_HOST_DIMS(self)[2], NUM_VECTOR_OP_THREADS_PER_BLOCK)
1150 );
1151 k_iop_3<<<n_blocks, n_threads>>>(
1152 CudaNdarray_HOST_DIMS(self)[0],
1153 CudaNdarray_HOST_DIMS(self)[1],
1154 CudaNdarray_HOST_DIMS(self)[2],
1155 CudaNdarray_DEV_DATA(self),
1156 CudaNdarray_HOST_STRIDES(self)[0],
1157 CudaNdarray_HOST_STRIDES(self)[1],
1158 CudaNdarray_HOST_STRIDES(self)[2],
1159 CudaNdarray_DEV_DATA(other),
1160 CudaNdarray_HOST_STRIDES(other)[0],
1161 CudaNdarray_HOST_STRIDES(other)[1],
1162 CudaNdarray_HOST_STRIDES(other)[2]);
1163 CNDA_THREAD_SYNC;
1164 cudaError_t err = cudaGetLastError();
1165 if( cudaSuccess != err)
1166 {
1167 PyErr_Format(PyExc_RuntimeError, "Cuda error: %s: %s.\n", "k_iop_3", cudaGetErrorString(err));
1168 return NULL;
1169 }
1170 Py_INCREF(py_self);
1171 return py_self;
1172 }
1173 case 4:
1174 {
1175 dim3 n_blocks(
1176 std::min(CudaNdarray_HOST_DIMS(self)[0], NUM_VECTOR_OP_BLOCKS),
1177 CudaNdarray_HOST_DIMS(self)[1]
1178 );
1179 while (n_blocks.x * n_blocks.y > NUM_VECTOR_OP_BLOCKS) n_blocks.y /= 2;
1180 dim3 n_threads(
1181 std::min(CudaNdarray_HOST_DIMS(self)[2], NUM_VECTOR_OP_THREADS_PER_BLOCK)
1182 );
1183 k_iop_4<<<n_blocks, n_threads>>>(
1184 CudaNdarray_HOST_DIMS(self)[0],
1185 CudaNdarray_HOST_DIMS(self)[1],
1186 CudaNdarray_HOST_DIMS(self)[2],
1187 CudaNdarray_HOST_DIMS(self)[3],
1188 CudaNdarray_DEV_DATA(self),
1189 CudaNdarray_HOST_STRIDES(self)[0],
1190 CudaNdarray_HOST_STRIDES(self)[1],
1191 CudaNdarray_HOST_STRIDES(self)[2],
1192 CudaNdarray_HOST_STRIDES(self)[3],
1193 CudaNdarray_DEV_DATA(other),
1194 CudaNdarray_HOST_STRIDES(other)[0],
1195 CudaNdarray_HOST_STRIDES(other)[1],
1196 CudaNdarray_HOST_STRIDES(other)[2],
1197 CudaNdarray_HOST_STRIDES(other)[3]);
1198 CNDA_THREAD_SYNC;
1199 cudaError_t err = cudaGetLastError();
1200 if( cudaSuccess != err)
1201 {
1202 PyErr_Format(PyExc_RuntimeError, "Cuda error: %s: %s.\n", "k_iop_4", cudaGetErrorString(err));
1203 return NULL;
1204 }
1205 Py_INCREF(py_self);
1206 return py_self;
1207 }
1208 case 5:
1209 {
1210 dim3 n_blocks(
1211 std::min(CudaNdarray_HOST_DIMS(self)[1], NUM_VECTOR_OP_BLOCKS),
1212 CudaNdarray_HOST_DIMS(self)[2]
1213 );
1214 while (n_blocks.x * n_blocks.y > NUM_VECTOR_OP_BLOCKS) n_blocks.y /= 2;
1215 dim3 n_threads(
1216 std::min(CudaNdarray_HOST_DIMS(self)[3], NUM_VECTOR_OP_THREADS_PER_BLOCK)
1217 );
1218 for (int i = 0; i < CudaNdarray_HOST_DIMS(self)[0]; ++i)
1219 {
1220 k_iop_4<<<n_blocks, n_threads>>>(
1221 CudaNdarray_HOST_DIMS(self)[1],
1222 CudaNdarray_HOST_DIMS(self)[2],
1223 CudaNdarray_HOST_DIMS(self)[3],
1224 CudaNdarray_HOST_DIMS(self)[4],
1225 CudaNdarray_DEV_DATA(self) + i * CudaNdarray_HOST_STRIDES(self)[0],
1226 CudaNdarray_HOST_STRIDES(self)[1],
1227 CudaNdarray_HOST_STRIDES(self)[2],
1228 CudaNdarray_HOST_STRIDES(self)[3],
1229 CudaNdarray_HOST_STRIDES(self)[4],
1230 CudaNdarray_DEV_DATA(other) + i * CudaNdarray_HOST_STRIDES(other)[0],
1231 CudaNdarray_HOST_STRIDES(other)[1],
1232 CudaNdarray_HOST_STRIDES(other)[2],
1233 CudaNdarray_HOST_STRIDES(other)[3],
1234 CudaNdarray_HOST_STRIDES(other)[4]);
1235 CNDA_THREAD_SYNC;
1236 cudaError_t err = cudaGetLastError();
1237 if( cudaSuccess != err)
1238 {
1239 PyErr_Format(PyExc_RuntimeError, "Cuda error: %s: %s.\n", "k_iop_4", cudaGetErrorString(err));
1240 return NULL;
1241 }
1242 }
1243 Py_INCREF(py_self);
1244 return py_self;
1245 }
1246 }
1247
1248 PyErr_Format(PyExc_NotImplementedError, "inplace_add w nd=%i\n", self->nd);
1249 return NULL;
1250 }
1251
1252 /*
1253 * We need this inplace Add to support IncSubTensor
1254 */
1255 // Will be called by __iadd__ in Python
1256 static PyObject *
1257 CudaNdarray_inplace_add(PyObject* py_self, PyObject * py_other){
1258 PyObject * rval = CudaNdarray_inplace_add_div(py_self, py_other, 0);
1259 //We should not increment the refcount as we are doing inplace operation
1260 //And in this syntax, their is no additional reference created!
1261 return rval;
1262 }
1263
1264 /*
1265 * We need this inplace div for cuda/tests/test_basic_ops.py:test_shared_options
1266 */
1267 // Will be called by __idiv__ in Python
1268 static PyObject *
1269 CudaNdarray_inplace_div(PyObject* py_self, PyObject * py_other){
1270 PyObject * rval = CudaNdarray_inplace_add_div(py_self, py_other, 1);
1271 //We should not increment the refcount as we are doing inplace operation
1272 //And in this syntax, their is no additional reference created!
1273 return rval;
1274 }
1275
1276 static PyNumberMethods CudaNdarrayNumberMethods =
1277 {
1278 (binaryfunc)CudaNdarray_add, //binaryfunc nb_add; __add__
1279 0, //binaryfunc nb_subtract; __sub__
1280 0, //binaryfunc nb_multiply; __mul__
1281 0, //binaryfunc nb_divide; __div__
1282 0, //binaryfunc nb_remainder; __mod__
1283 0, //binaryfunc nb_divmod; __divmod__
1284 0, //ternaryfunc nb_power; __pow__
1285 0, //unaryfunc nb_negative; __neg__
1286 0, //unaryfunc nb_positive; __pos__
1287 0, //unaryfunc nb_absolute; __abs__
1288 0, //inquiry nb_nonzero; __nonzero__ /* Used by PyObject_IsTrue */
1289 0, //unaryfunc nb_invert; __invert__
1290 0, //binaryfunc nb_lshift; __lshift__
1291 0, //binaryfunc nb_rshift; __rshift__
1292 0, //binaryfunc nb_and; __and__
1293 0, //binaryfunc nb_xor; __xor__
1294 0, //binaryfunc nb_or; __or__
1295 0, //coercion nb_coerce; __coerce__ /* Used by the coerce() function */
1296 0, //unaryfunc nb_int; __int__
1297 0, //unaryfunc nb_long; __long__
1298 0, //unaryfunc nb_float; __float__
1299 0, //unaryfunc nb_oct; __oct__
1300 0, //unaryfunc nb_hex; __hex__
1301
1302 /* Added in release 2.0 */
1303 (binaryfunc)CudaNdarray_inplace_add, //binaryfunc nb_inplace_add; __iadd__
1304 0, //binaryfunc nb_inplace_subtract; __isub__
1305 0, //binaryfunc nb_inplace_multiply; __imul__
1306 (binaryfunc)CudaNdarray_inplace_div, //binaryfunc nb_inplace_divide; __idiv__
1307 0, //binaryfunc nb_inplace_remainder; __imod__
1308 0, //ternaryfunc nb_inplace_power; __ipow__
1309 0, //binaryfunc nb_inplace_lshift; __ilshift__
1310 0, //binaryfunc nb_inplace_rshift; __irshift__
1311 0, //binaryfunc nb_inplace_and; __iand__
1312 0, //binaryfunc nb_inplace_xor; __ixor__
1313 0, //binaryfunc nb_inplace_or; __ior__
1314
1315 /* Added in release 2.2 */
1316 0, //binaryfunc nb_floor_divide; __floordiv__
1317 0, //binaryfunc nb_true_divide; __truediv__
1318 0, //binaryfunc nb_inplace_floor_divide; __ifloordiv__
1319 0, //binaryfunc nb_inplace_true_divide; __itruediv__
1320
1321 #if PY_MINOR_VERSION > 4
1322 /* Added in release 2.5 */
1323 0 //unaryfunc nb_index; __index__
1324 #endif
1325 };
1326
1327
1328 /////////////////////
1329 // Mapping protocol
1330 /////////////////////
1331
1332 // Will by called by __len__ in Python
1333 static Py_ssize_t
1334 CudaNdarray_len(PyObject * py_self)
1335 {
1336 CudaNdarray * self = (CudaNdarray*) py_self;
1337 if (self->nd <= 0)
1338 {
1339 return (Py_ssize_t) 0;
1340 }
1341 else
1342 {
1343 return (Py_ssize_t) CudaNdarray_HOST_DIMS(self)[0];
1344 }
1345 }
1346
1347 // Will by called by __getitem__ in Python
1348 static PyObject *
1349 CudaNdarray_Subscript(PyObject * py_self, PyObject * key)
1350 {
1351 int verbose = 0;
1352 if (verbose) fprintf(stderr, "Subscript .... \n");
1353 CudaNdarray * self = (CudaNdarray*) py_self;
1354 PyObject * py_rval = NULL;
1355 CudaNdarray * rval = NULL;
1356 PyObject * intobj = NULL;
1357
1358 //PyObject_Print(key, stderr, 0);
1359
1360 if (key == Py_Ellipsis)
1361 {
1362 Py_INCREF(py_self);
1363 return py_self;
1364 }
1365 if ((intobj=PyNumber_Int(key))) //INDEXING BY INTEGER
1366 //else if (PyInt_Check(key)) //INDEXING BY INTEGER
1367 {
1368 int d_idx = PyInt_AsLong(intobj);
1369 Py_DECREF(intobj); intobj=NULL;
1370 //int d_idx = PyInt_AsLong(key);
1371 if (self->nd == 0)
1372 {
1373 PyErr_SetString(PyExc_IndexError, "0-d arrays can't be indexed");
1374 return NULL;
1375 }
1376 int d_dim = CudaNdarray_HOST_DIMS(self)[0];
1377 int offset = 0;
1378
1379 if ((d_idx >= 0) && (d_idx < d_dim))
1380 {
1381 //normal indexing
1382 offset += d_idx * CudaNdarray_HOST_STRIDES(self)[0];
1383 }
1384 else if ((d_idx < 0) && (d_idx >= -d_dim))
1385 {
1386 //end-based indexing
1387 // d_idx is negative
1388 offset += (d_dim + d_idx) * CudaNdarray_HOST_STRIDES(self)[0];
1389 }
1390 else
1391 {
1392 PyErr_SetString(PyExc_IndexError, "index out of bounds");
1393 return NULL;
1394 }
1395
1396 //allocate our subtensor view
1397 py_rval = CudaNdarray_new_nd(self->nd - 1);
1398 rval = (CudaNdarray*) py_rval;
1399 if (!rval) return NULL;
1400 assert (0 == rval->data_allocated);
1401
1402 //initialize the view's data pointer to our own.
1403 if (CudaNdarray_set_device_data(rval, CudaNdarray_DEV_DATA(self) + offset, self))
1404 {
1405 Py_DECREF(rval);
1406 return NULL;
1407 }
1408 for (int d = 1; d < self->nd; ++d)
1409 {
1410 CudaNdarray_set_stride(rval, d-1, CudaNdarray_HOST_STRIDES(self)[d]);
1411 CudaNdarray_set_dim(rval, d-1, CudaNdarray_HOST_DIMS(self)[d]);
1412 }
1413 }
1414 else
1415 {
1416 PyErr_Clear();
1417 }
1418 if (PySlice_Check(key)) //INDEXING BY SLICE
1419 {
1420 if (self->nd == 0)
1421 {
1422 PyErr_SetString(PyExc_ValueError, "cannot slice a 0-d array");
1423 return NULL;
1424 }
1425
1426 int d_dim = CudaNdarray_HOST_DIMS(self)[0];
1427 Py_ssize_t start, stop, step, slen;
1428 if (PySlice_GetIndicesEx((PySliceObject*)key, d_dim, &start, &stop, &step, &slen))
1429 {
1430 return NULL;
1431 }
1432 if (verbose)
1433 {
1434 std::cerr << "start " << start << "\n";
1435 std::cerr << "stop " << stop << "\n";
1436 std::cerr << "step " << step << "\n";
1437 std::cerr << "slen " << slen << "\n";
1438 }
1439
1440 //allocate our subtensor view
1441 py_rval = CudaNdarray_new_nd(self->nd);
1442 rval = (CudaNdarray*) py_rval;
1443 if (!rval) return NULL;
1444 assert (0 == rval->data_allocated);
1445
1446
1447 //initialize the view's data pointer to our own.
1448 if (CudaNdarray_set_device_data(rval,
1449 CudaNdarray_DEV_DATA(self) + start * CudaNdarray_HOST_STRIDES(self)[0],
1450 self))
1451 {
1452 Py_DECREF(rval);
1453 return NULL;
1454 }
1455 //initialize dimension 0 of rval
1456 CudaNdarray_set_stride(rval, 0, step * CudaNdarray_HOST_STRIDES(self)[0]);
1457 CudaNdarray_set_dim(rval, 0, slen);
1458 if (verbose) std::cerr << "rval stride " << CudaNdarray_HOST_STRIDES(rval)[0] << "\n";
1459 // initialize dimensions > 0 of rval
1460 for (int d = 1; d < self->nd; ++d)
1461 {
1462 CudaNdarray_set_stride(rval, d, CudaNdarray_HOST_STRIDES(self)[d]);
1463 CudaNdarray_set_dim(rval, d, CudaNdarray_HOST_DIMS(self)[d]);
1464 }
1465 }
1466 if (PyTuple_Check(key)) //INDEXING BY TUPLE
1467 {
1468 //elements of the tuple can be either integers or slices
1469 //the dimensionality of the view we will return is diminished for each slice in the tuple
1470
1471 if (PyTuple_Size(key) > self->nd)
1472 {
1473 PyErr_SetString(PyExc_IndexError, "index error");
1474 return NULL;
1475 }
1476
1477 //calculate the number of dimensions in the return value
1478 int rval_nd = self->nd;
1479 for (int d = 0; d < PyTuple_Size(key); ++d)
1480 {
1481 //On some paltform PyInt_Check(<type 'numpy.int64'>) return true, other it return false.
1482 //So we use PyArray_IsAnyScalar that should covert everything.
1483 rval_nd -= PyArray_IsAnyScalar(PyTuple_GetItem(key, d));
1484 }
1485
1486 //allocate our subtensor view
1487 py_rval = CudaNdarray_new_nd(rval_nd);
1488 rval = (CudaNdarray*) py_rval;
1489 if (!rval) return NULL;
1490 assert (0 == rval->data_allocated);
1491
1492 //initialize the view's data pointer to our own.
1493 if (CudaNdarray_set_device_data(rval, CudaNdarray_DEV_DATA(self), self))
1494 {
1495 Py_DECREF(rval);
1496 return NULL;
1497 }
1498
1499 // rval_d will refer to the current dimension in the rval.
1500 // It will not be incremented for integer keys, but will be incremented for slice
1501 // keys
1502 int rval_d = 0;
1503
1504 for (int d = 0; d < self->nd; ++d)
1505 {
1506 // keys can be shorter than self->nd.
1507 // when that happens, it means that the remaining dimensions are "full slices"
1508 if (d >=PyTuple_Size(key))
1509 {
1510 CudaNdarray_set_stride(rval, rval_d, CudaNdarray_HOST_STRIDES(self)[d]);
1511 CudaNdarray_set_dim(rval, rval_d, CudaNdarray_HOST_DIMS(self)[d]);
1512 ++rval_d;
1513 }
1514 else
1515 {
1516 PyObject * key_d = PyTuple_GetItem(key, d);
1517
1518 if (PySlice_Check(key_d))
1519 {
1520 Py_ssize_t start, stop, step, slen;
1521 if (PySlice_GetIndicesEx((PySliceObject*)key_d, CudaNdarray_HOST_DIMS(self)[d], &start, &stop, &step, &slen))
1522 {
1523 Py_DECREF(rval);
1524 return NULL;
1525 }
1526 rval->devdata += start * CudaNdarray_HOST_STRIDES(self)[d];
1527 CudaNdarray_set_stride(rval, rval_d, step * CudaNdarray_HOST_STRIDES(self)[d]);
1528 CudaNdarray_set_dim(rval, rval_d, slen);
1529 if (0)
1530 {
1531 std::cerr << "start " << start << "\n";
1532 std::cerr << "stop " << stop << "\n";
1533 std::cerr << "step " << step << "\n";
1534 std::cerr << "slen " << slen << "\n";
1535 }
1536 ++rval_d;
1537 }
1538 else if ((intobj=PyNumber_Int(key_d)))
1539 {
1540 assert(PyArray_IsAnyScalar(key_d));
1541 int d_idx = PyInt_AsLong(intobj);
1542 Py_DECREF(intobj);
1543 intobj = NULL;
1544 int d_dim = CudaNdarray_HOST_DIMS(self)[d];
1545
1546 if ((d_idx >= 0) && (d_idx < d_dim))
1547 {
1548 //normal indexing
1549 rval->devdata += d_idx * CudaNdarray_HOST_STRIDES(self)[d];
1550 }
1551 else if ((d_idx < 0) && (d_idx >= -d_dim))
1552 {
1553 //end-based indexing
1554 rval->devdata += (d_dim + d_idx) * CudaNdarray_HOST_STRIDES(self)[d];
1555 }
1556 else
1557 {
1558 PyErr_SetString(PyExc_IndexError, "index out of bounds");
1559 Py_DECREF(rval);
1560 return NULL;
1561 }
1562 }
1563 else
1564 {
1565 PyErr_Clear(); // clear the error set by PyNumber_Int
1566 PyErr_SetString(PyExc_IndexError, "index must be either int or slice");
1567 Py_DECREF(rval);
1568 return NULL;
1569 }
1570 }
1571 }
1572 }
1573 if (py_rval)
1574 {
1575 if (verbose) fprint_CudaNdarray(stderr, self);
1576 if (verbose) fprint_CudaNdarray(stderr, rval);
1577 }
1578 else
1579 {
1580 PyErr_SetString(PyExc_NotImplementedError, "Unknown key type");
1581 return NULL;
1582 }
1583 return py_rval;
1584 }
1585
1586 // Will by called by __setitem__ in Python
1587 // See http://docs.python.org/dev/py3k/c-api/object.html#PyObject_SetItem
1588 // Doesn't handle broadcasting, e.g. a[:] = 5
1589 // Can only be assigned from a CudaNdarray on the right side
1590 // Or a ndarray when the left side part is c contiguous.
1591 static int
1592 CudaNdarray_setitem(PyObject *o, PyObject *key, PyObject *value)
1593 {
1594 if(CudaNdarray_Check(o) && PyArray_Check(value)){
1595 // We try to copy directly into this CudaNdarray from the ndarray
1596 CudaNdarray* rval = (CudaNdarray*)CudaNdarray_Subscript(o, key);
1597 int typenum = PyArray_TYPE(value);
1598
1599 if(!rval){
1600 // CudaNdarray_Subscript failed and set the error msg.
1601 Py_XDECREF(rval);
1602 return -1;
1603 }
1604 if (typenum != REAL_TYPENUM){
1605 PyErr_SetString(PyExc_TypeError, "CudaNdarray.__setitem__: can only copy from float32 arrays");
1606 Py_XDECREF(rval);
1607 return -1;
1608 }
1609 if(! CudaNdarray_is_c_contiguous(rval)){
1610 PyErr_SetString(PyExc_NotImplementedError, "CudaNdarray.__setitem__: When the new value is an ndarray the part where we copy it to must be c contiguous.");
1611 Py_XDECREF(rval);
1612 return -1;
1613 }
1614 if(rval->nd != ((PyArrayObject*)value)->nd){
1615 PyErr_Format(PyExc_NotImplementedError, "CudaNdarray.__setitem__: need same number of dims. destination nd=%d, source nd=%d. No broadcasting implemented.",
1616 rval->nd,((PyArrayObject*)value)->nd);
1617 Py_XDECREF(rval);
1618 return -1;
1619 }
1620 for(int i=0 ; i<rval->nd ; i++){
1621 if(CudaNdarray_HOST_DIMS(rval)[i] != ((PyArrayObject*)value)->dimensions[i]){
1622 PyErr_Format(PyExc_ValueError, "CudaNdarray.__setitem__: need same dimensions for dim %d, destination=%d, source=%ld",
1623 i,
1624 CudaNdarray_HOST_DIMS(rval)[i],
1625 (long int)(((PyArrayObject*)value)->dimensions[i]));
1626 Py_XDECREF(rval);
1627 return -1;
1628 }
1629 }
1630 PyArrayObject * py_v = (PyArrayObject*)PyArray_ContiguousFromAny((PyObject*)value, typenum,
1631 rval->nd, rval->nd);
1632 cublasSetVector(PyArray_SIZE(py_v),
1633 sizeof(real),
1634 PyArray_DATA(py_v), 1,
1635 rval->devdata, 1);
1636 CNDA_THREAD_SYNC;
1637 Py_XDECREF(py_v);
1638 Py_XDECREF(rval);
1639 if (CUBLAS_STATUS_SUCCESS != cublasGetError()){
1640 PyErr_SetString(PyExc_RuntimeError, "CudaNdarray.__setitem__: error copying ndarray data to device memory");
1641 return -1;
1642 }
1643 return 0;
1644 }
1645
1646
1647 if(!CudaNdarray_Check(o) || !CudaNdarray_Check(value))
1648 {
1649 PyErr_SetString(PyExc_TypeError, "CudaNdarray.__setitem__: left must be a CudaNdarrays and right must be a CudaNdarrays or ndarray");
1650 return -1;
1651 }
1652
1653 CudaNdarray* rval = (CudaNdarray*)CudaNdarray_Subscript(o, key);
1654
1655 if(rval == NULL)
1656 {
1657 // Actually error string was probably set if we get a NULL, so we leave it as it is
1658 //PyErr_SetString(PyExc_RuntimeError, "__getitem__ returned an error");
1659 return -1;
1660 }
1661 else if(rval != (CudaNdarray*)o &&
1662 (rval->data_allocated ||
1663 // The new array should have a base
1664 !(((CudaNdarray*)rval)->base) ||
1665 // If the original array has no base, the base of the new
1666 // array should be the original one
1667 (!((CudaNdarray*)o)->base && ((CudaNdarray*)rval)->base != o) ||
1668 // Else, the two arrays should have the same base
1669 (((CudaNdarray*)o)->base && ((CudaNdarray*)rval)->base != ((CudaNdarray*)o)->base)))
1670 {
1671 // This case shouldn't happen, based on what I see in Subscript
1672 // but just in case it happens sometime in the future
1673
1674 PyErr_Format(PyExc_RuntimeError, "__getitem__ must return a CudaNdarray that refers to the original CudaNdarray, not a copy. rval.base=%p o.base=%p o=%p",
1675 (((CudaNdarray*)rval)->base), ((CudaNdarray*)o)->base, o);
1676 Py_DECREF(rval);
1677 return -1;
1678 }
1679
1680 if (cnda_copy_structure_to_device(rval))
1681 {
1682 PyErr_SetString(PyExc_RuntimeError, "CudaNdarray.__setitem__: syncing structure to device failed");
1683 Py_DECREF(rval);
1684 return -1;
1685 }
1686
1687 PyObject *baseSavedForComparison = rval->base;
1688
1689 if(CudaNdarray_CopyFromCudaNdarray(rval, (CudaNdarray*)value, true))
1690 {
1691 Py_DECREF((PyObject*)rval);
1692 return -1;
1693 }
1694
1695 assert (rval->base == baseSavedForComparison);
1696 assert (rval->dev_structure_fresh);
1697
1698 // Clean up locally-created references
1699 Py_DECREF(rval);
1700
1701 return 0;
1702 }
1703
1704
1705 PyMappingMethods CudaNdarrayMappingMethods = {
1706 CudaNdarray_len, //lenfunc mp_length; __len__
1707 CudaNdarray_Subscript, //binaryfunc mp_subscript; __getitem__
1708 CudaNdarray_setitem //objobjargproc mp_ass_subscript; __setitem__
1709 };
1710
1711 ////////////////////
1712 //
1713 ////////////////////
1714
1715 static PyObject *
1716 CudaNdarray_get_shape(CudaNdarray *self, void *closure)
1717 {
1718 if (self->nd < 0)
1719 {
1720 PyErr_SetString(PyExc_ValueError, "CudaNdarray not initialized");
1721 return NULL;
1722 }
1723 PyObject * rval = PyTuple_New(self->nd);
1724 for (int i = 0; i < self->nd; ++i)
1725 {
1726 if (!rval || PyTuple_SetItem(rval, i, PyInt_FromLong(CudaNdarray_HOST_DIMS(self)[i])))
1727 {
1728 Py_XDECREF(rval);
1729 return NULL;
1730 }
1731
1732 }
1733 return rval;
1734 }
1735
1736 static int
1737 CudaNdarray_set_shape(CudaNdarray *self, PyObject *value, void *closure)
1738 {
1739 PyErr_SetString(PyExc_NotImplementedError, "TODO: call reshape");
1740 return -1;
1741 }
1742
1743 static PyObject *
1744 CudaNdarray_get_strides(CudaNdarray *self, void *closure)
1745 {
1746 if (self->nd < 0)
1747 {
1748 PyErr_SetString(PyExc_ValueError, "CudaNdarray not initialized");
1749 return NULL;
1750 }
1751 PyObject * rval = PyTuple_New(self->nd);
1752 for (int i = 0; i < self->nd; ++i)
1753 {
1754 if (!rval || PyTuple_SetItem(rval, i, PyInt_FromLong(CudaNdarray_HOST_STRIDES(self)[i])))
1755 {
1756 Py_XDECREF(rval);
1757 return NULL;
1758 }
1759
1760 }
1761 return rval;
1762 }
1763
1764 static int
1765 CudaNdarray_set_strides(CudaNdarray *self, PyObject *value, void *closure)
1766 {
1767 PyErr_SetString(PyExc_NotImplementedError, "");
1768 return -1;
1769 }
1770
1771 static PyObject *
1772 CudaNdarray_get_dev_data(CudaNdarray *self, void *closure)
1773 {
1774 float * p = CudaNdarray_DEV_DATA(self);
1775 //printf("get_dev_data %p %li \n", p, (long int)p );
1776 return PyInt_FromLong((long int) CudaNdarray_DEV_DATA(self));
1777 }
1778
1779 static int
1780 CudaNdarray_set_dev_data(CudaNdarray *self, PyObject *value, void *closure)
1781 {
1782 long int newdevdata = PyInt_AsLong(value);
1783 //printf("set_dev_data %p %li \n",(float*)newdevdata ,newdevdata);
1784 if (PyErr_Occurred())
1785 {
1786 return -1;
1787 }
1788 return CudaNdarray_set_device_data(self, (float*)newdevdata, (CudaNdarray*)self->base);
1789 }
1790
1791 static PyObject *
1792 CudaNdarray_get_dtype(CudaNdarray *self, void *closure)
1793 {
1794 return PyString_FromString("float32");
1795 }
1796
1797 static PyObject *
1798 CudaNdarray_get_ndim(CudaNdarray *self, void *closure)
1799 {
1800 return PyInt_FromLong(self->nd);
1801 }
1802
1803 static PyObject *
1804 CudaNdarray_get_base(CudaNdarray *self, void *closure)
1805 {
1806 PyObject * base = self->base;
1807 if (!base)
1808 {
1809 // We cannot return a NULL pointer, use None instead
1810 base = Py_None;
1811 }
1812 Py_INCREF(base);
1813 return base;
1814 }
1815
1816 static PyGetSetDef CudaNdarray_getset[] = {
1817 {"shape",
1818 (getter)CudaNdarray_get_shape,
1819 (setter)CudaNdarray_set_shape,
1820 "shape of this ndarray (tuple)",
1821 NULL},
1822 {"_strides",
1823 (getter)CudaNdarray_get_strides,
1824 (setter)CudaNdarray_set_strides,
1825 "data pointer strides (in elements)",
1826 NULL},
1827 //gpudata is needed to allow calling pycuda fct with CudaNdarray input.
1828 {"gpudata",
1829 (getter)CudaNdarray_get_dev_data,
1830 NULL,
1831 "device data pointer",
1832 NULL},
1833 {"_dev_data",
1834 (getter)CudaNdarray_get_dev_data,
1835 (setter)CudaNdarray_set_dev_data,
1836 "device data pointer",
1837 NULL},
1838 {"dtype",
1839 (getter)CudaNdarray_get_dtype,
1840 NULL,
1841 "The dtype of the element. Now always float32",
1842 NULL},
1843 {"size",
1844 (getter)CudaNdarray_SIZE_Object,
1845 NULL,
1846 "The number of elements in this object.",
1847 NULL},
1848 //mem_size is neede for pycuda.elementwise.ElementwiseKernel Why do they use size and mem_size of the same value?
1849 {"mem_size",
1850 (getter)CudaNdarray_SIZE_Object,
1851 NULL,
1852 "The number of elements in this object.",
1853 NULL},
1854 {"ndim",
1855 (getter)CudaNdarray_get_ndim,
1856 NULL,
1857 "The number of dimensions in this object.",
1858 NULL},
1859 {"base",
1860 (getter)CudaNdarray_get_base,
1861 NULL,
1862 "If this ndarray is a view, base is the original ndarray.",
1863 NULL},
1864
1865 {NULL, NULL, NULL, NULL} /* Sentinel */
1866 };
1867
1868
1869
1870 static PyTypeObject CudaNdarrayType =
1871 {
1872 PyObject_HEAD_INIT(NULL)
1873 0, /*ob_size*/
1874 "CudaNdarray", /*tp_name*/
1875 sizeof(CudaNdarray), /*tp_basicsize*/
1876 0, /*tp_itemsize*/
1877 (destructor)CudaNdarray_dealloc, /*tp_dealloc*/
1878 0, /*tp_print*/
1879 0, /*tp_getattr*/
1880 0, /*tp_setattr*/
1881 0, /*tp_compare*/
1882 0, /*tp_repr*/
1883 &CudaNdarrayNumberMethods, /*tp_as_number*/
1884 0, /*tp_as_sequence*/
1885 &CudaNdarrayMappingMethods,/*tp_as_mapping*/
1886 0, /*tp_hash */
1887 0, /*tp_call*/
1888 0, /*tp_str*/
1889 0, /*tp_getattro*/
1890 0, /*tp_setattro*/
1891 0, /*tp_as_buffer*/
1892 Py_TPFLAGS_DEFAULT | Py_TPFLAGS_BASETYPE | Py_TPFLAGS_CHECKTYPES, /*tp_flags*/
1893 "CudaNdarray objects", /* tp_doc */
1894 0, /* tp_traverse */
1895 0, /* tp_clear */
1896 0, /* tp_richcompare */
1897 0, /* tp_weaklistoffset */
1898 0, /* tp_iter */
1899 0, /* tp_iternext */
1900 CudaNdarray_methods, /* tp_methods */
1901 CudaNdarray_members, /* tp_members */
1902 CudaNdarray_getset, /* tp_getset */
1903 0, /* tp_base */
1904 0, /* tp_dict */
1905 0, /* tp_descr_get */
1906 0, /* tp_descr_set */
1907 0, /* tp_dictoffset */
1908 (initproc)CudaNdarray_init,/* tp_init */
1909 0, /* tp_alloc */
1910 CudaNdarray_new, /* tp_new */
1911 };
1912
1913 static __global__ void get_gpu_ptr_size(int* dst)
1914 {
1915 dst[0] = sizeof(float*);
1916 dst[1] = sizeof(int);
1917 }
1918
1919 PyObject *
1920 CudaNdarray_ptr_int_size(PyObject* _unused, PyObject* args)
1921 {
1922 int *gpu_data = (int*)device_malloc(sizeof(int)*2);
1923 if(gpu_data == NULL){
1924 return PyErr_Format(PyExc_MemoryError,
1925 "CudaNdarray_ptr_int_size: Can't allocate memory on the gpu.");
1926 }
1927 get_gpu_ptr_size<<<1,1>>>(gpu_data);
1928 if (cudaSuccess != cublasGetError()){
1929
1930 device_free(gpu_data);
1931 return PyErr_Format(PyExc_RuntimeError,
1932 "CudaNdarray_ptr_int_size: error when calling the gpu code.");
1933 }
1934
1935 // Transfer the result to cpu
1936 int gpu_sizes[] = {-1,-1};
1937 cublasGetVector(2, sizeof(int), gpu_data, 1, gpu_sizes, 1);
1938 device_free(gpu_data);
1939
1940 if (CUBLAS_STATUS_SUCCESS != cublasGetError()){
1941 PyErr_SetString(PyExc_RuntimeError, "error copying data to from memory");
1942 return NULL;
1943 }
1944 return Py_BuildValue("iiii", gpu_sizes[0], sizeof(float*), sizeof(int), gpu_sizes[1]);
1945 }
1946
1947 // Initialize the gpu.
1948 // Takes one optional parameter, the device number.
1949 // If provided, it sets that device to be the active device.
1950 // If not provided (usually just to test whether the gpu is available at all),
1951 // it does not set an active device.
1952 // Raises EnvironmentError or ValueError (as appropriate) if the initialization failed.
1953 PyObject *
1954 CudaNdarray_gpu_init(PyObject* _unused, PyObject* args)
1955 {
1956 int card_nb = 0;
1957 int card_number_provided = 1;
1958
1959 PyArg_ParseTuple(args, "|i", &card_nb); // if we're given something wildly invalid, this will throw a TypeError
1960
1961 if(PyTuple_Size(args) == 0) {
1962 card_number_provided = 0;
1963 card_nb = 0;
1964 }
1965
1966 int deviceCount;
1967 cudaError err = cudaGetDeviceCount(&deviceCount);
1968 if(cudaSuccess != err) {
1969 return PyErr_Format(PyExc_EnvironmentError,
1970 "Unable to get the number of gpus available: %s",
1971 cudaGetErrorString(cudaGetLastError()));
1972 }
1973
1974 // as soon as the first successful call to a cuda* function is made, a
1975 // gpu context has been created
1976 g_gpu_context_active = 1;
1977
1978 if(deviceCount <= 0) {
1979 return PyErr_Format(PyExc_EnvironmentError,
1980 "Can't use the GPU, no devices support CUDA");
1981 }
1982 if(card_number_provided && (card_nb < 0 || card_nb > (deviceCount - 1))) {
1983 return PyErr_Format(PyExc_ValueError,
1984 "Bad device number %d. Only %d devices available.",
1985 card_nb,
1986 deviceCount);
1987 }
1988
1989 cudaDeviceProp deviceProp;
1990 err = cudaGetDeviceProperties(&deviceProp, card_nb);
1991 if(cudaSuccess != err) {
1992 return PyErr_Format(PyExc_EnvironmentError,
1993 "Unable to get properties of gpu %i: %s",
1994 card_nb,
1995 cudaGetErrorString(cudaGetLastError()));
1996 }
1997
1998 if(deviceProp.major == 9999 && deviceProp.minor == 9999 ){
1999 return PyErr_Format(PyExc_EnvironmentError,
2000 "There is no device that supports CUDA");
2001 }
2002
2003 if(card_number_provided) {
2004 err = cudaSetDevice(card_nb);
2005 if(cudaSuccess != err) {
2006 return PyErr_Format(PyExc_EnvironmentError,
2007 "Unable to set device %i: %s",
2008 card_nb,
2009 cudaGetErrorString(cudaGetLastError()));
2010 }
2011 }
2012
2013 Py_INCREF(Py_None);
2014 return Py_None;
2015 }
2016
2017 PyObject *
2018 CudaNdarray_active_device_number(PyObject* _unused, PyObject* _unused_args) {
2019 // NB: No cuda error checking here; keeps things simple, and it's not
2020 // really necessary.
2021 int currentDevice;
2022 cudaGetDevice(&currentDevice);
2023 return PyInt_FromLong(currentDevice);
2024 }
2025
2026 PyObject *
2027 CudaNdarray_active_device_name(PyObject* _unused, PyObject* _unused_args) {
2028 // NB: No cuda error checking here; keeps things simple, and it's not
2029 // really necessary.
2030 int currentDevice;
2031 cudaGetDevice(&currentDevice);
2032
2033 cudaDeviceProp deviceProp;
2034 cudaGetDeviceProperties(&deviceProp, currentDevice);
2035 return PyString_FromString(deviceProp.name);
2036 }
2037
2038 PyObject *
2039 CudaNdarray_gpu_shutdown(PyObject* _unused, PyObject* _unused_args) {
2040 cudaThreadExit();
2041 g_gpu_context_active = 0; // context has now been closed down
2042 Py_INCREF(Py_None);
2043 return Py_None;
2044 }
2045
2046 /*
2047 * This function is tested in theano/misc/test_pycuda_theano_simple.py
2048 */
2049 PyObject *
2050 CudaNdarray_from_gpu_pointer(PyObject* _unused, PyObject* args)
2051 {
2052 PyObject *gpu_ptr = NULL;
2053 PyObject *shapes = NULL;
2054 PyObject *strides = NULL;
2055 PyObject *base = NULL;
2056 PyObject *rval = NULL;
2057
2058 //args should consist of 3 python objects
2059 //The first is the gpu ptr
2060 //The second if the shape
2061 //The third if the strides
2062 if (! PyArg_ParseTuple(args, "OOOO", &gpu_ptr, &shapes, &strides, &base))
2063 return NULL;
2064
2065 printf("In CudaNdarray_from_gpu_pointer\n");
2066 if (!PyLong_Check(gpu_ptr))
2067 {
2068 PyErr_Format(PyExc_Exception, "CudaNdarray_from_gpu_pointer: The gpu pointor is not an long");
2069 return NULL;
2070 }
2071
2072 Py_ssize_t nd = PyObject_Length(shapes);
2073 if (nd < 0)
2074 {
2075 PyErr_SetString(PyExc_TypeError, "CudaNdarray_from_gpu_pointer: Couldn't get length of second argument");
2076 return NULL;
2077 }
2078 Py_ssize_t nd_stride = PyObject_Length(strides);
2079 if (nd_stride < 0)
2080 {
2081 PyErr_SetString(PyExc_TypeError, "CudaNdarray_from_gpu_pointer: Couldn't get length of third argument");
2082 return NULL;
2083 }
2084
2085 if (nd != nd_stride)
2086 {
2087 PyErr_SetString(PyExc_TypeError, "CudaNdarray_from_gpu_pointer: We need the same number of shapes and strides");
2088 return NULL;
2089 }
2090
2091 rval = CudaNdarray_New();
2092
2093 if (CudaNdarray_set_nd((CudaNdarray *)rval, nd))
2094 {
2095 //CudaNdarray_set_nd set the error msg
2096 return NULL;
2097 }
2098 // set gpu pointeur
2099 assert(((CudaNdarray *)rval)->data_allocated == 0);
2100 if (CudaNdarray_set_device_data((CudaNdarray *)rval, (float *)PyInt_AsLong(gpu_ptr), base))
2101 {
2102 PyErr_SetString(PyExc_TypeError, "CudaNdarray_from_gpu_pointer: Error while setting the gpu pointor");
2103 return NULL;
2104
2105 }
2106
2107 // Set dims and strides
2108 for (int i = nd-1; i >= 0; --i)
2109 {
2110 PyObject * idx = PyLong_FromLong(i);
2111 if (idx == NULL)
2112 {
2113 PyErr_SetString(PyExc_Exception, "CudaNdarray_from_gpu_pointer: Couldn't make long object to loop over list/tuple");
2114 return NULL;
2115 }
2116 PyObject* dim_ = PyObject_GetItem(shapes, idx);
2117 PyObject* strd_ = PyObject_GetItem(strides, idx);
2118 if (!PyInt_Check(dim_))
2119 {
2120 PyErr_Format(PyExc_Exception, "CudaNdarray_from_gpu_pointer: shapes[%d] is not an int", i);
2121 return NULL;
2122 }
2123 if (!PyInt_Check(strd_))
2124 {
2125 PyErr_Format(PyExc_Exception, "CudaNdarray_from_gpu_pointer: strides[%d] is not an int", i);
2126 return NULL;
2127 }
2128 int dim = PyInt_AsLong(dim_);
2129 int strd = PyInt_AsLong(strd_);
2130 CudaNdarray_set_stride((CudaNdarray *)rval, i, strd);
2131 CudaNdarray_set_dim((CudaNdarray *)rval, i, dim);
2132 Py_DECREF(idx);
2133 Py_DECREF(dim_);
2134 Py_DECREF(strd_);
2135 }
2136 printf("CudaNdarray_from_gpu_pointer normal return\n");
2137 return rval;
2138 }
2139
2140 PyObject *
2141 CudaNdarray_Dot(PyObject* _unused, PyObject* args)
2142 {
2143 PyObject *l=NULL;
2144 PyObject *r=NULL;
2145 PyObject * rval = NULL;
2146
2147 //args should consist of two python objects ("OO")
2148 if (! PyArg_ParseTuple(args, "OO", &l, &r))
2149 return NULL;
2150
2151 if (!CudaNdarray_Check(l) || !CudaNdarray_Check(r))
2152 {
2153 PyErr_SetString(PyExc_TypeError, "CudaNdarray arguments required ");
2154 goto CudaNdarray_dot_fail;
2155 }
2156 if (((CudaNdarray*)l)->nd != 2)
2157 {
2158 PyErr_SetString(PyExc_TypeError, "need 2d CudaNdarray arg for now");
2159 goto CudaNdarray_dot_fail;
2160 }
2161 if (((CudaNdarray*)r)->nd != 2)
2162 {
2163 PyErr_SetString(PyExc_TypeError, "need 2d CudaNdarray arg for now");
2164 goto CudaNdarray_dot_fail;
2165 }
2166 rval = CudaNdarray_New();
2167 if (!rval)
2168 {
2169 goto CudaNdarray_dot_fail;
2170 }
2171 int dims[2];
2172 dims[0] = CudaNdarray_HOST_DIMS((CudaNdarray*)l)[0];
2173 dims[1] = CudaNdarray_HOST_DIMS((CudaNdarray*)r)[1];
2174 if (CudaNdarray_alloc_contiguous((CudaNdarray*)rval, 2, dims))
2175 {
2176 goto CudaNdarray_dot_fail;
2177 }
2178 if (CudaNdarray_gemm(1.0, (CudaNdarray*)l, (CudaNdarray*)r, 0.0, (CudaNdarray*)rval))
2179 {
2180 goto CudaNdarray_dot_fail;
2181 }
2182
2183 return rval;
2184
2185 CudaNdarray_dot_fail:
2186 Py_XDECREF(rval);
2187 return NULL;
2188 }
2189
2190 static PyObject *
2191 filter(PyObject* __unsed_self, PyObject *args) // args = (data, broadcastable, strict, storage)
2192 {
2193 /*
2194 * TODO: DOC what this function should do in the various cases of
2195 * What is 'strict' supposed to mean in the context of this function?
2196 * What do we do with input that could be interpreted as matching the broadcastable pattern in strict vs. non-strict cases?
2197 *
2198 */
2199 PyObject *py_data=NULL;
2200 PyArrayObject * data = NULL;
2201 int strict = 0;
2202 PyObject * broadcastable=NULL;
2203 PyObject * storage=NULL;
2204 CudaNdarray * rval=NULL;
2205
2206 //Python object references which are provided to the caller are borrowed references
2207 if (!PyArg_ParseTuple(args, "OOiO", &py_data, &broadcastable, &strict, &storage)) return NULL;
2208
2209 if (!PyTuple_Check(broadcastable)){
2210 PyErr_SetString(PyExc_TypeError, "broadcastable arg should be a tuple of int.");
2211 return NULL;
2212 }
2213 Py_INCREF(py_data);
2214 Py_INCREF(broadcastable);
2215
2216 CudaNdarray * cnda = (CudaNdarray*)py_data;
2217
2218 if (strict || CudaNdarray_Check(py_data))
2219 {
2220 //TODO: support non-strict "casting" from a vt to the broadcastable/type/size that we need.
2221 if (!CudaNdarray_Check(py_data))
2222 {
2223 Py_DECREF(py_data);
2224 Py_DECREF(broadcastable);
2225 PyErr_SetString(PyExc_TypeError, "strict mode requires CudaNdarray");
2226 return NULL;
2227 }
2228 if (cnda->nd != PyTuple_Size(broadcastable))
2229 {
2230 Py_DECREF(py_data);
2231 Py_DECREF(broadcastable);
2232 PyErr_Format(PyExc_TypeError, "Wrong rank: %i vs %li", cnda->nd, (long)PyTuple_Size(broadcastable));
2233 return NULL;
2234 }
2235 for (int i = 0; i < cnda->nd; ++i)
2236 {
2237 if ((CudaNdarray_HOST_DIMS(cnda)[i] > 1) && PyInt_AsLong(PyTuple_GetItem(broadcastable, Py_ssize_t(i))))
2238 {
2239 PyErr_Format(PyExc_TypeError, "Non-unit size in broadcastable vt dimension %i", i);
2240 Py_DECREF(py_data);
2241 Py_DECREF(broadcastable);
2242 return NULL;
2243 }
2244 }
2245 Py_DECREF(broadcastable);
2246 return py_data;
2247 }
2248 else
2249 {
2250 data = (PyArrayObject*)PyArray_FromObject(py_data, REAL_TYPENUM, PyTuple_Size(broadcastable), PyTuple_Size(broadcastable));
2251 if (!data)
2252 {
2253 //err message already defined
2254 Py_DECREF(py_data);
2255 Py_DECREF(broadcastable);
2256 return NULL;
2257 }
2258 for (int i = 0; i < data->nd; ++i)
2259 {
2260 if ((data->dimensions[i] > 1) && PyInt_AsLong(PyTuple_GetItem(broadcastable, Py_ssize_t(i))))
2261 {
2262 PyErr_Format(PyExc_TypeError, "Non-unit size in broadcastable dimension %i", i);
2263 Py_DECREF(data);
2264 Py_DECREF(py_data);
2265 Py_DECREF(broadcastable);
2266 return NULL;
2267 }
2268 }
2269 if (storage && CudaNdarray_Check(storage))
2270 {
2271 rval = (CudaNdarray*) storage;
2272 Py_INCREF(rval);
2273 }
2274 else
2275 {
2276 rval = (CudaNdarray*) CudaNdarray_New();
2277 }
2278 if (rval)
2279 {
2280 if (CudaNdarray_CopyFromArray(rval, data))
2281 {
2282 Py_DECREF(rval);
2283 rval = NULL;
2284 }
2285 }
2286 Py_DECREF(data);
2287 Py_DECREF(py_data);
2288 Py_DECREF(broadcastable);
2289 return (PyObject*)rval;
2290 }
2291 }
2292
2293 //TODO-- CudaNdarray_Dot and CudaNdarray_active_device_name are following different capitalization conventions.
2294 // Pick one and standardize it, this file is already annoying enough to grep through
2295 static PyMethodDef module_methods[] = {
2296 {"dimshuffle", CudaNdarray_Dimshuffle, METH_VARARGS, "Returns the dimshuffle of a CudaNdarray."},
2297 {"dot", CudaNdarray_Dot, METH_VARARGS, "Returns the matrix product of two CudaNdarray arguments."},
2298 {"gpu_init", CudaNdarray_gpu_init, METH_VARARGS, "Select the gpu card to use; also usable to test whether CUDA is available."},
2299 {"active_device_name", CudaNdarray_active_device_name, METH_VARARGS, "Get the name of the active device."},
2300 {"active_device_number", CudaNdarray_active_device_number, METH_VARARGS, "Get the number of the active device."},
2301 {"gpu_shutdown", CudaNdarray_gpu_shutdown, METH_VARARGS, "Shut down the gpu."},
2302 {"ptr_int_size", CudaNdarray_ptr_int_size, METH_VARARGS, "Return a tuple with the size of gpu pointer, cpu pointer and int in bytes."},
2303 {"filter", filter, METH_VARARGS, "filter(obj, broadcastable, strict, storage) returns a CudaNdarray initialized to obj if it matches the constraints of broadcastable. strict=True prevents any numeric casting. If storage is a CudaNdarray it may be overwritten and used as the return value."},
2304 {"outstanding_mallocs", outstanding_mallocs, METH_VARARGS, "how many more mallocs have been called than free's"},
2305 {"from_gpu_pointer", CudaNdarray_from_gpu_pointer, METH_VARARGS, "Used to create a CudaNdarray from already allocated memory on the gpu.(example by pycuda)"},
2306 {NULL, NULL, NULL, NULL} /* Sentinel */
2307 };
2308
2309 #ifndef PyMODINIT_FUNC /* declarations for DLL import/export */
2310 #define PyMODINIT_FUNC void
2311 #endif
2312 PyMODINIT_FUNC
2313 initcuda_ndarray(void)
2314 {
2315 import_array();
2316
2317 PyObject* m;
2318
2319 if (PyType_Ready(&CudaNdarrayType) < 0)
2320 return;
2321
2322 m = Py_InitModule3("cuda_ndarray", module_methods,
2323 "Example module that creates an extension type.");
2324
2325 if (m == NULL)
2326 return;
2327
2328 Py_INCREF(&CudaNdarrayType);
2329 PyModule_AddObject(m, "CudaNdarray", (PyObject *)&CudaNdarrayType);
2330 #if COMPUTE_GPU_MEM_USED
2331 for(int i=0;i<TABLE_SIZE;i++){
2332 _alloc_size_table[i].ptr=NULL;
2333 _alloc_size_table[i].size=0;
2334 }
2335 #endif
2336 // cublasInit();
2337 //if (0&&CUBLAS_STATUS_SUCCESS != cublasGetError())
2338 //{
2339 //std::cerr << "WARNING: initcuda_ndarray: error initializing device\n";
2340 //}
2341 if (0) //TODO: is this necessary?
2342 {
2343 int deviceId = 0; // TODO: what number goes here?
2344 cudaSetDevice(deviceId);
2345 cudaError_t err = cudaGetLastError();
2346 if( cudaSuccess != err)
2347 {
2348 std::cerr << "Error in SetDevice:" << cudaGetErrorString(err) << "\n";
2349 }
2350 }
2351 }
2352
2353
2354 //////////////////////////////////////
2355 //
2356 // C API FOR CudaNdarray
2357 //
2358 //////////////////////////////////////
2359
2360 int
2361 CudaNdarray_Check(const PyObject * ob)
2362 {
2363 //TODO: doesn't work with inheritance
2364 return CudaNdarray_CheckExact(ob);
2365 }
2366 int
2367 CudaNdarray_CheckExact(const PyObject * ob)
2368 {
2369 return ((ob->ob_type == &CudaNdarrayType) ? 1 : 0);
2370 }
2371
2372 PyObject *
2373 CudaNdarray_New(int nd)
2374 {
2375 CudaNdarray *self = (CudaNdarray *)CudaNdarrayType.tp_alloc(&CudaNdarrayType, 0);
2376 if (self == NULL)
2377 {
2378 PyErr_SetString(PyExc_RuntimeError, "CudaNdarray_New failed to allocate self");
2379 return NULL;
2380 }
2381 CudaNdarray_null_init(self);
2382
2383 if (nd == 0)
2384 {
2385 self->nd = 0;
2386 }
2387 else if (nd > 0)
2388 {
2389 if (CudaNdarray_set_nd(self, nd))
2390 {
2391 Py_DECREF(self);
2392 return NULL;
2393 }
2394 }
2395 ++_outstanding_mallocs[1];
2396 return (PyObject *)self;
2397 }
2398
2399
2400
2401 //////////////////////////////
2402 //
2403 // Published helper functions
2404 //
2405 //////////////////////////////
2406
2407 int
2408 cublas_init()
2409 {
2410 cublasInit();
2411 if (CUBLAS_STATUS_SUCCESS != cublasGetError())
2412 {
2413 PyErr_SetString(PyExc_RuntimeError, "error initializing device");
2414 return -1;
2415 }
2416 return 0;
2417 }
2418 int
2419 cublas_shutdown()
2420 {
2421 cublasShutdown();
2422 if (CUBLAS_STATUS_SUCCESS != cublasGetError())
2423 {
2424 PyErr_SetString(PyExc_RuntimeError, "error shutting down device");
2425 return -1;
2426 }
2427 return 0;
2428 }
2429
2430 int
2431 CudaNdarray_CopyFromArray(CudaNdarray * self, PyArrayObject*obj)
2432 {
2433 int err = CudaNdarray_alloc_contiguous(self, obj->nd, obj->dimensions);
2434 if (err) {
2435 return err;
2436 }
2437
2438 int typenum = PyArray_TYPE(obj);
2439 if (typenum != REAL_TYPENUM)
2440 {
2441 PyErr_SetString(PyExc_TypeError, "can only copy from float arrays");
2442 return -1;
2443 }
2444 assert( 4 == PyArray_ITEMSIZE(obj));
2445 PyObject * py_src = PyArray_ContiguousFromAny((PyObject*)obj, typenum, self->nd, self->nd);
2446 if (!py_src) {
2447 return -1;
2448 }
2449 cublasSetVector(PyArray_SIZE(py_src),
2450 sizeof(real),
2451 PyArray_DATA(py_src), 1,
2452 self->devdata, 1);
2453 CNDA_THREAD_SYNC;
2454 if (CUBLAS_STATUS_SUCCESS != cublasGetError())
2455 {
2456 PyErr_SetString(PyExc_RuntimeError, "error copying data to device memory");
2457 Py_DECREF(py_src);
2458 return -1;
2459 }
2460 Py_DECREF(py_src);
2461 return 0;
2462 }
2463 bool
2464 CudaNdarray_is_c_contiguous(const CudaNdarray * self)
2465 {
2466 bool c_contiguous = true;
2467 int size = 1;
2468 for (int i = self->nd-1; (i >= 0) && c_contiguous; --i)
2469 {
2470 if (CudaNdarray_HOST_DIMS(self)[i] == 1)
2471 continue;
2472 if (CudaNdarray_HOST_STRIDES(self)[i] != size)
2473 {
2474 c_contiguous = false;
2475 }
2476 size = size * CudaNdarray_HOST_DIMS(self)[i];
2477 }
2478 return c_contiguous;
2479 }
2480
2481 PyObject *
2482 CudaNdarray_new_nd(int nd)
2483 {
2484 CudaNdarray * rval = (CudaNdarray*) CudaNdarray_New();
2485 if (!rval || CudaNdarray_set_nd(rval, nd))
2486 {
2487 Py_XDECREF(rval);
2488 rval = NULL;
2489 }
2490 return (PyObject *) rval;
2491 }
2492
2493
2494 /**
2495 * Initialize 'self' as a view of 'base', with memory storage 'data'
2496 */
2497
2498 int CudaNdarray_set_device_data(CudaNdarray * self, float * data, PyObject * base)
2499 {
2500 if (self->data_allocated)
2501 {
2502 assert(self->devdata);
2503 if (device_free(self->devdata))
2504 {
2505 self->devdata = NULL;
2506 self->data_allocated = 0;
2507 return -1;
2508 }
2509 }
2510 // Get the original base object (base.base.base...)
2511 // TODO: check that base is indeed a CudaNdarray?
2512 PyObject * orig_base = base;
2513 while (((CudaNdarray*) orig_base)->base)
2514 {
2515 // base_base is itself a view
2516 orig_base = ((CudaNdarray*) orig_base)->base;
2517 }
2518 //N.B. XDECREF and XINCREF are no-ops for NULL pointers
2519 if (self->base != orig_base)
2520 {
2521 Py_XDECREF(self->base);
2522 self->base = orig_base;
2523 Py_XINCREF(self->base);
2524 }
2525 self->data_allocated = 0;
2526 self->devdata = data;
2527 return 0;
2528 }
2529
2530 static __global__ void k_copy_1d(const int N, const float * x, const int sx, float * y, const int sy)
2531 {
2532 for (int i = threadIdx.x + blockIdx.x * blockDim.x; i < N; i += gridDim.x*blockDim.x)
2533 {
2534 y[i*sy] = x[i*sx];
2535 }
2536 }
2537
2538 //copy from other into self
2539 int CudaNdarray_CopyFromCudaNdarray(CudaNdarray * self, CudaNdarray * other, bool unbroadcast)
2540 {
2541 int verbose = 0;
2542 if (verbose>1) fprintf(stderr, "CudaNdarray_CopyFromCudaNdarray\n");
2543
2544 //standard elemwise size checks
2545 if (self->nd == -1)
2546 {
2547 PyErr_SetString(PyExc_TypeError, "can't copy into un-initialized CudaNdarray");
2548 return -1;
2549 }
2550 if (self->nd != other->nd)
2551 {
2552 PyErr_Format(PyExc_NotImplementedError, "CudaNdarray_CopyFromCudaNdarray: need same number of dims. destination nd=%d, source nd=%d. No broadcasting implemented.", self->nd, other->nd);
2553 return -1;
2554 }
2555 //standard elemwise dim checks (also compute total size)
2556 unsigned int size = 1;
2557 unsigned int size_source = 1;
2558 for (int i = 0; i< self->nd; ++i)
2559 {
2560 if ((CudaNdarray_HOST_DIMS(self)[i] != CudaNdarray_HOST_DIMS(other)[i])
2561 && (1!=CudaNdarray_HOST_DIMS(other)[i] || !unbroadcast) )
2562 {
2563 PyErr_Format(PyExc_ValueError, "need same dimensions for dim %d, destination=%d, source=%d",
2564 i, CudaNdarray_HOST_DIMS(self)[i], CudaNdarray_HOST_DIMS(other)[i]);
2565 return -1;
2566 }
2567 size *= (unsigned int) CudaNdarray_HOST_DIMS(self)[i];
2568 size_source *= (unsigned int) CudaNdarray_HOST_DIMS(other)[i];
2569 }
2570 if (0 == size)
2571 {
2572 return 0; //nothing to copy, we're done.
2573 }
2574 if (CudaNdarray_is_c_contiguous(self) && CudaNdarray_is_c_contiguous(other) && size == size_source)
2575 {
2576 cublasScopy(size, CudaNdarray_DEV_DATA(other), 1, CudaNdarray_DEV_DATA(self), 1);
2577 CNDA_THREAD_SYNC;
2578 if (CUBLAS_STATUS_SUCCESS != cublasGetError())
2579 {
2580 PyErr_SetString(PyExc_RuntimeError, "Error copying memory");
2581 return -1;
2582 }
2583 return 0;
2584 }
2585 //TODO: rewrite these copy operations to be more efficient
2586 // See, for example the transpose example in the cuda_sdk.
2587 switch (self->nd)
2588 {
2589 case 0: // scalar
2590 {
2591 // THIS CASE SHOULD NEVER HAPPEN BECAUSE SCALARS ARE ALWAYS C CONTIGUOUS
2592 assert(0);
2593 }; break;
2594 case 1: // vector
2595 {
2596 if (verbose) fprintf(stderr, "Copying non-contiguous vector\n");
2597 if (verbose) fprint_CudaNdarray(stderr, other);
2598 unsigned int n_blocks = std::min(size, (unsigned int)NUM_VECTOR_OP_BLOCKS);
2599 unsigned int n_threads = std::min(ceil_intdiv(size, n_blocks), (unsigned int)NUM_VECTOR_OP_THREADS_PER_BLOCK);
2600 k_copy_1d<<<n_blocks, n_threads>>>(size,
2601 CudaNdarray_DEV_DATA(other), CudaNdarray_HOST_STRIDES(other)[0],
2602 CudaNdarray_DEV_DATA(self), CudaNdarray_HOST_STRIDES(self)[0]);
2603 CNDA_THREAD_SYNC;
2604 cudaError_t err = cudaGetLastError();
2605 if( cudaSuccess != err)
2606 {
2607 PyErr_Format(PyExc_RuntimeError, "Cuda error: %s: %s. (n_blocks=%i, n_threads_per_block=%i)\n", "k_copy_1d", cudaGetErrorString(err), n_blocks, n_threads);
2608 return -1;
2609 }
2610 }; break;
2611 default:
2612 {
2613 assert (cudaSuccess == cudaGetLastError());
2614 if (verbose) fprintf(stderr, "Copying with default version unbroadcast=%d\n", unbroadcast);
2615 // call worker routine
2616 unsigned int n_blocks = std::min(size, (unsigned int)NUM_VECTOR_OP_BLOCKS);
2617 unsigned int threads_per_block = std::min(ceil_intdiv(size, n_blocks), (unsigned int)NUM_VECTOR_OP_THREADS_PER_BLOCK);
2618 CudaNdarray * cuda_dims = other;
2619 if(unbroadcast)
2620 cuda_dims = self;
2621 //copy from other into self
2622 k_elemwise_unary_rowmajor_copy<<<n_blocks, threads_per_block>>>(
2623 size,
2624 (unsigned int)other->nd,
2625 (const int *)CudaNdarray_DEV_DIMS(cuda_dims),
2626 (const float*)CudaNdarray_DEV_DATA(other), (const int *)CudaNdarray_DEV_STRIDES(other),
2627 CudaNdarray_DEV_DATA(self), (const int *)CudaNdarray_DEV_STRIDES(self));
2628 CNDA_THREAD_SYNC;
2629 cudaError_t err = cudaGetLastError();
2630 if( cudaSuccess != err)
2631 {
2632 //fprint_CudaNdarray(stderr, self);
2633 //fprint_CudaNdarray(stderr, other);
2634 PyErr_Format(PyExc_RuntimeError, "Cuda error: %s: %s. (n_blocks=%i, n_threads_per_block=%i)\n", "k_elemwise_unary_rowmajor_copy", cudaGetErrorString(err), n_blocks, threads_per_block);
2635 return -1;
2636 }
2637 }
2638 };
2639 return 0;
2640 }
2641
2642 int CudaNdarray_gemm(float alpha, const CudaNdarray * A, const CudaNdarray * B, float beta, CudaNdarray * C)
2643 {
2644 if (A->nd != 2) { PyErr_SetString(PyExc_ValueError, "non-matrix arg to gemm"); return -1; }
2645 if (B->nd != 2) { PyErr_SetString(PyExc_ValueError, "non-matrix arg to gemm"); return -1; }
2646 if (C->nd != 2) { PyErr_SetString(PyExc_ValueError, "non-matrix arg to gemm"); return -1; }
2647
2648 if ((CudaNdarray_HOST_DIMS(A)[1] != CudaNdarray_HOST_DIMS(B)[0])
2649 || (CudaNdarray_HOST_DIMS(A)[0] != CudaNdarray_HOST_DIMS(C)[0])
2650 || (CudaNdarray_HOST_DIMS(B)[1] != CudaNdarray_HOST_DIMS(C)[1]))
2651 {
2652 PyErr_Format(PyExc_ValueError, "dimension mismatch in args to gemm (%i,%i)x(%i,%i)->(%i,%i)",
2653 CudaNdarray_HOST_DIMS(A)[0],
2654 CudaNdarray_HOST_DIMS(A)[1],
2655 CudaNdarray_HOST_DIMS(B)[0],
2656 CudaNdarray_HOST_DIMS(B)[1],
2657 CudaNdarray_HOST_DIMS(C)[0],
2658 CudaNdarray_HOST_DIMS(C)[1]);
2659 return -1;
2660 }
2661
2662 // a matrix has non-unit size and non-unit stride in both directions, we can't operate in-place
2663 // TODO: make a copy instead of returning in error
2664 if (((CudaNdarray_HOST_DIMS(A)[0] > 1) && (CudaNdarray_HOST_STRIDES(A)[0] != 1)) && ((CudaNdarray_HOST_DIMS(A)[1] > 1) && (CudaNdarray_HOST_STRIDES(A)[1] != 1)))
2665 { PyErr_SetString(PyExc_NotImplementedError, "non-unit stride in gemm arg"); return -1; }
2666 if (((CudaNdarray_HOST_DIMS(B)[0] > 1) && (CudaNdarray_HOST_STRIDES(B)[0] != 1)) && ((CudaNdarray_HOST_DIMS(B)[1] > 1) && (CudaNdarray_HOST_STRIDES(B)[1] != 1)))
2667 { PyErr_SetString(PyExc_NotImplementedError, "non-unit stride in gemm arg"); return -1; }
2668 if (((CudaNdarray_HOST_DIMS(C)[0] > 1) && (CudaNdarray_HOST_STRIDES(C)[0] != 1)) && ((CudaNdarray_HOST_DIMS(C)[1] > 1) && (CudaNdarray_HOST_STRIDES(C)[1] != 1)))
2669 { PyErr_SetString(PyExc_NotImplementedError, "non-unit stride in gemm arg"); return -1; }
2670
2671 // the unit integer is divided logically into three fields of 4 bits
2672 // the lowermost 4 bits encode the stride pattern of the output
2673 // the next higher 4 bits encode the B variable (or y)
2674 // the next higher 4 bits encode the C variable (or x)
2675 //
2676 // the stride pattern for each input is encoded as 0 for unit stride from col to col (Row major)
2677 // 1 for unit stride from row to row (Col major)
2678
2679 // a stride of 0 implies a dimension of 1 - so we can actually define
2680 // a stride of 0 as a 'unit' stride because gemm will never use it.
2681 int unit = 0;
2682 if (CudaNdarray_HOST_STRIDES(A)[1] == 1 || CudaNdarray_HOST_STRIDES(A)[1] == 0) {
2683 unit |= (0x0 << 8);
2684 } else if (CudaNdarray_HOST_STRIDES(A)[0] == 1 || CudaNdarray_HOST_STRIDES(A)[0] == 0) {
2685 unit |= (0x1 << 8);
2686 } else {
2687 unit |= (0x2 << 8);
2688 }
2689 if (CudaNdarray_HOST_STRIDES(B)[1] == 1 || CudaNdarray_HOST_STRIDES(B)[1] == 0) {
2690 unit |= (0x0 << 4);
2691 } else if (CudaNdarray_HOST_STRIDES(B)[0] == 1 || CudaNdarray_HOST_STRIDES(B)[0] == 0) {
2692 unit |= (0x1 << 4);
2693 } else {
2694 unit |= (0x2 << 4);
2695 }
2696 if (CudaNdarray_HOST_STRIDES(C)[1] == 1 || CudaNdarray_HOST_STRIDES(C)[1] == 0) {
2697 unit |= (0x0 << 0);
2698 } else if (CudaNdarray_HOST_STRIDES(C)[0] == 1 || CudaNdarray_HOST_STRIDES(C)[0] == 0) {
2699 unit |= (0x1 << 0);
2700 } else {
2701 unit |= (0x2 << 0);
2702 }
2703
2704 // I don't know if cudablas handles negative strides
2705 assert (CudaNdarray_HOST_STRIDES(A)[0] >= 0) ; // for now
2706 assert (CudaNdarray_HOST_STRIDES(A)[1] >= 0) ; // for now
2707 assert (CudaNdarray_HOST_STRIDES(B)[0] >= 0) ; // for now
2708 assert (CudaNdarray_HOST_STRIDES(B)[1] >= 0) ; // for now
2709 assert (CudaNdarray_HOST_STRIDES(C)[0] >= 0) ; // for now
2710 assert (CudaNdarray_HOST_STRIDES(C)[1] >= 0) ; // for now
2711
2712 /* create appropriate strides for malformed matrices that are row or column
2713 * vectors
2714 */
2715 int sa_0 = (CudaNdarray_HOST_DIMS(A)[0] > 1) ? CudaNdarray_HOST_STRIDES(A)[0] : CudaNdarray_HOST_DIMS(A)[1];
2716 int sa_1 = (CudaNdarray_HOST_DIMS(A)[1] > 1) ? CudaNdarray_HOST_STRIDES(A)[1] : CudaNdarray_HOST_DIMS(A)[0];
2717 int sb_0 = (CudaNdarray_HOST_DIMS(B)[0] > 1) ? CudaNdarray_HOST_STRIDES(B)[0] : CudaNdarray_HOST_DIMS(B)[1];
2718 int sb_1 = (CudaNdarray_HOST_DIMS(B)[1] > 1) ? CudaNdarray_HOST_STRIDES(B)[1] : CudaNdarray_HOST_DIMS(B)[0];
2719 int sc_0 = (CudaNdarray_HOST_DIMS(C)[0] > 1) ? CudaNdarray_HOST_STRIDES(C)[0] : CudaNdarray_HOST_DIMS(C)[1];
2720 int sc_1 = (CudaNdarray_HOST_DIMS(C)[1] > 1) ? CudaNdarray_HOST_STRIDES(C)[1] : CudaNdarray_HOST_DIMS(C)[0];
2721
2722 float* a = CudaNdarray_DEV_DATA(A);
2723 float* b = CudaNdarray_DEV_DATA(B);
2724 float* c = CudaNdarray_DEV_DATA(C);
2725 char N = 'N';
2726 char T = 'T';
2727 //std::cerr << (unit/256) MOD 16 << (unit / 16) MOD 16 << unit MOD 16<< '\\n';
2728 //TODO: recognize the negative stride and make a copy of the offending argument,
2729 //rather than aborting
2730 #define CHK_STRIDE_SGEMM(T0, T1, D0, D1, D2, a, x, sx, y, sy, b, z, sz) \
2731 if ((sx > 0) && (sy > 0) && (sz > 0)) { \
2732 cublasSgemm(T0, T1, D0, D1, D2, a, x, sx, y, sy, b, z, sz); \
2733 } else { \
2734 PyErr_SetString(PyExc_NotImplementedError, "negative stride to sGemm");\
2735 return -1; \
2736 }
2737
2738 switch(unit)
2739 {
2740 case 0x000: CHK_STRIDE_SGEMM(N, N, CudaNdarray_HOST_DIMS(C)[1], CudaNdarray_HOST_DIMS(C)[0], CudaNdarray_HOST_DIMS(A)[1], alpha, b, sb_0, a, sa_0, beta, c, sc_0); break;
2741 case 0x100: CHK_STRIDE_SGEMM(N, T, CudaNdarray_HOST_DIMS(C)[1], CudaNdarray_HOST_DIMS(C)[0], CudaNdarray_HOST_DIMS(A)[1], alpha, b, sb_0, a, sa_1, beta, c, sc_0); break;
2742 case 0x010: CHK_STRIDE_SGEMM(T, N, CudaNdarray_HOST_DIMS(C)[1], CudaNdarray_HOST_DIMS(C)[0], CudaNdarray_HOST_DIMS(A)[1], alpha, b, sb_1, a, sa_0, beta, c, sc_0); break;
2743 case 0x110: CHK_STRIDE_SGEMM(T, T, CudaNdarray_HOST_DIMS(C)[1], CudaNdarray_HOST_DIMS(C)[0], CudaNdarray_HOST_DIMS(A)[1], alpha, b, sb_1, a, sa_1, beta, c, sc_0); break;
2744 case 0x001: CHK_STRIDE_SGEMM(T, T, CudaNdarray_HOST_DIMS(C)[0], CudaNdarray_HOST_DIMS(C)[1], CudaNdarray_HOST_DIMS(A)[1], alpha, a, sa_0, b, sb_0, beta, c, sc_1); break;
2745 case 0x101: CHK_STRIDE_SGEMM(N, T, CudaNdarray_HOST_DIMS(C)[0], CudaNdarray_HOST_DIMS(C)[1], CudaNdarray_HOST_DIMS(A)[1], alpha, a, sa_1, b, sb_0, beta, c, sc_1); break;
2746 case 0x011: CHK_STRIDE_SGEMM(T, N, CudaNdarray_HOST_DIMS(C)[0], CudaNdarray_HOST_DIMS(C)[1], CudaNdarray_HOST_DIMS(A)[1], alpha, a, sa_0, b, sb_1, beta, c, sc_1); break;
2747 case 0x111: CHK_STRIDE_SGEMM(N, N, CudaNdarray_HOST_DIMS(C)[0], CudaNdarray_HOST_DIMS(C)[1], CudaNdarray_HOST_DIMS(A)[1], alpha, a, sa_1, b, sb_1, beta, c, sc_1); break;
2748 default: PyErr_Format(PyExc_ValueError, "some matrix has no unit stride (unit=%i)", unit);
2749 return -1;
2750 };
2751 CNDA_THREAD_SYNC;
2752 cudaError_t err = cudaGetLastError();
2753 if (CUBLAS_STATUS_SUCCESS != err)
2754 {
2755 PyErr_Format(PyExc_RuntimeError, "cublassGemm failed (%s)",cudaGetErrorString(err));
2756 return -1;
2757 }
2758 return 0;
2759 }
2760
2761 int CudaNdarray_sger(float alpha, CudaNdarray * x, CudaNdarray * y, CudaNdarray * A) {
2762 if (x->nd != 1) { PyErr_SetString(PyExc_ValueError, "non-vector arg x to sger"); return -1; }
2763 if (y->nd != 1) { PyErr_SetString(PyExc_ValueError, "non-vector arg y to sger"); return -1; }
2764 if (A->nd != 2) { PyErr_SetString(PyExc_ValueError, "non-matrix arg A to sger"); return -1; }
2765
2766 if ((CudaNdarray_HOST_DIMS(A)[0] != CudaNdarray_HOST_DIMS(x)[0])
2767 || (CudaNdarray_HOST_DIMS(A)[1] != CudaNdarray_HOST_DIMS(y)[0])) {
2768 PyErr_Format(PyExc_ValueError,
2769 "dimension mismatch in args to sger (%i)x(%i)->(%i,%i)",
2770 CudaNdarray_HOST_DIMS(x)[0],
2771 CudaNdarray_HOST_DIMS(y)[0],
2772 CudaNdarray_HOST_DIMS(A)[0],
2773 CudaNdarray_HOST_DIMS(A)[1]);
2774 return -1;
2775 }
2776
2777 // Maybe this could work, but be safe for now
2778 if (!CudaNdarray_is_c_contiguous(A)) {
2779 PyErr_SetString(PyExc_NotImplementedError, "non-c continugous A in sger");
2780 return -1;
2781 }
2782
2783 // Same for this, be safe
2784 assert (CudaNdarray_HOST_STRIDES(x)[0] >= 0);
2785 assert (CudaNdarray_HOST_STRIDES(y)[0] >= 0);
2786
2787 // Since Sger expects A in col-major, we invert x and y to fake this.
2788 cublasSger(CudaNdarray_HOST_DIMS(y)[0], CudaNdarray_HOST_DIMS(x)[0], alpha,
2789 CudaNdarray_DEV_DATA(y), CudaNdarray_HOST_STRIDES(y)[0],
2790 CudaNdarray_DEV_DATA(x), CudaNdarray_HOST_STRIDES(x)[0],
2791 CudaNdarray_DEV_DATA(A), CudaNdarray_HOST_DIMS(A)[1]);
2792 CNDA_THREAD_SYNC;
2793
2794 cudaError_t err = cudaGetLastError();
2795 if (CUBLAS_STATUS_SUCCESS != err)
2796 {
2797 PyErr_Format(PyExc_RuntimeError, "cublasSger failed (%s)",cudaGetErrorString(err));
2798 return -1;
2799 }
2800 return 0;
2801 }
2802
2803 /**
2804 *
2805 * Precondition:
2806 * a->dim[d] == (dims_a[d]==0) ? (1 << log2_dims_a[d]) : dims_a[d]
2807 * z->dim[d] == (z_str[d]==0) ? 1 : dims_a[d];
2808 *
2809 * TODO: templatize this function to support other reductions.
2810 * All that needs to change is the initial value for sum, and the reduction operator.
2811 */
2812
2813 static __global__ void kernel_reduce_sum(const unsigned int size_z,
2814 const unsigned int nd,
2815 const int * dims_a,
2816 const int * log2_dims_a,
2817 const int * a_str,
2818 const float * a_data,
2819 const int * z_str,
2820 float * z_data)
2821 {
2822 const unsigned int idx = blockIdx.x * blockDim.x + threadIdx.x;
2823 const unsigned int numThreads = blockDim.x * gridDim.x;
2824
2825 //structure data contains the strides and dimensions of both a and z
2826 // a_dim[0], a_dim[1], ... a_dim[nd-1],
2827 // a_log2dim[0], a_log2dim[1], ... a_log2dim[nd-1],
2828 // a_str[0], ... a_str[nd-1],
2829 // z_str[0], ... z_str[nd-1]
2830 extern __shared__ int structure_data[];
2831 for (unsigned int i = threadIdx.x; i < nd; i += blockDim.x)
2832 {
2833 structure_data[i+0*nd] = dims_a[i];
2834 structure_data[i+1*nd] = log2_dims_a[i];
2835 structure_data[i+2*nd] = a_str[i];
2836 structure_data[i+3*nd] = z_str[i];
2837 }
2838 dims_a = structure_data;
2839 log2_dims_a = structure_data + nd;
2840 a_str = structure_data + 2*nd;
2841 z_str = structure_data + 3*nd;
2842
2843 __syncthreads(); //wait for all the shared structure to be loaded
2844
2845 for (unsigned int i = idx; i < size_z; i += numThreads)
2846 {
2847 unsigned int ii = i;
2848 const float * a_data_i = a_data;
2849 float * z_data_i = z_data;
2850 unsigned int n_reduce_elements = 1;
2851 unsigned int n_reduce_dims = 0;
2852 unsigned int reduce_dim0 = nd-1;
2853
2854
2855 //In this loop, we locate the initial element of the slice that we'd like to reduce with this thread
2856 // At the same time, we [re]calculate the size of that slice (n_reduce_elements)
2857 for (unsigned int d = 0; d < nd; ++d)
2858 {
2859 if (a_str[d] && (!z_str[d])) // this means 'd' is a dimension we are reducing over
2860 {
2861 n_reduce_elements *= dims_a[d];
2862 n_reduce_dims += 1;
2863 reduce_dim0 = (d < reduce_dim0) ? d : reduce_dim0;
2864 }
2865 else //'d' is not a dimension that we are reducing over
2866 {
2867 unsigned int pos_d;
2868 if (log2_dims_a[d]==-1) //TODO: when things are working, use this switch
2869 {
2870 // this branch is not preferred,
2871 // because the manual said that integer mod and div operations are slow on gpu
2872 pos_d = (ii % dims_a[d]);
2873 ii = (ii / dims_a[d]);
2874 }
2875 else
2876 {
2877 pos_d = (ii & ((1 << log2_dims_a[d])-1)); //take the lower log2_dims bits
2878 ii = (ii >> log2_dims_a[d]); //shift those lower log2_dims bits off of ii
2879 }
2880 a_data_i += pos_d * a_str[d];
2881 z_data_i += pos_d * z_str[d];
2882 }
2883 }
2884 // now we've got pointers a_data_i and z_data_i into element 0 of the slice over which we are reducing
2885 // do a similar loop
2886
2887 float sum = 0.0f;
2888 switch(n_reduce_dims)
2889 {
2890 case 0:
2891 {
2892 sum = a_data_i[0];
2893 }
2894 break;
2895 case 1:
2896 {
2897 const int stride = a_str[reduce_dim0];
2898 const float * a_data_i_max = a_data_i + dims_a[reduce_dim0] * stride;
2899 while (a_data_i != a_data_i_max)
2900 {
2901 sum += a_data_i[0];
2902 a_data_i += stride;
2903 }
2904 }
2905 break;
2906 case 2:
2907 {
2908 int rd = reduce_dim0+1;
2909 for (; rd < nd; ++rd)
2910 {
2911 if (a_str[rd] && (!z_str[rd])) // this means 'rd' is a dimension we are reducing over
2912 break;
2913 }
2914 const int stride0 = a_str[reduce_dim0];
2915 const int stride1 = a_str[rd];
2916 for (int ii = 0; ii < dims_a[rd]; ++ii)
2917 {
2918 const float * a_data_ri = a_data_i + ii * stride1;
2919 const float * a_data_ri_max = a_data_ri + dims_a[reduce_dim0] * stride0;
2920 while (a_data_ri != a_data_ri_max)
2921 {
2922 sum += a_data_ri[0];
2923 a_data_ri += stride0;
2924 }
2925 }
2926 };
2927 break;
2928 default:
2929 {
2930 for (unsigned int reduce_i = 0; reduce_i < n_reduce_elements; ++reduce_i)
2931 {
2932 //TODO: optimize this loop to work more like theano's Elemwise. It's serial code.
2933 unsigned int reduce_ii = reduce_i;
2934 const float * a_data_ri = a_data_i;
2935
2936 //This loop finds the element in the a slice to add.
2937 for (unsigned int rd = reduce_dim0; rd < nd; ++rd)
2938 {
2939 unsigned int pos_d;
2940 if (a_str[rd] && (!z_str[rd])) // this means 'd' is a dimension we are reducing over
2941 {
2942 if (log2_dims_a[rd]==-1)
2943 {
2944 // this branch is not preferred,
2945 // because the manual said that integer mod and div operations are slow on gpu
2946 pos_d = (reduce_ii % dims_a[rd]);
2947 reduce_ii = (reduce_ii / dims_a[rd]);
2948 }
2949 else
2950 {
2951 pos_d = (reduce_ii & ((1 << log2_dims_a[rd])-1)); //take the lower log2_dims bits
2952 reduce_ii = (reduce_ii >> log2_dims_a[rd]); //shift those lower log2_dims bits off of ii
2953 }
2954 a_data_ri += pos_d * a_str[rd];
2955 }
2956 }
2957 sum += a_data_ri[0];
2958 }
2959 }
2960 }
2961 z_data_i[0] = sum;
2962 }
2963 }
2964
2965 static __global__ void kernel_reduce_sum_1011(
2966 const unsigned int d0,
2967 const unsigned int d1,
2968 const unsigned int d2,
2969 const unsigned int d3,
2970 const float *A, const int sA0, const int sA1, const int sA2, const int sA3,
2971 float * Z, const int sZ0)
2972 {
2973 const int threadCount = blockDim.x * blockDim.y * blockDim.z;
2974 const int threadNum = threadIdx.z * blockDim.x * blockDim.y + threadIdx.y * blockDim.x + threadIdx.x;
2975 extern __shared__ float buf[];
2976 float mysum = 0.0f;
2977
2978 if (warpSize != 32)
2979 {
2980 return; //TODO: set error code
2981 }
2982
2983 for (int i0 = threadIdx.z; i0 < d0; i0 += blockDim.z)
2984 {
2985 float Ai = A[i0 * sA0 + blockIdx.x * sA1 + threadIdx.y * sA2 + threadIdx.x * sA3];
2986 mysum += Ai;
2987 }
2988 buf[threadNum] = mysum;
2989 __syncthreads();
2990
2991 // rest of function is handled by one warp
2992 if (threadNum < warpSize)
2993 {
2994 for (int i = threadNum + warpSize; i < threadCount; i += warpSize)
2995 {
2996 mysum += buf[i];
2997 }
2998 buf[threadNum] = mysum;
2999 if (threadNum < 16)
3000 {
3001 //reduce so that threadNum 0 has the sum of everything
3002 if(threadNum + 16 < threadCount) buf[threadNum] += buf[threadNum+16];
3003 if(threadNum + 8 < threadCount) buf[threadNum] += buf[threadNum+8];
3004 if(threadNum + 4 < threadCount) buf[threadNum] += buf[threadNum+4];
3005 if(threadNum + 2 < threadCount) buf[threadNum] += buf[threadNum+2];
3006 if(threadNum + 1 < threadCount) buf[threadNum] += buf[threadNum+1];
3007 if (threadNum == 0)
3008 {
3009 Z[blockIdx.x*sZ0] = buf[0];
3010 }
3011 }
3012 }
3013 }
3014 /**
3015 * Dimensions in which the self has size 1 and A has size > 1 are considered summing dimensions
3016 * Dimensions in which self has size > 1 and A has size > 1 are considered non-summing dimensions, and in this case their sizes must be equal.
3017 */
3018 int
3019 CudaNdarray_reduce_sum(CudaNdarray * self, CudaNdarray * A)
3020 {
3021 int verbose = 0;
3022 //check input rank
3023 if (self->nd != A->nd)
3024 {
3025 PyErr_Format(PyExc_TypeError, "Rank mismatch in CudaNdarray_sum: %i vs %i", self->nd, A->nd);
3026 return -1;
3027 }
3028 for (int i = 0; i < self->nd; ++i)
3029 {
3030 if ((CudaNdarray_HOST_DIMS(self)[i] > 1) && (CudaNdarray_HOST_DIMS(self)[i] != CudaNdarray_HOST_DIMS(A)[i]))
3031 {
3032 PyErr_Format(PyExc_TypeError, "Dimension mismatch in CudaNdarray_sum: self->dim[%i] == %i , A->dim[%i] = %i",
3033 i, CudaNdarray_HOST_DIMS(self)[i], i, CudaNdarray_HOST_DIMS(A)[i]);
3034 return -1;
3035 }
3036 }
3037
3038 int n_summations = (unsigned int)CudaNdarray_SIZE(self);
3039 if (verbose)
3040 {
3041 std::cerr << "reduce_sum n_summations " << n_summations << '\n';
3042 std::cerr << "reduce_sum nd " << self->nd << '\n';
3043 fprint_CudaNdarray(stderr, A);
3044 fprint_CudaNdarray(stderr, self);
3045 }
3046 if (0 && (A->nd == 4) //check to see if kernel_reduce_sum_1011 applies
3047 && (CudaNdarray_HOST_DIMS(self)[0] == 1)
3048 && (CudaNdarray_HOST_DIMS(self)[2] == 1)
3049 && (CudaNdarray_HOST_DIMS(self)[3] == 1)
3050 )
3051 {
3052 dim3 n_threads(CudaNdarray_HOST_DIMS(A)[3], CudaNdarray_HOST_DIMS(A)[2]);
3053 dim3 n_blocks(CudaNdarray_HOST_DIMS(A)[1]);
3054 while (n_threads.x * n_threads.y * n_threads.z < NUM_VECTOR_OP_THREADS_PER_BLOCK) ++n_threads.z;
3055 n_threads.z -= 1;
3056 if (n_threads.z > 64) n_threads.z = 64;
3057 if (n_threads.z)
3058 {
3059 if (verbose) printf("trying kernel_reduce_sum_1011\n");
3060 int n_shared = sizeof(float) * n_threads.x * n_threads.y * n_threads.z;
3061 kernel_reduce_sum_1011<<<n_blocks, n_threads, n_shared>>>(
3062 CudaNdarray_HOST_DIMS(A)[0],
3063 CudaNdarray_HOST_DIMS(A)[1],
3064 CudaNdarray_HOST_DIMS(A)[2],
3065 CudaNdarray_HOST_DIMS(A)[3],
3066 CudaNdarray_DEV_DATA(A),
3067 CudaNdarray_HOST_STRIDES(A)[0],
3068 CudaNdarray_HOST_STRIDES(A)[1],
3069 CudaNdarray_HOST_STRIDES(A)[2],
3070 CudaNdarray_HOST_STRIDES(A)[3],
3071 CudaNdarray_DEV_DATA(self),
3072 CudaNdarray_HOST_STRIDES(self)[1]);
3073 CNDA_THREAD_SYNC;
3074 if (cudaSuccess == cudaGetLastError()) return 0;
3075 if (verbose) printf("failed, falling back to kernel_reduce_sum\n");
3076 }
3077 }
3078
3079 int n_threads_per_block = std::min(n_summations,
3080 NUM_VECTOR_OP_THREADS_PER_BLOCK);
3081 int n_blocks = std::min(ceil_intdiv(n_summations,n_threads_per_block),
3082 NUM_VECTOR_OP_BLOCKS);
3083 int n_structure_cache = self->nd * 4 * sizeof(int);
3084
3085 if (verbose)
3086 {
3087 std::cerr << "n_blocks, n_threads_per_block " << n_blocks << ' ' << n_threads_per_block << '\n';
3088 }
3089 assert (self->nd > 0);
3090 assert (self->nd == A->nd);
3091 kernel_reduce_sum<<<n_blocks, n_threads_per_block, n_structure_cache>>>(
3092 n_summations,
3093 self->nd,
3094 CudaNdarray_DEV_DIMS(A),
3095 CudaNdarray_DEV_LOG2DIMS(A),
3096 CudaNdarray_DEV_STRIDES(A),
3097 CudaNdarray_DEV_DATA(A),
3098 CudaNdarray_DEV_STRIDES(self),
3099 CudaNdarray_DEV_DATA(self));
3100 CNDA_THREAD_SYNC;
3101 cudaError_t err = cudaGetLastError();
3102 if (cudaSuccess != err)
3103 {
3104 PyErr_Format(PyExc_RuntimeError, "Cuda error: %s: %s.\n", "kernel_reduce_sum", cudaGetErrorString(err));
3105 return -1;
3106 }
3107 return 0;
3108 }
3109 int
3110 CudaNdarray_reduce_prod(CudaNdarray * self, const CudaNdarray * A)
3111 {
3112 PyErr_SetString(PyExc_NotImplementedError, "");
3113 return -1;
3114 }
3115 int
3116 CudaNdarray_reduce_min(CudaNdarray * self, const CudaNdarray * A)
3117 {
3118 PyErr_SetString(PyExc_NotImplementedError, "");
3119 return -1;
3120 }
3121 int
3122 CudaNdarray_reduce_max(CudaNdarray * self, const CudaNdarray * A)
3123 {
3124 PyErr_SetString(PyExc_NotImplementedError, "");
3125 return -1;
3126 }
3127
3128
3129 /**
3130 *
3131 * pattern is a permutation of [0, 1, ... self->nd-1] with the following twists:
3132 * - an element 'd' of the permutation can be dropped if CudaNdarray_HOST_DIMS(self)[d] == 1
3133 * - any number of '-1' elements can be in the pattern, and they will cause new ranks (with dim==1) to be inserted.
3134 *
3135 * For example, if CudaNdarray_HOST_DIMS(self) == [4, 5, 1, 6], and pattern = [0,3,-1,-1, 1], then CudaNdarray_HOST_DIMS(self) would be modified to become:
3136 * [4, 6, 1, 1, 5] (we dropped the original dim[2]==1, and inserted two singleton dimensions with the -1s.
3137 */
3138 int
3139 CudaNdarray_dimshuffle(CudaNdarray * self, unsigned int len, const int * pattern)
3140 {
3141 //TODO: pass a workspace pointer to avoid the internal malloc
3142 int * newdims = (int *)malloc(sizeof(int) * (len + len + self->nd)); //we tack on the taken buffer here for speed of not having to malloc twice.
3143 int * newstrides = newdims + len;
3144 int * dims_taken = newstrides + len;
3145 if (!newdims)
3146 {
3147 PyErr_SetString(PyExc_MemoryError, "CudaNdarray_dimshuffle: Failed to allocate temporary space");
3148 return -1;
3149 }
3150 for (int i = 0; i < self->nd; ++i)
3151 {
3152 dims_taken[i] = 0;
3153 }
3154 for (int i = 0; i < len; ++i)
3155 {
3156 if (pattern[i] < 0)
3157 {
3158 newdims[i] = 1;
3159 newstrides[i] = 0;
3160 }
3161 else if(dims_taken[pattern[i]])
3162 {
3163 PyErr_Format(PyExc_ValueError, "Cudandarray_dimshuffle: invalid pattern for Cudandarray_dimshuffle. You used the dimensions %d multiple time",
3164 pattern[i]);
3165 free(newdims);
3166 return -1;
3167 }
3168 else if (pattern[i]>= self->nd)
3169 {
3170 PyErr_Format(PyExc_ValueError, "Cudandarray_dimshuffle: invalid pattern for Cudandarray_dimshuffle. You asked for a dimensions that don't exist %d for a %d dims CudaNdarray",
3171 pattern[i], self->nd);
3172 free(newdims);
3173 return -1;
3174 }
3175 else
3176 {
3177 newdims[i] = CudaNdarray_HOST_DIMS(self)[pattern[i]];
3178 newstrides[i] = CudaNdarray_HOST_STRIDES(self)[pattern[i]];
3179 dims_taken[pattern[i]] = 1;
3180 }
3181 }
3182 //Check if we dropped not broadcastable dims
3183 for (int i = 0; i < self->nd; ++i)
3184 {
3185 if (dims_taken[i]==0 && CudaNdarray_HOST_DIMS(self)[i]!=1)
3186 {
3187 PyErr_SetString(PyExc_ValueError, "Cudandarray_dimshuffle: You cannot drop a non-broadcastable dimension.");
3188 free(newdims);
3189 return -1;
3190 }
3191 }
3192 //swap this structure in for the one in self, and sync to the card
3193 if (CudaNdarray_set_nd(self, len))
3194 {
3195 free(newdims);
3196 return -1;
3197 }
3198 for (int i = 0; i < len; ++i)
3199 {
3200 CudaNdarray_set_dim(self, i, newdims[i]);
3201 CudaNdarray_set_stride(self, i, newstrides[i]);
3202 }
3203 if (cnda_copy_structure_to_device(self))
3204 {
3205 free(newdims);
3206 return -1;
3207 }
3208 free(newdims);
3209 return 0;
3210 }
3211
3212
3213
3214 /**
3215 *
3216 * This is the function that bind to python.
3217 * See CudaNdarray_dimshuffle to call from C.
3218 * We use -1 to mean 'x' as in Tensor Dimshuffle.
3219 */
3220 PyObject *
3221 CudaNdarray_Dimshuffle(PyObject* _unused, PyObject* args)
3222 {
3223 PyObject * self = NULL;
3224 PyObject * pattern_object = NULL;
3225 int * pattern = NULL;
3226 PyObject * rval = NULL;
3227 int success = -1;
3228 //const int * dims = NULL;
3229
3230 //args should consist of two python objects ("OO")
3231 if (! PyArg_ParseTuple(args, "OO", &self, &pattern_object))
3232 return NULL;
3233
3234 if (!CudaNdarray_Check(self) )
3235 {
3236 PyErr_SetString(PyExc_TypeError, "First argument to cuda_ndarray.dimshuffle must be a CudaNdarray");
3237 return NULL;
3238 }
3239
3240 //parse pattern_object into int * pattern
3241
3242 Py_ssize_t pattern_dim = PyObject_Length(pattern_object);
3243
3244 if (pattern_dim < 0)
3245 {
3246 PyErr_SetString(PyExc_TypeError, "Couldn't get length of third argument to cuda_ndarray.dimshuffle");
3247 return NULL;
3248 }
3249
3250 pattern = (int *) malloc( pattern_dim * sizeof(int));
3251
3252 for (Py_ssize_t i = 0; i < pattern_dim; i++)
3253 {
3254 PyObject * idx = PyLong_FromLong(i);
3255
3256 if (idx == NULL)
3257 {
3258 PyErr_SetString(PyExc_Exception, "Couldn't make long object to loop over list/tuple");
3259 goto CudaNdarray_dimshuffle_fail;
3260 }
3261
3262 long elem_value = 0;
3263
3264 PyObject * elem = PyObject_GetItem(pattern_object, idx);
3265
3266 if (elem == NULL)
3267 {
3268 Py_XDECREF( elem);
3269 PyErr_SetString(PyExc_ValueError, "Third argument to dimshuffle must be list or tuple of integers");
3270 goto CudaNdarray_dimshuffle_fail;
3271 }
3272
3273 elem_value = PyInt_AsLong(elem);
3274
3275 if (elem_value == -1 && PyErr_Occurred() )
3276 {
3277 Py_XDECREF(elem);
3278 PyErr_SetString(PyExc_ValueError, "Third argument to dimshuffle must be list or tuple of integers");
3279 goto CudaNdarray_dimshuffle_fail;
3280 }
3281
3282 pattern[i] = elem_value;
3283
3284 Py_XDECREF( elem );
3285 Py_XDECREF( idx );
3286 }
3287
3288 //allocate rval
3289 rval = (PyObject *) CudaNdarray_View((CudaNdarray *) self);
3290
3291 if (rval == NULL)
3292 {
3293 //CudaNdarray_New should have set the exception string
3294 goto CudaNdarray_dimshuffle_fail;
3295 }
3296
3297
3298 //printf("pattern_dim: %d\n",pattern_dim);
3299 //printf("pattern: %d %d\n",pattern[0],pattern[1]);
3300 //dims = CudaNdarray_HOST_DIMS( (CudaNdarray *) self);
3301 //printf("dims before: %d %d\n",dims[0],dims[1]);
3302
3303 success = CudaNdarray_dimshuffle((CudaNdarray *) rval, pattern_dim, pattern);
3304
3305 if (success != 0)
3306 {
3307 //Exception string should already be set by CudaNdarray_dimshuffle
3308 goto CudaNdarray_dimshuffle_fail;
3309 }
3310
3311 free(pattern);
3312
3313 return rval;
3314
3315 CudaNdarray_dimshuffle_fail:
3316
3317 if (pattern != NULL)
3318 free(pattern);
3319
3320 Py_XDECREF(rval);
3321 return NULL;
3322 }
3323
3324 /*
3325 Local Variables:
3326 mode:c++
3327 c-basic-offset:4
3328 c-file-style:"stroustrup"
3329 c-file-offsets:((innamespace . 0)(inline-open . 0))
3330 indent-tabs-mode:nil
3331 fill-column:79
3332 End:
3333 */
3334 // vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=79 :
3335
ERROR (theano.sandbox.cuda): Failed to compile cuda_ndarray.cu: ('nvcc return status', 1, 'for cmd', 'nvcc -shared -g -m64 -Xcompiler -fPIC,-m64 -Xlinker -rpath,/usr/local/cuda/lib -I/Users/bayerj/devel/third-party/Theano/theano/sandbox/cuda -I/Library/Python/2.6/site-packages/numpy-2.0.0.dev_ac2c160_20110325-py2.6-macosx-10.6-universal.egg/numpy/core/include -I/System/Library/Frameworks/Python.framework/Versions/2.6/include/python2.6 -o /Users/bayerj/.theano/compiledir_Darwin-10.6.0-i386-64bit-i386-2.6.1/cuda_ndarray/cuda_ndarray.so mod.cu -L/usr/local/cuda/lib -lcublas -lcudart python version 2.6.1 can\'t run /usr/bin/python-config. Try the alternative(s): /usr/bin/python2.5-config (uses python 2.5) /usr/bin/python2.6-config (uses python 2.6) Run "man python" for more information about multiple version support in Mac OS X.')
WARNING (theano.sandbox.cuda): Cuda is disabled, cuda-based code will thus not be working properly
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment