· 9 years ago · Dec 05, 2016, 07:13 PM
1iMacDavid:~ david$ python3 test_theano.py
21 #define _CUDA_NDARRAY_C
32
43 #include <Python.h>
54 #include <structmember.h>
65 #include "theano_mod_helper.h"
76
87 #include <numpy/arrayobject.h>
98 #include <iostream>
109
1110 #include "cuda_ndarray.cuh"
1211
1312 #ifndef CNMEM_DLLEXPORT
1413 #define CNMEM_DLLEXPORT
1514 #endif
1615
1716 #include "cnmem.h"
1817 #include "cnmem.cpp"
1918
2019 //If true, when there is a gpu malloc or free error, we print the size of allocated memory on the device.
2120 #define COMPUTE_GPU_MEM_USED 0
2221
2322 //If true, we fill with NAN allocated device memory.
2423 #define ALLOC_MEMSET 0
2524
2625 //If true, we print out when we free a device pointer, uninitialize a
2726 //CudaNdarray, or allocate a device pointer
2827 #define PRINT_FREE_MALLOC 0
2928
3029 //If true, we do error checking at the start of functions, to make sure there
3130 //is not a pre-existing error when the function is called.
3231 //You probably need to set the environment variable
3332 //CUDA_LAUNCH_BLOCKING=1, and/or modify the CNDA_THREAD_SYNC
3433 //preprocessor macro in cuda_ndarray.cuh
3534 //if you want this to work.
3635 #define PRECHECK_ERROR 0
3736
3837 cublasHandle_t handle = NULL;
3938 int* err_var = NULL;
4039
4140 /////////////////////////
4241 // Alloc and Free
4342 /////////////////////////
4443
4544 static int g_gpu_context_active = 0;
4645
4746
4847 PyObject *
4948 CudaNdarray_Dimshuffle(PyObject* _unused, PyObject* args);
5049 static PyObject *CudaNdarray_get_shape(CudaNdarray *self, void *closure);
5150
5251
5352 /**
5453 *
5554 * In the test program I'm using, the _outstanding_mallocs decreases with every call.
5655 * This suggests there are more free() calls being made than alloc(), but I can't figure out why.
5756 *
5857 */
5958 int _outstanding_mallocs[] = {0,0};
6059
6160 #if COMPUTE_GPU_MEM_USED
6261 size_t _allocated_size = 0;
6362 size_t _max_allocated_size = 0;
6463
6564 const int TABLE_SIZE = 10000;
6665 struct table_struct{
6766 void* ptr;
6867 size_t size;
6968 };
7069 table_struct _alloc_size_table[TABLE_SIZE];
7170 #endif
7271
7372 void * device_malloc(size_t size)
7473 {
7574 return device_malloc(size, VERBOSE_DEVICE_MALLOC);
7675 }
7776
7877 static bool g_use_cnmem = false;
7978 static const int g_max_devices = 8;
8079 int initCnmem(int card_number_provided, int card_nb, size_t mem) {
8180 static bool cnmemInitialized = false;
8281 if(cnmemInitialized) {
8382 return 0;
8483 }
8584 // On stderr to be at the same place as "Using gpu device..."
8685 int numDevices = 0;
8786 cnmemDevice_t devices[g_max_devices];
8887 if(cudaGetDeviceCount(&numDevices) != cudaSuccess) {
8988 PyErr_Format(PyExc_RuntimeError,
9089 "initCnmem: 'cudaGetDeviceCount' failed! Reason=%s\n",
9190 cudaGetErrorString(cudaGetLastError()));
9291 return -1;
9392 }
9493 if(card_number_provided){
9594 numDevices = 1;
9695 int i = 0;
9796 devices[i].device = card_nb;
9897 devices[i].size = mem;
9998 ///@TODO: thejaswi: add support for multiple streams
10099 devices[i].numStreams = 0;
101100 devices[i].streams = NULL;
102101 devices[i].streamSizes = NULL;
103102 }else{
104103 for(int i=0;i<numDevices;++i) {
105104 devices[i].device = i;
106105 devices[i].size = mem;
107106 ///@TODO: thejaswi: add support for multiple streams
108107 devices[i].numStreams = 0;
109108 devices[i].streams = NULL;
110109 }
111110 }
112111
113112 ///@TODO: thejaswi: passing custom cnmem flags?
114113 cnmemStatus_t status = cnmemInit(numDevices, devices, CNMEM_FLAGS_DEFAULT);
115114 if(status != CNMEM_STATUS_SUCCESS) {
116115 PyErr_Format(PyExc_RuntimeError,
117116 "initCnmem: cnmemInit call failed! Reason=%s. numdev=%d\n",
118117 cnmemGetErrorString(status), numDevices);
119118 return -1;
120119 }
121120 cnmemInitialized = true;
122121 return 0;
123122 }
124123
125124 void * device_malloc(size_t size, int verbose)
126125 {
127126 #if PRECHECK_ERROR
128127 cudaThreadSynchronize();
129128 cudaError_t prevError = cudaGetLastError();
130129 if (cudaSuccess != prevError)
131130 {
132131 fprintf(stderr,
133132 "Error existed before calling device_malloc. %s\n",
134133 cudaGetErrorString(prevError)
135134 );
136135 }
137136 #endif
138137 void * rval=NULL;
139138 ///@TODO: thejaswi: support for multiple-streams?
140139 if(g_use_cnmem) {
141140 cnmemStatus_t status = CNMEM_STATUS_SUCCESS;
142141 status = cnmemMalloc(&rval, size, NULL);
143142 if(status != CNMEM_STATUS_SUCCESS) {
144143 PyErr_Format(PyExc_MemoryError,
145144 "Error allocating %llu bytes of device memory (%s).",
146145 (unsigned long long)size, cnmemGetErrorString(status));
147146 return NULL;
148147 }
149148 }
150149 else {
151150 cudaError_t err = cudaMalloc(&rval, size);
152151 if (cudaSuccess != err)
153152 {
154153 // Clear the error flag, cudaMalloc doesn't do it.
155154 // Currently this returns the same thing as err, but if in future
156155 // it returns something else I still don't see why we should ignore
157156 // it. All we want to do here is reset the flag.
158157 cudaGetLastError();
159158 if (verbose)
160159 {
161160 size_t free = 0, total = 0;
162161 cudaError_t err2 = cudaMemGetInfo(&free, &total);
163162 if (err2 != cudaSuccess){
164163 cudaGetLastError();
165164 fprintf(stderr,
166165 "Error when trying to find the memory information"
167166 " on the GPU: %s\n", cudaGetErrorString(err2));
168167 }
169168 #if COMPUTE_GPU_MEM_USED
170169 fprintf(stderr,
171170 "Error allocating %llu bytes of device memory (%s)."
172171 " new total bytes allocated: %llu."
173172 " Driver report %llu bytes free and %llu bytes total \n",
174173 (unsigned long long)size, cudaGetErrorString(err), (unsigned long long)_allocated_size,
175174 (unsigned long long)free, (unsigned long long)total);
176175 #else
177176 fprintf(stderr,
178177 "Error allocating %llu bytes of device memory (%s)."
179178 " Driver report %llu bytes free and %llu bytes total \n",
180179 (unsigned long long)size, cudaGetErrorString(err), (unsigned long long)free, (unsigned long long)total);
181180 #endif
182181 }
183182 PyErr_Format(PyExc_MemoryError,
184183 "Error allocating %llu bytes of device memory (%s).",
185184 (unsigned long long)size, cudaGetErrorString(err));
186185 return NULL;
187186 }
188187 }
189188 if (rval != NULL){
190189 // Can it happen that cudaMalloc return cudaSuccess, but return a NULL ptr?
191190 // Could this be what happen if size is 0?
192191 _outstanding_mallocs[0] += 1;
193192
194193 #if COMPUTE_GPU_MEM_USED
195194 _allocated_size += size;
196195 _max_allocated_size = std::max(_max_allocated_size, _allocated_size);
197196 int i = 0;
198197 for(;i<TABLE_SIZE;i++){
199198 if(NULL==_alloc_size_table[i].ptr){
200199 _alloc_size_table[i].ptr=rval;
201200 _alloc_size_table[i].size=size;
202201 break;
203202 }
204203 }
205204 if (i == TABLE_SIZE){
206205 fprintf(stderr,
207206 "When tracking GPU malloc, our table size wasn't big enough."
208207 " So we loose some tracking. Raise the value of TABLE_SIZE in the file cuda_ndarra.cu");
209208 }
210209 #endif
211210 }
212211 //fprintf(stderr,
213212 //"allocated %li bytes of device memory (%s). new total bytes allocated: %d. ptr: %p\n",
214213 //(long)size, cudaGetErrorString(err),_allocated_size,rval);
215214
216215 if(ALLOC_MEMSET){
217216 //We init them to nan to make sure we catch more debug case.
218217 cudaMemset(rval, 0xFF, size);
219218 //printf("MEMSET\n");
220219 }
221220 #if PRINT_FREE_MALLOC
222221 fprintf(stderr, "device malloc %p of size %d\n", rval, size);
223222 #endif
224223 return rval;
225224 }
226225
227226 int device_free(void *ptr)
228227 {
229228 #if PRECHECK_ERROR
230229 cudaThreadSynchronize();
231230 cudaError_t prevError = cudaGetLastError();
232231 if (cudaSuccess != prevError)
233232 {
234233 fprintf(stderr,
235234 "Error existed before calling device_free. %s\n",
236235 cudaGetErrorString(prevError)
237236 );
238237 }
239238 #endif
240239 #if PRINT_FREE_MALLOC
241240 size_t free = 0, total = 0;
242241 cudaError_t err2 = cudaMemGetInfo(&free, &total);
243242 if (err2 != cudaSuccess){
244243 cudaGetLastError();
245244 fprintf(stderr,
246245 "Error when tring to find the memory information"
247246 " on the GPU: %s\n", cudaGetErrorString(err2));
248247 }
249248 #if COMPUTE_GPU_MEM_USED
250249 {
251250 int i = 0;
252251 for(;i<TABLE_SIZE;i++)
253252 if(_alloc_size_table[i].ptr==ptr){
254253 break;
255254 }
256255 assert(i<TABLE_SIZE);
257256 fprintf(stderr, "device_free %p of size %d."
258257 " Driver report %d bytes free and %d bytes total \n",
259258 ptr, _alloc_size_table[i].size, free, total);
260259 }
261260 #else
262261 fprintf(stderr, "device_free %p."
263262 " Driver report %d bytes free and %d bytes total \n",
264263 ptr, free, total);
265264 #endif
266265 #endif
267266
268267 // if there is no gpu context, the call to cudaFree will fail; skip it entirely
269268 if(!g_gpu_context_active) {
270269 return 0;
271270 }
272271
273272 ///@TODO: thejaswi: multi-stream support
274273 if(g_use_cnmem) {
275274 cnmemStatus_t status = cnmemFree(ptr, NULL);
276275 if(status != CNMEM_STATUS_SUCCESS) {
277276 fprintf(stderr, "device_free: cnmemFree call failed! Reason=%s\n",
278277 cnmemGetErrorString(status));
279278 }
280279 }
281280 else {
282281 // We need sync as the Theano's GC could remove intermediate variable that
283282 // are still needed as the gpu kernel are running or in the queue.
284283 CNDA_BEGIN_ALLOW_THREADS
285284 cudaThreadSynchronize();
286285 CNDA_END_ALLOW_THREADS
287286
288287 cudaError_t err = cudaFree(ptr);
289288 if (cudaSuccess != err)
290289 {
291290 // Clear the error flag, cudaFree doesn't do it.
292291 // Currently this returns the same thing as err, but if in future
293292 // it returns something else I still don't see why we should ignore
294293 // it. All we want to do here is reset the flag.
295294 cudaGetLastError();
296295 size_t free = 0, total = 0;
297296 cudaError_t err2 = cudaMemGetInfo(&free, &total);
298297 if (err2 != cudaSuccess){
299298 cudaGetLastError();
300299 fprintf(stderr,
301300 "Error when tring to find the memory information"
302301 " on the GPU: %s\n", cudaGetErrorString(err2));
303302 }
304303 #if COMPUTE_GPU_MEM_USED
305304 {
306305 int i = 0;
307306 for(;i<TABLE_SIZE;i++)
308307 if(_alloc_size_table[i].ptr==ptr){
309308 break;
310309 }
311310 assert(i<TABLE_SIZE);
312311 fprintf(stderr,
313312 "Error freeing device pointer %p (%s) of size %llu. %llu byte already allocated."
314313 " Driver report %llu bytes free and %llu bytes total \n",
315314 ptr, cudaGetErrorString(err),
316315 (unsigned long long)_alloc_size_table[i].size, (unsigned long long)_allocated_size, (unsigned long long)free, (unsigned long long)total);
317316 }
318317 #else
319318 fprintf(stderr,
320319 "Error freeing device pointer %p (%s)."
321320 " Driver report %llu bytes free and %llu bytes total \n",
322321 ptr,
323322 cudaGetErrorString(err), (unsigned long long)free, (unsigned long long)total);
324323 #endif
325324 if (NULL != PyErr_Occurred()){
326325 fprintf(stderr,
327326 "device_free: cudaFree() returned an error, but there is already an"
328327 " Python error set. This happen during the clean up when there is a"
329328 " first error and the CUDA driver is in a so bad state that it don't"
330329 " work anymore. We keep the previous error set to help debugging it.");
331330 return -1;
332331 }
333332 PyErr_Format(PyExc_MemoryError,
334333 "error freeing device pointer %p (%s)",
335334 ptr,
336335 cudaGetErrorString(err));
337336 return -1;
338337 }
339338 }
340339 _outstanding_mallocs[0] -= (ptr != NULL);
341340 #if COMPUTE_GPU_MEM_USED
342341 int i=0;
343342 size_t total_freed = 0;
344343 for(;i<TABLE_SIZE;i++)
345344 if(_alloc_size_table[i].ptr==ptr){
346345 _allocated_size -= _alloc_size_table[i].size;
347346 total_freed += _alloc_size_table[i].size;
348347 _alloc_size_table[i].ptr=0;
349348 _alloc_size_table[i].size=0;
350349
351350 break;
352351 }
353352 //if(i==TABLE_SIZE)
354353 // printf("Unallocated unknow size!\n");
355354 //fprintf(stderr, "freed %li bytes of device memory (%s). %d already allocated, ptr=%p\n", (long)total_freed, cudaGetErrorString(err),_allocated_size,ptr);
356355 #endif
357356 return 0;
358357 }
359358
360359 static PyObject *
361360 outstanding_mallocs(PyObject* self, PyObject * args)
362361 {
363362 return PyInt_FromLong(_outstanding_mallocs[0]);
364363 }
365364
366365
367366 static void *work_mem = NULL;
368367 static size_t work_size = 0;
369368
370369 /*
371370 * Returns a chunk of memory for temporary work inside of an op. You can only
372371 * request a single chunk of memory at a time since it is reused.
373372 */
374373 void *get_work_mem(size_t sz) {
375374 if (sz <= work_size)
376375 return work_mem;
377376 device_free(work_mem);
378377 work_mem = device_malloc(sz);
379378 work_size = sz;
380379 if (work_mem == NULL)
381380 work_size = 0;
382381 return work_mem;
383382 }
384383
385384 /////////////////////////
386385 // Static helper methods
387386 /////////////////////////
388387
389388 static void
390389 CudaNdarray_null_init(CudaNdarray*self)
391390 {
392391 self->base = NULL;
393392 self->nd = -1;
394393 self->host_structure = NULL;
395394 self->data_allocated = 0;
396395 self->dev_structure_fresh = 1;
397396 self->dev_structure = NULL;
398397 self->devdata = NULL;
399398 }
400399
401400 static int
402401 CudaNdarray_uninit(CudaNdarray*self)
403402 {
404403 #if PRINT_FREE_MALLOC
405404 fprintf(stderr, "CudaNdarray_uninit %p\n", self);
406405 #endif
407406 int rval = 0;
408407 if (self->data_allocated) {
409408 assert(self->devdata);
410409 if (device_free(self->devdata))
411410 {
412411 fprintf(stderr,
413412 "CudaNdarray_uninit: error freeing self->devdata. (self=%p, self->devata=%p)\n",
414413 self, self->devdata);
415414 rval = -1;
416415 }
417416 self->devdata = NULL;
418417 self->data_allocated = 0;
419418 }
420419 if (self->dev_structure)
421420 {
422421 if (device_free(self->dev_structure))
423422 {
424423 fprintf(stderr,
425424 "CudaNdarray_uninit: error freeing dev_structure memory %p (self=%p)\n",
426425 self->dev_structure, self);
427426 rval = -1;
428427 }
429428 self->dev_structure = NULL;
430429 }
431430 if (self->host_structure)
432431 {
433432 free(self->host_structure);
434433 self->host_structure = NULL;
435434 }
436435 self->nd = -1;
437436 Py_XDECREF(self->base);
438437 self->base = NULL;
439438 return rval;
440439 }
441440
442441
443442 //make the rightmost coords change fastest
444443 //TODO: why does a downward for-loop not work????
445444 //TODO: use the log2_dims and driver code to remove / and %
446445 //TODO: skip the last division (when d == 0)
447446 #define decl_k_elemwise_unary_rowmajor(name, F) \
448447 __global__ void name (unsigned int numEls, \
449448 unsigned int nd, \
450449 const int * dim, \
451450 const float * a_data, const int * a_str, \
452451 float * z_data, const int * z_str) \
453452 { \
454453 const unsigned int idx = blockIdx.x * blockDim.x + threadIdx.x; \
455454 const unsigned int numThreads = blockDim.x * gridDim.x; \
456455 \
457456 for (unsigned int i = idx; i < numEls; i += numThreads) \
458457 { \
459458 unsigned int ii = i; \
460459 const float * a_i = a_data; \
461460 float * z_i = z_data; \
462461 for (unsigned int _d = 0; _d < nd; ++_d) \
463462 { \
464463 unsigned int d = nd - _d-1; \
465464 int i_d = ii % dim[d]; /* i_d is our position in the d'th dimension */ \
466465 ii = ii / dim[d]; \
467466 a_i += i_d * a_str[d]; /* increment our a and z pointers by i_d elements */ \
468467 z_i += i_d * z_str[d]; \
469468 } \
470469 z_i[0] = F(a_i[0]); \
471470 } \
472471 }
473472
474473 template<typename T> __device__ T unary_copy(T a) { return a; }
475474 decl_k_elemwise_unary_rowmajor(k_elemwise_unary_rowmajor_copy, unary_copy<float>)
476475
477476 template<typename T> __device__ T unary_exp(T a) { return exp(a); }
478477 decl_k_elemwise_unary_rowmajor(k_elemwise_unary_rowmajor_exp, unary_exp<float>)
479478
480479 /////////////////////////////
481480 // Satisfying reqs to be Type
482481 /////////////////////////////
483482
484483 //DON'T use directly(if their is other CudaNdarray that point to it, it will cause problem)! use Py_DECREF() instead
485484 static void
486485 CudaNdarray_dealloc(CudaNdarray* self)
487486 {
488487 if (0) std::cerr << "CudaNdarray dealloc " << self << " " << self->devdata << '\n';
489488 if(Py_REFCNT(self) > 1)
490489 printf("WARNING:CudaNdarray_dealloc called when there is still active reference to it.\n");
491490 CudaNdarray_uninit(self);
492491 Py_TYPE(self)->tp_free((PyObject*)self);
493492 --_outstanding_mallocs[1];
494493 if (0)
495494 {
496495 fprintf(stderr, "device_malloc_counts: (device) %i (obj) %i\n",
497496 _outstanding_mallocs[0],
498497 _outstanding_mallocs[1]);
499498 }
500499 }
501500
502501 static PyObject *
503502 CudaNdarray_new(PyTypeObject *type, PyObject *args, PyObject *kwds)
504503 {
505504 CudaNdarray *self;
506505
507506 self = (CudaNdarray *)type->tp_alloc(type, 0);
508507 if (self != NULL)
509508 {
510509 CudaNdarray_null_init(self);
511510 ++_outstanding_mallocs[1];
512511 }
513512 return (PyObject *)self;
514513 }
515514 static int
516515 CudaNdarray_init(CudaNdarray *self, PyObject *args, PyObject *kwds)
517516 {
518517 PyObject *arr=NULL;
519518
520519 if (! PyArg_ParseTuple(args, "O", &arr))
521520 return -1;
522521 if (! PyArray_Check(arr))
523522 {
524523 PyErr_SetString(PyExc_TypeError, "PyArray arg required");
525524 return -1;
526525 }
527526 int rval = CudaNdarray_CopyFromArray(self, (PyArrayObject*)arr);
528527 return rval;
529528 }
530529 static PyMemberDef CudaNdarray_members[] =
531530 {
532531 /*
533532 {"first", T_OBJECT_EX, offsetof(CudaNdarray, first), 0,
534533 "first name"},
535534 {"last", T_OBJECT_EX, offsetof(CudaNdarray, last), 0,
536535 "last name"},
537536 {"number", T_INT, offsetof(CudaNdarray, number), 0,
538537 "noddy number"},
539538 */
540539 {NULL} /* Sentinel */
541540 };
542541
543542 PyObject * CudaNdarray_CreateArrayObj(CudaNdarray * self, PyObject *args)
544543 {
545544 PyObject * dtype = NULL;
546545 if (args && !PyArg_ParseTuple(args, "|O", &dtype))
547546 return NULL;
548547 if (dtype) {
549548 PyArray_Descr* dtype2;
550549 // PyArray_DescrConverter try to convert anything to a PyArray_Descr.
551550 if(!PyArray_DescrConverter(dtype, &dtype2))
552551 {
553552 PyObject * str = PyObject_Repr(dtype);
554553 PyErr_Format(PyExc_TypeError,
555554 "CudaNdarray dtype parameter not understood: %s",
556555 PyString_AsString(str)
557556 );
558557 Py_CLEAR(str);
559558 return NULL;
560559 }
561560 int typeNum = dtype2->type_num;
562561 Py_DECREF(dtype2);
563562 if (typeNum != NPY_FLOAT32)
564563 {
565564 PyObject * str = PyObject_Repr(dtype);
566565 PyErr_Format(PyExc_TypeError,
567566 "CudaNdarray support only support float32 dtype, provided: %d",
568567 typeNum
569568 );
570569 Py_CLEAR(str);
571570 return NULL;
572571 }
573572 }
574573
575574 int verbose = 0;
576575 if(self->nd>=0 && CudaNdarray_SIZE(self)==0){
577576 npy_intp * npydims = (npy_intp*)malloc(self->nd * sizeof(npy_intp));
578577 assert (npydims);
579578 for (int i = 0; i < self->nd; ++i) npydims[i] = (npy_intp)(CudaNdarray_HOST_DIMS(self)[i]);
580579 PyObject * rval = PyArray_SimpleNew(self->nd, npydims, REAL_TYPENUM);
581580 free(npydims);
582581 if (!rval){
583582 return NULL;
584583 }
585584 assert (PyArray_ITEMSIZE((PyArrayObject *)rval) == sizeof(real));
586585 return rval;
587586 }
588587 if ((self->nd < 0) || (self->devdata == 0))
589588 {
590589 PyErr_SetString(PyExc_ValueError, "can't copy from un-initialized CudaNdarray");
591590 return NULL;
592591 }
593592 CudaNdarray * contiguous_self = NULL;
594593 if (CudaNdarray_is_c_contiguous(self))
595594 {
596595 contiguous_self = self;
597596 Py_INCREF(contiguous_self);
598597 if (verbose) std::cerr << "CreateArrayObj already contiguous" << contiguous_self << '\n';
599598 }
600599 else
601600 {
602601 contiguous_self = (CudaNdarray*)CudaNdarray_Copy(self);
603602 if (verbose) std::cerr << "CreateArrayObj created contiguous" << contiguous_self << '\n';
604603 }
605604 if (!contiguous_self)
606605 {
607606 return NULL;
608607 }
609608
610609 npy_intp * npydims = (npy_intp*)malloc(self->nd * sizeof(npy_intp));
611610 assert (npydims);
612611 for (int i = 0; i < self->nd; ++i)
613612 npydims[i] = (npy_intp)(CudaNdarray_HOST_DIMS(self)[i]);
614613 PyArrayObject * rval = (PyArrayObject *) PyArray_SimpleNew(self->nd,
615614 npydims,
616615 REAL_TYPENUM);
617616 free(npydims);
618617 if (!rval)
619618 {
620619 Py_DECREF(contiguous_self);
621620 return NULL;
622621 }
623622
624623 assert (PyArray_ITEMSIZE(rval) == sizeof(real));
625624
626625 npy_intp rval_size = PyArray_SIZE(rval);
627626 void *rval_data = PyArray_DATA(rval);
628627 cudaError_t err;
629628 CNDA_BEGIN_ALLOW_THREADS;
630629
631630 err = cudaMemcpy(rval_data, contiguous_self->devdata,
632631 rval_size * sizeof(real),
633632 cudaMemcpyDeviceToHost
634633 );
635634 //CNDA_THREAD_SYNC; // unneeded because cudaMemcpy is blocking anyway
636635 CNDA_END_ALLOW_THREADS;
637636
638637 if (cudaSuccess != err)
639638 {
640639 PyErr_Format(PyExc_RuntimeError, "error (%s)copying data to host",
641640 cudaGetErrorString(err));
642641 Py_DECREF(rval);
643642 rval = NULL;
644643 }
645644
646645 Py_DECREF(contiguous_self);
647646 return (PyObject *)rval;
648647 }
649648
650649 // TODO-- we have two functions here, ZEROS and Zeros.
651650 // ZEROS is meant to be called just from C code (you don't need to pass it PyObject * s)
652651 // but this naming is very weird, makes it look like a macro
653652 // we should figure out the correct convention and change to that
654653 PyObject* CudaNdarray_ZEROS(int n, int * dims)
655654 {
656655
657656 size_t total_elements = 1;
658657
659658 for(size_t i=0;i<n;i++){
660659 // Detect overflow on unsigned integer
661660 if (dims[i] != 0 && total_elements > (SIZE_MAX / dims[i])) {
662661 PyErr_Format(PyExc_RuntimeError,
663662 "Can't store in size_t for the bytes requested %llu * %llu",
664663 (unsigned long long)total_elements,
665664 (unsigned long long)dims[i]);
666665 return NULL;
667666 }
668667 total_elements*=dims[i];
669668 }
670669
671670 // total_elements now contains the size of the array, in reals
672671 if (total_elements > (SIZE_MAX / sizeof(real))){
673672 PyErr_Format(PyExc_RuntimeError,
674673 "Can't store in size_t for the bytes requested %llu * 4",
675674 (unsigned long long)total_elements);
676675 return NULL;
677676 }
678677 size_t total_size = total_elements * sizeof(real);
679678
680679 CudaNdarray* rval = (CudaNdarray*)CudaNdarray_New();
681680 if (!rval)
682681 {
683682 PyErr_SetString(PyExc_RuntimeError, "CudaNdarray_ZEROS: call to New failed");
684683 return NULL;
685684 }
686685
687686 if (CudaNdarray_alloc_contiguous(rval, n, dims))
688687 {
689688 PyErr_SetString(PyExc_RuntimeError, "CudaNdarray_ZEROS: allocation failed.");
690689 Py_DECREF(rval);
691690 return NULL;
692691 }
693692
694693 // Fill with zeros
695694 //fprintf(stdout, "Sizeof: %d\n", total_size);
696695 if (cudaSuccess != cudaMemset(rval->devdata, 0, total_size))
697696 {
698697 PyErr_Format(PyExc_MemoryError,
699698 "CudaNdarray_ZEROS: Error memsetting %llu bytes of device memory.",
700699 (unsigned long long)total_size);
701700 Py_DECREF(rval);
702701 return NULL;
703702 }
704703
705704 if (cnda_copy_structure_to_device(rval))
706705 {
707706 PyErr_SetString(PyExc_RuntimeError, "CudaNdarray_ZEROS: syncing structure to device failed");
708707 Py_DECREF(rval);
709708 return NULL;
710709 }
711710 return (PyObject*) rval;
712711 }
713712
714713 // declared as a static method (hence 1st parameter is not used)
715714 // Based on _Copy and _dimshuffle
716715 PyObject* CudaNdarray_Zeros(PyObject* _unused, PyObject* shape)
717716 {
718717 if(!shape)
719718 {
720719 PyErr_SetString(PyExc_TypeError, "CudaNdarray_Zeros: function takes at least 1 argument (0 given)");
721720 return NULL;
722721 }
723722 if(!PySequence_Check(shape))
724723 {
725724 PyErr_SetString(PyExc_TypeError, "shape argument must be a sequence");
726725 return NULL;
727726 }
728727
729728 int shplen = PySequence_Length(shape);
730729
731730 if (shplen == 0)
732731 {
733732 return CudaNdarray_ZEROS(0, NULL);
734733 }
735734
736735 int* newdims = (int *)malloc(sizeof(int) * shplen);
737736
738737 if (!newdims)
739738 {
740739 PyErr_SetString(PyExc_MemoryError,
741740 "CudaNdarray_Zeros: Failed to allocate temporary space");
742741 return NULL;
743742 }
744743
745744 // start from the end to compute strides
746745 for (int i = shplen-1; i >= 0; --i)
747746 {
748747 PyObject* shp_el_obj = PySequence_GetItem(shape, i);
749748 if(shp_el_obj == NULL)
750749 {
751750 // shouldn't happen since we checked length before...
752751 PyErr_SetString(PyExc_RuntimeError, "CudaNdarray_Zeros: Index out of bound in sequence");
753752 free(newdims);
754753 return NULL;
755754 }
756755
757756 int shp_el = PyInt_AsLong(shp_el_obj);
758757 Py_DECREF(shp_el_obj);
759758
760759 if (shp_el < 0)
761760 {
762761 PyErr_SetString(PyExc_ValueError, "CudaNdarray_Zeros: shape must contain only non-negative values for size of a dimension");
763762 free(newdims);
764763 return NULL;
765764 }
766765
767766 newdims[i] = shp_el;
768767 }
769768
770769 PyObject* rval = CudaNdarray_ZEROS(shplen,newdims);
771770
772771 free(newdims);
773772
774773 return (PyObject*)rval;
775774 }
776775
777776
778777
779778
780779
781780 PyObject * CudaNdarray_Copy(const CudaNdarray * self)
782781 {
783782 PyObject * rval = CudaNdarray_New();
784783 if ((!rval) || (-1 == self->nd))
785784 {
786785 return rval;
787786 }
788787 if (CudaNdarray_alloc_contiguous((CudaNdarray*)rval, self->nd, CudaNdarray_HOST_DIMS(self)))
789788 {
790789 Py_DECREF(rval);
791790 return NULL;
792791 }
793792 if (CudaNdarray_CopyFromCudaNdarray((CudaNdarray*)rval, self))
794793 {
795794 Py_DECREF(rval);
796795 return NULL;
797796 }
798797 return rval;
799798 }
800799 PyObject * CudaNdarray_DeepCopy(CudaNdarray * self, PyObject * memo)
801800 {
802801 assert(PyDict_Check(memo));
803802 PyObject * selfkey = PyInt_FromLong((long)self);
804803 assert(selfkey);
805804 if (PyDict_Contains(memo, selfkey))
806805 {
807806 PyObject * rval = PyDict_GetItem(memo, selfkey);
808807 Py_DECREF(selfkey);
809808 Py_XINCREF(rval);
810809 return rval;
811810 }
812811 else
813812 {
814813 PyObject * rval = CudaNdarray_Copy(self);
815814 if (0) std::cerr << "DeepCopy created " << rval << " devdata " << ((CudaNdarray*)rval)->devdata << "\n";
816815 if (NULL == rval)
817816 {
818817 Py_DECREF(selfkey);
819818 return NULL;
820819 }
821820 if (PyDict_SetItem(memo, selfkey, rval))
822821 {
823822 Py_DECREF(rval);
824823 Py_DECREF(selfkey);
825824 return NULL;
826825 }
827826 Py_DECREF(selfkey);
828827 return rval;
829828 }
830829 }
831830 PyObject * CudaNdarray_ReduceSum(CudaNdarray * self, PyObject * py_reduce_mask)
832831 {
833832 if (!PySequence_Check(py_reduce_mask))
834833 {
835834 PyErr_SetString(PyExc_TypeError, "reduce_mask must be sequence of ints");
836835 return NULL;
837836 }
838837 int len = PySequence_Length(py_reduce_mask);
839838 if (len != self->nd)
840839 {
841840 PyErr_SetString(PyExc_TypeError, "length of reduce_mask must match self->nd");
842841 return NULL;
843842 }
844843 CudaNdarray * self_sum = (CudaNdarray*)CudaNdarray_New();
845844 if (!self_sum)
846845 {
847846 return NULL;
848847 }
849848 //TODO: allocate a fixed size dimshuffle_pattern_cache on the stack,
850849 // and use it if it is big enough.
851850 int * dimshuffle_pattern = (int*)malloc(len * 2 * sizeof(int));
852851 int * sum_dims = dimshuffle_pattern + len;
853852 int n_remaining_dims = 0;
854853 if (!dimshuffle_pattern)
855854 {
856855 Py_DECREF(self_sum);
857856 PyErr_SetString(PyExc_MemoryError, "failed to alloc internal storage");
858857 return NULL;
859858 }
860859 for (int i = 0; i < len; ++i)
861860 {
862861 PyObject *o_i = PySequence_GetItem(py_reduce_mask, i);
863862 int o_i_int = PyInt_AsLong(o_i);
864863 Py_XDECREF(o_i);
865864 if (PyErr_Occurred())
866865 {
867866 Py_DECREF(self_sum);
868867 free(dimshuffle_pattern);
869868 return NULL;
870869 }
871870 if (o_i_int) // this is a dimension over which we are reducing
872871 {
873872 sum_dims[i] = 1;
874873 }
875874 else
876875 {
877876 sum_dims[i] = CudaNdarray_HOST_DIMS(self)[i];
878877 dimshuffle_pattern[n_remaining_dims++] = i;
879878 }
880879 }
881880 if (0 || CudaNdarray_alloc_contiguous(self_sum, len, sum_dims)
882881 || CudaNdarray_reduce_sum(self_sum, self)
883882 || CudaNdarray_dimshuffle(self_sum, n_remaining_dims, dimshuffle_pattern))
884883 {
885884 Py_DECREF(self_sum);
886885 free(dimshuffle_pattern);
887886 return NULL;
888887 }
889888 free(dimshuffle_pattern);
890889 return (PyObject*)self_sum;
891890 }
892891
893892 // Reshape self to the new shape gived by the tuple shape.
894893 //
895894 // If self is c contiguous, it return a view. Otherwise it always do a copy.
896895 // TODO: make it return a view when the strides allow it even if it is not
897896 // c contiguous
898897 PyObject * CudaNdarray_Reshape(CudaNdarray * self, PyObject * shape)
899898 {
900899 if(!CudaNdarray_is_c_contiguous(self))
901900 {
902901 // allocate new space
903902 //TODO: test to see if we can re-use old one and take a new param to
904903 // use this
905904 CudaNdarray* rval = (CudaNdarray*) CudaNdarray_Copy(self);
906905 if (!rval)
907906 {
908907 return NULL;
909908 }
910909
911910 CudaNdarray* ret = (CudaNdarray*) CudaNdarray_Reshape(rval, shape);
912911 Py_XDECREF(rval);
913912 return (PyObject*)ret;
914913 }
915914
916915 // check shape tuple
917916 unsigned int rval_nd;
918917 unsigned int * rval_dims;
919918 size_t rval_size = 1;
920919
921920 if (PyTuple_Check(shape)){
922921 // copy shape to integer array
923922 rval_nd = PyTuple_Size(shape);
924923 }else if (PyInt_Check(shape)){
925924 rval_nd = 1;
926925 }else{
927926 PyErr_SetString(PyExc_TypeError, "shape must be tuple of integers or an integer");
928927 return NULL;
929928 }
930929 rval_dims = (unsigned int*)malloc(rval_nd * sizeof(int));
931930
932931 if(PyTuple_Check(shape)){
933932 for (int i = 0; i < rval_nd; ++i)
934933 {
935934 rval_dims[i] = PyInt_AsLong(PyTuple_GetItem(shape, i)); //GetItem returns borrowed reference
936935 if (PyErr_Occurred()) //error in AsLong
937936 {
938937 free(rval_dims);
939938 return NULL;
940939 }
941940 if(rval_dims[i]<0){
942941 PyErr_Format(PyExc_ValueError, "Reshape has invalid dimension %i (must be >=0)",rval_dims[i]);
943942 free(rval_dims);
944943 return NULL;
945944 }
946945 rval_size = rval_size * rval_dims[i];
947946 }
948947 }else{
949948 rval_size = PyInt_AsLong(shape);
950949 rval_dims[0] = rval_size;
951950 }
952951 // calculate new size, assert same as old size
953952 if (rval_size != CudaNdarray_SIZE(self))
954953 {
955954 PyErr_Format(PyExc_ValueError, "size must remain unchanged, changed from %lld to %lld", CudaNdarray_SIZE(self), rval_size);
956955 free(rval_dims);
957956 return NULL;
958957 }
959958 if (rval_size==0)
960959 {
961960 PyObject * rval = CudaNdarray_NewDims(rval_nd, rval_dims);
962961 free(rval_dims);
963962 return rval;
964963 }
965964
966965 //return a view, not a copy
967966 //we can do this as we checked self is c_contiguous
968967 CudaNdarray * rval = (CudaNdarray * )CudaNdarray_New(rval_nd);
969968
970969 if (!rval || 0 != rval->data_allocated
971970 ||CudaNdarray_set_device_data(rval, CudaNdarray_DEV_DATA(self), self))
972971 {
973972 Py_XDECREF(rval);
974973 free(rval_dims);
975974 return NULL;
976975 }
977976 //set dim and stride
978977 int size = 1;
979978 for (int i = rval_nd-1; i >= 0; --i)
980979 {
981980 CudaNdarray_set_stride(rval, i, (rval_dims[i] == 1) ? 0 : size);
982981 CudaNdarray_set_dim(rval, i, rval_dims[i]);
983982 size = size * rval_dims[i];
984983 }
985984 free(rval_dims);
986985 return (PyObject*)rval;
987986 }
988987
989988 PyObject * CudaNdarray_View(const CudaNdarray * self)
990989 {
991990 CudaNdarray * rval = (CudaNdarray*)CudaNdarray_New(self->nd);
992991 if (!rval || CudaNdarray_set_device_data(rval, CudaNdarray_DEV_DATA(self), self))
993992 {
994993 Py_XDECREF(rval);
995994 rval = NULL;
996995 }
997996 else
998997 {
999998 for (int i = 0; i < self->nd; ++i)
1000999 {
10011000 CudaNdarray_set_dim(rval, i, CudaNdarray_HOST_DIMS(self)[i]);
10021001 CudaNdarray_set_stride(rval, i, CudaNdarray_HOST_STRIDES(self)[i]);
10031002 }
10041003 }
10051004 return (PyObject*)rval;
10061005 }
10071006
10081007 /*
10091008 * d0,... are the output dims
10101009 * indices are a list of index to operate on
10111010 * They are int32 viewed as float32.
10121011 * a is the output
10131012 * b is the input
10141013 * dB0, the source leading dimensions size
10151014 */
10161015 template <int operator_num>
10171016 __global__ void k_take_3(const int d0, const int d1, const int d2,
10181017 const npy_int64* indices,
10191018 float* a,
10201019 const int sA0, const int sA1, const int sA2,
10211020 const float* b, const int dB0,
10221021 const int sB0, const int sB1, const int sB2,
10231022 int* err){
10241023 for (int i0 = blockIdx.x; i0 < d0; i0 += gridDim.x){
10251024 npy_int64 idx = indices[i0];
10261025 if (idx<0)
10271026 idx += dB0; // To allow negative indexing.
10281027 if ((idx < 0) || (idx >= dB0)){
10291028 // Any value other the 0 probably work. But to be more safe, I want
10301029 // to change all bits to prevent problem with concurrent write that
10311030 // could cross cache line. But this should not happen with the
10321031 // current code and driver.
10331032 *err = 0xFFFF;
10341033 continue;
10351034 }
10361035 for (int i1 = threadIdx.x; i1 < d1; i1 += blockDim.x){
10371036 for (int i2 = threadIdx.y; i2 < d2; i2 += blockDim.y){
10381037 int a_idx = i0*sA0 + i1*sA1 + i2*sA2;
10391038 int b_idx = idx*sB0 + i1*sB1 + i2*sB2;
10401039 a[a_idx] = b[b_idx];
10411040 }
10421041 }
10431042 }
10441043 }
10451044
10461045 // We try to be similar to the PyArray_TakeFrom function
10471046 //http://docs.scipy.org/doc/numpy/reference/c-api.array.html
10481047 //TODO: support other clip mode then raise(clip, wrap)
10491048 //self is the input that we copy data from.
10501049 //The indices that we receive MUST be an CudaNdarray(float32)
10511050 // that is in fact a view to int64 indices
10521051 PyObject*
10531052 CudaNdarray_TakeFrom(CudaNdarray * self, PyObject *args){
10541053 int verbose = 0;
10551054 PyObject * indices_obj = NULL;
10561055 //int axis; Default None, that mean the flattened array.
10571056 PyObject * axis_obj = Py_None;
10581057 PyObject * out_obj = Py_None;
10591058 PyObject * clipmode_obj = NULL;
10601059 int max_threads = 1; // max threads per blocks
10611060
10621061 if (! PyArg_ParseTuple(args, "O|OOOi", &indices_obj, &axis_obj,
10631062 &out_obj, &clipmode_obj, &max_threads))
10641063 return NULL;
10651064
10661065 //Check argument indices
10671066 //TODO: if not a numpy.ndarray, convert to numpy.ndarray
10681067 //TODO: If a CudaNdarray, accept it and suppose the data is int32? is float32 number of int?
10691068 //TODO: Support ndarray of other dtype then int32
10701069 //TODO: support list of indices that are not c_contiguous
10711070 CudaNdarray * indices = NULL;
10721071 if (CudaNdarray_Check(indices_obj)) {
10731072 if (verbose) printf("cudandarray indices\n");
10741073 indices = (CudaNdarray*) indices_obj;
10751074 Py_INCREF(indices);
10761075 } else if (PyArray_Check(indices_obj)) {
10771076 if (verbose) printf("ndarray indices\n");
10781077 if (PyArray_TYPE((PyArrayObject *)indices_obj) != NPY_INT64) {
10791078 PyErr_SetString(PyExc_TypeError,
10801079 "CudaNdarray_TakeFrom: need a ndarray for indices"
10811080 " with dtype int64");
10821081 return NULL;
10831082 }
10841083 if (PyArray_NDIM(((PyArrayObject*)indices_obj)) != 1) {
10851084 PyErr_SetString(PyExc_TypeError,
10861085 "CudaNdarray_TakeFrom: need a CudaNdarray of"
10871086 " indices with only 1 dimensions");
10881087 return NULL;
10891088 }
10901089 // We need indices_obj to be contiguous, in order to take a view
10911090 // with a different dtype.
10921091 if (!PyArray_IS_C_CONTIGUOUS((PyArrayObject*) indices_obj)) {
10931092 PyObject* indices_obj_contig = PyArray_NewCopy((PyArrayObject*) indices_obj, NPY_CORDER);
10941093 if (!indices_obj_contig)
10951094 return NULL;
10961095 indices_obj = indices_obj_contig;
10971096 } else {
10981097 // Keep the refcount consistent
10991098 Py_INCREF(indices_obj);
11001099 }
11011100 PyArray_Descr* float32_descr = PyArray_DescrFromType(NPY_FLOAT32);
11021101 PyObject * indices_float32 = NULL;
11031102 indices_float32 = PyArray_View((PyArrayObject*)indices_obj,
11041103 float32_descr, NULL);
11051104 if (verbose) printf("ndarray indices\n");
11061105 if (!indices_float32) {
11071106 Py_DECREF(indices_obj);
11081107 return NULL;
11091108 }
11101109
11111110 indices = (CudaNdarray*) CudaNdarray_New();
11121111 if (verbose) printf("\nndarray after new\n");
11131112 if (! indices){
11141113 Py_DECREF(indices_obj);
11151114 Py_DECREF(indices_float32);
11161115 return NULL;
11171116 }
11181117 if (CudaNdarray_CopyFromArray(indices,
11191118 (PyArrayObject *)indices_float32)){
11201119 Py_DECREF(indices_obj);
11211120 Py_DECREF(indices_float32);
11221121 return NULL;
11231122 }
11241123 Py_DECREF(indices_obj);
11251124 Py_DECREF(indices_float32);
11261125 } else {
11271126 PyObject* py_s = PyObject_Str(indices_obj);
11281127 const char* s = PyString_AsString(py_s);
11291128 Py_DECREF(py_s);
11301129 PyErr_Format(PyExc_TypeError,
11311130 "CudaNdarray_TakeFrom: need an ndarray of int64 or a"
11321131 " CudaNdarray(float32) that is a view from int64 data"
11331132 " for indices. Got %s", s);
11341133 return NULL;
11351134 }
11361135
11371136 if (verbose) {
11381137 printf("indices used on the gpu\n");
11391138 fprint_CudaNdarray(stdout, indices);
11401139 PyObject * used_indices = CudaNdarray_CreateArrayObj(indices);
11411140 PyObject_Print(used_indices, stdout, 0);
11421141 Py_DECREF(used_indices);
11431142 }
11441143 if (verbose) printf("after print of object\n");
11451144 if(!CudaNdarray_is_c_contiguous(indices) != 0) {
11461145 PyErr_SetString(PyExc_NotImplementedError,
11471146 "CudaNdarray_TakeFrom: The indices must be contiguous in memory.");
11481147 Py_DECREF(indices);
11491148 return NULL;
11501149 }
11511150 int nb_indices = CudaNdarray_SIZE((CudaNdarray *)indices) / 2;// int64 are 8 bytes, float32 are 4 bytes
11521151
11531152 //Check argument axis
11541153 //TODO: implement the default and other axis
11551154 long axis = PyInt_AsLong(axis_obj);
11561155
11571156 if (axis != 0) {
11581157 PyErr_Format(PyExc_NotImplementedError,
11591158 "CudaNdarray_TakeFrom: only axis=0 is currently supported."
11601159 " Got %ld.", axis);
11611160 Py_DECREF(indices);
11621161 return NULL;
11631162 }
11641163
11651164 //Check argument out_obj
11661165 CudaNdarray * out = NULL;
11671166 if (out_obj && CudaNdarray_Check(out_obj))
11681167 out = (CudaNdarray*) out_obj;
11691168 if (out && (out->nd != self->nd ||
11701169 CudaNdarray_HOST_DIMS(out)[0] != nb_indices))
11711170 out = NULL;
11721171 int * dims = (int *)malloc(sizeof(int) * self->nd);
11731172 dims[0] = nb_indices;
11741173
11751174 for (int i=1 ; i<self->nd ; i++) {
11761175 dims[i] = CudaNdarray_HOST_DIMS(self)[i];
11771176 if (out && CudaNdarray_HOST_DIMS(out)[i] != dims[i]) {
11781177 out = NULL;
11791178 }
11801179 }
11811180 if (!out) {
11821181 out = (CudaNdarray*)CudaNdarray_New();
11831182 if (!out){
11841183 Py_DECREF(indices);
11851184 free(dims);
11861185 return NULL;
11871186 }
11881187 if (CudaNdarray_alloc_contiguous(out, self->nd, dims)) {
11891188 Py_DECREF(out);
11901189 Py_DECREF(indices);
11911190 free(dims);
11921191 return NULL;
11931192 }
11941193 }else {
11951194 Py_INCREF(out);
11961195 }
11971196
11981197 //Check argument clipmode
11991198 if (clipmode_obj) {
12001199 char * clipmode = PyString_AsString(clipmode_obj);
12011200 if (! clipmode){
12021201 Py_DECREF(indices);
12031202 Py_DECREF(out);
12041203 free(dims);
12051204 return NULL;
12061205 }
12071206 if (strcmp(clipmode, "raise") != 0) {
12081207 PyErr_Format(PyExc_NotImplementedError,
12091208 "CudaNdarray_TakeFrom: only the raise mode is currently supported. Got '%s'",
12101209 clipmode);
12111210 Py_DECREF(indices);
12121211 Py_DECREF(out);
12131212 free(dims);
12141213 return NULL;
12151214 }
12161215 }
12171216 void (*k3)(const int, const int, const int,
12181217 const npy_int64*,
12191218 float*, const int, const int, const int,
12201219 const float*, const int,
12211220 const int, const int, const int,
12221221 int*);
12231222 k3 = k_take_3<CPY>;
12241223
12251224 // Create the memory place that will store the error information.
12261225 if(init_err_var() != 0) return NULL;
12271226
12281227 dim3 n_blocks(std::min(CudaNdarray_HOST_DIMS(out)[0],65535),1,1);
12291228 if(CudaNdarray_HOST_DIMS(out)[0] == 0){
12301229 // We take 0 elements, so no need for the rest of the code.
12311230 // This speed up that case AND fix crash otherwise.
12321231 free(dims);
12331232 Py_DECREF(indices);
12341233 return (PyObject *)out;
12351234 }
12361235
12371236 switch (self->nd) {
12381237 case 1:
12391238 {
12401239 dim3 n_threads(1, 1, 1);
12411240 if (verbose)
12421241 printf("cudaGetLastError=%d, nd=%d"
12431242 " kernel config: (n_blocks.x=%d, n_blocks.y=%d,"
12441243 " n_threads.x=%i, n_threads.y=%i)\n",
12451244 cudaGetLastError(), self->nd,
12461245 n_blocks.x, n_blocks.y, n_threads.x, n_threads.y);
12471246 k3<<<n_blocks, n_threads>>>(
12481247 dims[0],
12491248 1,
12501249 1,
12511250 (npy_int64*) CudaNdarray_DEV_DATA(indices),
12521251 CudaNdarray_DEV_DATA(out),
12531252 CudaNdarray_HOST_STRIDES(out)[0], //strides
12541253 1,
12551254 1,
12561255 CudaNdarray_DEV_DATA(self),
12571256 CudaNdarray_HOST_DIMS(self)[0], //For indices check
12581257 CudaNdarray_HOST_STRIDES(self)[0], //strides
12591258 1,
12601259 1,
12611260 err_var);
12621261 }
12631262 break;
12641263 case 2:
12651264 {
12661265 dim3 n_threads(std::min(CudaNdarray_HOST_DIMS(out)[1], max_threads), 1, 1);
12671266
12681267 if (verbose)
12691268 printf("cudaGetLastError=%d, nd=%d"
12701269 " kernel config: (n_blocks.x=%d, n_blocks.y=%d,"
12711270 " n_threads.x=%i, n_threads.y=%i)\n",
12721271 cudaGetLastError(), self->nd,
12731272 n_blocks.x, n_blocks.y, n_threads.x, n_threads.y);
12741273
12751274 k3<<<n_blocks, n_threads>>>(
12761275 dims[0], //dimensions
12771276 dims[1],
12781277 1,
12791278 (npy_int64*) CudaNdarray_DEV_DATA(indices),
12801279 CudaNdarray_DEV_DATA(out),
12811280 CudaNdarray_HOST_STRIDES(out)[0], //strides
12821281 CudaNdarray_HOST_STRIDES(out)[1],
12831282 1,
12841283 CudaNdarray_DEV_DATA(self),
12851284 CudaNdarray_HOST_DIMS(self)[0], //For indices check
12861285 CudaNdarray_HOST_STRIDES(self)[0], //strides
12871286 CudaNdarray_HOST_STRIDES(self)[1],
12881287 1,
12891288 err_var);
12901289 }
12911290 break;
12921291 case 3:
12931292 {
12941293 int ty = std::min(CudaNdarray_HOST_DIMS(out)[2], max_threads);
12951294 int tx = std::min(CudaNdarray_HOST_DIMS(out)[1], max_threads / ty);
12961295 dim3 n_threads(tx, ty, 1);
12971296 if (verbose)
12981297 printf("cudaGetLastError=%d, nd=%d"
12991298 " kernel config: (n_blocks.x=%d, n_blocks.y=%d,"
13001299 " n_threads.x=%i, n_threads.y=%i)\n",
13011300 cudaGetLastError(), self->nd,
13021301 n_blocks.x, n_blocks.y, n_threads.x, n_threads.y);
13031302 k3<<<n_blocks, n_threads>>>(
13041303 dims[0], //dimensions
13051304 dims[1],
13061305 dims[2],
13071306 (npy_int64*) CudaNdarray_DEV_DATA(indices),
13081307 CudaNdarray_DEV_DATA(out),
13091308 CudaNdarray_HOST_STRIDES(out)[0], //strides
13101309 CudaNdarray_HOST_STRIDES(out)[1],
13111310 CudaNdarray_HOST_STRIDES(out)[2],
13121311 CudaNdarray_DEV_DATA(self),
13131312 CudaNdarray_HOST_DIMS(self)[0], //For indices check
13141313 CudaNdarray_HOST_STRIDES(self)[0], //strides
13151314 CudaNdarray_HOST_STRIDES(self)[1],
13161315 CudaNdarray_HOST_STRIDES(self)[2],
13171316 err_var);
13181317 }
13191318 break;
13201319 default:
13211320 PyErr_SetString(PyExc_NotImplementedError,
13221321 "CudaNdarray_TakeFrom: only input with 1, 2 or 3"
13231322 " dimensions are currently supported");
13241323
13251324 }
13261325 free(dims);
13271326 CNDA_THREAD_SYNC;
13281327 cudaError_t err = cudaGetLastError();
13291328 if (cudaSuccess != err) {
13301329 PyErr_Format(PyExc_RuntimeError,
13311330 "Cuda error: %s: %s.\n",
13321331 "CudaNdarray_TakeFrom",
13331332 cudaGetErrorString(err));
13341333 Py_DECREF(indices);
13351334 Py_DECREF(out);
13361335 return NULL;
13371336 }
13381337
13391338 int index_err = check_err_var();
13401339 Py_DECREF(indices);
13411340 if (index_err != 0) {
13421341 Py_DECREF(out);
13431342 return NULL;
13441343 }
13451344
13461345 if (verbose) printf("TAKE SUCCEDED\n");
13471346 return (PyObject *)out;
13481347 }
13491348
13501349
13511350 PyObject * CudaNdarray_SetStride(CudaNdarray * self, PyObject *args)
13521351 {
13531352 int pos, stride;
13541353 if (! PyArg_ParseTuple(args, "ii", &pos, &stride))
13551354 return NULL;
13561355 if ((pos < 0) || (pos >= self->nd))
13571356 {
13581357 PyErr_Format(PyExc_ValueError, "position argument out of legal range [0, %i)", self->nd);
13591358 return NULL;
13601359 }
13611360 CudaNdarray_set_stride(self, pos, stride);
13621361 if (cnda_copy_structure_to_device(self))
13631362 {
13641363 return NULL;
13651364 }
13661365 Py_INCREF(Py_None);
13671366 return Py_None;
13681367 }
13691368 PyObject * CudaNdarray_SetShapeI(CudaNdarray * self, PyObject *args)
13701369 {
13711370 int pos, dim;
13721371 if (! PyArg_ParseTuple(args, "ii", &pos, &dim))
13731372 return NULL;
13741373 if ((pos < 0) || (pos >= self->nd))
13751374 {
13761375 PyErr_Format(PyExc_ValueError, "position argument out of legal range [0, %i)", self->nd);
13771376 return NULL;
13781377 }
13791378 CudaNdarray_set_dim(self, pos, dim);
13801379 if (cnda_copy_structure_to_device(self))
13811380 {
13821381 return NULL;
13831382 }
13841383 Py_INCREF(Py_None);
13851384 return Py_None;
13861385 }
13871386
13881387 static PyObject *
13891388 CudaNdarray_exp(CudaNdarray* self)
13901389 {
13911390 CudaNdarray * rval = (CudaNdarray *)CudaNdarray_New();
13921391 if ((NULL == rval) || CudaNdarray_alloc_contiguous(rval, self->nd, CudaNdarray_HOST_DIMS(self)))
13931392 {
13941393 Py_XDECREF(rval);
13951394 return NULL;
13961395 }
13971396 unsigned int size = 1;
13981397 for (int i = 0; i < self->nd; i++)
13991398 {
14001399 size *= (unsigned int) CudaNdarray_HOST_DIMS(self)[i];
14011400 }
14021401 unsigned int threads_per_block = std::min(size, (unsigned int)NUM_VECTOR_OP_THREADS_PER_BLOCK);
14031402 unsigned int n_blocks = std::min(ceil_intdiv(size,threads_per_block), (unsigned int)NUM_VECTOR_OP_BLOCKS);
14041403 k_elemwise_unary_rowmajor_exp<<<n_blocks,threads_per_block>>>(size, self->nd, CudaNdarray_DEV_DIMS(self),
14051404 CudaNdarray_DEV_DATA(self), CudaNdarray_DEV_STRIDES(self),
14061405 CudaNdarray_DEV_DATA(rval), CudaNdarray_DEV_STRIDES(rval));
14071406
14081407 //TODO: don't do this right away, do it when we need the result
14091408 CNDA_THREAD_SYNC;
14101409 cudaError_t err = cudaGetLastError();
14111410 if( cudaSuccess != err)
14121411 {
14131412 Py_DECREF(rval);
14141413 PyErr_Format(PyExc_RuntimeError, "Cuda error: %s: %s.\n", "kExp", cudaGetErrorString(err));
14151414 return NULL;
14161415 }
14171416
14181417 return (PyObject*)rval;
14191418 }
14201419
14211420 static PyMethodDef CudaNdarray_methods[] =
14221421 {
14231422 {"__array__",
14241423 (PyCFunction)CudaNdarray_CreateArrayObj, METH_VARARGS,
14251424 "Copy from the device to a numpy ndarray"},
14261425 {"__copy__",
14271426 (PyCFunction)CudaNdarray_View, METH_NOARGS,
14281427 "Create a shallow copy of this object. used by module copy"},
14291428 {"__deepcopy__",
14301429 (PyCFunction)CudaNdarray_DeepCopy, METH_O,
14311430 "Create a copy of this object"},
14321431 {"zeros",
14331432 (PyCFunction)CudaNdarray_Zeros, METH_STATIC | METH_O,
14341433 "Create a new CudaNdarray with specified shape, filled with zeros."},
14351434 {"copy",
14361435 (PyCFunction)CudaNdarray_Copy, METH_NOARGS,
14371436 "Create a copy of this object"},
14381437 {"is_c_contiguous",
14391438 (PyCFunction)CudaNdarray_IS_C_Contiguous, METH_NOARGS,
14401439 "Return True is the object is c contiguous. False otherwise."},
14411440 {"reduce_sum",
14421441 (PyCFunction)CudaNdarray_ReduceSum, METH_O,
14431442 "Reduce over the given dimensions by summation"},
14441443 {"exp",
14451444 (PyCFunction)CudaNdarray_exp, METH_NOARGS,
14461445 "Return the exponential of all elements"},
14471446 {"reshape",
14481447 (PyCFunction)CudaNdarray_Reshape, METH_O,
14491448 "Return a reshaped view (or copy) of this ndarray\n\
14501449 The required argument is a tuple of integers specifying the shape of the new ndarray."},
14511450 {"view",
14521451 (PyCFunction)CudaNdarray_View, METH_NOARGS,
14531452 "Return an alias of this ndarray"},
14541453 {"_set_stride",
14551454 (PyCFunction)CudaNdarray_SetStride, METH_VARARGS,
14561455 "For integer arguments (i, s), set the 'i'th stride to 's'"},
14571456 {"take",
14581457 (PyCFunction)CudaNdarray_TakeFrom, METH_VARARGS,
14591458 "Equivalent of numpy.take"},
14601459 {"_set_shape_i",
14611460 (PyCFunction)CudaNdarray_SetShapeI, METH_VARARGS,
14621461 "For integer arguments (i, s), set the 'i'th shape to 's'"},
14631462 {NULL, NULL, NULL, NULL} /* Sentinel */
14641463 };
14651464
14661465
14671466 ////////////////////
14681467 // Number protocol
14691468 ////////////////////
14701469
14711470 __global__ void kAdd_contiguous(float* a, float* b, float* dest, unsigned int numEls) {
14721471 const unsigned int idx = blockIdx.x * blockDim.x + threadIdx.x;
14731472 const unsigned int numThreads = blockDim.x * gridDim.x;
14741473
14751474 for (unsigned int i = idx; i < numEls; i += numThreads) {
14761475 dest[i] = a[i] + b[i];
14771476 }
14781477 }
14791478
14801479 // Will be called by __add__ in Python
14811480 static PyObject *
14821481 CudaNdarray_add(PyObject* py_self, PyObject * py_other)
14831482 {
14841483 if (! CudaNdarray_Check(py_self)) {
14851484 PyErr_SetString(PyExc_TypeError, "need a CudaNdarray on left");
14861485 return NULL;
14871486 }
14881487 if (! CudaNdarray_Check(py_other)) {
14891488 PyErr_SetString(PyExc_TypeError, "need a CudaNdarray on right");
14901489 return NULL;
14911490 }
14921491 CudaNdarray * self = (CudaNdarray *)py_self;
14931492 CudaNdarray * other = (CudaNdarray *)py_other;
14941493 if(!CudaNdarray_is_c_contiguous(self) || !CudaNdarray_is_c_contiguous(other)){
14951494 PyErr_SetString(PyExc_TypeError, "We have implementet only the c_contiguous version for now.");
14961495 return NULL;
14971496 }
14981497
14991498 //standard elemwise size checks
15001499 if (self->nd != other->nd)
15011500 {
15021501 PyErr_SetString(PyExc_TypeError, "CudaNdarray_add: need same number of dims");
15031502 return NULL;
15041503 }
15051504 //standard elemwise dim checks
15061505 unsigned int size = 1;
15071506 for (int i = 0; i< self->nd; ++i)
15081507 {
15091508 if (CudaNdarray_HOST_DIMS(self)[i] != CudaNdarray_HOST_DIMS(other)[i])
15101509 {
15111510 PyErr_SetString(PyExc_TypeError, "need same dimensions");
15121511 return NULL;
15131512 }
15141513 size *= (unsigned int) CudaNdarray_HOST_DIMS(self)[i];
15151514 }
15161515 CudaNdarray * rval = (CudaNdarray *)CudaNdarray_New();
15171516 if (!rval || CudaNdarray_alloc_contiguous(rval, self->nd, CudaNdarray_HOST_DIMS(self)))
15181517 {
15191518 Py_XDECREF(rval);
15201519 return NULL;
15211520 }
15221521
15231522 if(CudaNdarray_SIZE((CudaNdarray *)py_self)==0 && CudaNdarray_SIZE((CudaNdarray *)py_other)==0){
15241523 return (PyObject *) rval;
15251524 }
15261525
15271526 int threads_per_block = std::min(size, (unsigned int)NUM_VECTOR_OP_THREADS_PER_BLOCK);
15281527 int n_blocks = std::min(ceil_intdiv(size,(unsigned int)threads_per_block), (unsigned int)NUM_VECTOR_OP_BLOCKS);
15291528 kAdd_contiguous<<<n_blocks,threads_per_block>>>(
15301529 self->devdata, other->devdata, rval->devdata, size);
15311530 CNDA_THREAD_SYNC;
15321531 cudaError_t err = cudaGetLastError();
15331532 if( cudaSuccess != err)
15341533 {
15351534 PyErr_Format(PyExc_RuntimeError, "Cuda error: %s: %s.\n", "kAdd", cudaGetErrorString(err));
15361535 Py_DECREF(rval);
15371536 return NULL;
15381537 }
15391538 return (PyObject *) rval;
15401539 }
15411540
15421541 template <int operator_num>
15431542 __global__ void k_ielem_3(const int d0, const int d1, const int d2,
15441543 float* a, const int sA0, const int sA1, const int sA2,
15451544 const float* b, const int sB0, const int sB1, const int sB2){
15461545 for (int i0 = blockIdx.x; i0 < d0; i0 += gridDim.x){
15471546 for (int i1 = blockIdx.y; i1 < d1; i1 += gridDim.y){
15481547 for (int i2 = threadIdx.x; i2 < d2; i2 += blockDim.x){
15491548 switch (operator_num)
15501549 {
15511550 case IADD:
15521551 a[i0*sA0 + i1*sA1 + i2*sA2] += b[i0*sB0 + i1*sB1 + i2*sB2];
15531552 break;
15541553 case IDIV:
15551554 a[i0*sA0 + i1*sA1 + i2*sA2] /= b[i0*sB0 + i1*sB1 + i2*sB2];
15561555 break;
15571556 case CPY:
15581557 a[i0*sA0 + i1*sA1 + i2*sA2] = b[i0*sB0 + i1*sB1 + i2*sB2];
15591558 break;
15601559 }
15611560 }
15621561 }
15631562 }
15641563 }
15651564
15661565 template <int operator_num>
15671566 __global__ void k_ielem_4(const int d0, const int d1, const int d2, const int d3,
15681567 float* a, const int sA0, const int sA1,
15691568 const int sA2, const int sA3,
15701569 const float* b, const int sB0, const int sB1,
15711570 const int sB2, const int sB3){
15721571 for (int i0 = blockIdx.x; i0 < d0; i0 += gridDim.x){
15731572 for (int i1 = blockIdx.y; i1 < d1; i1 += gridDim.y){
15741573 for (int i2 = threadIdx.x; i2 < d2; i2 += blockDim.x){
15751574 for (int i3 = threadIdx.y; i3 < d3; i3 += blockDim.y){
15761575 switch (operator_num) {
15771576 case IADD:
15781577 a[i0*sA0 + i1*sA1 + i2*sA2 + i3*sA3]
15791578 += b[i0*sB0 + i1*sB1 + i2*sB2 + i3*sB3];
15801579 break;
15811580 case IDIV:
15821581 a[i0*sA0 + i1*sA1 + i2*sA2 + i3*sA3]
15831582 /= b[i0*sB0 + i1*sB1 + i2*sB2 + i3*sB3];
15841583 break;
15851584 case CPY:
15861585 a[i0*sA0 + i1*sA1 + i2*sA2 + i3*sA3]
15871586 = b[i0*sB0 + i1*sB1 + i2*sB2 + i3*sB3];
15881587 break;
15891588 }
15901589 }
15911590 }
15921591 }
15931592 }
15941593 }
15951594
15961595 template <int operator_num>
15971596 __global__ void k_ielem_6(const int d0, const int d1,
15981597 const int d2, const int d3,
15991598 const int d4, const int d5,
16001599 float* a, const int sA0, const int sA1,
16011600 const int sA2, const int sA3,
16021601 const int sA4, const int sA5,
16031602 const float* b, const int sB0, const int sB1,
16041603 const int sB2, const int sB3,
16051604 const int sB4, const int sB5
16061605 ){
16071606 for (int i0 = blockIdx.x; i0 < d0; i0 += gridDim.x){
16081607 for (int i1 = blockIdx.y; i1 < d1; i1 += gridDim.y){
16091608 for (int i2 = blockIdx.z; i2 < d2; i2 += gridDim.z){
16101609 for (int i3 = threadIdx.x; i3 < d3; i3 += blockDim.x){
16111610 for (int i4 = threadIdx.y; i4 < d4; i4 += blockDim.y){
16121611 for (int i5 = threadIdx.z; i5 < d5; i5 += blockDim.z){
16131612 switch (operator_num) {
16141613 case IADD:
16151614 a[i0*sA0 + i1*sA1 + i2*sA2 + i3*sA3 + i4*sA4 + i5*sA5]
16161615 += b[i0*sB0 + i1*sB1 + i2*sB2 + i3*sB3 + i4*sB4 + i5*sB5];
16171616 break;
16181617 case IDIV:
16191618 a[i0*sA0 + i1*sA1 + i2*sA2 + i3*sA3 + i4*sA4 + i5*sA5]
16201619 /= b[i0*sB0 + i1*sB1 + i2*sB2 + i3*sB3 + i4*sB4 + i5*sB5];
16211620 break;
16221621 case CPY:
16231622 a[i0*sA0 + i1*sA1 + i2*sA2 + i3*sA3 + i4*sA4 + i5*sA5]
16241623 = b[i0*sB0 + i1*sB1 + i2*sB2 + i3*sB3 + i4*sB4 + i5*sB5];
16251624 break;
16261625 }
16271626 }
16281627 }
16291628 }
16301629 }
16311630 }
16321631 }
16331632 }
16341633
16351634 /*
16361635 CudaNdarray_inplace_elemwise
16371636 Compute elemwise, working inplace on A.
16381637 Currently implemented A / B, A + B and A = B
16391638 (the last is not tested and not used!)
16401639
16411640 py_self - the CudaNdarray that we'll modify (A)
16421641 py_other - the other argument (B)
16431642 fct_nb - which operation to perform (operator_t)
16441643
16451644 Returns 0 on success.
16461645 Returns -1 on failure, and sets Python exception.
16471646
16481647 */
16491648 int
16501649 CudaNdarray_inplace_elemwise(PyObject* py_self, PyObject * py_other, operator_t fct_nb)
16511650 {
16521651 int verbose = 0;
16531652 void (*k3)(const int, const int, const int,
16541653 float*, const int, const int, const int,
16551654 const float*, const int, const int, const int);
16561655 void (*k4)(const int, const int, const int, const int,
16571656 float*, const int, const int,
16581657 const int, const int,
16591658 const float*, const int, const int,
16601659 const int, const int);
16611660 void (*k6)(const int, const int,
16621661 const int, const int,
16631662 const int, const int,
16641663 float*, const int, const int,
16651664 const int, const int,
16661665 const int, const int,
16671666 const float*, const int, const int,
16681667 const int, const int,
16691668 const int, const int);
16701669 switch (fct_nb)
16711670 {
16721671 case IADD:
16731672 k3 = k_ielem_3<IADD>;
16741673 k4 = k_ielem_4<IADD>;
16751674 k6 = k_ielem_6<IADD>;
16761675 break;
16771676 case IDIV:
16781677 k3 = k_ielem_3<IDIV>;
16791678 k4 = k_ielem_4<IDIV>;
16801679 k6 = k_ielem_6<IDIV>;
16811680 break;
16821681 case CPY:
16831682 k3 = k_ielem_3<CPY>;
16841683 k4 = k_ielem_4<CPY>;
16851684 k6 = k_ielem_6<CPY>;
16861685 break;
16871686 default:
16881687 assert (0);
16891688 PyErr_Format(
16901689 PyExc_TypeError,
16911690 "CudaNdarray_inplace_elemwise invalid fct_nb (%i).",
16921691 (int)fct_nb);
16931692 return -1;
16941693 }
16951694 if (!CudaNdarray_Check(py_self)) {
16961695 PyErr_SetString(
16971696 PyExc_TypeError,
16981697 "CudaNdarray_inplace_elemwise need a CudaNdarray on left");
16991698 return -1;
17001699 }
17011700 CudaNdarray * new_other = NULL;
17021701 if (!CudaNdarray_Check(py_other)) {
17031702 new_other = (CudaNdarray*) CudaNdarray_New();
17041703 if(!new_other)
17051704 {
17061705 return -1;
17071706 }
17081707 if(CudaNdarray_CopyFromArray(new_other, (PyArrayObject *) py_other))
17091708 {
17101709 Py_XDECREF(new_other);
17111710 return -1;
17121711 }
17131712 py_other = (PyObject *) new_other;
17141713 }
17151714
17161715 CudaNdarray * self = (CudaNdarray *)py_self;
17171716 CudaNdarray * other = (CudaNdarray *)py_other;
17181717
17191718 if (verbose)
17201719 {
17211720 fprintf(stderr,
17221721 "INPLACE ADD/DIV for self->nd=%d other->nd=%d\n",
17231722 self->nd, other->nd);
17241723 }
17251724
17261725 //standard elemwise nb dim checks
17271726 if (self->nd < other->nd)
17281727 {
17291728 PyErr_Format(
17301729 PyExc_TypeError,
17311730 "CudaNdarray_inplace_elemwise: The destination need more or the"
17321731 " same number of dimensions then the source. Got %d and %d.",
17331732 self->nd, other->nd);
17341733 Py_XDECREF(new_other);
17351734 return -1;
17361735 }
17371736
17381737 //broadcast to the same number of dimensions.
17391738 int* other_dims = (int*) alloca(self->nd * sizeof(int));
17401739 int* other_strides = (int*) alloca(self->nd * sizeof(int));
17411740 int added_dims = self->nd - other->nd;
17421741 // Add the added broadcasted dimensions
17431742 for (int i = 0; i< added_dims; ++i)
17441743 {
17451744 other_dims[i] = 1;
17461745 other_strides[i] = 0;
17471746 }
17481747 // Copy the existing dimensions
17491748 for (int i = 0; i< other->nd; ++i)
17501749 {
17511750 other_dims[i+added_dims] = CudaNdarray_HOST_DIMS(other)[i];
17521751 other_strides[i+added_dims] = CudaNdarray_HOST_STRIDES(other)[i];
17531752 }
17541753
17551754 //standard elemwise dim checks
17561755 unsigned int size = 1;
17571756 for (int i = 0; i< self->nd; ++i)
17581757 {
17591758 if ((CudaNdarray_HOST_DIMS(self)[i] != other_dims[i])
17601759 && (other_dims[i] != 1))
17611760 {
17621761 PyErr_SetString(
17631762 PyExc_ValueError,
17641763 "CudaNdarray_inplace_elemwise need same dimensions (or broadcastable dimension)");
17651764 Py_XDECREF(new_other);
17661765 return -1;
17671766 }
17681767 // if we're broadcasting other, then make sure it has stride 0
17691768 assert ((CudaNdarray_HOST_DIMS(self)[i] == other_dims[i])
17701769 || (other_strides[i] == 0));
17711770 size *= (unsigned int) CudaNdarray_HOST_DIMS(self)[i];
17721771 }
17731772
17741773 if (size==0)
17751774 {
17761775 int other_size = CudaNdarray_SIZE((CudaNdarray *)py_other);
17771776 if (!(other_size == 0 || other_size == 1))
17781777 {
17791778 PyErr_SetString(
17801779 PyExc_ValueError,
17811780 "CudaNdarray_inplace_elemwise cannot work inplace on"
17821781 " un-initialized array when the new value have more than"
17831782 " 0 or 1 broadcastable dimensions");
17841783 Py_XDECREF(new_other);
17851784 return 0;
17861785 }
17871786 Py_XDECREF(new_other);
17881787 return 0;
17891788 }
17901789
17911790 switch(self->nd)
17921791 {
17931792 case 0:
17941793 {
17951794 dim3 n_blocks(1, 1, 1);
17961795 dim3 n_threads(1);
17971796 k3<<<n_blocks, n_threads>>>(
17981797 1, //d0
17991798 1, //d1
18001799 1, //d2
18011800 CudaNdarray_DEV_DATA(self),
18021801 1, //strides
18031802 1,
18041803 1,
18051804 CudaNdarray_DEV_DATA(other),
18061805 1, //strides
18071806 1,
18081807 1);
18091808 CNDA_THREAD_SYNC;
18101809 cudaError_t err = cudaGetLastError();
18111810 if (cudaSuccess != err)
18121811 {
18131812 PyErr_Format(
18141813 PyExc_RuntimeError,
18151814 "CudaNdarray_inplace_elemwise case0: Cuda error: %s: %s.\n",
18161815 "k3",
18171816 cudaGetErrorString(err));
18181817 Py_XDECREF(new_other);
18191818 return -1;
18201819 }
18211820 }
18221821 break;
18231822 case 1:
18241823 {
18251824 dim3 n_blocks(1, 1, 1);
18261825 dim3 n_threads(
18271826 std::min(
18281827 CudaNdarray_HOST_DIMS(self)[0],
18291828 NUM_VECTOR_OP_THREADS_PER_BLOCK));
18301829 k3<<<n_blocks, n_threads>>>(
18311830 1, //dimensions
18321831 1,
18331832 CudaNdarray_HOST_DIMS(self)[0],
18341833 CudaNdarray_DEV_DATA(self),
18351834 1, //strides
18361835 1,
18371836 CudaNdarray_HOST_STRIDES(self)[0],
18381837 CudaNdarray_DEV_DATA(other),
18391838 1, //strides
18401839 1,
18411840 other_strides[0]);
18421841 CNDA_THREAD_SYNC;
18431842 cudaError_t err = cudaGetLastError();
18441843 if (cudaSuccess != err)
18451844 {
18461845 PyErr_Format(
18471846 PyExc_RuntimeError,
18481847 "CudaNdarray_inplace_elemwise case1: Cuda error: %s: %s.\n",
18491848 "k3",
18501849 cudaGetErrorString(err));
18511850 Py_XDECREF(new_other);
18521851 return -1;
18531852 }
18541853 }
18551854 break;
18561855 case 2:
18571856 {
18581857 //TODO: if both self and other are f-contiguous
18591858 // Then flip the block and thread dimensions
18601859 // to make contiguous reads & writes
18611860 dim3 n_blocks(1,
18621861 std::min(
18631862 CudaNdarray_HOST_DIMS(self)[0],
18641863 NUM_VECTOR_OP_BLOCKS));
18651864 dim3 n_threads(
18661865 std::min(
18671866 CudaNdarray_HOST_DIMS(self)[1],
18681867 NUM_VECTOR_OP_THREADS_PER_BLOCK));
18691868 k3<<<n_blocks, n_threads>>>(1,
18701869 CudaNdarray_HOST_DIMS(self)[0],
18711870 CudaNdarray_HOST_DIMS(self)[1],
18721871 CudaNdarray_DEV_DATA(self),
18731872 1,
18741873 CudaNdarray_HOST_STRIDES(self)[0],
18751874 CudaNdarray_HOST_STRIDES(self)[1],
18761875 CudaNdarray_DEV_DATA(other),
18771876 1,
18781877 other_strides[0],
18791878 other_strides[1]);
18801879 CNDA_THREAD_SYNC;
18811880 cudaError_t err = cudaGetLastError();
18821881 if (cudaSuccess != err)
18831882 {
18841883 PyErr_Format(
18851884 PyExc_RuntimeError,
18861885 "CudaNdarray_inplace_elemwise case2: Cuda error: %s: %s.\n",
18871886 "k3",
18881887 cudaGetErrorString(err));
18891888 Py_XDECREF(new_other);
18901889 return -1;
18911890 }
18921891 }
18931892 break;
18941893 case 3:
18951894 {
18961895 //TODO: Dimshuffle so that at least one of the arrays
18971896 // has a contiguous dimension on the thread idx.
18981897 dim3 n_blocks(
18991898 std::min(
19001899 CudaNdarray_HOST_DIMS(self)[0],
19011900 NUM_VECTOR_OP_BLOCKS),
19021901 CudaNdarray_HOST_DIMS(self)[1]);
19031902 while (n_blocks.x * n_blocks.y > NUM_VECTOR_OP_BLOCKS)
19041903 n_blocks.y /= 2;
19051904 dim3 n_threads(
19061905 std::min(
19071906 CudaNdarray_HOST_DIMS(self)[2],
19081907 NUM_VECTOR_OP_THREADS_PER_BLOCK));
19091908 k3<<<n_blocks, n_threads>>>(
19101909 CudaNdarray_HOST_DIMS(self)[0],
19111910 CudaNdarray_HOST_DIMS(self)[1],
19121911 CudaNdarray_HOST_DIMS(self)[2],
19131912 CudaNdarray_DEV_DATA(self),
19141913 CudaNdarray_HOST_STRIDES(self)[0],
19151914 CudaNdarray_HOST_STRIDES(self)[1],
19161915 CudaNdarray_HOST_STRIDES(self)[2],
19171916 CudaNdarray_DEV_DATA(other),
19181917 other_strides[0],
19191918 other_strides[1],
19201919 other_strides[2]);
19211920 CNDA_THREAD_SYNC;
19221921 cudaError_t err = cudaGetLastError();
19231922 if (cudaSuccess != err)
19241923 {
19251924 PyErr_Format(
19261925 PyExc_RuntimeError,
19271926 "CudaNdarray_inplace_elemwise case3: Cuda error: %s: %s.\n",
19281927 "k3",
19291928 cudaGetErrorString(err));
19301929 Py_XDECREF(new_other);
19311930 return -1;
19321931 }
19331932 }
19341933 break;
19351934 case 4:
19361935 {
19371936 dim3 n_blocks(
19381937 std::min(
19391938 CudaNdarray_HOST_DIMS(self)[0],
19401939 NUM_VECTOR_OP_BLOCKS),
19411940 CudaNdarray_HOST_DIMS(self)[1]
19421941 );
19431942 while (n_blocks.x * n_blocks.y > NUM_VECTOR_OP_BLOCKS)
19441943 n_blocks.y /= 2;
19451944 dim3 n_threads(
19461945 std::min(
19471946 CudaNdarray_HOST_DIMS(self)[2],
19481947 NUM_VECTOR_OP_THREADS_PER_BLOCK)
19491948 //TODO: DON"T YOU NEED OT PUT DIMS[3] in here???
19501949 );
19511950 k4<<<n_blocks, n_threads>>>(
19521951 CudaNdarray_HOST_DIMS(self)[0],
19531952 CudaNdarray_HOST_DIMS(self)[1],
19541953 CudaNdarray_HOST_DIMS(self)[2],
19551954 CudaNdarray_HOST_DIMS(self)[3],
19561955 CudaNdarray_DEV_DATA(self),
19571956 CudaNdarray_HOST_STRIDES(self)[0],
19581957 CudaNdarray_HOST_STRIDES(self)[1],
19591958 CudaNdarray_HOST_STRIDES(self)[2],
19601959 CudaNdarray_HOST_STRIDES(self)[3],
19611960 CudaNdarray_DEV_DATA(other),
19621961 other_strides[0],
19631962 other_strides[1],
19641963 other_strides[2],
19651964 other_strides[3]);
19661965 CNDA_THREAD_SYNC;
19671966 cudaError_t err = cudaGetLastError();
19681967 if (cudaSuccess != err)
19691968 {
19701969 PyErr_Format(
19711970 PyExc_RuntimeError,
19721971 "CudaNdarray_inplace_elemwise case4: Cuda error: %s: %s.\n",
19731972 "k4",
19741973 cudaGetErrorString(err));
19751974 Py_XDECREF(new_other);
19761975 return -1;
19771976 }
19781977 }
19791978 break;
19801979 case 5:
19811980 {
19821981 dim3 n_blocks(
19831982 std::min(
19841983 CudaNdarray_HOST_DIMS(self)[1],
19851984 NUM_VECTOR_OP_BLOCKS),
19861985 CudaNdarray_HOST_DIMS(self)[2]);
19871986 while (n_blocks.x * n_blocks.y > NUM_VECTOR_OP_BLOCKS)
19881987 n_blocks.y /= 2;
19891988 dim3 n_threads(
19901989 std::min(
19911990 CudaNdarray_HOST_DIMS(self)[3],
19921991 NUM_VECTOR_OP_THREADS_PER_BLOCK)
19931992 //TODO: DON"T YOU NEED OT PUT DIMS[3] in here???
19941993 );
19951994 for (int i = 0; i < CudaNdarray_HOST_DIMS(self)[0]; ++i)
19961995 {
19971996 k4<<<n_blocks, n_threads>>>(
19981997 CudaNdarray_HOST_DIMS(self)[1],
19991998 CudaNdarray_HOST_DIMS(self)[2],
20001999 CudaNdarray_HOST_DIMS(self)[3],
20012000 CudaNdarray_HOST_DIMS(self)[4],
20022001 CudaNdarray_DEV_DATA(self) + i * CudaNdarray_HOST_STRIDES(self)[0],
20032002 CudaNdarray_HOST_STRIDES(self)[1],
20042003 CudaNdarray_HOST_STRIDES(self)[2],
20052004 CudaNdarray_HOST_STRIDES(self)[3],
20062005 CudaNdarray_HOST_STRIDES(self)[4],
20072006 CudaNdarray_DEV_DATA(other) + i * other_strides[0],
20082007 other_strides[1],
20092008 other_strides[2],
20102009 other_strides[3],
20112010 other_strides[4]);
20122011 CNDA_THREAD_SYNC;
20132012 cudaError_t err = cudaGetLastError();
20142013 if( cudaSuccess != err)
20152014 {
20162015 PyErr_Format(
20172016 PyExc_RuntimeError,
20182017 "CudaNdarray_inplace_elemwise case5: Cuda error: %s: %s. n_block=(%ld,%ld) n_threads=%ld\n",
20192018 "k5 with loop over k4",
20202019 cudaGetErrorString(err),
20212020 (long) n_blocks.x, (long) n_blocks.y, (long) n_threads.x);
20222021 Py_XDECREF(new_other);
20232022 return -1;
20242023 }
20252024 }
20262025 }
20272026 break;
20282027 case 6:
20292028 {
20302029 dim3 n_blocks(
20312030 std::min(
20322031 CudaNdarray_HOST_DIMS(self)[0],
20332032 NUM_VECTOR_OP_BLOCKS),
20342033 CudaNdarray_HOST_DIMS(self)[1],
20352034 CudaNdarray_HOST_DIMS(self)[2]
20362035 );
20372036 while (n_blocks.x * n_blocks.y > NUM_VECTOR_OP_BLOCKS)
20382037 n_blocks.y /= 2;
20392038 // GTX285(compute capabilities 1.3) don't support n_blocks.z > 1
20402039 // (compute capabilities 2.0) support 65535 for n_blocks.z
20412040 //while (n_blocks.x * n_blocks.y * n_blocks.z > NUM_VECTOR_OP_BLOCKS)
20422041 // n_blocks.z /= 2;
20432042 n_blocks.z = 1;
20442043 dim3 n_threads(
20452044 std::min(
20462045 CudaNdarray_HOST_DIMS(self)[3],
20472046 NUM_VECTOR_OP_THREADS_PER_BLOCK)
20482047 //TODO: DON'T YOU NEED TO PUT DIMS[4] in here???
20492048 //TODO: DON'T YOU NEED TO PUT DIMS[5] in here???
20502049 );
20512050 k6<<<n_blocks, n_threads>>>(
20522051 CudaNdarray_HOST_DIMS(self)[0],
20532052 CudaNdarray_HOST_DIMS(self)[1],
20542053 CudaNdarray_HOST_DIMS(self)[2],
20552054 CudaNdarray_HOST_DIMS(self)[3],
20562055 CudaNdarray_HOST_DIMS(self)[4],
20572056 CudaNdarray_HOST_DIMS(self)[5],
20582057 CudaNdarray_DEV_DATA(self),
20592058 CudaNdarray_HOST_STRIDES(self)[0],
20602059 CudaNdarray_HOST_STRIDES(self)[1],
20612060 CudaNdarray_HOST_STRIDES(self)[2],
20622061 CudaNdarray_HOST_STRIDES(self)[3],
20632062 CudaNdarray_HOST_STRIDES(self)[4],
20642063 CudaNdarray_HOST_STRIDES(self)[5],
20652064 CudaNdarray_DEV_DATA(other),
20662065 other_strides[0],
20672066 other_strides[1],
20682067 other_strides[2],
20692068 other_strides[3],
20702069 other_strides[4],
20712070 other_strides[5]);
20722071 CNDA_THREAD_SYNC;
20732072 cudaError_t err = cudaGetLastError();
20742073 if (cudaSuccess != err)
20752074 {
20762075 PyErr_Format(
20772076 PyExc_RuntimeError,
20782077 "CudaNdarray_inplace_elemwise case6: Cuda error: %s: %s. n_blocks=(%ld, %ld, %ld) n_threads=(%ld)\n",
20792078 "k6",
20802079 cudaGetErrorString(err),
20812080 (long) n_blocks.x, (long) n_blocks.y, (long) n_blocks.z,
20822081 (long) n_threads.x);
20832082 Py_XDECREF(new_other);
20842083 return -1;
20852084 }
20862085 }
20872086 break;
20882087 default:
20892088 {
20902089 PyErr_Format(
20912090 PyExc_NotImplementedError,
20922091 "inplace_elemwise w nd=%i\n",
20932092 self->nd);
20942093 Py_XDECREF(new_other);
20952094 return -1;
20962095 }
20972096 }
20982097 if (verbose)
20992098 fprintf(stderr, "INPLACE ADD/DIV end\n");
21002099 Py_XDECREF(new_other);
21012100 return 0;
21022101 }
21032102
21042103 /*
21052104 * We need this inplace Add to support IncSubTensor
21062105 * It returns py_self on success with an additional reference. Else NULL.
21072106 */
21082107 // Will be called by __iadd__ in Python
21092108 PyObject *
21102109 CudaNdarray_inplace_add(PyObject* py_self, PyObject * py_other)
21112110 {
21122111 if (CudaNdarray_inplace_elemwise(py_self, py_other, IADD))
21132112 {
21142113 return NULL;
21152114 }
21162115 Py_INCREF(py_self);
21172116 return py_self;
21182117 }
21192118
21202119 /*
21212120 * We need this inplace div for cuda/tests/test_basic_ops.py:test_shared_options
21222121 * It returns py_self on success with an additional reference. Else NULL.
21232122 */
21242123 // Will be called by __idiv__ in Python
21252124 static PyObject *
21262125 CudaNdarray_inplace_div(PyObject* py_self, PyObject * py_other)
21272126 {
21282127 if (CudaNdarray_inplace_elemwise(py_self, py_other, IDIV))
21292128 {
21302129 return NULL;
21312130 }
21322131 Py_INCREF(py_self);
21332132 return py_self;
21342133 }
21352134
21362135 // The PyNumberMethods struct layout changed in a non-trivial way from 2 to 3.
21372136 #if PY_MAJOR_VERSION == 3
21382137 static PyNumberMethods CudaNdarrayNumberMethods =
21392138 {
21402139 (binaryfunc)CudaNdarray_add, //binaryfunc nb_add; __add__
21412140 0, //binaryfunc nb_subtract;
21422141 0, //binaryfunc nb_multiply;
21432142 0, //binaryfunc nb_remainder;
21442143 0, //binaryfunc nb_divmod;
21452144 0, //ternaryfunc nb_power;
21462145 0, //unaryfunc nb_negative;
21472146 0, //unaryfunc nb_positive;
21482147 0, //unaryfunc nb_absolute;
21492148 0, //inquiry nb_bool;
21502149 0, //unaryfunc nb_invert;
21512150 0, //binaryfunc nb_lshift;
21522151 0, //binaryfunc nb_rshift;
21532152 0, //binaryfunc nb_and;
21542153 0, //binaryfunc nb_xor;
21552154 0, //binaryfunc nb_or;
21562155 0, //unaryfunc nb_int;
21572156 0, //void *nb_reserved;
21582157 0, //unaryfunc nb_float;
21592158
21602159 (binaryfunc)CudaNdarray_inplace_add, //binaryfunc nb_inplace_add; __iadd__
21612160 0, //binaryfunc nb_inplace_subtract;
21622161 0, //binaryfunc nb_inplace_multiply;
21632162 0, //binaryfunc nb_inplace_remainder;
21642163 0, //ternaryfunc nb_inplace_power;
21652164 0, //binaryfunc nb_inplace_lshift;
21662165 0, //binaryfunc nb_inplace_rshift;
21672166 0, //binaryfunc nb_inplace_and;
21682167 0, //binaryfunc nb_inplace_xor;
21692168 0, //binaryfunc nb_inplace_or;
21702169
21712170 0, //binaryfunc nb_floor_divide;
21722171 0, //binaryfunc nb_true_divide;
21732172 0, //binaryfunc nb_inplace_floor_divide;
21742173 (binaryfunc)CudaNdarray_inplace_div, //binaryfunc nb_inplace_true_divide; __idiv__
21752174
21762175 0, //unaryfunc nb_index
21772176 };
21782177 #else
21792178 static PyNumberMethods CudaNdarrayNumberMethods =
21802179 {
21812180 (binaryfunc)CudaNdarray_add, //binaryfunc nb_add; __add__
21822181 0, //binaryfunc nb_subtract; __sub__
21832182 0, //binaryfunc nb_multiply; __mul__
21842183 0, //binaryfunc nb_divide; __div__
21852184 0, //binaryfunc nb_remainder; __mod__
21862185 0, //binaryfunc nb_divmod; __divmod__
21872186 0, //ternaryfunc nb_power; __pow__
21882187 0, //unaryfunc nb_negative; __neg__
21892188 0, //unaryfunc nb_positive; __pos__
21902189 0, //unaryfunc nb_absolute; __abs__
21912190 0, //inquiry nb_nonzero; __nonzero__ /* Used by PyObject_IsTrue */
21922191 0, //unaryfunc nb_invert; __invert__
21932192 0, //binaryfunc nb_lshift; __lshift__
21942193 0, //binaryfunc nb_rshift; __rshift__
21952194 0, //binaryfunc nb_and; __and__
21962195 0, //binaryfunc nb_xor; __xor__
21972196 0, //binaryfunc nb_or; __or__
21982197 0, //coercion nb_coerce; __coerce__ /* Used by the coerce() function */
21992198 0, //unaryfunc nb_int; __int__
22002199 0, //unaryfunc nb_long; __long__
22012200 0, //unaryfunc nb_float; __float__
22022201 0, //unaryfunc nb_oct; __oct__
22032202 0, //unaryfunc nb_hex; __hex__
22042203
22052204 /* Added in release 2.0 */
22062205 (binaryfunc)CudaNdarray_inplace_add, //binaryfunc nb_inplace_add; __iadd__
22072206 0, //binaryfunc nb_inplace_subtract; __isub__
22082207 0, //binaryfunc nb_inplace_multiply; __imul__
22092208 (binaryfunc)CudaNdarray_inplace_div, //binaryfunc nb_inplace_divide; __idiv__
22102209 0, //binaryfunc nb_inplace_remainder; __imod__
22112210 0, //ternaryfunc nb_inplace_power; __ipow__
22122211 0, //binaryfunc nb_inplace_lshift; __ilshift__
22132212 0, //binaryfunc nb_inplace_rshift; __irshift__
22142213 0, //binaryfunc nb_inplace_and; __iand__
22152214 0, //binaryfunc nb_inplace_xor; __ixor__
22162215 0, //binaryfunc nb_inplace_or; __ior__
22172216
22182217 /* Added in release 2.2 */
22192218 0, //binaryfunc nb_floor_divide; __floordiv__
22202219 0, //binaryfunc nb_true_divide; __truediv__
22212220 0, //binaryfunc nb_inplace_floor_divide; __ifloordiv__
22222221 (binaryfunc)CudaNdarray_inplace_div, //binaryfunc nb_inplace_true_divide; __itruediv__
22232222
22242223 #if PY_MINOR_VERSION > 4
22252224 /* Added in release 2.5 */
22262225 0 //unaryfunc nb_index; __index__
22272226 #endif
22282227 };
22292228 #endif
22302229
22312230
22322231 /////////////////////
22332232 // Mapping protocol
22342233 /////////////////////
22352234
22362235 // Will by called by __len__ in Python
22372236 static Py_ssize_t
22382237 CudaNdarray_len(PyObject * py_self)
22392238 {
22402239 CudaNdarray * self = (CudaNdarray*) py_self;
22412240 if (self->nd <= 0)
22422241 {
22432242 return (Py_ssize_t) 0;
22442243 }
22452244 else
22462245 {
22472246 return (Py_ssize_t) CudaNdarray_HOST_DIMS(self)[0];
22482247 }
22492248 }
22502249
22512250 // Will by called by __getitem__ in Python
22522251 PyObject *
22532252 CudaNdarray_Subscript(PyObject * py_self, PyObject * key)
22542253 {
22552254 int verbose = 0;
22562255 if (verbose) fprintf(stderr, "Subscript .... \n");
22572256 CudaNdarray * self = (CudaNdarray*) py_self;
22582257 PyObject * py_rval = NULL;
22592258 CudaNdarray * rval = NULL;
22602259 PyObject * intobj = NULL;
22612260
22622261 //PyObject_Print(key, stderr, 0);
22632262
22642263 if (key == Py_Ellipsis)
22652264 {
22662265 Py_INCREF(py_self);
22672266 return py_self;
22682267 }
22692268 if ((intobj=PyNumber_Int(key))) //INDEXING BY INTEGER
22702269 //else if (PyInt_Check(key)) //INDEXING BY INTEGER
22712270 {
22722271 int d_idx = PyInt_AsLong(intobj);
22732272 Py_DECREF(intobj); intobj=NULL;
22742273 //int d_idx = PyInt_AsLong(key);
22752274 if (self->nd == 0)
22762275 {
22772276 PyErr_SetString(PyExc_IndexError, "0-d arrays can't be indexed");
22782277 return NULL;
22792278 }
22802279 int d_dim = CudaNdarray_HOST_DIMS(self)[0];
22812280 int offset = 0;
22822281
22832282 if ((d_idx >= 0) && (d_idx < d_dim))
22842283 {
22852284 //normal indexing
22862285 offset += d_idx * CudaNdarray_HOST_STRIDES(self)[0];
22872286 }
22882287 else if ((d_idx < 0) && (d_idx >= -d_dim))
22892288 {
22902289 //end-based indexing
22912290 // d_idx is negative
22922291 offset += (d_dim + d_idx) * CudaNdarray_HOST_STRIDES(self)[0];
22932292 }
22942293 else
22952294 {
22962295 PyErr_Format(PyExc_IndexError,
22972296 "index out of bounds. Asked %d, but size of %d",
22982297 d_idx, d_dim);
22992298 return NULL;
23002299 }
23012300
23022301 //allocate our subtensor view
23032302 py_rval = CudaNdarray_new_nd(self->nd - 1);
23042303 rval = (CudaNdarray*) py_rval;
23052304 if (!rval) return NULL;
23062305 assert (0 == rval->data_allocated);
23072306
23082307 //initialize the view's data pointer to our own.
23092308 if (CudaNdarray_set_device_data(rval, CudaNdarray_DEV_DATA(self) + offset, self))
23102309 {
23112310 Py_DECREF(rval);
23122311 return NULL;
23132312 }
23142313 for (int d = 1; d < self->nd; ++d)
23152314 {
23162315 CudaNdarray_set_stride(rval, d-1, CudaNdarray_HOST_STRIDES(self)[d]);
23172316 CudaNdarray_set_dim(rval, d-1, CudaNdarray_HOST_DIMS(self)[d]);
23182317 }
23192318 }
23202319 else
23212320 {
23222321 PyErr_Clear();
23232322 }
23242323 if (PySlice_Check(key)) //INDEXING BY SLICE
23252324 {
23262325 if (verbose) fprintf(stderr, "by slice\n");
23272326 if (self->nd == 0)
23282327 {
23292328 PyErr_SetString(PyExc_ValueError, "cannot slice a 0-d array");
23302329 return NULL;
23312330 }
23322331
23332332 int d_dim = CudaNdarray_HOST_DIMS(self)[0];
23342333 Py_ssize_t start, stop, step, slen;
23352334 if (PySlice_GetIndicesEx(SLICE_CAST(key), d_dim, &start, &stop, &step, &slen))
23362335 {
23372336 if (verbose)
23382337 fprintf(stderr, "PySlice_GetIndicesEx failed\n");
23392338 return NULL;
23402339 }
23412340 if (verbose)
23422341 {
23432342 std::cerr << "start " << start << "\n";
23442343 std::cerr << "stop " << stop << "\n";
23452344 std::cerr << "step " << step << "\n";
23462345 std::cerr << "slen " << slen << "\n";
23472346 }
23482347
23492348 //allocate our subtensor view
23502349 py_rval = CudaNdarray_new_nd(self->nd);
23512350 rval = (CudaNdarray*) py_rval;
23522351 if (!rval) return NULL;
23532352 assert (0 == rval->data_allocated);
23542353
23552354
23562355 //initialize the view's data pointer to our own.
23572356 if (CudaNdarray_set_device_data(rval,
23582357 CudaNdarray_DEV_DATA(self) + start * CudaNdarray_HOST_STRIDES(self)[0],
23592358 self))
23602359 {
23612360 Py_DECREF(rval);
23622361 return NULL;
23632362 }
23642363 //initialize dimension 0 of rval
23652364 CudaNdarray_set_stride(rval, 0,
23662365 (slen == 1) ? 0 : step * CudaNdarray_HOST_STRIDES(self)[0]);
23672366 CudaNdarray_set_dim(rval, 0, slen);
23682367 if (verbose) std::cerr << "rval stride " << CudaNdarray_HOST_STRIDES(rval)[0] << "\n";
23692368 // initialize dimensions > 0 of rval
23702369 for (int d = 1; d < self->nd; ++d)
23712370 {
23722371 CudaNdarray_set_stride(rval, d, CudaNdarray_HOST_STRIDES(self)[d]);
23732372 CudaNdarray_set_dim(rval, d, CudaNdarray_HOST_DIMS(self)[d]);
23742373 }
23752374 }
23762375 if (PyTuple_Check(key)) //INDEXING BY TUPLE
23772376 {
23782377 if (verbose) fprintf(stderr, "by tuple\n");
23792378 //elements of the tuple can be either integers or slices
23802379 //the dimensionality of the view we will return is diminished for each slice in the tuple
23812380
23822381 if (PyTuple_Size(key) > self->nd)
23832382 {
23842383 PyErr_SetString(PyExc_IndexError, "index error");
23852384 return NULL;
23862385 }
23872386
23882387 //calculate the number of dimensions in the return value
23892388 int rval_nd = self->nd;
23902389 for (int d = 0; d < PyTuple_Size(key); ++d)
23912390 {
23922391 //On some paltform PyInt_Check(<type 'numpy.int64'>) return true, other it return false.
23932392 //So we use PyArray_IsAnyScalar that should covert everything.
23942393 rval_nd -= PyArray_IsAnyScalar(PyTuple_GetItem(key, d));
23952394 }
23962395
23972396 //allocate our subtensor view
23982397 py_rval = CudaNdarray_new_nd(rval_nd);
23992398 rval = (CudaNdarray*) py_rval;
24002399 if (!rval) return NULL;
24012400 assert (0 == rval->data_allocated);
24022401
24032402 //initialize the view's data pointer to our own.
24042403 if (CudaNdarray_set_device_data(rval, CudaNdarray_DEV_DATA(self), self))
24052404 {
24062405 Py_DECREF(rval);
24072406 return NULL;
24082407 }
24092408
24102409 // rval_d will refer to the current dimension in the rval.
24112410 // It will not be incremented for integer keys, but will be incremented for slice
24122411 // keys
24132412 int rval_d = 0;
24142413
24152414 for (int d = 0; d < self->nd; ++d)
24162415 {
24172416 // keys can be shorter than self->nd.
24182417 // when that happens, it means that the remaining dimensions are "full slices"
24192418 if (d >=PyTuple_Size(key))
24202419 {
24212420 CudaNdarray_set_stride(rval, rval_d, CudaNdarray_HOST_STRIDES(self)[d]);
24222421 CudaNdarray_set_dim(rval, rval_d, CudaNdarray_HOST_DIMS(self)[d]);
24232422 ++rval_d;
24242423 }
24252424 else
24262425 {
24272426 PyObject * key_d = PyTuple_GetItem(key, d);
24282427
24292428 if (PySlice_Check(key_d))
24302429 {
24312430 Py_ssize_t start, stop, step, slen;
24322431 if (PySlice_GetIndicesEx(SLICE_CAST(key_d), CudaNdarray_HOST_DIMS(self)[d], &start, &stop, &step, &slen))
24332432 {
24342433 Py_DECREF(rval);
24352434 return NULL;
24362435 }
24372436 rval->devdata += start * CudaNdarray_HOST_STRIDES(self)[d];
24382437 CudaNdarray_set_stride(rval, rval_d,
24392438 (slen == 1) ? 0 : step * CudaNdarray_HOST_STRIDES(self)[d]);
24402439 CudaNdarray_set_dim(rval, rval_d, slen);
24412440 if (0)
24422441 {
24432442 std::cerr << "start " << start << "\n";
24442443 std::cerr << "stop " << stop << "\n";
24452444 std::cerr << "step " << step << "\n";
24462445 std::cerr << "slen " << slen << "\n";
24472446 }
24482447 ++rval_d;
24492448 }
24502449 else if ((intobj=PyNumber_Int(key_d)))
24512450 {
24522451 assert(PyArray_IsAnyScalar(key_d));
24532452 int d_idx = PyInt_AsLong(intobj);
24542453 Py_DECREF(intobj);
24552454 intobj = NULL;
24562455 int d_dim = CudaNdarray_HOST_DIMS(self)[d];
24572456
24582457 if ((d_idx >= 0) && (d_idx < d_dim))
24592458 {
24602459 //normal indexing
24612460 rval->devdata += d_idx * CudaNdarray_HOST_STRIDES(self)[d];
24622461 }
24632462 else if ((d_idx < 0) && (d_idx >= -d_dim))
24642463 {
24652464 //end-based indexing
24662465 rval->devdata += (d_dim + d_idx) * CudaNdarray_HOST_STRIDES(self)[d];
24672466 }
24682467 else
24692468 {
24702469 PyErr_Format(PyExc_IndexError,
24712470 "index out of bounds. Asked %d for dimensions %d, but size of %d",
24722471 d_idx, d, d_dim);
24732472 Py_DECREF(rval);
24742473 return NULL;
24752474 }
24762475 }
24772476 else
24782477 {
24792478 PyErr_Clear(); // clear the error set by PyNumber_Int
24802479 PyErr_SetString(PyExc_IndexError, "index must be either int or slice");
24812480 Py_DECREF(rval);
24822481 return NULL;
24832482 }
24842483 }
24852484 }
24862485 }
24872486 if (py_rval)
24882487 {
24892488 if (verbose) fprint_CudaNdarray(stderr, self);
24902489 if (verbose) fprint_CudaNdarray(stderr, rval);
24912490 }
24922491 else
24932492 {
24942493 PyErr_SetString(PyExc_NotImplementedError, "Unknown key type");
24952494 return NULL;
24962495 }
24972496 return py_rval;
24982497 }
24992498
25002499 // Will by called by __setitem__ in Python
25012500 // See http://docs.python.org/dev/py3k/c-api/object.html#PyObject_SetItem
25022501 // Doesn't handle broadcasting, e.g. a[:] = 5
25032502 // Can only be assigned from a CudaNdarray on the right side
25042503 // Or a ndarray
25052504 // Or a python scalar with value 0 when the left side part is c contiguous.
25062505 static int
25072506 CudaNdarray_setitem(PyObject *o, PyObject *key, PyObject *value)
25082507 {
25092508 int verbose = 0;
25102509 if (verbose) fprintf(stderr, "CudaNdarray_setitem start\n");
25112510 // We try to copy directly into this CudaNdarray from the ndarray
25122511 CudaNdarray* rval = (CudaNdarray*)CudaNdarray_Subscript(o, key);
25132512 CudaNdarray* new_value = NULL;
25142513
25152514 if(!rval){
25162515 // CudaNdarray_Subscript failed and set the error msg.
25172516 Py_XDECREF(rval);
25182517 return -1;
25192518 }
25202519
25212520 if(rval != (CudaNdarray*)o &&
25222521 (rval->data_allocated ||
25232522 // The new array should have a base
25242523 !(((CudaNdarray*)rval)->base) ||
25252524 // If the original array has no base, the base of the new
25262525 // array should be the original one
25272526 (!((CudaNdarray*)o)->base && ((CudaNdarray*)rval)->base != o) ||
25282527 // Else, the two arrays should have the same base
25292528 (((CudaNdarray*)o)->base && ((CudaNdarray*)rval)->base != ((CudaNdarray*)o)->base)))
25302529 {
25312530 // This case shouldn't happen, based on what I see in Subscript
25322531 // but just in case it happens sometime in the future
25332532
25342533 PyErr_Format(PyExc_RuntimeError,
25352534 "__getitem__ must return a CudaNdarray that refers to"
25362535 " the original CudaNdarray, not a copy. rval.base=%p"
25372536 " o.base=%p o=%p",
25382537 (((CudaNdarray*)rval)->base), ((CudaNdarray*)o)->base, o);
25392538 Py_DECREF(rval);
25402539 return -1;
25412540 }
25422541
25432542 PyObject * intobj = NULL;
25442543 if (CudaNdarray_Check(o) && PyArray_Check(value)){
25452544 if (verbose)
25462545 fprintf(stderr,
25472546 "CudaNdarray_setitem dest is a CudaNdarray and"
25482547 " value is a ndarray\n");
25492548 new_value = (CudaNdarray*) CudaNdarray_New();
25502549 if(!new_value)
25512550 {
25522551 return -1;
25532552 }
25542553 if (CudaNdarray_CopyFromArray(new_value, (PyArrayObject *) value))
25552554 {
25562555 Py_XDECREF(new_value);
25572556 Py_XDECREF(rval);
25582557 return -1;
25592558 }
25602559 value = (PyObject *) new_value;
25612560 }
25622561 else if ((intobj=PyNumber_Int(value)))
25632562 {
25642563 if (verbose)
25652564 fprintf(stderr,
25662565 "CudaNdarray_setitem dest and value is a python number\n");
25672566 if(! CudaNdarray_is_c_contiguous(rval)){
25682567 PyErr_SetString(PyExc_NotImplementedError,
25692568 "CudaNdarray.__setitem__: When the new value is a scalar"
25702569 " of value 0 the part where we copy to must be c contiguous.");
25712570 Py_XDECREF(rval);
25722571 return -1;
25732572 }
25742573
25752574 long val = PyInt_AsLong(intobj);
25762575 Py_DECREF(intobj); intobj=NULL;
25772576 if (val == 0)
25782577 {
25792578 cudaError_t err = cudaMemset(rval->devdata, 0,
25802579 CudaNdarray_SIZE(rval) * sizeof(real));
25812580 Py_XDECREF(rval);
25822581 if (err)
25832582 {
25842583 // Clear the error flag, cudaMemset doesn't do it.
25852584 // Currently this returns the same thing as err, but if in future
25862585 // it returns something else I still don't see why we should ignore
25872586 // it. All we want to do here is reset the flag.
25882587 cudaGetLastError();
25892588 PyErr_SetString(PyExc_RuntimeError,
25902589 "CudaNdarray.__setitem__: cudaMemset failed");
25912590 return -1;
25922591 }
25932592 return 0;
25942593 } else {
25952594 Py_XDECREF(rval);
25962595 PyErr_SetString(PyExc_NotImplementedError,
25972596 "CudaNdarray.__setitem__: we support setting only python"
25982597 " scalar of value 0, numpy nd array and CudaNdarray.");
25992598 return -1;
26002599 }
26012600 }
26022601
26032602 PyErr_Clear(); // clear PyNumber_Int error.
26042603
26052604 if(!CudaNdarray_Check(o) || !CudaNdarray_Check(value))
26062605 {
26072606 PyErr_SetString(PyExc_TypeError,
26082607 "CudaNdarray.__setitem__: left must be a CudaNdarrays and right"
26092608 " must be a CudaNdarrays, an ndarray or a python scalar of value 0.");
26102609 Py_XDECREF(new_value);
26112610 return -1;
26122611 }
26132612
26142613 if (verbose)
26152614 fprintf(stderr, "CudaNdarray_setitem dest and value are CudaNdarray\n");
26162615
26172616 if (cnda_copy_structure_to_device(rval))
26182617 {
26192618 PyErr_SetString(PyExc_RuntimeError,
26202619 "CudaNdarray.__setitem__: syncing structure to device failed");
26212620 Py_DECREF(rval);
26222621 Py_XDECREF(new_value);
26232622
26242623 if (verbose)
26252624 fprintf(stderr, "CudaNdarray_setitem error end\n");
26262625 return -1;
26272626 }
26282627
26292628 PyObject *baseSavedForComparison = rval->base;
26302629
26312630 if (CudaNdarray_CopyFromCudaNdarray(rval, (CudaNdarray*)value, true))
26322631 {
26332632 Py_DECREF((PyObject*)rval);
26342633 Py_XDECREF(new_value);
26352634
26362635 if (verbose)
26372636 fprintf(stderr, "CudaNdarray_setitem error end\n");
26382637 return -1;
26392638 }
26402639
26412640 assert (rval->base == baseSavedForComparison);
26422641 assert (rval->dev_structure_fresh);
26432642
26442643 // Clean up locally-created references
26452644 Py_DECREF(rval);
26462645 Py_XDECREF(new_value);
26472646
26482647 return 0;
26492648 }
26502649
26512650
26522651 PyMappingMethods CudaNdarrayMappingMethods = {
26532652 CudaNdarray_len, //lenfunc mp_length; __len__
26542653 CudaNdarray_Subscript, //binaryfunc mp_subscript; __getitem__
26552654 CudaNdarray_setitem //objobjargproc mp_ass_subscript; __setitem__
26562655 };
26572656
26582657 ////////////////////
26592658 //
26602659 ////////////////////
26612660
26622661 static PyObject *
26632662 CudaNdarray_get_shape(CudaNdarray *self, void *closure)
26642663 {
26652664 if (self->nd < 0)
26662665 {
26672666 PyErr_SetString(PyExc_ValueError, "CudaNdarray not initialized");
26682667 return NULL;
26692668 }
26702669 PyObject * rval = PyTuple_New(self->nd);
26712670 for (int i = 0; i < self->nd; ++i)
26722671 {
26732672 if (!rval || PyTuple_SetItem(rval, i, PyInt_FromLong(CudaNdarray_HOST_DIMS(self)[i])))
26742673 {
26752674 Py_XDECREF(rval);
26762675 return NULL;
26772676 }
26782677
26792678 }
26802679 return rval;
26812680 }
26822681
26832682 static int
26842683 CudaNdarray_set_shape(CudaNdarray *self, PyObject *value, void *closure)
26852684 {
26862685 PyErr_SetString(PyExc_NotImplementedError, "TODO: call reshape");
26872686 return -1;
26882687 }
26892688
26902689 static PyObject *
26912690 CudaNdarray_get_strides(CudaNdarray *self, void *closure)
26922691 {
26932692 if (self->nd < 0)
26942693 {
26952694 PyErr_SetString(PyExc_ValueError, "CudaNdarray not initialized");
26962695 return NULL;
26972696 }
26982697 PyObject * rval = PyTuple_New(self->nd);
26992698 for (int i = 0; i < self->nd; ++i)
27002699 {
27012700 if (!rval || PyTuple_SetItem(rval, i, PyInt_FromLong(CudaNdarray_HOST_STRIDES(self)[i])))
27022701 {
27032702 Py_XDECREF(rval);
27042703 return NULL;
27052704 }
27062705
27072706 }
27082707 return rval;
27092708 }
27102709
27112710 static int
27122711 CudaNdarray_set_strides(CudaNdarray *self, PyObject *value, void *closure)
27132712 {
27142713 //npy_intp newstrides_bytes[PyTuple_Size(value)];
27152714 if (PyTuple_Check(value)){
27162715 if (PyTuple_Size(value) != CudaNdarray_NDIM(self)){
27172716 PyErr_SetString(PyExc_ValueError,
27182717 "The new strides tuple must have the same length"
27192718 " as the number of dimensions");
27202719 return -1;
27212720 }
27222721 }else if (PyList_Check(value)){
27232722 if (PyList_Size(value) != CudaNdarray_NDIM(self)){
27242723 PyErr_SetString(PyExc_ValueError,
27252724 "The new strides list must have the same length"
27262725 " as the number of dimensions");
27272726 return -1;
27282727 }
27292728 }else{
27302729 PyErr_SetString(PyExc_ValueError,
27312730 "The new strides need to be encoded in a tuple or list");
27322731 return -1;
27332732 }
27342733 npy_intp* newstrides = (npy_intp*) alloca(CudaNdarray_NDIM(self) * sizeof(npy_intp));
27352734 if (PyTuple_Check(value)){
27362735 for(int i=0; i < CudaNdarray_NDIM(self); i++){
27372736 newstrides[i] = PyInt_AsLong(PyTuple_GetItem(value, Py_ssize_t(i)));
27382737 //newstrides_bytes[i] = newstrides[i] * 4;
27392738 }
27402739 }else if (PyList_Check(value)){
27412740 for(int i=0; i < CudaNdarray_NDIM(self); i++){
27422741 newstrides[i] = PyInt_AsLong(PyList_GetItem(value, Py_ssize_t(i)));
27432742 //newstrides_bytes[i] = newstrides[i] * 4;
27442743 }
27452744 }
27462745 /*
27472746 // Do not do this check, as ExtractDiag needs that, and NumPy does not seem
27482747 // to do it.
27492748 npy_intp dims[PyTuple_Size(value)];
27502749 for(int i=0; i < CudaNdarray_NDIM(self); i++){
27512750 dims[i] = CudaNdarray_HOST_DIMS(self)[i];
27522751 }
27532752 if (!PyArray_CheckStrides(4,
27542753 CudaNdarray_NDIM(self),
27552754 0, 0,
27562755 dims,
27572756 newstrides_bytes)){
27582757 PyErr_SetString(PyExc_ValueError, "bad new strides");
27592758 return -1;
27602759 }
27612760 */
27622761 for(int i=0; i < CudaNdarray_NDIM(self); i++){
27632762 CudaNdarray_set_stride(self, i, newstrides[i]);
27642763 }
27652764 return 0;
27662765 }
27672766
27682767 static PyObject *
27692768 CudaNdarray_get_dev_data(CudaNdarray *self, void *closure)
27702769 {
27712770 float * p = CudaNdarray_DEV_DATA(self);
27722771 //printf("get_dev_data %p %li \n", p, (long int)p );
27732772 return PyInt_FromSize_t((size_t) CudaNdarray_DEV_DATA(self));
27742773 }
27752774
27762775 static int
27772776 CudaNdarray_set_dev_data(CudaNdarray *self, PyObject *value, void *closure)
27782777 {
27792778 Py_ssize_t newdevdata = PyInt_AsSsize_t(value);
27802779 //printf("set_dev_data %p %li \n",(float*)newdevdata ,newdevdata);
27812780 if (PyErr_Occurred())
27822781 {
27832782 return -1;
27842783 }
27852784 return CudaNdarray_set_device_data(self, (float*)newdevdata, (CudaNdarray*)self->base);
27862785 }
27872786
27882787 static PyObject *
27892788 CudaNdarray_get_dtype(CudaNdarray *self, void *closure)
27902789 {
27912790 return PyString_FromString("float32");
27922791 }
27932792
27942793 static PyObject *
27952794 CudaNdarray_get_ndim(CudaNdarray *self, void *closure)
27962795 {
27972796 return PyInt_FromLong(self->nd);
27982797 }
27992798
28002799 static PyObject *
28012800 CudaNdarray_get_base(CudaNdarray *self, void *closure)
28022801 {
28032802 PyObject * base = self->base;
28042803 if (!base)
28052804 {
28062805 // We cannot return a NULL pointer, use None instead
28072806 base = Py_None;
28082807 }
28092808 Py_INCREF(base);
28102809 return base;
28112810 }
28122811
28132812 void put_in_dict(PyObject * dict, const char * key, int val)
28142813 {
28152814 PyObject * k = PyString_FromString(key);
28162815 PyObject * v = PyInt_FromLong(val);
28172816 PyDict_SetItem(dict, k, v);
28182817 Py_DECREF(k);
28192818 Py_DECREF(v);
28202819 }
28212820
28222821 PyObject *
28232822 GetDeviceProperties(PyObject* _unused, PyObject* args)
28242823 {
28252824 int dev_id = -1;
28262825 if (! PyArg_ParseTuple(args, "i", &dev_id))
28272826 return NULL;
28282827 cudaDeviceProp deviceProp;
28292828 cudaGetDeviceProperties(&deviceProp, dev_id);
28302829
28312830 PyObject * dict = PyDict_New();
28322831 PyObject * str= PyString_FromString("name");
28332832 PyObject * i = PyString_FromString(deviceProp.name);
28342833 PyDict_SetItem(dict, str, i);
28352834 Py_DECREF(str);
28362835 Py_DECREF(i);
28372836
28382837 put_in_dict(dict, "major", deviceProp.major);
28392838 put_in_dict(dict, "minor", deviceProp.minor);
28402839 #if CUDART_VERSION >= 2020
28412840 int driverVersion = 0, runtimeVersion = 0;
28422841 cudaDriverGetVersion(&driverVersion);
28432842 cudaRuntimeGetVersion(&runtimeVersion);
28442843 put_in_dict(dict, "driverVersion", driverVersion);
28452844 put_in_dict(dict, "runtimeVersion", runtimeVersion);
28462845 #endif
28472846 #if CUDART_VERSION >= 2000
28482847
28492848 put_in_dict(dict, "multiProcessorCount", deviceProp.multiProcessorCount);
28502849 //if ConvertSMVer2Cores is not defined in cuda_runtime_api.h, the run time is too old.
28512850 int sm_cores = -1;
28522851 if(deviceProp.major==1)
28532852 sm_cores = 32;
28542853 else if(deviceProp.major==2 && deviceProp.minor==0)
28552854 sm_cores = 32;
28562855 else if(deviceProp.major==2 && deviceProp.minor==1)
28572856 sm_cores = 48;
28582857 put_in_dict(dict, "coresCount", sm_cores * deviceProp.multiProcessorCount);
28592858 #endif
28602859 put_in_dict(dict, "totalConstMem", deviceProp.totalConstMem);
28612860 put_in_dict(dict, "sharedMemPerBlock", deviceProp.sharedMemPerBlock);
28622861 put_in_dict(dict, "regsPerBlock", deviceProp.regsPerBlock);
28632862 put_in_dict(dict, "warpSize", deviceProp.warpSize);
28642863 put_in_dict(dict, "maxThreadsPerBlock", deviceProp.maxThreadsPerBlock);
28652864 put_in_dict(dict, "maxThreadsDim0", deviceProp.maxThreadsDim[0]);
28662865 put_in_dict(dict, "maxThreadsDim1", deviceProp.maxThreadsDim[1]);
28672866 put_in_dict(dict, "maxThreadsDim2", deviceProp.maxThreadsDim[2]);
28682867 put_in_dict(dict, "maxGridSize0", deviceProp.maxGridSize[0]);
28692868 put_in_dict(dict, "maxGridSize1", deviceProp.maxGridSize[1]);
28702869 put_in_dict(dict, "maxGridSize2", deviceProp.maxGridSize[2]);
28712870 put_in_dict(dict, "memPitch", deviceProp.memPitch);
28722871 put_in_dict(dict, "textureAlignment", deviceProp.textureAlignment);
28732872 put_in_dict(dict, "clockRate", deviceProp.clockRate);
28742873 #if CUDART_VERSION >= 2000
28752874 put_in_dict(dict, "deviceOverlap", deviceProp.deviceOverlap);
28762875 #endif
28772876 #if CUDART_VERSION >= 2020
28782877 put_in_dict(dict, "kernelExecTimeoutEnabled", deviceProp.kernelExecTimeoutEnabled);
28792878 put_in_dict(dict, "integrated", deviceProp.integrated);
28802879 put_in_dict(dict, "canMapHostMemory", deviceProp.canMapHostMemory);
28812880 put_in_dict(dict, "computeMode", deviceProp.computeMode);
28822881 //in the doc of this fct tell that 0 - Normal mode, 1 - only 1 context, 2 - no context
28832882 #endif
28842883 #if CUDART_VERSION >= 3000
28852884 put_in_dict(dict, "concurrentKernels", deviceProp.concurrentKernels);
28862885 #endif
28872886 #if CUDART_VERSION >= 3010
28882887 put_in_dict(dict, "ECCEnabled", deviceProp.ECCEnabled);
28892888 #endif
28902889 #if CUDART_VERSION >= 3020
28912890 put_in_dict(dict, "tccDriver", deviceProp.tccDriver);
28922891 #endif
28932892
28942893 return dict;
28952894 }
28962895
28972896 /*
28982897 * Returns in *free and *total respectively, the free and total amount of memory available for allocation by the device in bytes.
28992898 */
29002899 PyObject *
29012900 GetDeviceMemInfo(PyObject* _unused, PyObject* dummy)
29022901 {
29032902 size_t free = 0, total = 0;
29042903 if(g_gpu_context_active == 0){
29052904 PyErr_Format(PyExc_RuntimeError, "No gpu device selected yet. Please make sure the gpu device was initialized by Theano before.");
29062905 return NULL;
29072906 }
29082907
29092908 cudaError_t err = cudaMemGetInfo(&free, &total);
29102909 if (err != cudaSuccess){
29112910 // Clear the error flag, cudaMemGetInfo doesn't do it.
29122911 // Currently this returns the same thing as err, but if in future
29132912 // it returns something else I still don't see why we should ignore
29142913 // it. All we want to do here is reset the flag.
29152914 cudaGetLastError();
29162915 PyErr_Format(PyExc_RuntimeError,
29172916 "Error while getting memory info about the gpu: %s",
29182917 cudaGetErrorString(err));
29192918 return NULL;
29202919 }
29212920 return PyTuple_Pack(2, PyLong_FromLong(free), PyLong_FromLong(total));
29222921 }
29232922
29242923 /*
29252924 * Synchronize with all the gpu device stream.
29262925 */
29272926 PyObject *
29282927 CudaNdarray_synchronize(PyObject* _unused, PyObject* dummy)
29292928 {
29302929 CNDA_BEGIN_ALLOW_THREADS
29312930 cudaThreadSynchronize();
29322931 CNDA_END_ALLOW_THREADS
29332932 Py_INCREF(Py_None);
29342933 return Py_None;
29352934 }
29362935
29372936 /*
29382937 * Exist and return true if we link with cublas v2.
29392938 */
29402939 PyObject *
29412940 CudaNdarray_cublasv2(PyObject* _unused, PyObject* dummy)
29422941 {
29432942 Py_INCREF(Py_True);
29442943 return Py_True;
29452944 }
29462945
29472946 PyObject *
29482947 CudaNdarray_select_a_gpu(PyObject* _unused, PyObject* dummy)
29492948 {
29502949 void * rval = NULL;
29512950 cudaError_t err;
29522951 int num_gpus = 0;
29532952
29542953 err = cudaGetDeviceCount(&num_gpus);
29552954 if (cudaSuccess != err){
29562955 printf("ERR!\\n");
29572956 PyErr_Format(PyExc_RuntimeError,
29582957 "Not able to get number of GPUs (%s).",
29592958 cudaGetErrorString(err));
29602959 return NULL;
29612960 }
29622961
29632962 for (int device = 0; device < num_gpus; device++) {
29642963 cudaSetDevice(device);
29652964 err = cudaDeviceSynchronize(); // << CUDA context gets created here.
29662965 cudaGetLastError(); // reset the error state
29672966 if (cudaSuccess == err)
29682967 break;
29692968 }
29702969
29712970 if (cudaSuccess != err){
29722971 printf("ERR!\\n");
29732972 PyErr_Format(PyExc_RuntimeError,
29742973 "Not able to select available GPU from %d cards (%s).",
29752974 num_gpus, cudaGetErrorString(err));
29762975 return NULL;
29772976 }
29782977
29792978 Py_INCREF(Py_None);
29802979 return Py_None;
29812980 }
29822981
29832982 #if COMPUTE_GPU_MEM_USED
29842983 /*
29852984 * Return the size in bytes that Theano currently have allocated on the gpu.
29862985 */
29872986 PyObject *
29882987 GetTheanoAllocInfo(PyObject* _unused, PyObject* dummy)
29892988 {
29902989 PyObject* a = PyLong_FromLong(_allocated_size);
29912990 PyObject* b = PyLong_FromLong(_max_allocated_size);
29922991
29932992 PyObject* tuple = PyTuple_New(2);
29942993 PyTuple_SetItem(tuple, 0, a);
29952994 PyTuple_SetItem(tuple, 1, b);
29962995 return tuple;
29972996 }
29982997 #endif
29992998
30002999 static PyGetSetDef CudaNdarray_getset[] = {
30013000 {"shape",
30023001 (getter)CudaNdarray_get_shape,
30033002 (setter)CudaNdarray_set_shape,
30043003 "shape of this ndarray (tuple)",
30053004 NULL},
30063005 {"_strides",
30073006 (getter)CudaNdarray_get_strides,
30083007 (setter)CudaNdarray_set_strides,
30093008 "data pointer strides (in elements)",
30103009 NULL},
30113010 {"strides",
30123011 (getter)CudaNdarray_get_strides,
30133012 (setter)CudaNdarray_set_strides,
30143013 "data pointer strides (in elements)",
30153014 NULL},
30163015 //gpudata is needed to allow calling pycuda fct with CudaNdarray input.
30173016 {"gpudata",
30183017 (getter)CudaNdarray_get_dev_data,
30193018 NULL,
30203019 "device data pointer",
30213020 NULL},
30223021 {"_dev_data",
30233022 (getter)CudaNdarray_get_dev_data,
30243023 (setter)CudaNdarray_set_dev_data,
30253024 "device data pointer",
30263025 NULL},
30273026 {"dtype",
30283027 (getter)CudaNdarray_get_dtype,
30293028 NULL,
30303029 "The dtype of the element. Now always float32",
30313030 NULL},
30323031 {"size",
30333032 (getter)CudaNdarray_SIZE_Object,
30343033 NULL,
30353034 "The number of elements in this object.",
30363035 NULL},
30373036 //mem_size is neede for pycuda.elementwise.ElementwiseKernel Why do they use size and mem_size of the same value?
30383037 {"mem_size",
30393038 (getter)CudaNdarray_SIZE_Object,
30403039 NULL,
30413040 "The number of elements in this object.",
30423041 NULL},
30433042 {"ndim",
30443043 (getter)CudaNdarray_get_ndim,
30453044 NULL,
30463045 "The number of dimensions in this object.",
30473046 NULL},
30483047 {"base",
30493048 (getter)CudaNdarray_get_base,
30503049 NULL,
30513050 "If this ndarray is a view, base is the original ndarray.",
30523051 NULL},
30533052
30543053 {NULL, NULL, NULL, NULL} /* Sentinel */
30553054 };
30563055
30573056 PyObject *CudaNdarray_repr(PyObject *self)
30583057 {
30593058 CudaNdarray *object = (CudaNdarray *)self;
30603059 PyObject * np_object = CudaNdarray_CreateArrayObj(object);
30613060 PyObject * str = PyObject_Str((PyObject *) np_object);
30623061 char * cstr = PyString_AsString(str);
30633062 PyObject * out = PyString_FromFormat("%s%s%s",
30643063 "CudaNdarray(",
30653064 cstr,
30663065 ")");
30673066 Py_DECREF(str);
30683067 Py_DECREF(np_object);
30693068 #if PY_MAJOR_VERSION >= 3
30703069 // In Python 3 PyString_FromFormat return a Bytes object
30713070 PyObject* out2 = PyObject_Str(out);
30723071 Py_DECREF(out);
30733072 return out2;
30743073 #endif
30753074 return out;
30763075 }
30773076
30783077 static PyTypeObject CudaNdarrayType =
30793078 {
30803079 #if PY_MAJOR_VERSION >= 3
30813080 PyVarObject_HEAD_INIT(NULL, 0)
30823081 #else
30833082 PyObject_HEAD_INIT(NULL)
30843083 0, /*ob_size*/
30853084 #endif
30863085 "CudaNdarray", /*tp_name*/
30873086 sizeof(CudaNdarray), /*tp_basicsize*/
30883087 0, /*tp_itemsize*/
30893088 (destructor)CudaNdarray_dealloc, /*tp_dealloc*/
30903089 0, /*tp_print*/
30913090 0, /*tp_getattr*/
30923091 0, /*tp_setattr*/
30933092 0, /*tp_compare*/
30943093 CudaNdarray_repr, /*tp_repr*/
30953094 &CudaNdarrayNumberMethods, /*tp_as_number*/
30963095 0, /*tp_as_sequence*/
30973096 &CudaNdarrayMappingMethods,/*tp_as_mapping*/
30983097 0, /*tp_hash */
30993098 0, /*tp_call*/
31003099 0, /*tp_str*/
31013100 0, /*tp_getattro*/
31023101 0, /*tp_setattro*/
31033102 0, /*tp_as_buffer*/
31043103 #if PY_MAJOR_VERSION >= 3
31053104 // Py_TPFLAGS_CHECKTYPES is always true and was removed in Python 3.
31063105 Py_TPFLAGS_DEFAULT | Py_TPFLAGS_BASETYPE, /*tp_flags*/
31073106 #else
31083107 Py_TPFLAGS_DEFAULT | Py_TPFLAGS_BASETYPE | Py_TPFLAGS_CHECKTYPES, /*tp_flags*/
31093108 #endif
31103109 "CudaNdarray objects", /* tp_doc */
31113110 0, /* tp_traverse */
31123111 0, /* tp_clear */
31133112 0, /* tp_richcompare */
31143113 0, /* tp_weaklistoffset */
31153114 0, /* tp_iter */
31163115 0, /* tp_iternext */
31173116 CudaNdarray_methods, /* tp_methods */
31183117 CudaNdarray_members, /* tp_members */
31193118 CudaNdarray_getset, /* tp_getset */
31203119 0, /* tp_base */
31213120 0, /* tp_dict */
31223121 0, /* tp_descr_get */
31233122 0, /* tp_descr_set */
31243123 0, /* tp_dictoffset */
31253124 (initproc)CudaNdarray_init,/* tp_init */
31263125 0, /* tp_alloc */
31273126 CudaNdarray_new, /* tp_new */
31283127 };
31293128
31303129 static __global__ void get_gpu_ptr_size(int* dst)
31313130 {
31323131 dst[0] = sizeof(float*);
31333132 dst[1] = sizeof(int);
31343133 }
31353134
31363135 PyObject *
31373136 CudaNdarray_ptr_int_size(PyObject* _unused, PyObject* args)
31383137 {
31393138 int *gpu_data = (int*)device_malloc(sizeof(int)*2);
31403139 if(gpu_data == NULL){
31413140 return NULL;
31423141 }
31433142 get_gpu_ptr_size<<<1,1>>>(gpu_data);
31443143
31453144 cudaError_t cudaErr = cudaGetLastError();
31463145 if (cudaSuccess != cudaErr){
31473146
31483147 device_free(gpu_data);
31493148 return PyErr_Format(PyExc_RuntimeError,
31503149 "CudaNdarray_ptr_int_size: error when calling the gpu code. (%s)",
31513150 cudaGetErrorString(cudaErr));
31523151 }
31533152
31543153 // Transfer the result to cpu
31553154 int gpu_sizes[] = {-1,-1};
31563155 cublasStatus_t err;
31573156 err = cublasGetVector(2, sizeof(int), gpu_data, 1, gpu_sizes, 1);
31583157 device_free(gpu_data);
31593158
31603159 if (CUBLAS_STATUS_SUCCESS != err){
31613160 PyErr_SetString(PyExc_RuntimeError, "error copying data to from memory");
31623161 return NULL;
31633162 }
31643163 return Py_BuildValue("iiii", (int) gpu_sizes[0], (int)sizeof(float*),
31653164 (int)sizeof(int), (int) gpu_sizes[1]);
31663165 }
31673166
31683167 static int cublas_init();
31693168 static void cublas_shutdown();
31703169 // Initialize the gpu.
31713170 // Takes two optional parameters, the device number and if we should use cnmem.
31723171 // If the device number is provided, it sets that device to be the active device.
31733172 // If not provided (usually just to test whether the gpu is available at all),
31743173 // it does not set an active device.
31753174 // Raises EnvironmentError or ValueError (as appropriate) if the initialization failed.
31763175 // cnmem is threaded like a bool. If converted to 0, don't use cnmem. Otherwise, use it.
31773176 PyObject *
31783177 CudaNdarray_gpu_init(PyObject* _unused, PyObject* args)
31793178 {
31803179 int card_nb = 0;
31813180 int card_number_provided = 1;
31823181 float cnmem = 0; // Theano flag lib.cnmem
31833182 // if we're given something wildly invalid, this will throw a TypeError
31843183 if(!PyArg_ParseTuple(args, "|if", &card_nb, &cnmem))
31853184 return NULL;
31863185 if(cnmem)
31873186 g_use_cnmem = true;
31883187
31893188 if(PyTuple_Size(args) == 0) {
31903189 card_number_provided = 0;
31913190 card_nb = 0;
31923191 }
31933192
31943193 int deviceCount;
31953194 cudaError err = cudaGetDeviceCount(&deviceCount);
31963195 if(cudaSuccess != err) {
31973196 return PyErr_Format(PyExc_EnvironmentError,
31983197 "Unable to get the number of gpus available: %s",
31993198 cudaGetErrorString(cudaGetLastError()));
32003199 }
32013200
32023201 // as soon as the first successful call to a cuda* function is made, a
32033202 // gpu context has been created
32043203 g_gpu_context_active = 1;
32053204
32063205 if(deviceCount <= 0) {
32073206 return PyErr_Format(PyExc_EnvironmentError,
32083207 "Can't use the GPU, no devices support CUDA");
32093208 }
32103209 if(card_number_provided && (card_nb < 0 || card_nb > (deviceCount - 1))) {
32113210 return PyErr_Format(PyExc_ValueError,
32123211 "Bad device number %d. Only %d devices available.",
32133212 card_nb,
32143213 deviceCount);
32153214 }
32163215
32173216 cudaDeviceProp deviceProp;
32183217 err = cudaGetDeviceProperties(&deviceProp, card_nb);
32193218 if(cudaSuccess != err) {
32203219 return PyErr_Format(PyExc_EnvironmentError,
32213220 "Unable to get properties of gpu %i: %s",
32223221 card_nb,
32233222 cudaGetErrorString(cudaGetLastError()));
32243223 }
32253224
32263225 if(deviceProp.major == 9999 && deviceProp.minor == 9999 ){
32273226 return PyErr_Format(PyExc_EnvironmentError,
32283227 "There is no device that supports CUDA");
32293228 }
32303229
32313230 if(card_number_provided) {
32323231 err = cudaSetDevice(card_nb);
32333232 if(cudaSuccess != err) {
32343233 return PyErr_Format(PyExc_EnvironmentError,
32353234 "Unable to set device %i: %s",
32363235 card_nb,
32373236 cudaGetErrorString(cudaGetLastError()));
32383237 }
32393238 if (cublas_init() == -1)
32403239 return NULL;
32413240 }
32423241 if(card_number_provided && g_use_cnmem) {
32433242 size_t mem = 0;
32443243 if (cnmem > 1)
32453244 mem = cnmem * 1024 * 1024;
32463245 else{
32473246 // Clip to 95% to let memory for the driver.
32483247 // 98% didn't worked in some cases.
32493248 if (cnmem > .95){
32503249 cnmem = .95;
32513250 }
32523251 size_t free = 0, total = 0;
32533252 cudaError_t err = cudaMemGetInfo(&free, &total);
32543253 if (err != cudaSuccess){
32553254 // Clear the error flag, cudaMemGetInfo doesn't do it.
32563255 // Currently this returns the same thing as err, but if in future
32573256 // it returns something else I still don't see why we should ignore
32583257 // it. All we want to do here is reset the flag.
32593258 cudaGetLastError();
32603259 PyErr_Format(PyExc_RuntimeError,
32613260 "Error while getting memory info about the gpu: %s",
32623261 cudaGetErrorString(err));
32633262 return NULL;
32643263 }
32653264 mem = total * cnmem;
32663265 }
32673266 if(initCnmem(card_number_provided, card_nb, mem) == -1){
32683267 return NULL;
32693268 }
32703269 }
32713270
32723271 Py_INCREF(Py_None);
32733272 return Py_None;
32743273 }
32753274
32763275 PyObject *
32773276 CudaNdarray_active_device_number(PyObject* _unused, PyObject* _unused_args) {
32783277 // NB: No cuda error checking here; keeps things simple, and it's not
32793278 // really necessary.
32803279 int currentDevice;
32813280 cudaGetDevice(¤tDevice);
32823281 return PyInt_FromLong(currentDevice);
32833282 }
32843283
32853284 PyObject *
32863285 CudaNdarray_active_device_name(PyObject* _unused, PyObject* _unused_args) {
32873286 // NB: No cuda error checking here; keeps things simple, and it's not
32883287 // really necessary.
32893288 int currentDevice;
32903289 cudaGetDevice(¤tDevice);
32913290
32923291 cudaDeviceProp deviceProp;
32933292 cudaGetDeviceProperties(&deviceProp, currentDevice);
32943293 return PyString_FromString(deviceProp.name);
32953294 }
32963295
32973296 PyObject *
32983297 CudaNdarray_gpu_shutdown(PyObject* _unused, PyObject* _unused_args) {
32993298 // Don't handle errors here
33003299 cublas_shutdown();
33013300 g_gpu_context_active = 0; // context has now been closed down
33023301 if(g_use_cnmem) {
33033302 cnmemStatus_t status = cnmemFinalize();
33043303 if(status != CNMEM_STATUS_SUCCESS) {
33053304 fprintf(stderr, "CudaNdarray_gpu_shutdown: cnmemFinalize failed! Reason=%s\n",
33063305 cnmemGetErrorString(status));
33073306 if(status == CNMEM_STATUS_CUDA_ERROR) {
33083307 fprintf(stderr, " Cuda-Reason=%s\n",
33093308 cudaGetErrorString(cudaGetLastError()));
33103309 }
33113310 }
33123311 }
33133312
33143313 Py_INCREF(Py_None);
33153314 return Py_None;
33163315 }
33173316
33183317 /*
33193318 * This function is tested in theano/misc/test_pycuda_theano_simple.py
33203319 */
33213320 PyObject *
33223321 CudaNdarray_from_gpu_pointer(PyObject* _unused, PyObject* args)
33233322 {
33243323 int verbose = 0;
33253324 PyObject *gpu_ptr = NULL;
33263325 PyObject *shapes = NULL;
33273326 PyObject *strides = NULL;
33283327 PyObject *base = NULL;
33293328 PyObject *rval = NULL;
33303329
33313330 //args should consist of 3 python objects
33323331 //The first is the gpu ptr
33333332 //The second if the shape
33343333 //The third if the strides
33353334 if (! PyArg_ParseTuple(args, "OOOO", &gpu_ptr, &shapes, &strides, &base))
33363335 return NULL;
33373336
33383337 if (verbose) printf("In CudaNdarray_from_gpu_pointer\n");
33393338 if (!PyLong_Check(gpu_ptr))
33403339 {
33413340 PyErr_Format(PyExc_Exception, "CudaNdarray_from_gpu_pointer: The gpu pointor is not an long");
33423341 return NULL;
33433342 }
33443343
33453344 Py_ssize_t nd = PyObject_Length(shapes);
33463345 if (nd < 0)
33473346 {
33483347 PyErr_SetString(PyExc_TypeError, "CudaNdarray_from_gpu_pointer: Couldn't get length of second argument");
33493348 return NULL;
33503349 }
33513350 Py_ssize_t nd_stride = PyObject_Length(strides);
33523351 if (nd_stride < 0)
33533352 {
33543353 PyErr_SetString(PyExc_TypeError, "CudaNdarray_from_gpu_pointer: Couldn't get length of third argument");
33553354 return NULL;
33563355 }
33573356
33583357 if (nd != nd_stride)
33593358 {
33603359 PyErr_SetString(PyExc_TypeError, "CudaNdarray_from_gpu_pointer: We need the same number of shapes and strides");
33613360 return NULL;
33623361 }
33633362
33643363 rval = CudaNdarray_New();
33653364
33663365 if (CudaNdarray_set_nd((CudaNdarray *)rval, nd))
33673366 {
33683367 //CudaNdarray_set_nd set the error msg
33693368 return NULL;
33703369 }
33713370 // set gpu pointeur
33723371 assert(((CudaNdarray *)rval)->data_allocated == 0);
33733372 if (CudaNdarray_set_device_data((CudaNdarray *)rval, (float *)PyInt_AsLong(gpu_ptr), base))
33743373 {
33753374 PyErr_SetString(PyExc_TypeError, "CudaNdarray_from_gpu_pointer: Error while setting the gpu pointor");
33763375 return NULL;
33773376
33783377 }
33793378
33803379 // Set dims and strides
33813380 for (int i = nd-1; i >= 0; --i)
33823381 {
33833382 PyObject * idx = PyLong_FromLong(i);
33843383 if (idx == NULL)
33853384 {
33863385 PyErr_SetString(PyExc_Exception, "CudaNdarray_from_gpu_pointer: Couldn't make long object to loop over list/tuple");
33873386 return NULL;
33883387 }
33893388 PyObject* dim_ = PyObject_GetItem(shapes, idx);
33903389 PyObject* strd_ = PyObject_GetItem(strides, idx);
33913390 if (!PyInt_Check(dim_))
33923391 {
33933392 PyErr_Format(PyExc_Exception, "CudaNdarray_from_gpu_pointer: shapes[%d] is not an int", i);
33943393 return NULL;
33953394 }
33963395 if (!PyInt_Check(strd_))
33973396 {
33983397 PyErr_Format(PyExc_Exception, "CudaNdarray_from_gpu_pointer: strides[%d] is not an int", i);
33993398 return NULL;
34003399 }
34013400 int dim = PyInt_AsLong(dim_);
34023401 int strd = PyInt_AsLong(strd_);
34033402 CudaNdarray_set_stride((CudaNdarray *)rval, i, strd);
34043403 CudaNdarray_set_dim((CudaNdarray *)rval, i, dim);
34053404 Py_DECREF(idx);
34063405 Py_DECREF(dim_);
34073406 Py_DECREF(strd_);
34083407 }
34093408 if (verbose) printf("CudaNdarray_from_gpu_pointer normal return\n");
34103409 return rval;
34113410 }
34123411
34133412 PyObject *
34143413 CudaNdarray_Dot(PyObject* _unused, PyObject* args)
34153414 {
34163415 PyObject *l=NULL;
34173416 PyObject *r=NULL;
34183417 PyObject * rval = NULL;
34193418
34203419 //args should consist of two python objects ("OO")
34213420 if (! PyArg_ParseTuple(args, "OO", &l, &r))
34223421 return NULL;
34233422
34243423 if (!CudaNdarray_Check(l) || !CudaNdarray_Check(r))
34253424 {
34263425 PyErr_SetString(PyExc_TypeError, "CudaNdarray arguments required ");
34273426 goto CudaNdarray_dot_fail;
34283427 }
34293428 if (((CudaNdarray*)l)->nd != 2)
34303429 {
34313430 PyErr_SetString(PyExc_TypeError, "need 2d CudaNdarray arg for now");
34323431 goto CudaNdarray_dot_fail;
34333432 }
34343433 if (((CudaNdarray*)r)->nd != 2)
34353434 {
34363435 PyErr_SetString(PyExc_TypeError, "need 2d CudaNdarray arg for now");
34373436 goto CudaNdarray_dot_fail;
34383437 }
34393438 rval = CudaNdarray_New();
34403439 if (!rval)
34413440 {
34423441 goto CudaNdarray_dot_fail;
34433442 }
34443443 int dims[2];
34453444 dims[0] = CudaNdarray_HOST_DIMS((CudaNdarray*)l)[0];
34463445 dims[1] = CudaNdarray_HOST_DIMS((CudaNdarray*)r)[1];
34473446 if (CudaNdarray_alloc_contiguous((CudaNdarray*)rval, 2, dims))
34483447 {
34493448 goto CudaNdarray_dot_fail;
34503449 }
34513450 if (CudaNdarray_gemm(1.0, (CudaNdarray*)l, (CudaNdarray*)r, 0.0, (CudaNdarray*)rval))
34523451 {
34533452 goto CudaNdarray_dot_fail;
34543453 }
34553454
34563455 return rval;
34573456
34583457 CudaNdarray_dot_fail:
34593458 Py_XDECREF(rval);
34603459 return NULL;
34613460 }
34623461
34633462 static PyObject *
34643463 filter(PyObject* __unsed_self, PyObject *args) // args = (data, broadcastable, strict, storage)
34653464 {
34663465 /*
34673466 * TODO: DOC what this function should do in the various cases of
34683467 * What is 'strict' supposed to mean in the context of this function?
34693468 * What do we do with input that could be interpreted as matching the broadcastable pattern in strict vs. non-strict cases?
34703469 *
34713470 */
34723471 PyObject *py_data=NULL;
34733472 PyArrayObject * data = NULL;
34743473 int strict = 0;
34753474 PyObject * broadcastable=NULL;
34763475 PyObject * storage=NULL;
34773476 CudaNdarray * rval=NULL;
34783477
34793478 //Python object references which are provided to the caller are borrowed references
34803479 if (!PyArg_ParseTuple(args, "OOiO", &py_data, &broadcastable, &strict, &storage)) return NULL;
34813480
34823481 if (!PyTuple_Check(broadcastable)){
34833482 PyErr_SetString(PyExc_TypeError, "broadcastable arg should be a tuple of int.");
34843483 return NULL;
34853484 }
34863485 Py_INCREF(py_data);
34873486 Py_INCREF(broadcastable);
34883487
34893488 CudaNdarray * cnda = (CudaNdarray*)py_data;
34903489
34913490 if (strict || CudaNdarray_Check(py_data))
34923491 {
34933492 //TODO: support non-strict "casting" from a vt to the broadcastable/type/size that we need.
34943493 if (!CudaNdarray_Check(py_data))
34953494 {
34963495 Py_DECREF(py_data);
34973496 Py_DECREF(broadcastable);
34983497 PyErr_SetString(PyExc_TypeError, "strict mode requires CudaNdarray");
34993498 return NULL;
35003499 }
35013500 if (cnda->nd != PyTuple_Size(broadcastable))
35023501 {
35033502 Py_DECREF(py_data);
35043503 Py_DECREF(broadcastable);
35053504 PyErr_Format(PyExc_TypeError, "Wrong rank: %i vs %li", cnda->nd, (long)PyTuple_Size(broadcastable));
35063505 return NULL;
35073506 }
35083507 for (int i = 0; i < cnda->nd; ++i)
35093508 {
35103509 if ((CudaNdarray_HOST_DIMS(cnda)[i] > 1) && PyInt_AsLong(PyTuple_GetItem(broadcastable, Py_ssize_t(i))))
35113510 {
35123511 PyErr_Format(PyExc_TypeError, "Non-unit size in broadcastable vt dimension %i", i);
35133512 Py_DECREF(py_data);
35143513 Py_DECREF(broadcastable);
35153514 return NULL;
35163515 }else if (CudaNdarray_HOST_DIMS(cnda)[i] == 1 && CudaNdarray_HOST_STRIDES(cnda)[i] != 0){
35173516 PyErr_Format(PyExc_TypeError, "Non-zeros strides(%d) on dimension %d of size 1",
35183517 CudaNdarray_HOST_STRIDES(cnda)[i], i);
35193518 Py_DECREF(py_data);
35203519 Py_DECREF(broadcastable);
35213520 return NULL;
35223521 }
35233522 }
35243523 Py_DECREF(broadcastable);
35253524 return py_data;
35263525 }
35273526 else
35283527 {
35293528 data = (PyArrayObject*)PyArray_FromObject(py_data, REAL_TYPENUM, PyTuple_Size(broadcastable), PyTuple_Size(broadcastable));
35303529 if (!data)
35313530 {
35323531 //err message already defined
35333532 Py_DECREF(py_data);
35343533 Py_DECREF(broadcastable);
35353534 return NULL;
35363535 }
35373536 for (int i = 0; i < PyArray_NDIM(data); ++i)
35383537 {
35393538 if ((PyArray_DIMS(data)[i] > 1) && PyInt_AsLong(PyTuple_GetItem(broadcastable, Py_ssize_t(i))))
35403539 {
35413540 PyErr_Format(PyExc_TypeError, "Non-unit size in broadcastable dimension %i", i);
35423541 Py_DECREF(data);
35433542 Py_DECREF(py_data);
35443543 Py_DECREF(broadcastable);
35453544 return NULL;
35463545 }
35473546 }
35483547 if (storage && CudaNdarray_Check(storage))
35493548 {
35503549 rval = (CudaNdarray*) storage;
35513550 Py_INCREF(rval);
35523551 }
35533552 else
35543553 {
35553554 rval = (CudaNdarray*) CudaNdarray_New();
35563555 }
35573556 if (rval)
35583557 {
35593558 if (CudaNdarray_CopyFromArray(rval, data))
35603559 {
35613560 Py_DECREF(rval);
35623561 rval = NULL;
35633562 }
35643563 }
35653564 Py_DECREF(data);
35663565 Py_DECREF(py_data);
35673566 Py_DECREF(broadcastable);
35683567 return (PyObject*)rval;
35693568 }
35703569 }
35713570
35723571 //TODO-- CudaNdarray_Dot and CudaNdarray_active_device_name are following different capitalization conventions.
35733572 // Pick one and standardize it, this file is already annoying enough to grep through
35743573 static PyMethodDef module_methods[] = {
35753574 {"dimshuffle", CudaNdarray_Dimshuffle, METH_VARARGS, "Returns the dimshuffle of a CudaNdarray."},
35763575 {"dot", CudaNdarray_Dot, METH_VARARGS, "Returns the matrix product of two CudaNdarray arguments."},
35773576 {"gpu_init", CudaNdarray_gpu_init, METH_VARARGS, "Select the gpu card to use; also usable to test whether CUDA is available."},
35783577 {"select_a_gpu", CudaNdarray_select_a_gpu, METH_NOARGS, "Call this method if you want to select a GPU before gpu_init call and let the driver choose the GPU."},
35793578 {"active_device_name", CudaNdarray_active_device_name, METH_VARARGS, "Get the name of the active device."},
35803579 {"active_device_number", CudaNdarray_active_device_number, METH_VARARGS, "Get the number of the active device."},
35813580 {"gpu_shutdown", CudaNdarray_gpu_shutdown, METH_VARARGS, "Shut down the gpu."},
35823581 {"device_properties", GetDeviceProperties, METH_VARARGS, "Return a dictionary with the device properties."},
35833582 {"mem_info", GetDeviceMemInfo, METH_NOARGS, "Return a tuple with the free and total memory on the gpu in bytes."},
35843583 #if COMPUTE_GPU_MEM_USED
35853584 {"theano_allocated", GetTheanoAllocInfo, METH_NOARGS, "Return the size in bytes of memory Theano currently have allocated on the gpu."},
35863585 #endif
35873586 {"ptr_int_size", CudaNdarray_ptr_int_size, METH_VARARGS, "Return a tuple with the size of gpu pointer, cpu pointer and int in bytes."},
35883587 {"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."},
35893588 {"outstanding_mallocs", outstanding_mallocs, METH_VARARGS, "how many more mallocs have been called than free's"},
35903589 {"from_gpu_pointer", CudaNdarray_from_gpu_pointer, METH_VARARGS, "Used to create a CudaNdarray from already allocated memory on the gpu.(example by pycuda)"},
35913590 {"synchronize", CudaNdarray_synchronize, METH_NOARGS, "Used to synchronize the device"},
35923591 {"cublas_v2", CudaNdarray_cublasv2, METH_NOARGS,
35933592 "Used to know if this version of cuda_ndarray is linked with cublas v2."},
35943593 {NULL, NULL, NULL, NULL} /* Sentinel */
35953594 };
35963595
35973596 #define CNDA_MOD_NAME "cuda_ndarray"
35983597 #define CNDA_DOCSTRING "CUDA implementation of a numpy ndarray-like object."
35993598
36003599 #if PY_MAJOR_VERSION == 3
36013600 static struct PyModuleDef cuda_ndarray_moduledef =
36023601 {
36033602 PyModuleDef_HEAD_INIT,
36043603 CNDA_MOD_NAME,
36053604 CNDA_DOCSTRING,
36063605 -1, /* size of per-interpreter state of the module,
36073606 or -1 if the module keeps state in global variables. */
36083607 module_methods
36093608 };
36103609
36113610 PyMODINIT_FUNC
36123611 PyInit_cuda_ndarray(void)
36133612 #else
36143613 PyMODINIT_FUNC
36153614 initcuda_ndarray(void)
36163615 #endif
36173616 {
36183617 import_array();
36193618
36203619 PyObject* m;
36213620
36223621 if (PyType_Ready(&CudaNdarrayType) < 0) {
36233622 #if PY_MAJOR_VERSION == 3
36243623 return NULL;
36253624 #else
36263625 return;
36273626 #endif
36283627 }
36293628
36303629 #if PY_MAJOR_VERSION == 3
36313630 m = PyModule_Create(&cuda_ndarray_moduledef);
36323631 #else
36333632 m = Py_InitModule3(CNDA_MOD_NAME, module_methods, CNDA_DOCSTRING);
36343633 #endif
36353634
36363635 if (m == NULL) {
36373636 #if PY_MAJOR_VERSION == 3
36383637 return NULL;
36393638 #else
36403639 return;
36413640 #endif
36423641 }
36433642
36443643 Py_INCREF(&CudaNdarrayType);
36453644 PyModule_AddObject(m, "CudaNdarray", (PyObject *)&CudaNdarrayType);
36463645 #if COMPUTE_GPU_MEM_USED
36473646 for(int i=0;i<TABLE_SIZE;i++){
36483647 _alloc_size_table[i].ptr=NULL;
36493648 _alloc_size_table[i].size=0;
36503649 }
36513650 #endif
36523651 // cublasInit();
36533652 //if (0&&CUBLAS_STATUS_SUCCESS != cublasGetError())
36543653 //{
36553654 //std::cerr << "WARNING: initcuda_ndarray: error initializing device\n";
36563655 //}
36573656 if (0) //TODO: is this necessary?
36583657 {
36593658 int deviceId = 0; // TODO: what number goes here?
36603659 cudaSetDevice(deviceId);
36613660 cudaError_t err = cudaGetLastError();
36623661 if( cudaSuccess != err)
36633662 {
36643663 std::cerr << "Error in SetDevice:" << cudaGetErrorString(err) << "\n";
36653664 }
36663665 }
36673666
36683667 #if PY_MAJOR_VERSION == 3
36693668 return m;
36703669 #endif
36713670 }
36723671
36733672
36743673 //////////////////////////////////////
36753674 //
36763675 // C API FOR CudaNdarray
36773676 //
36783677 //////////////////////////////////////
36793678
36803679 int
36813680 CudaNdarray_Check(const PyObject * ob)
36823681 {
36833682 //TODO: doesn't work with inheritance
36843683 return CudaNdarray_CheckExact(ob);
36853684 }
36863685 int
36873686 CudaNdarray_CheckExact(const PyObject * ob)
36883687 {
36893688 return ((Py_TYPE(ob) == &CudaNdarrayType) ? 1 : 0);
36903689 }
36913690
36923691 PyObject *
36933692 CudaNdarray_New(int nd)
36943693 {
36953694 CudaNdarray *self = (CudaNdarray *)CudaNdarrayType.tp_alloc(&CudaNdarrayType, 0);
36963695 if (self == NULL)
36973696 {
36983697 PyErr_SetString(PyExc_RuntimeError, "CudaNdarray_New failed to allocate self");
36993698 return NULL;
37003699 }
37013700 CudaNdarray_null_init(self);
37023701
37033702 if (nd == 0)
37043703 {
37053704 self->nd = 0;
37063705 }
37073706 else if (nd > 0)
37083707 {
37093708 if (CudaNdarray_set_nd(self, nd))
37103709 {
37113710 Py_DECREF(self);
37123711 return NULL;
37133712 }
37143713 }
37153714 ++_outstanding_mallocs[1];
37163715 return (PyObject *)self;
37173716 }
37183717
37193718
37203719
37213720 //////////////////////////////
37223721 //
37233722 // Published helper functions
37243723 //
37253724 //////////////////////////////
37263725
37273726 static int
37283727 cublas_init()
37293728 {
37303729 cublasStatus_t err;
37313730 err = cublasCreate(&handle);
37323731 if (CUBLAS_STATUS_SUCCESS != err)
37333732 {
37343733 if(CUBLAS_STATUS_NOT_INITIALIZED == err)
37353734 PyErr_SetString(PyExc_RuntimeError,
37363735 "cublasCreate() returned this error "
37373736 "'the CUDA Runtime initialization failed'");
37383737 else if(CUBLAS_STATUS_ALLOC_FAILED == err)
37393738 PyErr_SetString(PyExc_RuntimeError,
37403739 "cublasCreate() returned this error "
37413740 "'the resources could not be allocated'");
37423741 else
37433742 PyErr_SetString(PyExc_RuntimeError,
37443743 "unknow error during returned by cublasCreate()");
37453744 return -1;
37463745 }
37473746 // Set the default stream as the one to execute on (default)
37483747 cublasSetStream(handle, NULL);
37493748 // Pointer to scalars are on the host (also default)
37503749 cublasSetPointerMode(handle, CUBLAS_POINTER_MODE_HOST);
37513750 #if CUDA_VERSION >= 5000
37523751 // atomics can be used in kernels to speed up operations (not default)
37533752 // This may lead to a slight variance from run to run in some operations
37543753 cublasSetAtomicsMode(handle, CUBLAS_ATOMICS_ALLOWED);
37553754 #endif
37563755 return 0;
37573756 }
37583757
37593758 static void
37603759 cublas_shutdown()
37613760 {
37623761 if (handle != NULL)
37633762 cublasDestroy(handle);
37643763 // No point in handling any errors here
37653764 handle = NULL;
37663765 }
37673766
37683767 int
37693768 CudaNdarray_CopyFromArray(CudaNdarray * self, PyArrayObject*obj)
37703769 {
37713770 int err = CudaNdarray_alloc_contiguous(self, PyArray_NDIM(obj),
37723771 PyArray_DIMS(obj));
37733772 if (err) {
37743773 return err;
37753774 }
37763775
37773776 int typenum = PyArray_TYPE(obj);
37783777 if (typenum != REAL_TYPENUM)
37793778 {
37803779 PyErr_SetString(PyExc_TypeError, "can only copy from float arrays");
37813780 return -1;
37823781 }
37833782 assert( 4 == PyArray_ITEMSIZE(obj));
37843783 PyArrayObject * py_src = (PyArrayObject *)PyArray_ContiguousFromAny(
37853784 (PyObject*)obj, typenum, self->nd, self->nd);
37863785 if (!py_src) {
37873786 return -1;
37883787 }
37893788 npy_intp py_src_size = PyArray_SIZE(py_src);
37903789 void *py_src_data = PyArray_DATA(py_src);
37913790 cudaError_t cerr;
37923791 CNDA_BEGIN_ALLOW_THREADS;
37933792 cerr = cudaMemcpy(self->devdata, py_src_data,
37943793 py_src_size * sizeof(real),
37953794 cudaMemcpyHostToDevice);
37963795 //CNDA_THREAD_SYNC; // unneeded because cudaMemcpy is blocking anyway
37973796 CNDA_END_ALLOW_THREADS;
37983797 if (cudaSuccess != cerr)
37993798 {
38003799 PyErr_Format(PyExc_RuntimeError,
38013800 "Cuda error '%s' while copying %lli data element"
38023801 " to device memory. str ptr=%p. dst ptr=%p",
38033802 cudaGetErrorString(cerr),
38043803 (long long)py_src_size,
38053804 py_src_data,
38063805 self->devdata);
38073806 Py_DECREF(py_src);
38083807 return -1;
38093808 }
38103809 Py_DECREF(py_src);
38113810 return 0;
38123811 }
38133812
38143813 PyObject *
38153814 CudaNdarray_new_nd(int nd)
38163815 {
38173816 CudaNdarray * rval = (CudaNdarray*) CudaNdarray_New();
38183817 if (!rval || CudaNdarray_set_nd(rval, nd))
38193818 {
38203819 Py_XDECREF(rval);
38213820 rval = NULL;
38223821 }
38233822 return (PyObject *) rval;
38243823 }
38253824
38263825
38273826 /**
38283827 * Initialize 'self' as a view of 'base', with memory storage 'data'
38293828 */
38303829
38313830 int CudaNdarray_set_device_data(CudaNdarray * self, float * data, PyObject * base)
38323831 {
38333832 if (self->data_allocated)
38343833 {
38353834 assert(self->devdata);
38363835 if (device_free(self->devdata))
38373836 {
38383837 self->devdata = NULL;
38393838 self->data_allocated = 0;
38403839 return -1;
38413840 }
38423841 }
38433842 // Get the original base object (base.base.base...)
38443843 PyObject * orig_base = base;
38453844 // base is not always a CudaNdarray. It can be a GpuArray from pycuda, ...
38463845 while (orig_base && CudaNdarray_Check(orig_base) && ((CudaNdarray*) orig_base)->base)
38473846 {
38483847 // base_base is itself a view
38493848 orig_base = ((CudaNdarray*) orig_base)->base;
38503849 }
38513850 //N.B. XDECREF and XINCREF are no-ops for NULL pointers
38523851 if (self->base != orig_base)
38533852 {
38543853 Py_XDECREF(self->base);
38553854 self->base = orig_base;
38563855 Py_XINCREF(self->base);
38573856 }
38583857 self->data_allocated = 0;
38593858 self->devdata = data;
38603859 return 0;
38613860 }
38623861
38633862 static __global__ void k_copy_1d(const int N, const float * x, const int sx, float * y, const int sy)
38643863 {
38653864 for (int i = threadIdx.x + blockIdx.x * blockDim.x; i < N; i += gridDim.x*blockDim.x)
38663865 {
38673866 y[i*sy] = x[i*sx];
38683867 }
38693868 }
38703869
38713870 // N1 through N4 are the size of y
38723871 static __global__ void k_copy_4d(const int N1,
38733872 const int N2, const int N3, const int N4,
38743873 const float * x, const int sx1, const int sx2, const int sx3,
38753874 const int sx4, float * y, const int sy1, const int sy2,
38763875 const int sy3, const int sy4)
38773876 {
38783877 // These must be made int instead of unsigned int due to a bug in nvcc
38793878 int bx = blockIdx.x;
38803879 int by = blockIdx.y;
38813880
38823881 for (int i = bx; i < N1; i += gridDim.x)
38833882 {
38843883 for (int j = by; j < N2; j += gridDim.y)
38853884 {
38863885 for (int k = threadIdx.x; k < N3; k += (int) blockDim.x)
38873886 {
38883887 for (int l = threadIdx.y; l < N4; l += (int) blockDim.y)
38893888 {
38903889 y[i * sy1 + j * sy2 + k * sy3 + l * sy4] =
38913890 x[i * sx1 + j * sx2 + k * sx3 + l * sx4];
38923891 }
38933892 }
38943893 }
38953894 }
38963895 }
38973896
38983897 //copy from other into self
38993898 int CudaNdarray_CopyFromCudaNdarray(CudaNdarray * self,
39003899 const CudaNdarray * other,
39013900 bool unbroadcast)
39023901 {
39033902 int verbose = 0;
39043903 if (verbose>1) fprintf(stderr, "CudaNdarray_CopyFromCudaNdarray\n");
39053904
39063905 //standard elemwise size checks
39073906 if (self->nd == -1)
39083907 {
39093908 PyErr_SetString(PyExc_TypeError,
39103909 "can't copy into un-initialized CudaNdarray");
39113910 return -1;
39123911 }
39133912 CudaNdarray * new_other = NULL;
39143913
39153914 if (self->nd < other->nd)
39163915 {
39173916 PyErr_Format(PyExc_NotImplementedError,
39183917 "CudaNdarray_CopyFromCudaNdarray: The number of dimensions of the "
39193918 "destination needs to be >= the number of dimensions of the "
39203919 "source. Got %d and %d.", self->nd, other->nd);
39213920 return -1;
39223921 }
39233922 else if (self->nd != other->nd)
39243923 {
39253924 new_other = (CudaNdarray *) CudaNdarray_View(other);
39263925 int added_dims = self->nd - other->nd;
39273926 int* pattern = (int*) alloca(self->nd * sizeof(int));
39283927 for(int i = 0; i < added_dims; i++)
39293928 pattern[i] = -1;
39303929 for(int i = 0; i < other->nd; i++)
39313930 pattern[i + added_dims] = i;
39323931 CudaNdarray_dimshuffle(new_other, self->nd, pattern);
39333932 other = new_other;
39343933 }
39353934 assert(self->nd == other->nd);
39363935 //standard elemwise dim checks (also compute total size)
39373936 unsigned int size = 1;
39383937 unsigned int size_source = 1;
39393938 for (int i = 0; i< self->nd; ++i)
39403939 {
39413940 if ((CudaNdarray_HOST_DIMS(self)[i] != CudaNdarray_HOST_DIMS(other)[i])
39423941 && (1!=CudaNdarray_HOST_DIMS(other)[i] || !unbroadcast) )
39433942 {
39443943 PyErr_Format(PyExc_ValueError,
39453944 "CudaNdarray_CopyFromCudaNdarray:"
39463945 " need same dimensions for dim %d,"
39473946 " destination=%d, source=%d",
39483947 i, CudaNdarray_HOST_DIMS(self)[i],
39493948 CudaNdarray_HOST_DIMS(other)[i]);
39503949 Py_XDECREF(new_other);
39513950 return -1;
39523951 }
39533952 size *= (unsigned int) CudaNdarray_HOST_DIMS(self)[i];
39543953 size_source *= (unsigned int) CudaNdarray_HOST_DIMS(other)[i];
39553954 }
39563955 if (0 == size)
39573956 {
39583957 Py_XDECREF(new_other);
39593958 return 0; //nothing to copy, we're done.
39603959 }
39613960 if (CudaNdarray_is_c_contiguous(self) &&
39623961 CudaNdarray_is_c_contiguous(other) &&
39633962 size == size_source)
39643963 {
39653964 if (verbose)
39663965 fprintf(stderr, "Copying contiguous vector with cublasScopy\n");
39673966
39683967 cublasStatus_t err;
39693968 err = cublasScopy(handle, size, CudaNdarray_DEV_DATA(other), 1,
39703969 CudaNdarray_DEV_DATA(self), 1);
39713970 CNDA_THREAD_SYNC;
39723971 Py_XDECREF(new_other);
39733972 if (CUBLAS_STATUS_SUCCESS != err)
39743973 {
39753974 PyErr_SetString(PyExc_RuntimeError, "Error copying memory");
39763975 return -1;
39773976 }
39783977 return 0;
39793978 }
39803979 //TODO: rewrite these copy operations to be more efficient
39813980 // See, for example the transpose example in the cuda_sdk.
39823981 switch (self->nd)
39833982 {
39843983 case 0: // scalar
39853984 {
39863985 // THIS CASE SHOULD NEVER HAPPEN BECAUSE SCALARS ARE ALWAYS C CONTIGUOUS
39873986 assert(0);
39883987 }; break;
39893988 case 1: // vector
39903989 {
39913990 if (verbose) fprintf(stderr, "Copying non-contiguous vector\n");
39923991 if (verbose) fprint_CudaNdarray(stderr, other);
39933992 unsigned int n_blocks = std::min(size,
39943993 (unsigned int)NUM_VECTOR_OP_BLOCKS);
39953994 unsigned int n_threads = std::min(ceil_intdiv(size, n_blocks),
39963995 (unsigned int)NUM_VECTOR_OP_THREADS_PER_BLOCK);
39973996 k_copy_1d<<<n_blocks, n_threads>>>(size,
39983997 CudaNdarray_DEV_DATA(other),
39993998 CudaNdarray_HOST_STRIDES(other)[0],
40003999 CudaNdarray_DEV_DATA(self),
40014000 CudaNdarray_HOST_STRIDES(self)[0]);
40024001 CNDA_THREAD_SYNC;
40034002 cudaError_t err = cudaGetLastError();
40044003 if( cudaSuccess != err)
40054004 {
40064005 PyErr_Format(PyExc_RuntimeError,
40074006 "Cuda error: %s: %s. (n_blocks=%i,"
40084007 " n_threads_per_block=%i)\n", "k_copy_1d",
40094008 cudaGetErrorString(err), n_blocks, n_threads);
40104009 Py_XDECREF(new_other);
40114010 return -1;
40124011 }
40134012 }; break;
40144013 case 4: // 4-tensor
40154014 {
40164015 if (verbose)
40174016 {
40184017 if (0 != fprint_CudaNdarray(stderr, other))
40194018 {
40204019 Py_XDECREF(new_other);
40214020 return -1;
40224021 }
40234022 }
40244023
40254024 // The blocks implement the looping over the first two axes so
40264025 // this needs to be (N1, N2)
40274026 dim3 n_blocks( std::min(CudaNdarray_HOST_DIMS(self)[0],
40284027 NUM_VECTOR_OP_BLOCKS),
40294028 std::min(CudaNdarray_HOST_DIMS(self)[1],
40304029 NUM_VECTOR_OP_BLOCKS));
40314030 // For the threads, just make as many as possible
40324031 dim3 n_threads( std::min( (unsigned int) CudaNdarray_HOST_DIMS(self)[2],
40334032 (unsigned int) NUM_VECTOR_OP_THREADS_PER_BLOCK),
40344033 std::min( (unsigned int) CudaNdarray_HOST_DIMS(self)[3],
40354034 (unsigned int) NUM_VECTOR_OP_THREADS_PER_BLOCK));
40364035
40374036 n_threads.x = std::min( (unsigned int) 32, (unsigned int) n_threads.x);
40384037 n_threads.y = std::min( n_threads.y, NUM_VECTOR_OP_THREADS_PER_BLOCK / n_threads.x);
40394038
40404039 k_copy_4d<<<n_blocks, n_threads>>>(
40414040 // size of y
40424041 (unsigned int) CudaNdarray_HOST_DIMS(self)[0], // N1
40434042 (unsigned int) CudaNdarray_HOST_DIMS(self)[1], // N2
40444043 (unsigned int) CudaNdarray_HOST_DIMS(self)[2], // N3
40454044 (unsigned int) CudaNdarray_HOST_DIMS(self)[3], // N4
40464045 CudaNdarray_DEV_DATA(other), // x
40474046 // x strides
40484047 CudaNdarray_HOST_STRIDES(other)[0],
40494048 CudaNdarray_HOST_STRIDES(other)[1],
40504049 CudaNdarray_HOST_STRIDES(other)[2],
40514050 CudaNdarray_HOST_STRIDES(other)[3],
40524051 CudaNdarray_DEV_DATA(self), // y
40534052 // y strides
40544053 CudaNdarray_HOST_STRIDES(self)[0],
40554054 CudaNdarray_HOST_STRIDES(self)[1],
40564055 CudaNdarray_HOST_STRIDES(self)[2],
40574056 CudaNdarray_HOST_STRIDES(self)[3]
40584057 );
40594058 CNDA_THREAD_SYNC;
40604059 cudaError_t err = cudaGetLastError();
40614060 if( cudaSuccess != err)
40624061 {
40634062 PyErr_Format(PyExc_RuntimeError,
40644063 "Cuda error: %s: %s.",
40654064 "k_copy_4d",
40664065 cudaGetErrorString(err));
40674066 Py_XDECREF(new_other);
40684067 return -1;
40694068 }
40704069 }; break;
40714070 default:
40724071 {
40734072 cudaError_t err = cudaGetLastError();
40744073 if(cudaSuccess != err){
40754074 PyErr_Format(PyExc_RuntimeError,
40764075 "Unexpected Cuda error: %s: %s\n",
40774076 "CudaNdarray_CopyFromCudaNdarray",
40784077 cudaGetErrorString(err));
40794078 Py_XDECREF(new_other);
40804079 return -1;
40814080 }
40824081
40834082 if (verbose)
40844083 fprintf(stderr,
40854084 "Copying with default version unbroadcast=%d\n",
40864085 unbroadcast);
40874086 // call worker routine
40884087 unsigned int threads_per_block = std::min(size,
40894088 (unsigned int)NUM_VECTOR_OP_THREADS_PER_BLOCK);
40904089 unsigned int n_blocks = std::min(ceil_intdiv(size, threads_per_block),
40914090 (unsigned int)NUM_VECTOR_OP_BLOCKS);
40924091 const CudaNdarray * cuda_dims = other;
40934092 if(unbroadcast)
40944093 cuda_dims = self;
40954094 //copy from other into self
40964095 k_elemwise_unary_rowmajor_copy<<<n_blocks, threads_per_block>>>(
40974096 size,
40984097 (unsigned int)other->nd,
40994098 (const int *)CudaNdarray_DEV_DIMS(cuda_dims),
41004099 (const float*)CudaNdarray_DEV_DATA(other),
41014100 (const int *)CudaNdarray_DEV_STRIDES(other),
41024101 CudaNdarray_DEV_DATA(self),
41034102 (const int *)CudaNdarray_DEV_STRIDES(self));
41044103 CNDA_THREAD_SYNC;
41054104 err = cudaGetLastError();
41064105 if(verbose>1)
41074106 fprintf(stderr,
41084107 "INFO k_elemwise_unary_rowmaj (n_blocks=%i,"
41094108 " n_threads_per_block=%i)\n",
41104109 n_blocks, threads_per_block);
41114110 if( cudaSuccess != err)
41124111 {
41134112 //fprint_CudaNdarray(stderr, self);
41144113 //fprint_CudaNdarray(stderr, other);
41154114 PyErr_Format(PyExc_RuntimeError,
41164115 "Cuda error: %s: %s. (n_blocks=%i,"
41174116 " n_threads_per_block=%i)\n",
41184117 "k_elemwise_unary_rowmajor_copy",
41194118 cudaGetErrorString(err), n_blocks,
41204119 threads_per_block);
41214120 Py_XDECREF(new_other);
41224121 return -1;
41234122 }
41244123 }
41254124 };
41264125 Py_XDECREF(new_other);
41274126 return 0;
41284127 }
41294128
41304129 int CudaNdarray_gemm(float alpha, const CudaNdarray * A, const CudaNdarray * B, float beta, CudaNdarray * C)
41314130 {
41324131 if (A->nd != 2)
41334132 {
41344133 PyErr_SetString(PyExc_ValueError, "non-matrix arg A to gemm");
41354134 return -1;
41364135 }
41374136 if (B->nd != 2)
41384137 {
41394138 PyErr_SetString(PyExc_ValueError, "non-matrix arg B to gemm");
41404139 return -1;
41414140 }
41424141 if (C->nd != 2)
41434142 {
41444143 PyErr_SetString(PyExc_ValueError, "non-matrix arg C to gemm");
41454144 return -1;
41464145 }
41474146
41484147 // We must allow dimensions to be zeros.
41494148 if ((CudaNdarray_HOST_DIMS(A)[1] != CudaNdarray_HOST_DIMS(B)[0])
41504149 || (CudaNdarray_HOST_DIMS(A)[0] != CudaNdarray_HOST_DIMS(C)[0])
41514150 || (CudaNdarray_HOST_DIMS(B)[1] != CudaNdarray_HOST_DIMS(C)[1]))
41524151 {
41534152 PyErr_Format(PyExc_ValueError, "dimension mismatch in args to gemm (%i,%i)x(%i,%i)->(%i,%i)",
41544153 CudaNdarray_HOST_DIMS(A)[0],
41554154 CudaNdarray_HOST_DIMS(A)[1],
41564155 CudaNdarray_HOST_DIMS(B)[0],
41574156 CudaNdarray_HOST_DIMS(B)[1],
41584157 CudaNdarray_HOST_DIMS(C)[0],
41594158 CudaNdarray_HOST_DIMS(C)[1]);
41604159 return -1;
41614160 }
41624161
41634162 // If matrix A or B has non-unit size and non-unit stride in both
41644163 // dimensions, we can make a copy.
41654164 CudaNdarray * A_new = NULL;
41664165 CudaNdarray * B_new = NULL;
41674166 if (((CudaNdarray_HOST_DIMS(A)[0] > 1)
41684167 && (CudaNdarray_HOST_STRIDES(A)[0] != 1)
41694168 && (CudaNdarray_HOST_DIMS(A)[1] > 1)
41704169 && (CudaNdarray_HOST_STRIDES(A)[1] != 1))
41714170 || (CudaNdarray_HOST_STRIDES(A)[0] < 0)
41724171 || (CudaNdarray_HOST_STRIDES(A)[1] < 0))
41734172 {
41744173 A_new = (CudaNdarray*) CudaNdarray_Copy(A);
41754174 if (!A_new)
41764175 return -1;
41774176 A = A_new;
41784177 }
41794178
41804179 if (((CudaNdarray_HOST_DIMS(B)[0] > 1)
41814180 && (CudaNdarray_HOST_STRIDES(B)[0] != 1)
41824181 && (CudaNdarray_HOST_DIMS(B)[1] > 1)
41834182 && (CudaNdarray_HOST_STRIDES(B)[1] != 1))
41844183 || (CudaNdarray_HOST_STRIDES(B)[0] < 0)
41854184 || (CudaNdarray_HOST_STRIDES(B)[1] < 0))
41864185 {
41874186 B_new = (CudaNdarray*) CudaNdarray_Copy(B);
41884187 if (!B_new)
41894188 {
41904189 // If A_new is NULL, meaning A was not copied nothing happens
41914190 Py_XDECREF(A_new);
41924191 return -1;
41934192 }
41944193 B = B_new;
41954194 }
41964195
41974196 // If matrix C has non-unit size and non-unit stride in both
41984197 // dimensions, or negative strides, we can't operate. We cannot copy
41994198 // C either, because the calling code will expect the result to be
42004199 // in the original C container.
42014200 if (((CudaNdarray_HOST_DIMS(C)[0] > 1)
42024201 && (CudaNdarray_HOST_STRIDES(C)[0] != 1)
42034202 && (CudaNdarray_HOST_DIMS(C)[1] > 1)
42044203 && (CudaNdarray_HOST_STRIDES(C)[1] != 1))
42054204 || (CudaNdarray_HOST_STRIDES(C)[0] < 0)
42064205 || (CudaNdarray_HOST_STRIDES(C)[1] < 0))
42074206 {
42084207 PyErr_Format(PyExc_AssertionError,
42094208 "non-unit or negative stride in gemm arg C (%i,%i) of shape (%i,%i)",
42104209 CudaNdarray_HOST_STRIDES(C)[0],
42114210 CudaNdarray_HOST_STRIDES(C)[1],
42124211 CudaNdarray_HOST_DIMS(C)[0],
42134212 CudaNdarray_HOST_DIMS(C)[1]);
42144213 Py_XDECREF(A_new);
42154214 Py_XDECREF(B_new);
42164215 return -1;
42174216 }
42184217
42194218 // the unit integer is divided logically into three fields of 4 bits
42204219 // the lowermost 4 bits encode the stride pattern of the output
42214220 // the next higher 4 bits encode the B variable (or y)
42224221 // the next higher 4 bits encode the C variable (or x)
42234222 //
42244223 // the stride pattern for each input is encoded as 0 for unit stride from col to col (Row major)
42254224 // 1 for unit stride from row to row (Col major)
42264225
42274226 // a stride of 0 implies a dimension of 1 - so we can actually define
42284227 // a stride of 0 as a 'unit' stride because gemm will never use it.
42294228 // If a dimension is 0, its stride will not be used either, so we can
42304229 // consider it a 'unit' stride too.
42314230 int unit = 0;
42324231 if (CudaNdarray_HOST_STRIDES(A)[1] == 1 || CudaNdarray_HOST_DIMS(A)[1] <= 1) {
42334232 unit |= (0x0 << 8);
42344233 } else if (CudaNdarray_HOST_STRIDES(A)[0] == 1 || CudaNdarray_HOST_DIMS(A)[0] <= 1) {
42354234 unit |= (0x1 << 8);
42364235 } else {
42374236 unit |= (0x2 << 8);
42384237 }
42394238 if (CudaNdarray_HOST_STRIDES(B)[1] == 1 || CudaNdarray_HOST_DIMS(B)[1] <= 1) {
42404239 unit |= (0x0 << 4);
42414240 } else if (CudaNdarray_HOST_STRIDES(B)[0] == 1 || CudaNdarray_HOST_DIMS(B)[0] <= 1) {
42424241 unit |= (0x1 << 4);
42434242 } else {
42444243 unit |= (0x2 << 4);
42454244 }
42464245 if (CudaNdarray_HOST_STRIDES(C)[1] == 1 || CudaNdarray_HOST_DIMS(C)[1] <= 1) {
42474246 unit |= (0x0 << 0);
42484247 } else if (CudaNdarray_HOST_STRIDES(C)[0] == 1 || CudaNdarray_HOST_DIMS(C)[0] <= 1) {
42494248 unit |= (0x1 << 0);
42504249 } else {
42514250 unit |= (0x2 << 0);
42524251 }
42534252
42544253 /* create appropriate strides for malformed matrices that are row or column
42554254 * vectors
42564255 */
42574256 int sa_0 = (CudaNdarray_HOST_DIMS(A)[0] > 1) ? CudaNdarray_HOST_STRIDES(A)[0] : CudaNdarray_HOST_DIMS(A)[1];
42584257 int sa_1 = (CudaNdarray_HOST_DIMS(A)[1] > 1) ? CudaNdarray_HOST_STRIDES(A)[1] : CudaNdarray_HOST_DIMS(A)[0];
42594258 int sb_0 = (CudaNdarray_HOST_DIMS(B)[0] > 1) ? CudaNdarray_HOST_STRIDES(B)[0] : CudaNdarray_HOST_DIMS(B)[1];
42604259 int sb_1 = (CudaNdarray_HOST_DIMS(B)[1] > 1) ? CudaNdarray_HOST_STRIDES(B)[1] : CudaNdarray_HOST_DIMS(B)[0];
42614260 int sc_0 = (CudaNdarray_HOST_DIMS(C)[0] > 1) ? CudaNdarray_HOST_STRIDES(C)[0] : CudaNdarray_HOST_DIMS(C)[1];
42624261 int sc_1 = (CudaNdarray_HOST_DIMS(C)[1] > 1) ? CudaNdarray_HOST_STRIDES(C)[1] : CudaNdarray_HOST_DIMS(C)[0];
42634262
42644263 float* a = CudaNdarray_DEV_DATA(A);
42654264 float* b = CudaNdarray_DEV_DATA(B);
42664265 float* c = CudaNdarray_DEV_DATA(C);
42674266 cublasOperation_t N = CUBLAS_OP_N;
42684267 cublasOperation_t T = CUBLAS_OP_T;
42694268 //std::cerr << (unit/256) MOD 16 << (unit / 16) MOD 16 << unit MOD 16<< '\\n';
42704269 // There should be no negative stride at that point
42714270 #define CHK_STRIDE_SGEMM(T0, T1, D0, D1, D2, a, x, sx, y, sy, b, z, sz) \
42724271 if (sx == 0){sx = 1;}\
42734272 if (sy == 0){sy = 1;}\
42744273 if (sz == 0){sz = 1;}\
42754274 if ((sx > 0) && (sy > 0) && (sz > 0)) { \
42764275 err = cublasSgemm(handle, T0, T1, D0, D1, D2, &a, x, sx, y, sy, &b, z, sz); \
42774276 } else { \
42784277 PyErr_SetString(PyExc_AssertionError, "negative stride to sGemm");\
42794278 Py_XDECREF(A_new);\
42804279 Py_XDECREF(B_new);\
42814280 return -1; \
42824281 }
42834282
42844283 cublasStatus_t err;
42854284 switch(unit)
42864285 {
42874286 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;
42884287 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;
42894288 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;
42904289 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;
42914290 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;
42924291 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;
42934292 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;
42944293 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;
42954294 default: PyErr_Format(PyExc_ValueError, "some matrix has no unit stride (unit=%x)", unit);
42964295 return -1;
42974296 };
42984297 CNDA_THREAD_SYNC;
42994298 Py_XDECREF(A_new);
43004299 Py_XDECREF(B_new);
43014300
43024301 if (CUBLAS_STATUS_SUCCESS != err)
43034302 {
43044303 PyErr_Format(PyExc_RuntimeError,
43054304 "cublasSgemm failed (%i) %s\n"
43064305 " unit=%x N=%d, c.dims=[%d %d], a.dim=[%d %d], alpha=%f, beta=%f, a=%p, b=%p, c=%p"
43074306 " sa_0=%d, sa_1=%d, sb_0=%d, sb_1=%d, sc_0=%d, sc_1=%d",
43084307 err, cublasGetErrorString(err),
43094308 unit, N,
43104309 CudaNdarray_HOST_DIMS(C)[0],
43114310 CudaNdarray_HOST_DIMS(C)[1],
43124311 CudaNdarray_HOST_DIMS(A)[0], CudaNdarray_HOST_DIMS(A)[1],
43134312 alpha, beta, a, b, c, sa_0, sa_1, sb_0, sb_1, sc_0, sc_1);
43144313
43154314 return -1;
43164315 }
43174316 return 0;
43184317 }
43194318
43204319 int CudaNdarray_sgemv(float alpha, const CudaNdarray * A, const CudaNdarray * B, float beta, CudaNdarray * C)
43214320 {
43224321 /**
43234322 * C <- alpha A B + beta C
43244323 * A : matrix
43254324 * B, C: vector
43264325 * alpha, beta: scalars
43274326 */
43284327 if (A->nd != 2) { PyErr_SetString(PyExc_ValueError, "non-matrix arg to gemv"); return -1; }
43294328 if (B->nd != 1) { PyErr_SetString(PyExc_ValueError, "non-vector arg to gemv"); return -1; }
43304329 if (C->nd != 1) { PyErr_SetString(PyExc_ValueError, "non-vector arg to gemv"); return -1; }
43314330
43324331 // We must allow dimensions to be zeros.
43334332 if ((CudaNdarray_HOST_DIMS(A)[1] != CudaNdarray_HOST_DIMS(B)[0])
43344333 || (CudaNdarray_HOST_DIMS(A)[0] != CudaNdarray_HOST_DIMS(C)[0]))
43354334 {
43364335 PyErr_Format(PyExc_ValueError, "dimension mismatch in args to gemv (%i,%i)x(%i)->(%i)",
43374336 CudaNdarray_HOST_DIMS(A)[0],
43384337 CudaNdarray_HOST_DIMS(A)[1],
43394338 CudaNdarray_HOST_DIMS(B)[0],
43404339 CudaNdarray_HOST_DIMS(C)[0]);
43414340 return -1;
43424341 }
43434342
43444343 // If matrix A has non-unit size and non-unit stride in both
43454344 // dimensions, or negative strides, we cannot operate, but we can
43464345 // make a copy.
43474346 CudaNdarray * A_new = NULL;
43484347 CudaNdarray * B_new = NULL;
43494348 if (((CudaNdarray_HOST_DIMS(A)[0] > 1)
43504349 && (CudaNdarray_HOST_STRIDES(A)[0] != 1)
43514350 && (CudaNdarray_HOST_DIMS(A)[1] > 1)
43524351 && (CudaNdarray_HOST_STRIDES(A)[1] != 1))
43534352 || (CudaNdarray_HOST_STRIDES(A)[0] < 0)
43544353 || (CudaNdarray_HOST_STRIDES(A)[1] < 0))
43554354 {
43564355 A_new = (CudaNdarray*) CudaNdarray_Copy(A);
43574356 if (!A_new)
43584357 return -1;
43594358 A = A_new;
43604359 }
43614360
43624361 // If vector B as a negative stride, we also have to make a copy.
43634362 if (CudaNdarray_HOST_STRIDES(B)[0] < 0)
43644363 {
43654364 B_new = (CudaNdarray*) CudaNdarray_Copy(B);
43664365 if (!B_new)
43674366 {
43684367 // If A was not copied, A_new is NULL, and Py_XDECREF does not
43694368 // do anything
43704369 Py_XDECREF(A_new);
43714370 return -1;
43724371 }
43734372 B = B_new;
43744373 }
43754374
43764375 // cudablas does not handle negative strides as expected
43774376 if ( (CudaNdarray_HOST_STRIDES(A)[0] < 0)
43784377 || (CudaNdarray_HOST_STRIDES(A)[1] < 0))
43794378 {
43804379 PyErr_Format(PyExc_ValueError, "illegal strides in args to gemv (%i,%i)",
43814380 CudaNdarray_HOST_STRIDES(A)[0],
43824381 CudaNdarray_HOST_STRIDES(A)[1]);
43834382 Py_XDECREF(A_new);
43844383 Py_XDECREF(B_new);
43854384 return -1;
43864385 }
43874386
43884387 /* create appropriate strides for malformed matrices that are row or column
43894388 * vectors
43904389 */
43914390 int sa_0 = (CudaNdarray_HOST_DIMS(A)[0] > 1) ? CudaNdarray_HOST_STRIDES(A)[0] : CudaNdarray_HOST_DIMS(A)[1];
43924391 int sa_1 = (CudaNdarray_HOST_DIMS(A)[1] > 1) ? CudaNdarray_HOST_STRIDES(A)[1] : CudaNdarray_HOST_DIMS(A)[0];
43934392 int sb_0 = (CudaNdarray_HOST_DIMS(B)[0] > 1) ? CudaNdarray_HOST_STRIDES(B)[0] : 1;
43944393 int sc_0 = (CudaNdarray_HOST_DIMS(C)[0] > 1) ? CudaNdarray_HOST_STRIDES(C)[0] : 1;
43954394
43964395 if (sa_0 == 0) sa_0 = 1;
43974396 if (sa_1 == 0) sa_1 = 1;
43984397
43994398 int used_dot = 0;
44004399
44014400 // This is important because we can end up not calling Sgemv at all
44024401 cublasStatus_t err = CUBLAS_STATUS_SUCCESS;
44034402 if (CudaNdarray_SIZE(C)) {
44044403 // A is row vector & alpha==1 & beta==0 -> use cublasSdot
44054404 if (CudaNdarray_HOST_DIMS(A)[0] == 1 && alpha==1.f && beta==0.f) {
44064405 //replace this with custom inner product kernel with alpha and beta parameter?
44074406 cublasPointerMode_t pmode;
44084407 //set pointer mode to make sure cublas not storing on host pointer
44094408 cublasGetPointerMode(handle, &pmode);
44104409 cublasSetPointerMode(handle, CUBLAS_POINTER_MODE_DEVICE);
44114410 err = cublasSdot(
44124411 handle, CudaNdarray_HOST_DIMS(A)[1],
44134412 CudaNdarray_DEV_DATA(A), sa_1,
44144413 CudaNdarray_DEV_DATA(B), sb_0,
44154414 CudaNdarray_DEV_DATA(C));
44164415 cublasSetPointerMode(handle, pmode);
44174416 used_dot = 1;
44184417 }
44194418 // A is row-contiguous | row vector
44204419 else if ((CudaNdarray_HOST_DIMS(A)[0] <= 1)
44214420 || ((CudaNdarray_HOST_STRIDES(A)[0] == 1)
44224421 && (CudaNdarray_HOST_STRIDES(A)[1] > 0)))
44234422 {
44244423 err = cublasSgemv(handle, CUBLAS_OP_N,
44254424 CudaNdarray_HOST_DIMS(A)[0], CudaNdarray_HOST_DIMS(A)[1],
44264425 &alpha,
44274426 CudaNdarray_DEV_DATA(A), sa_1,
44284427 CudaNdarray_DEV_DATA(B), sb_0,
44294428 &beta,
44304429 CudaNdarray_DEV_DATA(C), sc_0);
44314430 }
44324431 // A is column-contiguous | column vector
44334432 else if ((CudaNdarray_HOST_DIMS(A)[1] <= 1)
44344433 || ((CudaNdarray_HOST_STRIDES(A)[1] == 1)
44354434 && (CudaNdarray_HOST_STRIDES(A)[0] > 0)))
44364435 {
44374436 err = cublasSgemv(handle, CUBLAS_OP_T,
44384437 CudaNdarray_HOST_DIMS(A)[1], CudaNdarray_HOST_DIMS(A)[0],
44394438 &alpha,
44404439 CudaNdarray_DEV_DATA(A), sa_0,
44414440 CudaNdarray_DEV_DATA(B), sb_0,
44424441 &beta,
44434442 CudaNdarray_DEV_DATA(C), sc_0);
44444443 }
44454444 // A is non vector and have malformed strides
44464445 else
44474446 {
44484447 PyErr_Format(PyExc_AssertionError,
44494448 "Unexpected stride pattern in gemv: (%i, %i) x %i -> %i.\n"
44504449 "Shapes are: (%i, %i) x %i -> %i\n",
44514450 CudaNdarray_HOST_STRIDES(A)[0],
44524451 CudaNdarray_HOST_STRIDES(A)[1],
44534452 CudaNdarray_HOST_STRIDES(B)[0],
44544453 CudaNdarray_HOST_STRIDES(C)[0],
44554454 CudaNdarray_HOST_DIMS(A)[0],
44564455 CudaNdarray_HOST_DIMS(A)[1],
44574456 CudaNdarray_HOST_DIMS(B)[0],
44584457 CudaNdarray_HOST_DIMS(C)[0]);
44594458 Py_XDECREF(A_new);
44604459 Py_XDECREF(B_new);
44614460 return -1;
44624461 }
44634462 }
44644463
44654464 CNDA_THREAD_SYNC;
44664465 Py_XDECREF(A_new);
44674466 Py_XDECREF(B_new);
44684467
44694468 if (CUBLAS_STATUS_SUCCESS != err)
44704469 {
44714470 if (!used_dot)
44724471 {
44734472 PyErr_Format(PyExc_RuntimeError,
44744473 "cublasSgemv failed (%i)",
44754474 err);
44764475 } else {
44774476 PyErr_Format(PyExc_RuntimeError,
44784477 "cublasSdot failed (%i)",
44794478 err);
44804479 }
44814480 return -1;
44824481 }
44834482 return 0;
44844483 }
44854484
44864485 int CudaNdarray_sger(float alpha, const CudaNdarray * x, const CudaNdarray * y, CudaNdarray * A) {
44874486 if (x->nd != 1) { PyErr_SetString(PyExc_ValueError, "non-vector arg x to sger"); return -1; }
44884487 if (y->nd != 1) { PyErr_SetString(PyExc_ValueError, "non-vector arg y to sger"); return -1; }
44894488 if (A->nd != 2) { PyErr_SetString(PyExc_ValueError, "non-matrix arg A to sger"); return -1; }
44904489
44914490 if ((CudaNdarray_HOST_DIMS(A)[0] != CudaNdarray_HOST_DIMS(x)[0])
44924491 || (CudaNdarray_HOST_DIMS(A)[1] != CudaNdarray_HOST_DIMS(y)[0])) {
44934492 PyErr_Format(PyExc_ValueError,
44944493 "dimension mismatch in args to sger (%i)x(%i)->(%i,%i)",
44954494 CudaNdarray_HOST_DIMS(x)[0],
44964495 CudaNdarray_HOST_DIMS(y)[0],
44974496 CudaNdarray_HOST_DIMS(A)[0],
44984497 CudaNdarray_HOST_DIMS(A)[1]);
44994498 return -1;
45004499 }
45014500
45024501 int x_strides = CudaNdarray_HOST_STRIDES(x)[0];
45034502 CudaNdarray * x_new = NULL;
45044503 if(x_strides == 0){
45054504 if(CudaNdarray_HOST_DIMS(x)[0] != 1){
45064505 PyErr_Format(PyExc_RuntimeError,
45074506 "CudaNdarray_sger: Invalid input x (should not happen)."
45084507 " We received a CudaNdarray vector with a stride of 0"
45094508 " that has more than 1 element!");
45104509 return -1;
45114510 }
45124511 x_strides = 1;
45134512 } else if(x_strides < 0){
45144513 x_new = (CudaNdarray*) CudaNdarray_Copy(x);
45154514 x = x_new;
45164515 x_strides = CudaNdarray_HOST_STRIDES(x)[0];
45174516 }
45184517
45194518 int y_strides = CudaNdarray_HOST_STRIDES(y)[0];
45204519 CudaNdarray * y_new = NULL;
45214520 if(y_strides == 0){
45224521 if(CudaNdarray_HOST_DIMS(y)[0] != 1){
45234522 PyErr_Format(PyExc_RuntimeError,
45244523 "CudaNdarray_sger: Invalid input y (should not happen)."
45254524 " We received a CudaNdarray vector with a stride of 0"
45264525 " that has more than 1 elements!");
45274526 Py_XDECREF(x_new);
45284527 return -1;
45294528 }
45304529 y_strides = 1;
45314530 } else if(y_strides < 0){
45324531 y_new = (CudaNdarray*) CudaNdarray_Copy(y);
45334532 y = y_new;
45344533 y_strides = CudaNdarray_HOST_STRIDES(y)[0];
45354534 }
45364535
45374536 // Create appropriate strides if A is a row or column vector
45384537 int sa_0 = (CudaNdarray_HOST_DIMS(A)[0] > 1) ? CudaNdarray_HOST_STRIDES(A)[0]
45394538 : CudaNdarray_HOST_DIMS(A)[1];
45404539 int sa_1 = (CudaNdarray_HOST_DIMS(A)[1] > 1) ? CudaNdarray_HOST_STRIDES(A)[1]
45414540 : CudaNdarray_HOST_DIMS(A)[0];
45424541
45434542 // This is important because we can end up not calling Sger at all
45444543 cublasStatus_t err = CUBLAS_STATUS_SUCCESS;
45454544 if(CudaNdarray_SIZE(A)){
45464545 // If A is in col-major
45474546 if ((CudaNdarray_HOST_DIMS(A)[0] <= 1)
45484547 || ((CudaNdarray_HOST_STRIDES(A)[0] == 1)
45494548 && (CudaNdarray_HOST_STRIDES(A)[1] > 0)))
45504549 {
45514550 err = cublasSger(handle, CudaNdarray_HOST_DIMS(x)[0], CudaNdarray_HOST_DIMS(y)[0], &alpha,
45524551 CudaNdarray_DEV_DATA(x), x_strides,
45534552 CudaNdarray_DEV_DATA(y), y_strides,
45544553 CudaNdarray_DEV_DATA(A), sa_1);
45554554 }
45564555 // Since Sger expects A in col-major, we invert x and y to fake this.
45574556 else if ((CudaNdarray_HOST_DIMS(A)[1] <= 1)
45584557 || ((CudaNdarray_HOST_STRIDES(A)[1] == 1)
45594558 && (CudaNdarray_HOST_STRIDES(A)[0] > 0)))
45604559 {
45614560 err = cublasSger(handle, CudaNdarray_HOST_DIMS(y)[0], CudaNdarray_HOST_DIMS(x)[0], &alpha,
45624561 CudaNdarray_DEV_DATA(y), y_strides,
45634562 CudaNdarray_DEV_DATA(x), x_strides,
45644563 CudaNdarray_DEV_DATA(A), sa_0);
45654564 }
45664565 // A has to be either c- or f-contiguous, with no negative strides
45674566 else
45684567 {
45694568 PyErr_SetString(PyExc_NotImplementedError,
45704569 "non-contiguous A, or negative strides, in sger");
45714570 Py_XDECREF(x_new);
45724571 Py_XDECREF(y_new);
45734572 return -1;
45744573 }
45754574 }
45764575 CNDA_THREAD_SYNC;
45774576 Py_XDECREF(x_new);
45784577 Py_XDECREF(y_new);
45794578
45804579 if (CUBLAS_STATUS_SUCCESS != err)
45814580 {
45824581 PyErr_Format(PyExc_RuntimeError,
45834582 "cublasSger failed (%i)",
45844583 err);
45854584 return -1;
45864585 }
45874586
45884587 return 0;
45894588 }
45904589
45914590 /**
45924591 *
45934592 * Precondition:
45944593 * a->dim[d] == (dims_a[d]==0) ? (1 << log2_dims_a[d]) : dims_a[d]
45954594 * z->dim[d] == (z_str[d]==0) ? 1 : dims_a[d];
45964595 *
45974596 * TODO: templatize this function to support other reductions.
45984597 * All that needs to change is the initial value for sum, and the reduction operator.
45994598 */
46004599
46014600 static __global__ void kernel_reduce_sum(const unsigned int size_z,
46024601 const unsigned int nd,
46034602 const int * dims_a,
46044603 const int * log2_dims_a,
46054604 const int * a_str,
46064605 const float * a_data,
46074606 const int * z_str,
46084607 float * z_data)
46094608 {
46104609 const unsigned int idx = blockIdx.x * blockDim.x + threadIdx.x;
46114610 const unsigned int numThreads = blockDim.x * gridDim.x;
46124611
46134612 //structure data contains the strides and dimensions of both a and z
46144613 // a_dim[0], a_dim[1], ... a_dim[nd-1],
46154614 // a_log2dim[0], a_log2dim[1], ... a_log2dim[nd-1],
46164615 // a_str[0], ... a_str[nd-1],
46174616 // z_str[0], ... z_str[nd-1]
46184617 extern __shared__ int structure_data[];
46194618 for (unsigned int i = threadIdx.x; i < nd; i += blockDim.x)
46204619 {
46214620 structure_data[i+0*nd] = dims_a[i];
46224621 structure_data[i+1*nd] = log2_dims_a[i];
46234622 structure_data[i+2*nd] = a_str[i];
46244623 structure_data[i+3*nd] = z_str[i];
46254624 }
46264625 dims_a = structure_data;
46274626 log2_dims_a = structure_data + nd;
46284627 a_str = structure_data + 2*nd;
46294628 z_str = structure_data + 3*nd;
46304629
46314630 __syncthreads(); //wait for all the shared structure to be loaded
46324631
46334632 for (unsigned int i = idx; i < size_z; i += numThreads)
46344633 {
46354634 unsigned int ii = i;
46364635 const float * a_data_i = a_data;
46374636 float * z_data_i = z_data;
46384637 unsigned int n_reduce_elements = 1;
46394638 unsigned int n_reduce_dims = 0;
46404639 unsigned int reduce_dim0 = nd-1;
46414640
46424641
46434642 //In this loop, we locate the initial element of the slice that we'd like to reduce with this thread
46444643 // At the same time, we [re]calculate the size of that slice (n_reduce_elements)
46454644 for (unsigned int d = 0; d < nd; ++d)
46464645 {
46474646 if (a_str[d] && (!z_str[d])) // this means 'd' is a dimension we are reducing over
46484647 {
46494648 n_reduce_elements *= dims_a[d];
46504649 n_reduce_dims += 1;
46514650 reduce_dim0 = (d < reduce_dim0) ? d : reduce_dim0;
46524651 }
46534652 else //'d' is not a dimension that we are reducing over
46544653 {
46554654 unsigned int pos_d;
46564655 if (log2_dims_a[d]==-1) //TODO: when things are working, use this switch
46574656 {
46584657 // this branch is not preferred,
46594658 // because the manual said that integer mod and div operations are slow on gpu
46604659 pos_d = (ii % dims_a[d]);
46614660 ii = (ii / dims_a[d]);
46624661 }
46634662 else
46644663 {
46654664 pos_d = (ii & ((1 << log2_dims_a[d])-1)); //take the lower log2_dims bits
46664665 ii = (ii >> log2_dims_a[d]); //shift those lower log2_dims bits off of ii
46674666 }
46684667 a_data_i += pos_d * a_str[d];
46694668 z_data_i += pos_d * z_str[d];
46704669 }
46714670 }
46724671 // now we've got pointers a_data_i and z_data_i into element 0 of the slice over which we are reducing
46734672 // do a similar loop
46744673
46754674 float sum = 0.0f;
46764675 switch(n_reduce_dims)
46774676 {
46784677 case 0:
46794678 {
46804679 sum = a_data_i[0];
46814680 }
46824681 break;
46834682 case 1:
46844683 {
46854684 const int stride = a_str[reduce_dim0];
46864685 const float * a_data_i_max = a_data_i + dims_a[reduce_dim0] * stride;
46874686 while (a_data_i != a_data_i_max)
46884687 {
46894688 sum += a_data_i[0];
46904689 a_data_i += stride;
46914690 }
46924691 }
46934692 break;
46944693 case 2:
46954694 {
46964695 int rd = reduce_dim0+1;
46974696 for (; rd < nd; ++rd)
46984697 {
46994698 if (a_str[rd] && (!z_str[rd])) // this means 'rd' is a dimension we are reducing over
47004699 break;
47014700 }
47024701 const int stride0 = a_str[reduce_dim0];
47034702 const int stride1 = a_str[rd];
47044703 for (int ii = 0; ii < dims_a[rd]; ++ii)
47054704 {
47064705 const float * a_data_ri = a_data_i + ii * stride1;
47074706 const float * a_data_ri_max = a_data_ri + dims_a[reduce_dim0] * stride0;
47084707 while (a_data_ri != a_data_ri_max)
47094708 {
47104709 sum += a_data_ri[0];
47114710 a_data_ri += stride0;
47124711 }
47134712 }
47144713 };
47154714 break;
47164715 default:
47174716 {
47184717 for (unsigned int reduce_i = 0; reduce_i < n_reduce_elements; ++reduce_i)
47194718 {
47204719 //TODO: optimize this loop to work more like theano's Elemwise. It's serial code.
47214720 unsigned int reduce_ii = reduce_i;
47224721 const float * a_data_ri = a_data_i;
47234722
47244723 //This loop finds the element in the a slice to add.
47254724 for (unsigned int rd = reduce_dim0; rd < nd; ++rd)
47264725 {
47274726 unsigned int pos_d;
47284727 if (a_str[rd] && (!z_str[rd])) // this means 'd' is a dimension we are reducing over
47294728 {
47304729 if (log2_dims_a[rd]==-1)
47314730 {
47324731 // this branch is not preferred,
47334732 // because the manual said that integer mod and div operations are slow on gpu
47344733 pos_d = (reduce_ii % dims_a[rd]);
47354734 reduce_ii = (reduce_ii / dims_a[rd]);
47364735 }
47374736 else
47384737 {
47394738 pos_d = (reduce_ii & ((1 << log2_dims_a[rd])-1)); //take the lower log2_dims bits
47404739 reduce_ii = (reduce_ii >> log2_dims_a[rd]); //shift those lower log2_dims bits off of ii
47414740 }
47424741 a_data_ri += pos_d * a_str[rd];
47434742 }
47444743 }
47454744 sum += a_data_ri[0];
47464745 }
47474746 }
47484747 }
47494748 z_data_i[0] = sum;
47504749 }
47514750 }
47524751
47534752 static __global__ void kernel_reduce_sum_1011(
47544753 const unsigned int d0,
47554754 const unsigned int d1,
47564755 const unsigned int d2,
47574756 const unsigned int d3,
47584757 const float *A, const int sA0, const int sA1, const int sA2, const int sA3,
47594758 float * Z, const int sZ0)
47604759 {
47614760 const int threadCount = blockDim.x * blockDim.y * blockDim.z;
47624761 const int threadNum = threadIdx.z * blockDim.x * blockDim.y + threadIdx.y * blockDim.x + threadIdx.x;
47634762 extern __shared__ float buf[];
47644763 float mysum = 0.0f;
47654764
47664765 if (warpSize != 32)
47674766 {
47684767 return; //TODO: set error code
47694768 }
47704769
47714770 for (int i0 = threadIdx.z; i0 < d0; i0 += blockDim.z)
47724771 {
47734772 float Ai = A[i0 * sA0 + blockIdx.x * sA1 + threadIdx.y * sA2 + threadIdx.x * sA3];
47744773 mysum += Ai;
47754774 }
47764775 buf[threadNum] = mysum;
47774776 __syncthreads();
47784777
47794778 // rest of function is handled by one warp
47804779 if (threadNum < warpSize)
47814780 {
47824781 for (int i = threadNum + warpSize; i < threadCount; i += warpSize)
47834782 {
47844783 mysum += buf[i];
47854784 }
47864785 buf[threadNum] = mysum;
47874786 if (threadNum < 16)
47884787 {
47894788 //reduce so that threadNum 0 has the sum of everything
47904789 if(threadNum + 16 < threadCount) buf[threadNum] += buf[threadNum+16];
47914790 if(threadNum + 8 < threadCount) buf[threadNum] += buf[threadNum+8];
47924791 if(threadNum + 4 < threadCount) buf[threadNum] += buf[threadNum+4];
47934792 if(threadNum + 2 < threadCount) buf[threadNum] += buf[threadNum+2];
47944793 if(threadNum + 1 < threadCount) buf[threadNum] += buf[threadNum+1];
47954794 if (threadNum == 0)
47964795 {
47974796 Z[blockIdx.x*sZ0] = buf[0];
47984797 }
47994798 }
48004799 }
48014800 }
48024801 /**
48034802 * Dimensions in which the self has size 1 and A has size > 1 are considered summing dimensions
48044803 * 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.
48054804 */
48064805 int
48074806 CudaNdarray_reduce_sum(CudaNdarray * self, CudaNdarray * A)
48084807 {
48094808 int verbose = 0;
48104809 //check input rank
48114810 if (self->nd != A->nd)
48124811 {
48134812 PyErr_Format(PyExc_TypeError, "Rank mismatch in CudaNdarray_sum: %i vs %i", self->nd, A->nd);
48144813 return -1;
48154814 }
48164815 for (int i = 0; i < self->nd; ++i)
48174816 {
48184817 if ((CudaNdarray_HOST_DIMS(self)[i] > 1) && (CudaNdarray_HOST_DIMS(self)[i] != CudaNdarray_HOST_DIMS(A)[i]))
48194818 {
48204819 PyErr_Format(PyExc_TypeError, "Dimension mismatch in CudaNdarray_sum: self->dim[%i] == %i , A->dim[%i] = %i",
48214820 i, CudaNdarray_HOST_DIMS(self)[i], i, CudaNdarray_HOST_DIMS(A)[i]);
48224821 return -1;
48234822 }
48244823 }
48254824
48264825 int n_summations = (unsigned int)CudaNdarray_SIZE(self);
48274826 if (verbose)
48284827 {
48294828 std::cerr << "reduce_sum n_summations " << n_summations << '\n';
48304829 std::cerr << "reduce_sum nd " << self->nd << '\n';
48314830 fprint_CudaNdarray(stderr, A);
48324831 fprint_CudaNdarray(stderr, self);
48334832 }
48344833 if (0 && (A->nd == 4) //check to see if kernel_reduce_sum_1011 applies
48354834 && (CudaNdarray_HOST_DIMS(self)[0] == 1)
48364835 && (CudaNdarray_HOST_DIMS(self)[2] == 1)
48374836 && (CudaNdarray_HOST_DIMS(self)[3] == 1)
48384837 )
48394838 {
48404839 dim3 n_threads(CudaNdarray_HOST_DIMS(A)[3], CudaNdarray_HOST_DIMS(A)[2]);
48414840 dim3 n_blocks(CudaNdarray_HOST_DIMS(A)[1]);
48424841 while (n_threads.x * n_threads.y * n_threads.z < NUM_VECTOR_OP_THREADS_PER_BLOCK) ++n_threads.z;
48434842 n_threads.z -= 1;
48444843 if (n_threads.z > 64) n_threads.z = 64;
48454844 if (n_threads.z)
48464845 {
48474846 if (verbose) printf("trying kernel_reduce_sum_1011\n");
48484847 int n_shared = sizeof(float) * n_threads.x * n_threads.y * n_threads.z;
48494848 kernel_reduce_sum_1011<<<n_blocks, n_threads, n_shared>>>(
48504849 CudaNdarray_HOST_DIMS(A)[0],
48514850 CudaNdarray_HOST_DIMS(A)[1],
48524851 CudaNdarray_HOST_DIMS(A)[2],
48534852 CudaNdarray_HOST_DIMS(A)[3],
48544853 CudaNdarray_DEV_DATA(A),
48554854 CudaNdarray_HOST_STRIDES(A)[0],
48564855 CudaNdarray_HOST_STRIDES(A)[1],
48574856 CudaNdarray_HOST_STRIDES(A)[2],
48584857 CudaNdarray_HOST_STRIDES(A)[3],
48594858 CudaNdarray_DEV_DATA(self),
48604859 CudaNdarray_HOST_STRIDES(self)[1]);
48614860 CNDA_THREAD_SYNC;
48624861 if (cudaSuccess == cudaGetLastError()) return 0;
48634862 if (verbose) printf("failed, falling back to kernel_reduce_sum\n");
48644863 }
48654864 }
48664865
48674866 int n_threads_per_block = std::min(n_summations,
48684867 NUM_VECTOR_OP_THREADS_PER_BLOCK);
48694868 int n_blocks = std::min(ceil_intdiv(n_summations,n_threads_per_block),
48704869 NUM_VECTOR_OP_BLOCKS);
48714870 int n_structure_cache = self->nd * 4 * sizeof(int);
48724871
48734872 if (verbose)
48744873 {
48754874 std::cerr << "n_blocks, n_threads_per_block " << n_blocks << ' ' << n_threads_per_block << '\n';
48764875 }
48774876 assert (self->nd > 0);
48784877 assert (self->nd == A->nd);
48794878 kernel_reduce_sum<<<n_blocks, n_threads_per_block, n_structure_cache>>>(
48804879 n_summations,
48814880 self->nd,
48824881 CudaNdarray_DEV_DIMS(A),
48834882 CudaNdarray_DEV_LOG2DIMS(A),
48844883 CudaNdarray_DEV_STRIDES(A),
48854884 CudaNdarray_DEV_DATA(A),
48864885 CudaNdarray_DEV_STRIDES(self),
48874886 CudaNdarray_DEV_DATA(self));
48884887 CNDA_THREAD_SYNC;
48894888 cudaError_t err = cudaGetLastError();
48904889 if (cudaSuccess != err)
48914890 {
48924891 PyErr_Format(PyExc_RuntimeError, "Cuda error: %s: %s.\n", "kernel_reduce_sum", cudaGetErrorString(err));
48934892 return -1;
48944893 }
48954894 return 0;
48964895 }
48974896 int
48984897 CudaNdarray_reduce_prod(CudaNdarray * self, const CudaNdarray * A)
48994898 {
49004899 PyErr_SetString(PyExc_NotImplementedError, "");
49014900 return -1;
49024901 }
49034902 int
49044903 CudaNdarray_reduce_min(CudaNdarray * self, const CudaNdarray * A)
49054904 {
49064905 PyErr_SetString(PyExc_NotImplementedError, "");
49074906 return -1;
49084907 }
49094908 int
49104909 CudaNdarray_reduce_max(CudaNdarray * self, const CudaNdarray * A)
49114910 {
49124911 PyErr_SetString(PyExc_NotImplementedError, "");
49134912 return -1;
49144913 }
49154914
49164915
49174916 /**
49184917 *
49194918 * pattern is a permutation of [0, 1, ... self->nd-1] with the following twists:
49204919 * - an element 'd' of the permutation can be dropped if CudaNdarray_HOST_DIMS(self)[d] == 1
49214920 * - any number of '-1' elements can be in the pattern, and they will cause new ranks (with dim==1) to be inserted.
49224921 *
49234922 * 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:
49244923 * [4, 6, 1, 1, 5] (we dropped the original dim[2]==1, and inserted two singleton dimensions with the -1s.
49254924 */
49264925 int
49274926 CudaNdarray_dimshuffle(CudaNdarray * self, unsigned int len, const int * pattern)
49284927 {
49294928 //TODO: pass a workspace pointer to avoid the internal malloc
49304929 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.
49314930 int * newstrides = newdims + len;
49324931 int * dims_taken = newstrides + len;
49334932 if (!newdims)
49344933 {
49354934 PyErr_SetString(PyExc_MemoryError, "CudaNdarray_dimshuffle: Failed to allocate temporary space");
49364935 return -1;
49374936 }
49384937 for (int i = 0; i < self->nd; ++i)
49394938 {
49404939 dims_taken[i] = 0;
49414940 }
49424941 for (int i = 0; i < len; ++i)
49434942 {
49444943 if (pattern[i] < 0)
49454944 {
49464945 newdims[i] = 1;
49474946 newstrides[i] = 0;
49484947 }
49494948 else if(dims_taken[pattern[i]])
49504949 {
49514950 PyErr_Format(PyExc_ValueError, "Cudandarray_dimshuffle: invalid pattern for Cudandarray_dimshuffle. You used the dimensions %d multiple time",
49524951 pattern[i]);
49534952 free(newdims);
49544953 return -1;
49554954 }
49564955 else if (pattern[i]>= self->nd)
49574956 {
49584957 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",
49594958 pattern[i], self->nd);
49604959 free(newdims);
49614960 return -1;
49624961 }
49634962 else
49644963 {
49654964 newdims[i] = CudaNdarray_HOST_DIMS(self)[pattern[i]];
49664965 newstrides[i] = CudaNdarray_HOST_STRIDES(self)[pattern[i]];
49674966 dims_taken[pattern[i]] = 1;
49684967 }
49694968 }
49704969 //Check if we dropped not broadcastable dims
49714970 for (int i = 0; i < self->nd; ++i)
49724971 {
49734972 if (dims_taken[i]==0 && CudaNdarray_HOST_DIMS(self)[i]!=1)
49744973 {
49754974 PyErr_SetString(PyExc_ValueError, "Cudandarray_dimshuffle: You cannot drop a non-broadcastable dimension.");
49764975 free(newdims);
49774976 return -1;
49784977 }
49794978 }
49804979 //swap this structure in for the one in self, and sync to the card
49814980 if (CudaNdarray_set_nd(self, len))
49824981 {
49834982 free(newdims);
49844983 return -1;
49854984 }
49864985 for (int i = 0; i < len; ++i)
49874986 {
49884987 CudaNdarray_set_dim(self, i, newdims[i]);
49894988 CudaNdarray_set_stride(self, i, newstrides[i]);
49904989 }
49914990 if (cnda_copy_structure_to_device(self))
49924991 {
49934992 free(newdims);
49944993 return -1;
49954994 }
49964995 free(newdims);
49974996 return 0;
49984997 }
49994998
50004999
50015000
50025001 /**
50035002 *
50045003 * This is the function that bind to python.
50055004 * See CudaNdarray_dimshuffle to call from C.
50065005 * We use -1 to mean 'x' as in Tensor Dimshuffle.
50075006 */
50085007 PyObject *
50095008 CudaNdarray_Dimshuffle(PyObject* _unused, PyObject* args)
50105009 {
50115010 PyObject * self = NULL;
50125011 PyObject * pattern_object = NULL;
50135012 int * pattern = NULL;
50145013 PyObject * rval = NULL;
50155014 int success = -1;
50165015 //const int * dims = NULL;
50175016
50185017 //args should consist of two python objects ("OO")
50195018 if (! PyArg_ParseTuple(args, "OO", &self, &pattern_object))
50205019 return NULL;
50215020
50225021 if (!CudaNdarray_Check(self) )
50235022 {
50245023 PyErr_SetString(PyExc_TypeError, "First argument to cuda_ndarray.dimshuffle must be a CudaNdarray");
50255024 return NULL;
50265025 }
50275026
50285027 //parse pattern_object into int * pattern
50295028
50305029 Py_ssize_t pattern_dim = PyObject_Length(pattern_object);
50315030
50325031 if (pattern_dim < 0)
50335032 {
50345033 PyErr_SetString(PyExc_TypeError, "Couldn't get length of third argument to cuda_ndarray.dimshuffle");
50355034 return NULL;
50365035 }
50375036
50385037 pattern = (int *) malloc( pattern_dim * sizeof(int));
50395038
50405039 for (Py_ssize_t i = 0; i < pattern_dim; i++)
50415040 {
50425041 PyObject * idx = PyLong_FromLong(i);
50435042
50445043 if (idx == NULL)
50455044 {
50465045 PyErr_SetString(PyExc_Exception, "Couldn't make long object to loop over list/tuple");
50475046 goto CudaNdarray_dimshuffle_fail;
50485047 }
50495048
50505049 long elem_value = 0;
50515050
50525051 PyObject * elem = PyObject_GetItem(pattern_object, idx);
50535052
50545053 if (elem == NULL)
50555054 {
50565055 Py_XDECREF( elem);
50575056 PyErr_SetString(PyExc_ValueError, "Third argument to dimshuffle must be list or tuple of integers");
50585057 goto CudaNdarray_dimshuffle_fail;
50595058 }
50605059
50615060 elem_value = PyInt_AsLong(elem);
50625061
50635062 if (elem_value == -1 && PyErr_Occurred() )
50645063 {
50655064 Py_XDECREF(elem);
50665065 PyErr_SetString(PyExc_ValueError, "Third argument to dimshuffle must be list or tuple of integers");
50675066 goto CudaNdarray_dimshuffle_fail;
50685067 }
50695068
50705069 pattern[i] = elem_value;
50715070
50725071 Py_XDECREF( elem );
50735072 Py_XDECREF( idx );
50745073 }
50755074
50765075 //allocate rval
50775076 rval = (PyObject *) CudaNdarray_View((CudaNdarray *) self);
50785077
50795078 if (rval == NULL)
50805079 {
50815080 //CudaNdarray_New should have set the exception string
50825081 goto CudaNdarray_dimshuffle_fail;
50835082 }
50845083
50855084
50865085 //printf("pattern_dim: %d\n",pattern_dim);
50875086 //printf("pattern: %d %d\n",pattern[0],pattern[1]);
50885087 //dims = CudaNdarray_HOST_DIMS( (CudaNdarray *) self);
50895088 //printf("dims before: %d %d\n",dims[0],dims[1]);
50905089
50915090 success = CudaNdarray_dimshuffle((CudaNdarray *) rval, pattern_dim, pattern);
50925091
50935092 if (success != 0)
50945093 {
50955094 //Exception string should already be set by CudaNdarray_dimshuffle
50965095 goto CudaNdarray_dimshuffle_fail;
50975096 }
50985097
50995098 free(pattern);
51005099
51015100 return rval;
51025101
51035102 CudaNdarray_dimshuffle_fail:
51045103
51055104 if (pattern != NULL)
51065105 free(pattern);
51075106
51085107 Py_XDECREF(rval);
51095108 return NULL;
51105109 }
51115110
51125111
51135112 int
51145113 cnda_structure_size(int nd)
51155114 {
51165115 // dim0, dim1, ...
51175116 // str0, str1, ...
51185117 // log2(dim0), log2(dim1), ...
51195118 return nd + nd + nd;
51205119 }
51215120
51225121 const int *
51235122 CudaNdarray_HOST_DIMS(const CudaNdarray * self)
51245123 {
51255124 return self->host_structure;
51265125 }
51275126
51285127 const int *
51295128 CudaNdarray_HOST_STRIDES(const CudaNdarray * self)
51305129 {
51315130 return self->host_structure + self->nd;
51325131 }
51335132 const int *
51345133 CudaNdarray_HOST_LOG2DIMS(const CudaNdarray * self)
51355134 {
51365135 return self->host_structure + 2*self->nd;
51375136 }
51385137
51395138 int
51405139 CudaNdarray_EqualAndIgnore(CudaNdarray *cnda1, CudaNdarray *cnda2, int ignoreSync, int ignoreBase)
51415140 {
51425141 int verbose = 0;
51435142
51445143 if (!ignoreSync && cnda1->dev_structure_fresh != cnda2->dev_structure_fresh)
51455144 {
51465145 if(verbose) fprintf(stdout, "CUDANDARRAY_EQUAL FAILED : 1\n");
51475146 return 0;
51485147 }
51495148
51505149 if (cnda1->nd != cnda2->nd)
51515150 {
51525151 if(verbose) fprintf(stdout, "CUDANDARRAY_EQUAL FAILED : 2\n");
51535152 return 0;
51545153 }
51555154
51565155 for (int i=0; i < 2*cnda1->nd; i++)
51575156 {
51585157 if (cnda1->host_structure[i] != cnda2->host_structure[i])
51595158 {
51605159 if(verbose)
51615160 fprintf(stdout, "CUDANDARRAY_EQUAL : host_structure : %d, %d, %d\n", i, cnda1->host_structure[i], cnda2->host_structure[i]);
51625161 return 0;
51635162 }
51645163 }
51655164
51665165 if (!ignoreBase && cnda1->base != cnda2->base)
51675166 {
51685167 if(verbose) fprintf(stdout, "CUDANDARRAY_EQUAL FAILED : 4");
51695168 return 0;
51705169 }
51715170 else if (cnda1->data_allocated != cnda2->data_allocated)
51725171 {
51735172 if(verbose) fprintf(stdout, "CUDANDARRAY_EQUAL FAILED : 5");
51745173 return 0;
51755174 }
51765175 else if (cnda1->data_allocated && cnda1->devdata != cnda2->devdata)
51775176 {
51785177 if(verbose) fprintf(stdout, "CUDANDARRAY_EQUAL FAILED : 6");
51795178 // no need to check devdata if data is not allocated
51805179 return 0;
51815180 }
51825181
51835182 return 1;
51845183 }
51855184
51865185
51875186 int
51885187 CudaNdarray_Equal(CudaNdarray *cnda1, CudaNdarray *cnda2)
51895188 {
51905189 return CudaNdarray_EqualAndIgnore(cnda1, cnda2, 0, 0);
51915190 }
51925191
51935192 int
51945193 cnda_copy_structure_to_device(const CudaNdarray * self)
51955194 {
51965195 //If the device structure do not exists, create it.
51975196 //We allocate it here as we do not need it often.
51985197 //In fact, we need it so infrequently that we expect
51995198 //that most object won't need it. Not allocating it
52005199 //save a significant when creating object.
52015200 //This speed up a benchmark by 8% with the gc.
52025201 if (!self->dev_structure)
52035202 {
52045203 int struct_size = cnda_structure_size(self->nd);
52055204 if (struct_size)
52065205 {
52075206 self->dev_structure = (int*)device_malloc(struct_size* sizeof(int));
52085207 if (NULL == self->dev_structure)
52095208 {
52105209 return -1;
52115210 }
52125211 }
52135212 }
52145213 if (cublasSetVector(cnda_structure_size(self->nd),
52155214 sizeof(int),
52165215 self->host_structure,
52175216 1,
52185217 self->dev_structure,
52195218 1) != CUBLAS_STATUS_SUCCESS)
52205219 {
52215220 PyErr_SetString(PyExc_RuntimeError, "error copying structure to device memory");
52225221 return -1;
52235222 }
52245223 self->dev_structure_fresh = 1;
52255224 return 0;
52265225 }
52275226
52285227 const int *
52295228 CudaNdarray_DEV_DIMS(const CudaNdarray * self)
52305229 {
52315230 if (!self->dev_structure_fresh)
52325231 {
52335232 if (cnda_copy_structure_to_device(self))
52345233 return NULL;
52355234 }
52365235 return self->dev_structure;
52375236 }
52385237 const int *
52395238 CudaNdarray_DEV_STRIDES(const CudaNdarray * self)
52405239 {
52415240 if (!self->dev_structure_fresh)
52425241 {
52435242 if (cnda_copy_structure_to_device(self))
52445243 return NULL;
52455244 }
52465245 return self->dev_structure + self->nd;
52475246 }
52485247 const int *
52495248 CudaNdarray_DEV_LOG2DIMS(const CudaNdarray * self)
52505249 {
52515250 if (!self->dev_structure_fresh)
52525251 {
52535252 if (cnda_copy_structure_to_device(self))
52545253 return NULL;
52555254 }
52565255 return self->dev_structure + 2*self->nd;
52575256 }
52585257 float *
52595258 CudaNdarray_DEV_DATA(const CudaNdarray * self)
52605259 {
52615260 return self->devdata;
52625261 }
52635262
52645263 /**
52655264 * Return the number of elements in the ndarray (product of the dimensions)
52665265 */
52675266 size_t
52685267 CudaNdarray_SIZE(const CudaNdarray *self)
52695268 {
52705269 if (self->nd == -1) return 0;
52715270 size_t size = 1;
52725271 for (int i = 0; i < self->nd; ++i)
52735272 {
52745273 size *= CudaNdarray_HOST_DIMS(self)[i];
52755274 }
52765275 return size;
52775276 }
52785277
52795278 PyObject *
52805279 CudaNdarray_SIZE_Object(const CudaNdarray *self, void *closure)
52815280 {
52825281 return PyInt_FromLong(CudaNdarray_SIZE(self));
52835282 }
52845283
52855284 int CudaNdarray_set_device_data(CudaNdarray * self, float * data, const CudaNdarray * base)
52865285 {
52875286 return CudaNdarray_set_device_data(self, data, (PyObject *) base);
52885287 }
52895288
52905289 PyObject * CudaNdarray_IS_C_Contiguous(CudaNdarray * self)
52915290 {
52925291 return PyBool_FromLong(CudaNdarray_is_c_contiguous(self));
52935292 }
52945293
52955294 int fprint_CudaNdarray(FILE * fd, const CudaNdarray *self)
52965295 {
52975296 cudaError_t err = cudaGetLastError();
52985297 if( cudaSuccess != err)
52995298 {
53005299 PyErr_Format(PyExc_RuntimeError,
53015300 "Cuda error: %s: %s.",
53025301 "fprint_CudaNdarray was called with an uncleared error",
53035302 cudaGetErrorString(err));
53045303 return -1;
53055304 }
53065305 fprintf(fd, "CudaNdarray <%p, %p> nd=%i dev_structure_fresh=%d data_allocated=%d\n",
53075306 self, self->devdata, self->nd, self->dev_structure_fresh, self->data_allocated);
53085307 fprintf(fd, "\tHOST_DIMS: ");
53095308 for (int i = 0; i < self->nd; ++i)
53105309 {
53115310 fprintf(fd, "%i\t", CudaNdarray_HOST_DIMS(self)[i]);
53125311 }
53135312 fprintf(fd, "\n\tHOST_STRIDES: ");
53145313 for (int i = 0; i < self->nd; ++i)
53155314 {
53165315 fprintf(fd, "%i\t", CudaNdarray_HOST_STRIDES(self)[i]);
53175316 }
53185317
53195318 if (self->dev_structure)
53205319 {
53215320 int data=0;
53225321 fprintf(fd, "\n\tDEV_DIMS: ");
53235322 for (int i = 0; i < self->nd; ++i)
53245323 {
53255324 cublasGetVector(1, sizeof(int),
53265325 self->dev_structure+i, 1,
53275326 &data, 1);
53285327 fprintf(fd, "%i\t", data);
53295328 }
53305329 fprintf(fd, "\n\tDEV_STRIDES: ");
53315330 for (int i = 0; i < self->nd; ++i)
53325331 {
53335332 cublasGetVector(1, sizeof(int),
53345333 self->dev_structure + self->nd+i, 1,
53355334 &data, 1);
53365335 fprintf(fd, "%i \t", data);
53375336 }
53385337 fprintf(fd, "\n");
53395338 }
53405339 else
53415340 {
53425341 fprintf(fd, "\n\tdev_structure not allocated\n");
53435342 }
53445343
53455344 err = cudaGetLastError();
53465345 if( cudaSuccess != err)
53475346 {
53485347 PyErr_Format(PyExc_RuntimeError,
53495348 "Cuda error: %s: %s.",
53505349 "fprint_CudaNdarray",
53515350 cudaGetErrorString(err));
53525351 return -1;
53535352 }
53545353 return 0;
53555354 }
53565355
53575356
53585357 int CudaNdarray_prep_output(CudaNdarray ** arr, int nd,
53595358 const int * dims, int fortran)
53605359 {
53615360 bool allocated = false;
53625361 if (*arr == NULL)
53635362 {
53645363 // This allocates the metadata but not the data
53655364 *arr = (CudaNdarray *) CudaNdarray_new_nd(nd);
53665365 if (*arr == NULL)
53675366 return -1;
53685367 allocated = true;
53695368 }
53705369
53715370 if (CudaNdarray_alloc_contiguous(*arr, nd, dims, fortran))
53725371 {
53735372 if (allocated)
53745373 {
53755374 Py_DECREF(*arr);
53765375 *arr = NULL;
53775376 }
53785377 return -1;
53795378 }
53805379 return 0;
53815380 }
53825381
53835382
53845383 /*
53855384 Local Variables:
53865385 mode:c++
53875386 c-basic-offset:4
53885387 c-file-style:"stroustrup"
53895388 indent-tabs-mode:nil
53905389 fill-column:79
53915390 End:
53925391 */
53935392 // vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:textwidth=79 :
53945393
5395===============================
5396nvcc warning : The 'compute_20', 'sm_20', and 'sm_21' architectures are deprecated, and may be removed in a future release (Use -Wno-deprecated-gpu-targets to suppress warning).
5397nvcc fatal : The version ('80000') of the host compiler ('Apple clang') is not supported
5398
5399['nvcc', '-shared', '-O3', '-use_fast_math', '--compiler-bindir', '/usr/bin/clang', '-m64', '-Xcompiler', '-DCUDA_NDARRAY_CUH=mc72d035fdf91890f3b36710688069b2e,-DNPY_NO_DEPRECATED_API=NPY_1_7_API_VERSION,-fPIC,-fvisibility=hidden', '-Xlinker', '-rpath,/Users/david/.theano/compiledir_Darwin-16.0.0-x86_64-i386-64bit-i386-3.5.2-64/cuda_ndarray', '-Xlinker', '-rpath,/usr/local/cuda/lib', '-I/usr/local/lib/python3.5/site-packages/theano/sandbox/cuda', '-I/usr/local/lib/python3.5/site-packages/numpy/core/include', '-I/usr/local/Cellar/python3/3.5.2_3/Frameworks/Python.framework/Versions/3.5/include/python3.5m', '-I/usr/local/lib/python3.5/site-packages/theano/gof', '-L/usr/local/Cellar/python3/3.5.2_3/Frameworks/Python.framework/Versions/3.5/lib', '-o', '/Users/david/.theano/compiledir_Darwin-16.0.0-x86_64-i386-64bit-i386-3.5.2-64/cuda_ndarray/cuda_ndarray.so', 'mod.cu', '-lcublas', '-lcudart', '-Xcompiler', '-undefined,dynamic_lookup', '-Xlinker', '-pie']
5400ERROR (theano.sandbox.cuda): Failed to compile cuda_ndarray.cu: ('nvcc return status', 1, 'for cmd', 'nvcc -shared -O3 -use_fast_math --compiler-bindir /usr/bin/clang -m64 -Xcompiler -DCUDA_NDARRAY_CUH=mc72d035fdf91890f3b36710688069b2e,-DNPY_NO_DEPRECATED_API=NPY_1_7_API_VERSION,-fPIC,-fvisibility=hidden -Xlinker -rpath,/Users/david/.theano/compiledir_Darwin-16.0.0-x86_64-i386-64bit-i386-3.5.2-64/cuda_ndarray -Xlinker -rpath,/usr/local/cuda/lib -I/usr/local/lib/python3.5/site-packages/theano/sandbox/cuda -I/usr/local/lib/python3.5/site-packages/numpy/core/include -I/usr/local/Cellar/python3/3.5.2_3/Frameworks/Python.framework/Versions/3.5/include/python3.5m -I/usr/local/lib/python3.5/site-packages/theano/gof -L/usr/local/Cellar/python3/3.5.2_3/Frameworks/Python.framework/Versions/3.5/lib -o /Users/david/.theano/compiledir_Darwin-16.0.0-x86_64-i386-64bit-i386-3.5.2-64/cuda_ndarray/cuda_ndarray.so mod.cu -lcublas -lcudart -Xcompiler -undefined,dynamic_lookup -Xlinker -pie')
5401Traceback (most recent call last):
5402 File "test_theano.py", line 1, in <module>
5403 from theano import function, config, shared, sandbox
5404 File "/usr/local/lib/python3.5/site-packages/theano/__init__.py", line 107, in <module>
5405 import theano.sandbox.cuda
5406 File "/usr/local/lib/python3.5/site-packages/theano/sandbox/cuda/__init__.py", line 713, in <module>
5407 use(device=config.device, force=config.force_device, test_driver=False)
5408 File "/usr/local/lib/python3.5/site-packages/theano/sandbox/cuda/__init__.py", line 503, in use
5409 cuda_initialization_error_message))
5410OSError: You forced the use of gpu device gpu, but CUDA initialization failed with error:
5411cuda unavailable