· 9 years ago · Nov 28, 2016, 07:18 PM
1aetros start burgalon/digit-convolution
2
3...
4Training '1xv2KAPRD' created and started. Open http://aetros.com/trainer/app?training=1xv2KAPRD to monitor the training.
5start network ...
6Using Theano backend.
71 #define _CUDA_NDARRAY_C
82
93 #include <Python.h>
104 #include <structmember.h>
115 #include "theano_mod_helper.h"
126
137 #include <numpy/arrayobject.h>
148 #include <iostream>
159
1610 #include "cuda_ndarray.cuh"
1711
1812 #ifndef CNMEM_DLLEXPORT
1913 #define CNMEM_DLLEXPORT
2014 #endif
2115
2216 #include "cnmem.h"
2317 #include "cnmem.cpp"
2418
2519 //If true, when there is a gpu malloc or free error, we print the size of allocated memory on the device.
2620 #define COMPUTE_GPU_MEM_USED 0
2721
2822 //If true, we fill with NAN allocated device memory.
2923 #define ALLOC_MEMSET 0
3024
3125 //If true, we print out when we free a device pointer, uninitialize a
3226 //CudaNdarray, or allocate a device pointer
3327 #define PRINT_FREE_MALLOC 0
3428
3529 //If true, we do error checking at the start of functions, to make sure there
3630 //is not a pre-existing error when the function is called.
3731 //You probably need to set the environment variable
3832 //CUDA_LAUNCH_BLOCKING=1, and/or modify the CNDA_THREAD_SYNC
3933 //preprocessor macro in cuda_ndarray.cuh
4034 //if you want this to work.
4135 #define PRECHECK_ERROR 0
4236
4337 cublasHandle_t handle = NULL;
4438 int* err_var = NULL;
4539
4640 /////////////////////////
4741 // Alloc and Free
4842 /////////////////////////
4943
5044 static int g_gpu_context_active = 0;
5145
5246
5347 PyObject *
5448 CudaNdarray_Dimshuffle(PyObject* _unused, PyObject* args);
5549 static PyObject *CudaNdarray_get_shape(CudaNdarray *self, void *closure);
5650
5751
5852 /**
5953 *
6054 * In the test program I'm using, the _outstanding_mallocs decreases with every call.
6155 * This suggests there are more free() calls being made than alloc(), but I can't figure out why.
6256 *
6357 */
6458 int _outstanding_mallocs[] = {0,0};
6559
6660 #if COMPUTE_GPU_MEM_USED
6761 size_t _allocated_size = 0;
6862 size_t _max_allocated_size = 0;
6963
7064 const int TABLE_SIZE = 10000;
7165 struct table_struct{
7266 void* ptr;
7367 size_t size;
7468 };
7569 table_struct _alloc_size_table[TABLE_SIZE];
7670 #endif
7771
7872 void * device_malloc(size_t size)
7973 {
8074 return device_malloc(size, VERBOSE_DEVICE_MALLOC);
8175 }
8276
8377 ///@TODO: thejaswi: link this option to a theano config variable?
8478 static bool g_use_cnmem = false;
8579 static const int g_max_devices = 8;
8680 int initCnmem(int card_number_provided, int card_nb, size_t mem) {
8781 static bool cnmemInitialized = false;
8882 if(cnmemInitialized) {
8983 return 0;
9084 }
9185 // On stderr to be at the same place as "Using gpu device..."
9286 int numDevices = 0;
9387 cnmemDevice_t devices[g_max_devices];
9488 if(cudaGetDeviceCount(&numDevices) != cudaSuccess) {
9589 PyErr_Format(PyExc_RuntimeError,
9690 "initCnmem: 'cudaGetDeviceCount' failed! Reason=%s\n",
9791 cudaGetErrorString(cudaGetLastError()));
9892 return -1;
9993 }
10094 if(card_number_provided){
10195 numDevices = 1;
10296 int i = 0;
10397 devices[i].device = card_nb;
10498 devices[i].size = mem;
10599 ///@TODO: thejaswi: add support for multiple streams
106100 devices[i].numStreams = 0;
107101 devices[i].streams = NULL;
108102 devices[i].streamSizes = NULL;
109103 }else{
110104 for(int i=0;i<numDevices;++i) {
111105 devices[i].device = i;
112106 devices[i].size = mem;
113107 ///@TODO: thejaswi: add support for multiple streams
114108 devices[i].numStreams = 0;
115109 devices[i].streams = NULL;
116110 }
117111 }
118112
119113 ///@TODO: thejaswi: passing custom cnmem flags?
120114 cnmemStatus_t status = cnmemInit(numDevices, devices, CNMEM_FLAGS_DEFAULT);
121115 if(status != CNMEM_STATUS_SUCCESS) {
122116 PyErr_Format(PyExc_RuntimeError,
123117 "initCnmem: cnmemInit call failed! Reason=%s. numdev=%d\n",
124118 cnmemGetErrorString(status), numDevices);
125119 return -1;
126120 }
127121 cnmemInitialized = true;
128122 return 0;
129123 }
130124
131125 void * device_malloc(size_t size, int verbose)
132126 {
133127 #if PRECHECK_ERROR
134128 cudaThreadSynchronize();
135129 cudaError_t prevError = cudaGetLastError();
136130 if (cudaSuccess != prevError)
137131 {
138132 fprintf(stderr,
139133 "Error existed before calling device_malloc. %s\n",
140134 cudaGetErrorString(prevError)
141135 );
142136 }
143137 #endif
144138 void * rval=NULL;
145139 ///@TODO: thejaswi: support for multiple-streams?
146140 if(g_use_cnmem) {
147141 cnmemStatus_t status = CNMEM_STATUS_SUCCESS;
148142 status = cnmemMalloc(&rval, size, NULL);
149143 if(status != CNMEM_STATUS_SUCCESS) {
150144 PyErr_Format(PyExc_MemoryError,
151145 "Error allocating %llu bytes of device memory (%s).",
152146 (unsigned long long)size, cnmemGetErrorString(status));
153147 return NULL;
154148 }
155149 }
156150 else {
157151 cudaError_t err = cudaMalloc(&rval, size);
158152 if (cudaSuccess != err)
159153 {
160154 // Clear the error flag, cudaMalloc doesn't do it.
161155 // Currently this returns the same thing as err, but if in future
162156 // it returns something else I still don't see why we should ignore
163157 // it. All we want to do here is reset the flag.
164158 cudaGetLastError();
165159 if (verbose)
166160 {
167161 size_t free = 0, total = 0;
168162 cudaError_t err2 = cudaMemGetInfo(&free, &total);
169163 if (err2 != cudaSuccess){
170164 cudaGetLastError();
171165 fprintf(stderr,
172166 "Error when trying to find the memory information"
173167 " on the GPU: %s\n", cudaGetErrorString(err2));
174168 }
175169 #if COMPUTE_GPU_MEM_USED
176170 fprintf(stderr,
177171 "Error allocating %llu bytes of device memory (%s)."
178172 " new total bytes allocated: %llu."
179173 " Driver report %llu bytes free and %llu bytes total \n",
180174 (unsigned long long)size, cudaGetErrorString(err), (unsigned long long)_allocated_size,
181175 (unsigned long long)free, (unsigned long long)total);
182176 #else
183177 fprintf(stderr,
184178 "Error allocating %llu bytes of device memory (%s)."
185179 " Driver report %llu bytes free and %llu bytes total \n",
186180 (unsigned long long)size, cudaGetErrorString(err), (unsigned long long)free, (unsigned long long)total);
187181 #endif
188182 }
189183 PyErr_Format(PyExc_MemoryError,
190184 "Error allocating %llu bytes of device memory (%s).",
191185 (unsigned long long)size, cudaGetErrorString(err));
192186 return NULL;
193187 }
194188 }
195189 if (rval != NULL){
196190 // Can it happen that cudaMalloc return cudaSuccess, but return a NULL ptr?
197191 // Could this be what happen if size is 0?
198192 _outstanding_mallocs[0] += 1;
199193
200194 #if COMPUTE_GPU_MEM_USED
201195 _allocated_size += size;
202196 _max_allocated_size = std::max(_max_allocated_size, _allocated_size);
203197 int i = 0;
204198 for(;i<TABLE_SIZE;i++){
205199 if(NULL==_alloc_size_table[i].ptr){
206200 _alloc_size_table[i].ptr=rval;
207201 _alloc_size_table[i].size=size;
208202 break;
209203 }
210204 }
211205 if (i == TABLE_SIZE){
212206 fprintf(stderr,
213207 "When tracking GPU malloc, our table size wasn't big enough."
214208 " So we loose some tracking. Raise the value of TABLE_SIZE in the file cuda_ndarra.cu");
215209 }
216210 #endif
217211 }
218212 //fprintf(stderr,
219213 //"allocated %li bytes of device memory (%s). new total bytes allocated: %d. ptr: %p\n",
220214 //(long)size, cudaGetErrorString(err),_allocated_size,rval);
221215
222216 if(ALLOC_MEMSET){
223217 //We init them to nan to make sure we catch more debug case.
224218 cudaMemset(rval, 0xFF, size);
225219 //printf("MEMSET\n");
226220 }
227221 #if PRINT_FREE_MALLOC
228222 fprintf(stderr, "device malloc %p of size %d\n", rval, size);
229223 #endif
230224 return rval;
231225 }
232226
233227 int device_free(void *ptr)
234228 {
235229 #if PRECHECK_ERROR
236230 cudaThreadSynchronize();
237231 cudaError_t prevError = cudaGetLastError();
238232 if (cudaSuccess != prevError)
239233 {
240234 fprintf(stderr,
241235 "Error existed before calling device_free. %s\n",
242236 cudaGetErrorString(prevError)
243237 );
244238 }
245239 #endif
246240 #if PRINT_FREE_MALLOC
247241 size_t free = 0, total = 0;
248242 cudaError_t err2 = cudaMemGetInfo(&free, &total);
249243 if (err2 != cudaSuccess){
250244 cudaGetLastError();
251245 fprintf(stderr,
252246 "Error when tring to find the memory information"
253247 " on the GPU: %s\n", cudaGetErrorString(err2));
254248 }
255249 #if COMPUTE_GPU_MEM_USED
256250 {
257251 int i = 0;
258252 for(;i<TABLE_SIZE;i++)
259253 if(_alloc_size_table[i].ptr==ptr){
260254 break;
261255 }
262256 assert(i<TABLE_SIZE);
263257 fprintf(stderr, "device_free %p of size %d."
264258 " Driver report %d bytes free and %d bytes total \n",
265259 ptr, _alloc_size_table[i].size, free, total);
266260 }
267261 #else
268262 fprintf(stderr, "device_free %p."
269263 " Driver report %d bytes free and %d bytes total \n",
270264 ptr, free, total);
271265 #endif
272266 #endif
273267
274268 // if there is no gpu context, the call to cudaFree will fail; skip it entirely
275269 if(!g_gpu_context_active) {
276270 return 0;
277271 }
278272
279273 ///@TODO: thejaswi: multi-stream support
280274 if(g_use_cnmem) {
281275 cnmemStatus_t status = cnmemFree(ptr, NULL);
282276 if(status != CNMEM_STATUS_SUCCESS) {
283277 fprintf(stderr, "device_free: cnmemFree call failed! Reason=%s\n",
284278 cnmemGetErrorString(status));
285279 }
286280 }
287281 else {
288282 // We need sync as the Theano's GC could remove intermediate variable that
289283 // are still needed as the gpu kernel are running or in the queue.
290284 CNDA_BEGIN_ALLOW_THREADS
291285 cudaThreadSynchronize();
292286 CNDA_END_ALLOW_THREADS
293287
294288 cudaError_t err = cudaFree(ptr);
295289 if (cudaSuccess != err)
296290 {
297291 // Clear the error flag, cudaFree doesn't do it.
298292 // Currently this returns the same thing as err, but if in future
299293 // it returns something else I still don't see why we should ignore
300294 // it. All we want to do here is reset the flag.
301295 cudaGetLastError();
302296 size_t free = 0, total = 0;
303297 cudaError_t err2 = cudaMemGetInfo(&free, &total);
304298 if (err2 != cudaSuccess){
305299 cudaGetLastError();
306300 fprintf(stderr,
307301 "Error when tring to find the memory information"
308302 " on the GPU: %s\n", cudaGetErrorString(err2));
309303 }
310304 #if COMPUTE_GPU_MEM_USED
311305 {
312306 int i = 0;
313307 for(;i<TABLE_SIZE;i++)
314308 if(_alloc_size_table[i].ptr==ptr){
315309 break;
316310 }
317311 assert(i<TABLE_SIZE);
318312 fprintf(stderr,
319313 "Error freeing device pointer %p (%s) of size %llu. %llu byte already allocated."
320314 " Driver report %llu bytes free and %llu bytes total \n",
321315 ptr, cudaGetErrorString(err),
322316 (unsigned long long)_alloc_size_table[i].size, (unsigned long long)_allocated_size, (unsigned long long)free, (unsigned long long)total);
323317 }
324318 #else
325319 fprintf(stderr,
326320 "Error freeing device pointer %p (%s)."
327321 " Driver report %llu bytes free and %llu bytes total \n",
328322 ptr,
329323 cudaGetErrorString(err), (unsigned long long)free, (unsigned long long)total);
330324 #endif
331325 if (NULL != PyErr_Occurred()){
332326 fprintf(stderr,
333327 "device_free: cudaFree() returned an error, but there is already an"
334328 " Python error set. This happen during the clean up when there is a"
335329 " first error and the CUDA driver is in a so bad state that it don't"
336330 " work anymore. We keep the previous error set to help debugging it.");
337331 return -1;
338332 }
339333 PyErr_Format(PyExc_MemoryError,
340334 "error freeing device pointer %p (%s)",
341335 ptr,
342336 cudaGetErrorString(err));
343337 return -1;
344338 }
345339 }
346340 _outstanding_mallocs[0] -= (ptr != NULL);
347341 #if COMPUTE_GPU_MEM_USED
348342 int i=0;
349343 size_t total_freed = 0;
350344 for(;i<TABLE_SIZE;i++)
351345 if(_alloc_size_table[i].ptr==ptr){
352346 _allocated_size -= _alloc_size_table[i].size;
353347 total_freed += _alloc_size_table[i].size;
354348 _alloc_size_table[i].ptr=0;
355349 _alloc_size_table[i].size=0;
356350
357351 break;
358352 }
359353 //if(i==TABLE_SIZE)
360354 // printf("Unallocated unknow size!\n");
361355 //fprintf(stderr, "freed %li bytes of device memory (%s). %d already allocated, ptr=%p\n", (long)total_freed, cudaGetErrorString(err),_allocated_size,ptr);
362356 #endif
363357 return 0;
364358 }
365359
366360 static PyObject *
367361 outstanding_mallocs(PyObject* self, PyObject * args)
368362 {
369363 return PyInt_FromLong(_outstanding_mallocs[0]);
370364 }
371365
372366
373367 static void *work_mem = NULL;
374368 static size_t work_size = 0;
375369
376370 /*
377371 * Returns a chunk of memory for temporary work inside of an op. You can only
378372 * request a single chunk of memory at a time since it is reused.
379373 */
380374 void *get_work_mem(size_t sz) {
381375 if (sz <= work_size)
382376 return work_mem;
383377 device_free(work_mem);
384378 work_mem = device_malloc(sz);
385379 work_size = sz;
386380 if (work_mem == NULL)
387381 work_size = 0;
388382 return work_mem;
389383 }
390384
391385 /////////////////////////
392386 // Static helper methods
393387 /////////////////////////
394388
395389 static void
396390 CudaNdarray_null_init(CudaNdarray*self)
397391 {
398392 self->base = NULL;
399393 self->nd = -1;
400394 self->host_structure = NULL;
401395 self->data_allocated = 0;
402396 self->dev_structure_fresh = 1;
403397 self->dev_structure = NULL;
404398 self->devdata = NULL;
405399 }
406400
407401 static int
408402 CudaNdarray_uninit(CudaNdarray*self)
409403 {
410404 #if PRINT_FREE_MALLOC
411405 fprintf(stderr, "CudaNdarray_uninit %p\n", self);
412406 #endif
413407 int rval = 0;
414408 if (self->data_allocated) {
415409 assert(self->devdata);
416410 if (device_free(self->devdata))
417411 {
418412 fprintf(stderr,
419413 "CudaNdarray_uninit: error freeing self->devdata. (self=%p, self->devata=%p)\n",
420414 self, self->devdata);
421415 rval = -1;
422416 }
423417 self->devdata = NULL;
424418 self->data_allocated = 0;
425419 }
426420 if (self->dev_structure)
427421 {
428422 if (device_free(self->dev_structure))
429423 {
430424 fprintf(stderr,
431425 "CudaNdarray_uninit: error freeing dev_structure memory %p (self=%p)\n",
432426 self->dev_structure, self);
433427 rval = -1;
434428 }
435429 self->dev_structure = NULL;
436430 }
437431 if (self->host_structure)
438432 {
439433 free(self->host_structure);
440434 self->host_structure = NULL;
441435 }
442436 self->nd = -1;
443437 Py_XDECREF(self->base);
444438 self->base = NULL;
445439 return rval;
446440 }
447441
448442
449443 //make the rightmost coords change fastest
450444 //TODO: why does a downward for-loop not work????
451445 //TODO: use the log2_dims and driver code to remove / and %
452446 //TODO: skip the last division (when d == 0)
453447 #define decl_k_elemwise_unary_rowmajor(name, F) \
454448 __global__ void name (unsigned int numEls, \
455449 unsigned int nd, \
456450 const int * dim, \
457451 const float * a_data, const int * a_str, \
458452 float * z_data, const int * z_str) \
459453 { \
460454 const unsigned int idx = blockIdx.x * blockDim.x + threadIdx.x; \
461455 const unsigned int numThreads = blockDim.x * gridDim.x; \
462456 \
463457 for (unsigned int i = idx; i < numEls; i += numThreads) \
464458 { \
465459 unsigned int ii = i; \
466460 const float * a_i = a_data; \
467461 float * z_i = z_data; \
468462 for (unsigned int _d = 0; _d < nd; ++_d) \
469463 { \
470464 unsigned int d = nd - _d-1; \
471465 int i_d = ii % dim[d]; /* i_d is our position in the d'th dimension */ \
472466 ii = ii / dim[d]; \
473467 a_i += i_d * a_str[d]; /* increment our a and z pointers by i_d elements */ \
474468 z_i += i_d * z_str[d]; \
475469 } \
476470 z_i[0] = F(a_i[0]); \
477471 } \
478472 }
479473
480474 template<typename T> __device__ T unary_copy(T a) { return a; }
481475 decl_k_elemwise_unary_rowmajor(k_elemwise_unary_rowmajor_copy, unary_copy<float>)
482476
483477 template<typename T> __device__ T unary_exp(T a) { return exp(a); }
484478 decl_k_elemwise_unary_rowmajor(k_elemwise_unary_rowmajor_exp, unary_exp<float>)
485479
486480 /////////////////////////////
487481 // Satisfying reqs to be Type
488482 /////////////////////////////
489483
490484 //DON'T use directly(if their is other CudaNdarray that point to it, it will cause problem)! use Py_DECREF() instead
491485 static void
492486 CudaNdarray_dealloc(CudaNdarray* self)
493487 {
494488 if (0) std::cerr << "CudaNdarray dealloc " << self << " " << self->devdata << '\n';
495489 if(Py_REFCNT(self) > 1)
496490 printf("WARNING:CudaNdarray_dealloc called when there is still active reference to it.\n");
497491 CudaNdarray_uninit(self);
498492 Py_TYPE(self)->tp_free((PyObject*)self);
499493 --_outstanding_mallocs[1];
500494 if (0)
501495 {
502496 fprintf(stderr, "device_malloc_counts: (device) %i (obj) %i\n",
503497 _outstanding_mallocs[0],
504498 _outstanding_mallocs[1]);
505499 }
506500 }
507501
508502 static PyObject *
509503 CudaNdarray_new(PyTypeObject *type, PyObject *args, PyObject *kwds)
510504 {
511505 CudaNdarray *self;
512506
513507 self = (CudaNdarray *)type->tp_alloc(type, 0);
514508 if (self != NULL)
515509 {
516510 CudaNdarray_null_init(self);
517511 ++_outstanding_mallocs[1];
518512 }
519513 return (PyObject *)self;
520514 }
521515 static int
522516 CudaNdarray_init(CudaNdarray *self, PyObject *args, PyObject *kwds)
523517 {
524518 PyObject *arr=NULL;
525519
526520 if (! PyArg_ParseTuple(args, "O", &arr))
527521 return -1;
528522 if (! PyArray_Check(arr))
529523 {
530524 PyErr_SetString(PyExc_TypeError, "PyArray arg required");
531525 return -1;
532526 }
533527 int rval = CudaNdarray_CopyFromArray(self, (PyArrayObject*)arr);
534528 return rval;
535529 }
536530 static PyMemberDef CudaNdarray_members[] =
537531 {
538532 /*
539533 {"first", T_OBJECT_EX, offsetof(CudaNdarray, first), 0,
540534 "first name"},
541535 {"last", T_OBJECT_EX, offsetof(CudaNdarray, last), 0,
542536 "last name"},
543537 {"number", T_INT, offsetof(CudaNdarray, number), 0,
544538 "noddy number"},
545539 */
546540 {NULL} /* Sentinel */
547541 };
548542
549543 PyObject * CudaNdarray_CreateArrayObj(CudaNdarray * self, PyObject *args)
550544 {
551545 PyObject * dtype = NULL;
552546 if (args && !PyArg_ParseTuple(args, "|O", &dtype))
553547 return NULL;
554548 if (dtype) {
555549 PyArray_Descr* dtype2;
556550 // PyArray_DescrConverter try to convert anything to a PyArray_Descr.
557551 if(!PyArray_DescrConverter(dtype, &dtype2))
558552 {
559553 PyObject * str = PyObject_Repr(dtype);
560554 PyErr_Format(PyExc_TypeError,
561555 "CudaNdarray dtype parameter not understood: %s",
562556 PyString_AsString(str)
563557 );
564558 Py_CLEAR(str);
565559 return NULL;
566560 }
567561 int typeNum = dtype2->type_num;
568562 Py_DECREF(dtype2);
569563 if (typeNum != NPY_FLOAT32)
570564 {
571565 PyObject * str = PyObject_Repr(dtype);
572566 PyErr_Format(PyExc_TypeError,
573567 "CudaNdarray support only support float32 dtype, provided: %d",
574568 typeNum
575569 );
576570 Py_CLEAR(str);
577571 return NULL;
578572 }
579573 }
580574
581575 int verbose = 0;
582576 if(self->nd>=0 && CudaNdarray_SIZE(self)==0){
583577 npy_intp * npydims = (npy_intp*)malloc(self->nd * sizeof(npy_intp));
584578 assert (npydims);
585579 for (int i = 0; i < self->nd; ++i) npydims[i] = (npy_intp)(CudaNdarray_HOST_DIMS(self)[i]);
586580 PyObject * rval = PyArray_SimpleNew(self->nd, npydims, REAL_TYPENUM);
587581 free(npydims);
588582 if (!rval){
589583 return NULL;
590584 }
591585 assert (PyArray_ITEMSIZE((PyArrayObject *)rval) == sizeof(real));
592586 return rval;
593587 }
594588 if ((self->nd < 0) || (self->devdata == 0))
595589 {
596590 PyErr_SetString(PyExc_ValueError, "can't copy from un-initialized CudaNdarray");
597591 return NULL;
598592 }
599593 CudaNdarray * contiguous_self = NULL;
600594 if (CudaNdarray_is_c_contiguous(self))
601595 {
602596 contiguous_self = self;
603597 Py_INCREF(contiguous_self);
604598 if (verbose) std::cerr << "CreateArrayObj already contiguous" << contiguous_self << '\n';
605599 }
606600 else
607601 {
608602 contiguous_self = (CudaNdarray*)CudaNdarray_Copy(self);
609603 if (verbose) std::cerr << "CreateArrayObj created contiguous" << contiguous_self << '\n';
610604 }
611605 if (!contiguous_self)
612606 {
613607 return NULL;
614608 }
615609
616610 npy_intp * npydims = (npy_intp*)malloc(self->nd * sizeof(npy_intp));
617611 assert (npydims);
618612 for (int i = 0; i < self->nd; ++i)
619613 npydims[i] = (npy_intp)(CudaNdarray_HOST_DIMS(self)[i]);
620614 PyArrayObject * rval = (PyArrayObject *) PyArray_SimpleNew(self->nd,
621615 npydims,
622616 REAL_TYPENUM);
623617 free(npydims);
624618 if (!rval)
625619 {
626620 Py_DECREF(contiguous_self);
627621 return NULL;
628622 }
629623
630624 assert (PyArray_ITEMSIZE(rval) == sizeof(real));
631625
632626 npy_intp rval_size = PyArray_SIZE(rval);
633627 void *rval_data = PyArray_DATA(rval);
634628 cudaError_t err;
635629 CNDA_BEGIN_ALLOW_THREADS;
636630
637631 err = cudaMemcpy(rval_data, contiguous_self->devdata,
638632 rval_size * sizeof(real),
639633 cudaMemcpyDeviceToHost
640634 );
641635 //CNDA_THREAD_SYNC; // unneeded because cudaMemcpy is blocking anyway
642636 CNDA_END_ALLOW_THREADS;
643637
644638 if (cudaSuccess != err)
645639 {
646640 PyErr_Format(PyExc_RuntimeError, "error (%s)copying data to host",
647641 cudaGetErrorString(err));
648642 Py_DECREF(rval);
649643 rval = NULL;
650644 }
651645
652646 Py_DECREF(contiguous_self);
653647 return (PyObject *)rval;
654648 }
655649
656650 // TODO-- we have two functions here, ZEROS and Zeros.
657651 // ZEROS is meant to be called just from C code (you don't need to pass it PyObject * s)
658652 // but this naming is very weird, makes it look like a macro
659653 // we should figure out the correct convention and change to that
660654 PyObject* CudaNdarray_ZEROS(int n, int * dims)
661655 {
662656
663657 size_t total_elements = 1;
664658
665659 for(size_t i=0;i<n;i++){
666660 // Detect overflow on unsigned integer
667661 if (dims[i] != 0 && total_elements > (SIZE_MAX / dims[i])) {
668662 PyErr_Format(PyExc_RuntimeError,
669663 "Can't store in size_t for the bytes requested %llu * %llu",
670664 (unsigned long long)total_elements,
671665 (unsigned long long)dims[i]);
672666 return NULL;
673667 }
674668 total_elements*=dims[i];
675669 }
676670
677671 // total_elements now contains the size of the array, in reals
678672 if (total_elements > (SIZE_MAX / sizeof(real))){
679673 PyErr_Format(PyExc_RuntimeError,
680674 "Can't store in size_t for the bytes requested %llu * 4",
681675 (unsigned long long)total_elements);
682676 return NULL;
683677 }
684678 size_t total_size = total_elements * sizeof(real);
685679
686680 CudaNdarray* rval = (CudaNdarray*)CudaNdarray_New();
687681 if (!rval)
688682 {
689683 PyErr_SetString(PyExc_RuntimeError, "CudaNdarray_ZEROS: call to New failed");
690684 return NULL;
691685 }
692686
693687 if (CudaNdarray_alloc_contiguous(rval, n, dims))
694688 {
695689 PyErr_SetString(PyExc_RuntimeError, "CudaNdarray_ZEROS: allocation failed.");
696690 Py_DECREF(rval);
697691 return NULL;
698692 }
699693
700694 // Fill with zeros
701695 //fprintf(stdout, "Sizeof: %d\n", total_size);
702696 if (cudaSuccess != cudaMemset(rval->devdata, 0, total_size))
703697 {
704698 PyErr_Format(PyExc_MemoryError,
705699 "CudaNdarray_ZEROS: Error memsetting %llu bytes of device memory.",
706700 (unsigned long long)total_size);
707701 Py_DECREF(rval);
708702 return NULL;
709703 }
710704
711705 if (cnda_copy_structure_to_device(rval))
712706 {
713707 PyErr_SetString(PyExc_RuntimeError, "CudaNdarray_ZEROS: syncing structure to device failed");
714708 Py_DECREF(rval);
715709 return NULL;
716710 }
717711 return (PyObject*) rval;
718712 }
719713
720714 // declared as a static method (hence 1st parameter is not used)
721715 // Based on _Copy and _dimshuffle
722716 PyObject* CudaNdarray_Zeros(PyObject* _unused, PyObject* shape)
723717 {
724718 if(!shape)
725719 {
726720 PyErr_SetString(PyExc_TypeError, "CudaNdarray_Zeros: function takes at least 1 argument (0 given)");
727721 return NULL;
728722 }
729723 if(!PySequence_Check(shape))
730724 {
731725 PyErr_SetString(PyExc_TypeError, "shape argument must be a sequence");
732726 return NULL;
733727 }
734728
735729 int shplen = PySequence_Length(shape);
736730
737731 if (shplen == 0)
738732 {
739733 return CudaNdarray_ZEROS(0, NULL);
740734 }
741735
742736 int* newdims = (int *)malloc(sizeof(int) * shplen);
743737
744738 if (!newdims)
745739 {
746740 PyErr_SetString(PyExc_MemoryError,
747741 "CudaNdarray_Zeros: Failed to allocate temporary space");
748742 return NULL;
749743 }
750744
751745 // start from the end to compute strides
752746 for (int i = shplen-1; i >= 0; --i)
753747 {
754748 PyObject* shp_el_obj = PySequence_GetItem(shape, i);
755749 if(shp_el_obj == NULL)
756750 {
757751 // shouldn't happen since we checked length before...
758752 PyErr_SetString(PyExc_RuntimeError, "CudaNdarray_Zeros: Index out of bound in sequence");
759753 free(newdims);
760754 return NULL;
761755 }
762756
763757 int shp_el = PyInt_AsLong(shp_el_obj);
764758 Py_DECREF(shp_el_obj);
765759
766760 if (shp_el < 0)
767761 {
768762 PyErr_SetString(PyExc_ValueError, "CudaNdarray_Zeros: shape must contain only non-negative values for size of a dimension");
769763 free(newdims);
770764 return NULL;
771765 }
772766
773767 newdims[i] = shp_el;
774768 }
775769
776770 PyObject* rval = CudaNdarray_ZEROS(shplen,newdims);
777771
778772 free(newdims);
779773
780774 return (PyObject*)rval;
781775 }
782776
783777
784778
785779
786780
787781 PyObject * CudaNdarray_Copy(const CudaNdarray * self)
788782 {
789783 PyObject * rval = CudaNdarray_New();
790784 if ((!rval) || (-1 == self->nd))
791785 {
792786 return rval;
793787 }
794788 if (CudaNdarray_alloc_contiguous((CudaNdarray*)rval, self->nd, CudaNdarray_HOST_DIMS(self)))
795789 {
796790 Py_DECREF(rval);
797791 return NULL;
798792 }
799793 if (CudaNdarray_CopyFromCudaNdarray((CudaNdarray*)rval, self))
800794 {
801795 Py_DECREF(rval);
802796 return NULL;
803797 }
804798 return rval;
805799 }
806800 PyObject * CudaNdarray_DeepCopy(CudaNdarray * self, PyObject * memo)
807801 {
808802 assert(PyDict_Check(memo));
809803 PyObject * selfkey = PyInt_FromLong((long)self);
810804 assert(selfkey);
811805 if (PyDict_Contains(memo, selfkey))
812806 {
813807 PyObject * rval = PyDict_GetItem(memo, selfkey);
814808 Py_DECREF(selfkey);
815809 Py_XINCREF(rval);
816810 return rval;
817811 }
818812 else
819813 {
820814 PyObject * rval = CudaNdarray_Copy(self);
821815 if (0) std::cerr << "DeepCopy created " << rval << " devdata " << ((CudaNdarray*)rval)->devdata << "\n";
822816 if (NULL == rval)
823817 {
824818 Py_DECREF(selfkey);
825819 return NULL;
826820 }
827821 if (PyDict_SetItem(memo, selfkey, rval))
828822 {
829823 Py_DECREF(rval);
830824 Py_DECREF(selfkey);
831825 return NULL;
832826 }
833827 Py_DECREF(selfkey);
834828 return rval;
835829 }
836830 }
837831 PyObject * CudaNdarray_ReduceSum(CudaNdarray * self, PyObject * py_reduce_mask)
838832 {
839833 if (!PySequence_Check(py_reduce_mask))
840834 {
841835 PyErr_SetString(PyExc_TypeError, "reduce_mask must be sequence of ints");
842836 return NULL;
843837 }
844838 int len = PySequence_Length(py_reduce_mask);
845839 if (len != self->nd)
846840 {
847841 PyErr_SetString(PyExc_TypeError, "length of reduce_mask must match self->nd");
848842 return NULL;
849843 }
850844 CudaNdarray * self_sum = (CudaNdarray*)CudaNdarray_New();
851845 if (!self_sum)
852846 {
853847 return NULL;
854848 }
855849 //TODO: allocate a fixed size dimshuffle_pattern_cache on the stack,
856850 // and use it if it is big enough.
857851 int * dimshuffle_pattern = (int*)malloc(len * 2 * sizeof(int));
858852 int * sum_dims = dimshuffle_pattern + len;
859853 int n_remaining_dims = 0;
860854 if (!dimshuffle_pattern)
861855 {
862856 Py_DECREF(self_sum);
863857 PyErr_SetString(PyExc_MemoryError, "failed to alloc internal storage");
864858 return NULL;
865859 }
866860 for (int i = 0; i < len; ++i)
867861 {
868862 PyObject *o_i = PySequence_GetItem(py_reduce_mask, i);
869863 int o_i_int = PyInt_AsLong(o_i);
870864 Py_XDECREF(o_i);
871865 if (PyErr_Occurred())
872866 {
873867 Py_DECREF(self_sum);
874868 free(dimshuffle_pattern);
875869 return NULL;
876870 }
877871 if (o_i_int) // this is a dimension over which we are reducing
878872 {
879873 sum_dims[i] = 1;
880874 }
881875 else
882876 {
883877 sum_dims[i] = CudaNdarray_HOST_DIMS(self)[i];
884878 dimshuffle_pattern[n_remaining_dims++] = i;
885879 }
886880 }
887881 if (0 || CudaNdarray_alloc_contiguous(self_sum, len, sum_dims)
888882 || CudaNdarray_reduce_sum(self_sum, self)
889883 || CudaNdarray_dimshuffle(self_sum, n_remaining_dims, dimshuffle_pattern))
890884 {
891885 Py_DECREF(self_sum);
892886 free(dimshuffle_pattern);
893887 return NULL;
894888 }
895889 free(dimshuffle_pattern);
896890 return (PyObject*)self_sum;
897891 }
898892
899893 // Reshape self to the new shape gived by the tuple shape.
900894 //
901895 // If self is c contiguous, it return a view. Otherwise it always do a copy.
902896 // TODO: make it return a view when the strides allow it even if it is not
903897 // c contiguous
904898 PyObject * CudaNdarray_Reshape(CudaNdarray * self, PyObject * shape)
905899 {
906900 if(!CudaNdarray_is_c_contiguous(self))
907901 {
908902 // allocate new space
909903 //TODO: test to see if we can re-use old one and take a new param to
910904 // use this
911905 CudaNdarray* rval = (CudaNdarray*) CudaNdarray_Copy(self);
912906 if (!rval)
913907 {
914908 return NULL;
915909 }
916910
917911 CudaNdarray* ret = (CudaNdarray*) CudaNdarray_Reshape(rval, shape);
918912 Py_XDECREF(rval);
919913 return (PyObject*)ret;
920914 }
921915
922916 // check shape tuple
923917 unsigned int rval_nd;
924918 unsigned int * rval_dims;
925919 size_t rval_size = 1;
926920
927921 if (PyTuple_Check(shape)){
928922 // copy shape to integer array
929923 rval_nd = PyTuple_Size(shape);
930924 }else if (PyInt_Check(shape)){
931925 rval_nd = 1;
932926 }else{
933927 PyErr_SetString(PyExc_TypeError, "shape must be tuple of integers or an integer");
934928 return NULL;
935929 }
936930 rval_dims = (unsigned int*)malloc(rval_nd * sizeof(int));
937931
938932 if(PyTuple_Check(shape)){
939933 for (int i = 0; i < rval_nd; ++i)
940934 {
941935 rval_dims[i] = PyInt_AsLong(PyTuple_GetItem(shape, i)); //GetItem returns borrowed reference
942936 if (PyErr_Occurred()) //error in AsLong
943937 {
944938 free(rval_dims);
945939 return NULL;
946940 }
947941 if(rval_dims[i]<0){
948942 PyErr_Format(PyExc_ValueError, "Reshape has invalid dimension %i (must be >=0)",rval_dims[i]);
949943 free(rval_dims);
950944 return NULL;
951945 }
952946 rval_size = rval_size * rval_dims[i];
953947 }
954948 }else{
955949 rval_size = PyInt_AsLong(shape);
956950 rval_dims[0] = rval_size;
957951 }
958952 // calculate new size, assert same as old size
959953 if (rval_size != CudaNdarray_SIZE(self))
960954 {
961955 PyErr_Format(PyExc_ValueError, "size must remain unchanged, changed from %lld to %lld", CudaNdarray_SIZE(self), rval_size);
962956 free(rval_dims);
963957 return NULL;
964958 }
965959 if (rval_size==0)
966960 {
967961 PyObject * rval = CudaNdarray_NewDims(rval_nd, rval_dims);
968962 free(rval_dims);
969963 return rval;
970964 }
971965
972966 //return a view, not a copy
973967 //we can do this as we checked self is c_contiguous
974968 CudaNdarray * rval = (CudaNdarray * )CudaNdarray_New(rval_nd);
975969
976970 if (!rval || 0 != rval->data_allocated
977971 ||CudaNdarray_set_device_data(rval, CudaNdarray_DEV_DATA(self), self))
978972 {
979973 Py_XDECREF(rval);
980974 free(rval_dims);
981975 return NULL;
982976 }
983977 //set dim and stride
984978 int size = 1;
985979 for (int i = rval_nd-1; i >= 0; --i)
986980 {
987981 CudaNdarray_set_stride(rval, i, (rval_dims[i] == 1) ? 0 : size);
988982 CudaNdarray_set_dim(rval, i, rval_dims[i]);
989983 size = size * rval_dims[i];
990984 }
991985 free(rval_dims);
992986 return (PyObject*)rval;
993987 }
994988
995989 PyObject * CudaNdarray_View(const CudaNdarray * self)
996990 {
997991 CudaNdarray * rval = (CudaNdarray*)CudaNdarray_New(self->nd);
998992 if (!rval || CudaNdarray_set_device_data(rval, CudaNdarray_DEV_DATA(self), self))
999993 {
1000994 Py_XDECREF(rval);
1001995 rval = NULL;
1002996 }
1003997 else
1004998 {
1005999 for (int i = 0; i < self->nd; ++i)
10061000 {
10071001 CudaNdarray_set_dim(rval, i, CudaNdarray_HOST_DIMS(self)[i]);
10081002 CudaNdarray_set_stride(rval, i, CudaNdarray_HOST_STRIDES(self)[i]);
10091003 }
10101004 }
10111005 return (PyObject*)rval;
10121006 }
10131007
10141008 /*
10151009 * d0,... are the output dims
10161010 * indices are a list of index to operate on
10171011 * They are int32 viewed as float32.
10181012 * a is the output
10191013 * b is the input
10201014 * dB0, the source leading dimensions size
10211015 */
10221016 template <int operator_num>
10231017 __global__ void k_take_3(const int d0, const int d1, const int d2,
10241018 const npy_int64* indices,
10251019 float* a,
10261020 const int sA0, const int sA1, const int sA2,
10271021 const float* b, const int dB0,
10281022 const int sB0, const int sB1, const int sB2,
10291023 int* err){
10301024 for (int i0 = blockIdx.x; i0 < d0; i0 += gridDim.x){
10311025 npy_int64 idx = indices[i0];
10321026 if (idx<0)
10331027 idx += dB0; // To allow negative indexing.
10341028 if ((idx < 0) || (idx >= dB0)){
10351029 // Any value other the 0 probably work. But to be more safe, I want
10361030 // to change all bits to prevent problem with concurrent write that
10371031 // could cross cache line. But this should not happen with the
10381032 // current code and driver.
10391033 *err = 0xFFFF;
10401034 continue;
10411035 }
10421036 for (int i1 = threadIdx.x; i1 < d1; i1 += blockDim.x){
10431037 for (int i2 = threadIdx.y; i2 < d2; i2 += blockDim.y){
10441038 int a_idx = i0*sA0 + i1*sA1 + i2*sA2;
10451039 int b_idx = idx*sB0 + i1*sB1 + i2*sB2;
10461040 a[a_idx] = b[b_idx];
10471041 }
10481042 }
10491043 }
10501044 }
10511045
10521046 // We try to be similar to the PyArray_TakeFrom function
10531047 //http://docs.scipy.org/doc/numpy/reference/c-api.array.html
10541048 //TODO: support other clip mode then raise(clip, wrap)
10551049 //self is the input that we copy data from.
10561050 //The indices that we receive MUST be an CudaNdarray(float32)
10571051 // that is in fact a view to int64 indices
10581052 PyObject*
10591053 CudaNdarray_TakeFrom(CudaNdarray * self, PyObject *args){
10601054 int verbose = 0;
10611055 PyObject * indices_obj = NULL;
10621056 //int axis; Default None, that mean the flattened array.
10631057 PyObject * axis_obj = Py_None;
10641058 PyObject * out_obj = Py_None;
10651059 PyObject * clipmode_obj = NULL;
10661060 int max_threads = 1; // max threads per blocks
10671061
10681062 if (! PyArg_ParseTuple(args, "O|OOOi", &indices_obj, &axis_obj,
10691063 &out_obj, &clipmode_obj, &max_threads))
10701064 return NULL;
10711065
10721066 //Check argument indices
10731067 //TODO: if not a numpy.ndarray, convert to numpy.ndarray
10741068 //TODO: If a CudaNdarray, accept it and suppose the data is int32? is float32 number of int?
10751069 //TODO: Support ndarray of other dtype then int32
10761070 //TODO: support list of indices that are not c_contiguous
10771071 CudaNdarray * indices = NULL;
10781072 if (CudaNdarray_Check(indices_obj)) {
10791073 if (verbose) printf("cudandarray indices\n");
10801074 indices = (CudaNdarray*) indices_obj;
10811075 Py_INCREF(indices);
10821076 } else if (PyArray_Check(indices_obj)) {
10831077 if (verbose) printf("ndarray indices\n");
10841078 if (PyArray_TYPE((PyArrayObject *)indices_obj) != NPY_INT64) {
10851079 PyErr_SetString(PyExc_TypeError,
10861080 "CudaNdarray_TakeFrom: need a ndarray for indices"
10871081 " with dtype int64");
10881082 return NULL;
10891083 }
10901084 if (PyArray_NDIM(((PyArrayObject*)indices_obj)) != 1) {
10911085 PyErr_SetString(PyExc_TypeError,
10921086 "CudaNdarray_TakeFrom: need a CudaNdarray of"
10931087 " indices with only 1 dimensions");
10941088 return NULL;
10951089 }
10961090 // We need indices_obj to be contiguous, in order to take a view
10971091 // with a different dtype.
10981092 if (!PyArray_IS_C_CONTIGUOUS((PyArrayObject*) indices_obj)) {
10991093 PyObject* indices_obj_contig = PyArray_NewCopy((PyArrayObject*) indices_obj, NPY_CORDER);
11001094 if (!indices_obj_contig)
11011095 return NULL;
11021096 indices_obj = indices_obj_contig;
11031097 } else {
11041098 // Keep the refcount consistent
11051099 Py_INCREF(indices_obj);
11061100 }
11071101 PyArray_Descr* float32_descr = PyArray_DescrFromType(NPY_FLOAT32);
11081102 PyObject * indices_float32 = NULL;
11091103 indices_float32 = PyArray_View((PyArrayObject*)indices_obj,
11101104 float32_descr, NULL);
11111105 if (verbose) printf("ndarray indices\n");
11121106 if (!indices_float32) {
11131107 Py_DECREF(indices_obj);
11141108 return NULL;
11151109 }
11161110
11171111 indices = (CudaNdarray*) CudaNdarray_New();
11181112 if (verbose) printf("\nndarray after new\n");
11191113 if (! indices){
11201114 Py_DECREF(indices_obj);
11211115 Py_DECREF(indices_float32);
11221116 return NULL;
11231117 }
11241118 if (CudaNdarray_CopyFromArray(indices,
11251119 (PyArrayObject *)indices_float32)){
11261120 Py_DECREF(indices_obj);
11271121 Py_DECREF(indices_float32);
11281122 return NULL;
11291123 }
11301124 Py_DECREF(indices_obj);
11311125 Py_DECREF(indices_float32);
11321126 } else {
11331127 PyObject* py_s = PyObject_Str(indices_obj);
11341128 const char* s = PyString_AsString(py_s);
11351129 Py_DECREF(py_s);
11361130 PyErr_Format(PyExc_TypeError,
11371131 "CudaNdarray_TakeFrom: need an ndarray of int64 or a"
11381132 " CudaNdarray(float32) that is a view from int64 data"
11391133 " for indices. Got %s", s);
11401134 return NULL;
11411135 }
11421136
11431137 if (verbose) {
11441138 printf("indices used on the gpu\n");
11451139 fprint_CudaNdarray(stdout, indices);
11461140 PyObject * used_indices = CudaNdarray_CreateArrayObj(indices);
11471141 PyObject_Print(used_indices, stdout, 0);
11481142 Py_DECREF(used_indices);
11491143 }
11501144 if (verbose) printf("after print of object\n");
11511145 if(!CudaNdarray_is_c_contiguous(indices) != 0) {
11521146 PyErr_SetString(PyExc_NotImplementedError,
11531147 "CudaNdarray_TakeFrom: The indices must be contiguous in memory.");
11541148 Py_DECREF(indices);
11551149 return NULL;
11561150 }
11571151 int nb_indices = CudaNdarray_SIZE((CudaNdarray *)indices) / 2;// int64 are 8 bytes, float32 are 4 bytes
11581152
11591153 //Check argument axis
11601154 //TODO: implement the default and other axis
11611155 long axis = PyInt_AsLong(axis_obj);
11621156
11631157 if (axis != 0) {
11641158 PyErr_Format(PyExc_NotImplementedError,
11651159 "CudaNdarray_TakeFrom: only axis=0 is currently supported."
11661160 " Got %ld.", axis);
11671161 Py_DECREF(indices);
11681162 return NULL;
11691163 }
11701164
11711165 //Check argument out_obj
11721166 CudaNdarray * out = NULL;
11731167 if (out_obj && CudaNdarray_Check(out_obj))
11741168 out = (CudaNdarray*) out_obj;
11751169 if (out && (out->nd != self->nd ||
11761170 CudaNdarray_HOST_DIMS(out)[0] != nb_indices))
11771171 out = NULL;
11781172 int * dims = (int *)malloc(sizeof(int) * self->nd);
11791173 dims[0] = nb_indices;
11801174
11811175 for (int i=1 ; i<self->nd ; i++) {
11821176 dims[i] = CudaNdarray_HOST_DIMS(self)[i];
11831177 if (out && CudaNdarray_HOST_DIMS(out)[i] != dims[i]) {
11841178 out = NULL;
11851179 }
11861180 }
11871181 if (!out) {
11881182 out = (CudaNdarray*)CudaNdarray_New();
11891183 if (!out){
11901184 Py_DECREF(indices);
11911185 free(dims);
11921186 return NULL;
11931187 }
11941188 if (CudaNdarray_alloc_contiguous(out, self->nd, dims)) {
11951189 Py_DECREF(out);
11961190 Py_DECREF(indices);
11971191 free(dims);
11981192 return NULL;
11991193 }
12001194 }else {
12011195 Py_INCREF(out);
12021196 }
12031197
12041198 //Check argument clipmode
12051199 if (clipmode_obj) {
12061200 char * clipmode = PyString_AsString(clipmode_obj);
12071201 if (! clipmode){
12081202 Py_DECREF(indices);
12091203 Py_DECREF(out);
12101204 free(dims);
12111205 return NULL;
12121206 }
12131207 if (strcmp(clipmode, "raise") != 0) {
12141208 PyErr_Format(PyExc_NotImplementedError,
12151209 "CudaNdarray_TakeFrom: only the raise mode is currently supported. Got '%s'",
12161210 clipmode);
12171211 Py_DECREF(indices);
12181212 Py_DECREF(out);
12191213 free(dims);
12201214 return NULL;
12211215 }
12221216 }
12231217 void (*k3)(const int, const int, const int,
12241218 const npy_int64*,
12251219 float*, const int, const int, const int,
12261220 const float*, const int,
12271221 const int, const int, const int,
12281222 int*);
12291223 k3 = k_take_3<CPY>;
12301224
12311225 // Create the memory place that will store the error information.
12321226 if(init_err_var() != 0) return NULL;
12331227
12341228 dim3 n_blocks(std::min(CudaNdarray_HOST_DIMS(out)[0],65535),1,1);
12351229 if(CudaNdarray_HOST_DIMS(out)[0] == 0){
12361230 // We take 0 elements, so no need for the rest of the code.
12371231 // This speed up that case AND fix crash otherwise.
12381232 free(dims);
12391233 Py_DECREF(indices);
12401234 return (PyObject *)out;
12411235 }
12421236
12431237 switch (self->nd) {
12441238 case 1:
12451239 {
12461240 dim3 n_threads(1, 1, 1);
12471241 if (verbose)
12481242 printf("cudaGetLastError=%d, nd=%d"
12491243 " kernel config: (n_blocks.x=%d, n_blocks.y=%d,"
12501244 " n_threads.x=%i, n_threads.y=%i)\n",
12511245 cudaGetLastError(), self->nd,
12521246 n_blocks.x, n_blocks.y, n_threads.x, n_threads.y);
12531247 k3<<<n_blocks, n_threads>>>(
12541248 dims[0],
12551249 1,
12561250 1,
12571251 (npy_int64*) CudaNdarray_DEV_DATA(indices),
12581252 CudaNdarray_DEV_DATA(out),
12591253 CudaNdarray_HOST_STRIDES(out)[0], //strides
12601254 1,
12611255 1,
12621256 CudaNdarray_DEV_DATA(self),
12631257 CudaNdarray_HOST_DIMS(self)[0], //For indices check
12641258 CudaNdarray_HOST_STRIDES(self)[0], //strides
12651259 1,
12661260 1,
12671261 err_var);
12681262 }
12691263 break;
12701264 case 2:
12711265 {
12721266 dim3 n_threads(std::min(CudaNdarray_HOST_DIMS(out)[1], max_threads), 1, 1);
12731267
12741268 if (verbose)
12751269 printf("cudaGetLastError=%d, nd=%d"
12761270 " kernel config: (n_blocks.x=%d, n_blocks.y=%d,"
12771271 " n_threads.x=%i, n_threads.y=%i)\n",
12781272 cudaGetLastError(), self->nd,
12791273 n_blocks.x, n_blocks.y, n_threads.x, n_threads.y);
12801274
12811275 k3<<<n_blocks, n_threads>>>(
12821276 dims[0], //dimensions
12831277 dims[1],
12841278 1,
12851279 (npy_int64*) CudaNdarray_DEV_DATA(indices),
12861280 CudaNdarray_DEV_DATA(out),
12871281 CudaNdarray_HOST_STRIDES(out)[0], //strides
12881282 CudaNdarray_HOST_STRIDES(out)[1],
12891283 1,
12901284 CudaNdarray_DEV_DATA(self),
12911285 CudaNdarray_HOST_DIMS(self)[0], //For indices check
12921286 CudaNdarray_HOST_STRIDES(self)[0], //strides
12931287 CudaNdarray_HOST_STRIDES(self)[1],
12941288 1,
12951289 err_var);
12961290 }
12971291 break;
12981292 case 3:
12991293 {
13001294 int ty = std::min(CudaNdarray_HOST_DIMS(out)[2], max_threads);
13011295 int tx = std::min(CudaNdarray_HOST_DIMS(out)[1], max_threads / ty);
13021296 dim3 n_threads(tx, ty, 1);
13031297 if (verbose)
13041298 printf("cudaGetLastError=%d, nd=%d"
13051299 " kernel config: (n_blocks.x=%d, n_blocks.y=%d,"
13061300 " n_threads.x=%i, n_threads.y=%i)\n",
13071301 cudaGetLastError(), self->nd,
13081302 n_blocks.x, n_blocks.y, n_threads.x, n_threads.y);
13091303 k3<<<n_blocks, n_threads>>>(
13101304 dims[0], //dimensions
13111305 dims[1],
13121306 dims[2],
13131307 (npy_int64*) CudaNdarray_DEV_DATA(indices),
13141308 CudaNdarray_DEV_DATA(out),
13151309 CudaNdarray_HOST_STRIDES(out)[0], //strides
13161310 CudaNdarray_HOST_STRIDES(out)[1],
13171311 CudaNdarray_HOST_STRIDES(out)[2],
13181312 CudaNdarray_DEV_DATA(self),
13191313 CudaNdarray_HOST_DIMS(self)[0], //For indices check
13201314 CudaNdarray_HOST_STRIDES(self)[0], //strides
13211315 CudaNdarray_HOST_STRIDES(self)[1],
13221316 CudaNdarray_HOST_STRIDES(self)[2],
13231317 err_var);
13241318 }
13251319 break;
13261320 default:
13271321 PyErr_SetString(PyExc_NotImplementedError,
13281322 "CudaNdarray_TakeFrom: only input with 1, 2 or 3"
13291323 " dimensions are currently supported");
13301324
13311325 }
13321326 free(dims);
13331327 CNDA_THREAD_SYNC;
13341328 cudaError_t err = cudaGetLastError();
13351329 if (cudaSuccess != err) {
13361330 PyErr_Format(PyExc_RuntimeError,
13371331 "Cuda error: %s: %s.\n",
13381332 "CudaNdarray_TakeFrom",
13391333 cudaGetErrorString(err));
13401334 Py_DECREF(indices);
13411335 Py_DECREF(out);
13421336 return NULL;
13431337 }
13441338
13451339 int index_err = check_err_var();
13461340 Py_DECREF(indices);
13471341 if (index_err != 0) {
13481342 Py_DECREF(out);
13491343 return NULL;
13501344 }
13511345
13521346 if (verbose) printf("TAKE SUCCEDED\n");
13531347 return (PyObject *)out;
13541348 }
13551349
13561350
13571351 PyObject * CudaNdarray_SetStride(CudaNdarray * self, PyObject *args)
13581352 {
13591353 int pos, stride;
13601354 if (! PyArg_ParseTuple(args, "ii", &pos, &stride))
13611355 return NULL;
13621356 if ((pos < 0) || (pos >= self->nd))
13631357 {
13641358 PyErr_Format(PyExc_ValueError, "position argument out of legal range [0, %i)", self->nd);
13651359 return NULL;
13661360 }
13671361 CudaNdarray_set_stride(self, pos, stride);
13681362 if (cnda_copy_structure_to_device(self))
13691363 {
13701364 return NULL;
13711365 }
13721366 Py_INCREF(Py_None);
13731367 return Py_None;
13741368 }
13751369 PyObject * CudaNdarray_SetShapeI(CudaNdarray * self, PyObject *args)
13761370 {
13771371 int pos, dim;
13781372 if (! PyArg_ParseTuple(args, "ii", &pos, &dim))
13791373 return NULL;
13801374 if ((pos < 0) || (pos >= self->nd))
13811375 {
13821376 PyErr_Format(PyExc_ValueError, "position argument out of legal range [0, %i)", self->nd);
13831377 return NULL;
13841378 }
13851379 CudaNdarray_set_dim(self, pos, dim);
13861380 if (cnda_copy_structure_to_device(self))
13871381 {
13881382 return NULL;
13891383 }
13901384 Py_INCREF(Py_None);
13911385 return Py_None;
13921386 }
13931387
13941388 static PyObject *
13951389 CudaNdarray_exp(CudaNdarray* self)
13961390 {
13971391 CudaNdarray * rval = (CudaNdarray *)CudaNdarray_New();
13981392 if ((NULL == rval) || CudaNdarray_alloc_contiguous(rval, self->nd, CudaNdarray_HOST_DIMS(self)))
13991393 {
14001394 Py_XDECREF(rval);
14011395 return NULL;
14021396 }
14031397 unsigned int size = 1;
14041398 for (int i = 0; i < self->nd; i++)
14051399 {
14061400 size *= (unsigned int) CudaNdarray_HOST_DIMS(self)[i];
14071401 }
14081402 unsigned int threads_per_block = std::min(size, (unsigned int)NUM_VECTOR_OP_THREADS_PER_BLOCK);
14091403 unsigned int n_blocks = std::min(ceil_intdiv(size,threads_per_block), (unsigned int)NUM_VECTOR_OP_BLOCKS);
14101404 k_elemwise_unary_rowmajor_exp<<<n_blocks,threads_per_block>>>(size, self->nd, CudaNdarray_DEV_DIMS(self),
14111405 CudaNdarray_DEV_DATA(self), CudaNdarray_DEV_STRIDES(self),
14121406 CudaNdarray_DEV_DATA(rval), CudaNdarray_DEV_STRIDES(rval));
14131407
14141408 //TODO: don't do this right away, do it when we need the result
14151409 CNDA_THREAD_SYNC;
14161410 cudaError_t err = cudaGetLastError();
14171411 if( cudaSuccess != err)
14181412 {
14191413 Py_DECREF(rval);
14201414 PyErr_Format(PyExc_RuntimeError, "Cuda error: %s: %s.\n", "kExp", cudaGetErrorString(err));
14211415 return NULL;
14221416 }
14231417
14241418 return (PyObject*)rval;
14251419 }
14261420
14271421 static PyMethodDef CudaNdarray_methods[] =
14281422 {
14291423 {"__array__",
14301424 (PyCFunction)CudaNdarray_CreateArrayObj, METH_VARARGS,
14311425 "Copy from the device to a numpy ndarray"},
14321426 {"__copy__",
14331427 (PyCFunction)CudaNdarray_View, METH_NOARGS,
14341428 "Create a shallow copy of this object. used by module copy"},
14351429 {"__deepcopy__",
14361430 (PyCFunction)CudaNdarray_DeepCopy, METH_O,
14371431 "Create a copy of this object"},
14381432 {"zeros",
14391433 (PyCFunction)CudaNdarray_Zeros, METH_STATIC | METH_O,
14401434 "Create a new CudaNdarray with specified shape, filled with zeros."},
14411435 {"copy",
14421436 (PyCFunction)CudaNdarray_Copy, METH_NOARGS,
14431437 "Create a copy of this object"},
14441438 {"is_c_contiguous",
14451439 (PyCFunction)CudaNdarray_IS_C_Contiguous, METH_NOARGS,
14461440 "Return True is the object is c contiguous. False otherwise."},
14471441 {"reduce_sum",
14481442 (PyCFunction)CudaNdarray_ReduceSum, METH_O,
14491443 "Reduce over the given dimensions by summation"},
14501444 {"exp",
14511445 (PyCFunction)CudaNdarray_exp, METH_NOARGS,
14521446 "Return the exponential of all elements"},
14531447 {"reshape",
14541448 (PyCFunction)CudaNdarray_Reshape, METH_O,
14551449 "Return a reshaped view (or copy) of this ndarray\n\
14561450 The required argument is a tuple of integers specifying the shape of the new ndarray."},
14571451 {"view",
14581452 (PyCFunction)CudaNdarray_View, METH_NOARGS,
14591453 "Return an alias of this ndarray"},
14601454 {"_set_stride",
14611455 (PyCFunction)CudaNdarray_SetStride, METH_VARARGS,
14621456 "For integer arguments (i, s), set the 'i'th stride to 's'"},
14631457 {"take",
14641458 (PyCFunction)CudaNdarray_TakeFrom, METH_VARARGS,
14651459 "Equivalent of numpy.take"},
14661460 {"_set_shape_i",
14671461 (PyCFunction)CudaNdarray_SetShapeI, METH_VARARGS,
14681462 "For integer arguments (i, s), set the 'i'th shape to 's'"},
14691463 {NULL, NULL, NULL, NULL} /* Sentinel */
14701464 };
14711465
14721466
14731467 ////////////////////
14741468 // Number protocol
14751469 ////////////////////
14761470
14771471 __global__ void kAdd_contiguous(float* a, float* b, float* dest, unsigned int numEls) {
14781472 const unsigned int idx = blockIdx.x * blockDim.x + threadIdx.x;
14791473 const unsigned int numThreads = blockDim.x * gridDim.x;
14801474
14811475 for (unsigned int i = idx; i < numEls; i += numThreads) {
14821476 dest[i] = a[i] + b[i];
14831477 }
14841478 }
14851479
14861480 // Will be called by __add__ in Python
14871481 static PyObject *
14881482 CudaNdarray_add(PyObject* py_self, PyObject * py_other)
14891483 {
14901484 if (! CudaNdarray_Check(py_self)) {
14911485 PyErr_SetString(PyExc_TypeError, "need a CudaNdarray on left");
14921486 return NULL;
14931487 }
14941488 if (! CudaNdarray_Check(py_other)) {
14951489 PyErr_SetString(PyExc_TypeError, "need a CudaNdarray on right");
14961490 return NULL;
14971491 }
14981492 CudaNdarray * self = (CudaNdarray *)py_self;
14991493 CudaNdarray * other = (CudaNdarray *)py_other;
15001494 if(!CudaNdarray_is_c_contiguous(self) || !CudaNdarray_is_c_contiguous(other)){
15011495 PyErr_SetString(PyExc_TypeError, "We have implementet only the c_contiguous version for now.");
15021496 return NULL;
15031497 }
15041498
15051499 //standard elemwise size checks
15061500 if (self->nd != other->nd)
15071501 {
15081502 PyErr_SetString(PyExc_TypeError, "CudaNdarray_add: need same number of dims");
15091503 return NULL;
15101504 }
15111505 //standard elemwise dim checks
15121506 unsigned int size = 1;
15131507 for (int i = 0; i< self->nd; ++i)
15141508 {
15151509 if (CudaNdarray_HOST_DIMS(self)[i] != CudaNdarray_HOST_DIMS(other)[i])
15161510 {
15171511 PyErr_SetString(PyExc_TypeError, "need same dimensions");
15181512 return NULL;
15191513 }
15201514 size *= (unsigned int) CudaNdarray_HOST_DIMS(self)[i];
15211515 }
15221516 CudaNdarray * rval = (CudaNdarray *)CudaNdarray_New();
15231517 if (!rval || CudaNdarray_alloc_contiguous(rval, self->nd, CudaNdarray_HOST_DIMS(self)))
15241518 {
15251519 Py_XDECREF(rval);
15261520 return NULL;
15271521 }
15281522
15291523 if(CudaNdarray_SIZE((CudaNdarray *)py_self)==0 && CudaNdarray_SIZE((CudaNdarray *)py_other)==0){
15301524 return (PyObject *) rval;
15311525 }
15321526
15331527 int threads_per_block = std::min(size, (unsigned int)NUM_VECTOR_OP_THREADS_PER_BLOCK);
15341528 int n_blocks = std::min(ceil_intdiv(size,(unsigned int)threads_per_block), (unsigned int)NUM_VECTOR_OP_BLOCKS);
15351529 kAdd_contiguous<<<n_blocks,threads_per_block>>>(
15361530 self->devdata, other->devdata, rval->devdata, size);
15371531 CNDA_THREAD_SYNC;
15381532 cudaError_t err = cudaGetLastError();
15391533 if( cudaSuccess != err)
15401534 {
15411535 PyErr_Format(PyExc_RuntimeError, "Cuda error: %s: %s.\n", "kAdd", cudaGetErrorString(err));
15421536 Py_DECREF(rval);
15431537 return NULL;
15441538 }
15451539 return (PyObject *) rval;
15461540 }
15471541
15481542 template <int operator_num>
15491543 __global__ void k_ielem_3(const int d0, const int d1, const int d2,
15501544 float* a, const int sA0, const int sA1, const int sA2,
15511545 const float* b, const int sB0, const int sB1, const int sB2){
15521546 for (int i0 = blockIdx.x; i0 < d0; i0 += gridDim.x){
15531547 for (int i1 = blockIdx.y; i1 < d1; i1 += gridDim.y){
15541548 for (int i2 = threadIdx.x; i2 < d2; i2 += blockDim.x){
15551549 switch (operator_num)
15561550 {
15571551 case IADD:
15581552 a[i0*sA0 + i1*sA1 + i2*sA2] += b[i0*sB0 + i1*sB1 + i2*sB2];
15591553 break;
15601554 case IDIV:
15611555 a[i0*sA0 + i1*sA1 + i2*sA2] /= b[i0*sB0 + i1*sB1 + i2*sB2];
15621556 break;
15631557 case CPY:
15641558 a[i0*sA0 + i1*sA1 + i2*sA2] = b[i0*sB0 + i1*sB1 + i2*sB2];
15651559 break;
15661560 }
15671561 }
15681562 }
15691563 }
15701564 }
15711565
15721566 template <int operator_num>
15731567 __global__ void k_ielem_4(const int d0, const int d1, const int d2, const int d3,
15741568 float* a, const int sA0, const int sA1,
15751569 const int sA2, const int sA3,
15761570 const float* b, const int sB0, const int sB1,
15771571 const int sB2, const int sB3){
15781572 for (int i0 = blockIdx.x; i0 < d0; i0 += gridDim.x){
15791573 for (int i1 = blockIdx.y; i1 < d1; i1 += gridDim.y){
15801574 for (int i2 = threadIdx.x; i2 < d2; i2 += blockDim.x){
15811575 for (int i3 = threadIdx.y; i3 < d3; i3 += blockDim.y){
15821576 switch (operator_num) {
15831577 case IADD:
15841578 a[i0*sA0 + i1*sA1 + i2*sA2 + i3*sA3]
15851579 += b[i0*sB0 + i1*sB1 + i2*sB2 + i3*sB3];
15861580 break;
15871581 case IDIV:
15881582 a[i0*sA0 + i1*sA1 + i2*sA2 + i3*sA3]
15891583 /= b[i0*sB0 + i1*sB1 + i2*sB2 + i3*sB3];
15901584 break;
15911585 case CPY:
15921586 a[i0*sA0 + i1*sA1 + i2*sA2 + i3*sA3]
15931587 = b[i0*sB0 + i1*sB1 + i2*sB2 + i3*sB3];
15941588 break;
15951589 }
15961590 }
15971591 }
15981592 }
15991593 }
16001594 }
16011595
16021596 template <int operator_num>
16031597 __global__ void k_ielem_6(const int d0, const int d1,
16041598 const int d2, const int d3,
16051599 const int d4, const int d5,
16061600 float* a, const int sA0, const int sA1,
16071601 const int sA2, const int sA3,
16081602 const int sA4, const int sA5,
16091603 const float* b, const int sB0, const int sB1,
16101604 const int sB2, const int sB3,
16111605 const int sB4, const int sB5
16121606 ){
16131607 for (int i0 = blockIdx.x; i0 < d0; i0 += gridDim.x){
16141608 for (int i1 = blockIdx.y; i1 < d1; i1 += gridDim.y){
16151609 for (int i2 = blockIdx.z; i2 < d2; i2 += gridDim.z){
16161610 for (int i3 = threadIdx.x; i3 < d3; i3 += blockDim.x){
16171611 for (int i4 = threadIdx.y; i4 < d4; i4 += blockDim.y){
16181612 for (int i5 = threadIdx.z; i5 < d5; i5 += blockDim.z){
16191613 switch (operator_num) {
16201614 case IADD:
16211615 a[i0*sA0 + i1*sA1 + i2*sA2 + i3*sA3 + i4*sA4 + i5*sA5]
16221616 += b[i0*sB0 + i1*sB1 + i2*sB2 + i3*sB3 + i4*sB4 + i5*sB5];
16231617 break;
16241618 case IDIV:
16251619 a[i0*sA0 + i1*sA1 + i2*sA2 + i3*sA3 + i4*sA4 + i5*sA5]
16261620 /= b[i0*sB0 + i1*sB1 + i2*sB2 + i3*sB3 + i4*sB4 + i5*sB5];
16271621 break;
16281622 case CPY:
16291623 a[i0*sA0 + i1*sA1 + i2*sA2 + i3*sA3 + i4*sA4 + i5*sA5]
16301624 = b[i0*sB0 + i1*sB1 + i2*sB2 + i3*sB3 + i4*sB4 + i5*sB5];
16311625 break;
16321626 }
16331627 }
16341628 }
16351629 }
16361630 }
16371631 }
16381632 }
16391633 }
16401634
16411635 /*
16421636 CudaNdarray_inplace_elemwise
16431637 Compute elemwise, working inplace on A.
16441638 Currently implemented A / B, A + B and A = B
16451639 (the last is not tested and not used!)
16461640
16471641 py_self - the CudaNdarray that we'll modify (A)
16481642 py_other - the other argument (B)
16491643 fct_nb - which operation to perform (operator_t)
16501644
16511645 Returns 0 on success.
16521646 Returns -1 on failure, and sets Python exception.
16531647
16541648 */
16551649 int
16561650 CudaNdarray_inplace_elemwise(PyObject* py_self, PyObject * py_other, operator_t fct_nb)
16571651 {
16581652 int verbose = 0;
16591653 void (*k3)(const int, const int, const int,
16601654 float*, const int, const int, const int,
16611655 const float*, const int, const int, const int);
16621656 void (*k4)(const int, const int, const int, const int,
16631657 float*, const int, const int,
16641658 const int, const int,
16651659 const float*, const int, const int,
16661660 const int, const int);
16671661 void (*k6)(const int, const int,
16681662 const int, const int,
16691663 const int, const int,
16701664 float*, const int, const int,
16711665 const int, const int,
16721666 const int, const int,
16731667 const float*, const int, const int,
16741668 const int, const int,
16751669 const int, const int);
16761670 switch (fct_nb)
16771671 {
16781672 case IADD:
16791673 k3 = k_ielem_3<IADD>;
16801674 k4 = k_ielem_4<IADD>;
16811675 k6 = k_ielem_6<IADD>;
16821676 break;
16831677 case IDIV:
16841678 k3 = k_ielem_3<IDIV>;
16851679 k4 = k_ielem_4<IDIV>;
16861680 k6 = k_ielem_6<IDIV>;
16871681 break;
16881682 case CPY:
16891683 k3 = k_ielem_3<CPY>;
16901684 k4 = k_ielem_4<CPY>;
16911685 k6 = k_ielem_6<CPY>;
16921686 break;
16931687 default:
16941688 assert (0);
16951689 PyErr_Format(
16961690 PyExc_TypeError,
16971691 "CudaNdarray_inplace_elemwise invalid fct_nb (%i).",
16981692 (int)fct_nb);
16991693 return -1;
17001694 }
17011695 if (!CudaNdarray_Check(py_self)) {
17021696 PyErr_SetString(
17031697 PyExc_TypeError,
17041698 "CudaNdarray_inplace_elemwise need a CudaNdarray on left");
17051699 return -1;
17061700 }
17071701 CudaNdarray * new_other = NULL;
17081702 if (!CudaNdarray_Check(py_other)) {
17091703 new_other = (CudaNdarray*) CudaNdarray_New();
17101704 if(!new_other)
17111705 {
17121706 return -1;
17131707 }
17141708 if(CudaNdarray_CopyFromArray(new_other, (PyArrayObject *) py_other))
17151709 {
17161710 Py_XDECREF(new_other);
17171711 return -1;
17181712 }
17191713 py_other = (PyObject *) new_other;
17201714 }
17211715
17221716 CudaNdarray * self = (CudaNdarray *)py_self;
17231717 CudaNdarray * other = (CudaNdarray *)py_other;
17241718
17251719 if (verbose)
17261720 {
17271721 fprintf(stderr,
17281722 "INPLACE ADD/DIV for self->nd=%d other->nd=%d\n",
17291723 self->nd, other->nd);
17301724 }
17311725
17321726 //standard elemwise nb dim checks
17331727 if (self->nd < other->nd)
17341728 {
17351729 PyErr_Format(
17361730 PyExc_TypeError,
17371731 "CudaNdarray_inplace_elemwise: The destination need more or the"
17381732 " same number of dimensions then the source. Got %d and %d.",
17391733 self->nd, other->nd);
17401734 Py_XDECREF(new_other);
17411735 return -1;
17421736 }
17431737
17441738 //broadcast to the same number of dimensions.
17451739 int* other_dims = (int*) alloca(self->nd * sizeof(int));
17461740 int* other_strides = (int*) alloca(self->nd * sizeof(int));
17471741 int added_dims = self->nd - other->nd;
17481742 // Add the added broadcasted dimensions
17491743 for (int i = 0; i< added_dims; ++i)
17501744 {
17511745 other_dims[i] = 1;
17521746 other_strides[i] = 0;
17531747 }
17541748 // Copy the existing dimensions
17551749 for (int i = 0; i< other->nd; ++i)
17561750 {
17571751 other_dims[i+added_dims] = CudaNdarray_HOST_DIMS(other)[i];
17581752 other_strides[i+added_dims] = CudaNdarray_HOST_STRIDES(other)[i];
17591753 }
17601754
17611755 //standard elemwise dim checks
17621756 unsigned int size = 1;
17631757 for (int i = 0; i< self->nd; ++i)
17641758 {
17651759 if ((CudaNdarray_HOST_DIMS(self)[i] != other_dims[i])
17661760 && (other_dims[i] != 1))
17671761 {
17681762 PyErr_SetString(
17691763 PyExc_ValueError,
17701764 "CudaNdarray_inplace_elemwise need same dimensions (or broadcastable dimension)");
17711765 Py_XDECREF(new_other);
17721766 return -1;
17731767 }
17741768 // if we're broadcasting other, then make sure it has stride 0
17751769 assert ((CudaNdarray_HOST_DIMS(self)[i] == other_dims[i])
17761770 || (other_strides[i] == 0));
17771771 size *= (unsigned int) CudaNdarray_HOST_DIMS(self)[i];
17781772 }
17791773
17801774 if (size==0)
17811775 {
17821776 int other_size = CudaNdarray_SIZE((CudaNdarray *)py_other);
17831777 if (!(other_size == 0 || other_size == 1))
17841778 {
17851779 PyErr_SetString(
17861780 PyExc_ValueError,
17871781 "CudaNdarray_inplace_elemwise cannot work inplace on"
17881782 " un-initialized array when the new value have more than"
17891783 " 0 or 1 broadcastable dimensions");
17901784 Py_XDECREF(new_other);
17911785 return 0;
17921786 }
17931787 Py_XDECREF(new_other);
17941788 return 0;
17951789 }
17961790
17971791 switch(self->nd)
17981792 {
17991793 case 0:
18001794 {
18011795 dim3 n_blocks(1, 1, 1);
18021796 dim3 n_threads(1);
18031797 k3<<<n_blocks, n_threads>>>(
18041798 1, //d0
18051799 1, //d1
18061800 1, //d2
18071801 CudaNdarray_DEV_DATA(self),
18081802 1, //strides
18091803 1,
18101804 1,
18111805 CudaNdarray_DEV_DATA(other),
18121806 1, //strides
18131807 1,
18141808 1);
18151809 CNDA_THREAD_SYNC;
18161810 cudaError_t err = cudaGetLastError();
18171811 if (cudaSuccess != err)
18181812 {
18191813 PyErr_Format(
18201814 PyExc_RuntimeError,
18211815 "CudaNdarray_inplace_elemwise case0: Cuda error: %s: %s.\n",
18221816 "k3",
18231817 cudaGetErrorString(err));
18241818 Py_XDECREF(new_other);
18251819 return -1;
18261820 }
18271821 }
18281822 break;
18291823 case 1:
18301824 {
18311825 dim3 n_blocks(1, 1, 1);
18321826 dim3 n_threads(
18331827 std::min(
18341828 CudaNdarray_HOST_DIMS(self)[0],
18351829 NUM_VECTOR_OP_THREADS_PER_BLOCK));
18361830 k3<<<n_blocks, n_threads>>>(
18371831 1, //dimensions
18381832 1,
18391833 CudaNdarray_HOST_DIMS(self)[0],
18401834 CudaNdarray_DEV_DATA(self),
18411835 1, //strides
18421836 1,
18431837 CudaNdarray_HOST_STRIDES(self)[0],
18441838 CudaNdarray_DEV_DATA(other),
18451839 1, //strides
18461840 1,
18471841 other_strides[0]);
18481842 CNDA_THREAD_SYNC;
18491843 cudaError_t err = cudaGetLastError();
18501844 if (cudaSuccess != err)
18511845 {
18521846 PyErr_Format(
18531847 PyExc_RuntimeError,
18541848 "CudaNdarray_inplace_elemwise case1: Cuda error: %s: %s.\n",
18551849 "k3",
18561850 cudaGetErrorString(err));
18571851 Py_XDECREF(new_other);
18581852 return -1;
18591853 }
18601854 }
18611855 break;
18621856 case 2:
18631857 {
18641858 //TODO: if both self and other are f-contiguous
18651859 // Then flip the block and thread dimensions
18661860 // to make contiguous reads & writes
18671861 dim3 n_blocks(1,
18681862 std::min(
18691863 CudaNdarray_HOST_DIMS(self)[0],
18701864 NUM_VECTOR_OP_BLOCKS));
18711865 dim3 n_threads(
18721866 std::min(
18731867 CudaNdarray_HOST_DIMS(self)[1],
18741868 NUM_VECTOR_OP_THREADS_PER_BLOCK));
18751869 k3<<<n_blocks, n_threads>>>(1,
18761870 CudaNdarray_HOST_DIMS(self)[0],
18771871 CudaNdarray_HOST_DIMS(self)[1],
18781872 CudaNdarray_DEV_DATA(self),
18791873 1,
18801874 CudaNdarray_HOST_STRIDES(self)[0],
18811875 CudaNdarray_HOST_STRIDES(self)[1],
18821876 CudaNdarray_DEV_DATA(other),
18831877 1,
18841878 other_strides[0],
18851879 other_strides[1]);
18861880 CNDA_THREAD_SYNC;
18871881 cudaError_t err = cudaGetLastError();
18881882 if (cudaSuccess != err)
18891883 {
18901884 PyErr_Format(
18911885 PyExc_RuntimeError,
18921886 "CudaNdarray_inplace_elemwise case2: Cuda error: %s: %s.\n",
18931887 "k3",
18941888 cudaGetErrorString(err));
18951889 Py_XDECREF(new_other);
18961890 return -1;
18971891 }
18981892 }
18991893 break;
19001894 case 3:
19011895 {
19021896 //TODO: Dimshuffle so that at least one of the arrays
19031897 // has a contiguous dimension on the thread idx.
19041898 dim3 n_blocks(
19051899 std::min(
19061900 CudaNdarray_HOST_DIMS(self)[0],
19071901 NUM_VECTOR_OP_BLOCKS),
19081902 CudaNdarray_HOST_DIMS(self)[1]);
19091903 while (n_blocks.x * n_blocks.y > NUM_VECTOR_OP_BLOCKS)
19101904 n_blocks.y /= 2;
19111905 dim3 n_threads(
19121906 std::min(
19131907 CudaNdarray_HOST_DIMS(self)[2],
19141908 NUM_VECTOR_OP_THREADS_PER_BLOCK));
19151909 k3<<<n_blocks, n_threads>>>(
19161910 CudaNdarray_HOST_DIMS(self)[0],
19171911 CudaNdarray_HOST_DIMS(self)[1],
19181912 CudaNdarray_HOST_DIMS(self)[2],
19191913 CudaNdarray_DEV_DATA(self),
19201914 CudaNdarray_HOST_STRIDES(self)[0],
19211915 CudaNdarray_HOST_STRIDES(self)[1],
19221916 CudaNdarray_HOST_STRIDES(self)[2],
19231917 CudaNdarray_DEV_DATA(other),
19241918 other_strides[0],
19251919 other_strides[1],
19261920 other_strides[2]);
19271921 CNDA_THREAD_SYNC;
19281922 cudaError_t err = cudaGetLastError();
19291923 if (cudaSuccess != err)
19301924 {
19311925 PyErr_Format(
19321926 PyExc_RuntimeError,
19331927 "CudaNdarray_inplace_elemwise case3: Cuda error: %s: %s.\n",
19341928 "k3",
19351929 cudaGetErrorString(err));
19361930 Py_XDECREF(new_other);
19371931 return -1;
19381932 }
19391933 }
19401934 break;
19411935 case 4:
19421936 {
19431937 dim3 n_blocks(
19441938 std::min(
19451939 CudaNdarray_HOST_DIMS(self)[0],
19461940 NUM_VECTOR_OP_BLOCKS),
19471941 CudaNdarray_HOST_DIMS(self)[1]
19481942 );
19491943 while (n_blocks.x * n_blocks.y > NUM_VECTOR_OP_BLOCKS)
19501944 n_blocks.y /= 2;
19511945 dim3 n_threads(
19521946 std::min(
19531947 CudaNdarray_HOST_DIMS(self)[2],
19541948 NUM_VECTOR_OP_THREADS_PER_BLOCK)
19551949 //TODO: DON"T YOU NEED OT PUT DIMS[3] in here???
19561950 );
19571951 k4<<<n_blocks, n_threads>>>(
19581952 CudaNdarray_HOST_DIMS(self)[0],
19591953 CudaNdarray_HOST_DIMS(self)[1],
19601954 CudaNdarray_HOST_DIMS(self)[2],
19611955 CudaNdarray_HOST_DIMS(self)[3],
19621956 CudaNdarray_DEV_DATA(self),
19631957 CudaNdarray_HOST_STRIDES(self)[0],
19641958 CudaNdarray_HOST_STRIDES(self)[1],
19651959 CudaNdarray_HOST_STRIDES(self)[2],
19661960 CudaNdarray_HOST_STRIDES(self)[3],
19671961 CudaNdarray_DEV_DATA(other),
19681962 other_strides[0],
19691963 other_strides[1],
19701964 other_strides[2],
19711965 other_strides[3]);
19721966 CNDA_THREAD_SYNC;
19731967 cudaError_t err = cudaGetLastError();
19741968 if (cudaSuccess != err)
19751969 {
19761970 PyErr_Format(
19771971 PyExc_RuntimeError,
19781972 "CudaNdarray_inplace_elemwise case4: Cuda error: %s: %s.\n",
19791973 "k4",
19801974 cudaGetErrorString(err));
19811975 Py_XDECREF(new_other);
19821976 return -1;
19831977 }
19841978 }
19851979 break;
19861980 case 5:
19871981 {
19881982 dim3 n_blocks(
19891983 std::min(
19901984 CudaNdarray_HOST_DIMS(self)[1],
19911985 NUM_VECTOR_OP_BLOCKS),
19921986 CudaNdarray_HOST_DIMS(self)[2]);
19931987 while (n_blocks.x * n_blocks.y > NUM_VECTOR_OP_BLOCKS)
19941988 n_blocks.y /= 2;
19951989 dim3 n_threads(
19961990 std::min(
19971991 CudaNdarray_HOST_DIMS(self)[3],
19981992 NUM_VECTOR_OP_THREADS_PER_BLOCK)
19991993 //TODO: DON"T YOU NEED OT PUT DIMS[3] in here???
20001994 );
20011995 for (int i = 0; i < CudaNdarray_HOST_DIMS(self)[0]; ++i)
20021996 {
20031997 k4<<<n_blocks, n_threads>>>(
20041998 CudaNdarray_HOST_DIMS(self)[1],
20051999 CudaNdarray_HOST_DIMS(self)[2],
20062000 CudaNdarray_HOST_DIMS(self)[3],
20072001 CudaNdarray_HOST_DIMS(self)[4],
20082002 CudaNdarray_DEV_DATA(self) + i * CudaNdarray_HOST_STRIDES(self)[0],
20092003 CudaNdarray_HOST_STRIDES(self)[1],
20102004 CudaNdarray_HOST_STRIDES(self)[2],
20112005 CudaNdarray_HOST_STRIDES(self)[3],
20122006 CudaNdarray_HOST_STRIDES(self)[4],
20132007 CudaNdarray_DEV_DATA(other) + i * other_strides[0],
20142008 other_strides[1],
20152009 other_strides[2],
20162010 other_strides[3],
20172011 other_strides[4]);
20182012 CNDA_THREAD_SYNC;
20192013 cudaError_t err = cudaGetLastError();
20202014 if( cudaSuccess != err)
20212015 {
20222016 PyErr_Format(
20232017 PyExc_RuntimeError,
20242018 "CudaNdarray_inplace_elemwise case5: Cuda error: %s: %s. n_block=(%ld,%ld) n_threads=%ld\n",
20252019 "k5 with loop over k4",
20262020 cudaGetErrorString(err),
20272021 (long) n_blocks.x, (long) n_blocks.y, (long) n_threads.x);
20282022 Py_XDECREF(new_other);
20292023 return -1;
20302024 }
20312025 }
20322026 }
20332027 break;
20342028 case 6:
20352029 {
20362030 dim3 n_blocks(
20372031 std::min(
20382032 CudaNdarray_HOST_DIMS(self)[0],
20392033 NUM_VECTOR_OP_BLOCKS),
20402034 CudaNdarray_HOST_DIMS(self)[1],
20412035 CudaNdarray_HOST_DIMS(self)[2]
20422036 );
20432037 while (n_blocks.x * n_blocks.y > NUM_VECTOR_OP_BLOCKS)
20442038 n_blocks.y /= 2;
20452039 // GTX285(compute capabilities 1.3) don't support n_blocks.z > 1
20462040 // (compute capabilities 2.0) support 65535 for n_blocks.z
20472041 //while (n_blocks.x * n_blocks.y * n_blocks.z > NUM_VECTOR_OP_BLOCKS)
20482042 // n_blocks.z /= 2;
20492043 n_blocks.z = 1;
20502044 dim3 n_threads(
20512045 std::min(
20522046 CudaNdarray_HOST_DIMS(self)[3],
20532047 NUM_VECTOR_OP_THREADS_PER_BLOCK)
20542048 //TODO: DON'T YOU NEED TO PUT DIMS[4] in here???
20552049 //TODO: DON'T YOU NEED TO PUT DIMS[5] in here???
20562050 );
20572051 k6<<<n_blocks, n_threads>>>(
20582052 CudaNdarray_HOST_DIMS(self)[0],
20592053 CudaNdarray_HOST_DIMS(self)[1],
20602054 CudaNdarray_HOST_DIMS(self)[2],
20612055 CudaNdarray_HOST_DIMS(self)[3],
20622056 CudaNdarray_HOST_DIMS(self)[4],
20632057 CudaNdarray_HOST_DIMS(self)[5],
20642058 CudaNdarray_DEV_DATA(self),
20652059 CudaNdarray_HOST_STRIDES(self)[0],
20662060 CudaNdarray_HOST_STRIDES(self)[1],
20672061 CudaNdarray_HOST_STRIDES(self)[2],
20682062 CudaNdarray_HOST_STRIDES(self)[3],
20692063 CudaNdarray_HOST_STRIDES(self)[4],
20702064 CudaNdarray_HOST_STRIDES(self)[5],
20712065 CudaNdarray_DEV_DATA(other),
20722066 other_strides[0],
20732067 other_strides[1],
20742068 other_strides[2],
20752069 other_strides[3],
20762070 other_strides[4],
20772071 other_strides[5]);
20782072 CNDA_THREAD_SYNC;
20792073 cudaError_t err = cudaGetLastError();
20802074 if (cudaSuccess != err)
20812075 {
20822076 PyErr_Format(
20832077 PyExc_RuntimeError,
20842078 "CudaNdarray_inplace_elemwise case6: Cuda error: %s: %s. n_blocks=(%ld, %ld, %ld) n_threads=(%ld)\n",
20852079 "k6",
20862080 cudaGetErrorString(err),
20872081 (long) n_blocks.x, (long) n_blocks.y, (long) n_blocks.z,
20882082 (long) n_threads.x);
20892083 Py_XDECREF(new_other);
20902084 return -1;
20912085 }
20922086 }
20932087 break;
20942088 default:
20952089 {
20962090 PyErr_Format(
20972091 PyExc_NotImplementedError,
20982092 "inplace_elemwise w nd=%i\n",
20992093 self->nd);
21002094 Py_XDECREF(new_other);
21012095 return -1;
21022096 }
21032097 }
21042098 if (verbose)
21052099 fprintf(stderr, "INPLACE ADD/DIV end\n");
21062100 Py_XDECREF(new_other);
21072101 return 0;
21082102 }
21092103
21102104 /*
21112105 * We need this inplace Add to support IncSubTensor
21122106 * It returns py_self on success with an additional reference. Else NULL.
21132107 */
21142108 // Will be called by __iadd__ in Python
21152109 PyObject *
21162110 CudaNdarray_inplace_add(PyObject* py_self, PyObject * py_other)
21172111 {
21182112 if (CudaNdarray_inplace_elemwise(py_self, py_other, IADD))
21192113 {
21202114 return NULL;
21212115 }
21222116 Py_INCREF(py_self);
21232117 return py_self;
21242118 }
21252119
21262120 /*
21272121 * We need this inplace div for cuda/tests/test_basic_ops.py:test_shared_options
21282122 * It returns py_self on success with an additional reference. Else NULL.
21292123 */
21302124 // Will be called by __idiv__ in Python
21312125 static PyObject *
21322126 CudaNdarray_inplace_div(PyObject* py_self, PyObject * py_other)
21332127 {
21342128 if (CudaNdarray_inplace_elemwise(py_self, py_other, IDIV))
21352129 {
21362130 return NULL;
21372131 }
21382132 Py_INCREF(py_self);
21392133 return py_self;
21402134 }
21412135
21422136 // The PyNumberMethods struct layout changed in a non-trivial way from 2 to 3.
21432137 #if PY_MAJOR_VERSION == 3
21442138 static PyNumberMethods CudaNdarrayNumberMethods =
21452139 {
21462140 (binaryfunc)CudaNdarray_add, //binaryfunc nb_add; __add__
21472141 0, //binaryfunc nb_subtract;
21482142 0, //binaryfunc nb_multiply;
21492143 0, //binaryfunc nb_remainder;
21502144 0, //binaryfunc nb_divmod;
21512145 0, //ternaryfunc nb_power;
21522146 0, //unaryfunc nb_negative;
21532147 0, //unaryfunc nb_positive;
21542148 0, //unaryfunc nb_absolute;
21552149 0, //inquiry nb_bool;
21562150 0, //unaryfunc nb_invert;
21572151 0, //binaryfunc nb_lshift;
21582152 0, //binaryfunc nb_rshift;
21592153 0, //binaryfunc nb_and;
21602154 0, //binaryfunc nb_xor;
21612155 0, //binaryfunc nb_or;
21622156 0, //unaryfunc nb_int;
21632157 0, //void *nb_reserved;
21642158 0, //unaryfunc nb_float;
21652159
21662160 (binaryfunc)CudaNdarray_inplace_add, //binaryfunc nb_inplace_add; __iadd__
21672161 0, //binaryfunc nb_inplace_subtract;
21682162 0, //binaryfunc nb_inplace_multiply;
21692163 0, //binaryfunc nb_inplace_remainder;
21702164 0, //ternaryfunc nb_inplace_power;
21712165 0, //binaryfunc nb_inplace_lshift;
21722166 0, //binaryfunc nb_inplace_rshift;
21732167 0, //binaryfunc nb_inplace_and;
21742168 0, //binaryfunc nb_inplace_xor;
21752169 0, //binaryfunc nb_inplace_or;
21762170
21772171 0, //binaryfunc nb_floor_divide;
21782172 0, //binaryfunc nb_true_divide;
21792173 0, //binaryfunc nb_inplace_floor_divide;
21802174 (binaryfunc)CudaNdarray_inplace_div, //binaryfunc nb_inplace_true_divide; __idiv__
21812175
21822176 0, //unaryfunc nb_index
21832177 };
21842178 #else
21852179 static PyNumberMethods CudaNdarrayNumberMethods =
21862180 {
21872181 (binaryfunc)CudaNdarray_add, //binaryfunc nb_add; __add__
21882182 0, //binaryfunc nb_subtract; __sub__
21892183 0, //binaryfunc nb_multiply; __mul__
21902184 0, //binaryfunc nb_divide; __div__
21912185 0, //binaryfunc nb_remainder; __mod__
21922186 0, //binaryfunc nb_divmod; __divmod__
21932187 0, //ternaryfunc nb_power; __pow__
21942188 0, //unaryfunc nb_negative; __neg__
21952189 0, //unaryfunc nb_positive; __pos__
21962190 0, //unaryfunc nb_absolute; __abs__
21972191 0, //inquiry nb_nonzero; __nonzero__ /* Used by PyObject_IsTrue */
21982192 0, //unaryfunc nb_invert; __invert__
21992193 0, //binaryfunc nb_lshift; __lshift__
22002194 0, //binaryfunc nb_rshift; __rshift__
22012195 0, //binaryfunc nb_and; __and__
22022196 0, //binaryfunc nb_xor; __xor__
22032197 0, //binaryfunc nb_or; __or__
22042198 0, //coercion nb_coerce; __coerce__ /* Used by the coerce() function */
22052199 0, //unaryfunc nb_int; __int__
22062200 0, //unaryfunc nb_long; __long__
22072201 0, //unaryfunc nb_float; __float__
22082202 0, //unaryfunc nb_oct; __oct__
22092203 0, //unaryfunc nb_hex; __hex__
22102204
22112205 /* Added in release 2.0 */
22122206 (binaryfunc)CudaNdarray_inplace_add, //binaryfunc nb_inplace_add; __iadd__
22132207 0, //binaryfunc nb_inplace_subtract; __isub__
22142208 0, //binaryfunc nb_inplace_multiply; __imul__
22152209 (binaryfunc)CudaNdarray_inplace_div, //binaryfunc nb_inplace_divide; __idiv__
22162210 0, //binaryfunc nb_inplace_remainder; __imod__
22172211 0, //ternaryfunc nb_inplace_power; __ipow__
22182212 0, //binaryfunc nb_inplace_lshift; __ilshift__
22192213 0, //binaryfunc nb_inplace_rshift; __irshift__
22202214 0, //binaryfunc nb_inplace_and; __iand__
22212215 0, //binaryfunc nb_inplace_xor; __ixor__
22222216 0, //binaryfunc nb_inplace_or; __ior__
22232217
22242218 /* Added in release 2.2 */
22252219 0, //binaryfunc nb_floor_divide; __floordiv__
22262220 0, //binaryfunc nb_true_divide; __truediv__
22272221 0, //binaryfunc nb_inplace_floor_divide; __ifloordiv__
22282222 0, //binaryfunc nb_inplace_true_divide; __itruediv__
22292223
22302224 #if PY_MINOR_VERSION > 4
22312225 /* Added in release 2.5 */
22322226 0 //unaryfunc nb_index; __index__
22332227 #endif
22342228 };
22352229 #endif
22362230
22372231
22382232 /////////////////////
22392233 // Mapping protocol
22402234 /////////////////////
22412235
22422236 // Will by called by __len__ in Python
22432237 static Py_ssize_t
22442238 CudaNdarray_len(PyObject * py_self)
22452239 {
22462240 CudaNdarray * self = (CudaNdarray*) py_self;
22472241 if (self->nd <= 0)
22482242 {
22492243 return (Py_ssize_t) 0;
22502244 }
22512245 else
22522246 {
22532247 return (Py_ssize_t) CudaNdarray_HOST_DIMS(self)[0];
22542248 }
22552249 }
22562250
22572251 // Will by called by __getitem__ in Python
22582252 PyObject *
22592253 CudaNdarray_Subscript(PyObject * py_self, PyObject * key)
22602254 {
22612255 int verbose = 0;
22622256 if (verbose) fprintf(stderr, "Subscript .... \n");
22632257 CudaNdarray * self = (CudaNdarray*) py_self;
22642258 PyObject * py_rval = NULL;
22652259 CudaNdarray * rval = NULL;
22662260 PyObject * intobj = NULL;
22672261
22682262 //PyObject_Print(key, stderr, 0);
22692263
22702264 if (key == Py_Ellipsis)
22712265 {
22722266 Py_INCREF(py_self);
22732267 return py_self;
22742268 }
22752269 if ((intobj=PyNumber_Int(key))) //INDEXING BY INTEGER
22762270 //else if (PyInt_Check(key)) //INDEXING BY INTEGER
22772271 {
22782272 int d_idx = PyInt_AsLong(intobj);
22792273 Py_DECREF(intobj); intobj=NULL;
22802274 //int d_idx = PyInt_AsLong(key);
22812275 if (self->nd == 0)
22822276 {
22832277 PyErr_SetString(PyExc_IndexError, "0-d arrays can't be indexed");
22842278 return NULL;
22852279 }
22862280 int d_dim = CudaNdarray_HOST_DIMS(self)[0];
22872281 int offset = 0;
22882282
22892283 if ((d_idx >= 0) && (d_idx < d_dim))
22902284 {
22912285 //normal indexing
22922286 offset += d_idx * CudaNdarray_HOST_STRIDES(self)[0];
22932287 }
22942288 else if ((d_idx < 0) && (d_idx >= -d_dim))
22952289 {
22962290 //end-based indexing
22972291 // d_idx is negative
22982292 offset += (d_dim + d_idx) * CudaNdarray_HOST_STRIDES(self)[0];
22992293 }
23002294 else
23012295 {
23022296 PyErr_Format(PyExc_IndexError,
23032297 "index out of bounds. Asked %d, but size of %d",
23042298 d_idx, d_dim);
23052299 return NULL;
23062300 }
23072301
23082302 //allocate our subtensor view
23092303 py_rval = CudaNdarray_new_nd(self->nd - 1);
23102304 rval = (CudaNdarray*) py_rval;
23112305 if (!rval) return NULL;
23122306 assert (0 == rval->data_allocated);
23132307
23142308 //initialize the view's data pointer to our own.
23152309 if (CudaNdarray_set_device_data(rval, CudaNdarray_DEV_DATA(self) + offset, self))
23162310 {
23172311 Py_DECREF(rval);
23182312 return NULL;
23192313 }
23202314 for (int d = 1; d < self->nd; ++d)
23212315 {
23222316 CudaNdarray_set_stride(rval, d-1, CudaNdarray_HOST_STRIDES(self)[d]);
23232317 CudaNdarray_set_dim(rval, d-1, CudaNdarray_HOST_DIMS(self)[d]);
23242318 }
23252319 }
23262320 else
23272321 {
23282322 PyErr_Clear();
23292323 }
23302324 if (PySlice_Check(key)) //INDEXING BY SLICE
23312325 {
23322326 if (verbose) fprintf(stderr, "by slice\n");
23332327 if (self->nd == 0)
23342328 {
23352329 PyErr_SetString(PyExc_ValueError, "cannot slice a 0-d array");
23362330 return NULL;
23372331 }
23382332
23392333 int d_dim = CudaNdarray_HOST_DIMS(self)[0];
23402334 Py_ssize_t start, stop, step, slen;
23412335 if (PySlice_GetIndicesEx(SLICE_CAST(key), d_dim, &start, &stop, &step, &slen))
23422336 {
23432337 if (verbose)
23442338 fprintf(stderr, "PySlice_GetIndicesEx failed\n");
23452339 return NULL;
23462340 }
23472341 if (verbose)
23482342 {
23492343 std::cerr << "start " << start << "\n";
23502344 std::cerr << "stop " << stop << "\n";
23512345 std::cerr << "step " << step << "\n";
23522346 std::cerr << "slen " << slen << "\n";
23532347 }
23542348
23552349 //allocate our subtensor view
23562350 py_rval = CudaNdarray_new_nd(self->nd);
23572351 rval = (CudaNdarray*) py_rval;
23582352 if (!rval) return NULL;
23592353 assert (0 == rval->data_allocated);
23602354
23612355
23622356 //initialize the view's data pointer to our own.
23632357 if (CudaNdarray_set_device_data(rval,
23642358 CudaNdarray_DEV_DATA(self) + start * CudaNdarray_HOST_STRIDES(self)[0],
23652359 self))
23662360 {
23672361 Py_DECREF(rval);
23682362 return NULL;
23692363 }
23702364 //initialize dimension 0 of rval
23712365 CudaNdarray_set_stride(rval, 0,
23722366 (slen == 1) ? 0 : step * CudaNdarray_HOST_STRIDES(self)[0]);
23732367 CudaNdarray_set_dim(rval, 0, slen);
23742368 if (verbose) std::cerr << "rval stride " << CudaNdarray_HOST_STRIDES(rval)[0] << "\n";
23752369 // initialize dimensions > 0 of rval
23762370 for (int d = 1; d < self->nd; ++d)
23772371 {
23782372 CudaNdarray_set_stride(rval, d, CudaNdarray_HOST_STRIDES(self)[d]);
23792373 CudaNdarray_set_dim(rval, d, CudaNdarray_HOST_DIMS(self)[d]);
23802374 }
23812375 }
23822376 if (PyTuple_Check(key)) //INDEXING BY TUPLE
23832377 {
23842378 if (verbose) fprintf(stderr, "by tuple\n");
23852379 //elements of the tuple can be either integers or slices
23862380 //the dimensionality of the view we will return is diminished for each slice in the tuple
23872381
23882382 if (PyTuple_Size(key) > self->nd)
23892383 {
23902384 PyErr_SetString(PyExc_IndexError, "index error");
23912385 return NULL;
23922386 }
23932387
23942388 //calculate the number of dimensions in the return value
23952389 int rval_nd = self->nd;
23962390 for (int d = 0; d < PyTuple_Size(key); ++d)
23972391 {
23982392 //On some paltform PyInt_Check(<type 'numpy.int64'>) return true, other it return false.
23992393 //So we use PyArray_IsAnyScalar that should covert everything.
24002394 rval_nd -= PyArray_IsAnyScalar(PyTuple_GetItem(key, d));
24012395 }
24022396
24032397 //allocate our subtensor view
24042398 py_rval = CudaNdarray_new_nd(rval_nd);
24052399 rval = (CudaNdarray*) py_rval;
24062400 if (!rval) return NULL;
24072401 assert (0 == rval->data_allocated);
24082402
24092403 //initialize the view's data pointer to our own.
24102404 if (CudaNdarray_set_device_data(rval, CudaNdarray_DEV_DATA(self), self))
24112405 {
24122406 Py_DECREF(rval);
24132407 return NULL;
24142408 }
24152409
24162410 // rval_d will refer to the current dimension in the rval.
24172411 // It will not be incremented for integer keys, but will be incremented for slice
24182412 // keys
24192413 int rval_d = 0;
24202414
24212415 for (int d = 0; d < self->nd; ++d)
24222416 {
24232417 // keys can be shorter than self->nd.
24242418 // when that happens, it means that the remaining dimensions are "full slices"
24252419 if (d >=PyTuple_Size(key))
24262420 {
24272421 CudaNdarray_set_stride(rval, rval_d, CudaNdarray_HOST_STRIDES(self)[d]);
24282422 CudaNdarray_set_dim(rval, rval_d, CudaNdarray_HOST_DIMS(self)[d]);
24292423 ++rval_d;
24302424 }
24312425 else
24322426 {
24332427 PyObject * key_d = PyTuple_GetItem(key, d);
24342428
24352429 if (PySlice_Check(key_d))
24362430 {
24372431 Py_ssize_t start, stop, step, slen;
24382432 if (PySlice_GetIndicesEx(SLICE_CAST(key_d), CudaNdarray_HOST_DIMS(self)[d], &start, &stop, &step, &slen))
24392433 {
24402434 Py_DECREF(rval);
24412435 return NULL;
24422436 }
24432437 rval->devdata += start * CudaNdarray_HOST_STRIDES(self)[d];
24442438 CudaNdarray_set_stride(rval, rval_d,
24452439 (slen == 1) ? 0 : step * CudaNdarray_HOST_STRIDES(self)[d]);
24462440 CudaNdarray_set_dim(rval, rval_d, slen);
24472441 if (0)
24482442 {
24492443 std::cerr << "start " << start << "\n";
24502444 std::cerr << "stop " << stop << "\n";
24512445 std::cerr << "step " << step << "\n";
24522446 std::cerr << "slen " << slen << "\n";
24532447 }
24542448 ++rval_d;
24552449 }
24562450 else if ((intobj=PyNumber_Int(key_d)))
24572451 {
24582452 assert(PyArray_IsAnyScalar(key_d));
24592453 int d_idx = PyInt_AsLong(intobj);
24602454 Py_DECREF(intobj);
24612455 intobj = NULL;
24622456 int d_dim = CudaNdarray_HOST_DIMS(self)[d];
24632457
24642458 if ((d_idx >= 0) && (d_idx < d_dim))
24652459 {
24662460 //normal indexing
24672461 rval->devdata += d_idx * CudaNdarray_HOST_STRIDES(self)[d];
24682462 }
24692463 else if ((d_idx < 0) && (d_idx >= -d_dim))
24702464 {
24712465 //end-based indexing
24722466 rval->devdata += (d_dim + d_idx) * CudaNdarray_HOST_STRIDES(self)[d];
24732467 }
24742468 else
24752469 {
24762470 PyErr_Format(PyExc_IndexError,
24772471 "index out of bounds. Asked %d for dimensions %d, but size of %d",
24782472 d_idx, d, d_dim);
24792473 Py_DECREF(rval);
24802474 return NULL;
24812475 }
24822476 }
24832477 else
24842478 {
24852479 PyErr_Clear(); // clear the error set by PyNumber_Int
24862480 PyErr_SetString(PyExc_IndexError, "index must be either int or slice");
24872481 Py_DECREF(rval);
24882482 return NULL;
24892483 }
24902484 }
24912485 }
24922486 }
24932487 if (py_rval)
24942488 {
24952489 if (verbose) fprint_CudaNdarray(stderr, self);
24962490 if (verbose) fprint_CudaNdarray(stderr, rval);
24972491 }
24982492 else
24992493 {
25002494 PyErr_SetString(PyExc_NotImplementedError, "Unknown key type");
25012495 return NULL;
25022496 }
25032497 return py_rval;
25042498 }
25052499
25062500 // Will by called by __setitem__ in Python
25072501 // See http://docs.python.org/dev/py3k/c-api/object.html#PyObject_SetItem
25082502 // Doesn't handle broadcasting, e.g. a[:] = 5
25092503 // Can only be assigned from a CudaNdarray on the right side
25102504 // Or a ndarray
25112505 // Or a python scalar with value 0 when the left side part is c contiguous.
25122506 static int
25132507 CudaNdarray_setitem(PyObject *o, PyObject *key, PyObject *value)
25142508 {
25152509 int verbose = 0;
25162510 if (verbose) fprintf(stderr, "CudaNdarray_setitem start\n");
25172511 // We try to copy directly into this CudaNdarray from the ndarray
25182512 CudaNdarray* rval = (CudaNdarray*)CudaNdarray_Subscript(o, key);
25192513 CudaNdarray* new_value = NULL;
25202514
25212515 if(!rval){
25222516 // CudaNdarray_Subscript failed and set the error msg.
25232517 Py_XDECREF(rval);
25242518 return -1;
25252519 }
25262520
25272521 if(rval != (CudaNdarray*)o &&
25282522 (rval->data_allocated ||
25292523 // The new array should have a base
25302524 !(((CudaNdarray*)rval)->base) ||
25312525 // If the original array has no base, the base of the new
25322526 // array should be the original one
25332527 (!((CudaNdarray*)o)->base && ((CudaNdarray*)rval)->base != o) ||
25342528 // Else, the two arrays should have the same base
25352529 (((CudaNdarray*)o)->base && ((CudaNdarray*)rval)->base != ((CudaNdarray*)o)->base)))
25362530 {
25372531 // This case shouldn't happen, based on what I see in Subscript
25382532 // but just in case it happens sometime in the future
25392533
25402534 PyErr_Format(PyExc_RuntimeError,
25412535 "__getitem__ must return a CudaNdarray that refers to"
25422536 " the original CudaNdarray, not a copy. rval.base=%p"
25432537 " o.base=%p o=%p",
25442538 (((CudaNdarray*)rval)->base), ((CudaNdarray*)o)->base, o);
25452539 Py_DECREF(rval);
25462540 return -1;
25472541 }
25482542
25492543 PyObject * intobj = NULL;
25502544 if (CudaNdarray_Check(o) && PyArray_Check(value)){
25512545 if (verbose)
25522546 fprintf(stderr,
25532547 "CudaNdarray_setitem dest is a CudaNdarray and"
25542548 " value is a ndarray\n");
25552549 new_value = (CudaNdarray*) CudaNdarray_New();
25562550 if(!new_value)
25572551 {
25582552 return -1;
25592553 }
25602554 if (CudaNdarray_CopyFromArray(new_value, (PyArrayObject *) value))
25612555 {
25622556 Py_XDECREF(new_value);
25632557 Py_XDECREF(rval);
25642558 return -1;
25652559 }
25662560 value = (PyObject *) new_value;
25672561 }
25682562 else if ((intobj=PyNumber_Int(value)))
25692563 {
25702564 if (verbose)
25712565 fprintf(stderr,
25722566 "CudaNdarray_setitem dest and value is a python number\n");
25732567 if(! CudaNdarray_is_c_contiguous(rval)){
25742568 PyErr_SetString(PyExc_NotImplementedError,
25752569 "CudaNdarray.__setitem__: When the new value is a scalar"
25762570 " of value 0 the part where we copy to must be c contiguous.");
25772571 Py_XDECREF(rval);
25782572 return -1;
25792573 }
25802574
25812575 long val = PyInt_AsLong(intobj);
25822576 Py_DECREF(intobj); intobj=NULL;
25832577 if (val == 0)
25842578 {
25852579 cudaError_t err = cudaMemset(rval->devdata, 0,
25862580 CudaNdarray_SIZE(rval) * sizeof(real));
25872581 Py_XDECREF(rval);
25882582 if (err)
25892583 {
25902584 // Clear the error flag, cudaMemset doesn't do it.
25912585 // Currently this returns the same thing as err, but if in future
25922586 // it returns something else I still don't see why we should ignore
25932587 // it. All we want to do here is reset the flag.
25942588 cudaGetLastError();
25952589 PyErr_SetString(PyExc_RuntimeError,
25962590 "CudaNdarray.__setitem__: cudaMemset failed");
25972591 return -1;
25982592 }
25992593 return 0;
26002594 } else {
26012595 Py_XDECREF(rval);
26022596 PyErr_SetString(PyExc_NotImplementedError,
26032597 "CudaNdarray.__setitem__: we support setting only python"
26042598 " scalar of value 0, numpy nd array and CudaNdarray.");
26052599 return -1;
26062600 }
26072601 }
26082602
26092603 PyErr_Clear(); // clear PyNumber_Int error.
26102604
26112605 if(!CudaNdarray_Check(o) || !CudaNdarray_Check(value))
26122606 {
26132607 PyErr_SetString(PyExc_TypeError,
26142608 "CudaNdarray.__setitem__: left must be a CudaNdarrays and right"
26152609 " must be a CudaNdarrays, an ndarray or a python scalar of value 0.");
26162610 Py_XDECREF(new_value);
26172611 return -1;
26182612 }
26192613
26202614 if (verbose)
26212615 fprintf(stderr, "CudaNdarray_setitem dest and value are CudaNdarray\n");
26222616
26232617 if (cnda_copy_structure_to_device(rval))
26242618 {
26252619 PyErr_SetString(PyExc_RuntimeError,
26262620 "CudaNdarray.__setitem__: syncing structure to device failed");
26272621 Py_DECREF(rval);
26282622 Py_XDECREF(new_value);
26292623
26302624 if (verbose)
26312625 fprintf(stderr, "CudaNdarray_setitem error end\n");
26322626 return -1;
26332627 }
26342628
26352629 PyObject *baseSavedForComparison = rval->base;
26362630
26372631 if (CudaNdarray_CopyFromCudaNdarray(rval, (CudaNdarray*)value, true))
26382632 {
26392633 Py_DECREF((PyObject*)rval);
26402634 Py_XDECREF(new_value);
26412635
26422636 if (verbose)
26432637 fprintf(stderr, "CudaNdarray_setitem error end\n");
26442638 return -1;
26452639 }
26462640
26472641 assert (rval->base == baseSavedForComparison);
26482642 assert (rval->dev_structure_fresh);
26492643
26502644 // Clean up locally-created references
26512645 Py_DECREF(rval);
26522646 Py_XDECREF(new_value);
26532647
26542648 return 0;
26552649 }
26562650
26572651
26582652 PyMappingMethods CudaNdarrayMappingMethods = {
26592653 CudaNdarray_len, //lenfunc mp_length; __len__
26602654 CudaNdarray_Subscript, //binaryfunc mp_subscript; __getitem__
26612655 CudaNdarray_setitem //objobjargproc mp_ass_subscript; __setitem__
26622656 };
26632657
26642658 ////////////////////
26652659 //
26662660 ////////////////////
26672661
26682662 static PyObject *
26692663 CudaNdarray_get_shape(CudaNdarray *self, void *closure)
26702664 {
26712665 if (self->nd < 0)
26722666 {
26732667 PyErr_SetString(PyExc_ValueError, "CudaNdarray not initialized");
26742668 return NULL;
26752669 }
26762670 PyObject * rval = PyTuple_New(self->nd);
26772671 for (int i = 0; i < self->nd; ++i)
26782672 {
26792673 if (!rval || PyTuple_SetItem(rval, i, PyInt_FromLong(CudaNdarray_HOST_DIMS(self)[i])))
26802674 {
26812675 Py_XDECREF(rval);
26822676 return NULL;
26832677 }
26842678
26852679 }
26862680 return rval;
26872681 }
26882682
26892683 static int
26902684 CudaNdarray_set_shape(CudaNdarray *self, PyObject *value, void *closure)
26912685 {
26922686 PyErr_SetString(PyExc_NotImplementedError, "TODO: call reshape");
26932687 return -1;
26942688 }
26952689
26962690 static PyObject *
26972691 CudaNdarray_get_strides(CudaNdarray *self, void *closure)
26982692 {
26992693 if (self->nd < 0)
27002694 {
27012695 PyErr_SetString(PyExc_ValueError, "CudaNdarray not initialized");
27022696 return NULL;
27032697 }
27042698 PyObject * rval = PyTuple_New(self->nd);
27052699 for (int i = 0; i < self->nd; ++i)
27062700 {
27072701 if (!rval || PyTuple_SetItem(rval, i, PyInt_FromLong(CudaNdarray_HOST_STRIDES(self)[i])))
27082702 {
27092703 Py_XDECREF(rval);
27102704 return NULL;
27112705 }
27122706
27132707 }
27142708 return rval;
27152709 }
27162710
27172711 static int
27182712 CudaNdarray_set_strides(CudaNdarray *self, PyObject *value, void *closure)
27192713 {
27202714 //npy_intp newstrides_bytes[PyTuple_Size(value)];
27212715 if (PyTuple_Check(value)){
27222716 if (PyTuple_Size(value) != CudaNdarray_NDIM(self)){
27232717 PyErr_SetString(PyExc_ValueError,
27242718 "The new strides tuple must have the same length"
27252719 " as the number of dimensions");
27262720 return -1;
27272721 }
27282722 }else if (PyList_Check(value)){
27292723 if (PyList_Size(value) != CudaNdarray_NDIM(self)){
27302724 PyErr_SetString(PyExc_ValueError,
27312725 "The new strides list must have the same length"
27322726 " as the number of dimensions");
27332727 return -1;
27342728 }
27352729 }else{
27362730 PyErr_SetString(PyExc_ValueError,
27372731 "The new strides need to be encoded in a tuple or list");
27382732 return -1;
27392733 }
27402734 npy_intp* newstrides = (npy_intp*) alloca(CudaNdarray_NDIM(self) * sizeof(npy_intp));
27412735 if (PyTuple_Check(value)){
27422736 for(int i=0; i < CudaNdarray_NDIM(self); i++){
27432737 newstrides[i] = PyInt_AsLong(PyTuple_GetItem(value, Py_ssize_t(i)));
27442738 //newstrides_bytes[i] = newstrides[i] * 4;
27452739 }
27462740 }else if (PyList_Check(value)){
27472741 for(int i=0; i < CudaNdarray_NDIM(self); i++){
27482742 newstrides[i] = PyInt_AsLong(PyList_GetItem(value, Py_ssize_t(i)));
27492743 //newstrides_bytes[i] = newstrides[i] * 4;
27502744 }
27512745 }
27522746 /*
27532747 // Do not do this check, as ExtractDiag needs that, and NumPy does not seem
27542748 // to do it.
27552749 npy_intp dims[PyTuple_Size(value)];
27562750 for(int i=0; i < CudaNdarray_NDIM(self); i++){
27572751 dims[i] = CudaNdarray_HOST_DIMS(self)[i];
27582752 }
27592753 if (!PyArray_CheckStrides(4,
27602754 CudaNdarray_NDIM(self),
27612755 0, 0,
27622756 dims,
27632757 newstrides_bytes)){
27642758 PyErr_SetString(PyExc_ValueError, "bad new strides");
27652759 return -1;
27662760 }
27672761 */
27682762 for(int i=0; i < CudaNdarray_NDIM(self); i++){
27692763 CudaNdarray_set_stride(self, i, newstrides[i]);
27702764 }
27712765 return 0;
27722766 }
27732767
27742768 static PyObject *
27752769 CudaNdarray_get_dev_data(CudaNdarray *self, void *closure)
27762770 {
27772771 float * p = CudaNdarray_DEV_DATA(self);
27782772 //printf("get_dev_data %p %li \n", p, (long int)p );
27792773 return PyInt_FromSize_t((size_t) CudaNdarray_DEV_DATA(self));
27802774 }
27812775
27822776 static int
27832777 CudaNdarray_set_dev_data(CudaNdarray *self, PyObject *value, void *closure)
27842778 {
27852779 Py_ssize_t newdevdata = PyInt_AsSsize_t(value);
27862780 //printf("set_dev_data %p %li \n",(float*)newdevdata ,newdevdata);
27872781 if (PyErr_Occurred())
27882782 {
27892783 return -1;
27902784 }
27912785 return CudaNdarray_set_device_data(self, (float*)newdevdata, (CudaNdarray*)self->base);
27922786 }
27932787
27942788 static PyObject *
27952789 CudaNdarray_get_dtype(CudaNdarray *self, void *closure)
27962790 {
27972791 return PyString_FromString("float32");
27982792 }
27992793
28002794 static PyObject *
28012795 CudaNdarray_get_ndim(CudaNdarray *self, void *closure)
28022796 {
28032797 return PyInt_FromLong(self->nd);
28042798 }
28052799
28062800 static PyObject *
28072801 CudaNdarray_get_base(CudaNdarray *self, void *closure)
28082802 {
28092803 PyObject * base = self->base;
28102804 if (!base)
28112805 {
28122806 // We cannot return a NULL pointer, use None instead
28132807 base = Py_None;
28142808 }
28152809 Py_INCREF(base);
28162810 return base;
28172811 }
28182812
28192813 void put_in_dict(PyObject * dict, const char * key, int val)
28202814 {
28212815 PyObject * k = PyString_FromString(key);
28222816 PyObject * v = PyInt_FromLong(val);
28232817 PyDict_SetItem(dict, k, v);
28242818 Py_DECREF(k);
28252819 Py_DECREF(v);
28262820 }
28272821
28282822 PyObject *
28292823 GetDeviceProperties(PyObject* _unused, PyObject* args)
28302824 {
28312825 int dev_id = -1;
28322826 if (! PyArg_ParseTuple(args, "i", &dev_id))
28332827 return NULL;
28342828 cudaDeviceProp deviceProp;
28352829 cudaGetDeviceProperties(&deviceProp, dev_id);
28362830
28372831 PyObject * dict = PyDict_New();
28382832 PyObject * str= PyString_FromString("name");
28392833 PyObject * i = PyString_FromString(deviceProp.name);
28402834 PyDict_SetItem(dict, str, i);
28412835 Py_DECREF(str);
28422836 Py_DECREF(i);
28432837
28442838 put_in_dict(dict, "major", deviceProp.major);
28452839 put_in_dict(dict, "minor", deviceProp.minor);
28462840 #if CUDART_VERSION >= 2020
28472841 int driverVersion = 0, runtimeVersion = 0;
28482842 cudaDriverGetVersion(&driverVersion);
28492843 cudaRuntimeGetVersion(&runtimeVersion);
28502844 put_in_dict(dict, "driverVersion", driverVersion);
28512845 put_in_dict(dict, "runtimeVersion", runtimeVersion);
28522846 #endif
28532847 #if CUDART_VERSION >= 2000
28542848
28552849 put_in_dict(dict, "multiProcessorCount", deviceProp.multiProcessorCount);
28562850 //if ConvertSMVer2Cores is not defined in cuda_runtime_api.h, the run time is too old.
28572851 int sm_cores = -1;
28582852 if(deviceProp.major==1)
28592853 sm_cores = 32;
28602854 else if(deviceProp.major==2 && deviceProp.minor==0)
28612855 sm_cores = 32;
28622856 else if(deviceProp.major==2 && deviceProp.minor==1)
28632857 sm_cores = 48;
28642858 put_in_dict(dict, "coresCount", sm_cores * deviceProp.multiProcessorCount);
28652859 #endif
28662860 put_in_dict(dict, "totalConstMem", deviceProp.totalConstMem);
28672861 put_in_dict(dict, "sharedMemPerBlock", deviceProp.sharedMemPerBlock);
28682862 put_in_dict(dict, "regsPerBlock", deviceProp.regsPerBlock);
28692863 put_in_dict(dict, "warpSize", deviceProp.warpSize);
28702864 put_in_dict(dict, "maxThreadsPerBlock", deviceProp.maxThreadsPerBlock);
28712865 put_in_dict(dict, "maxThreadsDim0", deviceProp.maxThreadsDim[0]);
28722866 put_in_dict(dict, "maxThreadsDim1", deviceProp.maxThreadsDim[1]);
28732867 put_in_dict(dict, "maxThreadsDim2", deviceProp.maxThreadsDim[2]);
28742868 put_in_dict(dict, "maxGridSize0", deviceProp.maxGridSize[0]);
28752869 put_in_dict(dict, "maxGridSize1", deviceProp.maxGridSize[1]);
28762870 put_in_dict(dict, "maxGridSize2", deviceProp.maxGridSize[2]);
28772871 put_in_dict(dict, "memPitch", deviceProp.memPitch);
28782872 put_in_dict(dict, "textureAlignment", deviceProp.textureAlignment);
28792873 put_in_dict(dict, "clockRate", deviceProp.clockRate);
28802874 #if CUDART_VERSION >= 2000
28812875 put_in_dict(dict, "deviceOverlap", deviceProp.deviceOverlap);
28822876 #endif
28832877 #if CUDART_VERSION >= 2020
28842878 put_in_dict(dict, "kernelExecTimeoutEnabled", deviceProp.kernelExecTimeoutEnabled);
28852879 put_in_dict(dict, "integrated", deviceProp.integrated);
28862880 put_in_dict(dict, "canMapHostMemory", deviceProp.canMapHostMemory);
28872881 put_in_dict(dict, "computeMode", deviceProp.computeMode);
28882882 //in the doc of this fct tell that 0 - Normal mode, 1 - only 1 context, 2 - no context
28892883 #endif
28902884 #if CUDART_VERSION >= 3000
28912885 put_in_dict(dict, "concurrentKernels", deviceProp.concurrentKernels);
28922886 #endif
28932887 #if CUDART_VERSION >= 3010
28942888 put_in_dict(dict, "ECCEnabled", deviceProp.ECCEnabled);
28952889 #endif
28962890 #if CUDART_VERSION >= 3020
28972891 put_in_dict(dict, "tccDriver", deviceProp.tccDriver);
28982892 #endif
28992893
29002894 return dict;
29012895 }
29022896
29032897 /*
29042898 * Returns in *free and *total respectively, the free and total amount of memory available for allocation by the device in bytes.
29052899 */
29062900 PyObject *
29072901 GetDeviceMemInfo(PyObject* _unused, PyObject* dummy)
29082902 {
29092903 size_t free = 0, total = 0;
29102904 if(g_gpu_context_active == 0){
29112905 PyErr_Format(PyExc_RuntimeError, "No gpu device selected yet. Please make sure the gpu device was initialized by Theano before.");
29122906 return NULL;
29132907 }
29142908
29152909 cudaError_t err = cudaMemGetInfo(&free, &total);
29162910 if (err != cudaSuccess){
29172911 // Clear the error flag, cudaMemGetInfo doesn't do it.
29182912 // Currently this returns the same thing as err, but if in future
29192913 // it returns something else I still don't see why we should ignore
29202914 // it. All we want to do here is reset the flag.
29212915 cudaGetLastError();
29222916 PyErr_Format(PyExc_RuntimeError,
29232917 "Error while getting memory info about the gpu: %s",
29242918 cudaGetErrorString(err));
29252919 return NULL;
29262920 }
29272921 return PyTuple_Pack(2, PyLong_FromLong(free), PyLong_FromLong(total));
29282922 }
29292923
29302924 /*
29312925 * Synchronize with all the gpu device stream.
29322926 */
29332927 PyObject *
29342928 CudaNdarray_synchronize(PyObject* _unused, PyObject* dummy)
29352929 {
29362930 CNDA_BEGIN_ALLOW_THREADS
29372931 cudaThreadSynchronize();
29382932 CNDA_END_ALLOW_THREADS
29392933 Py_INCREF(Py_None);
29402934 return Py_None;
29412935 }
29422936
29432937 /*
29442938 * Exist and return true if we link with cublas v2.
29452939 */
29462940 PyObject *
29472941 CudaNdarray_cublasv2(PyObject* _unused, PyObject* dummy)
29482942 {
29492943 Py_INCREF(Py_True);
29502944 return Py_True;
29512945 }
29522946
29532947 PyObject *
29542948 CudaNdarray_select_a_gpu(PyObject* _unused, PyObject* dummy)
29552949 {
29562950 void * rval = NULL;
29572951 cudaError_t err;
29582952 int num_gpus = 0;
29592953
29602954 err = cudaGetDeviceCount(&num_gpus);
29612955 if (cudaSuccess != err){
29622956 printf("ERR!\\n");
29632957 PyErr_Format(PyExc_RuntimeError,
29642958 "Not able to get number of GPUs (%s).",
29652959 cudaGetErrorString(err));
29662960 return NULL;
29672961 }
29682962
29692963 for (int device = 0; device < num_gpus; device++) {
29702964 cudaSetDevice(device);
29712965 err = cudaDeviceSynchronize(); // << CUDA context gets created here.
29722966 cudaGetLastError(); // reset the error state
29732967 if (cudaSuccess == err)
29742968 break;
29752969 }
29762970
29772971 if (cudaSuccess != err){
29782972 printf("ERR!\\n");
29792973 PyErr_Format(PyExc_RuntimeError,
29802974 "Not able to select available GPU from %d cards (%s).",
29812975 num_gpus, cudaGetErrorString(err));
29822976 return NULL;
29832977 }
29842978
29852979 Py_INCREF(Py_None);
29862980 return Py_None;
29872981 }
29882982
29892983 #if COMPUTE_GPU_MEM_USED
29902984 /*
29912985 * Return the size in bytes that Theano currently have allocated on the gpu.
29922986 */
29932987 PyObject *
29942988 GetTheanoAllocInfo(PyObject* _unused, PyObject* dummy)
29952989 {
29962990 PyObject* a = PyLong_FromLong(_allocated_size);
29972991 PyObject* b = PyLong_FromLong(_max_allocated_size);
29982992
29992993 PyObject* tuple = PyTuple_New(2);
30002994 PyTuple_SetItem(tuple, 0, a);
30012995 PyTuple_SetItem(tuple, 1, b);
30022996 return tuple;
30032997 }
30042998 #endif
30052999
30063000 static PyGetSetDef CudaNdarray_getset[] = {
30073001 {"shape",
30083002 (getter)CudaNdarray_get_shape,
30093003 (setter)CudaNdarray_set_shape,
30103004 "shape of this ndarray (tuple)",
30113005 NULL},
30123006 {"_strides",
30133007 (getter)CudaNdarray_get_strides,
30143008 (setter)CudaNdarray_set_strides,
30153009 "data pointer strides (in elements)",
30163010 NULL},
30173011 {"strides",
30183012 (getter)CudaNdarray_get_strides,
30193013 (setter)CudaNdarray_set_strides,
30203014 "data pointer strides (in elements)",
30213015 NULL},
30223016 //gpudata is needed to allow calling pycuda fct with CudaNdarray input.
30233017 {"gpudata",
30243018 (getter)CudaNdarray_get_dev_data,
30253019 NULL,
30263020 "device data pointer",
30273021 NULL},
30283022 {"_dev_data",
30293023 (getter)CudaNdarray_get_dev_data,
30303024 (setter)CudaNdarray_set_dev_data,
30313025 "device data pointer",
30323026 NULL},
30333027 {"dtype",
30343028 (getter)CudaNdarray_get_dtype,
30353029 NULL,
30363030 "The dtype of the element. Now always float32",
30373031 NULL},
30383032 {"size",
30393033 (getter)CudaNdarray_SIZE_Object,
30403034 NULL,
30413035 "The number of elements in this object.",
30423036 NULL},
30433037 //mem_size is neede for pycuda.elementwise.ElementwiseKernel Why do they use size and mem_size of the same value?
30443038 {"mem_size",
30453039 (getter)CudaNdarray_SIZE_Object,
30463040 NULL,
30473041 "The number of elements in this object.",
30483042 NULL},
30493043 {"ndim",
30503044 (getter)CudaNdarray_get_ndim,
30513045 NULL,
30523046 "The number of dimensions in this object.",
30533047 NULL},
30543048 {"base",
30553049 (getter)CudaNdarray_get_base,
30563050 NULL,
30573051 "If this ndarray is a view, base is the original ndarray.",
30583052 NULL},
30593053
30603054 {NULL, NULL, NULL, NULL} /* Sentinel */
30613055 };
30623056
30633057 PyObject *CudaNdarray_repr(PyObject *self)
30643058 {
30653059 CudaNdarray *object = (CudaNdarray *)self;
30663060 PyObject * np_object = CudaNdarray_CreateArrayObj(object);
30673061 PyObject * str = PyObject_Str((PyObject *) np_object);
30683062 char * cstr = PyString_AsString(str);
30693063 PyObject * out = PyString_FromFormat("%s%s%s",
30703064 "CudaNdarray(",
30713065 cstr,
30723066 ")");
30733067 Py_DECREF(str);
30743068 Py_DECREF(np_object);
30753069 #if PY_MAJOR_VERSION >= 3
30763070 // In Python 3 PyString_FromFormat return a Bytes object
30773071 PyObject* out2 = PyObject_Str(out);
30783072 Py_DECREF(out);
30793073 return out2;
30803074 #endif
30813075 return out;
30823076 }
30833077
30843078 static PyTypeObject CudaNdarrayType =
30853079 {
30863080 #if PY_MAJOR_VERSION >= 3
30873081 PyVarObject_HEAD_INIT(NULL, 0)
30883082 #else
30893083 PyObject_HEAD_INIT(NULL)
30903084 0, /*ob_size*/
30913085 #endif
30923086 "CudaNdarray", /*tp_name*/
30933087 sizeof(CudaNdarray), /*tp_basicsize*/
30943088 0, /*tp_itemsize*/
30953089 (destructor)CudaNdarray_dealloc, /*tp_dealloc*/
30963090 0, /*tp_print*/
30973091 0, /*tp_getattr*/
30983092 0, /*tp_setattr*/
30993093 0, /*tp_compare*/
31003094 CudaNdarray_repr, /*tp_repr*/
31013095 &CudaNdarrayNumberMethods, /*tp_as_number*/
31023096 0, /*tp_as_sequence*/
31033097 &CudaNdarrayMappingMethods,/*tp_as_mapping*/
31043098 0, /*tp_hash */
31053099 0, /*tp_call*/
31063100 0, /*tp_str*/
31073101 0, /*tp_getattro*/
31083102 0, /*tp_setattro*/
31093103 0, /*tp_as_buffer*/
31103104 #if PY_MAJOR_VERSION >= 3
31113105 // Py_TPFLAGS_CHECKTYPES is always true and was removed in Python 3.
31123106 Py_TPFLAGS_DEFAULT | Py_TPFLAGS_BASETYPE, /*tp_flags*/
31133107 #else
31143108 Py_TPFLAGS_DEFAULT | Py_TPFLAGS_BASETYPE | Py_TPFLAGS_CHECKTYPES, /*tp_flags*/
31153109 #endif
31163110 "CudaNdarray objects", /* tp_doc */
31173111 0, /* tp_traverse */
31183112 0, /* tp_clear */
31193113 0, /* tp_richcompare */
31203114 0, /* tp_weaklistoffset */
31213115 0, /* tp_iter */
31223116 0, /* tp_iternext */
31233117 CudaNdarray_methods, /* tp_methods */
31243118 CudaNdarray_members, /* tp_members */
31253119 CudaNdarray_getset, /* tp_getset */
31263120 0, /* tp_base */
31273121 0, /* tp_dict */
31283122 0, /* tp_descr_get */
31293123 0, /* tp_descr_set */
31303124 0, /* tp_dictoffset */
31313125 (initproc)CudaNdarray_init,/* tp_init */
31323126 0, /* tp_alloc */
31333127 CudaNdarray_new, /* tp_new */
31343128 };
31353129
31363130 static __global__ void get_gpu_ptr_size(int* dst)
31373131 {
31383132 dst[0] = sizeof(float*);
31393133 dst[1] = sizeof(int);
31403134 }
31413135
31423136 PyObject *
31433137 CudaNdarray_ptr_int_size(PyObject* _unused, PyObject* args)
31443138 {
31453139 int *gpu_data = (int*)device_malloc(sizeof(int)*2);
31463140 if(gpu_data == NULL){
31473141 return NULL;
31483142 }
31493143 get_gpu_ptr_size<<<1,1>>>(gpu_data);
31503144
31513145 cudaError_t cudaErr = cudaGetLastError();
31523146 if (cudaSuccess != cudaErr){
31533147
31543148 device_free(gpu_data);
31553149 return PyErr_Format(PyExc_RuntimeError,
31563150 "CudaNdarray_ptr_int_size: error when calling the gpu code. (%s)",
31573151 cudaGetErrorString(cudaErr));
31583152 }
31593153
31603154 // Transfer the result to cpu
31613155 int gpu_sizes[] = {-1,-1};
31623156 cublasStatus_t err;
31633157 err = cublasGetVector(2, sizeof(int), gpu_data, 1, gpu_sizes, 1);
31643158 device_free(gpu_data);
31653159
31663160 if (CUBLAS_STATUS_SUCCESS != err){
31673161 PyErr_SetString(PyExc_RuntimeError, "error copying data to from memory");
31683162 return NULL;
31693163 }
31703164 return Py_BuildValue("iiii", (int) gpu_sizes[0], (int)sizeof(float*),
31713165 (int)sizeof(int), (int) gpu_sizes[1]);
31723166 }
31733167
31743168 static int cublas_init();
31753169 static void cublas_shutdown();
31763170 // Initialize the gpu.
31773171 // Takes two optional parameters, the device number and if we should use cnmem.
31783172 // If the device number is provided, it sets that device to be the active device.
31793173 // If not provided (usually just to test whether the gpu is available at all),
31803174 // it does not set an active device.
31813175 // Raises EnvironmentError or ValueError (as appropriate) if the initialization failed.
31823176 // cnmem is threaded like a bool. If converted to 0, don't use cnmem. Otherwise, use it.
31833177 PyObject *
31843178 CudaNdarray_gpu_init(PyObject* _unused, PyObject* args)
31853179 {
31863180 int card_nb = 0;
31873181 int card_number_provided = 1;
31883182 float cnmem = 0; // Theano flag lib.cnmem
31893183 // if we're given something wildly invalid, this will throw a TypeError
31903184 if(!PyArg_ParseTuple(args, "|if", &card_nb, &cnmem))
31913185 return NULL;
31923186 if(cnmem)
31933187 g_use_cnmem = true;
31943188
31953189 if(PyTuple_Size(args) == 0) {
31963190 card_number_provided = 0;
31973191 card_nb = 0;
31983192 }
31993193
32003194 int deviceCount;
32013195 cudaError err = cudaGetDeviceCount(&deviceCount);
32023196 if(cudaSuccess != err) {
32033197 return PyErr_Format(PyExc_EnvironmentError,
32043198 "Unable to get the number of gpus available: %s",
32053199 cudaGetErrorString(cudaGetLastError()));
32063200 }
32073201
32083202 // as soon as the first successful call to a cuda* function is made, a
32093203 // gpu context has been created
32103204 g_gpu_context_active = 1;
32113205
32123206 if(deviceCount <= 0) {
32133207 return PyErr_Format(PyExc_EnvironmentError,
32143208 "Can't use the GPU, no devices support CUDA");
32153209 }
32163210 if(card_number_provided && (card_nb < 0 || card_nb > (deviceCount - 1))) {
32173211 return PyErr_Format(PyExc_ValueError,
32183212 "Bad device number %d. Only %d devices available.",
32193213 card_nb,
32203214 deviceCount);
32213215 }
32223216
32233217 cudaDeviceProp deviceProp;
32243218 err = cudaGetDeviceProperties(&deviceProp, card_nb);
32253219 if(cudaSuccess != err) {
32263220 return PyErr_Format(PyExc_EnvironmentError,
32273221 "Unable to get properties of gpu %i: %s",
32283222 card_nb,
32293223 cudaGetErrorString(cudaGetLastError()));
32303224 }
32313225
32323226 if(deviceProp.major == 9999 && deviceProp.minor == 9999 ){
32333227 return PyErr_Format(PyExc_EnvironmentError,
32343228 "There is no device that supports CUDA");
32353229 }
32363230
32373231 if(card_number_provided) {
32383232 err = cudaSetDevice(card_nb);
32393233 if(cudaSuccess != err) {
32403234 return PyErr_Format(PyExc_EnvironmentError,
32413235 "Unable to set device %i: %s",
32423236 card_nb,
32433237 cudaGetErrorString(cudaGetLastError()));
32443238 }
32453239 if (cublas_init() == -1)
32463240 return NULL;
32473241 }
32483242 if(card_number_provided && g_use_cnmem) {
32493243 size_t mem = 0;
32503244 if (cnmem > 1)
32513245 mem = cnmem * 1024 * 1024;
32523246 else{
32533247 // Clip to 95% to let memory for the driver.
32543248 // 98% didn't worked in some cases.
32553249 if (cnmem > .95){
32563250 cnmem = .95;
32573251 }
32583252 size_t free = 0, total = 0;
32593253 cudaError_t err = cudaMemGetInfo(&free, &total);
32603254 if (err != cudaSuccess){
32613255 // Clear the error flag, cudaMemGetInfo doesn't do it.
32623256 // Currently this returns the same thing as err, but if in future
32633257 // it returns something else I still don't see why we should ignore
32643258 // it. All we want to do here is reset the flag.
32653259 cudaGetLastError();
32663260 PyErr_Format(PyExc_RuntimeError,
32673261 "Error while getting memory info about the gpu: %s",
32683262 cudaGetErrorString(err));
32693263 return NULL;
32703264 }
32713265 mem = total * cnmem;
32723266 }
32733267 if(initCnmem(card_number_provided, card_nb, mem) == -1){
32743268 return NULL;
32753269 }
32763270 }
32773271
32783272 Py_INCREF(Py_None);
32793273 return Py_None;
32803274 }
32813275
32823276 PyObject *
32833277 CudaNdarray_active_device_number(PyObject* _unused, PyObject* _unused_args) {
32843278 // NB: No cuda error checking here; keeps things simple, and it's not
32853279 // really necessary.
32863280 int currentDevice;
32873281 cudaGetDevice(¤tDevice);
32883282 return PyInt_FromLong(currentDevice);
32893283 }
32903284
32913285 PyObject *
32923286 CudaNdarray_active_device_name(PyObject* _unused, PyObject* _unused_args) {
32933287 // NB: No cuda error checking here; keeps things simple, and it's not
32943288 // really necessary.
32953289 int currentDevice;
32963290 cudaGetDevice(¤tDevice);
32973291
32983292 cudaDeviceProp deviceProp;
32993293 cudaGetDeviceProperties(&deviceProp, currentDevice);
33003294 return PyString_FromString(deviceProp.name);
33013295 }
33023296
33033297 PyObject *
33043298 CudaNdarray_gpu_shutdown(PyObject* _unused, PyObject* _unused_args) {
33053299 // Don't handle errors here
33063300 cublas_shutdown();
33073301 g_gpu_context_active = 0; // context has now been closed down
33083302 if(g_use_cnmem) {
33093303 cnmemStatus_t status = cnmemFinalize();
33103304 if(status != CNMEM_STATUS_SUCCESS) {
33113305 fprintf(stderr, "CudaNdarray_gpu_shutdown: cnmemFinalize failed! Reason=%s\n",
33123306 cnmemGetErrorString(status));
33133307 if(status == CNMEM_STATUS_CUDA_ERROR) {
33143308 fprintf(stderr, " Cuda-Reason=%s\n",
33153309 cudaGetErrorString(cudaGetLastError()));
33163310 }
33173311 }
33183312 }
33193313 cudaThreadExit();
33203314
33213315 Py_INCREF(Py_None);
33223316 return Py_None;
33233317 }
33243318
33253319 /*
33263320 * This function is tested in theano/misc/test_pycuda_theano_simple.py
33273321 */
33283322 PyObject *
33293323 CudaNdarray_from_gpu_pointer(PyObject* _unused, PyObject* args)
33303324 {
33313325 int verbose = 0;
33323326 PyObject *gpu_ptr = NULL;
33333327 PyObject *shapes = NULL;
33343328 PyObject *strides = NULL;
33353329 PyObject *base = NULL;
33363330 PyObject *rval = NULL;
33373331
33383332 //args should consist of 3 python objects
33393333 //The first is the gpu ptr
33403334 //The second if the shape
33413335 //The third if the strides
33423336 if (! PyArg_ParseTuple(args, "OOOO", &gpu_ptr, &shapes, &strides, &base))
33433337 return NULL;
33443338
33453339 if (verbose) printf("In CudaNdarray_from_gpu_pointer\n");
33463340 if (!PyLong_Check(gpu_ptr))
33473341 {
33483342 PyErr_Format(PyExc_Exception, "CudaNdarray_from_gpu_pointer: The gpu pointor is not an long");
33493343 return NULL;
33503344 }
33513345
33523346 Py_ssize_t nd = PyObject_Length(shapes);
33533347 if (nd < 0)
33543348 {
33553349 PyErr_SetString(PyExc_TypeError, "CudaNdarray_from_gpu_pointer: Couldn't get length of second argument");
33563350 return NULL;
33573351 }
33583352 Py_ssize_t nd_stride = PyObject_Length(strides);
33593353 if (nd_stride < 0)
33603354 {
33613355 PyErr_SetString(PyExc_TypeError, "CudaNdarray_from_gpu_pointer: Couldn't get length of third argument");
33623356 return NULL;
33633357 }
33643358
33653359 if (nd != nd_stride)
33663360 {
33673361 PyErr_SetString(PyExc_TypeError, "CudaNdarray_from_gpu_pointer: We need the same number of shapes and strides");
33683362 return NULL;
33693363 }
33703364
33713365 rval = CudaNdarray_New();
33723366
33733367 if (CudaNdarray_set_nd((CudaNdarray *)rval, nd))
33743368 {
33753369 //CudaNdarray_set_nd set the error msg
33763370 return NULL;
33773371 }
33783372 // set gpu pointeur
33793373 assert(((CudaNdarray *)rval)->data_allocated == 0);
33803374 if (CudaNdarray_set_device_data((CudaNdarray *)rval, (float *)PyInt_AsLong(gpu_ptr), base))
33813375 {
33823376 PyErr_SetString(PyExc_TypeError, "CudaNdarray_from_gpu_pointer: Error while setting the gpu pointor");
33833377 return NULL;
33843378
33853379 }
33863380
33873381 // Set dims and strides
33883382 for (int i = nd-1; i >= 0; --i)
33893383 {
33903384 PyObject * idx = PyLong_FromLong(i);
33913385 if (idx == NULL)
33923386 {
33933387 PyErr_SetString(PyExc_Exception, "CudaNdarray_from_gpu_pointer: Couldn't make long object to loop over list/tuple");
33943388 return NULL;
33953389 }
33963390 PyObject* dim_ = PyObject_GetItem(shapes, idx);
33973391 PyObject* strd_ = PyObject_GetItem(strides, idx);
33983392 if (!PyInt_Check(dim_))
33993393 {
34003394 PyErr_Format(PyExc_Exception, "CudaNdarray_from_gpu_pointer: shapes[%d] is not an int", i);
34013395 return NULL;
34023396 }
34033397 if (!PyInt_Check(strd_))
34043398 {
34053399 PyErr_Format(PyExc_Exception, "CudaNdarray_from_gpu_pointer: strides[%d] is not an int", i);
34063400 return NULL;
34073401 }
34083402 int dim = PyInt_AsLong(dim_);
34093403 int strd = PyInt_AsLong(strd_);
34103404 CudaNdarray_set_stride((CudaNdarray *)rval, i, strd);
34113405 CudaNdarray_set_dim((CudaNdarray *)rval, i, dim);
34123406 Py_DECREF(idx);
34133407 Py_DECREF(dim_);
34143408 Py_DECREF(strd_);
34153409 }
34163410 if (verbose) printf("CudaNdarray_from_gpu_pointer normal return\n");
34173411 return rval;
34183412 }
34193413
34203414 PyObject *
34213415 CudaNdarray_Dot(PyObject* _unused, PyObject* args)
34223416 {
34233417 PyObject *l=NULL;
34243418 PyObject *r=NULL;
34253419 PyObject * rval = NULL;
34263420
34273421 //args should consist of two python objects ("OO")
34283422 if (! PyArg_ParseTuple(args, "OO", &l, &r))
34293423 return NULL;
34303424
34313425 if (!CudaNdarray_Check(l) || !CudaNdarray_Check(r))
34323426 {
34333427 PyErr_SetString(PyExc_TypeError, "CudaNdarray arguments required ");
34343428 goto CudaNdarray_dot_fail;
34353429 }
34363430 if (((CudaNdarray*)l)->nd != 2)
34373431 {
34383432 PyErr_SetString(PyExc_TypeError, "need 2d CudaNdarray arg for now");
34393433 goto CudaNdarray_dot_fail;
34403434 }
34413435 if (((CudaNdarray*)r)->nd != 2)
34423436 {
34433437 PyErr_SetString(PyExc_TypeError, "need 2d CudaNdarray arg for now");
34443438 goto CudaNdarray_dot_fail;
34453439 }
34463440 rval = CudaNdarray_New();
34473441 if (!rval)
34483442 {
34493443 goto CudaNdarray_dot_fail;
34503444 }
34513445 int dims[2];
34523446 dims[0] = CudaNdarray_HOST_DIMS((CudaNdarray*)l)[0];
34533447 dims[1] = CudaNdarray_HOST_DIMS((CudaNdarray*)r)[1];
34543448 if (CudaNdarray_alloc_contiguous((CudaNdarray*)rval, 2, dims))
34553449 {
34563450 goto CudaNdarray_dot_fail;
34573451 }
34583452 if (CudaNdarray_gemm(1.0, (CudaNdarray*)l, (CudaNdarray*)r, 0.0, (CudaNdarray*)rval))
34593453 {
34603454 goto CudaNdarray_dot_fail;
34613455 }
34623456
34633457 return rval;
34643458
34653459 CudaNdarray_dot_fail:
34663460 Py_XDECREF(rval);
34673461 return NULL;
34683462 }
34693463
34703464 static PyObject *
34713465 filter(PyObject* __unsed_self, PyObject *args) // args = (data, broadcastable, strict, storage)
34723466 {
34733467 /*
34743468 * TODO: DOC what this function should do in the various cases of
34753469 * What is 'strict' supposed to mean in the context of this function?
34763470 * What do we do with input that could be interpreted as matching the broadcastable pattern in strict vs. non-strict cases?
34773471 *
34783472 */
34793473 PyObject *py_data=NULL;
34803474 PyArrayObject * data = NULL;
34813475 int strict = 0;
34823476 PyObject * broadcastable=NULL;
34833477 PyObject * storage=NULL;
34843478 CudaNdarray * rval=NULL;
34853479
34863480 //Python object references which are provided to the caller are borrowed references
34873481 if (!PyArg_ParseTuple(args, "OOiO", &py_data, &broadcastable, &strict, &storage)) return NULL;
34883482
34893483 if (!PyTuple_Check(broadcastable)){
34903484 PyErr_SetString(PyExc_TypeError, "broadcastable arg should be a tuple of int.");
34913485 return NULL;
34923486 }
34933487 Py_INCREF(py_data);
34943488 Py_INCREF(broadcastable);
34953489
34963490 CudaNdarray * cnda = (CudaNdarray*)py_data;
34973491
34983492 if (strict || CudaNdarray_Check(py_data))
34993493 {
35003494 //TODO: support non-strict "casting" from a vt to the broadcastable/type/size that we need.
35013495 if (!CudaNdarray_Check(py_data))
35023496 {
35033497 Py_DECREF(py_data);
35043498 Py_DECREF(broadcastable);
35053499 PyErr_SetString(PyExc_TypeError, "strict mode requires CudaNdarray");
35063500 return NULL;
35073501 }
35083502 if (cnda->nd != PyTuple_Size(broadcastable))
35093503 {
35103504 Py_DECREF(py_data);
35113505 Py_DECREF(broadcastable);
35123506 PyErr_Format(PyExc_TypeError, "Wrong rank: %i vs %li", cnda->nd, (long)PyTuple_Size(broadcastable));
35133507 return NULL;
35143508 }
35153509 for (int i = 0; i < cnda->nd; ++i)
35163510 {
35173511 if ((CudaNdarray_HOST_DIMS(cnda)[i] > 1) && PyInt_AsLong(PyTuple_GetItem(broadcastable, Py_ssize_t(i))))
35183512 {
35193513 PyErr_Format(PyExc_TypeError, "Non-unit size in broadcastable vt dimension %i", i);
35203514 Py_DECREF(py_data);
35213515 Py_DECREF(broadcastable);
35223516 return NULL;
35233517 }else if (CudaNdarray_HOST_DIMS(cnda)[i] == 1 && CudaNdarray_HOST_STRIDES(cnda)[i] != 0){
35243518 PyErr_Format(PyExc_TypeError, "Non-zeros strides(%d) on dimension %d of size 1",
35253519 CudaNdarray_HOST_STRIDES(cnda)[i], i);
35263520 Py_DECREF(py_data);
35273521 Py_DECREF(broadcastable);
35283522 return NULL;
35293523 }
35303524 }
35313525 Py_DECREF(broadcastable);
35323526 return py_data;
35333527 }
35343528 else
35353529 {
35363530 data = (PyArrayObject*)PyArray_FromObject(py_data, REAL_TYPENUM, PyTuple_Size(broadcastable), PyTuple_Size(broadcastable));
35373531 if (!data)
35383532 {
35393533 //err message already defined
35403534 Py_DECREF(py_data);
35413535 Py_DECREF(broadcastable);
35423536 return NULL;
35433537 }
35443538 for (int i = 0; i < PyArray_NDIM(data); ++i)
35453539 {
35463540 if ((PyArray_DIMS(data)[i] > 1) && PyInt_AsLong(PyTuple_GetItem(broadcastable, Py_ssize_t(i))))
35473541 {
35483542 PyErr_Format(PyExc_TypeError, "Non-unit size in broadcastable dimension %i", i);
35493543 Py_DECREF(data);
35503544 Py_DECREF(py_data);
35513545 Py_DECREF(broadcastable);
35523546 return NULL;
35533547 }
35543548 }
35553549 if (storage && CudaNdarray_Check(storage))
35563550 {
35573551 rval = (CudaNdarray*) storage;
35583552 Py_INCREF(rval);
35593553 }
35603554 else
35613555 {
35623556 rval = (CudaNdarray*) CudaNdarray_New();
35633557 }
35643558 if (rval)
35653559 {
35663560 if (CudaNdarray_CopyFromArray(rval, data))
35673561 {
35683562 Py_DECREF(rval);
35693563 rval = NULL;
35703564 }
35713565 }
35723566 Py_DECREF(data);
35733567 Py_DECREF(py_data);
35743568 Py_DECREF(broadcastable);
35753569 return (PyObject*)rval;
35763570 }
35773571 }
35783572
35793573 //TODO-- CudaNdarray_Dot and CudaNdarray_active_device_name are following different capitalization conventions.
35803574 // Pick one and standardize it, this file is already annoying enough to grep through
35813575 static PyMethodDef module_methods[] = {
35823576 {"dimshuffle", CudaNdarray_Dimshuffle, METH_VARARGS, "Returns the dimshuffle of a CudaNdarray."},
35833577 {"dot", CudaNdarray_Dot, METH_VARARGS, "Returns the matrix product of two CudaNdarray arguments."},
35843578 {"gpu_init", CudaNdarray_gpu_init, METH_VARARGS, "Select the gpu card to use; also usable to test whether CUDA is available."},
35853579 {"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."},
35863580 {"active_device_name", CudaNdarray_active_device_name, METH_VARARGS, "Get the name of the active device."},
35873581 {"active_device_number", CudaNdarray_active_device_number, METH_VARARGS, "Get the number of the active device."},
35883582 {"gpu_shutdown", CudaNdarray_gpu_shutdown, METH_VARARGS, "Shut down the gpu."},
35893583 {"device_properties", GetDeviceProperties, METH_VARARGS, "Return a dictionary with the device properties."},
35903584 {"mem_info", GetDeviceMemInfo, METH_NOARGS, "Return a tuple with the free and total memory on the gpu in bytes."},
35913585 #if COMPUTE_GPU_MEM_USED
35923586 {"theano_allocated", GetTheanoAllocInfo, METH_NOARGS, "Return the size in bytes of memory Theano currently have allocated on the gpu."},
35933587 #endif
35943588 {"ptr_int_size", CudaNdarray_ptr_int_size, METH_VARARGS, "Return a tuple with the size of gpu pointer, cpu pointer and int in bytes."},
35953589 {"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."},
35963590 {"outstanding_mallocs", outstanding_mallocs, METH_VARARGS, "how many more mallocs have been called than free's"},
35973591 {"from_gpu_pointer", CudaNdarray_from_gpu_pointer, METH_VARARGS, "Used to create a CudaNdarray from already allocated memory on the gpu.(example by pycuda)"},
35983592 {"synchronize", CudaNdarray_synchronize, METH_NOARGS, "Used to synchronize the device"},
35993593 {"cublas_v2", CudaNdarray_cublasv2, METH_NOARGS,
36003594 "Used to know if this version of cuda_ndarray is linked with cublas v2."},
36013595 {NULL, NULL, NULL, NULL} /* Sentinel */
36023596 };
36033597
36043598 #define CNDA_MOD_NAME "cuda_ndarray"
36053599 #define CNDA_DOCSTRING "CUDA implementation of a numpy ndarray-like object."
36063600
36073601 #if PY_MAJOR_VERSION == 3
36083602 static struct PyModuleDef cuda_ndarray_moduledef =
36093603 {
36103604 PyModuleDef_HEAD_INIT,
36113605 CNDA_MOD_NAME,
36123606 CNDA_DOCSTRING,
36133607 -1, /* size of per-interpreter state of the module,
36143608 or -1 if the module keeps state in global variables. */
36153609 module_methods
36163610 };
36173611
36183612 PyMODINIT_FUNC
36193613 PyInit_cuda_ndarray(void)
36203614 #else
36213615 PyMODINIT_FUNC
36223616 initcuda_ndarray(void)
36233617 #endif
36243618 {
36253619 import_array();
36263620
36273621 PyObject* m;
36283622
36293623 if (PyType_Ready(&CudaNdarrayType) < 0) {
36303624 #if PY_MAJOR_VERSION == 3
36313625 return NULL;
36323626 #else
36333627 return;
36343628 #endif
36353629 }
36363630
36373631 #if PY_MAJOR_VERSION == 3
36383632 m = PyModule_Create(&cuda_ndarray_moduledef);
36393633 #else
36403634 m = Py_InitModule3(CNDA_MOD_NAME, module_methods, CNDA_DOCSTRING);
36413635 #endif
36423636
36433637 if (m == NULL) {
36443638 #if PY_MAJOR_VERSION == 3
36453639 return NULL;
36463640 #else
36473641 return;
36483642 #endif
36493643 }
36503644
36513645 Py_INCREF(&CudaNdarrayType);
36523646 PyModule_AddObject(m, "CudaNdarray", (PyObject *)&CudaNdarrayType);
36533647 #if COMPUTE_GPU_MEM_USED
36543648 for(int i=0;i<TABLE_SIZE;i++){
36553649 _alloc_size_table[i].ptr=NULL;
36563650 _alloc_size_table[i].size=0;
36573651 }
36583652 #endif
36593653 // cublasInit();
36603654 //if (0&&CUBLAS_STATUS_SUCCESS != cublasGetError())
36613655 //{
36623656 //std::cerr << "WARNING: initcuda_ndarray: error initializing device\n";
36633657 //}
36643658 if (0) //TODO: is this necessary?
36653659 {
36663660 int deviceId = 0; // TODO: what number goes here?
36673661 cudaSetDevice(deviceId);
36683662 cudaError_t err = cudaGetLastError();
36693663 if( cudaSuccess != err)
36703664 {
36713665 std::cerr << "Error in SetDevice:" << cudaGetErrorString(err) << "\n";
36723666 }
36733667 }
36743668
36753669 #if PY_MAJOR_VERSION == 3
36763670 return m;
36773671 #endif
36783672 }
36793673
36803674
36813675 //////////////////////////////////////
36823676 //
36833677 // C API FOR CudaNdarray
36843678 //
36853679 //////////////////////////////////////
36863680
36873681 int
36883682 CudaNdarray_Check(const PyObject * ob)
36893683 {
36903684 //TODO: doesn't work with inheritance
36913685 return CudaNdarray_CheckExact(ob);
36923686 }
36933687 int
36943688 CudaNdarray_CheckExact(const PyObject * ob)
36953689 {
36963690 return ((Py_TYPE(ob) == &CudaNdarrayType) ? 1 : 0);
36973691 }
36983692
36993693 PyObject *
37003694 CudaNdarray_New(int nd)
37013695 {
37023696 CudaNdarray *self = (CudaNdarray *)CudaNdarrayType.tp_alloc(&CudaNdarrayType, 0);
37033697 if (self == NULL)
37043698 {
37053699 PyErr_SetString(PyExc_RuntimeError, "CudaNdarray_New failed to allocate self");
37063700 return NULL;
37073701 }
37083702 CudaNdarray_null_init(self);
37093703
37103704 if (nd == 0)
37113705 {
37123706 self->nd = 0;
37133707 }
37143708 else if (nd > 0)
37153709 {
37163710 if (CudaNdarray_set_nd(self, nd))
37173711 {
37183712 Py_DECREF(self);
37193713 return NULL;
37203714 }
37213715 }
37223716 ++_outstanding_mallocs[1];
37233717 return (PyObject *)self;
37243718 }
37253719
37263720
37273721
37283722 //////////////////////////////
37293723 //
37303724 // Published helper functions
37313725 //
37323726 //////////////////////////////
37333727
37343728 static int
37353729 cublas_init()
37363730 {
37373731 cublasStatus_t err;
37383732 err = cublasCreate(&handle);
37393733 if (CUBLAS_STATUS_SUCCESS != err)
37403734 {
37413735 if(CUBLAS_STATUS_NOT_INITIALIZED == err)
37423736 PyErr_SetString(PyExc_RuntimeError,
37433737 "cublasCreate() returned this error "
37443738 "'the CUDA Runtime initialization failed'");
37453739 else if(CUBLAS_STATUS_ALLOC_FAILED == err)
37463740 PyErr_SetString(PyExc_RuntimeError,
37473741 "cublasCreate() returned this error "
37483742 "'the resources could not be allocated'");
37493743 else
37503744 PyErr_SetString(PyExc_RuntimeError,
37513745 "unknow error during returned by cublasCreate()");
37523746 return -1;
37533747 }
37543748 // Set the default stream as the one to execute on (default)
37553749 cublasSetStream(handle, NULL);
37563750 // Pointer to scalars are on the host (also default)
37573751 cublasSetPointerMode(handle, CUBLAS_POINTER_MODE_HOST);
37583752 #if CUDA_VERSION >= 5000
37593753 // atomics can be used in kernels to speed up operations (not default)
37603754 // This may lead to a slight variance from run to run in some operations
37613755 cublasSetAtomicsMode(handle, CUBLAS_ATOMICS_ALLOWED);
37623756 #endif
37633757 return 0;
37643758 }
37653759
37663760 static void
37673761 cublas_shutdown()
37683762 {
37693763 if (handle != NULL)
37703764 cublasDestroy(handle);
37713765 // No point in handling any errors here
37723766 handle = NULL;
37733767 }
37743768
37753769 int
37763770 CudaNdarray_CopyFromArray(CudaNdarray * self, PyArrayObject*obj)
37773771 {
37783772 int err = CudaNdarray_alloc_contiguous(self, PyArray_NDIM(obj),
37793773 PyArray_DIMS(obj));
37803774 if (err) {
37813775 return err;
37823776 }
37833777
37843778 int typenum = PyArray_TYPE(obj);
37853779 if (typenum != REAL_TYPENUM)
37863780 {
37873781 PyErr_SetString(PyExc_TypeError, "can only copy from float arrays");
37883782 return -1;
37893783 }
37903784 assert( 4 == PyArray_ITEMSIZE(obj));
37913785 PyArrayObject * py_src = (PyArrayObject *)PyArray_ContiguousFromAny(
37923786 (PyObject*)obj, typenum, self->nd, self->nd);
37933787 if (!py_src) {
37943788 return -1;
37953789 }
37963790 npy_intp py_src_size = PyArray_SIZE(py_src);
37973791 void *py_src_data = PyArray_DATA(py_src);
37983792 cudaError_t cerr;
37993793 CNDA_BEGIN_ALLOW_THREADS;
38003794 cerr = cudaMemcpy(self->devdata, py_src_data,
38013795 py_src_size * sizeof(real),
38023796 cudaMemcpyHostToDevice);
38033797 //CNDA_THREAD_SYNC; // unneeded because cudaMemcpy is blocking anyway
38043798 CNDA_END_ALLOW_THREADS;
38053799 if (cudaSuccess != cerr)
38063800 {
38073801 PyErr_Format(PyExc_RuntimeError,
38083802 "Cuda error '%s' while copying %lli data element"
38093803 " to device memory",
38103804 cudaGetErrorString(cerr),
38113805 (long long)py_src_size);
38123806 Py_DECREF(py_src);
38133807 return -1;
38143808 }
38153809 Py_DECREF(py_src);
38163810 return 0;
38173811 }
38183812
38193813 PyObject *
38203814 CudaNdarray_new_nd(int nd)
38213815 {
38223816 CudaNdarray * rval = (CudaNdarray*) CudaNdarray_New();
38233817 if (!rval || CudaNdarray_set_nd(rval, nd))
38243818 {
38253819 Py_XDECREF(rval);
38263820 rval = NULL;
38273821 }
38283822 return (PyObject *) rval;
38293823 }
38303824
38313825
38323826 /**
38333827 * Initialize 'self' as a view of 'base', with memory storage 'data'
38343828 */
38353829
38363830 int CudaNdarray_set_device_data(CudaNdarray * self, float * data, PyObject * base)
38373831 {
38383832 if (self->data_allocated)
38393833 {
38403834 assert(self->devdata);
38413835 if (device_free(self->devdata))
38423836 {
38433837 self->devdata = NULL;
38443838 self->data_allocated = 0;
38453839 return -1;
38463840 }
38473841 }
38483842 // Get the original base object (base.base.base...)
38493843 PyObject * orig_base = base;
38503844 // base is not always a CudaNdarray. It can be a GpuArray from pycuda, ...
38513845 while (orig_base && CudaNdarray_Check(orig_base) && ((CudaNdarray*) orig_base)->base)
38523846 {
38533847 // base_base is itself a view
38543848 orig_base = ((CudaNdarray*) orig_base)->base;
38553849 }
38563850 //N.B. XDECREF and XINCREF are no-ops for NULL pointers
38573851 if (self->base != orig_base)
38583852 {
38593853 Py_XDECREF(self->base);
38603854 self->base = orig_base;
38613855 Py_XINCREF(self->base);
38623856 }
38633857 self->data_allocated = 0;
38643858 self->devdata = data;
38653859 return 0;
38663860 }
38673861
38683862 static __global__ void k_copy_1d(const int N, const float * x, const int sx, float * y, const int sy)
38693863 {
38703864 for (int i = threadIdx.x + blockIdx.x * blockDim.x; i < N; i += gridDim.x*blockDim.x)
38713865 {
38723866 y[i*sy] = x[i*sx];
38733867 }
38743868 }
38753869
38763870 // N1 through N4 are the size of y
38773871 static __global__ void k_copy_4d(const int N1,
38783872 const int N2, const int N3, const int N4,
38793873 const float * x, const int sx1, const int sx2, const int sx3,
38803874 const int sx4, float * y, const int sy1, const int sy2,
38813875 const int sy3, const int sy4)
38823876 {
38833877 // These must be made int instead of unsigned int due to a bug in nvcc
38843878 int bx = blockIdx.x;
38853879 int by = blockIdx.y;
38863880
38873881 for (int i = bx; i < N1; i += gridDim.x)
38883882 {
38893883 for (int j = by; j < N2; j += gridDim.y)
38903884 {
38913885 for (int k = threadIdx.x; k < N3; k += (int) blockDim.x)
38923886 {
38933887 for (int l = threadIdx.y; l < N4; l += (int) blockDim.y)
38943888 {
38953889 y[i * sy1 + j * sy2 + k * sy3 + l * sy4] =
38963890 x[i * sx1 + j * sx2 + k * sx3 + l * sx4];
38973891 }
38983892 }
38993893 }
39003894 }
39013895 }
39023896
39033897 //copy from other into self
39043898 int CudaNdarray_CopyFromCudaNdarray(CudaNdarray * self,
39053899 const CudaNdarray * other,
39063900 bool unbroadcast)
39073901 {
39083902 int verbose = 0;
39093903 if (verbose>1) fprintf(stderr, "CudaNdarray_CopyFromCudaNdarray\n");
39103904
39113905 //standard elemwise size checks
39123906 if (self->nd == -1)
39133907 {
39143908 PyErr_SetString(PyExc_TypeError,
39153909 "can't copy into un-initialized CudaNdarray");
39163910 return -1;
39173911 }
39183912 CudaNdarray * new_other = NULL;
39193913
39203914 if (self->nd < other->nd)
39213915 {
39223916 PyErr_Format(PyExc_NotImplementedError,
39233917 "CudaNdarray_CopyFromCudaNdarray: The number of dimensions of the "
39243918 "destination needs to be >= the number of dimensions of the "
39253919 "source. Got %d and %d.", self->nd, other->nd);
39263920 return -1;
39273921 }
39283922 else if (self->nd != other->nd)
39293923 {
39303924 new_other = (CudaNdarray *) CudaNdarray_View(other);
39313925 int added_dims = self->nd - other->nd;
39323926 int* pattern = (int*) alloca(self->nd * sizeof(int));
39333927 for(int i = 0; i < added_dims; i++)
39343928 pattern[i] = -1;
39353929 for(int i = 0; i < other->nd; i++)
39363930 pattern[i + added_dims] = i;
39373931 CudaNdarray_dimshuffle(new_other, self->nd, pattern);
39383932 other = new_other;
39393933 }
39403934 assert(self->nd == other->nd);
39413935 //standard elemwise dim checks (also compute total size)
39423936 unsigned int size = 1;
39433937 unsigned int size_source = 1;
39443938 for (int i = 0; i< self->nd; ++i)
39453939 {
39463940 if ((CudaNdarray_HOST_DIMS(self)[i] != CudaNdarray_HOST_DIMS(other)[i])
39473941 && (1!=CudaNdarray_HOST_DIMS(other)[i] || !unbroadcast) )
39483942 {
39493943 PyErr_Format(PyExc_ValueError,
39503944 "CudaNdarray_CopyFromCudaNdarray:"
39513945 " need same dimensions for dim %d,"
39523946 " destination=%d, source=%d",
39533947 i, CudaNdarray_HOST_DIMS(self)[i],
39543948 CudaNdarray_HOST_DIMS(other)[i]);
39553949 Py_XDECREF(new_other);
39563950 return -1;
39573951 }
39583952 size *= (unsigned int) CudaNdarray_HOST_DIMS(self)[i];
39593953 size_source *= (unsigned int) CudaNdarray_HOST_DIMS(other)[i];
39603954 }
39613955 if (0 == size)
39623956 {
39633957 Py_XDECREF(new_other);
39643958 return 0; //nothing to copy, we're done.
39653959 }
39663960 if (CudaNdarray_is_c_contiguous(self) &&
39673961 CudaNdarray_is_c_contiguous(other) &&
39683962 size == size_source)
39693963 {
39703964 if (verbose)
39713965 fprintf(stderr, "Copying contiguous vector with cublasScopy\n");
39723966
39733967 cublasStatus_t err;
39743968 err = cublasScopy(handle, size, CudaNdarray_DEV_DATA(other), 1,
39753969 CudaNdarray_DEV_DATA(self), 1);
39763970 CNDA_THREAD_SYNC;
39773971 Py_XDECREF(new_other);
39783972 if (CUBLAS_STATUS_SUCCESS != err)
39793973 {
39803974 PyErr_SetString(PyExc_RuntimeError, "Error copying memory");
39813975 return -1;
39823976 }
39833977 return 0;
39843978 }
39853979 //TODO: rewrite these copy operations to be more efficient
39863980 // See, for example the transpose example in the cuda_sdk.
39873981 switch (self->nd)
39883982 {
39893983 case 0: // scalar
39903984 {
39913985 // THIS CASE SHOULD NEVER HAPPEN BECAUSE SCALARS ARE ALWAYS C CONTIGUOUS
39923986 assert(0);
39933987 }; break;
39943988 case 1: // vector
39953989 {
39963990 if (verbose) fprintf(stderr, "Copying non-contiguous vector\n");
39973991 if (verbose) fprint_CudaNdarray(stderr, other);
39983992 unsigned int n_blocks = std::min(size,
39993993 (unsigned int)NUM_VECTOR_OP_BLOCKS);
40003994 unsigned int n_threads = std::min(ceil_intdiv(size, n_blocks),
40013995 (unsigned int)NUM_VECTOR_OP_THREADS_PER_BLOCK);
40023996 k_copy_1d<<<n_blocks, n_threads>>>(size,
40033997 CudaNdarray_DEV_DATA(other),
40043998 CudaNdarray_HOST_STRIDES(other)[0],
40053999 CudaNdarray_DEV_DATA(self),
40064000 CudaNdarray_HOST_STRIDES(self)[0]);
40074001 CNDA_THREAD_SYNC;
40084002 cudaError_t err = cudaGetLastError();
40094003 if( cudaSuccess != err)
40104004 {
40114005 PyErr_Format(PyExc_RuntimeError,
40124006 "Cuda error: %s: %s. (n_blocks=%i,"
40134007 " n_threads_per_block=%i)\n", "k_copy_1d",
40144008 cudaGetErrorString(err), n_blocks, n_threads);
40154009 Py_XDECREF(new_other);
40164010 return -1;
40174011 }
40184012 }; break;
40194013 case 4: // 4-tensor
40204014 {
40214015 if (verbose)
40224016 {
40234017 if (0 != fprint_CudaNdarray(stderr, other))
40244018 {
40254019 Py_XDECREF(new_other);
40264020 return -1;
40274021 }
40284022 }
40294023
40304024 // The blocks implement the looping over the first two axes so
40314025 // this needs to be (N1, N2)
40324026 dim3 n_blocks( std::min(CudaNdarray_HOST_DIMS(self)[0],
40334027 NUM_VECTOR_OP_BLOCKS),
40344028 std::min(CudaNdarray_HOST_DIMS(self)[1],
40354029 NUM_VECTOR_OP_BLOCKS));
40364030 // For the threads, just make as many as possible
40374031 dim3 n_threads( std::min( (unsigned int) CudaNdarray_HOST_DIMS(self)[2],
40384032 (unsigned int) NUM_VECTOR_OP_THREADS_PER_BLOCK),
40394033 std::min( (unsigned int) CudaNdarray_HOST_DIMS(self)[3],
40404034 (unsigned int) NUM_VECTOR_OP_THREADS_PER_BLOCK));
40414035
40424036 n_threads.x = std::min( (unsigned int) 32, (unsigned int) n_threads.x);
40434037 n_threads.y = std::min( n_threads.y, NUM_VECTOR_OP_THREADS_PER_BLOCK / n_threads.x);
40444038
40454039 k_copy_4d<<<n_blocks, n_threads>>>(
40464040 // size of y
40474041 (unsigned int) CudaNdarray_HOST_DIMS(self)[0], // N1
40484042 (unsigned int) CudaNdarray_HOST_DIMS(self)[1], // N2
40494043 (unsigned int) CudaNdarray_HOST_DIMS(self)[2], // N3
40504044 (unsigned int) CudaNdarray_HOST_DIMS(self)[3], // N4
40514045 CudaNdarray_DEV_DATA(other), // x
40524046 // x strides
40534047 CudaNdarray_HOST_STRIDES(other)[0],
40544048 CudaNdarray_HOST_STRIDES(other)[1],
40554049 CudaNdarray_HOST_STRIDES(other)[2],
40564050 CudaNdarray_HOST_STRIDES(other)[3],
40574051 CudaNdarray_DEV_DATA(self), // y
40584052 // y strides
40594053 CudaNdarray_HOST_STRIDES(self)[0],
40604054 CudaNdarray_HOST_STRIDES(self)[1],
40614055 CudaNdarray_HOST_STRIDES(self)[2],
40624056 CudaNdarray_HOST_STRIDES(self)[3]
40634057 );
40644058 CNDA_THREAD_SYNC;
40654059 cudaError_t err = cudaGetLastError();
40664060 if( cudaSuccess != err)
40674061 {
40684062 PyErr_Format(PyExc_RuntimeError,
40694063 "Cuda error: %s: %s.",
40704064 "k_copy_4d",
40714065 cudaGetErrorString(err));
40724066 Py_XDECREF(new_other);
40734067 return -1;
40744068 }
40754069 }; break;
40764070 default:
40774071 {
40784072 cudaError_t err = cudaGetLastError();
40794073 if(cudaSuccess != err){
40804074 PyErr_Format(PyExc_RuntimeError,
40814075 "Unexpected Cuda error: %s: %s\n",
40824076 "CudaNdarray_CopyFromCudaNdarray",
40834077 cudaGetErrorString(err));
40844078 Py_XDECREF(new_other);
40854079 return -1;
40864080 }
40874081
40884082 if (verbose)
40894083 fprintf(stderr,
40904084 "Copying with default version unbroadcast=%d\n",
40914085 unbroadcast);
40924086 // call worker routine
40934087 unsigned int threads_per_block = std::min(size,
40944088 (unsigned int)NUM_VECTOR_OP_THREADS_PER_BLOCK);
40954089 unsigned int n_blocks = std::min(ceil_intdiv(size, threads_per_block),
40964090 (unsigned int)NUM_VECTOR_OP_BLOCKS);
40974091 const CudaNdarray * cuda_dims = other;
40984092 if(unbroadcast)
40994093 cuda_dims = self;
41004094 //copy from other into self
41014095 k_elemwise_unary_rowmajor_copy<<<n_blocks, threads_per_block>>>(
41024096 size,
41034097 (unsigned int)other->nd,
41044098 (const int *)CudaNdarray_DEV_DIMS(cuda_dims),
41054099 (const float*)CudaNdarray_DEV_DATA(other),
41064100 (const int *)CudaNdarray_DEV_STRIDES(other),
41074101 CudaNdarray_DEV_DATA(self),
41084102 (const int *)CudaNdarray_DEV_STRIDES(self));
41094103 CNDA_THREAD_SYNC;
41104104 err = cudaGetLastError();
41114105 if(verbose>1)
41124106 fprintf(stderr,
41134107 "INFO k_elemwise_unary_rowmaj (n_blocks=%i,"
41144108 " n_threads_per_block=%i)\n",
41154109 n_blocks, threads_per_block);
41164110 if( cudaSuccess != err)
41174111 {
41184112 //fprint_CudaNdarray(stderr, self);
41194113 //fprint_CudaNdarray(stderr, other);
41204114 PyErr_Format(PyExc_RuntimeError,
41214115 "Cuda error: %s: %s. (n_blocks=%i,"
41224116 " n_threads_per_block=%i)\n",
41234117 "k_elemwise_unary_rowmajor_copy",
41244118 cudaGetErrorString(err), n_blocks,
41254119 threads_per_block);
41264120 Py_XDECREF(new_other);
41274121 return -1;
41284122 }
41294123 }
41304124 };
41314125 Py_XDECREF(new_other);
41324126 return 0;
41334127 }
41344128
41354129 int CudaNdarray_gemm(float alpha, const CudaNdarray * A, const CudaNdarray * B, float beta, CudaNdarray * C)
41364130 {
41374131 if (A->nd != 2)
41384132 {
41394133 PyErr_SetString(PyExc_ValueError, "non-matrix arg A to gemm");
41404134 return -1;
41414135 }
41424136 if (B->nd != 2)
41434137 {
41444138 PyErr_SetString(PyExc_ValueError, "non-matrix arg B to gemm");
41454139 return -1;
41464140 }
41474141 if (C->nd != 2)
41484142 {
41494143 PyErr_SetString(PyExc_ValueError, "non-matrix arg C to gemm");
41504144 return -1;
41514145 }
41524146
41534147 // We must allow dimensions to be zeros.
41544148 if ((CudaNdarray_HOST_DIMS(A)[1] != CudaNdarray_HOST_DIMS(B)[0])
41554149 || (CudaNdarray_HOST_DIMS(A)[0] != CudaNdarray_HOST_DIMS(C)[0])
41564150 || (CudaNdarray_HOST_DIMS(B)[1] != CudaNdarray_HOST_DIMS(C)[1]))
41574151 {
41584152 PyErr_Format(PyExc_ValueError, "dimension mismatch in args to gemm (%i,%i)x(%i,%i)->(%i,%i)",
41594153 CudaNdarray_HOST_DIMS(A)[0],
41604154 CudaNdarray_HOST_DIMS(A)[1],
41614155 CudaNdarray_HOST_DIMS(B)[0],
41624156 CudaNdarray_HOST_DIMS(B)[1],
41634157 CudaNdarray_HOST_DIMS(C)[0],
41644158 CudaNdarray_HOST_DIMS(C)[1]);
41654159 return -1;
41664160 }
41674161
41684162 // If matrix A or B has non-unit size and non-unit stride in both
41694163 // dimensions, we can make a copy.
41704164 CudaNdarray * A_new = NULL;
41714165 CudaNdarray * B_new = NULL;
41724166 if (((CudaNdarray_HOST_DIMS(A)[0] > 1)
41734167 && (CudaNdarray_HOST_STRIDES(A)[0] != 1)
41744168 && (CudaNdarray_HOST_DIMS(A)[1] > 1)
41754169 && (CudaNdarray_HOST_STRIDES(A)[1] != 1))
41764170 || (CudaNdarray_HOST_STRIDES(A)[0] < 0)
41774171 || (CudaNdarray_HOST_STRIDES(A)[1] < 0))
41784172 {
41794173 A_new = (CudaNdarray*) CudaNdarray_Copy(A);
41804174 if (!A_new)
41814175 return -1;
41824176 A = A_new;
41834177 }
41844178
41854179 if (((CudaNdarray_HOST_DIMS(B)[0] > 1)
41864180 && (CudaNdarray_HOST_STRIDES(B)[0] != 1)
41874181 && (CudaNdarray_HOST_DIMS(B)[1] > 1)
41884182 && (CudaNdarray_HOST_STRIDES(B)[1] != 1))
41894183 || (CudaNdarray_HOST_STRIDES(B)[0] < 0)
41904184 || (CudaNdarray_HOST_STRIDES(B)[1] < 0))
41914185 {
41924186 B_new = (CudaNdarray*) CudaNdarray_Copy(B);
41934187 if (!B_new)
41944188 {
41954189 // If A_new is NULL, meaning A was not copied nothing happens
41964190 Py_XDECREF(A_new);
41974191 return -1;
41984192 }
41994193 B = B_new;
42004194 }
42014195
42024196 // If matrix C has non-unit size and non-unit stride in both
42034197 // dimensions, or negative strides, we can't operate. We cannot copy
42044198 // C either, because the calling code will expect the result to be
42054199 // in the original C container.
42064200 if (((CudaNdarray_HOST_DIMS(C)[0] > 1)
42074201 && (CudaNdarray_HOST_STRIDES(C)[0] != 1)
42084202 && (CudaNdarray_HOST_DIMS(C)[1] > 1)
42094203 && (CudaNdarray_HOST_STRIDES(C)[1] != 1))
42104204 || (CudaNdarray_HOST_STRIDES(C)[0] < 0)
42114205 || (CudaNdarray_HOST_STRIDES(C)[1] < 0))
42124206 {
42134207 PyErr_Format(PyExc_AssertionError,
42144208 "non-unit or negative stride in gemm arg C (%i,%i) of shape (%i,%i)",
42154209 CudaNdarray_HOST_STRIDES(C)[0],
42164210 CudaNdarray_HOST_STRIDES(C)[1],
42174211 CudaNdarray_HOST_DIMS(C)[0],
42184212 CudaNdarray_HOST_DIMS(C)[1]);
42194213 Py_XDECREF(A_new);
42204214 Py_XDECREF(B_new);
42214215 return -1;
42224216 }
42234217
42244218 // the unit integer is divided logically into three fields of 4 bits
42254219 // the lowermost 4 bits encode the stride pattern of the output
42264220 // the next higher 4 bits encode the B variable (or y)
42274221 // the next higher 4 bits encode the C variable (or x)
42284222 //
42294223 // the stride pattern for each input is encoded as 0 for unit stride from col to col (Row major)
42304224 // 1 for unit stride from row to row (Col major)
42314225
42324226 // a stride of 0 implies a dimension of 1 - so we can actually define
42334227 // a stride of 0 as a 'unit' stride because gemm will never use it.
42344228 // If a dimension is 0, its stride will not be used either, so we can
42354229 // consider it a 'unit' stride too.
42364230 int unit = 0;
42374231 if (CudaNdarray_HOST_STRIDES(A)[1] == 1 || CudaNdarray_HOST_DIMS(A)[1] <= 1) {
42384232 unit |= (0x0 << 8);
42394233 } else if (CudaNdarray_HOST_STRIDES(A)[0] == 1 || CudaNdarray_HOST_DIMS(A)[0] <= 1) {
42404234 unit |= (0x1 << 8);
42414235 } else {
42424236 unit |= (0x2 << 8);
42434237 }
42444238 if (CudaNdarray_HOST_STRIDES(B)[1] == 1 || CudaNdarray_HOST_DIMS(B)[1] <= 1) {
42454239 unit |= (0x0 << 4);
42464240 } else if (CudaNdarray_HOST_STRIDES(B)[0] == 1 || CudaNdarray_HOST_DIMS(B)[0] <= 1) {
42474241 unit |= (0x1 << 4);
42484242 } else {
42494243 unit |= (0x2 << 4);
42504244 }
42514245 if (CudaNdarray_HOST_STRIDES(C)[1] == 1 || CudaNdarray_HOST_DIMS(C)[1] <= 1) {
42524246 unit |= (0x0 << 0);
42534247 } else if (CudaNdarray_HOST_STRIDES(C)[0] == 1 || CudaNdarray_HOST_DIMS(C)[0] <= 1) {
42544248 unit |= (0x1 << 0);
42554249 } else {
42564250 unit |= (0x2 << 0);
42574251 }
42584252
42594253 /* create appropriate strides for malformed matrices that are row or column
42604254 * vectors
42614255 */
42624256 int sa_0 = (CudaNdarray_HOST_DIMS(A)[0] > 1) ? CudaNdarray_HOST_STRIDES(A)[0] : CudaNdarray_HOST_DIMS(A)[1];
42634257 int sa_1 = (CudaNdarray_HOST_DIMS(A)[1] > 1) ? CudaNdarray_HOST_STRIDES(A)[1] : CudaNdarray_HOST_DIMS(A)[0];
42644258 int sb_0 = (CudaNdarray_HOST_DIMS(B)[0] > 1) ? CudaNdarray_HOST_STRIDES(B)[0] : CudaNdarray_HOST_DIMS(B)[1];
42654259 int sb_1 = (CudaNdarray_HOST_DIMS(B)[1] > 1) ? CudaNdarray_HOST_STRIDES(B)[1] : CudaNdarray_HOST_DIMS(B)[0];
42664260 int sc_0 = (CudaNdarray_HOST_DIMS(C)[0] > 1) ? CudaNdarray_HOST_STRIDES(C)[0] : CudaNdarray_HOST_DIMS(C)[1];
42674261 int sc_1 = (CudaNdarray_HOST_DIMS(C)[1] > 1) ? CudaNdarray_HOST_STRIDES(C)[1] : CudaNdarray_HOST_DIMS(C)[0];
42684262
42694263 float* a = CudaNdarray_DEV_DATA(A);
42704264 float* b = CudaNdarray_DEV_DATA(B);
42714265 float* c = CudaNdarray_DEV_DATA(C);
42724266 cublasOperation_t N = CUBLAS_OP_N;
42734267 cublasOperation_t T = CUBLAS_OP_T;
42744268 //std::cerr << (unit/256) MOD 16 << (unit / 16) MOD 16 << unit MOD 16<< '\\n';
42754269 // There should be no negative stride at that point
42764270 #define CHK_STRIDE_SGEMM(T0, T1, D0, D1, D2, a, x, sx, y, sy, b, z, sz) \
42774271 if (sx == 0){sx = 1;}\
42784272 if (sy == 0){sy = 1;}\
42794273 if (sz == 0){sz = 1;}\
42804274 if ((sx > 0) && (sy > 0) && (sz > 0)) { \
42814275 err = cublasSgemm(handle, T0, T1, D0, D1, D2, &a, x, sx, y, sy, &b, z, sz); \
42824276 } else { \
42834277 PyErr_SetString(PyExc_AssertionError, "negative stride to sGemm");\
42844278 Py_XDECREF(A_new);\
42854279 Py_XDECREF(B_new);\
42864280 return -1; \
42874281 }
42884282
42894283 cublasStatus_t err;
42904284 switch(unit)
42914285 {
42924286 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;
42934287 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;
42944288 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;
42954289 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;
42964290 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;
42974291 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;
42984292 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;
42994293 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;
43004294 default: PyErr_Format(PyExc_ValueError, "some matrix has no unit stride (unit=%x)", unit);
43014295 return -1;
43024296 };
43034297 CNDA_THREAD_SYNC;
43044298 Py_XDECREF(A_new);
43054299 Py_XDECREF(B_new);
43064300
43074301 if (CUBLAS_STATUS_SUCCESS != err)
43084302 {
43094303 PyErr_Format(PyExc_RuntimeError,
43104304 "cublasSgemm failed (%i) %s\n"
43114305 " unit=%x N=%d, c.dims=[%d %d], a.dim=[%d %d], alpha=%f, beta=%f, a=%p, b=%p, c=%p"
43124306 " sa_0=%d, sa_1=%d, sb_0=%d, sb_1=%d, sc_0=%d, sc_1=%d",
43134307 err, cublasGetErrorString(err),
43144308 unit, N,
43154309 CudaNdarray_HOST_DIMS(C)[0],
43164310 CudaNdarray_HOST_DIMS(C)[1],
43174311 CudaNdarray_HOST_DIMS(A)[0], CudaNdarray_HOST_DIMS(A)[1],
43184312 alpha, beta, a, b, c, sa_0, sa_1, sb_0, sb_1, sc_0, sc_1);
43194313
43204314 return -1;
43214315 }
43224316 return 0;
43234317 }
43244318
43254319 int CudaNdarray_sgemv(float alpha, const CudaNdarray * A, const CudaNdarray * B, float beta, CudaNdarray * C)
43264320 {
43274321 /**
43284322 * C <- alpha A B + beta C
43294323 * A : matrix
43304324 * B, C: vector
43314325 * alpha, beta: scalars
43324326 */
43334327 if (A->nd != 2) { PyErr_SetString(PyExc_ValueError, "non-matrix arg to gemv"); return -1; }
43344328 if (B->nd != 1) { PyErr_SetString(PyExc_ValueError, "non-vector arg to gemv"); return -1; }
43354329 if (C->nd != 1) { PyErr_SetString(PyExc_ValueError, "non-vector arg to gemv"); return -1; }
43364330
43374331 // We must allow dimensions to be zeros.
43384332 if ((CudaNdarray_HOST_DIMS(A)[1] != CudaNdarray_HOST_DIMS(B)[0])
43394333 || (CudaNdarray_HOST_DIMS(A)[0] != CudaNdarray_HOST_DIMS(C)[0]))
43404334 {
43414335 PyErr_Format(PyExc_ValueError, "dimension mismatch in args to gemv (%i,%i)x(%i)->(%i)",
43424336 CudaNdarray_HOST_DIMS(A)[0],
43434337 CudaNdarray_HOST_DIMS(A)[1],
43444338 CudaNdarray_HOST_DIMS(B)[0],
43454339 CudaNdarray_HOST_DIMS(C)[0]);
43464340 return -1;
43474341 }
43484342
43494343 // If matrix A has non-unit size and non-unit stride in both
43504344 // dimensions, or negative strides, we cannot operate, but we can
43514345 // make a copy.
43524346 CudaNdarray * A_new = NULL;
43534347 CudaNdarray * B_new = NULL;
43544348 if (((CudaNdarray_HOST_DIMS(A)[0] > 1)
43554349 && (CudaNdarray_HOST_STRIDES(A)[0] != 1)
43564350 && (CudaNdarray_HOST_DIMS(A)[1] > 1)
43574351 && (CudaNdarray_HOST_STRIDES(A)[1] != 1))
43584352 || (CudaNdarray_HOST_STRIDES(A)[0] < 0)
43594353 || (CudaNdarray_HOST_STRIDES(A)[1] < 0))
43604354 {
43614355 A_new = (CudaNdarray*) CudaNdarray_Copy(A);
43624356 if (!A_new)
43634357 return -1;
43644358 A = A_new;
43654359 }
43664360
43674361 // If vector B as a negative stride, we also have to make a copy.
43684362 if (CudaNdarray_HOST_STRIDES(B)[0] < 0)
43694363 {
43704364 B_new = (CudaNdarray*) CudaNdarray_Copy(B);
43714365 if (!B_new)
43724366 {
43734367 // If A was not copied, A_new is NULL, and Py_XDECREF does not
43744368 // do anything
43754369 Py_XDECREF(A_new);
43764370 return -1;
43774371 }
43784372 B = B_new;
43794373 }
43804374
43814375 // cudablas does not handle negative strides as expected
43824376 if ( (CudaNdarray_HOST_STRIDES(A)[0] < 0)
43834377 || (CudaNdarray_HOST_STRIDES(A)[1] < 0))
43844378 {
43854379 PyErr_Format(PyExc_ValueError, "illegal strides in args to gemv (%i,%i)",
43864380 CudaNdarray_HOST_STRIDES(A)[0],
43874381 CudaNdarray_HOST_STRIDES(A)[1]);
43884382 Py_XDECREF(A_new);
43894383 Py_XDECREF(B_new);
43904384 return -1;
43914385 }
43924386
43934387 /* create appropriate strides for malformed matrices that are row or column
43944388 * vectors
43954389 */
43964390 int sa_0 = (CudaNdarray_HOST_DIMS(A)[0] > 1) ? CudaNdarray_HOST_STRIDES(A)[0] : CudaNdarray_HOST_DIMS(A)[1];
43974391 int sa_1 = (CudaNdarray_HOST_DIMS(A)[1] > 1) ? CudaNdarray_HOST_STRIDES(A)[1] : CudaNdarray_HOST_DIMS(A)[0];
43984392 int sb_0 = (CudaNdarray_HOST_DIMS(B)[0] > 1) ? CudaNdarray_HOST_STRIDES(B)[0] : 1;
43994393 int sc_0 = (CudaNdarray_HOST_DIMS(C)[0] > 1) ? CudaNdarray_HOST_STRIDES(C)[0] : 1;
44004394
44014395 if (sa_0 == 0)
44024396 sa_0 = 1;
44034397 if (sa_1 == 0)
44044398 sa_1 = 1;
44054399
44064400 // This is important because we can end up not calling Sgemv at all
44074401 cublasStatus_t err = CUBLAS_STATUS_SUCCESS;
44084402 if (CudaNdarray_SIZE(C)) {
44094403 if ((CudaNdarray_HOST_DIMS(A)[0] <= 1)
44104404 || ((CudaNdarray_HOST_STRIDES(A)[0] == 1)
44114405 && (CudaNdarray_HOST_STRIDES(A)[1] > 0)))
44124406 {
44134407 err = cublasSgemv(handle, CUBLAS_OP_N,
44144408 CudaNdarray_HOST_DIMS(A)[0], CudaNdarray_HOST_DIMS(A)[1],
44154409 &alpha,
44164410 CudaNdarray_DEV_DATA(A), sa_1,
44174411 CudaNdarray_DEV_DATA(B), sb_0,
44184412 &beta,
44194413 CudaNdarray_DEV_DATA(C), sc_0);
44204414 }
44214415 else if ((CudaNdarray_HOST_DIMS(A)[1] <= 1)
44224416 || ((CudaNdarray_HOST_STRIDES(A)[1] == 1)
44234417 && (CudaNdarray_HOST_STRIDES(A)[0] > 0)))
44244418 {
44254419 err = cublasSgemv(handle, CUBLAS_OP_T,
44264420 CudaNdarray_HOST_DIMS(A)[1], CudaNdarray_HOST_DIMS(A)[0],
44274421 &alpha,
44284422 CudaNdarray_DEV_DATA(A), sa_0,
44294423 CudaNdarray_DEV_DATA(B), sb_0,
44304424 &beta,
44314425 CudaNdarray_DEV_DATA(C), sc_0);
44324426 }
44334427 else
44344428 {
44354429 PyErr_Format(PyExc_AssertionError,
44364430 "Unexpected stride pattern in gemv: (%i, %i) x %i -> %i.\n"
44374431 "Shapes are: (%i, %i) x %i -> %i\n",
44384432 CudaNdarray_HOST_STRIDES(A)[0],
44394433 CudaNdarray_HOST_STRIDES(A)[1],
44404434 CudaNdarray_HOST_STRIDES(B)[0],
44414435 CudaNdarray_HOST_STRIDES(C)[0],
44424436 CudaNdarray_HOST_DIMS(A)[0],
44434437 CudaNdarray_HOST_DIMS(A)[1],
44444438 CudaNdarray_HOST_DIMS(B)[0],
44454439 CudaNdarray_HOST_DIMS(C)[0]);
44464440 Py_XDECREF(A_new);
44474441 Py_XDECREF(B_new);
44484442 return -1;
44494443 }
44504444 }
44514445
44524446 CNDA_THREAD_SYNC;
44534447 Py_XDECREF(A_new);
44544448 Py_XDECREF(B_new);
44554449
44564450 if (CUBLAS_STATUS_SUCCESS != err)
44574451 {
44584452 PyErr_Format(PyExc_RuntimeError,
44594453 "cublasSgemv failed (%i)",
44604454 err);
44614455 return -1;
44624456 }
44634457 return 0;
44644458 }
44654459
44664460 int CudaNdarray_sger(float alpha, const CudaNdarray * x, const CudaNdarray * y, CudaNdarray * A) {
44674461 if (x->nd != 1) { PyErr_SetString(PyExc_ValueError, "non-vector arg x to sger"); return -1; }
44684462 if (y->nd != 1) { PyErr_SetString(PyExc_ValueError, "non-vector arg y to sger"); return -1; }
44694463 if (A->nd != 2) { PyErr_SetString(PyExc_ValueError, "non-matrix arg A to sger"); return -1; }
44704464
44714465 if ((CudaNdarray_HOST_DIMS(A)[0] != CudaNdarray_HOST_DIMS(x)[0])
44724466 || (CudaNdarray_HOST_DIMS(A)[1] != CudaNdarray_HOST_DIMS(y)[0])) {
44734467 PyErr_Format(PyExc_ValueError,
44744468 "dimension mismatch in args to sger (%i)x(%i)->(%i,%i)",
44754469 CudaNdarray_HOST_DIMS(x)[0],
44764470 CudaNdarray_HOST_DIMS(y)[0],
44774471 CudaNdarray_HOST_DIMS(A)[0],
44784472 CudaNdarray_HOST_DIMS(A)[1]);
44794473 return -1;
44804474 }
44814475
44824476 int x_strides = CudaNdarray_HOST_STRIDES(x)[0];
44834477 CudaNdarray * x_new = NULL;
44844478 if(x_strides == 0){
44854479 if(CudaNdarray_HOST_DIMS(x)[0] != 1){
44864480 PyErr_Format(PyExc_RuntimeError,
44874481 "CudaNdarray_sger: Invalid input x (should not happen)."
44884482 " We received a CudaNdarray vector with a stride of 0"
44894483 " that has more than 1 element!");
44904484 return -1;
44914485 }
44924486 x_strides = 1;
44934487 } else if(x_strides < 0){
44944488 x_new = (CudaNdarray*) CudaNdarray_Copy(x);
44954489 x = x_new;
44964490 x_strides = CudaNdarray_HOST_STRIDES(x)[0];
44974491 }
44984492
44994493 int y_strides = CudaNdarray_HOST_STRIDES(y)[0];
45004494 CudaNdarray * y_new = NULL;
45014495 if(y_strides == 0){
45024496 if(CudaNdarray_HOST_DIMS(y)[0] != 1){
45034497 PyErr_Format(PyExc_RuntimeError,
45044498 "CudaNdarray_sger: Invalid input y (should not happen)."
45054499 " We received a CudaNdarray vector with a stride of 0"
45064500 " that has more than 1 elements!");
45074501 Py_XDECREF(x_new);
45084502 return -1;
45094503 }
45104504 y_strides = 1;
45114505 } else if(y_strides < 0){
45124506 y_new = (CudaNdarray*) CudaNdarray_Copy(y);
45134507 y = y_new;
45144508 y_strides = CudaNdarray_HOST_STRIDES(y)[0];
45154509 }
45164510
45174511 // Create appropriate strides if A is a row or column vector
45184512 int sa_0 = (CudaNdarray_HOST_DIMS(A)[0] > 1) ? CudaNdarray_HOST_STRIDES(A)[0]
45194513 : CudaNdarray_HOST_DIMS(A)[1];
45204514 int sa_1 = (CudaNdarray_HOST_DIMS(A)[1] > 1) ? CudaNdarray_HOST_STRIDES(A)[1]
45214515 : CudaNdarray_HOST_DIMS(A)[0];
45224516
45234517 // This is important because we can end up not calling Sger at all
45244518 cublasStatus_t err = CUBLAS_STATUS_SUCCESS;
45254519 if(CudaNdarray_SIZE(A)){
45264520 // If A is in col-major
45274521 if ((CudaNdarray_HOST_DIMS(A)[0] <= 1)
45284522 || ((CudaNdarray_HOST_STRIDES(A)[0] == 1)
45294523 && (CudaNdarray_HOST_STRIDES(A)[1] > 0)))
45304524 {
45314525 err = cublasSger(handle, CudaNdarray_HOST_DIMS(x)[0], CudaNdarray_HOST_DIMS(y)[0], &alpha,
45324526 CudaNdarray_DEV_DATA(x), x_strides,
45334527 CudaNdarray_DEV_DATA(y), y_strides,
45344528 CudaNdarray_DEV_DATA(A), sa_1);
45354529 }
45364530 // Since Sger expects A in col-major, we invert x and y to fake this.
45374531 else if ((CudaNdarray_HOST_DIMS(A)[1] <= 1)
45384532 || ((CudaNdarray_HOST_STRIDES(A)[1] == 1)
45394533 && (CudaNdarray_HOST_STRIDES(A)[0] > 0)))
45404534 {
45414535 err = cublasSger(handle, CudaNdarray_HOST_DIMS(y)[0], CudaNdarray_HOST_DIMS(x)[0], &alpha,
45424536 CudaNdarray_DEV_DATA(y), y_strides,
45434537 CudaNdarray_DEV_DATA(x), x_strides,
45444538 CudaNdarray_DEV_DATA(A), sa_0);
45454539 }
45464540 // A has to be either c- or f-contiguous, with no negative strides
45474541 else
45484542 {
45494543 PyErr_SetString(PyExc_NotImplementedError,
45504544 "non-contiguous A, or negative strides, in sger");
45514545 Py_XDECREF(x_new);
45524546 Py_XDECREF(y_new);
45534547 return -1;
45544548 }
45554549 }
45564550 CNDA_THREAD_SYNC;
45574551 Py_XDECREF(x_new);
45584552 Py_XDECREF(y_new);
45594553
45604554 if (CUBLAS_STATUS_SUCCESS != err)
45614555 {
45624556 PyErr_Format(PyExc_RuntimeError,
45634557 "cublasSger failed (%i)",
45644558 err);
45654559 return -1;
45664560 }
45674561
45684562 return 0;
45694563 }
45704564
45714565 /**
45724566 *
45734567 * Precondition:
45744568 * a->dim[d] == (dims_a[d]==0) ? (1 << log2_dims_a[d]) : dims_a[d]
45754569 * z->dim[d] == (z_str[d]==0) ? 1 : dims_a[d];
45764570 *
45774571 * TODO: templatize this function to support other reductions.
45784572 * All that needs to change is the initial value for sum, and the reduction operator.
45794573 */
45804574
45814575 static __global__ void kernel_reduce_sum(const unsigned int size_z,
45824576 const unsigned int nd,
45834577 const int * dims_a,
45844578 const int * log2_dims_a,
45854579 const int * a_str,
45864580 const float * a_data,
45874581 const int * z_str,
45884582 float * z_data)
45894583 {
45904584 const unsigned int idx = blockIdx.x * blockDim.x + threadIdx.x;
45914585 const unsigned int numThreads = blockDim.x * gridDim.x;
45924586
45934587 //structure data contains the strides and dimensions of both a and z
45944588 // a_dim[0], a_dim[1], ... a_dim[nd-1],
45954589 // a_log2dim[0], a_log2dim[1], ... a_log2dim[nd-1],
45964590 // a_str[0], ... a_str[nd-1],
45974591 // z_str[0], ... z_str[nd-1]
45984592 extern __shared__ int structure_data[];
45994593 for (unsigned int i = threadIdx.x; i < nd; i += blockDim.x)
46004594 {
46014595 structure_data[i+0*nd] = dims_a[i];
46024596 structure_data[i+1*nd] = log2_dims_a[i];
46034597 structure_data[i+2*nd] = a_str[i];
46044598 structure_data[i+3*nd] = z_str[i];
46054599 }
46064600 dims_a = structure_data;
46074601 log2_dims_a = structure_data + nd;
46084602 a_str = structure_data + 2*nd;
46094603 z_str = structure_data + 3*nd;
46104604
46114605 __syncthreads(); //wait for all the shared structure to be loaded
46124606
46134607 for (unsigned int i = idx; i < size_z; i += numThreads)
46144608 {
46154609 unsigned int ii = i;
46164610 const float * a_data_i = a_data;
46174611 float * z_data_i = z_data;
46184612 unsigned int n_reduce_elements = 1;
46194613 unsigned int n_reduce_dims = 0;
46204614 unsigned int reduce_dim0 = nd-1;
46214615
46224616
46234617 //In this loop, we locate the initial element of the slice that we'd like to reduce with this thread
46244618 // At the same time, we [re]calculate the size of that slice (n_reduce_elements)
46254619 for (unsigned int d = 0; d < nd; ++d)
46264620 {
46274621 if (a_str[d] && (!z_str[d])) // this means 'd' is a dimension we are reducing over
46284622 {
46294623 n_reduce_elements *= dims_a[d];
46304624 n_reduce_dims += 1;
46314625 reduce_dim0 = (d < reduce_dim0) ? d : reduce_dim0;
46324626 }
46334627 else //'d' is not a dimension that we are reducing over
46344628 {
46354629 unsigned int pos_d;
46364630 if (log2_dims_a[d]==-1) //TODO: when things are working, use this switch
46374631 {
46384632 // this branch is not preferred,
46394633 // because the manual said that integer mod and div operations are slow on gpu
46404634 pos_d = (ii % dims_a[d]);
46414635 ii = (ii / dims_a[d]);
46424636 }
46434637 else
46444638 {
46454639 pos_d = (ii & ((1 << log2_dims_a[d])-1)); //take the lower log2_dims bits
46464640 ii = (ii >> log2_dims_a[d]); //shift those lower log2_dims bits off of ii
46474641 }
46484642 a_data_i += pos_d * a_str[d];
46494643 z_data_i += pos_d * z_str[d];
46504644 }
46514645 }
46524646 // now we've got pointers a_data_i and z_data_i into element 0 of the slice over which we are reducing
46534647 // do a similar loop
46544648
46554649 float sum = 0.0f;
46564650 switch(n_reduce_dims)
46574651 {
46584652 case 0:
46594653 {
46604654 sum = a_data_i[0];
46614655 }
46624656 break;
46634657 case 1:
46644658 {
46654659 const int stride = a_str[reduce_dim0];
46664660 const float * a_data_i_max = a_data_i + dims_a[reduce_dim0] * stride;
46674661 while (a_data_i != a_data_i_max)
46684662 {
46694663 sum += a_data_i[0];
46704664 a_data_i += stride;
46714665 }
46724666 }
46734667 break;
46744668 case 2:
46754669 {
46764670 int rd = reduce_dim0+1;
46774671 for (; rd < nd; ++rd)
46784672 {
46794673 if (a_str[rd] && (!z_str[rd])) // this means 'rd' is a dimension we are reducing over
46804674 break;
46814675 }
46824676 const int stride0 = a_str[reduce_dim0];
46834677 const int stride1 = a_str[rd];
46844678 for (int ii = 0; ii < dims_a[rd]; ++ii)
46854679 {
46864680 const float * a_data_ri = a_data_i + ii * stride1;
46874681 const float * a_data_ri_max = a_data_ri + dims_a[reduce_dim0] * stride0;
46884682 while (a_data_ri != a_data_ri_max)
46894683 {
46904684 sum += a_data_ri[0];
46914685 a_data_ri += stride0;
46924686 }
46934687 }
46944688 };
46954689 break;
46964690 default:
46974691 {
46984692 for (unsigned int reduce_i = 0; reduce_i < n_reduce_elements; ++reduce_i)
46994693 {
47004694 //TODO: optimize this loop to work more like theano's Elemwise. It's serial code.
47014695 unsigned int reduce_ii = reduce_i;
47024696 const float * a_data_ri = a_data_i;
47034697
47044698 //This loop finds the element in the a slice to add.
47054699 for (unsigned int rd = reduce_dim0; rd < nd; ++rd)
47064700 {
47074701 unsigned int pos_d;
47084702 if (a_str[rd] && (!z_str[rd])) // this means 'd' is a dimension we are reducing over
47094703 {
47104704 if (log2_dims_a[rd]==-1)
47114705 {
47124706 // this branch is not preferred,
47134707 // because the manual said that integer mod and div operations are slow on gpu
47144708 pos_d = (reduce_ii % dims_a[rd]);
47154709 reduce_ii = (reduce_ii / dims_a[rd]);
47164710 }
47174711 else
47184712 {
47194713 pos_d = (reduce_ii & ((1 << log2_dims_a[rd])-1)); //take the lower log2_dims bits
47204714 reduce_ii = (reduce_ii >> log2_dims_a[rd]); //shift those lower log2_dims bits off of ii
47214715 }
47224716 a_data_ri += pos_d * a_str[rd];
47234717 }
47244718 }
47254719 sum += a_data_ri[0];
47264720 }
47274721 }
47284722 }
47294723 z_data_i[0] = sum;
47304724 }
47314725 }
47324726
47334727 static __global__ void kernel_reduce_sum_1011(
47344728 const unsigned int d0,
47354729 const unsigned int d1,
47364730 const unsigned int d2,
47374731 const unsigned int d3,
47384732 const float *A, const int sA0, const int sA1, const int sA2, const int sA3,
47394733 float * Z, const int sZ0)
47404734 {
47414735 const int threadCount = blockDim.x * blockDim.y * blockDim.z;
47424736 const int threadNum = threadIdx.z * blockDim.x * blockDim.y + threadIdx.y * blockDim.x + threadIdx.x;
47434737 extern __shared__ float buf[];
47444738 float mysum = 0.0f;
47454739
47464740 if (warpSize != 32)
47474741 {
47484742 return; //TODO: set error code
47494743 }
47504744
47514745 for (int i0 = threadIdx.z; i0 < d0; i0 += blockDim.z)
47524746 {
47534747 float Ai = A[i0 * sA0 + blockIdx.x * sA1 + threadIdx.y * sA2 + threadIdx.x * sA3];
47544748 mysum += Ai;
47554749 }
47564750 buf[threadNum] = mysum;
47574751 __syncthreads();
47584752
47594753 // rest of function is handled by one warp
47604754 if (threadNum < warpSize)
47614755 {
47624756 for (int i = threadNum + warpSize; i < threadCount; i += warpSize)
47634757 {
47644758 mysum += buf[i];
47654759 }
47664760 buf[threadNum] = mysum;
47674761 if (threadNum < 16)
47684762 {
47694763 //reduce so that threadNum 0 has the sum of everything
47704764 if(threadNum + 16 < threadCount) buf[threadNum] += buf[threadNum+16];
47714765 if(threadNum + 8 < threadCount) buf[threadNum] += buf[threadNum+8];
47724766 if(threadNum + 4 < threadCount) buf[threadNum] += buf[threadNum+4];
47734767 if(threadNum + 2 < threadCount) buf[threadNum] += buf[threadNum+2];
47744768 if(threadNum + 1 < threadCount) buf[threadNum] += buf[threadNum+1];
47754769 if (threadNum == 0)
47764770 {
47774771 Z[blockIdx.x*sZ0] = buf[0];
47784772 }
47794773 }
47804774 }
47814775 }
47824776 /**
47834777 * Dimensions in which the self has size 1 and A has size > 1 are considered summing dimensions
47844778 * 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.
47854779 */
47864780 int
47874781 CudaNdarray_reduce_sum(CudaNdarray * self, CudaNdarray * A)
47884782 {
47894783 int verbose = 0;
47904784 //check input rank
47914785 if (self->nd != A->nd)
47924786 {
47934787 PyErr_Format(PyExc_TypeError, "Rank mismatch in CudaNdarray_sum: %i vs %i", self->nd, A->nd);
47944788 return -1;
47954789 }
47964790 for (int i = 0; i < self->nd; ++i)
47974791 {
47984792 if ((CudaNdarray_HOST_DIMS(self)[i] > 1) && (CudaNdarray_HOST_DIMS(self)[i] != CudaNdarray_HOST_DIMS(A)[i]))
47994793 {
48004794 PyErr_Format(PyExc_TypeError, "Dimension mismatch in CudaNdarray_sum: self->dim[%i] == %i , A->dim[%i] = %i",
48014795 i, CudaNdarray_HOST_DIMS(self)[i], i, CudaNdarray_HOST_DIMS(A)[i]);
48024796 return -1;
48034797 }
48044798 }
48054799
48064800 int n_summations = (unsigned int)CudaNdarray_SIZE(self);
48074801 if (verbose)
48084802 {
48094803 std::cerr << "reduce_sum n_summations " << n_summations << '\n';
48104804 std::cerr << "reduce_sum nd " << self->nd << '\n';
48114805 fprint_CudaNdarray(stderr, A);
48124806 fprint_CudaNdarray(stderr, self);
48134807 }
48144808 if (0 && (A->nd == 4) //check to see if kernel_reduce_sum_1011 applies
48154809 && (CudaNdarray_HOST_DIMS(self)[0] == 1)
48164810 && (CudaNdarray_HOST_DIMS(self)[2] == 1)
48174811 && (CudaNdarray_HOST_DIMS(self)[3] == 1)
48184812 )
48194813 {
48204814 dim3 n_threads(CudaNdarray_HOST_DIMS(A)[3], CudaNdarray_HOST_DIMS(A)[2]);
48214815 dim3 n_blocks(CudaNdarray_HOST_DIMS(A)[1]);
48224816 while (n_threads.x * n_threads.y * n_threads.z < NUM_VECTOR_OP_THREADS_PER_BLOCK) ++n_threads.z;
48234817 n_threads.z -= 1;
48244818 if (n_threads.z > 64) n_threads.z = 64;
48254819 if (n_threads.z)
48264820 {
48274821 if (verbose) printf("trying kernel_reduce_sum_1011\n");
48284822 int n_shared = sizeof(float) * n_threads.x * n_threads.y * n_threads.z;
48294823 kernel_reduce_sum_1011<<<n_blocks, n_threads, n_shared>>>(
48304824 CudaNdarray_HOST_DIMS(A)[0],
48314825 CudaNdarray_HOST_DIMS(A)[1],
48324826 CudaNdarray_HOST_DIMS(A)[2],
48334827 CudaNdarray_HOST_DIMS(A)[3],
48344828 CudaNdarray_DEV_DATA(A),
48354829 CudaNdarray_HOST_STRIDES(A)[0],
48364830 CudaNdarray_HOST_STRIDES(A)[1],
48374831 CudaNdarray_HOST_STRIDES(A)[2],
48384832 CudaNdarray_HOST_STRIDES(A)[3],
48394833 CudaNdarray_DEV_DATA(self),
48404834 CudaNdarray_HOST_STRIDES(self)[1]);
48414835 CNDA_THREAD_SYNC;
48424836 if (cudaSuccess == cudaGetLastError()) return 0;
48434837 if (verbose) printf("failed, falling back to kernel_reduce_sum\n");
48444838 }
48454839 }
48464840
48474841 int n_threads_per_block = std::min(n_summations,
48484842 NUM_VECTOR_OP_THREADS_PER_BLOCK);
48494843 int n_blocks = std::min(ceil_intdiv(n_summations,n_threads_per_block),
48504844 NUM_VECTOR_OP_BLOCKS);
48514845 int n_structure_cache = self->nd * 4 * sizeof(int);
48524846
48534847 if (verbose)
48544848 {
48554849 std::cerr << "n_blocks, n_threads_per_block " << n_blocks << ' ' << n_threads_per_block << '\n';
48564850 }
48574851 assert (self->nd > 0);
48584852 assert (self->nd == A->nd);
48594853 kernel_reduce_sum<<<n_blocks, n_threads_per_block, n_structure_cache>>>(
48604854 n_summations,
48614855 self->nd,
48624856 CudaNdarray_DEV_DIMS(A),
48634857 CudaNdarray_DEV_LOG2DIMS(A),
48644858 CudaNdarray_DEV_STRIDES(A),
48654859 CudaNdarray_DEV_DATA(A),
48664860 CudaNdarray_DEV_STRIDES(self),
48674861 CudaNdarray_DEV_DATA(self));
48684862 CNDA_THREAD_SYNC;
48694863 cudaError_t err = cudaGetLastError();
48704864 if (cudaSuccess != err)
48714865 {
48724866 PyErr_Format(PyExc_RuntimeError, "Cuda error: %s: %s.\n", "kernel_reduce_sum", cudaGetErrorString(err));
48734867 return -1;
48744868 }
48754869 return 0;
48764870 }
48774871 int
48784872 CudaNdarray_reduce_prod(CudaNdarray * self, const CudaNdarray * A)
48794873 {
48804874 PyErr_SetString(PyExc_NotImplementedError, "");
48814875 return -1;
48824876 }
48834877 int
48844878 CudaNdarray_reduce_min(CudaNdarray * self, const CudaNdarray * A)
48854879 {
48864880 PyErr_SetString(PyExc_NotImplementedError, "");
48874881 return -1;
48884882 }
48894883 int
48904884 CudaNdarray_reduce_max(CudaNdarray * self, const CudaNdarray * A)
48914885 {
48924886 PyErr_SetString(PyExc_NotImplementedError, "");
48934887 return -1;
48944888 }
48954889
48964890
48974891 /**
48984892 *
48994893 * pattern is a permutation of [0, 1, ... self->nd-1] with the following twists:
49004894 * - an element 'd' of the permutation can be dropped if CudaNdarray_HOST_DIMS(self)[d] == 1
49014895 * - any number of '-1' elements can be in the pattern, and they will cause new ranks (with dim==1) to be inserted.
49024896 *
49034897 * 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:
49044898 * [4, 6, 1, 1, 5] (we dropped the original dim[2]==1, and inserted two singleton dimensions with the -1s.
49054899 */
49064900 int
49074901 CudaNdarray_dimshuffle(CudaNdarray * self, unsigned int len, const int * pattern)
49084902 {
49094903 //TODO: pass a workspace pointer to avoid the internal malloc
49104904 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.
49114905 int * newstrides = newdims + len;
49124906 int * dims_taken = newstrides + len;
49134907 if (!newdims)
49144908 {
49154909 PyErr_SetString(PyExc_MemoryError, "CudaNdarray_dimshuffle: Failed to allocate temporary space");
49164910 return -1;
49174911 }
49184912 for (int i = 0; i < self->nd; ++i)
49194913 {
49204914 dims_taken[i] = 0;
49214915 }
49224916 for (int i = 0; i < len; ++i)
49234917 {
49244918 if (pattern[i] < 0)
49254919 {
49264920 newdims[i] = 1;
49274921 newstrides[i] = 0;
49284922 }
49294923 else if(dims_taken[pattern[i]])
49304924 {
49314925 PyErr_Format(PyExc_ValueError, "Cudandarray_dimshuffle: invalid pattern for Cudandarray_dimshuffle. You used the dimensions %d multiple time",
49324926 pattern[i]);
49334927 free(newdims);
49344928 return -1;
49354929 }
49364930 else if (pattern[i]>= self->nd)
49374931 {
49384932 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",
49394933 pattern[i], self->nd);
49404934 free(newdims);
49414935 return -1;
49424936 }
49434937 else
49444938 {
49454939 newdims[i] = CudaNdarray_HOST_DIMS(self)[pattern[i]];
49464940 newstrides[i] = CudaNdarray_HOST_STRIDES(self)[pattern[i]];
49474941 dims_taken[pattern[i]] = 1;
49484942 }
49494943 }
49504944 //Check if we dropped not broadcastable dims
49514945 for (int i = 0; i < self->nd; ++i)
49524946 {
49534947 if (dims_taken[i]==0 && CudaNdarray_HOST_DIMS(self)[i]!=1)
49544948 {
49554949 PyErr_SetString(PyExc_ValueError, "Cudandarray_dimshuffle: You cannot drop a non-broadcastable dimension.");
49564950 free(newdims);
49574951 return -1;
49584952 }
49594953 }
49604954 //swap this structure in for the one in self, and sync to the card
49614955 if (CudaNdarray_set_nd(self, len))
49624956 {
49634957 free(newdims);
49644958 return -1;
49654959 }
49664960 for (int i = 0; i < len; ++i)
49674961 {
49684962 CudaNdarray_set_dim(self, i, newdims[i]);
49694963 CudaNdarray_set_stride(self, i, newstrides[i]);
49704964 }
49714965 if (cnda_copy_structure_to_device(self))
49724966 {
49734967 free(newdims);
49744968 return -1;
49754969 }
49764970 free(newdims);
49774971 return 0;
49784972 }
49794973
49804974
49814975
49824976 /**
49834977 *
49844978 * This is the function that bind to python.
49854979 * See CudaNdarray_dimshuffle to call from C.
49864980 * We use -1 to mean 'x' as in Tensor Dimshuffle.
49874981 */
49884982 PyObject *
49894983 CudaNdarray_Dimshuffle(PyObject* _unused, PyObject* args)
49904984 {
49914985 PyObject * self = NULL;
49924986 PyObject * pattern_object = NULL;
49934987 int * pattern = NULL;
49944988 PyObject * rval = NULL;
49954989 int success = -1;
49964990 //const int * dims = NULL;
49974991
49984992 //args should consist of two python objects ("OO")
49994993 if (! PyArg_ParseTuple(args, "OO", &self, &pattern_object))
50004994 return NULL;
50014995
50024996 if (!CudaNdarray_Check(self) )
50034997 {
50044998 PyErr_SetString(PyExc_TypeError, "First argument to cuda_ndarray.dimshuffle must be a CudaNdarray");
50054999 return NULL;
50065000 }
50075001
50085002 //parse pattern_object into int * pattern
50095003
50105004 Py_ssize_t pattern_dim = PyObject_Length(pattern_object);
50115005
50125006 if (pattern_dim < 0)
50135007 {
50145008 PyErr_SetString(PyExc_TypeError, "Couldn't get length of third argument to cuda_ndarray.dimshuffle");
50155009 return NULL;
50165010 }
50175011
50185012 pattern = (int *) malloc( pattern_dim * sizeof(int));
50195013
50205014 for (Py_ssize_t i = 0; i < pattern_dim; i++)
50215015 {
50225016 PyObject * idx = PyLong_FromLong(i);
50235017
50245018 if (idx == NULL)
50255019 {
50265020 PyErr_SetString(PyExc_Exception, "Couldn't make long object to loop over list/tuple");
50275021 goto CudaNdarray_dimshuffle_fail;
50285022 }
50295023
50305024 long elem_value = 0;
50315025
50325026 PyObject * elem = PyObject_GetItem(pattern_object, idx);
50335027
50345028 if (elem == NULL)
50355029 {
50365030 Py_XDECREF( elem);
50375031 PyErr_SetString(PyExc_ValueError, "Third argument to dimshuffle must be list or tuple of integers");
50385032 goto CudaNdarray_dimshuffle_fail;
50395033 }
50405034
50415035 elem_value = PyInt_AsLong(elem);
50425036
50435037 if (elem_value == -1 && PyErr_Occurred() )
50445038 {
50455039 Py_XDECREF(elem);
50465040 PyErr_SetString(PyExc_ValueError, "Third argument to dimshuffle must be list or tuple of integers");
50475041 goto CudaNdarray_dimshuffle_fail;
50485042 }
50495043
50505044 pattern[i] = elem_value;
50515045
50525046 Py_XDECREF( elem );
50535047 Py_XDECREF( idx );
50545048 }
50555049
50565050 //allocate rval
50575051 rval = (PyObject *) CudaNdarray_View((CudaNdarray *) self);
50585052
50595053 if (rval == NULL)
50605054 {
50615055 //CudaNdarray_New should have set the exception string
50625056 goto CudaNdarray_dimshuffle_fail;
50635057 }
50645058
50655059
50665060 //printf("pattern_dim: %d\n",pattern_dim);
50675061 //printf("pattern: %d %d\n",pattern[0],pattern[1]);
50685062 //dims = CudaNdarray_HOST_DIMS( (CudaNdarray *) self);
50695063 //printf("dims before: %d %d\n",dims[0],dims[1]);
50705064
50715065 success = CudaNdarray_dimshuffle((CudaNdarray *) rval, pattern_dim, pattern);
50725066
50735067 if (success != 0)
50745068 {
50755069 //Exception string should already be set by CudaNdarray_dimshuffle
50765070 goto CudaNdarray_dimshuffle_fail;
50775071 }
50785072
50795073 free(pattern);
50805074
50815075 return rval;
50825076
50835077 CudaNdarray_dimshuffle_fail:
50845078
50855079 if (pattern != NULL)
50865080 free(pattern);
50875081
50885082 Py_XDECREF(rval);
50895083 return NULL;
50905084 }
50915085
50925086
50935087 int
50945088 cnda_structure_size(int nd)
50955089 {
50965090 // dim0, dim1, ...
50975091 // str0, str1, ...
50985092 // log2(dim0), log2(dim1), ...
50995093 return nd + nd + nd;
51005094 }
51015095
51025096 const int *
51035097 CudaNdarray_HOST_DIMS(const CudaNdarray * self)
51045098 {
51055099 return self->host_structure;
51065100 }
51075101
51085102 const int *
51095103 CudaNdarray_HOST_STRIDES(const CudaNdarray * self)
51105104 {
51115105 return self->host_structure + self->nd;
51125106 }
51135107 const int *
51145108 CudaNdarray_HOST_LOG2DIMS(const CudaNdarray * self)
51155109 {
51165110 return self->host_structure + 2*self->nd;
51175111 }
51185112
51195113 int
51205114 CudaNdarray_EqualAndIgnore(CudaNdarray *cnda1, CudaNdarray *cnda2, int ignoreSync, int ignoreBase)
51215115 {
51225116 int verbose = 0;
51235117
51245118 if (!ignoreSync && cnda1->dev_structure_fresh != cnda2->dev_structure_fresh)
51255119 {
51265120 if(verbose) fprintf(stdout, "CUDANDARRAY_EQUAL FAILED : 1\n");
51275121 return 0;
51285122 }
51295123
51305124 if (cnda1->nd != cnda2->nd)
51315125 {
51325126 if(verbose) fprintf(stdout, "CUDANDARRAY_EQUAL FAILED : 2\n");
51335127 return 0;
51345128 }
51355129
51365130 for (int i=0; i < 2*cnda1->nd; i++)
51375131 {
51385132 if (cnda1->host_structure[i] != cnda2->host_structure[i])
51395133 {
51405134 if(verbose)
51415135 fprintf(stdout, "CUDANDARRAY_EQUAL : host_structure : %d, %d, %d\n", i, cnda1->host_structure[i], cnda2->host_structure[i]);
51425136 return 0;
51435137 }
51445138 }
51455139
51465140 if (!ignoreBase && cnda1->base != cnda2->base)
51475141 {
51485142 if(verbose) fprintf(stdout, "CUDANDARRAY_EQUAL FAILED : 4");
51495143 return 0;
51505144 }
51515145 else if (cnda1->data_allocated != cnda2->data_allocated)
51525146 {
51535147 if(verbose) fprintf(stdout, "CUDANDARRAY_EQUAL FAILED : 5");
51545148 return 0;
51555149 }
51565150 else if (cnda1->data_allocated && cnda1->devdata != cnda2->devdata)
51575151 {
51585152 if(verbose) fprintf(stdout, "CUDANDARRAY_EQUAL FAILED : 6");
51595153 // no need to check devdata if data is not allocated
51605154 return 0;
51615155 }
51625156
51635157 return 1;
51645158 }
51655159
51665160
51675161 int
51685162 CudaNdarray_Equal(CudaNdarray *cnda1, CudaNdarray *cnda2)
51695163 {
51705164 return CudaNdarray_EqualAndIgnore(cnda1, cnda2, 0, 0);
51715165 }
51725166
51735167 int
51745168 cnda_copy_structure_to_device(const CudaNdarray * self)
51755169 {
51765170 //If the device structure do not exists, create it.
51775171 //We allocate it here as we do not need it often.
51785172 //In fact, we need it so infrequently that we expect
51795173 //that most object won't need it. Not allocating it
51805174 //save a significant when creating object.
51815175 //This speed up a benchmark by 8% with the gc.
51825176 if (!self->dev_structure)
51835177 {
51845178 int struct_size = cnda_structure_size(self->nd);
51855179 if (struct_size)
51865180 {
51875181 self->dev_structure = (int*)device_malloc(struct_size* sizeof(int));
51885182 if (NULL == self->dev_structure)
51895183 {
51905184 return -1;
51915185 }
51925186 }
51935187 }
51945188 if (cublasSetVector(cnda_structure_size(self->nd),
51955189 sizeof(int),
51965190 self->host_structure,
51975191 1,
51985192 self->dev_structure,
51995193 1) != CUBLAS_STATUS_SUCCESS)
52005194 {
52015195 PyErr_SetString(PyExc_RuntimeError, "error copying structure to device memory");
52025196 return -1;
52035197 }
52045198 self->dev_structure_fresh = 1;
52055199 return 0;
52065200 }
52075201
52085202 const int *
52095203 CudaNdarray_DEV_DIMS(const CudaNdarray * self)
52105204 {
52115205 if (!self->dev_structure_fresh)
52125206 {
52135207 if (cnda_copy_structure_to_device(self))
52145208 return NULL;
52155209 }
52165210 return self->dev_structure;
52175211 }
52185212 const int *
52195213 CudaNdarray_DEV_STRIDES(const CudaNdarray * self)
52205214 {
52215215 if (!self->dev_structure_fresh)
52225216 {
52235217 if (cnda_copy_structure_to_device(self))
52245218 return NULL;
52255219 }
52265220 return self->dev_structure + self->nd;
52275221 }
52285222 const int *
52295223 CudaNdarray_DEV_LOG2DIMS(const CudaNdarray * self)
52305224 {
52315225 if (!self->dev_structure_fresh)
52325226 {
52335227 if (cnda_copy_structure_to_device(self))
52345228 return NULL;
52355229 }
52365230 return self->dev_structure + 2*self->nd;
52375231 }
52385232 float *
52395233 CudaNdarray_DEV_DATA(const CudaNdarray * self)
52405234 {
52415235 return self->devdata;
52425236 }
52435237
52445238 /**
52455239 * Return the number of elements in the ndarray (product of the dimensions)
52465240 */
52475241 size_t
52485242 CudaNdarray_SIZE(const CudaNdarray *self)
52495243 {
52505244 if (self->nd == -1) return 0;
52515245 size_t size = 1;
52525246 for (int i = 0; i < self->nd; ++i)
52535247 {
52545248 size *= CudaNdarray_HOST_DIMS(self)[i];
52555249 }
52565250 return size;
52575251 }
52585252
52595253 PyObject *
52605254 CudaNdarray_SIZE_Object(const CudaNdarray *self, void *closure)
52615255 {
52625256 return PyInt_FromLong(CudaNdarray_SIZE(self));
52635257 }
52645258
52655259 int CudaNdarray_set_device_data(CudaNdarray * self, float * data, const CudaNdarray * base)
52665260 {
52675261 return CudaNdarray_set_device_data(self, data, (PyObject *) base);
52685262 }
52695263
52705264 PyObject * CudaNdarray_IS_C_Contiguous(CudaNdarray * self)
52715265 {
52725266 return PyBool_FromLong(CudaNdarray_is_c_contiguous(self));
52735267 }
52745268
52755269 int fprint_CudaNdarray(FILE * fd, const CudaNdarray *self)
52765270 {
52775271 cudaError_t err = cudaGetLastError();
52785272 if( cudaSuccess != err)
52795273 {
52805274 PyErr_Format(PyExc_RuntimeError,
52815275 "Cuda error: %s: %s.",
52825276 "fprint_CudaNdarray was called with an uncleared error",
52835277 cudaGetErrorString(err));
52845278 return -1;
52855279 }
52865280 fprintf(fd, "CudaNdarray <%p, %p> nd=%i dev_structure_fresh=%d data_allocated=%d\n",
52875281 self, self->devdata, self->nd, self->dev_structure_fresh, self->data_allocated);
52885282 fprintf(fd, "\tHOST_DIMS: ");
52895283 for (int i = 0; i < self->nd; ++i)
52905284 {
52915285 fprintf(fd, "%i\t", CudaNdarray_HOST_DIMS(self)[i]);
52925286 }
52935287 fprintf(fd, "\n\tHOST_STRIDES: ");
52945288 for (int i = 0; i < self->nd; ++i)
52955289 {
52965290 fprintf(fd, "%i\t", CudaNdarray_HOST_STRIDES(self)[i]);
52975291 }
52985292
52995293 if (self->dev_structure)
53005294 {
53015295 int data=0;
53025296 fprintf(fd, "\n\tDEV_DIMS: ");
53035297 for (int i = 0; i < self->nd; ++i)
53045298 {
53055299 cublasGetVector(1, sizeof(int),
53065300 self->dev_structure+i, 1,
53075301 &data, 1);
53085302 fprintf(fd, "%i\t", data);
53095303 }
53105304 fprintf(fd, "\n\tDEV_STRIDES: ");
53115305 for (int i = 0; i < self->nd; ++i)
53125306 {
53135307 cublasGetVector(1, sizeof(int),
53145308 self->dev_structure + self->nd+i, 1,
53155309 &data, 1);
53165310 fprintf(fd, "%i \t", data);
53175311 }
53185312 fprintf(fd, "\n");
53195313 }
53205314 else
53215315 {
53225316 fprintf(fd, "\n\tdev_structure not allocated\n");
53235317 }
53245318
53255319 err = cudaGetLastError();
53265320 if( cudaSuccess != err)
53275321 {
53285322 PyErr_Format(PyExc_RuntimeError,
53295323 "Cuda error: %s: %s.",
53305324 "fprint_CudaNdarray",
53315325 cudaGetErrorString(err));
53325326 return -1;
53335327 }
53345328 return 0;
53355329 }
53365330
53375331
53385332 int CudaNdarray_prep_output(CudaNdarray ** arr, int nd,
53395333 const int * dims, int fortran)
53405334 {
53415335 bool allocated = false;
53425336 if (*arr == NULL)
53435337 {
53445338 // This allocates the metadata but not the data
53455339 *arr = (CudaNdarray *) CudaNdarray_new_nd(nd);
53465340 if (*arr == NULL)
53475341 return -1;
53485342 allocated = true;
53495343 }
53505344
53515345 if (CudaNdarray_alloc_contiguous(*arr, nd, dims, fortran))
53525346 {
53535347 if (allocated)
53545348 {
53555349 Py_DECREF(*arr);
53565350 *arr = NULL;
53575351 }
53585352 return -1;
53595353 }
53605354 return 0;
53615355 }
53625356
53635357
53645358 /*
53655359 Local Variables:
53665360 mode:c++
53675361 c-basic-offset:4
53685362 c-file-style:"stroustrup"
53695363 indent-tabs-mode:nil
53705364 fill-column:79
53715365 End:
53725366 */
53735367 // vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:textwidth=79 :
53745368
5375===============================
5376nvcc 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).
5377nvcc fatal : The version ('80000') of the host compiler ('Apple clang') is not supported
5378
5379['nvcc', '-shared', '-O3', '-m64', '-Xcompiler', '-DCUDA_NDARRAY_CUH=mc72d035fdf91890f3b36710688069b2e,-DNPY_NO_DEPRECATED_API=NPY_1_7_API_VERSION,-fPIC,-fvisibility=hidden', '-Xlinker', '-rpath,/Users/burg/.theano/compiledir_Darwin-16.1.0-x86_64-i386-64bit-i386-3.5.2-64/cuda_ndarray', '-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_1/Frameworks/Python.framework/Versions/3.5/include/python3.5m', '-I/usr/local/lib/python3.5/site-packages/theano/gof', '-o', '/Users/burg/.theano/compiledir_Darwin-16.1.0-x86_64-i386-64bit-i386-3.5.2-64/cuda_ndarray/cuda_ndarray.so', 'mod.cu', '-L/usr/local/Cellar/python3/3.5.2_1/Frameworks/Python.framework/Versions/3.5/lib', '-lcublas', '-lcudart', '-Xcompiler', '-undefined,dynamic_lookup', '-Xlinker', '-pie']
5380ERROR (theano.sandbox.cuda): Failed to compile cuda_ndarray.cu: ('nvcc return status', 1, 'for cmd', 'nvcc -shared -O3 -m64 -Xcompiler -DCUDA_NDARRAY_CUH=mc72d035fdf91890f3b36710688069b2e,-DNPY_NO_DEPRECATED_API=NPY_1_7_API_VERSION,-fPIC,-fvisibility=hidden -Xlinker -rpath,/Users/burg/.theano/compiledir_Darwin-16.1.0-x86_64-i386-64bit-i386-3.5.2-64/cuda_ndarray -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_1/Frameworks/Python.framework/Versions/3.5/include/python3.5m -I/usr/local/lib/python3.5/site-packages/theano/gof -o /Users/burg/.theano/compiledir_Darwin-16.1.0-x86_64-i386-64bit-i386-3.5.2-64/cuda_ndarray/cuda_ndarray.so mod.cu -L/usr/local/Cellar/python3/3.5.2_1/Frameworks/Python.framework/Versions/3.5/lib -lcublas -lcudart -Xcompiler -undefined,dynamic_lookup -Xlinker -pie')
5381Setup training
5382Start training
5383Training status changed to STARTING
5384Imported model_provider in /Users/burg/git/nexar/keras/aetros-cli-data/networks/burgalon/digit-convolution/1xv2KAPRD/model_provider.py
5385Training status changed to LOAD DATA
5386Imported dataset provider in /Users/burg/git/nexar/keras/aetros-cli-data/networks/burgalon/digit-convolution/1xv2KAPRD/datasets/burgalon__dataset__mnist-digits.py
5387X_train shape: (60000, 1, 28, 28)
538860000 train samples
538910000 test samples
5390trainer.input_shape = []
5391trainer.classes = []
5392Possible data keys 'burgalon/dataset/mnist-digits'
5393Training status changed to CONSTRUCT
5394input_shapes [(None, 1, 28, 28)]
5395input_shapes [(None, 64, 26, 26)]
5396input_shapes [(None, 64, 13, 13)]
5397input_shapes [(None, 64, 11, 11)]
5398input_shapes [(None, 64, 5, 5)]
5399input_shapes [(None, 80, 3, 3)]
5400input_shapes [(None, 80, 1, 1)]
5401input_shapes [(None, 80)]
5402input_shapes [(None, 64)]
5403input_shapes [(None, 64)]
5404Training status changed to COMPILING
5405____________________________________________________________________________________________________
5406Layer (type) Output Shape Param # Connected to
5407====================================================================================================
5408Input (InputLayer) (None, 1, 28, 28) 0
5409____________________________________________________________________________________________________
5410conv1 (Convolution2D) (None, 64, 26, 26) 640 Input[0][0]
5411____________________________________________________________________________________________________
5412pool1 (MaxPooling2D) (None, 64, 13, 13) 0 conv1[0][0]
5413____________________________________________________________________________________________________
5414Node_9 (Convolution2D) (None, 64, 11, 11) 36928 pool1[0][0]
5415____________________________________________________________________________________________________
5416Node_10 (MaxPooling2D) (None, 64, 5, 5) 0 Node_9[0][0]
5417____________________________________________________________________________________________________
5418conv2 (Convolution2D) (None, 80, 3, 3) 46160 Node_10[0][0]
5419____________________________________________________________________________________________________
5420pool2 (MaxPooling2D) (None, 80, 1, 1) 0 conv2[0][0]
5421____________________________________________________________________________________________________
5422flatten (Flatten) (None, 80) 0 pool2[0][0]
5423____________________________________________________________________________________________________
5424dense (Dense) (None, 64) 5184 flatten[0][0]
5425____________________________________________________________________________________________________
5426dense_dropout (Dropout) (None, 64) 0 dense[0][0]
5427____________________________________________________________________________________________________
5428out (Dense) (None, 10) 650 dense_dropout[0][0]
5429====================================================================================================
5430Total params: 89562
5431____________________________________________________________________________________________________
5432Training status changed to TRAINING
5433Epoch 1: loss=0.332328, acc=0.897267, val_loss=0.109936, val_acc=0.964800
5434
5435Epoch 2: loss=0.141836, acc=0.959633, val_loss=0.077582, val_acc=0.976300
5436
5437
5438
5439Epoch 3: loss=0.107202, acc=0.969700, val_loss=0.058112, val_acc=0.982500
5440
5441
5442
5443
5444
5445Epoch 4: loss=0.088381, acc=0.974533, val_loss=0.059488, val_acc=0.982600
5446Epoch 5: loss=0.076564, acc=0.978117, val_loss=0.055403, val_acc=0.985100
5447Epoch 6: loss=0.068382, acc=0.980733, val_loss=0.052204, val_acc=0.984900
5448Epoch 7: loss=0.061393, acc=0.982967, val_loss=0.062597, val_acc=0.984100